.. currentmodule:: fretbursts.burstlib Burst Search in FRETBursts ========================== This section describes details and conventions used to implement burst search in FRETBursts. For a more general explanation of burst search concepts see `(Ingargiola PLOS ONE 2016) `__. For usage examples see the `μs-ALEX notebook `__. An analysis of implementation performances of the :ref:`low-level burst search ` can be found in this blog post: `Optimizing burst search in python `__. Defining the rate estimator --------------------------- Before describing FRETBursts implementation let me introduce an expression for computing rates of random events that will be used later on. A general expression, used by FRETBursts (since version 0.5.6), for estimating the rate using *m* consecutive timestamps is: .. math:: \hat{\lambda} = \frac{m - 1 - c}{t_{i + m - 1} - t_{i}} :label: ratedef where :math:`c` is a parameter that can be passed to all FRETBursts functions that deal with photon rates. Note that :math:`m` is the number of photons and :math:`m - 1` is the number of inter-photon delays. For example, using :math:`c=1`, yields an unbiased estimator of the rate for events generated by a stationary Poisson process. See `this notebook `__ for a discussion of the different estimator properties as a function of :math:`c`. In practice, the choice of :math:`c` is just a convention and it is provided for flexibility and to match results of other software that may use a different definition. In FRETBursts version 0.5.5 or earlier, there is no `c` parameter and the rate is always computed as :math:`\hat{\lambda} = m / (t_{i + m - 1} - t_{i})` (equivalent to :math:`c=-1`). Conventions in burst search --------------------------- Burst search is mainly performed calling the method :meth:`Data.burst_search`. The AND-gate burst search function (:meth:`fretbursts.burstlib_ext.burst_search_and`) calls :meth:`Data.burst_search` under the hood, so all the considerations below are also valid for the AND-gate version. With :meth:`Data.burst_search`, you can perform burst search by setting a "rate threshold" *F* times larger than the background rate (argument `F`), or you can just set a single fixed rate for the full measurement (argument `min_rate_cps`). In both cases the real burst search is performed by the low-level function :func:`.phtools.burstsearch.bsearch_py`, which takes as input parameters *m* and *T*. This function finds bursts when a group of *m* consecutive photons lies within a time window *T*. You can find an analysis of the algorithm implementation and performance considerations in this `blog post `__. When using the `F` argument, FRETBursts will choose the appropriate `T` for each background period in order to obtain a "rate threshold" *F times larger than background rate*. In this case, FRETBursts uses the following expression to compute *T* (derived from :eq:`ratedef`): .. math:: T(t) = \frac{m - 1 - c}{F \cdot \hat{\lambda}_{bg}(t)} where :math:`\hat{\lambda}_{bg}(t)` is the estimated background rate as a function of time (:math:`t`). Conversely, when directly fixing a rate with the argument `min_rate_cps` (:math:`\lambda_{th}`), FRETBursts computes *T* using the expression: .. math:: T = \frac{m - 1 - c}{\lambda_{th}} The parameter :math:`c` can be specified when performing burst search. When not specified, the default value of :math:`c=-1` is used. This choice preserves backward compatibility with results obtained with FRETBursts 0.5.5 or earlier. .. _burstsearch_core: The Core Algorithm ------------------ The different types of burst search described in the previous sections are implemented calling the same low-level burst search function which implements the core "sliding window" algorithm. Here we explain in details this core algorithm. The low-level burst search takes as an input the array of (monotonically increasing) photon timestamps, as well as two other arguments *m* (the number of timestamps) and *T* (the time window duration). Starting from the the first element of the array, we consider all the m-tuple of timestamps [0..m-1], [1..m], etc. **Point 1.** For each m-tuple if the timestamps are contained in a time window smaller or equal to *T* we mark a burst start at the position of the first timestamp in the current m-tuple. Otherwise we take the next m-tuple and repeat the check. Once a burst starts, we keep "sliding" the m-tuple one timestamp a time. If the current m-tuple is still contained in a window of duration *T* the burst continues. When the current m-tuple is contained in a window larger than *T* the burst ends. When this happens, the last timestamp in a burst is the (m-1)-th timestamp of current m-tuple (i.e. the last timestamps of the **previous** m-tuple which was still contained in a window *T*). After the burst ends, we continue as in point 1 checking the next m-tuple. This is shifted by only one timestamp (i.e. there is no jump when the burst ends). At this point, it can happen that the current m-tuple is contained in *T* and a new burst starts right away. In this situation the new bursts will have m-2 timestamps overlapping with the previous one. At the end of the timestamp array, if a burst is currently started we end it by marking the last timestamp as burst stop. The set of bursts obtained in this way has the minimum-rate property, i.e. all the m-tuple of consecutive timestamps in any burst are guaranteed to be contained in a windows *T* or smaller. Conversely, a few bursts will overlap and thus share some timestamps. If the user wants to avoid overlapping bursts a burst fusion steps must be applied as described in next section. Note, however, that after fusing overlapping bursts at least one m-tuple inside each fused burst will not have the minimum-rate property, i.e. the m-tuple is contained in a window larger than *T*. .. py:currentmodule:: fretbursts The previous function is implemented in :func:`phtools.burstsearch.bsearch_py` (pure python version) and in :func:`phtools.burstsearch_c.bsearch_c` (optimized cython version). Several tests make sure that the two functions return numerically identical results. An analysis of performance of of different implementations can be found in this blog post: `Optimizing burst search in python `__. Burst Fusion ------------ Burst fusion is an operation which fuses consecutive bursts if the start of the second bursts minus the end of the first burst (called burst separation) is <= of a fusion time :math:`t_f`. When bursts are overlapping (see previous section) the burst separation is negative. Therefore, to avoid overlapping bursts, we need to apply fusion with separation of 0. Note that with this condition, if a bursts ends on a timestamp which is the start of the next burst (i.e. 1 overlapping photon) the two bursts will be fused. Conversely if one burst ends and the next burst starts one photon later (0 overlapping photons) the two bursts will be kept separated. In the latter case, there will be no timestamp between the end of the previous burst and the start of the next one. .. py:currentmodule:: fretbursts.burstlib To perform burst fusion use the method :meth:`Data.fuse_bursts`. Low-level burst search functions -------------------------------- .. automodule:: fretbursts.phtools.burstsearch :members: