%pylab inline
import pandas as pd
import csv
from threshold_functions import *
from equivalent_ellipse import *
from scipy.interpolate import SmoothBivariateSpline, bisplrep
Populating the interactive namespace from numpy and matplotlib
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])
cutout_dimensions = pd.DataFrame.from_csv('../data/cutout_dimensions.csv')
#cutout_dimensions
boundBox = cutout_dimensions['box_bounds'].values
#ref = boundBox < 11
width = cutout_dimensions['width'].values
length = cutout_dimensions['length'].values
cutoutRef = arange(num_cutouts)
ratio = width/length
# ref = True
# while any(ref):
# sortingArray = array([width,ratio,cutoutRef]).transpose()
# sortingTable = pd.DataFrame(sortingArray)
# sortingTable = sortingTable.sort([0, 1])
# width = sortingTable[0].values
# ratio = sortingTable[1].values
# cutoutRef = sortingTable[2].values
# ref = (width[0:-1] == width[1:]) & (ratio[0:-1] == ratio[1:])
# width = delete(width,find(ref))
# ratio = delete(ratio,find(ref))
# cutoutRef = delete(cutoutRef,find(ref))
plot(append(width,width),append(log2(ratio),log2(1/ratio)),'.')
[<matplotlib.lines.Line2D at 0x7fb87d324748>]
def give_calc(xGrid,log2yGrid,width,ratio,w):
w = w/max(w)
ref = (w == 0)
width = delete(width,find(ref))
ratio = delete(ratio,find(ref))
w = delete(w,find(ref))
return fit_give(xGrid,log2yGrid,
append(width,width),
append(log2(ratio),log2(1/ratio)),
zeros(len(width)*2),append(w,w))
def gap_calc(xGrid,log2yGrid,width,ratio):
return angle_gap(xGrid,log2yGrid,append(width,width),
append(log2(ratio),log2(1/ratio)),1,2)
resltn = 11 # Must be odd
log2ratio = log2(ratio)
mirr_log2ratio = append(log2ratio,-log2ratio)
w = ones(len(width))
xRange = [mean(width)-3.5*std(width),mean(width)+3.5*std(width)]
yRange = [-3.5*std(mirr_log2ratio),0]
xVec = linspace(xRange[0],xRange[1],resltn)
yVec = linspace(yRange[0],0,(resltn//2+1))
xMeshHlf, yMeshHlf = meshgrid(xVec, yVec)
gapMeshHlf = gap_calc(xMeshHlf,yMeshHlf,width,ratio)
giveMeshHlf = give_calc(xMeshHlf,yMeshHlf,width,ratio,w)
xMesh = vstack((xMeshHlf,flipud(xMeshHlf)[1:,:]))
yMesh = vstack((yMeshHlf,-flipud(yMeshHlf)[1:,:]))
gapMesh = vstack((gapMeshHlf,flipud(gapMeshHlf)[1:,:]))
giveMesh = vstack((giveMeshHlf,flipud(giveMeshHlf)[1:,:]))
contour(xMesh,yMesh,gapMesh, levels=[130], colors='r',alpha=0.7)
contour(xMesh,yMesh,giveMesh, levels=[0.15], colors='g',alpha=0.7)
<matplotlib.contour.QuadContourSet at 0x7fb87b8e09e8>
fig = plt.figure()
ax = fig.add_subplot(111)
cf = ax.contourf(xMesh,yMesh,giveMesh,60)
give_contour = ax.contour(xMesh,yMesh,giveMesh,levels=[0.15],colors='k')
ax.scatter(append(width,width),append(log2ratio,-log2ratio),color='black')
colorbar(cf, label='Fit give')
ax.set_xlabel('Width')
ax.set_ylabel('log2(Ratio)')
ax.set_title('Contour plot of the fit give')
<matplotlib.text.Text at 0x7fb87b766e48>
def give_area(w):
giveMeshHlf = give_calc(xMeshHlf,yMeshHlf,width,ratio,w)
return sum(giveMeshHlf <= 0.15) / size(giveMeshHlf)
def give_sigmoid_bias(w):
giveMeshHlf = give_calc(xMeshHlf,yMeshHlf,width,ratio,w)
return mean(200/(1+exp(-(giveMeshHlf-0.15)*100)) - 100)
init_area = give_area(ones(len(width)))
d = rand(len(width))*0.2+0.8
give_area(d)/give_area(ones(len(width)))
0.94444444444444453
give_sigmoid_bias(append(ones(len(width)-1),0))
nan
w = append(ones(len(width)*2-1),0.0001)
SmoothBivariateSpline(append(width,width),append(log2ratio,-log2ratio),zeros(len(width)*2),w=w).ev(1,2)
array(0.0)
give_sigmoid_bias(d)
48.919278811184718
# rand_weights = zeros(100)
# for i in range(100):
# rand_weights[i] = give_sigmoid_bias(rand(len(width)))
# mean(rand_weights)
65.355469460037824
def toMinimise(w):
if (sum(w!=0) < 13):
bias = 100
else:
t = (give_area(w) / init_area)
bias = exp((-(t-0.95)*100))/10
return bias + mean(w)
bounds = ([0,1],)*len(width)
minimizer_kwargs = {"method": 'L-BFGS-B', "bounds": bounds }
x0 = np.ones(len(width))
output = basinhopping(toMinimise,x0,
niter=1,minimizer_kwargs=minimizer_kwargs)
output
nit: 1 x: array([ 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861944, 0.25861947, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041, 0.25862041]) message: ['requested number of basinhopping iterations completed successfully'] minimization_failures: 0 fun: 0.25929367093269146 nfev: 1298
t = linspace(0.9,1,1000)
y = exp((-(t-0.95)*100))/10
plot(t,y)
ylim([0,1])
(0, 1)
0.7189271775207503
sum(giveMeshHlf <= 0.15) / size(giveMeshHlf)
0.28253968253968254
sum(giveMesh <= 0.15) / size(giveMesh)
0.2742857142857143
width_scale = std(width)
log2ratio_scale = std(append(log2ratio,-log2ratio))
give_xylist = give_contour.collections[0].get_paths()[0].vertices
shapely_give = sh.Polygon([(i[0], i[1]) for i in
zip(give_xylist[:,0]/width_scale,
give_xylist[:,1]/log2ratio_scale)])
shapely_give
fig = plt.figure()
ax = fig.add_subplot(111)
cf = ax.contourf(xMesh,yMesh,gapMesh,60)
gap_contour = ax.contour(xMesh,yMesh,gapMesh,levels=[130],colors='k')
ax.scatter(append(width,width),append(log2ratio,-log2ratio),color='black')
colorbar(cf, label='Angle gap')
ax.set_xlabel('Width')
ax.set_ylabel('log2(Ratio)')
ax.set_title('Contour plot of the angle gap')
<matplotlib.text.Text at 0x7fb87c46ad30>
gap_xylist = gap_contour.collections[0].get_paths()[0].vertices
shapely_gap = sh.Polygon([(i[0], i[1]) for i in
zip(gap_xylist[:,0]/width_scale,
gap_xylist[:,1]/log2ratio_scale)])
shapely_gap
shapely_valid = shapely_gap.intersection(shapely_give)
shapely_valid
shapely_valid.area
13.429001009871161
# def valid_area()
p = give_contour.collections[0].get_paths()[0]
v = p.vertices
x = v[:,0]
y = v[:,1]
# with open('../data/optimised_cutout_x.csv', 'w') as x_csvfile:
# with open('../data/optimised_cutout_y.csv', 'w') as y_csvfile:
# x_writer = csv.writer(x_csvfile, delimiter=',', lineterminator='\n')
# y_writer = csv.writer(y_csvfile, delimiter=',', lineterminator='\n')
# for i in range(len(x_list_optimised)):
# x_writer.writerow(x_list_optimised[i])
# y_writer.writerow(y_list_optimised[i])