In [1]:
%pylab inline
import pandas as pd
import csv
from scipy.interpolate import SmoothBivariateSpline, bisplrep
from threshold_functions import *
from equivalent_ellipse import *
In [2]:
x_list = list()
y_list = list()
with open('../data/cutout_x.csv', 'r') as x_csvfile:
with open('../data/cutout_y.csv', 'r') as y_csvfile:
x_reader = csv.reader(x_csvfile, delimiter=',', lineterminator='\n')
y_reader = csv.reader(y_csvfile, delimiter=',', lineterminator='\n')
for row in x_reader:
x_list += [row]
for row in y_reader:
y_list += [row]
num_cutouts = len(x_list)
x_array = [0,]*num_cutouts
y_array = [0,]*num_cutouts
for i in range(num_cutouts):
x_array[i] = array(x_list[i], dtype='float')
y_array[i] = array(y_list[i], dtype='float')
cutout = [0,]*num_cutouts
for i in range(num_cutouts):
cutout[i] = shapely_cutout(x_array[i],y_array[i])
In [3]:
cutout_dimensions = pd.DataFrame.from_csv('../data/cutout_dimensions.csv')
In [4]:
ellipse_def = transpose(array([cutout_dimensions['a'].values,
cutout_dimensions['b'].values,
cutout_dimensions['width'].values,
cutout_dimensions['length'].values,
cutout_dimensions['theta'].values]))
In [5]:
width_init = cutout_dimensions['width'].values
cutoutRef_init = arange(len(width_init))
sorted_ref = cutoutRef_init[argsort(width_init)]
In [6]:
similar_list = diff(width_init[sorted_ref]) < 0.0001
duplicates_removed_ref = sorted_ref[~similar_list]
In [7]:
boundBox = cutout_dimensions['box_bounds'].values
too_big = boundBox[duplicates_removed_ref] > 11
too_big_removed_ref = duplicates_removed_ref[~too_big]
In [8]:
final_ref = too_big_removed_ref
In [9]:
width = cutout_dimensions['width'].values[final_ref]
length = cutout_dimensions['length'].values[final_ref]
cutoutRef = final_ref
ratio = width/length
In [10]:
# plot(append(width,width),append(log2(ratio),log2(1/ratio)),'.')
In [11]:
def give_calc(xGrid,yGrid,width,ratio,w):
w = w/max(w)
ref = (w < 0.01)
width = delete(width,find(ref))
ratio = delete(ratio,find(ref))
w = delete(w,find(ref))
return fit_give(xGrid,log2(yGrid),
append(width,width),
append(log2(ratio),log2(1/ratio)),
zeros(len(width)*2),append(w,w))
def gap_calc(xGrid,yGrid,width,ratio):
return angle_gap(xGrid,log2(yGrid),append(width,width),
append(log2(ratio),log2(1/ratio)),1,2)
In [12]:
# w = ones(len(width))
# valid = (give_calc(width,ratio,width,ratio,w) < 0.15) & (gap_calc(width,ratio,width,ratio) < 130)
# curr_give = give_calc(width[valid],ratio[valid],width,ratio,w)
# curr_gap = gap_calc(width[valid],ratio[valid],width,ratio)
# curr_gap
In [13]:
# curr_give = give_calc(width[valid],ratio[valid],width,ratio,w)
In [14]:
# curr_give[2] = 0.2
# sum((curr_give/0.15)**15)
In [15]:
# curr_gap[2] = 100
# sum((curr_gap/130)**30)
In [16]:
# sum(w)
In [17]:
# sort(rand(100))[0:-5]
In [18]:
global keeping_track
def toMinimise(w):
if (sum(w) > 13):
curr_give = give_calc(width,ratio,width,ratio,w)
else:
curr_give = ones(len(w))
weights_bias = sum((2*exp(-((w-0.5)*2)**2) + 3*w - 1)/3)
too_few_bias = exp(-(sum(w)-15))
give_bias_init = 4/(1 + exp(-(curr_give-0.15)*200))
give_bias = sum(sort(give_bias_init)[0:-5])
output = weights_bias + give_bias + too_few_bias
return output
In [19]:
# t = linspace(0,20,1000)
# y = exp(-(t-13))
# plot(t,y)
# ylim([0,1])
In [20]:
# t = linspace(0,1,1000)
# y = (2*exp(-((t-0.5)*2)**2) + 3*t - 1)/3
# plot(t,y)
In [21]:
# t = linspace(0,1,1000)
# y = 4/(1 + exp(-(t-0.15)*200))
# t5 = 0.15
# y5 = 4/(1 + exp(-(t5-0.15)*200))
# print(y5)
# plot(t,y)
# scatter(t5,y5)
In [22]:
def print_fun(w, f, accepted):
global successCount_basinhopping
global minVals_basinhopping
global nSuccess_basinhopping
print("Output %.4f" % (f))
w_show = reshape(w,[3,23])
imshow(w_show, interpolation='nearest', cmap=cm.gray, vmin=0, vmax=1)
show()
In [23]:
bounds = ([0,1],)*len(width)
In [24]:
# x0 = ones(len(width))
x0 = loadtxt('x0')
In [25]:
keeping_track = 0
global nSuccess_basinhopping
nSuccess_basinhopping = 5
global minVals_basinhopping
minVals_basinhopping = np.empty(nSuccess_basinhopping)
minVals_basinhopping[:] = np.nan
global successCount_basinhopping
successCount_basinhopping = 0
minimizer_kwargs = {"tol":0.01,"method": 'SLSQP', "bounds": bounds}
In [26]:
# curr_give = give_calc(width[valid],ratio[valid],width,ratio,ones(len(width)))
# curr_give
In [27]:
# give_calc(width[valid],ratio[valid],width,ratio,x0)
In [28]:
# (give_calc(width[valid],ratio[valid],width,ratio,ones(len(width))) -
# give_calc(width[valid],ratio[valid],width,ratio,x0))
In [ ]:
resltn = 25 # Must be odd
xRange = [min(width),max(width)]
yRange = [min(ratio),max(ratio)]
xVec = linspace(xRange[0],xRange[1],resltn)
yVec = linspace(yRange[0],yRange[1],resltn)
xMesh, yMesh = meshgrid(xVec, yVec)
giveMesh = give_calc(xMesh,yMesh,width,ratio,ones(len(width)))
In [ ]:
while 1:
output = basinhopping(toMinimise,x0,
niter=20,minimizer_kwargs=minimizer_kwargs,
callback=print_fun)
print("Final %.4f" % (output.fun))
w_show = reshape(output.x,[3,23])
imshow(w_show, interpolation='nearest', cmap=cm.gray, vmin=0, vmax=1)
show()
w = output.x/max(output.x)
# give_pic = give_calc(width,ratio,width,ratio,w)
# print(around(give_pic,4))
saveStuffs = w > 0.5
widthSave = width[saveStuffs]
ratioSave = ratio[saveStuffs]
plot(width[~saveStuffs],ratio[~saveStuffs],'rx')
scatter(widthSave,ratioSave)
print("Currently %d cutouts require measurement:" % (sum(saveStuffs)))
newGiveMesh = give_calc(xMesh,yMesh,width,ratio,w)
contour(xMesh,yMesh,giveMesh, levels=[0.05,0.1,0.15], colors='g',alpha=0.7)
contour(xMesh,yMesh,newGiveMesh, levels=[0.05,0.1,0.15], colors='r',alpha=0.7)
title('Optimisation of where to measure')
xlabel('Width (cm)')
ylabel('Ratio')
x0 = w
# print(output)
savetxt('cutouts_to_be_measured',cutoutRef[saveStuffs])
savetxt('x0',x0)
show()
In [ ]:
# x0 = output.x
# oldF = 200
In [ ]:
# savetxt('oldF',oldF)
In [ ]:
# newGiveMesh = give_calc(xMesh,yMesh,width,ratio,x0)
# contour(xMesh,yMesh,giveMesh, levels=[0.05,0.1,0.15], colors='g',alpha=0.7)
# contour(xMesh,yMesh,newGiveMesh, levels=[0.05,0.1,0.15], colors='r',alpha=0.7)
# len(width) - sum(w)
In [ ]:
# saveStuffs = x0 > 0.5
# widthSave = width[saveStuffs]
# ratioSave = ratio[saveStuffs]
# plot(width[~saveStuffs],ratio[~saveStuffs],'rx')
# scatter(widthSave,ratioSave)
# sum(saveStuffs)
In [ ]:
# oldF
In [ ]:
# output.fun
In [ ]: