I have a data output file from a LabVIEW program (that I cannot modify), and I need to plot the data and fit to a theoretical curve. I accomplished this with Python (with which I am entirely self-taught).
My code runs and produces the desired output just fine, but I feel that it is much too messy, unreadable, and hacked together, not to mention that it was time-consuming to produce (especially for a task as common as reading and plotting data). What code smells are present, and is there a faster/preferred way to do this?
What needs to be plotted?
From the data file (below), there are 4 distinct values of HWP B
: 0, 45, 90, and 135. With this angle held constant, HWP A
(which varies from 0 to 180), must be plotted on the x-axis against a function of the columns A
, B
, and AB
(namely, AB - 2*A*B*dt
, where dt
is a constant). The following plot is produced, as well as 4 lines which are used later by me.
(For those wondering, this data follows from an experiment outlined by Dehlinger and Mitchell in a paper on entanglement available on the arXiv.)
Code:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
''' Read, analyze and plot data from EPR state quantification data file '''
# time interval for data collection
dt = 5e-9
# read in data from file
fname_epr_quant = "qie-epr-quant.dat"
with open(fname_epr_quant,'r') as f:
eq_dat = np.array([ line.split(', ') for line in list(filter(None,f.read().split('\n'))) ])
# create dictionary for the data so that it can be referenced by name
varnames = eq_dat[0]
eq_n2i = { name:index for index,name in enumerate(varnames)} #name-to-index
eq_dict = {}
for l_line in eq_dat[1:]:
seriesname = l_line[ eq_n2i['Comments:'] ]
if seriesname not in eq_dict.keys():
eq_dict[seriesname] = { varname:np.array([]) for varname in varnames }
for l_val,l_varname in zip(l_line,varnames):
try:
value = float(l_val)
except:
value = l_val
eq_dict[seriesname][l_varname] = np.append(eq_dict[seriesname][l_varname], value)
# compute 2 new quantities from the data ('acc' and 'AB-acc')
# given that accidentals = 2*A*B*dt
for l_seriesname,l_series in eq_dict.items():
accidentals = 2*l_series['A']*l_series['B']*dt
eq_dict[l_seriesname]['acc'] = accidentals
eq_dict[l_seriesname]['AB-acc'] = l_series['AB'] - accidentals
# theoretical N for plotting and fitting data points
def N_theo(a,b,A,C,th,phi):
'''Accepts a,b,th,phi in DEGREES'''
a,b,th,phi = np.pi/180 * np.array([a,b,th,phi])
return C + A*( (np.sin(a)*np.sin(b)*np.cos(th))**2 +
(np.cos(a)*np.cos(b)*np.sin(th))**2 +
1/4*np.sin(2*a)*np.sin(2*b)*np.sin(2*th)*np.cos(phi) )
# range of values for the x-axis
deg = np.linspace(0,180)
# plotting
plt.figure(fname_epr_quant)
plt.cla()
for (l_seriesname,l_vars),color in zip(eq_dict.items(),'rgbk'):
x = l_vars['HWP A']
y = l_vars['AB-acc']
b_plt = l_vars['HWP B'][0]
N_plt = lambda a,A,C,th,phi: N_theo(a,b_plt,A,C,th,phi)
popt, pcov = curve_fit( N_plt, x, y, bounds=([100,1,30,15],[1000,100,60,35]) )
print(b_plt,popt) # print fitting parameters for later analysis
plt.plot( x, y, '.--'+color, label=l_seriesname+' data' )
plt.plot( deg, N_theo(deg, b_plt, *popt), '-'+color, label=l_seriesname+' fit' )
plt.legend()
plt.show(block=False)
(Note: I use l_
for loop variables due to my relatively recent discovery that Python doesn't keep a namespace for loops, so I may be overwriting important variables.)
Data file:
(Note: I have delimited the data with ,<space>
, since \t
doesn't survive the copy/paste into a text editor. :%s/, /\t/g
will trivially recover tabs, if you wish.)
A, B, A', B', A'B, AB, AB', A'B', HWP A, HWP B, E, Comments:, Duration (s), Sub. Acc.?
19346.20, 23218.60, 0.00, 0.00, 0.00, 442.70, 0.00, 0.00, 0.00, 0.00, 1.00, B=0, 10.0, N
19356.10, 23186.20, 0.00, 0.00, 0.00, 426.50, 0.00, 0.00, 8.00, 0.00, 1.00, B=0, 10.0, N
19348.30, 23126.20, 0.00, 0.00, 0.00, 401.90, 0.00, 0.00, 16.00, 0.00, 1.00, B=0, 10.0, N
19075.30, 23145.20, 0.00, 0.00, 0.00, 349.60, 0.00, 0.00, 24.00, 0.00, 1.00, B=0, 10.0, N
18851.40, 23115.40, 0.00, 0.00, 0.00, 299.80, 0.00, 0.00, 32.00, 0.00, 1.00, B=0, 10.0, N
18791.60, 23106.60, 0.00, 0.00, 0.00, 259.00, 0.00, 0.00, 40.00, 0.00, 1.00, B=0, 10.0, N
18561.20, 23129.80, 0.00, 0.00, 0.00, 193.90, 0.00, 0.00, 48.00, 0.00, 1.00, B=0, 10.0, N
18335.30, 23102.50, 0.00, 0.00, 0.00, 136.50, 0.00, 0.00, 56.00, 0.00, 1.00, B=0, 10.0, N
18107.90, 23129.20, 0.00, 0.00, 0.00, 92.70, 0.00, 0.00, 64.00, 0.00, 1.00, B=0, 10.0, N
18042.30, 23135.50, 0.00, 0.00, 0.00, 46.70, 0.00, 0.00, 72.00, 0.00, 1.00, B=0, 10.0, N
17909.20, 23042.20, 0.00, 0.00, 0.00, 22.60, 0.00, 0.00, 80.00, 0.00, 1.00, B=0, 10.0, N
17920.30, 23020.20, 0.00, 0.00, 0.00, 13.40, 0.00, 0.00, 88.00, 0.00, 1.00, B=0, 10.0, N
17945.50, 23048.60, 0.00, 0.00, 0.00, 18.30, 0.00, 0.00, 96.00, 0.00, 1.00, B=0, 10.0, N
17944.50, 23076.20, 0.00, 0.00, 0.00, 44.50, 0.00, 0.00, 104.00, 0.00, 1.00, B=0, 10.0, N
18032.40, 23111.60, 0.00, 0.00, 0.00, 80.50, 0.00, 0.00, 112.00, 0.00, 1.00, B=0, 10.0, N
18254.20, 23131.00, 0.00, 0.00, 0.00, 128.90, 0.00, 0.00, 120.00, 0.00, 1.00, B=0, 10.0, N
18370.90, 23088.10, 0.00, 0.00, 0.00, 174.80, 0.00, 0.00, 128.00, 0.00, 1.00, B=0, 10.0, N
18598.20, 23112.80, 0.00, 0.00, 0.00, 238.80, 0.00, 0.00, 136.00, 0.00, 1.00, B=0, 10.0, N
18780.40, 23076.10, 0.00, 0.00, 0.00, 293.10, 0.00, 0.00, 144.00, 0.00, 1.00, B=0, 10.0, N
18965.40, 23039.90, 0.00, 0.00, 0.00, 338.60, 0.00, 0.00, 152.00, 0.00, 1.00, B=0, 10.0, N
19112.00, 23097.40, 0.00, 0.00, 0.00, 387.60, 0.00, 0.00, 160.00, 0.00, 1.00, B=0, 10.0, N
19185.00, 22957.50, 0.00, 0.00, 0.00, 417.90, 0.00, 0.00, 168.00, 0.00, 1.00, B=0, 10.0, N
19230.40, 23019.70, 0.00, 0.00, 0.00, 423.50, 0.00, 0.00, 176.00, 0.00, 1.00, B=0, 10.0, N
19255.50, 23201.20, 0.00, 0.00, 0.00, 428.50, 0.00, 0.00, 180.00, 0.00, 1.00, B=0, 10.0, N
18910.00, 22827.90, 0.00, 0.00, 0.00, 248.50, 0.00, 0.00, 0.00, 45.00, 1.00, B=45, 10.0, N
18828.00, 22829.80, 0.00, 0.00, 0.00, 321.40, 0.00, 0.00, 12.00, 45.00, 1.00, B=45, 10.0, N
18554.50, 22812.60, 0.00, 0.00, 0.00, 357.30, 0.00, 0.00, 24.00, 45.00, 1.00, B=45, 10.0, N
19159.90, 24043.50, 0.00, 0.00, 0.00, 386.70, 0.00, 0.00, 36.00, 45.00, 1.00, B=45, 10.0, N
18439.10, 23410.90, 0.00, 0.00, 0.00, 382.70, 0.00, 0.00, 48.00, 45.00, 1.00, B=45, 10.0, N
17666.20, 22784.00, 0.00, 0.00, 0.00, 342.20, 0.00, 0.00, 60.00, 45.00, 1.00, B=45, 10.0, N
17479.40, 22887.70, 0.00, 0.00, 0.00, 276.30, 0.00, 0.00, 72.00, 45.00, 1.00, B=45, 10.0, N
17353.90, 22847.90, 0.00, 0.00, 0.00, 213.30, 0.00, 0.00, 84.00, 45.00, 1.00, B=45, 10.0, N
17328.40, 22803.70, 0.00, 0.00, 0.00, 132.30, 0.00, 0.00, 96.00, 45.00, 1.00, B=45, 10.0, N
17613.40, 22917.50, 0.00, 0.00, 0.00, 74.30, 0.00, 0.00, 108.00, 45.00, 1.00, B=45, 10.0, N
17809.60, 22779.80, 0.00, 0.00, 0.00, 41.10, 0.00, 0.00, 120.00, 45.00, 1.00, B=45, 10.0, N
18226.40, 22855.30, 0.00, 0.00, 0.00, 35.10, 0.00, 0.00, 132.00, 45.00, 1.00, B=45, 10.0, N
18498.50, 22838.80, 0.00, 0.00, 0.00, 59.00, 0.00, 0.00, 144.00, 45.00, 1.00, B=45, 10.0, N
18939.00, 22883.70, 0.00, 0.00, 0.00, 107.50, 0.00, 0.00, 156.00, 45.00, 1.00, B=45, 10.0, N
19016.30, 22821.90, 0.00, 0.00, 0.00, 175.40, 0.00, 0.00, 168.00, 45.00, 1.00, B=45, 10.0, N
19077.90, 22873.60, 0.00, 0.00, 0.00, 255.10, 0.00, 0.00, 180.00, 45.00, 1.00, B=45, 10.0, N
19041.60, 21765.30, 0.00, 0.00, 0.00, 20.90, 0.00, 0.00, 0.00, 90.00, 1.00, B=90, 10.0, N
19085.70, 21803.40, 0.00, 0.00, 0.00, 43.20, 0.00, 0.00, 12.00, 90.00, 1.00, B=90, 10.0, N
18822.70, 21839.40, 0.00, 0.00, 0.00, 97.10, 0.00, 0.00, 24.00, 90.00, 1.00, B=90, 10.0, N
18511.70, 21795.40, 0.00, 0.00, 0.00, 163.40, 0.00, 0.00, 36.00, 90.00, 1.00, B=90, 10.0, N
18051.10, 21839.20, 0.00, 0.00, 0.00, 250.30, 0.00, 0.00, 48.00, 90.00, 1.00, B=90, 10.0, N
17727.90, 21838.70, 0.00, 0.00, 0.00, 312.30, 0.00, 0.00, 60.00, 90.00, 1.00, B=90, 10.0, N
17489.30, 21813.20, 0.00, 0.00, 0.00, 360.50, 0.00, 0.00, 72.00, 90.00, 1.00, B=90, 10.0, N
17371.70, 21847.80, 0.00, 0.00, 0.00, 383.50, 0.00, 0.00, 84.00, 90.00, 1.00, B=90, 10.0, N
17249.30, 21918.90, 0.00, 0.00, 0.00, 361.90, 0.00, 0.00, 96.00, 90.00, 1.00, B=90, 10.0, N
17475.90, 21818.20, 0.00, 0.00, 0.00, 314.60, 0.00, 0.00, 108.00, 90.00, 1.00, B=90, 10.0, N
17643.20, 21747.90, 0.00, 0.00, 0.00, 249.60, 0.00, 0.00, 120.00, 90.00, 1.00, B=90, 10.0, N
18132.30, 21869.30, 0.00, 0.00, 0.00, 181.90, 0.00, 0.00, 132.00, 90.00, 1.00, B=90, 10.0, N
18301.00, 21805.40, 0.00, 0.00, 0.00, 111.80, 0.00, 0.00, 144.00, 90.00, 1.00, B=90, 10.0, N
18835.90, 22196.70, 0.00, 0.00, 0.00, 46.00, 0.00, 0.00, 156.00, 90.00, 1.00, B=90, 10.0, N
18534.50, 21750.50, 0.00, 0.00, 0.00, 14.30, 0.00, 0.00, 168.00, 90.00, 1.00, B=90, 10.0, N
18568.90, 21831.90, 0.00, 0.00, 0.00, 14.80, 0.00, 0.00, 180.00, 90.00, 1.00, B=90, 10.0, N
18544.40, 22340.60, 0.00, 0.00, 0.00, 173.10, 0.00, 0.00, 0.00, 135.00, 1.00, B=135, 10.0, N
18485.30, 22349.00, 0.00, 0.00, 0.00, 105.70, 0.00, 0.00, 12.00, 135.00, 1.00, B=135, 10.0, N
18310.80, 22352.90, 0.00, 0.00, 0.00, 62.10, 0.00, 0.00, 24.00, 135.00, 1.00, B=135, 10.0, N
17917.60, 22278.50, 0.00, 0.00, 0.00, 35.50, 0.00, 0.00, 36.00, 135.00, 1.00, B=135, 10.0, N
17576.00, 22315.20, 0.00, 0.00, 0.00, 38.70, 0.00, 0.00, 48.00, 135.00, 1.00, B=135, 10.0, N
17259.60, 22274.40, 0.00, 0.00, 0.00, 62.10, 0.00, 0.00, 60.00, 135.00, 1.00, B=135, 10.0, N
16985.30, 22214.70, 0.00, 0.00, 0.00, 112.50, 0.00, 0.00, 72.00, 135.00, 1.00, B=135, 10.0, N
17070.90, 22526.70, 0.00, 0.00, 0.00, 179.80, 0.00, 0.00, 84.00, 135.00, 1.00, B=135, 10.0, N
17139.60, 22546.70, 0.00, 0.00, 0.00, 238.10, 0.00, 0.00, 96.00, 135.00, 1.00, B=135, 10.0, N
17320.80, 22484.00, 0.00, 0.00, 0.00, 299.10, 0.00, 0.00, 108.00, 135.00, 1.00, B=135, 10.0, N
17398.80, 22275.50, 0.00, 0.00, 0.00, 330.90, 0.00, 0.00, 120.00, 135.00, 1.00, B=135, 10.0, N
17798.80, 22231.60, 0.00, 0.00, 0.00, 342.80, 0.00, 0.00, 132.00, 135.00, 1.00, B=135, 10.0, N
18224.00, 22234.00, 0.00, 0.00, 0.00, 339.50, 0.00, 0.00, 144.00, 135.00, 1.00, B=135, 10.0, N
18391.90, 22246.30, 0.00, 0.00, 0.00, 292.10, 0.00, 0.00, 156.00, 135.00, 1.00, B=135, 10.0, N
18690.00, 22221.70, 0.00, 0.00, 0.00, 233.00, 0.00, 0.00, 168.00, 135.00, 1.00, B=135, 10.0, N
18667.20, 22349.00, 0.00, 0.00, 0.00, 169.00, 0.00, 0.00, 180.00, 135.00, 1.00, B=135, 10.0, N
TypeError: leastsq() got an unexpected keyword argument 'bounds'
fromcurve_fit()
\$\endgroup\$curve_fit
method takes abounds
parameter since scipy 0.17 according to the docs. \$\endgroup\$