{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial about analyzing blink statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SMLM depends critically on the fluorescence intermittency or, in other words, the blinking of fluorescence dyes. To characterize blinking properties you can compute on- and off-periods from clustered localizations assuming that they originate from the same fluorophore."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as stats\n",
    "\n",
    "import locan as lc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc.show_versions(system=False, dependencies=False, verbose=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Synthetic data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use synthetic data that represents localizations from a single fluorophore being normal distributed in space and emitting at a constant intensity. We assume that the on- and off-times in units of frames are distributed like a geometric distribution with a mean on_period `mean_on` and a mean off_period `mean_off`. Typically a geometric distribution is parameterized by a variable `p` with `p = 1 / mean`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(seed=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samples = 10_000\n",
    "mean_on = 5\n",
    "mean_off = 20\n",
    "\n",
    "on_periods = stats.geom.rvs(p=1/mean_on, size=n_samples, random_state=rng)\n",
    "off_periods = stats.geom.rvs(p=1/mean_off, size=n_samples, random_state=rng)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On- and off-times are converted in a series of frame numbers at which a localization was detected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def periods_to_frames(on_periods, off_periods):\n",
    "    \"\"\"\n",
    "    Convert on- and off-periods into a series of increasing frame values.\n",
    "    \"\"\"\n",
    "    on_frames = np.arange(np.sum(on_periods))\n",
    "    cumsums = np.r_[0, np.cumsum(off_periods)[:-1]]\n",
    "    add_on = np.repeat(cumsums, on_periods)\n",
    "    frames = on_frames + add_on\n",
    "    return frames[:len(on_periods)]\n",
    "\n",
    "frames = periods_to_frames(on_periods, off_periods)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "offspring = [rng.normal(loc=0, scale=10, size=(n_samples, 2))]\n",
    "locdata = lc.simulate_cluster(centers=[(50, 50)], region=[(0, 100), (0, 100)], offspring=offspring, clip=False, shuffle=False, seed=rng)\n",
    "locdata.dataframe['intensity'] = 1\n",
    "locdata.dataframe['frame'] = frames\n",
    "\n",
    "locdata = lc.LocData.from_dataframe(dataframe=locdata.data)\n",
    "\n",
    "print('Data head:')\n",
    "print(locdata.data.head(), '\\n')\n",
    "print('Summary:')\n",
    "locdata.print_summary()\n",
    "print('Properties:')\n",
    "print(locdata.properties)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc.render_2d(locdata, bin_size=5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Blinking statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To determine on- and off-times for the observed blink events use the analysis class `BlinkStatistics`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bs = lc.BlinkStatistics(memory=0, remove_heading_off_periods=False).compute(locdata)\n",
    "bs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bs.results.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When plotting the histogram an exponential distribution is fitted by default."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bs.hist(data_identifier='on_periods');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bs.hist(data_identifier='off_periods');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bs.distribution_statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fit results provide `loc` and `scale` parameter (see `scipy.stats` documentation). For `loc = 0`, `scale` describes the mean of the distribution.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bs.distribution_statistics['on_periods'].parameter_dict()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bs.distribution_statistics['off_periods'].parameter_dict()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to the default setting for the scaling parameter `loc` the mean on_period is `on_periods_scale + on_periods_loc` in agreement with our input value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geometric distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare this with a geometric distribution that is estimated from the observed mean on_period `on_periods_mean`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "on_periods_mean = bs.results['on_periods'].mean()\n",
    "on_periods_mean.round(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "off_periods_mean = bs.results['off_periods'].mean()\n",
    "off_periods_mean.round(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test result\n",
    "x = np.arange(stats.geom.ppf(0.01, 1/on_periods_mean), stats.geom.ppf(0.9999, 1/on_periods_mean))\n",
    "y = stats.geom.pmf(x, 1/on_periods_mean)\n",
    "fig, ax = plt.subplots()\n",
    "bs.hist(data_identifier='on_periods', fit=False, label='data')\n",
    "bs.distribution_statistics['on_periods'].plot(label='exponential')\n",
    "ax.plot(x, y, '-go', label='geometric')\n",
    "ax.set_yscale('log')\n",
    "ax.legend(loc='best')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}