In [451]:
%pylab inline
import csv
import pandas as pd
from scipy.interpolate import SmoothBivariateSpline
from scipy.optimize import minimize
from scipy.stats import shapiro
from mpl_toolkits.mplot3d import Axes3D
from threshold_functions import *
In [452]:
pylab.rcParams['savefig.dpi'] = 130
In [453]:
# Define thresholds
giveThreshold = 0.5
gapThreshold = 130
In [454]:
with open('cutout_factors','r') as f:
loaded_factor = eval(f.read())
dictArray = array(list(loaded_factor.items()),float)
cutoutRef = dictArray[:,0].astype(int)
factor = dictArray[:,1]
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 [455]:
# loaded_factor.items()
In [456]:
cutout_dimensions = pd.DataFrame.from_csv('../data/cutout_dimensions.csv')
width = cutout_dimensions['width'].values[cutoutRef]
length = cutout_dimensions['length'].values[cutoutRef]
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 [457]:
shuffleArray = arange(len(width),dtype=int)
shuffle(shuffleArray)
numToKeep = 11
toKeepRef = zeros(len(width))
toKeepRef[shuffleArray[0:numToKeep]] = 1
toKeepRef = toKeepRef.astype(bool)
In [497]:
arrayPicker = array([29,28,27,26,34,40,41,0,1,2,4,5,7,37,32],dtype=int) #,30,25,24
toKeepRef = zeros(len(width))
toKeepRef[arrayPicker] = 1
toKeepRef = toKeepRef.astype(bool)
print(len(arrayPicker))
scatter(width[~toKeepRef],ratio[~toKeepRef],c="grey")
scatter(width[toKeepRef],ratio[toKeepRef],c="red")
title('Data this time round')
xlabel('Width (cm)')
ylabel('Ratio')
Out[497]:
In [498]:
# cutoutRef,width, length, factor, ratio
In [498]:
In [499]:
# Spline fitting
widthMirrored = np.append(width[toKeepRef],width[toKeepRef])
ratioLogMirrored = np.log2(np.append(width[toKeepRef]/length[toKeepRef],length[toKeepRef]/width[toKeepRef]))
factorMirrored = np.append(factor[toKeepRef],factor[toKeepRef])
weightMirrored = append(weight[toKeepRef],weight[toKeepRef])
splineFit = SmoothBivariateSpline(widthMirrored,ratioLogMirrored,factorMirrored,w=weightMirrored,s=len(width[toKeepRef]))
In [500]:
# Create surface mesh
xVec2 = linspace(width.min(),width.max(),100)
yVec2 = linspace(ratio.min(),ratio.max(),100)
xMesh2, yMesh2 = meshgrid(xVec2, yVec2)
gapMesh2 = angle_gap(xMesh2,log2(yMesh2),widthMirrored,ratioLogMirrored,1,2)
giveMesh2 = fit_give(xMesh2,log2(yMesh2),widthMirrored,ratioLogMirrored,factorMirrored,s=len(width))
zMesh2 = splineFit.ev(xMesh2, log2(yMesh2))
ref2 = (gapMesh2>gapThreshold) | (giveMesh2>giveThreshold)
zMesh2[ref2] = nan
In [516]:
fig = plt.figure()
ax = fig.add_subplot(111)
cf = ax.contourf(xMesh2,yMesh2,giveMesh2,60)
colorbar(cf, label='Fit give')
ax.contour(xMesh2,yMesh2,giveMesh2,levels=[giveThreshold],colors='k')
ax.scatter(width[~toKeepRef],ratio[~toKeepRef],c="grey")
ax.scatter(width[toKeepRef],ratio[toKeepRef],c="red")
ax.set_xlabel('Width')
ax.set_ylabel('Ratio')
ax.set_title('Contour plot of the fit give')
Out[516]:
In [519]:
fig = plt.figure()
ax = fig.add_subplot(111)
cf = ax.contourf(xMesh2,yMesh2,gapMesh2,60)
colorbar(cf, label='Angle gap')
ax.contour(xMesh2,yMesh2,gapMesh2,levels=[gapThreshold],colors='k')
ax.scatter(width[~toKeepRef],ratio[~toKeepRef],c="grey")
ax.scatter(width[toKeepRef],ratio[toKeepRef],c="red")
ax.set_xlabel('Width')
ax.set_ylabel('Ratio')
ax.set_title('Contour plot of the angle gap')
Out[519]:
In [502]:
# 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()],[ratio.min(),ratio.min()],[bot,bot],'r-',alpha=0.3)
proxyGive, = ax.plot([width.min(),width.min()],[ratio.min(),ratio.min()],[bot,bot],'g-',alpha=0.3)
ax.plot_surface(xMesh2,yMesh2,zMesh2,alpha=0.3,rstride=4,cstride=4)
ax.scatter(width[~toKeepRef],ratio[~toKeepRef],factor[~toKeepRef],c="grey")
ax.scatter(width[toKeepRef],ratio[toKeepRef],factor[toKeepRef],c="red")
ax.set_title('Fit displayed against width and aspect ratio')
ax.set_xlabel('Width')
ax.set_ylabel('Ratio')
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 [503]:
# 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 = log2(width[i]/length[i])
adjToKeepRef = toKeepRef.copy()
adjToKeepRef[i] = 0
adjWidth = width[adjToKeepRef]
adjLength = length[adjToKeepRef]
adjFactor = factor[adjToKeepRef]
adjWeight = weight[adjToKeepRef]
adjWidthMirrored = np.append(adjWidth,adjWidth)
adjRatioLogMirrored = np.log2(np.append(adjWidth/adjLength,adjLength/adjWidth))
adjFactorMirrored = np.append(adjFactor,adjFactor)
adjWeightMirrored = append(adjWeight,adjWeight)
adjSplineFit = SmoothBivariateSpline(adjWidthMirrored,adjRatioLogMirrored,adjFactorMirrored,w=adjWeightMirrored,s=len(adjWidth))
predictedFactor[i] = adjSplineFit.ev(x,y)
predictionAngleGap[i] = angle_gap(x,y,adjWidthMirrored,adjRatioLogMirrored,1,2)
predictionFitGive[i] = fit_give(x,y,adjWidthMirrored,adjRatioLogMirrored,adjFactorMirrored,s=len(adjWidth))
difference = factor - predictedFactor
withinThresholds = ((predictionAngleGap<gapThreshold) &
(predictionFitGive<giveThreshold)
)
In [504]:
# 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[504]:
In [505]:
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 [506]:
# 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[506]:
In [507]:
W, pvalue = shapiro(difference[withinThresholds])
pvalue
Out[507]:
In [508]:
# 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,ratio,factor,c=diffColour, vmin=-colourRange, vmax=colourRange)
colorbar(sc)
ax.set_title('Differences')
ax.set_xlabel('Width')
ax.set_ylabel('Ratio')
ax.set_zlabel('Factor')
ax.set_ylim(right,1)
ax.set_zlim(bot,top)
ax.view_init(elev=10, azim=-80)
In [509]:
diffTest = abs(difference)
diffTest[~withinThresholds] = 0
testRef = argmax(diffTest)
cutoutRef[testRef]
Out[509]:
In [510]:
factor[testRef]
Out[510]:
In [511]:
scatter(width,ratio)
scatter(width[testRef],ratio[testRef],c='red')
Out[511]:
In [512]:
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,log2(test_ratio)))**2
In [513]:
output = minimize(to_minimise,0,args=(testRef,))
output
Out[513]:
In [514]:
# scatter(width,ratio)
In [150]:
In [ ]: