Write your desired energy and applicator in the cell immediately below then from the "Cell" menu above select "Run All"
In [1]:
energyRequest = 15
applicatorRequest = 25
In [2]:
path2datafile = 'S:/Dosimetry/Elekta_EFacs/measured_cutout_factors.csv'
path2outputfile = 'S:/Dosimetry/Elekta_EFacs/%dMeV%dcm_report.html' %(energyRequest, applicatorRequest)
# path2datafile = 'archive/measured_cutout_factors.csv'
# path2outputfile = '%dMeV%dcm_report.html' %(energyRequest, applicatorRequest)
In [3]:
try:
import plotly as py
except:
!pip install plotly
import plotly as py
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import pylab
%matplotlib inline
from scipy.interpolate import UnivariateSpline
from scipy.interpolate import SmoothBivariateSpline
from scipy.stats import shapiro, ttest_1samp
import plotly.plotly as py
import plotly.tools as tls
from plotly.graph_objs import *
import os
from IPython.display import HTML
pylab.rcParams['savefig.dpi'] = 120
In [4]:
py.sign_in('SimonBiggs', 'oas177ahva')
if (os.path.isfile('key.txt')):
f = open('key.txt','r')
key = f.read()
f.close()
else:
key = "%04d-%04d-%04d-%04d-%04d" %tuple((np.floor(np.random.rand(5)*1e4)).astype(int))
f = open('key.txt','w')
f.write(key)
f.close()
In [5]:
measuredRawData = pd.DataFrame.from_csv(path2datafile,index_col=False)
In [6]:
energy = measuredRawData['energy'].values
applicator = measuredRawData['applicator'].values
label = measuredRawData['label'].values
width = measuredRawData['width'].values
length = measuredRawData['length'].values
cutoutFactor = measuredRawData['cutout factor'].values
wrongWayRef = width>length
widthTemp = width.copy()
width[wrongWayRef] = length[wrongWayRef]
length[wrongWayRef] = widthTemp[wrongWayRef]
ratio = width/length
eqPonA = 2*( 3*(ratio+1) - np.sqrt( (3*ratio+1)*(ratio+3) ) ) / width
In [7]:
ref = (energy == energyRequest) & (applicator == applicatorRequest)
In [8]:
marker_labels = [0,]*sum(ref)
for i in range(sum(ref)):
marker_labels[i] = ('Width: %0.1f cm'
+'<br>Length: %0.1f cm'
+'<br>Cutout factor: %0.3f '
+'<br>Label: %s') %(width[ref][i],
length[ref][i],
cutoutFactor[ref][i],
label[ref][i])
In [9]:
plotScatter = Scatter(x=width[ref], y=cutoutFactor[ref],
text=marker_labels,
mode='markers',
marker=Marker(size=15),
name='Measured data')
dataList = [plotScatter]
xSpline = np.linspace(width[ref].min(),width[ref].max(),100)
xSpline = np.unique(np.append(xSpline,width[ref]))
try:
spline1 = UnivariateSpline(width[ref],cutoutFactor[ref],k=1)
ySpline1 = spline1(xSpline)
plotSpline1 = Scatter(x=xSpline, y=ySpline1,
line=Line(width=2,opacity=0.5),
name='First order')
except:
print("Cannot create first order fit")
else:
dataList = dataList + [plotSpline1]
try:
spline2 = UnivariateSpline(width[ref],cutoutFactor[ref],k=2)
ySpline2 = spline2(xSpline)
plotSpline2 = Scatter(x=xSpline, y=ySpline2,
line=Line(width=2,opacity=0.5),
name='Second order')
except:
print("Cannot create second order fit")
else:
dataList = dataList + [plotSpline2]
data = Data(dataList)
layout = Layout(title=('First and second order univarite spline fits for '
+'the %d MeV and %d cm applicator') %(energyRequest, applicatorRequest),
xaxis=XAxis(title='Width (cm)'),
yaxis=YAxis(title='Cutout factor'),
hovermode='closest')
fig_width = Figure(data=data, layout=layout)
py.iplot(fig_width, filename='%s_electron-fits_widthalone_%dMeV_%dcm' %(key,energyRequest, applicatorRequest))
Out[9]:
In [10]:
def single_angle_gap(xTest,yTest,xData,yData,xScale,yScale):
xi = xScale * (xData - xData.min()) / xData.ptp()
yi = yScale * (yData - yData.min()) / yData.ptp()
xVal = xScale * (xTest - xData.min()) / xData.ptp()
yVal = yScale * (yTest - yData.min()) / yData.ptp()
dx = xi - xVal
dy = yi - yVal
if any((dx == 0) & (dy == 0)):
gap = 0
return gap
theta = np.arctan(dy/dx)
theta[dx<0] = theta[dx<0] + np.pi
theta[(dx>0) & (dy<0)] = theta[(dx>0) & (dy<0)] + 2*np.pi
theta[(dx==0) & (dy>0)] = np.pi/2
theta[(dx==0) & (dy<0)] = 3*np.pi/2
test = np.sort(theta)
test = np.append(test,test[0] + 2*np.pi)
gap = np.max(np.diff(test))*180/np.pi
return gap
def angle_gap(xTest,yTest,xData,yData,xScale,yScale):
dim = np.core.fromnumeric.shape(xTest)
if np.size(dim) == 0:
gap = single_angle_gap(xTest,yTest,xData,yData,xScale,yScale)
return gap
gap = np.zeros(dim)
if np.size(dim) == 1:
for i in range(dim[0]):
gap[i] = single_angle_gap(xTest[i],yTest[i],xData,yData,xScale,yScale)
return gap
for i in range(dim[0]):
for j in range (dim[1]):
gap[i,j] = single_angle_gap(xTest[i,j],yTest[i,j],xData,yData,xScale,yScale)
return gap
In [11]:
def single_fit_give(xTest,yTest,xData,yData,zData,s=None):
initialFitReturn = SmoothBivariateSpline(xData,yData,zData,kx=2,ky=2,s=s).ev(xTest,yTest)
adjXData = np.append(xData,xTest)
adjYData = np.append(yData,yTest)
posAdjZData = np.append(zData,initialFitReturn + 1)
negAdjZData = np.append(zData,initialFitReturn - 1)
posFitReturn = SmoothBivariateSpline(adjXData,adjYData,posAdjZData,kx=2,ky=2,s=s).ev(xTest,yTest)
negFitReturn = SmoothBivariateSpline(adjXData,adjYData,negAdjZData,kx=2,ky=2,s=s).ev(xTest,yTest)
posGive = posFitReturn - initialFitReturn
negGive = initialFitReturn - negFitReturn
give = np.mean([posGive, negGive])
return give
def fit_give(xTest,yTest,xData,yData,zData,s=None):
dim = np.core.fromnumeric.shape(xTest)
if np.size(dim) == 0:
give = single_fit_give(xTest,yTest,xData,yData,zData,s=s)
return give
give = np.zeros(dim)
if np.size(dim) == 1:
for i in range(dim[0]):
give[i] = single_fit_give(xTest[i],yTest[i],xData,yData,zData,s=s)
return give
for i in range(dim[0]):
for j in range (dim[1]):
give[i,j] = single_fit_give(xTest[i,j],yTest[i,j],xData,yData,zData,s=s)
return give
In [12]:
try:
bivariateSpline = SmoothBivariateSpline(width[ref],eqPonA[ref],cutoutFactor[ref],kx=2,ky=2)
except:
splinefail = True
else:
splinefail = False
In [13]:
x = np.arange(2,14,0.1)
y = x
xx, yy = np.meshgrid(x,y)
if not(splinefail):
mesh_width = xx.copy()
mesh_width[xx>yy] = yy[xx>yy]
mesh_length = xx.copy()
mesh_length[yy>xx] = yy[yy>xx]
mesh_ratio = mesh_width / mesh_length
mesh_eqPonA = 2*( 3*(mesh_ratio+1) - np.sqrt( (3*mesh_ratio+1)*(mesh_ratio+3) ) ) / mesh_width
zz = bivariateSpline.ev(mesh_width, mesh_eqPonA)
give = fit_give(mesh_width,mesh_eqPonA,width[ref],eqPonA[ref],cutoutFactor[ref])
gap = angle_gap(mesh_width,mesh_eqPonA,width[ref],eqPonA[ref],1,1)
outOfTolerance = (give > 0.5) | (gap > 180)
zz[outOfTolerance] = np.nan
else:
zz = np.empty(np.shape(xx))
zz[:] = np.nan
zz = np.floor(zz*1e4)/1e4
In [14]:
if not(splinefail):
c = plt.contourf(xx,yy,zz,100)
plt.colorbar(c)
plt.contour(xx,yy,give,levels=[0.5],colors='g')
plt.contour(xx,yy,gap,levels=[180],colors='r')
plt.scatter(width[ref],length[ref],s=30)
plt.scatter(length[ref],width[ref],s=30,alpha=0.2)
plt.xlabel('Width (cm)')
plt.ylabel('Length (cm)')
plt.title('Bivariate spline fit for %d MeV and %d cm applicator' %(energyRequest, applicatorRequest))
plt.axis('equal')
In [15]:
dtick = 0.02
colourWashMin = int(np.floor(cutoutFactor.min()*100))/100
colourWashMax = int(np.ceil(cutoutFactor.max()*100))/100
nticks = int(np.floor((colourWashMax - colourWashMin)*100 + 1))
tick0 = colourWashMin
disp = 0
if not(splinefail):
measured1 = Scatter(x = width[ref],
y = length[ref],
text=marker_labels,
mode='markers',
marker=Marker(size=10),
name='Measured data',
xaxis='x1')
measured2 = Scatter(x = length[ref],
y = width[ref],
text=marker_labels,
mode='markers',
marker=Marker(size=10,opacity=0.3),
name='Measured data',
xaxis='x1')
colorbar = ColorBar(title='Cutout factor',
titleside='right',
autotick=False,
nticks=nticks,
dtick=dtick,
tick0=tick0)
contour1 = Contour(z=zz,
x=x,
y=y,
opacity=0.8,
colorscale='Jet',
ncontours='100',
zauto=False,
zmin=colourWashMin,
zmax=colourWashMax,
autocontour=True,
contours=Contours(coloring='heatmap'),
line=Line(width=0),
name='Spline fit',
xaxis='x1',
colorbar=colorbar)
contour2 = Contour(z=zz,
x=x,
y=y,
opacity=0.8,
colorscale='Jet',
ncontours='100',
zauto=False,
zmin=colourWashMin,
zmax=colourWashMax,
autocontour=True,
contours=Contours(coloring='heatmap'),
line=Line(width=0),
name='Spline fit',
xaxis='x2',
colorbar=colorbar)
data = Data([measured1,measured2,contour1,contour2])
layout = Layout(xaxis1=XAxis(domain=[0,0.45],title='Width (cm)'),
xaxis2=XAxis(domain=[0.55,1],title='Width (cm)'),
yaxis=YAxis(title='Length (cm)'),
title='Bivariate spline fit for %d MeV and %d cm applicator' %(energyRequest, applicatorRequest),
hovermode='closest',
showlegend=False)
fig_widthlength = Figure(data=data, layout=layout)
disp = py.iplot(fig_widthlength,
filename='%s_electron-fits_widthlength_%dMeV_%dcm' %(key,energyRequest, applicatorRequest))
disp
Out[15]:
In [16]:
if not(splinefail):
residuals_lengthwidth = cutoutFactor[ref] - bivariateSpline.ev(width[ref], eqPonA[ref])
binSize = 0.0015
binStart = np.floor(residuals_lengthwidth.min()/binSize)*binSize
binEnd = np.ceil(residuals_lengthwidth.max()/binSize)*binSize
bins = np.arange(binStart,binEnd+binSize,binSize)
dbins = bins[1] - bins[0]
binsTrans = bins - dbins/2
binsTrans = binsTrans.reshape(-1,1)
binNum = np.argmin(abs(binsTrans - residuals_lengthwidth),0)
representative_height = np.zeros(len(binNum))
for i in range(len(bins)):
binRef = (binNum == i)
representative_height[binRef] = np.arange(sum(binRef)) + 1
plt.hist(residuals_lengthwidth,bins,alpha=0.5)
plt.scatter(residuals_lengthwidth,representative_height,zorder=2,s=50,)
In [17]:
disp = 0
if not(splinefail):
residualsScatter = Scatter(x=residuals_lengthwidth,
y=representative_height,
text=marker_labels,
mode='markers',
marker=Marker(size=15),
name='Measured data')
residualsHistogram = Histogram(x=residuals_lengthwidth,
autobinx=False,
xbins=XBins(start=binStart,end=binEnd,size=binSize),
marker=Marker(color='rgb(165,231,255,0.5)',
line=Line(color='black',width=2)),
name='Residuals histogram')
# residualsHistogram = Histogram(x=residuals_lengthwidth)
data = Data([residualsScatter,residualsHistogram])
layout = Layout(title='Histogram of bivariate spline fit residuals',
xaxis=XAxis(title='Cutout factor residual error (measured factor - fitted factor)', zeroline=False),
yaxis=YAxis(title='Frequency', zeroline=False),
hovermode='closest',
showlegend=False)
fig_residualsHist = Figure(data=data, layout=layout)
disp = py.iplot(fig_residualsHist, filename='%s_electron-fits_widthlength_residual_%dMeV_%dcm' %(key,energyRequest, applicatorRequest))
disp
Out[17]:
In [18]:
if not(splinefail):
t, ttestProb = ttest_1samp(residuals_lengthwidth,0)
res_ttest = "%0.4f" %(ttestProb)
print(res_ttest)
W, shapiroProb = shapiro(residuals_lengthwidth)
res_norm = "%0.4f" %(shapiroProb)
print(res_norm)
res_std = "%0.4f" %(np.std(residuals_lengthwidth))
print(res_std)
res_mean = "%0.4f" %(np.mean(residuals_lengthwidth))
print(res_mean)
In [19]:
width_url = py.plot(fig_width,
filename='%s_electron-fits_widthalone_%dMeV_%dcm' %(key,energyRequest, applicatorRequest),
auto_open=False,)
if not(splinefail):
widthlength_url = py.plot(fig_widthlength,
filename='%s_electron-fits_widthlength_%dMeV_%dcm' %(key,energyRequest, applicatorRequest),
auto_open=False,)
residual_widthlength_url = py.plot(fig_residualsHist,
filename='%s_electron-fits_widthlength_residual_%dMeV_%dcm' %(key,energyRequest, applicatorRequest),
auto_open=False,)
In [20]:
html_string = '''
<html>
<head>
<link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.1/css/bootstrap.min.css">
<style>
body {
margin:60 130;
background:whitesmoke;
}
section.left {
width:65%;
float:left;
}
section {
padding:20px
}
aside {
padding:20px
}
td {
padding:20px;
}
</style>
</head>
<body>
<section>
<h1>Electron spline modelling report</h1>
<h2>Energy: ''' + str(energyRequest) + ''' MeV   
Applicator: ''' + str(applicatorRequest) + ''' cm   
<i>Include reference conditions here</i></h2>
</section>
<section class="left">
<h3>Width alone fits for lookup purposes</h3>
<p>Your measured data is given as blue dots. Hovering your mouse over these dots will provide
you with the relevant information for that datapoint.</p>
<p>The univariate fits provided, in the figure just below, are intented to be used
as a guide for lookup purposes.</p>
</section>
<section>
<iframe width="100%" height="400px" frameborder="0" seamless="seamless" scrolling="no" \
src="''' + width_url + '''.embed?width=100%&height=100%"></iframe>
</section>'''
if not(splinefail):
html_string += '''
<section class="left">
<h3>Bivariate spline fit for interpolation purposes</h3>
<p>Below is the bivariate spline fit formed using the parameters <i>width</i> and <i>equivalent perimeter on area</i>.
For ease of use the results are plotted with respect to <i>length</i> and <i>width</i>. The valid fitting region is given as a
colourwash. The colour of the colourwash is representative of the fitted cutout factor.</p>
<p>Hovering your mouse over the righthand plot will give the cutout factor results in the <i>z</i> parameter.
Where the <i>z</i> paramater returns null, and there is no colour wash, is not a valid fitting region. This
may be due either to not having enough points in the direct vicinity, or the region being outside the external data
boundary.</p>
<p>Hovering your mouse over the lefthand plot will provide you with a way to lookup individual measured data points.
Your measured data is given in blue. The equivalent shape with <i>width</i> and <i>length</i> inverted is given in orange.</p>
</section>
<aside><i>Diagram of length and width using an ellipse goes here...
maybe have a semi irregular cutout and an equivalent ellipse overlay to represent length and width</i></aside>
<section>
<iframe width="100%" height="600px" frameborder="0" seamless="seamless" scrolling="no" \
src="''' + widthlength_url + '''.embed?width=100%&height=100%"></iframe>
</section>
<section class="left">
<h4>Histogram of differences between measured cutout factors and fitted factors</h4>
<p>The blue data points are your measured data grouped according to which histogram bin
it was placed within. Just as before hovering your mouse over these data points will allow you
to lookup the relevant datapoint.</p>
</section>
<section>
<iframe width="100%" height="600px" frameborder="0" seamless="seamless" scrolling="no" \
src="''' + residual_widthlength_url + '''.embed?width=100%&height=100%"></iframe>
</section>
<section>
<h4>Statistics with respect to the cutout factor residual error given in the prior histogram</h4>
<table>
<tr>
<td>Mean:</td> <td>''' + res_mean + '''</td>
</tr>
<tr>
<td>Individual T-test for zero population mean p-value:</td> <td>''' + res_ttest + '''</td>
</tr>
<tr>
<td>Standard deviation:</td> <td>''' + res_std + '''</td>
</tr>
<tr>
<td>Shapiro normality p-value:</td> <td>''' + res_norm + '''</td>
</tr>
</table>
</section>
'''
else:
html_string += '''
<section>
<h3>Bivariate spline fit for interpolation purposes</h3>
<p>Bivariate spline failed. Likely not enough data.</p>
<section>
'''
html_string += '''
</body>
</html>'''
f = open(path2outputfile,'w')
f.write(html_string)
f.close()