In [2]:
%pylab inline

import csv
import pandas as pd

from scipy.interpolate import SmoothBivariateSpline
from scipy.optimize import minimize
from mpl_toolkits.mplot3d import Axes3D

from threshold_functions import *


Populating the interactive namespace from numpy and matplotlib

In [3]:
pylab.rcParams['savefig.dpi'] = 130

In [4]:
# Define thresholds
giveThreshold = 0.25
gapThreshold = 130

In [5]:
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 [6]:
# loaded_factor.items()

In [7]:
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 [8]:
circleDiameter


Out[8]:
array([ 9.,  4.,  3.,  7.,  5.,  6.,  8.])

In [9]:
scatter(circleDiameter,circleFactor)


Out[9]:
<matplotlib.collections.PathCollection at 0x9090940>

In [10]:
scatter(width,ratio)
title('Data so far')
xlabel('Width (cm)')
ylabel('Ratio')


Out[10]:
<matplotlib.text.Text at 0x947b860>

In [11]:
# Spline fitting

widthMirrored = np.append(width,width)
ratioLogMirrored = np.log2(np.append(width/length,length/width))
factorMirrored = np.append(factor,factor)
weightMirrored = append(weight,weight)

splineFit = SmoothBivariateSpline(widthMirrored,ratioLogMirrored,factorMirrored,w=weightMirrored,s=len(width))

In [12]:
# 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 [13]:
# 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)
sc = ax.scatter(width,ratio,factor)

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 [14]:
# 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])

    adjWidth = np.delete(width,i)
    adjLength = np.delete(length,i)
    adjFactor = np.delete(factor,i)
    adjWeight = delete(weight,i)
    
    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(width))
    
    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(width))
    
difference = factor - predictedFactor
withinThresholds = ((predictionAngleGap<gapThreshold) & 
                    (predictionFitGive<giveThreshold)
                    )

In [15]:
# 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[15]:
Reference width length factor Predicted Factor Difference Angle Gap Fit Give Within thresholds
0 6 7.245140 10.062915 1.001302 0.994350 0.006952 124.219522 0.148445 1
1 38 4.571007 5.636698 1.027190 1.032128 -0.004938 25.676602 0.075467 1
2 19 5.804896 10.176211 1.001089 1.005055 -0.003966 141.439273 0.180585 0
3 43 6.103849 7.736637 1.007783 1.006442 0.001340 48.859845 0.082024 1
4 41 6.882961 10.087318 0.994829 0.997270 -0.002441 139.513440 0.159651 0
5 3 4.536055 6.035738 1.032704 1.031558 0.001147 33.823656 0.087725 1
6 73 7.253007 9.401196 0.995472 0.996086 -0.000613 56.076632 0.134567 1
7 32 7.728437 10.890306 0.993330 0.993562 -0.000231 183.630498 0.607070 0
8 106 6.042542 9.643383 1.000871 1.003816 -0.002945 68.723386 0.127254 1
9 83 5.676411 7.903129 1.010552 1.011009 -0.000457 45.037387 0.078091 1
10 58 5.968616 7.249113 1.008999 1.008610 0.000389 28.441917 0.070398 1
11 70 7.461885 8.181920 0.997190 0.995502 0.001688 85.806264 0.138505 1
12 14 4.507535 4.917415 1.043027 1.032378 0.010648 31.086185 0.063607 1
13 104 5.578054 8.115864 1.006601 1.012040 -0.005439 37.862942 0.079914 1
14 109 6.478522 7.636965 1.006979 1.002617 0.004363 36.628409 0.089597 1
15 33 5.933829 6.141369 1.006789 1.010254 -0.003464 16.824736 0.050451 1
16 57 3.541777 4.655761 1.059704 1.054385 0.005319 125.010990 0.086374 1
17 34 7.095174 9.471365 0.993330 0.997021 -0.003691 71.785309 0.109425 1
18 16 6.176973 10.235482 0.999567 1.001357 -0.001790 145.702397 0.234242 0
19 22 5.216353 9.919577 1.011713 1.011827 -0.000114 125.617963 0.127913 1
20 18 7.637566 9.560354 0.998054 0.993977 0.004077 113.373515 0.249033 1
21 82 6.276570 7.584045 1.000867 1.005346 -0.004479 41.925218 0.083059 1
22 20 4.265897 5.041729 1.037967 1.038378 -0.000412 32.578995 0.071180 1
23 112 3.510689 4.141420 1.053904 1.056729 -0.002825 123.077411 0.074717 1
24 NaN 9.000000 9.000000 0.991661 0.995775 -0.004114 274.790017 0.747814 0
25 NaN 4.000000 4.000000 1.045826 1.044633 0.001193 54.104885 0.056470 1
26 NaN 3.000000 3.000000 1.075771 1.067972 0.007800 180.000000 0.210307 0
27 NaN 7.000000 7.000000 0.996776 0.999399 -0.002622 39.217156 0.092973 1
28 NaN 5.000000 5.000000 1.030396 1.023590 0.006806 24.415841 0.055600 1
29 NaN 6.000000 6.000000 1.014436 1.008536 0.005900 14.837542 0.050962 1
30 NaN 8.000000 8.000000 0.993360 0.993709 -0.000349 100.951895 0.132215 1
31 NaN 4.000000 10.000000 1.029933 1.032597 -0.002664 80.882662 0.238387 1
32 NaN 3.000000 9.000000 1.053516 1.057177 -0.003661 180.000000 0.298737 0
33 NaN 5.000000 8.000000 1.017153 1.019998 -0.002844 48.068865 0.093269 1
34 NaN 3.000000 5.000000 1.069714 1.066858 0.002856 180.000000 0.263698 0
35 NaN 3.000000 13.000000 1.049271 1.046207 0.003063 320.414144 0.886623 0
36 NaN 5.000000 10.000000 1.012277 1.014866 -0.002589 117.706883 0.138994 1
37 NaN 4.000000 6.500000 1.040668 1.041220 -0.000552 68.073634 0.139233 1
38 NaN 5.000000 13.000000 1.013374 0.998744 0.014630 176.308101 0.766304 0
39 NaN 4.000000 13.000000 1.024896 1.021441 0.003455 172.746432 0.799694 0
40 NaN 4.000000 8.000000 1.030055 1.038281 -0.008227 54.858700 0.175037 1
41 NaN 3.000000 6.500000 1.054423 1.065631 -0.011208 180.000000 0.313602 0

In [16]:
print("There are", sum(withinThresholds), "points within give and gap thresholds."+
      "\nThe standard deviation of the prediction differences for these is", std(difference[withinThresholds]))


There are 30 points within give and gap thresholds.
The standard deviation of the prediction differences for these is 0.00419634920895

In [17]:
# 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[17]:
<matplotlib.text.Text at 0x9420a20>

In [18]:
# 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 [19]:
diffTest = abs(difference)
diffTest[~withinThresholds] = 0

testRef = argmax(diffTest)
cutoutRef[testRef]


Out[19]:
14.0

In [20]:
factor[testRef]


Out[20]:
1.0430265009664661

In [21]:
scatter(width,ratio)
scatter(width[testRef],ratio[testRef],c='red')


Out[21]:
<matplotlib.collections.PathCollection at 0x947bef0>

In [22]:
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 [23]:
output = minimize(to_minimise,0,args=(testRef,))
output


Out[23]:
   status: 0
      jac: array([  7.06460591e-07])
     njev: 17
  success: True
  message: 'Optimization terminated successfully.'
      fun: 2.4108718840179435e-10
        x: array([-0.43756466])
 hess_inv: array([[ 1001.21558256]])
     nfev: 51

In [24]:
# scatter(width,ratio)

In [24]: