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

# # Working with timestamps and bursts
# 
# *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 show how to access different streams of timestamps,
# > burst data (i.e. start and stop times and indexes, counts, etc...).
# > These operations are useful for users wanting to access or process bursts data 
# > and timestamps in custom ways.
# > For a complete tutorial on burst analysis see 
# > [FRETBursts - us-ALEX smFRET burst analysis](FRETBursts - us-ALEX smFRET burst analysis.ipynb).

# # Load data
# 
# We start by loading the data, computing background and performing a standard burst search:

# In[1]:


from fretbursts import *
sns = init_notebook()


# In[2]:


filename = "./data/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5"
d = loader.photon_hdf5(filename)
loader.alex_apply_period(d)
d.calc_bg(bg.exp_fit, time_s=30, tail_min_us='auto', F_bg=1.7)
d.burst_search()


# # Getting the timestamps
# 
# To get the timestamps arrays for a given [photon stream](http://fretbursts.readthedocs.org/en/latest/ph_sel.html) we use [Data.get_ph_times](http://fretbursts.readthedocs.org/en/latest/data_class.html?highlight=get_ph_times#fretbursts.burstlib.Data.get_ph_times). Here a few example:

# In[3]:


ph = d.get_ph_times()                              # all the recorded photons
ph_dd = d.get_ph_times(ph_sel=Ph_sel(Dex='Dem'))   # donor excitation, donor emission
ph_d = d.get_ph_times(ph_sel=Ph_sel(Dex='DAem'))   # donor excitation, donor+acceptor emission
ph_aa = d.get_ph_times(ph_sel=Ph_sel(Aex='Aem'))   # acceptor excitation, acceptor emission


# This are streams of all timestamps (both inside and outside the bursts).
# Similarly, we can get "masks" of photons for each photon stream using 
# [Data.get_ph_mask](http://fretbursts.readthedocs.org/en/latest/data_class.html?highlight=get_ph_mask#fretbursts.burstlib.Data.get_ph_mask):

# In[4]:


mask_dd = d.get_ph_mask(ph_sel=Ph_sel(Dex='Dem'))   # donor excitation, donor emission
mask_d = d.get_ph_mask(ph_sel=Ph_sel(Dex='DAem'))   # donor excitation, donor+acceptor emission
mask_aa = d.get_ph_mask(ph_sel=Ph_sel(Aex='Aem'))   # acceptor excitation, acceptor emission


# Masks are arrays of booleans (True or False values) which are True
# when the corresponding photon is in the stream. Note that all masks
# have same number of elements as the all-photons timestamps array:

# In[5]:


ph.size, mask_dd.size, mask_d.size, mask_aa.size


# Masks can be used to count photons in one stream:

# In[6]:


mask_d.sum()


# and to obtain the timestamps for one stream:

# In[7]:


ph[mask_d]


# Note that the arrays `ph[mask_d]` and `ph_d` are equal. This is an important point to understand.

# # Burst data
# 
# There are several fields containing burst data:
# 
# **Start-stop**:
# 
# - `Data.mburst`: start-stop information for each burst
# 
# **Counts**:
# - `Data.nd`: donor detector counts during donor excitation
# - `Data.na`: acceptor detector counts during donor excitation
# - `Data.naa`: acceptor detector counts during acceptor excitation (ALEX only)
# - `Data.nda`: donor detector counts during acceptor excitation
# 
# **FRET**:
# - `Data.E`: Proximity Ratio (when uncorrected) or FRET efficiency (when corrected)
# - `Data.S`: "Stoichiometry" (ALEX only)
# 
# All previous fields are **lists** containing one element per excitation spot.
# In single-spot data, these lists have only one element which is accessed
# using the `[0]` notation:

# In[8]:


bursts = d.mburst[0]
nd = d.nd[0]
na = d.na[0]
naa = d.naa[0]
E = d.E[0]
S = d.S[0]


# All previous variables are numpy arrays, except for `bursts` which is
# a `Bursts` object (see next section).

# ## Burst start and stop
# 
# The start-stop burst data is in `bursts` (a variable of type [Bursts](), plural):

# In[9]:


bursts


# Indexing `bursts` we can access a single burst:

# In[10]:


firstburst = bursts[0]
firstburst


# The first two "columns" (both in `bursts` or `firstburst`) are the index of 
# **first** and **last** timestamps (relative to the all-photons timestamps).
# The last two columns (`start` and `stop`) are the actual times of burst 
# start and stop. To access any of these fields we type:

# In[11]:


bursts.istart


# In[12]:


firstburst.istart


# Note that `ph[firstburst.istart]` is equal to `firstburst.start`:

# In[13]:


ph[firstburst.istart], firstburst.start


# The same holds for stop:

# In[14]:


ph[firstburst.istop], firstburst.stop


# Note that `bursts` is a `Bursts` object (plural, a bursts-set) 
# and `firstburst` is a `Burst` object (singular, only one burst).
# You can find more info on these objects in the documentation:
# 
# - [Low-level burst search functions](http://fretbursts.readthedocs.org/en/latest/burstsearch.html)

# ## Burst photon-counts
# 
# The variables `nd`, `na`, `naa` contains the number of photon in different photon streams.
# By default these values are background corrected and, if the correction coefficients 
# are specified, are also corrected for leakage, direct excitation and gamma factor.
# 
# To get the raw counts before correction we can redo the burst search as follow:

# In[15]:


d.burst_search(computefret=False)
d.calc_fret(count_ph=True, corrections=False)


# Note that if you select bursts, you also need to use `computefret=False`
# to avoid recomputing E and S (which by default applies the corrections):

# In[16]:


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


# In[17]:


nd = ds.nd[0]      # Donor-detector counts during donor excitation
na = ds.na[0]      # Acceptor-detector counts during donor excitation
naa = ds.naa[0]    # Acceptor-detector counts during acceptor excitation
E = ds.E[0]        # FRET efficiency or Proximity Ratio
S = ds.S[0]        # Stoichiometry, as defined in µs-ALEX experiments


# In[18]:


nd


# Note that the burst counts are integer values, confirming that the background
# correction was not applied.

# # Slice bursts in time bins
# 
# Here we slice each burst in fixed time bins.

# In[19]:


from fretbursts.phtools.burstsearch import Burst, Bursts


# In[20]:


times = d.ph_times_m[0]  # timestamps array


# We start fusing bursts with separation <= 0 milliseconds,
# to avoid having overlapping bursts:

# In[21]:


ds_fused = ds.fuse_bursts(ms=0)
bursts = ds_fused.mburst[0]
print('\nNumber of bursts:', bursts.num_bursts)


# Now we can slice each burst using a constant time bin:

# In[22]:


time_bin = 0.5e-3  # 0.5 ms


# In[23]:


time_bin_clk = time_bin / ds.clk_p

sub_bursts_list = []
for burst in bursts:
    # Compute binning of current bursts
    bins = np.arange(burst.start, burst.stop + time_bin_clk, time_bin_clk)
    counts, _ = np.histogram(times[burst.istart:burst.istop+1], bins)
    
    # From `counts` in each bin, find start-stop times and indexes (sub-burst).
    # Note that start and stop are the min and max timestamps in the bin,
    # therefore they are not on the bin edges. Also the burst width is not
    # exactly equal to the bin width.
    sub_bursts_l = []
    sub_start = burst.start
    sub_istart = burst.istart
    for count in counts:
        # Let's skip bins with 0 photons
        if count == 0:
            continue
            
        sub_istop = sub_istart + count - 1
        sub_bursts_l.append(Burst(istart=sub_istart, istop=sub_istop,
                                  start=sub_start, stop=times[sub_istop]))
        
        sub_istart += count 
        sub_start = times[sub_istart]
    
    sub_bursts = Bursts.from_list(sub_bursts_l)
    assert sub_bursts.num_bursts > 0
    assert sub_bursts.width.max() < time_bin_clk
    sub_bursts_list.append(sub_bursts)


# The list `sub_bursts_list` has one set of sub-bursts per each original burst:

# In[24]:


len(sub_bursts_list)


# In[25]:


ds_fused.num_bursts


# Each set of sub-bursts is a usual [Bursts](http://fretbursts.readthedocs.org/en/latest/burstsearch.html?highlight=burstsearchlib.bursts#fretbursts.burstsearch.burstsearchlib.Bursts) object:

# In[26]:


print('Sub-bursts from burst 0:')
sub_bursts_list[0]


# In[27]:


iburst = 10
print('Sub-bursts from burst %d:' % iburst)
sub_bursts_list[iburst]


# # Photon counts in custom bursts
# 
# When performing a burst search,
# FRETBursts automatically computes donor and acceptor counts (in both 
# excitation periods). These quantities are available as `Data` attributes:
# `nd`, `na`, `naa` and `nda` (as described in [Burst data](#Burst-data)).
# 
# When a custom bursts-set is created, like in the previous section in which we
# sliced bursts in sub-bursts, we may want to computed the photon counts 
# in the various photon streams. Let consider, as an example, the following `Bursts` object:

# In[28]:


bursts = sub_bursts_list[0]
bursts


# <p class="lead">How do we count the <b>donor</b> and <b>acceptor</b> photons in these bursts?<p>
# 
# First we need to prepare the masks for the various photon streams 
# (as explained [before](#Getting-the-timestamps)):

# In[29]:


mask_dd = d.get_ph_mask(ph_sel=Ph_sel(Dex='Dem'))   # donor excitation, donor emission
mask_ad = d.get_ph_mask(ph_sel=Ph_sel(Dex='Aem'))   # donor excitation, acceptor emission
mask_aa = d.get_ph_mask(ph_sel=Ph_sel(Aex='Aem'))   # acceptor excitation, acceptor emission


# Next, we use the function `counts_ph_in_bursts`:

# In[30]:


from fretbursts.phtools.burstsearch import count_ph_in_bursts


# In[31]:


counts_dd = count_ph_in_bursts(bursts, mask_dd)
counts_dd


# `counts_dd` contains the raw counts in each burst (in `bursts`)
# in the Donor-emission during Donor-acceptor stream. Similarly,
# we can compute counts for the other photon streams:

# In[32]:


counts_ad = count_ph_in_bursts(bursts, mask_ad)
counts_aa = count_ph_in_bursts(bursts, mask_aa)


# With these values, we can compute, for example, the uncorrected
# Proximity Ratio (PR):

# In[33]:


counts_ad / (counts_dd + counts_ad)


# Note that the first value is `nan` (not-a-number) because there
# is a zero in the numerator, i.e. the first bursts has 0 donor counts
# during donor excitation.
# 
# Finally, remember that no correction (background, leakage, etc.) has been 
# applied so far.

# # See also
# 
# This section of the documentation has more info how to use `Bursts` objects:
# 
# - [Low-level burst search functions](http://fretbursts.readthedocs.org/en/latest/burstsearch.html)
# 
# In particular, these 3 methods, allow transforming bursts to refer a new timestamps arrays:
# 
# - [Bursts.recompute_times](http://fretbursts.readthedocs.org/en/latest/burstsearch.html?highlight=recompute_times#fretbursts.burstsearch.burstsearchlib.Bursts.recompute_times)
# - [Bursts.recompute_index_expand](http://fretbursts.readthedocs.org/en/latest/burstsearch.html?highlight=recompute_index_expand#fretbursts.burstsearch.burstsearchlib.Bursts.recompute_index_expand)
# - [Bursts.recompute_index_reduce](http://fretbursts.readthedocs.org/en/latest/burstsearch.html?highlight=recompute_index_reduce#fretbursts.burstsearch.burstsearchlib.Bursts.recompute_index_reduce)
# 
#