In [ ]:
from fretbursts import *
sns = init_notebook()
In [ ]:
import seaborn
In [ ]:
import os
import pandas as pd
from IPython.display import display
from IPython.display import display, Math
In [ ]:
data_dir = './data/multispot/'
In [ ]:
data_dir = os.path.abspath(data_dir) + '/'
assert os.path.exists(data_dir), "Path '%s' does not exist." % data_dir
In [ ]:
from glob import glob
file_list = sorted(glob(data_dir + '*.hdf5'))
In [ ]:
labels = ['7d', '12d', '17d', '22d', '27d', 'DO']
files_dict = {lab: fname for lab, fname in zip(sorted(labels), file_list)}
files_dict
In [ ]:
PLOT_DIR = './figure/'
In [ ]:
import matplotlib as mpl
from cycler import cycler
bmap = sns.color_palette("Set1", 9)
colors = np.array(bmap)[(1,0,2,3,4,8,6,7), :]
colors_labels = ['blue', 'red', 'green', 'violet', 'orange', 'gray', 'brown', 'pink', ]
for c, cl in zip(colors, colors_labels):
locals()[cl] = tuple(c) # assign variables with color names
mpl.rcParams['axes.prop_cycle'] = cycler('color', colors)
sns.palplot(colors)
In [ ]:
## Background fit
bg_kwargs_auto = dict(fun=bg.exp_fit,
time_s = 30,
tail_min_us = 'auto',
F_bg=1.7,
)
## Burst search
F = 6
ph_sel = Ph_sel(Dex='Dem')
#ph_sel = Ph_sel('all')
size_min = 80
## D-only peak fit with KDE
bandwidth = 0.03
binwidth = 0.025
E_range_do = (-0.05, 0.1)
weights = 'size'
In [ ]:
df = pd.DataFrame(index=['7d', '12d', '17d', 'DO'], columns=range(8), dtype=float)
df.index.name = 'Sample'
df.columns.name = 'Channel'
E_do = df.copy()
E_do_g = df.copy()
nbursts = df.copy()
In [ ]:
def print_fit_report(E_pr, gamma=1, leakage=0, dir_ex_t=0, math=True):
"""Print fit and standard deviation for both corrected and uncorrected E
Returns d.E_fit.
"""
E_corr = fretmath.correct_E_gamma_leak_dir(E_pr, gamma=gamma, leakage=leakage, dir_ex_t=dir_ex_t)
E_pr_mean = E_pr.mean()*100
E_pr_delta = (E_pr.max() - E_pr.min())*100
E_corr_mean = E_corr.mean()*100
E_corr_delta = (E_corr.max() - E_corr.min())*100
if math:
display(Math(r'\text{Pre}\;\gamma\quad\langle{E}_{fit}\rangle = %.1f\%% \qquad'
'\Delta E_{fit} = %.2f \%%' % \
(E_pr_mean, E_pr_delta)))
display(Math(r'\text{Post}\;\gamma\quad\langle{E}_{fit}\rangle = %.1f\%% \qquad'
'\Delta E_{fit} = %.2f \%%' % \
(E_corr_mean, E_corr_delta)))
else:
print('Pre-gamma E (delta, mean): %.2f %.2f' % (E_pr_mean, E_pr_delta))
print('Post-gamma E (delta, mean): %.2f %.2f' % (E_corr_mean, E_corr_delta))
In [ ]:
str(ph_sel)
In [ ]:
data_id = '7d'
d7 = loader.photon_hdf5(files_dict[data_id])
d7.calc_bg(**bg_kwargs_auto)
In [ ]:
d7.burst_search(m=10, F=F, ph_sel=ph_sel)
In [ ]:
ds7 = d7.select_bursts(select_bursts.nd, th1=size_min)
dx = ds7
In [ ]:
## KDE Fit
E_do.loc[data_id] = bext.fit_bursts_kde_peak(dx, bandwidth=bandwidth, x_range=E_range_do,
weights=weights)
## Gaussian fit
dx.E_fitter.histogram(binwidth=binwidth, weights=weights)
dx.E_fitter.fit_histogram(mfit.factory_gaussian())
E_do_g.loc[data_id] = dx.E_fitter.params['center']
## D-only selection
do_s = dx.select_bursts(select_bursts.E, E2=0.1)
nbursts.loc[data_id] = do_s.num_bursts
In [ ]:
dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_kde=True, show_kde_peak=True, show_fit_value=True);
plt.xlim(xmax = 0.49)
print_fit_report(E_do.loc[data_id])
In [ ]:
dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_model=True, show_fit_value=True, fit_from='center');
plt.xlim(xmax = 0.49)
print_fit_report(E_do_g.loc[data_id])
In [ ]:
fig, axes = plt.subplots(4, 2, figsize=(12, 8), sharex=True, sharey=True)
fig.subplots_adjust(left=0.08, right=0.96, top=0.93, bottom=0.07,
wspace=0.06, hspace=0.18)
for ich, ax in enumerate(axes.ravel()):
mfit.plot_mfit(dx.E_fitter, ich=ich, ax=ax, plot_model=False, plot_kde=True)
plt.xlim(-0.2, 0.3)
print_fit_report(E_do.loc[data_id])
In [ ]:
fig, axes = plt.subplots(4, 2, figsize=(12, 8), sharex=True, sharey=True)
fig.subplots_adjust(left=0.08, right=0.96, top=0.93, bottom=0.07,
wspace=0.06, hspace=0.18)
for ich, ax in enumerate(axes.ravel()):
mfit.plot_mfit(dx.E_fitter, ich=ich, ax=ax)
plt.xlim(-0.2, 0.3)
print_fit_report(E_do_g.loc[data_id])
In [ ]:
data_id = '12d'
d12 = loader.photon_hdf5(files_dict[data_id])
d12.calc_bg(**bg_kwargs_auto)
In [ ]:
d12.burst_search(m=10, F=F, ph_sel=ph_sel)
In [ ]:
ds12 = d12.select_bursts(select_bursts.nd, th1=size_min)
dx = ds12
In [ ]:
## KDE Fit
E_do.loc[data_id] = bext.fit_bursts_kde_peak(dx, bandwidth=bandwidth, x_range=E_range_do,
weights=weights)
## Gaussian fit
dx.E_fitter.histogram(binwidth=binwidth, weights=weights)
dx.E_fitter.fit_histogram(mfit.factory_gaussian())
E_do_g.loc[data_id] = dx.E_fitter.params['center']
## D-only selection
do_s = dx.select_bursts(select_bursts.E, E2=0.1)
nbursts.loc[data_id] = do_s.num_bursts
In [ ]:
dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_kde=True, show_kde_peak=True, show_fit_value=True);
plt.xlim(xmax = 0.49)
print_fit_report(E_do.loc[data_id])
In [ ]:
dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_model=True, show_fit_value=True, fit_from='center');
plt.xlim(xmax = 0.49)
print_fit_report(E_do_g.loc[data_id])
In [ ]:
data_id = '17d'
d17 = loader.photon_hdf5(files_dict[data_id])
d17.calc_bg(**bg_kwargs_auto)
In [ ]:
d17.burst_search(m=10, F=F, ph_sel=ph_sel)
In [ ]:
ds17 = d17.select_bursts(select_bursts.nd, th1=size_min)
dx = ds17
In [ ]:
## KDE Fit
E_do.loc[data_id] = bext.fit_bursts_kde_peak(dx, bandwidth=bandwidth, x_range=E_range_do,
weights=weights)
## Gaussian fit
dx.E_fitter.histogram(binwidth=binwidth, weights=weights)
dx.E_fitter.fit_histogram(mfit.factory_two_gaussians(p1_center=0.03, p2_center=0.25))
E_do_g.loc[data_id] = dx.E_fitter.params['p1_center']
## D-only selection
do_s = Sel(dx, select_bursts.E, E2=0.1)
nbursts.loc[data_id] = do_s.num_bursts
In [ ]:
dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_kde=True, show_kde_peak=True, show_fit_value=True);
plt.xlim(xmax = 0.49)
print_fit_report(E_do.loc[data_id])
In [ ]:
dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_model=True, show_fit_value=True, fit_from='p1_center');
plt.xlim(xmax = 0.49)
print_fit_report(E_do_g.loc[data_id])
In [ ]:
data_id = 'DO'
do = loader.photon_hdf5(files_dict[data_id])
do.calc_bg(**bg_kwargs_auto)
In [ ]:
do.burst_search(L=10, m=10, F=F, ph_sel=ph_sel)
In [ ]:
dos = do.select_bursts(select_bursts.nd, th1=size_min)
dx = dos
In [ ]:
## KDE Fit
E_do.loc[data_id] = bext.fit_bursts_kde_peak(dx, bandwidth=bandwidth, x_range=E_range_do,
weights=weights)
## Gaussian fit
dx.E_fitter.histogram(binwidth=binwidth, weights=weights)
dx.E_fitter.fit_histogram(mfit.factory_gaussian())
E_do_g.loc[data_id] = dx.E_fitter.params['center']
## D-only selection
nbursts.loc[data_id] = dx.num_bursts
In [ ]:
dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_kde=True, show_kde_peak=True, show_fit_value=True);
plt.xlim(xmax = 0.49)
print_fit_report(E_do.loc[data_id])
In [ ]:
dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_model=True, show_fit_value=True, fit_from='center');
plt.xlim(xmax = 0.49)
print_fit_report(E_do_g.loc[data_id])
In [ ]:
ph_sel
In [ ]:
nbursts = nbursts.astype(int)
nbursts
In [ ]:
E_do_kde = E_do.copy()
In [ ]:
leakage_kde = (E_do_kde / (1 - E_do_kde)).round(6)
leakage_kde
In [ ]:
leakage_gauss = (E_do_g / (1 - E_do_g)).round(6)
leakage_gauss
In [ ]:
sns.set(style='ticks', font_scale=1.4)
In [ ]:
colors = sns.color_palette('Paired', 8)
mpl.rcParams['axes.prop_cycle'] = cycler('color', colors)
sns.palplot(colors)
In [ ]:
fig, ax = plt.subplots()
ax2 = ax.twinx()
kws = dict(lw=2, marker='o', ms=8)
for i, did in enumerate(('7d', '12d', '17d', 'DO')):
(100*leakage_kde).loc[did].plot(label='%s KDE' % did, ax=ax, color=colors[1+i*2], **kws)
nbursts.loc[did].plot(ax=ax2, ls='--', lw=2.5, color=colors[1+i*2])
for i, did in enumerate(('7d', '12d', '17d', 'DO')):
(100*leakage_gauss).loc[did].plot(label='%s Gauss' % did, ax=ax, color=colors[i*2], **kws)
handles, lab = ax.get_legend_handles_labels()
h = handles#[1::2] + handles[::2]
l = lab[1::2] + lab[::2]
ax.legend(ncol=2, loc=1, bbox_to_anchor=(1, 0.5), borderaxespad=0.)
ax.set_ylim(0)
ax2.set_ylim(0, 3200)
plt.xlim(-0.25, 7.25)
plt.xlabel('Channel')
ax.set_ylabel('Leakage %')
ax2.set_ylabel('# Bursts')
sns.despine(offset=10, trim=True, right=False)
In [ ]:
leakage_kde.to_csv('results/Multi-spot - leakage coefficient all values KDE %s.csv' % str(ph_sel))
leakage_gauss.to_csv('results/Multi-spot - leakage coefficient all values gauss %s.csv' % str(ph_sel))
nbursts.to_csv('results/Multi-spot - leakage coefficient all values nbursts %s.csv' % str(ph_sel))
In [ ]:
lk_s = pd.DataFrame(index=['mean', 'std'], columns=E_do.index)
lk_s.loc['mean'] = leakage_kde.mean(1)*100
lk_s.loc['std'] = leakage_kde.std(1)*100
lk_s['mean'] = lk_s.mean(1)
lk_s
In [ ]:
lk_s.T[['std']]
In [ ]:
lk_s.T[['mean']].plot(table=True)
plt.gca().xaxis.set_visible(False) # Hide Ticks
In [ ]:
nbursts
In [ ]:
lk_sw = pd.DataFrame(index=['mean', 'std'], columns=E_do.index)
lk_sw.loc['mean'] = (nbursts*leakage_kde).sum(1)/nbursts.sum(1)*100
lk_sw.loc['std'] = (nbursts*leakage_kde).std(1)/nbursts.sum(1)*100
lk_sw['mean'] = (nbursts.sum(1)*lk_sw).sum(1)/nbursts.sum(1).sum()
lk_sw
In [ ]:
lk_sw.loc['mean'].plot()
ylim(2, 4)
In [ ]:
lk_c = pd.DataFrame(index=['mean', 'std'], columns=E_do.columns)
lk_c.loc['mean'] = leakage_kde.mean()*100
lk_c.loc['std'] = leakage_kde.std()*100
lk_c['mean'] = lk_c.mean(1)
lk_c
In [ ]:
lk_c.loc['mean'].plot()
ylim(2, 4)
In [ ]:
lk_cw = pd.DataFrame(index=['mean', 'std'], columns=E_do.columns)
lk_cw.loc['mean'] = (nbursts*leakage_kde).sum()/nbursts.sum()*100
lk_cw.loc['std'] = (nbursts*leakage_kde).std()/nbursts.sum()*100
lk_cw['mean'] = lk_cw.mean(1)
lk_cw
In [ ]:
lk_cw.loc['mean'].plot()
ylim(2, 4)
In [ ]:
mch_plot_bg(d7)
NOTE: There is a per-channel trend that cannot be ascribed to the background because we performend a D-emission burst search and selection and the leakage vs ch does not resemble the D-background vs channel curve.
The effect is probably due to slight PDE variations (detectors + optics) that slightly change $\gamma$ on a per-spot basis.
In [ ]:
lk = (lk_cw.ix['mean', :-1]*nbursts.mean()).sum()/(nbursts.mean().sum())/100
lk_2 = (lk_sw.ix['mean', :-1]*nbursts.mean(1)).sum()/(nbursts.mean(1).sum())/100
assert np.allclose(lk, lk_2)
print('Mean leakage: %.6f' % lk)