In [24]:
%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

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

In [26]:
num_cutouts = len(x_list)

In [27]:
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')

In [28]:
cutout = [0,]*num_cutouts

for i in range(num_cutouts):
    cutout[i] = shapely_cutout(x_array[i],y_array[i])

In [29]:
cutout_dimensions = pd.DataFrame.from_csv('../data/cutout_dimensions.csv')
#cutout_dimensions

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

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

In [32]:
plot(append(width,width),append(log2(ratio),log2(1/ratio)),'.')


Out[32]:
[<matplotlib.lines.Line2D at 0x7fb87d324748>]

In [334]:
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)

In [313]:
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:,:]))

In [314]:
contour(xMesh,yMesh,gapMesh, levels=[130], colors='r',alpha=0.7)
contour(xMesh,yMesh,giveMesh, levels=[0.15], colors='g',alpha=0.7)


Out[314]:
<matplotlib.contour.QuadContourSet at 0x7fb87b8e09e8>

In [315]:
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')


Out[315]:
<matplotlib.text.Text at 0x7fb87b766e48>

In [316]:
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)

In [317]:
init_area = give_area(ones(len(width)))

In [318]:
d = rand(len(width))*0.2+0.8
give_area(d)/give_area(ones(len(width)))


Out[318]:
0.94444444444444453

In [330]:
give_sigmoid_bias(append(ones(len(width)-1),0))


Out[330]:
nan

In [333]:
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)


Out[333]:
array(0.0)

In [320]:
give_sigmoid_bias(d)


Out[320]:
48.919278811184718

In [321]:
# rand_weights = zeros(100)

# for i in range(100):
#     rand_weights[i] = give_sigmoid_bias(rand(len(width)))
    
# mean(rand_weights)


Out[321]:
65.355469460037824

In [336]:
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)

In [326]:
bounds = ([0,1],)*len(width)

In [ ]:
minimizer_kwargs = {"method": 'L-BFGS-B', "bounds": bounds  }
x0 = np.ones(len(width))

output = basinhopping(toMinimise,x0,
                      niter=1,minimizer_kwargs=minimizer_kwargs)

In [329]:
output


Out[329]:
                   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

In [307]:
t = linspace(0.9,1,1000)
y = exp((-(t-0.95)*100))/10
plot(t,y)
ylim([0,1])


Out[307]:
(0, 1)

In [206]:



Out[206]:
0.7189271775207503

In [118]:
sum(giveMeshHlf <= 0.15) / size(giveMeshHlf)


Out[118]:
0.28253968253968254

In [119]:
sum(giveMesh <= 0.15) / size(giveMesh)


Out[119]:
0.2742857142857143

In [ ]:


In [104]:
width_scale = std(width)

In [105]:
log2ratio_scale = std(append(log2ratio,-log2ratio))

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


Out[106]:

In [107]:
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')


Out[107]:
<matplotlib.text.Text at 0x7fb87c46ad30>

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


Out[108]:

In [109]:
shapely_valid = shapely_gap.intersection(shapely_give)
shapely_valid


Out[109]:

In [110]:
shapely_valid.area


Out[110]:
13.429001009871161

In [111]:
# def valid_area()

In [112]:
p = give_contour.collections[0].get_paths()[0]
v = p.vertices
x = v[:,0]
y = v[:,1]

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