In [1]:
%pylab
import csv
import pandas as pd
from scipy.interpolate import SmoothBivariateSpline
from scipy.optimize import minimize
from mpl_toolkits.mplot3d import Axes3D
In [2]:
def single_angle_gap(xTest,yTest,xData,yData,xScale,yScale):
xi = xScale * (xData - xData.min()) / xData.ptp()
yi = yScale * (yData - yData.min()) / yData.ptp()
xVal = xScale * (xTest - xData.min()) / xData.ptp()
yVal = yScale * (yTest - yData.min()) / yData.ptp()
dx = xi - xVal
dy = yi - yVal
if any((dx == 0) & (dy == 0)):
gap = 0
return gap
theta = np.arctan(dy/dx)
theta[dx<0] = theta[dx<0] + np.pi
theta[(dx>0) & (dy<0)] = theta[(dx>0) & (dy<0)] + 2*np.pi
theta[(dx==0) & (dy>0)] = np.pi/2
theta[(dx==0) & (dy<0)] = 3*np.pi/2
test = np.sort(theta)
test = np.append(test,test[0] + 2*np.pi)
gap = np.max(np.diff(test))*180/np.pi
return gap
def angle_gap(xTest,yTest,xData,yData,xScale,yScale):
dim = np.core.fromnumeric.shape(xTest)
if np.size(dim) == 0:
gap = single_angle_gap(xTest,yTest,xData,yData,xScale,yScale)
return gap
gap = np.zeros(dim)
if np.size(dim) == 1:
for i in range(dim[0]):
gap[i] = single_angle_gap(xTest[i],yTest[i],xData,yData,xScale,yScale)
return gap
for i in range(dim[0]):
for j in range (dim[1]):
gap[i,j] = single_angle_gap(xTest[i,j],yTest[i,j],xData,yData,xScale,yScale)
return gap
In [3]:
def single_fit_give(xTest,yTest,xData,yData,zData,s=None):
initialFitReturn = SmoothBivariateSpline(xData,yData,zData,kx=3,ky=3,s=s).ev(xTest,yTest)
adjXData = np.append(xData,xTest)
adjYData = np.append(yData,yTest)
posAdjZData = np.append(zData,initialFitReturn + 1)
negAdjZData = np.append(zData,initialFitReturn - 1)
posFitReturn = SmoothBivariateSpline(adjXData,adjYData,posAdjZData,kx=3,ky=3,s=s).ev(xTest,yTest)
negFitReturn = SmoothBivariateSpline(adjXData,adjYData,negAdjZData,kx=3,ky=3,s=s).ev(xTest,yTest)
posGive = posFitReturn - initialFitReturn
negGive = initialFitReturn - negFitReturn
give = np.mean([posGive, negGive])
return give
def fit_give(xTest,yTest,xData,yData,zData,s=None):
dim = np.core.fromnumeric.shape(xTest)
if np.size(dim) == 0:
give = single_fit_give(xTest,yTest,xData,yData,zData,s=s)
return give
give = np.zeros(dim)
if np.size(dim) == 1:
for i in range(dim[0]):
give[i] = single_fit_give(xTest[i],yTest[i],xData,yData,zData,s=s)
return give
for i in range(dim[0]):
for j in range (dim[1]):
give[i,j] = single_fit_give(xTest[i,j],yTest[i,j],xData,yData,zData,s=s)
return give
In [4]:
pylab.rcParams['savefig.dpi'] = 130
In [5]:
# Define thresholds
giveThreshold = 0.25
gapThreshold = 130
In [6]:
with open('cutout_factors','r') as f:
loaded_factor = eval(f.read())
dictArray = array(list(loaded_factor.items()),float)
cutoutRef = dictArray[:,0].astype(int)
ref = argsort(cutoutRef)
factor = dictArray[:,1]
cutoutRef = cutoutRef[ref]
factor = factor[ref]
with open('circle_cutout_factors','r') as f:
loaded_circle_factor = eval(f.read())
circleDictArray = array(list(loaded_circle_factor.items()),float)
circleDiameter = circleDictArray[:,0]
circleFactor = circleDictArray[:,1]
factor = append(factor,circleFactor)
with open('custom_cutout_factors','r') as f:
loaded_custom_factor = eval(f.read())
customDimensionsDict = {'5x13':[5,13],
'5x10':[5,10],
'5x8':[5,8],
'4x13':[4,13],
'4x10':[4,10],
'4x8':[4,8],
'4x6.5':[4,6.5],
'3x13':[3,13],
'3x9':[3,9],
'3x6.5':[3,6.5],
'3x5':[3,5]}
customDataLabelList = list(customDimensionsDict.keys())
custom_factors = [loaded_custom_factor[x] for x in customDataLabelList]
custom_widths = [customDimensionsDict[x][0] for x in customDataLabelList]
custom_lengths = [customDimensionsDict[x][1] for x in customDataLabelList]
factor = append(factor,custom_factors)
In [7]:
# loaded_factor.items()
In [8]:
cutoutRef
Out[8]:
In [8]:
In [9]:
cutout_dimensions = pd.DataFrame.from_csv('new_cutoutdata.csv',index_col =None)
ref = argsort(cutout_dimensions['cutoutRef'].values.astype(int))
width = cutout_dimensions['new_width'].values[ref]
length = cutout_dimensions['new_length'].values[ref]
weight = concatenate((ones(shape(width)), 1*ones(shape(circleDiameter)), ones(shape(custom_widths))))
width = append(width,circleDiameter)
length = append(length,circleDiameter)
nanFill = empty(shape(circleDiameter))
nanFill.fill(nan)
cutoutRef = append(cutoutRef,nanFill)
width = append(width,custom_widths)
length = append(length,custom_lengths)
nanFill = empty(shape(custom_widths))
nanFill.fill(nan)
cutoutRef = append(cutoutRef,nanFill)
# editRef = cutoutRef == 14
# width[editRef] = 4.3
# length[editRef] = 4.6
ratio = width / length
In [10]:
circleDiameter
Out[10]:
In [11]:
scatter(circleDiameter,circleFactor)
Out[11]:
In [12]:
scatter(width,ratio)
title('Data so far')
xlabel('Width (cm)')
ylabel('Ratio')
Out[12]:
In [13]:
# Spline fitting
extendedWidth = np.append(width,length)
extendedLength = np.append(length,width)
extendedFactor = np.append(factor,factor)
splineFit = SmoothBivariateSpline(extendedWidth,extendedLength,extendedFactor,kx=3,ky=3,s=len(width))
In [14]:
# Create surface mesh
xVec2 = linspace(width.min(),length.max(),100)
yVec2 = xVec2
xMesh2, yMesh2 = meshgrid(xVec2, yVec2)
gapMesh2 = angle_gap(xMesh2,yMesh2,extendedWidth,extendedLength,1,1)
giveMesh2 = fit_give(xMesh2,yMesh2,extendedWidth,extendedLength,extendedFactor,s=len(width))
zMesh2 = splineFit.ev(xMesh2, yMesh2)
ref2 = (gapMesh2>gapThreshold) | (giveMesh2>giveThreshold)
zMesh2[ref2] = nan
In [15]:
# Plot the fit
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
bot = np.amin(factor) - 0.04 * np.ptp(factor)
top = np.amax(factor) + 0.04 * np.ptp(factor)
ax.contour(xMesh2,yMesh2,gapMesh2, offset=bot, levels=[gapThreshold], colors='r',alpha=0.3)
ax.contour(xMesh2,yMesh2,giveMesh2, offset=bot, levels=[giveThreshold], colors='g',alpha=0.3)
proxyGap, = ax.plot([width.min(),width.min()],[width.min(),width.min()],[bot,bot],'r-',alpha=0.3)
proxyGive, = ax.plot([width.min(),width.min()],[width.min(),width.min()],[bot,bot],'g-',alpha=0.3)
ax.plot_surface(xMesh2,yMesh2,zMesh2,alpha=0.3,rstride=4,cstride=4)
sc = ax.scatter(width,length,factor)
ax.set_title('Fit displayed against width and aspect ratio')
ax.set_xlabel('Width')
ax.set_ylabel('Length')
ax.set_zlabel('Factor')
legend([proxyGap,proxyGive],['Gap threshold', 'Give threshold'],loc=7,prop={'size':8})
# right = ratio.min()-0.01*ratio.ptp()
# ax.set_ylim(right,1)
ax.set_zlim(bot,top)
ax.view_init(elev=10, azim=-80)
In [16]:
# Remove data and predict
predictedFactor = zeros(shape(width))
predictionAngleGap = zeros(shape(width))
predictionFitGive = zeros(shape(width))
for i in range(len(width)):
x = width[i]
y = length[i]
adjWidth = np.delete(width,i)
adjLength = np.delete(length,i)
adjFactor = np.delete(factor,i)
adjExtendedWidth = np.append(adjWidth,adjLength)
adjExtendedLength = np.append(adjLength,adjWidth)
adjExtendedFactor = np.append(adjFactor,adjFactor)
adjSplineFit = SmoothBivariateSpline(adjExtendedWidth,adjExtendedLength,adjExtendedFactor,kx=3,ky=3,s=len(width))
predictedFactor[i] = adjSplineFit.ev(x,y)
predictionAngleGap[i] = angle_gap(x,y,adjExtendedWidth,adjExtendedLength,1,1)
predictionFitGive[i] = fit_give(x,y,adjExtendedWidth,adjExtendedLength,adjExtendedFactor,s=len(width))
difference = factor - predictedFactor
withinThresholds = ((predictionAngleGap<gapThreshold) &
(predictionFitGive<giveThreshold)
)
In [17]:
# Display table
df = pd.DataFrame(transpose([cutoutRef,width, length, factor, predictedFactor, difference,
predictionAngleGap, predictionFitGive, withinThresholds]))
df.columns = ['Reference','width', 'length', 'factor', 'Predicted Factor', 'Difference',
'Angle Gap', 'Fit Give', 'Within thresholds']
# df.to_csv(path_or_buf = "cutoutdata.csv", index=False)
df
Out[17]:
In [18]:
print("There are", sum(withinThresholds), "points within give and gap thresholds."+
"\nThe standard deviation of the prediction differences for these is", std(difference[withinThresholds]))
In [19]:
# Histogram
fig = plt.figure()
ax = fig.add_subplot(111)
n, bins, patches = hist(difference[withinThresholds])
setp(patches, 'facecolor', 'g', 'alpha', 0.5)
ax.set_xlabel('Prediction differences')
ax.set_ylabel('Frequency')
ax.set_title('Histogram of prediction differences')
Out[19]:
In [20]:
# Visual representation of differences
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(xMesh2,yMesh2,zMesh2,alpha=0.3,rstride=4,cstride=4)
diffColour = copy(difference)
diffColour[~withinThresholds] = nan
colourRange = nanmax(abs(diffColour))
sc = ax.scatter(width,length,factor,c=diffColour, vmin=-colourRange, vmax=colourRange,s=50)
colorbar(sc)
ax.set_title('Differences')
ax.set_xlabel('Width')
ax.set_ylabel('Length')
ax.set_zlabel('Factor')
# ax.set_ylim(right,1)
ax.set_zlim(bot,top)
ax.view_init(elev=10, azim=-80)
In [21]:
diffTest = abs(difference)
diffTest[~withinThresholds] = 0
testRef = argmax(diffTest)
cutoutRef[testRef]
Out[21]:
In [22]:
factor[testRef]
Out[22]:
In [23]:
scatter(width,length)
scatter(width[testRef],length[testRef],c='red')
Out[23]:
In [24]:
def to_minimise(cut_increase,ref):
test_width = width[ref] + cut_increase
test_length = length[ref] + cut_increase
test_ratio = test_width / test_length
return (factor[ref] - splineFit.ev(test_width,test_length))
In [25]:
output = minimize(to_minimise,0,args=(testRef,))
output
Out[25]:
In [ ]: