{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [PyBroMo](http://tritemio.github.io/PyBroMo/) 4. Two-state dynamics - Static smFRET simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<small><i>\n",
    "This notebook is part of <a href=\"http://tritemio.github.io/PyBroMo\" target=\"_blank\">PyBroMo</a> a \n",
    "python-based single-molecule Brownian motion diffusion simulator \n",
    "that simulates confocal smFRET\n",
    "experiments.\n",
    "</i></small>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## *Overview*\n",
    "\n",
    "*In this notebook we generate single-polulation static smFRET data files \n",
    "from the same diffusion trajectories. These files are needed by the [next notebook](PyBroMo - 5. Two-state dynamics - Dynamic smFRET simulation)\n",
    "to generate smFRET data with 2-state dynamics.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading the software\n",
    "\n",
    "Import all the relevant libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from pathlib import Path\n",
    "import numpy as np\n",
    "import tables\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import pybromo as pbm\n",
    "import phconvert as phc\n",
    "print('Numpy version:', np.__version__)\n",
    "print('PyTables version:', tables.__version__)\n",
    "print('PyBroMo version:', pbm.__version__)\n",
    "print('phconvert version:', phc.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIM_PATH = 'data/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define populations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We assume a $\\gamma = 0.7$ and two populations, one with $E_{PR}=0.75$ and the other $E_{PR}=0.4$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Epr1 = 0.75\n",
    "Epr2 = 0.4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The corrected $E$ for the two populations are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gamma = 0.7\n",
    "E1 = Epr1 /(Epr1 * (1 - gamma) + gamma)\n",
    "E2 = Epr2 /(Epr2 * (1 - gamma) + gamma)\n",
    "E1, E2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulation takes the uncorrected $E_{PR}$ as input. \n",
    "\n",
    "We want to simulate a second population that has measured brightness scaling as if the difference was\n",
    "only due to the $\\gamma$ factor. Using the definitions $\\Lambda$ and $\\Lambda_\\gamma$ from \n",
    "[(Ingargiola 2017)](https://doi.org/10.1101/156182),\n",
    "we can use the relation:\n",
    "\n",
    "$$\\frac{E}{E_{PR}} = \\frac{\\Lambda}{\\Lambda_\\gamma}$$\n",
    "\n",
    "Solving for $\\Lambda_\\gamma$ or $\\Lambda$ we get:\n",
    "\n",
    "$$ \\Lambda_\\gamma  = \\Lambda\\frac{E_{PR}}{E}$$\n",
    "\n",
    "$${\\Lambda} = {\\Lambda_\\gamma}\\frac{E}{E_{PR}} $$\n",
    "\n",
    "Since $\\Lambda_\\gamma$ is gamma corrected does not depend on $E$. We can compute it from \n",
    "the parameters of population 1 and then use it for finding $\\Lambda$ for population 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Λ1 = 200e3  # kcps, the detected (i.e. uncorrected) peak emission rate for population 1\n",
    "Λγ = Λ1 * Epr1 / E1\n",
    "Λ2 = Λγ * E2 / Epr2\n",
    "Λ2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Λ2 = np.round(Λ2, -3)\n",
    "Λ2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create smFRET data-files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a file for storing timestamps\n",
    "\n",
    "Here we load a diffusion simulation opening a file to save\n",
    "timstamps in *write* mode. Use `'a'` (i.e. *append*) to keep \n",
    "previously simulated timestamps for the given diffusion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = pbm.ParticlesSimulation.from_datafile('0eb9', mode='a', path=SIM_PATH)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S.particles.diffusion_coeff_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#S = pbm.ParticlesSimulation.from_datafile('44dc', mode='a', path=SIM_PATH)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate timestamps of smFRET\n",
    "\n",
    "We want to simulate two separate smFRET files representing two static populations.\n",
    "\n",
    "We start definint the simulation parameters for population 1 with the following syntax:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params1 = dict(\n",
    "    em_rates = (Λ1,),       # Peak emission rates (cps) for each population (D+A)\n",
    "    E_values = (Epr1,),     # FRET efficiency for each population\n",
    "    num_particles = (35,),  # Number of particles in each population\n",
    "    bg_rate_d = 900,        # Poisson background rate (cps) Donor channel\n",
    "    bg_rate_a = 600,        # Poisson background rate (cps) Acceptor channel\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now define population 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params2 = dict(\n",
    "    em_rates = (Λ2,),       # Peak emission rates (cps) for each population (D+A)\n",
    "    E_values = (Epr2,),     # FRET efficiency for each population\n",
    "    num_particles = (35,),  # Number of particles in each population\n",
    "    bg_rate_d = 900,        # Poisson background rate (cps) Donor channel\n",
    "    bg_rate_a = 600,        # Poisson background rate (cps) Acceptor channel\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we also define a static mixture of the two populations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params_mix = dict(\n",
    "    em_rates = (Λ1, Λ2),       # Peak emission rates (cps) for each population (D+A)\n",
    "    E_values = (Epr1, Epr2),   # FRET efficiency for each population\n",
    "    num_particles = (20, 15),  # Number of particles in each population\n",
    "    bg_rate_d = 900,           # Poisson background rate (cps) Donor channel\n",
    "    bg_rate_a = 600,           # Poisson background rate (cps) Acceptor channel\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate static population 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Population 1**: Create the object that will run the simulation and print a summary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mix_sim = pbm.TimestapSimulation(S, **params1)\n",
    "mix_sim.summarize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the simualtion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rs = np.random.RandomState(1234)\n",
    "mix_sim.run(rs=rs, overwrite=False, skip_existing=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save simulation to a smFRET [Photon-HDF5](http://photon-hdf5.org) file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mix_sim.save_photon_hdf5(identity=dict(author='Antonino Ingargiola', \n",
    "                                       author_affiliation='UCLA'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate static population 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Population 2**: Create the object that will run the simulation and print a summary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mix_sim = pbm.TimestapSimulation(S, **params2)\n",
    "mix_sim.summarize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rs = np.random.RandomState(1234)\n",
    "mix_sim.run(rs=rs, overwrite=False, skip_existing=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mix_sim.save_photon_hdf5(identity=dict(author='Antonino Ingargiola', \n",
    "                                       author_affiliation='UCLA'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate static mixture"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Static mixture**: Create the object that will run the simulation and print a summary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mix_sim = pbm.TimestapSimulation(S, **params_mix)\n",
    "mix_sim.summarize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rs = np.random.RandomState(1234)\n",
    "mix_sim.run(rs=rs, overwrite=False, skip_existing=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mix_sim.save_photon_hdf5(identity=dict(author='Antonino Ingargiola', \n",
    "                                       author_affiliation='UCLA'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!rsync -av --exclude 'pybromo_*.hdf5' /mnt/archive/Antonio/pybromo /mnt/wAntonio/dd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Burst analysis\n",
    "\n",
    "The generated Photon-HDF5 files can be analyzed by any smFRET burst\n",
    "analysis program. Here we show an example using the opensource\n",
    "[FRETBursts](https://github.com/tritemio/FRETBursts/) program:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fretbursts as fb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filepath = list(Path(SIM_PATH).glob('smFRET_*'))\n",
    "filepath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = fb.loader.photon_hdf5(str(filepath[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d.A_em"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fb.dplot(d, fb.timetrace);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d.calc_bg(fun=fb.bg.exp_fit, tail_min_us='auto', F_bg=1.7, time_s=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fb.dplot(d, fb.timetrace_bg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d.burst_search(F=7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d.num_bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = d.select_bursts(fb.select_bursts.size, th1=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.num_bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with plt.rc_context({#'font.size': 10, \n",
    "                     #'savefig.dpi': 200, \n",
    "                     'figure.dpi': 150}):\n",
    "    for i in range(3):\n",
    "        fig, ax = plt.subplots(figsize=(100, 3))\n",
    "        fb.dplot(d, fb.timetrace, binwidth=0.5e-3, tmin=i*10, tmax=(i+1)*10, bursts=True, \n",
    "                 plot_style=dict(lw=1), ax=ax);\n",
    "        ax.set_xlim(i*10, (i+1)*10);\n",
    "        display(fig)\n",
    "        plt.close(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fb.dplot(ds, fb.hist_fret, pdf=False)\n",
    "plt.axvline(0.4, color='k', ls='--');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fb.bext.burst_data(ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fb.dplot(d, fb.hist_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fb.dplot(d, fb.hist_width)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.5.2"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}