#!/usr/bin/env python
# coding: utf-8

# # FRETBursts - μs-ALEX smFRET burst analysis
# 
# *This notebook is part of a [tutorial series](https://github.com/OpenSMFS/FRETBursts_notebooks) for the [FRETBursts](http://opensmfs.github.io/FRETBursts/) burst analysis software.*

# > In this notebook, we present a typical [FRETBursts](http://opensmfs.github.io/FRETBursts/) 
# > workflow for **μs-ALEX smFRET burst analysis**. 
# > Briefly, we show how to perform background estimation, burst search, burst selection, 
# > compute FRET histograms and ALEX histograms, do sub-population selection and finally, FRET efficiency fit.
# 
# 
# <br>
# <div class="alert alert-success">
# If you are new to the notebook interface see 
# <a href="1. First Steps - Start here if new to Jupyter Notebooks.ipynb">First Steps</a>
# before continuing.
# </div>
# 
# Before running the notebook, you can click on menu *Cell* -> *All Output* -> *Clear* 
# to clear all previous output. This will avoid mixing output from current execution and 
# the previously saved one.

# # Loading the software

# We start by loading **`FRETBursts`**:

# In[1]:


from fretbursts import *


# <div class="alert alert-info">
# Thanks in advance for remembering to <b>cite</b> FRETBursts in publications or presentations!
# </div>
# 
# Note that FRETBursts **version string** tells you the exact FRETBursts version (and revision) in use.
# Storing the version in the notebook helps with reproducibility and 
# tracking software regressions.
# 
# Next, we initialize the default plot style for the notebook 
# (under the hood it uses [seaborn](http://stanford.edu/~mwaskom/software/seaborn/)):

# In[2]:


sns = init_notebook()


# Note that the previous command has no output. Finally, we print the version of some dependencies:

# In[3]:


import lmfit; lmfit.__version__


# In[4]:


import phconvert; phconvert.__version__


# # Downloading the data file

# The full list of smFRET measurements used in the [FRETBursts tutorials](https://github.com/OpenSMFS/FRETBursts_notebooks) 
# can be found on [Figshare](http://dx.doi.org/10.6084/m9.figshare.1456362).
# 
# This is the file we will download:

# In[5]:


url = 'http://files.figshare.com/2182601/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5'


# <div class="alert alert-success">
# You can change the <code>url</code> variable above to download your own data file.
# This is useful if you are executing FRETBursts online and you want to use
# your own data file. See 
# <a href="1. First Steps - Start here if new to Jupyter Notebooks.ipynb">First Steps</a>.
# </div>
# 
# Here, we download the data file and put it in a folder named `data`, 
# inside the notebook folder:

# In[6]:


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


# > **NOTE**: If you modified the `url` variable providing an invalid URL
# > the previous command will fail. In this case, edit the cell containing 
# > the `url` and re-try the download.

# # Selecting the data file

# Use one of the following 2 methods to select a data file.

# ## Option 1: Paste the file-name

# Here, we can directly define the file name to be loaded:

# In[7]:


filename = "./data/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5"
filename


# Now `filename` contains the path of the file you just selected. 
# Run again the previous cell to select a new file. In a following cell
# we will check if the file actually exists.

# <div class="alert alert-success">
# When running the notebook online and using your own data file,
# make sure to modify the previous cell. 
# <br>
# See
# <a href="1. First Steps - Start here if new to Jupyter Notebooks.ipynb">First Steps</a>.
# </div>

# ## Option 2: Use an "Open File" dialog

# Alternatively, you can select a data file with an "Open File" windows.
# Note that, since this only works in a local installation, the next commands
# are commented (so nothing will happen when running the cell). 
# 
# If you want to try the File Dialog, you need to remove the `#` signs:

# In[8]:


# filename = OpenFileDialog()
# filename


# ## Check that the data file exists

# Let's check that the file exists:

# In[9]:


import os
if os.path.isfile(filename):
    print("Perfect, file found!")
else:
    print("Sorry, file:\n%s not found" % filename)


# # Load the selected file

# We can finally load the data and store it in a variable called `d`:

# In[10]:


d = loader.photon_hdf5(filename)


# If you don't get any message, the file is loaded successfully.
# 
# We can also set the 3 correction coefficients: 
# 
# * leakage or bleed-through: `leakage`
# * direct excitation: `dir_ex` (ALEX-only)
# * gamma-factor `gamma`

# In[11]:


d.leakage = 0.11
d.dir_ex = 0.04
d.gamma = 1.


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

# # μs-ALEX parameters
# 
# At this point, timestamps and detectors numbers are contained in the `ph_times_t` and `det_t` attributes of `d`. Let's print them:

# In[12]:


d.ph_times_t, d.det_t


# We need to define some ALEX parameters: 

# In[13]:


d.add(det_donor_accept = (0, 1), 
      alex_period = 4000,
      offset = 700,
      D_ON = (2180, 3900), 
      A_ON = (200, 1800))


# Here the parameters are:
# 
# - `det_donor_accept`: donor and acceptor channels
# - `alex_period`: length of excitation period (in timestamps units)
# - `D_ON` and `A_ON`: donor and acceptor excitation windows
# - `offset`: the offset between the start of alternation and start of timestamping 
#   (see also [Definition of alternation periods](http://photon-hdf5.readthedocs.org/en/latest/phdata.html#definition-of-alternation-periods)).
# 
# To check that the above parameters are correct, we need to plot the histogram of timestamps (modulo the alternation period) and superimpose the two excitation period definitions to it:

# In[14]:


bpl.plot_alternation_hist(d)


# If the previous alternation histogram looks correct, 
# the corresponding definitions of the excitation periods can be applied to the data using the following command:

# In[15]:


loader.alex_apply_period(d)


# If the previous histogram does not look right, the parameters in the `d.add(...)` cell can be modified and checked by running the histogram plot cell until everything looks fine. Don't forget to apply the 
# parameters with `loader.usalex_apply_period(d)` as a last step.
# 
# > **NOTE:** After applying the ALEX parameters a new array of 
# > timestamps containing only photons inside the excitation periods 
# > is created (name `d.ph_times_m`). To save memory, by default, 
# > the old timestamps array (`d.ph_times_t`) is deleted. Therefore, 
# > in the following, when we talk about all-photon selection we always 
# > refer to all photons inside both excitation periods.

# ## Measurement infos

# The entire measurement data is now stored in the variable `d`. Printing it
# will give a compact representation containing the file-name and additional parameters:

# In[16]:


d


# To check the **measurement duration** (in seconds) run:

# In[17]:


d.time_max


# # Plotting basics
# 
# In this section basic concepts of plotting with FRETBursts using the
# timetrace plot as an example.
# 
# To plot a timetrace of the measurement we use:

# In[18]:


dplot(d, timetrace);


# Here, `dplot` is a generic wrapper (the same for all plots)
# that takes care of setting up the figure, title and axis 
# (in the multispot case `dplot` creates multi-panel plot).
# The second argument, `timetrace`, is the actual plot function. 
# All the eventual additional arguments passed to `dplot` are, 
# in turn, passed to the plot function (e.g. `timetrace`).
# 
# If we look at the documentation for [`timetrace`](http://fretbursts.readthedocs.org/en/latest/plots.html#fretbursts.burst_plot.timetrace)
# function we notice that it accepts a long list of arguments.
# In python, when an argument is not specified, it will take the default
# value specified in the function definition (see previus link).
# 
# As an example, to change the bin size (i.e. duration) of the timetrace histogram,
# we can look up in the [`timetrace` documentation](http://fretbursts.readthedocs.org/en/latest/plots.html#fretbursts.burst_plot.timetrace)
# and find that the argument we need to modify is `binwidth`
# (we can also see that the default value is `0.001` seconds).
# We can then re-plot the timetrace using a bin of 0.5 ms:

# In[19]:


dplot(d, timetrace, binwidth=0.5e-3);


# The timetrace is **computed** between `tmin` and `tmax` (by default 0 and 200s),
# but as you can see is displayed only between 0 an 1 second, just because these
# are the default x-axis limits. The axis limits can be changes by using the
# standard matplotlib command `plt.xlim()`.
# On the other hand, to change the range where the timetrace is **computed**,
# we pass the additional arguments `tmin` and `tmax` as follows:

# In[20]:


dplot(d, timetrace, binwidth=0.5e-3, tmin=50, tmax=150)
plt.xlim(51, 52);


# When using FRETBursts in a notebook, all plots are static by default.
# This is because we use the so called `inline` backend of matplotlib.
# If you want to manipulate figures interactively, you can switch
# to the interactive `notebook` backend with:
# 
# ```
# %matplotlib notebook
# ```
# 
# to go back to inline use:
# 
# 
# ```
# %matplotlib inline
# ```
# 
# **NOTE**: Currently, the `notebook` backend is incompatible with the QT backend. 
# If in a session you activate the `notebook` backend, then switching to the QT backend requires 
# restarting the notebook. Conversely, you can switch between `inline` and `notebook` 
# or between `inline` and `qt4` backends in the same session wihtou issues.
# 
# > ### See also:
# >
# > - [bpl.timetrace](http://fretbursts.readthedocs.org/en/latest/plots.html#fretbursts.burst_plot.timetrace)
# >    function documentation
# > - [bpl.ratetrace](http://fretbursts.readthedocs.org/en/latest/plots.html#fretbursts.burst_plot.ratetrace)
# >   function documentation
# > - [Intensity timetrace and Rate-timetrace](#Intensity-timetrace-and-Rate-timetrace), a later section in this tutorial.

# # Background estimation
# 
# As a first step of the analysis, we need to estimate the background. 
# Here we will compute the background using the recommended approach of using
# the auto-threshold.
# 
# For more details see [Background estimation](Example - Background estimation.ipynb).
# 

# ## Automatic threshold
# 
# It is a good practice to monitor background rates as a function of time.
# Here, we compute background in adjacent 30s windows (called *background periods*)
# and plot the estimated rates as a function of time.

# In[21]:


d.calc_bg(bg.exp_fit, time_s=30, tail_min_us='auto')


# In[22]:


dplot(d, timetrace_bg)


# > **NOTE**: All background data is stored in `d.bg`. For details on how to
# > to export it see the [Background estimation](Example - Background estimation.ipynb) notebook.

# # Burst analysis

# The first step of burst analysis is the burst search. 
# 
# We will use the sliding-window algorithm on all photons. Note 
# that "all photons", as mentioned before, means all photons selected in the 
# alternation histogram.
# An important variation compared to the classical sliding-windows 
# is that the threshold-rate for burst start is computed as 
# a function of the background and changes when the background
# changes during the measurement.
# 
# To perform a burst search evaluating the photon rate with 
# 10 photons (`m=10`), and selecting a minimum rate 6 times larger than 
# the background rate (F=6) calculated with all photons (default):

# In[23]:


d.burst_search(L=10, m=10, F=6)


# The previous command performs the burst search, corrects 
# the bursts sizes for background, spectral leakage and direct excitation,
# and computes $\gamma$-corrected FRET and Stoichiometry. 
# 
# See the 
# [`burst_search` documentation](http://fretbursts.readthedocs.org/en/latest/data_class.html#fretbursts.burstlib.Data.burst_search) for more details.

# We can plot the resulting FRET histogram using the following command:

# In[24]:


dplot(d, hist_fret);


# All pre-defined plots follow this pattern: 
# call the generic `dplot()` function, passing 2 parameters:
# 
# - the measurement data (`d` in this case) 
# - the plot function (`hist_fret`)
# 
# In some case we can add other optional parameters to tweak the plot.
# 
# All plot functions start with `hist_` for histograms, 
# `scatter_` for scatter-plots  or `timetrace_` for plots as a function 
# of measurement time. You can use autocompletion to find all 
# plot function or you can look in `bursts_plot.py` where 
# all plot functions are defined.

# Instead of `hist_fret` we can use `hist_fret_kde` to add a [KDE](http://en.wikipedia.org/wiki/Kernel_density_estimation) overlay. Also, we can plot a **weighted histogram** by passing an additional parameter `weights`:

# In[25]:


dplot(d, hist_fret, show_kde=True);
dplot(d, hist_fret, show_kde=True, weights='size');


# You can experiment with different weighting schema (for all 
# supported weights see `get_weigths()` function in `fret_fit.py`).

# ## Burst selection

# When we performed the burst search, we specified `L=10` without 
# explaining what this parameter means. *L* is traditionally the minimum size 
# (number of photons) for a burst: smaller bursts will be rejected. 
# By setting L=m (10 in this case) we are deciding to not discard 
# any burst (because the smallest detected burst has at least *m* counts). 
# 
# Selecting the bursts in a second step, by applying a minimum burst size criterion, 
# results in a more accurate and unbiased selection.
# 
# For example, we can select bursts with more than 30 photons (after 
# background, gamma, leakage and direct excitation corrections)
# and store the result in a new
# `Data()` variable `ds`:

# In[26]:


ds = d.select_bursts(select_bursts.size, th1=30)


# By defaults the burst size includes donor and acceptor photons 
# during donor excitation. To add acceptor photons during 
# acceptor excitation (`naa`), we add the parameter `add_naa=True`:

# In[27]:


ds = d.select_bursts(select_bursts.size, add_naa=True, th1=30)


# Similar to plot functions, all selection functions 
# are defined in `select_bursts.py` and you can access them by typing 
# `select_bursts.` and using the TAB key for autocompletion. 
# 
# > **See also:** 
# > * [Burst selection](http://fretbursts.readthedocs.org/en/latest/burst_selection.html) in the documentation. 
# > In particular the function [`select_bursts.size`](http://fretbursts.readthedocs.org/en/latest/burst_selection.html#fretbursts.select_bursts.size) and [`Data.select_bursts`](http://fretbursts.readthedocs.org/en/latest/data_class.html#fretbursts.burstlib.Data.select_bursts).
# 
# To replot the FRET histogram after selection (note that now 
# we are passing `ds` to the plot function):

# In[28]:


dplot(ds, hist_fret);


# Note how the histogram exhibits much more clearly defined peaks after burst selection. 

# ## Histogram fitting and plotting style
# 
# Under the hood the previous `hist_fret` plot creates a `MultiFitter` 
# object for $E$ values. This object, stored as `ds.E_fitter`, operates
# on multi-channel data and computes the histogram, KDE and can fit 
# the histogram with a model ([lmfit.Model](http://lmfit.github.io/lmfit-py/model.html)).
# 
# Now, just for illustration purposes, we fit the previous histogram with 3 Gaussians, using the already created `ds.E_fitter` object:

# In[29]:


ds.E_fitter.fit_histogram(mfit.factory_three_gaussians(), verbose=False)


# In[30]:


dplot(ds, hist_fret, show_model=True);


# The bin width can be changed with `binwidth` argument. Alternatively,
# an arbitrary array of bin edges can be passed in `bins` 
# (overriding `binwidth`).
# 
# We can customize the appearance of this plot (type 
# `hist_fret?` for the complete set of arguments). 
# 
# For example to change from a bar plot to a line-plot 
# we use the `hist_style` argument:

# In[31]:


dplot(ds, hist_fret, show_model=True, hist_style='line')


# We can customize the line-plot, bar-plot, the model 
# plot and the KDE plot by passing dictionaries with matplotlib 
# style. The name of the arguments are:
# 
# - `hist_plot_style`: style for the histogram line-plot
# - `hist_bar_style`: style for the histogram bar-plot
# - `model_plot_style`: style for the model plot
# - `kde_plot_style`: style for the KDE plot
# 
# As an example:

# In[32]:


dplot(ds, hist_fret, show_model=True, hist_style='bar', show_kde=True,
      kde_plot_style = dict(linewidth=5, color='orange', alpha=0.6),
      hist_plot_style = dict(linewidth=3, markersize=8, color='b', alpha=0.6))
plt.legend();


# # Other plots

# Similarly, we can plot the burst size using all photons 
# (type `hist_size?` to learn about all plot options):

# In[33]:


dplot(ds, hist_size, add_naa=True);


# Or plot the burst size histogram for the different components:

# In[34]:


dplot(ds, hist_size_all);


# > **NOTE:** The previous plot may generate a benign warning 
# > due to the presence of zeroes when switching to log scale. Just ignore it.

# A scatterplot of Size *vs* FRET is created by:

# In[35]:


dplot(ds, scatter_fret_nd_na)
xlim(-1, 2)


# # Study of different populations
# 
# ## Removing multi-mers

# We can further select only bursts smaller than 300 photons 
# to get rid of multi-molecule events:

# In[36]:


ds2 = ds.select_bursts(select_bursts.size, th2=300)


# and superimpose the two histograms before and after selection to see the difference:

# In[37]:


ax = dplot(ds2, hist_fret, hist_style='bar', show_kde=True, 
              hist_bar_style = dict(facecolor='r', alpha=0.5, 
                                    label='Hist. no large bursts'),
              kde_plot_style = dict(lw=3, color='m', 
                                    label='KDE no large bursts'))
dplot(ds, hist_fret, ax=ax, hist_style='bar', show_kde=True,
      hist_bar_style = dict(label='Hist. with large bursts'),
      kde_plot_style = dict(lw=3, label='KDE with large bursts'))
plt.legend();


# > **NOTE:** It is not necessarily true that bursts with more that 300 photons 
# > represents multiple molecules. 
# > To asses the valididty of this assumption it can be useful to 
# > plot the peak count rates in each burst. See `hist_burst_phrate` 
# > for this kind of plot.

# ## Fit and plot peak positions
# 
# We can find the KDE peak position in a range (let say 0.2 ... 0.6):

# In[38]:


ds.E_fitter.find_kde_max(np.r_[0:1:0.0002], xmin=0.2, xmax=0.6)


# and plot it with `show_kde_peak=True`, we also use `show_fit_value=True` to show a box with the fitted value:

# In[39]:


dplot(ds, hist_fret, hist_style='line', 
      show_fit_value=True, 
      show_kde=True, show_kde_peak=True);


# Instead of using the KDE, we can use the peak position as fitted from a gaussian model. 

# In[40]:


ds.E_fitter.fit_histogram(mfit.factory_three_gaussians(), verbose=False)


# To select which peak to show we use `fit_from='p1_center'`:

# In[41]:


dplot(ds, hist_fret, hist_style='line', 
      show_fit_value=True, fit_from='p2_center', 
      show_model=True);


# The string `'p2_center'` is the name of the parameter of the 
# gaussian fit that we want to show in the text box. To see all 
# the parameters of the model we look in:

# In[42]:


ds.E_fitter.params  # <-- pandas DataFrame, one row per channel


# ## ALEX plots

# We can create a simple E-S scatter plot with `scatter_alex`:

# In[43]:


dplot(ds, scatter_alex, figsize=(4,4), mew=1, ms=4, mec='black', color='purple');


# We can also plot the ALEX histogram with a scatterplot overlay using `hist2d_alex`:

# In[44]:


dplot(ds, hist2d_alex);


# Finally we can also plot an ALEX histogram and marginals 
# (joint plots) as follow (for more options see: 
# [Example - usALEX histogram](Example - usALEX histogram.ipynb)):

# In[45]:


alex_jointplot(ds);


# To get rid of the large donor-only population, we can simply 
# select bursts with at least 10 photons in the acceptor channel 
# (during acceptor excitation). At the same time,
# with a burst size selection using Dex photons we get rid
# of the A-only population:

# In[46]:


ds = d.select_bursts(select_bursts.size, th1=20)
ds2 = ds.select_bursts(select_bursts.naa, th1=10)
alex_jointplot(ds2);


# As you can see, we remained with only the FRET sub-populations.
# 
# See next sections show how to select a region on the E/S 
# histogram to isolate a subpopulation.

# ## Graphical selection of an E-S region
# <br>
# <div class="alert alert-warning">
# The graphical selection of E-S regions requires a local FRETBursts installation.
# Therefore the commands below are commented by default.
# <br><br>
# If you have a local installation and you want to try commands below,
# please uncomment (i.e. remove the inital <code>#</code>) from the lines
# containing the <code>%matplotlib</code> command.
# </div>
# 
# 
# To select bursts graphically, we need to open the ALEX histogram 
# in a new (QT) window, drag the mouse to define a selection and 
# have it printed here in the notebook.
# 
# Here we describe how to do it in 3 steps.
# 
# **Step 1** Switch the plot modality to QT (i.e. opens graphs in a separate window):

# In[47]:


# Switches to open plot in external window
#%matplotlib qt


# **Step 2** Plot the E-S histogram in an external windows. There you can
# use the mouse to select a region:

# In[48]:


# ALEX histogram with GUI selection enabled
dplot(ds, hist2d_alex, gui_sel=True)


# The region selection is printed above.
# 
# **Step 3** Restore the normal inline plotting (no more external windows).

# In[49]:


# Switch back to show plots inline in the notebook
#%matplotlib inline


# ## Selection bursts by E-S values

# To apply a selection based on E-S values, 
# we can paste the values obtained from the previous plot
# (or we can type them in manually).
# 
# The selection function used here is `select_bursts.ES`.
# The same *E* and *S* boundaries can define either a rectangular 
# or an elliptical selection by using 
# respectively `rect=True` or `rect=False`. 
# Here we use the elliptical selection:

# In[50]:


roi = dict(E1=-0.07, E2=1.17, S1=0.18, S2=0.70, rect=False)
d_fret_mix = ds.select_bursts(select_bursts.ES, **roi)


# By plotting the FRET histogram we can double check that 
# the selection has been applied:

# In[51]:


g = alex_jointplot(d_fret_mix)
bpl.plot_ES_selection(g.axes[0], **roi);


# Now we can further separate high- and low FRET sub-populations. 

# In[52]:


roi_high_fret = dict(E1=0.65, E2=1.09, S1=-0.13, S2=0.96)
d_high_fret = d_fret_mix.select_bursts(select_bursts.ES, **roi_high_fret)
roi_low_fret = dict(E1=-0.19, E2=0.64, S1=-0.05, S2=0.92)
d_low_fret = d_fret_mix.select_bursts(select_bursts.ES, **roi_low_fret)


# In[53]:


alex_jointplot(d_high_fret);
alex_jointplot(d_low_fret);


# We can for example compute the ratio of high- and low-fret bursts:

# In[54]:


d_low_fret.num_bursts / d_high_fret.num_bursts


# ## Burst Width analysis

# To plot a burst-width histogram, we use `hist_width` instead of `hist_fret`:

# In[55]:


dplot(d_low_fret, hist_width)


# To use a larger bin size, plots two sub-populations (in different color) and add a legend:

# In[56]:


ax = dplot(d_high_fret, hist_width, bins=(0, 10, 0.2))
dplot(d_low_fret, hist_width, bins=(0, 10, 0.2), ax=ax)
plt.legend(['High-FRET population', 'Low-FRET population']);


# Finally, we compute the mean burst width for each subpopulation:

# In[57]:


millisec = d.clk_p * 1e3
mean_b_width_low_fret  = d_low_fret.mburst[0].width.mean() * millisec
mean_b_width_high_fret = d_high_fret.mburst[0].width.mean() * millisec

print('Mean burst width: %.1f ms (high-fret), %.1f (low-fret)' % 
      (mean_b_width_high_fret, mean_b_width_low_fret))


# # FRET fit: in-depth example
# 
# We can fit a FRET distribution to any model. For example, 
# we will fit the FRET selection (2 FRET sub-populations) with 2 Gaussians. 
# The fitting model is a [`Model` object](http://lmfit.github.io/lmfit-py/model.html) 
# from the [`lmfit` library](http://lmfit.github.io/lmfit-py/).
# 
# The first step, previously performed implicitely by the `hist_fret()`
# plot function, is to create a [MultiFitter](http://fretbursts.readthedocs.org/en/latest/mfit.html#the-multifitter-class) 
# object for $E$ and/or $S$, by calling [bext.bursts_fitter()](http://fretbursts.readthedocs.org/en/latest/plugins.html?highlight=bursts_fitter#fretbursts.burstlib_ext.bursts_fitter). 
# With `MultiFitter` objects, we can compute histograms, 
# KDEs and fit the histogram in one single step, as in the following example:

# In[58]:


E_fitter = bext.bursts_fitter(d_fret_mix, 'E', binwidth=0.03, bandwidth=0.03, 
                              model=mfit.factory_two_gaussians())
S_fitter = bext.bursts_fitter(d_fret_mix, 'S', binwidth=0.03, bandwidth=0.03, 
                              model=mfit.factory_gaussian())


# However if we want to modify the model (for example to add a 
# constrain) we need to perform the fit in a second step. 
# To skip the fitting, we simply avoid passing a `model`:

# In[59]:


E_fitter = bext.bursts_fitter(d_fret_mix, 'E', binwidth=0.03, bandwidth=0.03)
S_fitter = bext.bursts_fitter(d_fret_mix, 'S', binwidth=0.03, bandwidth=0.03) 


# Now we create a model and initialize the parameters 
# using [`mfit.factory_two_gaussians()`](http://fretbursts.readthedocs.org/en/latest/mfit.html#fretbursts.mfit.factory_two_gaussians) 
# (see [Model factory functions](http://fretbursts.readthedocs.org/en/latest/mfit.html#model-factory-functions) 
# for a list of pre-defined models in FRETBursts):

# In[60]:


model = mfit.factory_two_gaussians(add_bridge=True)


# We can see the list of parameters, initial values and constraints:

# In[61]:


model.print_param_hints()


# We can change initial values and/or constrains (bounds):

# In[62]:


model.set_param_hint('p1_center', value=0.3, min=-0.1, max=0.6)
model.set_param_hint('p2_center', value=0.85, min=0.5, max=1.1)
model.print_param_hints()


# Finally, we fit the histogram with one 
# of the [supported minimization methods](http://lmfit.github.io/lmfit-py/fitting.html#fit-engines-label) 
# (the default is least-squares):

# In[63]:


E_fitter.fit_histogram(model=model, verbose=False)     # default method is 'leastsq'
#E_fitter.fit_histogram(model=model, method='nelder')  # example using simplex method


# To plot the model with the fitted parameters on top of the FRET 
# histogram, we add `show_model=True` as seen before:

# In[64]:


dplot(d_fret_mix, hist_fret, show_model=True);


# The fitting results of `lmfit` models are stored in [lmfit.ModelResult](https://lmfit.github.io/lmfit-py/model.html#the-modelresult-class) objects. In FRETBursts, these objects are available in `MultiFitter.fit_res`:

# In[65]:


E_fitter.fit_res[0]


# Here the `[0]` indicates CH=0. This index is used in multispot measurements, in which 
# there are several channels.
# 
# The values of the fitted parameters are available in the `best_values` 
# dictionary:

# In[66]:


results = E_fitter.fit_res[0]


# In[67]:


results.params.pretty_print(columns=['value'])


# We can also print a complete report of fitted parameters 
# including reduced $\chi^2$, error ranges 
# ($\pm 1 \sigma$) and correlations:

# In[68]:


print(results.fit_report(min_correl=0.5))


# To customize the printed report, see [lmfit.fit_report()](https://lmfit.github.io/lmfit-py/fitting.html?highlight=fit_report#getting-and-printing-fit-reports) docs.

# We can also take a look at the initial parameters:

# In[69]:


results.init_params.pretty_print()


# We can compute more accurate confidence-intervals (note that 
# it can take several seconds, 
# see [lmfit docs](https://lmfit.github.io/lmfit-py/confidence.html) 
# for the method details):

# In[70]:


ci = results.conf_interval()
lmfit.report_ci(ci)


# Finally, FRETBursts's 
# [MultiFitter object](http://fretbursts.readthedocs.org/en/latest/mfit.html#the-multifitter-class) 
# (e.g. `E_fitter`) stores the fitted parameters in `E_fitter.params` as pandas DataFrame:

# In[71]:


E_fitter.params


# For more information on fitting see:
# 
# - [FRETBursts: Fit Framework](http://fretbursts.readthedocs.org/en/latest/fit.html)
# - [lmfit Documentation](http://cars9.uchicago.edu/software/python/lmfit/index.html).

# # Exporting Data
# 
# To export data computed by FRETBursts, you need to find first where 
# the data is stored or, for data computed on fly, which function/method 
# is used to computed it. 
# 
# Most data computed by FRETBursts is stored in the `Data` object
# (see the [Data docs](http://fretbursts.readthedocs.org/en/latest/data_class.html) for the list of attributes).
# For example:
# 
# - corrected burst counts are stored in `Data.nd`, `Data.na`, `Data.naa`
# - timestamps are stored in `Data.ph_times_m`
# 
# All these attributes are list of arrays, one per excitation spot. 
# This means that, even for single-spot data, we need to use indexing (`[0]`)
# to obtain the array. For example:

# In[72]:


d.nd[0]


# Any numpy array can be saved to disk with the method `.tofile()`. For example:

# In[73]:


d.nd[0].tofile('n_dd.csv', sep=',')


# The background data is stored in `Data.bg` (and `Data.bg_mean`) attribute (see [Getting the background rates](#Getting-the-background-rates) in this notebook). This attribute is dict containing lists containing arrays (scalars) and ca be saved to disk similarly.
# 

# ## Exporting burst data
# 
# Part of burst data is stored in `Data` attributes (`mburst`, `nd`, `na`, `naa`)
# and part is computed on fly (duration, raw counts, etc.).
# 
# To put all burst data in a single "table" (a `pandas.DataFrame`), one row
# per burst, we can use `bext.burst_data`:

# In[74]:


bursts = bext.burst_data(ds)
bursts


# > **NOTE**: For multi-spot data use the `ich` argument to get burst data
# > for the different spots.
# 
# This table (as any `pandas.DataFrame`) can be saved in a CSV text file with:

# In[75]:


bursts.to_csv('bursts.csv')


# # Intensity timetrace and Rate timetrace
# 
# For an initial inspection of a data file, it is common to compute an intensity *timetrace* plot (or timetrace for short). 
# We also can compute a similar plot called *ratetrace*, which does not bin counts 
# but shows the instantaneous count rate. In both cases, it is convenient to scroll 
# along the plot interactively.
# 
# In FRETBursts, the two *timetrace* and *ratetrace* plots support 
# interactive scrolling. We just need to switch from the inline backend 
# to QT, as we did before. 
# 
# <br>
# <div class="alert alert-warning">
# The graphical scrolling of timetraces requires a local FRETBursts installation.
# Therefore the commands below are commented by default.
# <br><br>
# If you have a local installation and you want to try commands below,
# please uncomment them by removing the inital <code>#</code>.
# </div>
# 
# 
# Here are the corresponding commands:

# In[76]:


#%matplotlib qt


# In[77]:


#dplot(ds, ratetrace, scroll=True, bursts=True)


# In[78]:


#dplot(ds, timetrace, tmax=600, scroll=True, bursts=True)


# In[79]:


#%matplotlib inline


# Note that, to save memory, the previous plots stops will show the first 200 s
# of measurement.
# For both `timetrace` and `ratetrace`, you can extent (or reduce) the plotted 
# region passing the `tmin` and `tmax` arguments (values in seconds).
# 
# > ### Related documentation:
# >
# > - [bpl.timetrace](http://fretbursts.readthedocs.org/en/latest/plots.html#fretbursts.burst_plot.timetrace)
# > - [bpl.ratetrace](http://fretbursts.readthedocs.org/en/latest/plots.html#fretbursts.burst_plot.ratetrace)

# <br>
# <div class="alert alert-success">
# This ends of the tutorial.
# </div>
# 
# ----
# 
# *Please send feedback and/or question by opening a new [GitHub Issue](https://github.com/OpenSMFS/FRETBursts/issues).*