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

# # Example - Burst Variance Analysis
# 
# *This notebook is part of smFRET burst analysis software [FRETBursts](http://tritemio.github.io/FRETBursts/).*

# > This notebook shows how to implement Burst Variance Analysis (BVA) ([Torella 2011](http://dx.doi.org/10.1016/j.bpj.2011.01.066)) using FRETBursts.
# 
# > For a complete tutorial on burst analysis see 
# > [FRETBursts - us-ALEX smFRET burst analysis](FRETBursts - us-ALEX smFRET burst analysis.ipynb).

# # Loading the software

# We start loading the **`FRETBursts`** software:

# In[1]:


from fretbursts import *
sns = init_notebook()


# # Data file

# In[2]:


url = 'http://files.figshare.com/2182601/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5'
download_file(url, save_dir='./data')


# In[3]:


file_name = "0023uLRpitc_NTP_20dT_0.5GndCl.hdf5"

# Here the folder is the subfolder "data" of current notebook folder
folder_name = './data/'
full_fname = folder_name + file_name
full_fname


# In[4]:


import os
if os.path.isfile(full_fname):
    print ("Perfect, I found the file!")
else:
    print ("Sorry, I can't find the file:\n", full_fname)


# # Load the selected file

# In[5]:


d = loader.photon_hdf5(full_fname)


# # us-ALEX parameters
# 
# At this point, in `d`, we only have the timestamps (`ph_times_t`) and the detector numbers (`det_t`):

# In[6]:


d.add(det_donor_accept=(0, 1), 
      alex_period=4000, 
      D_ON=(2100, 3900), 
      A_ON=(150, 1900),
      offset=700)
bpl.plot_alternation_hist (d)
loader.usalex_apply_period(d)


# # Burst Search and Selection

# Here we perform a standard burst search using the donor-excitation photon stream 
# (`Ph_sel(Dex='DAem')`) as required for BVA.
# Then, we apply a few selection filters to discard D-only and A-only bursts.

# In[7]:


d.calc_bg(bg.exp_fit, time_s=50.1, tail_min_us='auto', F_bg=1.7)


# In[8]:


d.burst_search(m=10, computefret=False, ph_sel=Ph_sel(Dex='DAem'))
d.calc_fret(count_ph=True, corrections=False)


# In[9]:


ds = d.select_bursts(select_bursts.naa, th1=30, computefret=False)
ds1 = ds.select_bursts(select_bursts.size, th1=30, computefret=False)
ds_FRET = ds1.select_bursts(select_bursts.S, S1=0.25, S2=0.85, computefret=False)


# In[10]:


dx=ds_FRET
alex_jointplot(dx)


# # Burst Variance Analysis

# We define a function to compute $s_E$:

# In[11]:


def bva_sigma_E(n, bursts, DexAem_mask, out=None):
    """
    Perform BVA analysis computing std.dev. of E for sub-bursts in each burst.
    
    Split each burst in n-photons chunks (sub-bursts), compute E for each sub-burst,
    then compute std.dev. of E across the sub-bursts.

    For details on BVA see:

    - Torella et al. (2011) Biophys. J. doi.org/10.1016/j.bpj.2011.01.066
    - Ingargiola et al. (2016) bioRxiv, doi.org/10.1101/039198

    Arguments:
        n (int): number of photons in each sub-burst
        bursts (Bursts object): burst-data object with indexes relative 
            to the Dex photon stream.
        DexAem_mask (bool array): mask of A-emitted photons during D-excitation 
            periods. It is a boolean array indexing the array of Dex timestamps 
            (`Ph_sel(Dex='DAem')`).
        out (None or list): append the result to the passed list. If None,
            creates a new list. This is useful to accumulate data from
            different spots in a single list.

    Returns:
        E_sub_std (1D array): contains for each burst, the standard deviation of 
        sub-bursts FRET efficiency. Same length of input argument `bursts`.
    """
    E_sub_std = [] if out is None else out
    
    for burst in bursts:
        E_sub_bursts = []
        startlist = range(burst.istart, burst.istop + 2 - n, n)
        stoplist = [i + n for i in startlist]
        for start, stop in zip(startlist, stoplist):
            A_D = DexAem_mask[start:stop].sum()
            assert stop - start == n
            E = A_D / n
            E_sub_bursts.append(E)
        E_sub_std.append(np.std(E_sub_bursts))
        
    return E_sub_std


# Next we prepare the data for BVA:

# In[12]:


ph_d = ds_FRET.get_ph_times(ph_sel=Ph_sel(Dex='DAem'))
bursts = ds_FRET.mburst[0]
bursts_d = bursts.recompute_index_reduce(ph_d)


# In[13]:


Dex_mask = ds_FRET.get_ph_mask(ph_sel=Ph_sel(Dex='DAem'))   
DexAem_mask = ds_FRET.get_ph_mask(ph_sel=Ph_sel(Dex='Aem')) 
DexAem_mask_d = DexAem_mask[Dex_mask]


# and call the `bva_sigma_E` function:

# In[14]:


n = 7
E_sub_std = bva_sigma_E(n, bursts_d, DexAem_mask_d)


# Finally, we make a KDE plot of the 2D distribution `E_sub_std` versus the burst FRET efficiency:

# In[15]:


plt.figure(figsize=(6,6))
x = np.arange(0,1.01,0.01)
y = np.sqrt((x*(1-x))/n)
plt.plot(x,y, lw=3, color='red')
im = sns.kdeplot(ds_FRET.E[0], np.asfarray(E_sub_std), shade=True, cmap='viridis', shade_lowest=False)
plt.xlim(0,1)
plt.ylim(0,0.4)
plt.xlabel('E', fontsize=14)
plt.ylabel(r'$s_E$', fontsize=24);


# ---
# **Executed:** Tue Jul 11 21:51:51 2017
# 
# **Duration:** 7 seconds.
# 
# **Autogenerated from:** [Example - Burst Variance Analysis.ipynb](out/Example - Burst Variance Analysis.ipynb)