from fretbursts import *

sns = init_notebook()

import lmfit; lmfit.__version__

import phconvert; phconvert.__version__

"""
Explanation: FRETBursts - 8-spot smFRET burst analysis
This notebook is part of a tutorial series for the FRETBursts burst analysis software.

For a step-by-step introduction to FRETBursts usage please refer to 
us-ALEX smFRET burst analysis.
In this notebook we present a typical FRETBursts
workflow for multi-spot smFRET burst analysis. 
Briefly, we show how to perform background estimation, burst search, burst selection, 
FRET histograms, and FRET efficiency fit using different methods. 

Loading the software
End of explanation
"""


url = 'http://files.figshare.com/2182604/12d_New_30p_320mW_steer_3.hdf5'

download_file(url, save_dir='./data')

"""
Explanation: Downloading the sample data file
The complete example dataset can be downloaded 
from here.
Here we download an 8-spot smFRET measurement file using 
the download_file in FRETBursts:
End of explanation
"""


filename = "./data/12d_New_30p_320mW_steer_3.hdf5"

import os
assert os.path.exists(filename)

"""
Explanation: Selecting a data file
End of explanation
"""


d = loader.photon_hdf5(filename)

"""
Explanation: Data load and Burst search
Load and process the data:
End of explanation
"""


d.leakage = 0.038
d.gamma = 0.43

"""
Explanation: For convenience we can set the correction coefficients right away 
so that they will be used in the subsequent analysis. 
The correction coefficients are: 

leakage or bleed-through: leakage
direct excitation: dir_ex (ALEX-only)
gamma-factor gamma

The direct excitation cannot be applied to non-ALEX (single-laser) 
smFRET measurements (like the current one).
End of explanation
"""


d.calc_bg(bg.exp_fit, time_s=30, tail_min_us='auto', F_bg=1.7)
d.burst_search(L=10, m=10, F=7)

"""
Explanation: NOTE: at any later moment, after burst search, a simple 
reassignment of these coefficient will update the burst data 
with the new correction values.

Compute background and burst search:
End of explanation
"""


mch_plot_bg(d)

"""
Explanation: Perform a background plot as a function of the channel:
End of explanation
"""


dplot(d, hist_bg);

"""
Explanation: Let's take a look at the photon waiting times histograms and at the fitted background rates:
End of explanation
"""


dplot(d, timetrace_bg);

"""
Explanation: 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 significant variations during the measurement:
End of explanation
"""


dplot(d, timetrace)
xlim(2, 3); ylim(-100, 100);

"""
Explanation: We can look at the timetrace of the photon stream (binning):
End of explanation
"""


#%matplotlib qt

#dplot(d, timetrace, scroll=True);

#ylim(-100, 100)

#%matplotlib inline

"""
Explanation: We can also open the same plot in an interactive window that allows scrolling (uncomment the following lines):
End of explanation
"""


gamma = d.gamma
gamma

d.gamma = 1
ds = d.select_bursts(select_bursts.size, th1=30, gamma=1)
dplot(ds, hist_fret);

ds = d.select_bursts(select_bursts.size, th1=25, gamma=gamma, donor_ref=False)
dplot(ds, hist_fret);

ds = d.select_bursts(select_bursts.size, th1=25, gamma=gamma)
dplot(ds, hist_fret, weights='size', gamma=gamma);

dplot(ds, scatter_fret_nd_na); ylim(0,200);

"""
Explanation: Burst selection and FRET
Selecting bursts by burst size (select_bursts.size)
End of explanation
"""


ds.gamma = 1.
bext.bursts_fitter(ds, weights=None)
ds.E_fitter.fit_histogram(mfit.factory_two_gaussians(), verbose=False)

"""
Explanation: FRET Fitting
2-Gaussian mixture
Let's fit the $E$ histogram with a 2-Gaussians model:
End of explanation
"""


ds.E_fitter.params

dplot(ds, hist_fret, weights=None, show_model=True,
      show_fit_stats=True, fit_from='p2_center');

"""
Explanation: The fitted parameters are stored in a pandas DataFrame:
End of explanation
"""


# bl.two_gaussian_fit_EM??

EM_results = ds.fit_E_two_gauss_EM(weights=None, gamma=1.)
EM_results

"""
Explanation: Weighted Expectation Maximization
The expectation maximization 
(EM) method is particularly suited to resolve population 
mixtures. Note that the EM algorithm does not fit the histogram 
but the $E$ distribution with no binning.
FRETBursts include a weighted version of the EM algorithm that 
can take into account the burst size.
The algorithm and benchmarks with the 2-Gaussian histogram fit 
are reported here.
You can find the EM algorithm in fretbursts/fit/gaussian_fit.py or typing:
bl.two_gaussian_fit_EM??
End of explanation
"""


ds.fit_E_name, ds.fit_E_res

"""
Explanation: The fitted parameters for each channel are stored in the fit_E_res attribute:
End of explanation
"""


ds.fit_E_model

"""
Explanation: The model function is stored in:
End of explanation
"""


AX = dplot(ds, hist_fret, weights=None)

x = np.r_[-0.2: 1.2 : 0.01]
for ich, (ax, E_fit) in enumerate(zip(AX.ravel(), EM_results)):
    ax.axvline(E_fit, ls='--', color='r')
    ax.plot(x, ds.fit_E_model(x, ds.fit_E_res[ich]))

print('E mean: %.2f%%   E delta: %.2f%%' %\
      (EM_results.mean()*100, (EM_results.max() - EM_results.min())*100))

"""
Explanation: Let's plot the histogram and the model with parameters from the EM fit:
End of explanation
"""


import pandas as pd

EM_results = pd.DataFrame(ds.fit_E_res, columns=['p1_center', 'p1_sigma', 'p2_center', 'p2_sigma', 'p1_amplitude'])
EM_results * 100

ds.E_fitter.params * 100

"""
Explanation: Comparing 2-Gaussian and EM fit
To quickly compare the 2-Gaussians with the EM fit we convert the EM fit results in a DataFrame:
End of explanation
"""


(ds.E_fitter.params - EM_results) * 100

"""
Explanation: And we compute the difference between the two sets of parameters:
End of explanation
"""


