Executed: Wed Jan 18 21:31:13 2017
Duration: 697 seconds.
In [1]:
measurement_id = 0 # possible values: 0, 1, 2
windows = (5, 30, 60)
In [2]:
# Cell inserted during automated execution.
windows = (5, 30, 60)
measurement_id = 0
In [3]:
import time
from pathlib import Path
import pandas as pd
from scipy.stats import linregress
from IPython.display import display
In [4]:
from fretbursts import *
In [5]:
sns = init_notebook(fs=14)
In [6]:
import lmfit; lmfit.__version__
Out[6]:
In [7]:
import phconvert; phconvert.__version__
Out[7]:
In [8]:
dir_ = 'data/multispot_'
filenames = [
dir_+'2015-07-31_bubble-bubble-run-off-kinetics-800mW-steer110_12.hdf5',
dir_+'2015-07-29_bubble-bubble-open-complex-run-off-kinetics-600mW-steer110_7.hdf5',
dir_+'2015-07-30_bubble-bubble-run-off-kinetics-800mW-steer110_8.hdf5']
start_times = [900, 600, 900] # time of NTP injection and start of kinetics
In [9]:
filename = filenames[measurement_id]
start_time = start_times[measurement_id]
In [10]:
filename
Out[10]:
In [11]:
import os
assert os.path.exists(filename)
Load and process the data:
In [12]:
d = loader.photon_hdf5(filename)
In [13]:
d.time_max
Out[13]:
Compute background and burst search:
In [14]:
d.calc_bg(bg.exp_fit, time_s=10, tail_min_us='auto', F_bg=1.7)
Perform a background plot as a function of the channel:
In [15]:
mch_plot_bg(d)
Let's take a look at the photon waiting times histograms and at the fitted background rates:
In [16]:
dplot(d, hist_bg);
Using dplot exactly in the same way as for the single-spot data has now generated 8 subplots, one for each channel.
Let's plot a timetrace for the background to see is there are significat variations during the measurement:
In [17]:
dplot(d, timetrace_bg);
xlim(start_time - 150, start_time + 150)
Out[17]:
We can look at the timetrace of the photon stream (binning):
In [18]:
#dplot(d, timetrace)
#xlim(2, 3); ylim(-100, 100);
In [19]:
d.burst_search(m=10, F=5)
In [20]:
ds = d.select_bursts(select_bursts.size, th1=30)
In [21]:
ds0 = ds.select_bursts(select_bursts.time, time_s1=0, time_s2=start_time-10)
In [22]:
dplot(ds0, hist_fret, pdf=False);
In [23]:
dm0 = ds0.collapse()
In [24]:
dplot(dm0, hist_fret, pdf=False);
In [25]:
weights = 'size'
bext.bursts_fitter(dm0, weights=weights)
dm0.E_fitter.fit_histogram(mfit.factory_three_gaussians(p1_center=0.05, p2_center=0.6, p3_center=0.9), verbose=False)
dplot(dm0, hist_fret, show_model=True, weights=weights);
dm0.E_fitter.params
Out[25]:
In [26]:
weights = None
bext.bursts_fitter(dm0, weights=weights)
dm0.E_fitter.fit_histogram(mfit.factory_three_gaussians(p1_center=0.05, p2_center=0.6, p3_center=0.9), verbose=False)
dplot(dm0, hist_fret, show_model=True, weights=weights);
dm0.E_fitter.params
Out[26]:
In [27]:
def gauss3(**params0):
peak1 = lmfit.models.GaussianModel(prefix='p1_')
peak3 = lmfit.models.GaussianModel(prefix='p3_')
peak2 = lmfit.models.GaussianModel(prefix='p2_')
model = peak1 + peak2 + peak3
model.set_param_hint('p1_center', **{'value': 0.0, 'min': 0.0, 'max': 0.2, **params0.get('p1_center', {})})
model.set_param_hint('p2_center', **{'value': 0.5, 'min': 0.0, 'max': 1.0, **params0.get('p2_center', {})})
model.set_param_hint('p3_center', **{'value': 0.9, 'min': 0.8, 'max': 1.0, **params0.get('p3_center', {})})
for sigma in ['p%d_sigma' % i for i in (1, 2, 3)]:
model.set_param_hint(sigma, **{'value': 0.02, 'min': 0.01, **params0.get(sigma, {})})
for ampl in ['p%d_amplitude' % i for i in (1, 2, 3)]:
model.set_param_hint(ampl, **{'value': 0.333, 'min': 0.01, **params0.get(ampl, {})})
model.name = '3 gauss peaks'
return model
In [28]:
#%matplotlib notebook
In [29]:
#fig, ax = plt.subplots(figsize=(12, 8))
#dplot(dm0, scatter_fret_size, ax=ax)
In [30]:
bext.bursts_fitter(dm0, weights=None)
dm0.E_fitter.fit_histogram(gauss3(), verbose=False)
mfit.plot_mfit(dm0.E_fitter)
params_3gauss = dm0.E_fitter.params
plt.xlabel('E')
plt.ylabel('PDF')
plt.title('')
#dir_ = r'C:\Data\Antonio\docs\conferences\Seaborg2015\figures/'
#plt.savefig(dir_+'Realtime kinetics FRET hist', dpi=200, bbox_inches='tight')
params_3gauss
Out[30]:
In [31]:
dsc = ds.collapse()
In [32]:
dm_final = dsc.select_bursts(select_bursts.time, time_s1=start_time+300, time_s2=ds.time_max + 1)
dm_final.num_bursts
Out[32]:
In [33]:
dm_final1 = dsc.select_bursts(select_bursts.time, time_s1=start_time+100, time_s2=start_time+1600)
dm_final1.num_bursts
Out[33]:
In [34]:
dm_final2 = dsc.select_bursts(select_bursts.time, time_s1=start_time + 2100, time_s2=ds.time_max + 1)
dm_final2.num_bursts
Out[34]:
In [35]:
bext.bursts_fitter(dm_final1, weights=None)
model = gauss3()
model.set_param_hint('p2_center', value=params_3gauss.p2_center[0], vary=False)
dm_final1.E_fitter.fit_histogram(model, verbose=False)
fig, ax = plt.subplots(figsize=(12, 6))
mfit.plot_mfit(dm_final1.E_fitter, ax=ax)
params_3gauss1 = dm_final1.E_fitter.params
params_3gauss1
Out[35]:
In [36]:
bext.bursts_fitter(dm_final2, weights=None)
model = gauss3()
model.set_param_hint('p2_center', value=params_3gauss.p2_center[0], vary=False)
dm_final2.E_fitter.fit_histogram(model, verbose=False)
fig, ax = plt.subplots(figsize=(12, 6))
mfit.plot_mfit(dm_final2.E_fitter, ax=ax)
params_3gauss1 = dm_final2.E_fitter.params
params_3gauss1
Out[36]:
In [37]:
bext.bursts_fitter(dm_final, weights=None)
model = gauss3()
model.set_param_hint('p2_center', value=params_3gauss.p2_center[0], vary=False)
dm_final.E_fitter.fit_histogram(model, verbose=False)
fig, ax = plt.subplots(figsize=(12, 6))
mfit.plot_mfit(dm_final.E_fitter, ax=ax)
params_3gauss1 = dm_final.E_fitter.params
params_3gauss1
Out[37]:
In [38]:
#del params_3gauss0
In [39]:
if 'params_3gauss0' not in locals():
params_3gauss0 = params_3gauss.copy()
params_3gauss0.p3_center = params_3gauss1.p3_center
In [40]:
params_3gauss0.p1_amplitude + params_3gauss0.p2_amplitude + params_3gauss0.p3_amplitude
Out[40]:
In [41]:
'params_3gauss0' in locals()
Out[41]:
In [42]:
from scipy import optimize
In [43]:
params_fixed = dict(
mu1=float(params_3gauss0.p1_center),
mu2=float(params_3gauss0.p2_center),
mu3=float(params_3gauss0.p3_center),
sig1=float(params_3gauss0.p1_sigma),
sig2=float(params_3gauss0.p2_sigma),
sig3=float(params_3gauss0.p3_sigma),
)
In [44]:
def em_weights_3gauss(x, a2, a3, mu1, mu2, mu3, sig1, sig2, sig3):
"""Responsibility function for a 3-Gaussian model.
Returns 3 arrays of size = x.size: the responsibility of
each Gaussian population.
"""
a1 = 1 - a2 - a3
assert np.abs(a1 + a2 + a3 - 1) < 1e-3
f1 = a1 * gauss_pdf(x, mu1, sig1)
f2 = a2 * gauss_pdf(x, mu2, sig2)
f3 = a3 * gauss_pdf(x, mu3, sig3)
γ1 = f1 / (f1 + f2 + f3)
γ2 = f2 / (f1 + f2 + f3)
γ3 = f3 / (f1 + f2 + f3)
return γ1, γ2, γ3
def em_fit_3gauss(x, a2_0, a3_0, params_fixed, print_every=10, max_iter=100, rtol=1e-3):
"""Fit amplitude of 3_Gaussian model using Expectation-Maximization.
Only 2 amplitudes are fitted (a2, a3), the first peak is derived imposing
that the PDF sums to 1.
"""
a2_new, a3_new = a2_0, a3_0
rel_change = 1
i = 0
while rel_change > rtol and i < max_iter:
# E-step
γ1, γ2, γ3 = em_weights_3gauss(x, a2_new, a3_new, **params_fixed)
assert np.allclose(γ1.sum() + γ2.sum() + γ3.sum(), x.size)
# M-step
a2_old, a3_old = a2_new, a3_new
a2_new = γ2.sum()/γ2.size
a3_new = γ3.sum()/γ3.size
# Convergence
rel_change = (np.abs((a2_old - a2_new)/a2_new)
+ np.abs((a3_old - a3_new)/a3_new))
i += 1
if (i % print_every) == 0:
print(i, a2_new, a3_new, rel_change)
return a2_new, a3_new, i
In [45]:
from matplotlib.pylab import normpdf as gauss_pdf
# Model PDF to be maximized
def model_pdf(x, a2, a3, mu1, mu2, mu3, sig1, sig2, sig3):
a1 = 1 - a2 - a3
#assert np.abs(a1 + a2 + a3 - 1) < 1e-3
return (a1 * gauss_pdf(x, mu1, sig1) +
a2 * gauss_pdf(x, mu2, sig2) +
a3 * gauss_pdf(x, mu3, sig3))
# Function to be minimized by lmfit
def func2min_lmfit(params, x):
a2 = params['a2'].value
a3 = params['a3'].value
mu1 = params['mu1'].value
mu2 = params['mu2'].value
mu3 = params['mu3'].value
sig1 = params['sig1'].value
sig2 = params['sig2'].value
sig3 = params['sig3'].value
return -np.sqrt(np.log(model_pdf(x, a2, a3, mu1, mu2, mu3, sig1, sig2, sig3)))
# Function to be minimized by scipy
def func2min_scipy(params_fit, params_fixed, x):
a2, a3 = params_fit
mu1 = params_fixed['mu1']
mu2 = params_fixed['mu2']
mu3 = params_fixed['mu3']
sig1 = params_fixed['sig1']
sig2 = params_fixed['sig2']
sig3 = params_fixed['sig3']
return -np.log(model_pdf(x, a2, a3, mu1, mu2, mu3, sig1, sig2, sig3)).sum()
# create a set of Parameters
params = lmfit.Parameters()
params.add('a2', value=0.33, min=0)
params.add('a3', value=0.33, min=0)
for k, v in params_fixed.items():
params.add(k, value=v, vary=False)
In [46]:
x = dm0.E_
x
#result = lmfit.minimize(func2min_lmfit, params, args=(x,), method='nelder')
#lmfit.report_fit(result.params)
Out[46]:
In [47]:
#optimize.brute(func2min_scipy, ranges=((0.01, 0.99), (0.01, 0.99)), Ns=101, args=(params, x))
In [48]:
res = optimize.minimize(func2min_scipy, x0=[0.33, 0.33], args=(params_fixed, x), method='Nelder-Mead')
res
Out[48]:
In [49]:
res = optimize.minimize(func2min_scipy, x0=[0.33, 0.33], args=(params_fixed, x), bounds=((0,1), (0,1)), method='SLSQP')
res
Out[49]:
In [50]:
res = optimize.minimize(func2min_scipy, x0=[0.33, 0.33], args=(params_fixed, x), bounds=((0,1), (0,1)), method='TNC')
res
Out[50]:
In [51]:
bins = np.arange(-0.1, 1.1, 0.025)
plt.hist(x, bins, histtype='step', lw=2, normed=True);
xx = np.arange(-0.1, 1.1, 0.005)
#plt.plot(xx, model_pdf(xx, params))
plt.plot(xx, model_pdf(xx, a2=res.x[0], a3=res.x[1], **params_fixed))
Out[51]:
In [52]:
def _kinetics_fit_em(dx, a2_0, a3_0, params_fixed, **kwargs):
kwargs = {'max_iter': 200, 'print_every': 201, **kwargs}
a2, a3, i = em_fit_3gauss(dx.E_, a2_0, a3_0, params_fixed, **kwargs)
return a2, a3, i < kwargs['max_iter']
def _kinetics_fit_ll(dx, a2_0, a3_0, params_fixed, **kwargs):
kwargs = {'method':'Nelder-Mead', **kwargs}
res = optimize.minimize(func2min_scipy, x0=[a2_0, a3_0], args=(params_fixed, dx.E_),
**kwargs)
return res.x[0], res.x[1], res.success
def _kinetics_fit_hist(dx, a2_0, a3_0, params_fixed):
E_fitter = bext.bursts_fitter(dx)
model = gauss3(p1_amplitude={'value': 1 - a2_0 - a3_0},
p2_amplitude={'value': a2_0},
p3_amplitude={'value': a3_0})
model.set_param_hint('p1_center', value=params_fixed['mu1'], vary=False)
model.set_param_hint('p2_center', value=params_fixed['mu2'], vary=False)
model.set_param_hint('p3_center', value=params_fixed['mu3'], vary=False)
model.set_param_hint('p1_sigma', value=params_fixed['sig1'], vary=False)
model.set_param_hint('p2_sigma', value=params_fixed['sig2'], vary=False)
model.set_param_hint('p3_sigma', value=params_fixed['sig3'], vary=False)
E_fitter.fit_histogram(model, verbose=False)
return (float(E_fitter.params.p2_amplitude),
float(E_fitter.params.p3_amplitude),
dx.E_fitter.fit_res[0].success)
def kinetics_fit(ds_slices, params_fixed, a2_0=0.33, a3_0=0.33, method='em', **method_kws):
fit_func = {
'em': _kinetics_fit_em,
'll': _kinetics_fit_ll,
'hist': _kinetics_fit_hist}
fit_list = []
for dx in ds_slices:
a2, a3, success = fit_func[method](dx, a2_0, a3_0, params_fixed, **method_kws)
df_i = pd.DataFrame(data=dict(p2_amplitude=a2, p3_amplitude=a3,
p1_center=params_fixed['mu1'], p2_center=params_fixed['mu2'],
p3_center=params_fixed['mu3'], p1_sigma=params_fixed['sig1'],
p2_sigma=params_fixed['sig2'], p3_sigma=params_fixed['sig3'],
tstart=dx.slice_tstart, tstop=dx.slice_tstop,
tmean=0.5*(dx.slice_tstart + dx.slice_tstop)),
index=[0.5*(dx.slice_tstart + dx.slice_tstop)])
if not success:
print('* ', end='', flush=True)
continue
fit_list.append(df_i)
return pd.concat(fit_list)
In [53]:
start_time/60
Out[53]:
In [54]:
def print_slices(moving_window_params):
msg = ' - Slicing measurement:'
for name in ('start', 'stop', 'step', 'window'):
msg += ' %s = %.1fs' % (name, moving_window_params[name])
print(msg, flush=True)
num_slices = len(bext.moving_window_startstop(**moving_window_params))
print(' Number of slices %d' % num_slices, flush=True)
In [55]:
t1 = time.time()
time.ctime()
Out[55]:
In [56]:
dsc = ds.collapse()
dsc.calc_max_rate(m=10)
dsc_high = dsc.select_bursts(select_bursts.E, E1=0.88)
In [57]:
step = 1
params = {}
for window in windows:
moving_window_params = dict(start=0, stop=dsc.time_max, step=step, window=window)
print_slices(moving_window_params)
ds_slices = bext.moving_window_chunks(dsc, time_zero=start_time, **moving_window_params)
for meth in ['em', 'll', 'hist']:
print(' >>> Fitting method %s ' % meth, end='', flush=True)
p = kinetics_fit(ds_slices, params_fixed, method=meth)
print(flush=True)
p['kinetics'] = p.p3_amplitude / (p.p2_amplitude + p.p3_amplitude)
p = p.round(dict(p1_center=3, p1_sigma=4, p2_amplitude=4, p2_center=3, p2_sigma=4, kinetics=4,
p3_amplitude=4, p3_center=3, p3_sigma=4))
params[meth, window, step] = p
In [58]:
print('Moving-window processing duration: %d seconds.' % (time.time() - t1))
In [59]:
moving_window_params['window'] = 30
moving_window_params
Out[59]:
In [60]:
ds_slices = bext.moving_window_chunks(dsc, **moving_window_params)
ds_slices_high = bext.moving_window_chunks(dsc_high, **moving_window_params)
df = bext.moving_window_dataframe(**moving_window_params) - start_time
df['size_mean'] = [di.nt_.mean() for di in ds_slices]
df['size_max'] = [di.nt_.max() for di in ds_slices]
df['num_bursts'] = [di.num_bursts[0] for di in ds_slices]
df['burst_width'] = [di.mburst_.width.mean()*di.clk_p*1e3 for di in ds_slices]
df['burst_width_high'] = [di.mburst_.width.mean()*di.clk_p*1e3 for di in ds_slices_high]
df['phrate_mean'] = [di.max_rate_.mean() for di in ds_slices]
In [61]:
df = df.round(dict(tmean=1, tstart=1, tstop=1, size_mean=2, size_max=1,
burst_width=2, burst_width_high=2, phrate_mean=1))
df
Out[61]:
In [62]:
labels = ('num_bursts', 'burst_width', 'phrate_mean')
fig, axes = plt.subplots(len(labels), 1, figsize=(12, 3*len(labels)))
for ax, label in zip(axes, labels):
ax.plot(label, data=df)
ax.legend(loc='best')
#ax.set_ylim(0)
In [63]:
# %%timeit -n1 -r1
# meth = 'em'
# print(' >>> Fitting method %s' % meth, flush=True)
# p = kinetics_fit(ds_slices, params_fixed, method=meth)
In [64]:
# %%timeit -n1 -r1
# meth = 'hist'
# print(' >>> Fitting method %s' % meth, flush=True)
# p = kinetics_fit(ds_slices, params_fixed, method=meth)
In [65]:
# %%timeit -n1 -r1
# meth = 'll'
# print(' >>> Fitting method %s' % meth, flush=True)
# p = kinetics_fit(ds_slices, params_fixed, method=meth)
In [66]:
out_fname = 'results/%s_burst_data_vs_time__window%ds_step%ds.csv' % (
Path(filename).stem, moving_window_params['window'], moving_window_params['step'])
out_fname
Out[66]:
In [67]:
df.to_csv(out_fname)
In [68]:
# np.abs((params['em', 30, 1] - params['ll', 30, 1]).p2_amplitude).max()
In [69]:
methods = ('em', 'll', 'hist')
In [70]:
for meth in methods:
plt.figure(figsize=(14, 3))
plt.plot(params[meth, 5, 1].index, params[meth, 5, 1].kinetics, 'h', color='gray', alpha=0.2)
plt.plot(params[meth, 30, 1].index, params[meth, 30, 1].kinetics, 'h', alpha=0.3)
In [71]:
# (params['em', 5, 1].kinetics - params['ll', 5, 1].kinetics).plot()
In [72]:
step = 1
for window in windows:
for meth in methods:
out_fname = ('results/' + Path(filename).stem +
'_%sfit_ampl_only__window%ds_step%ds.csv' % (meth, window, step))
print('- Saving: ', out_fname)
params[meth, window, step].to_csv(out_fname)
In [73]:
d
Out[73]:
In [74]: