from fretbursts import *
sns = init_notebook()

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()

"""
Explanation: Working with timestamps and bursts
This notebook is part of a tutorial series for the 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.

Load data
We start by loading the data, computing background and performing a standard burst search:
End of explanation
"""


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

"""
Explanation: Getting the timestamps
To get the timestamps arrays for a given photon stream we use Data.get_ph_times. Here a few example:
End of explanation
"""


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

"""
Explanation: 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:
End of explanation
"""


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

"""
Explanation: 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:
End of explanation
"""


mask_d.sum()

"""
Explanation: Masks can be used to count photons in one stream:
End of explanation
"""


ph[mask_d]

"""
Explanation: and to obtain the timestamps for one stream:
End of explanation
"""


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

"""
Explanation: 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:
End of explanation
"""


bursts

"""
Explanation: 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):
End of explanation
"""


firstburst = bursts[0]
firstburst

"""
Explanation: Indexing bursts we can access a single burst:
End of explanation
"""


bursts.istart

firstburst.istart

"""
Explanation: 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:
End of explanation
"""


ph[firstburst.istart], firstburst.start

"""
Explanation: Note that ph[firstburst.istart] is equal to firstburst.start:
End of explanation
"""


ph[firstburst.istop], firstburst.stop

"""
Explanation: The same holds for stop:
End of explanation
"""


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

"""
Explanation: 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

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:
End of explanation
"""


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

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

nd

"""
Explanation: 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):
End of explanation
"""


from fretbursts.phtools.burstsearch import Burst, Bursts

times = d.ph_times_m[0]  # timestamps array

"""
Explanation: 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.
End of explanation
"""


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

"""
Explanation: We start fusing bursts with separation <= 0 milliseconds,
to avoid having overlapping bursts:
End of explanation
"""


time_bin = 0.5e-3  # 0.5 ms

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)

"""
Explanation: Now we can slice each burst using a constant time bin:
End of explanation
"""


len(sub_bursts_list)

ds_fused.num_bursts

"""
Explanation: The list sub_bursts_list has one set of sub-bursts per each original burst:
End of explanation
"""


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

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

"""
Explanation: Each set of sub-bursts is a usual Bursts object:
End of explanation
"""


bursts = sub_bursts_list[0]
bursts

"""
Explanation: 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).
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:
End of explanation
"""


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

"""
Explanation: <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)):
End of explanation
"""


from fretbursts.phtools.burstsearch import count_ph_in_bursts

counts_dd = count_ph_in_bursts(bursts, mask_dd)
counts_dd

"""
Explanation: Next, we use the function counts_ph_in_bursts:
End of explanation
"""


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

"""
Explanation: 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:
End of explanation
"""


counts_ad / (counts_dd + counts_ad)

"""
Explanation: With these values, we can compute, for example, the uncorrected
Proximity Ratio (PR):
End of explanation
"""


