notebooks/PyBroMo - 1. Simulate 3D trajectories - single core.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [PyBroMo](http://opensmfs.github.io/PyBroMo/) - 1. Simulate 3D trajectories - single core"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<small><i>\n",
    "This notebook is part of <a href=\"http://opensmfs.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 show how to perform a 3-D trajectories simulation of a set of freely diffusing molecules. The simulation computes (and saves!) 3-D trajectories and emission rates due to a confocal excitation PSF for each single molecule. Depending on the simulation length, the required disk space can be significant (~ 750MB per minute of simulated diffusion).*\n",
    "\n",
    "*For more info see [PyBroMo Homepage](http://opensmfs.github.io/PyBroMo/)*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Together with a few standard python libraries we import **PyBroMo** using the short name `pbm`. \n",
    "All **PyBroMo** functions will be available as `pbm.`*something*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\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",
    "print('Numpy version:', np.__version__)\n",
    "print('PyTables version:', tables.__version__)\n",
    "print('PyBroMo version:', pbm.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we define the simulation parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the random state\n",
    "rs = np.random.RandomState(seed=1)\n",
    "print('Initial random state:', pbm.hash_(rs.get_state()))\n",
    "\n",
    "# Simulation time step (seconds)\n",
    "t_step = 0.5e-6\n",
    "\n",
    "# Time duration of the simulation (seconds)\n",
    "t_max = 1\n",
    "\n",
    "# Diffusion coefficient\n",
    "Du = 12.0            # um^2 / s\n",
    "D1 = Du*(1e-6)**2    # m^2 / s\n",
    "D2 = D1/2\n",
    "\n",
    "# Simulation box definition\n",
    "box = pbm.Box(x1=-4.e-6, x2=4.e-6, y1=-4.e-6, y2=4.e-6, z1=-6e-6, z2=6e-6)\n",
    "\n",
    "# Particles definition\n",
    "P = pbm.Particles.from_specs(\n",
    "    num_particles=(3, 1),\n",
    "    D=(D1, D2),\n",
    "    box=box, rs=rs)\n",
    "\n",
    "# PSF definition\n",
    "psf = pbm.NumericPSF()\n",
    "\n",
    "# Particle simulation definition\n",
    "S = pbm.ParticlesSimulation(t_step=t_step, t_max=t_max, \n",
    "                            particles=P, box=box, psf=psf)\n",
    "\n",
    "print('Current random state:', pbm.hash_(rs.get_state()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most important line is the last line which creates an object `S` \n",
    "that contains all the simulation parameters (it also contains methods to run \n",
    "the simulation). You can print `S` and check the current parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other useful simulation info:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S.compact_name()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S.particles.diffusion_coeff_counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arrays sizes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S.print_sizes()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **NOTE:** This is the maximum in-memory array size when using a single chunk. \n",
    "> In the following, we simulate the diffusion in smaller time windows (chunks), \n",
    "> thus requiring only a few tens MB of RAM, regardless of the simulated duration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Brownian motion simulation\n",
    "\n",
    "In the brownian motion simulation we keep using the same random state object `rs`. \n",
    "Initial and final state are saved so the same simulation can be reproduced. \n",
    "See [PyBroMo - A1. Reference - Data format and internals.ipynb](PyBroMo - A1. Reference - Data format and internals.ipynb) \n",
    "for more info on the random state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Current random state:', pbm.hash_(rs.get_state()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S.simulate_diffusion(total_emission=False, save_pos=True, verbose=True, \n",
    "                     rs=rs, chunksize=2**19, chunkslice='times')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Current random state:', pbm.hash_(rs.get_state()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The normalized emission rate (peaks at 1) for each particle is stored \n",
    "in a 2D pytable array and accessible through the `emission` attribute:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Simulation file size: %d MB' % (S.store.h5file.get_filesize()/1024./1024.))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S.store.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load trajectories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S.hash()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = pbm.ParticlesSimulation.from_datafile(S.hash()[:4])  # Read-only by default"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting the emission"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This are basic debug plots. For more advanged interactive \n",
    "exploration of trajectory and emission arrays see the next notebook:\n",
    "\n",
    "- [PyBroMo - GUI Trajectory explorer](PyBroMo - GUI Trajectory explorer.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_emission(S, s=0, size=2e6, slice_=None, em_th=0.01, save=False, figsize=(9, 4.5)):\n",
    "    if slice_ is None:\n",
    "        slice_ = (s*size, (s+1)*size)\n",
    "    slice_ = slice(*slice_)\n",
    "    em = S.emission[:, slice_]\n",
    "    dec = 1 if slice_.step is None else slice_.step\n",
    "    t_step = S.t_step*dec\n",
    "    t = np.arange(em.shape[1])*(t_step*1e3)\n",
    "    fig, ax = plt.subplots(figsize=figsize)\n",
    "    for ip, em_ip in enumerate(em):\n",
    "        if em_ip.max() < em_th: continue\n",
    "        plt.plot(t, em_ip, label='P%d' % ip)\n",
    "    ax.set_xlabel('Time (ms)')\n",
    "    \n",
    "    rs_hash = pbm.hash_(S.traj_group._v_attrs['init_random_state'])[:3]\n",
    "    ax.set_title('%ds ID-EID: %d-%d, sim rs = %s, part rs = %s' %\\\n",
    "              (s, S.ID, S.EID, rs_hash, S.particles.rs_hash[:3]))\n",
    "    ax.legend(bbox_to_anchor=(1.03, 1), loc=2, borderaxespad=0.)\n",
    "    if save:\n",
    "        plt.savefig('em %ds ID-EID %d-%d, rs=%s' %\\\n",
    "                (s, S.ID, S.EID, rs_hash), \n",
    "                dpi=200, bbox_inches='tight')\n",
    "    #plt.close(fig)\n",
    "    #display(fig)\n",
    "    #fig.clear()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_emission(S, slice_=(0, 2e6, 10), em_th=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_tracks(S, slice_=None, particles=None):\n",
    "    if slice_ is None:\n",
    "        slice_ = (0, 100e3, 100)\n",
    "    duration = (slice_[1] - slice_[0])*S.t_step\n",
    "    slice_ = slice(*slice_)\n",
    "    \n",
    "    if particles is None:\n",
    "        particles = range(S.num_particles)\n",
    "    \n",
    "    fig, AX = plt.subplots(1, 2, figsize=(11, 5), sharey=True)\n",
    "    plt.subplots_adjust(left=0.05, right=0.93, top=0.95, bottom=0.09,\n",
    "                        wspace=0.05)\n",
    "    plt.suptitle(\"Total: %.1f s, Visualized: %.2f ms\" % (\n",
    "                 S.t_step*S.n_samples, duration*1e3))\n",
    "\n",
    "    for ip in particles:\n",
    "        x, y, z = S.position[ip, :, slice_]\n",
    "        x0, y0, z0 = S.particles[ip].r0\n",
    "        plot_kwargs = dict(ls='', marker='o', mew=0, ms=2, alpha=0.5, \n",
    "                           label='P%d' % ip)\n",
    "        l, = AX[0].plot(x*1e6, y*1e6, **plot_kwargs)\n",
    "        AX[1].plot(z*1e6, y*1e6, color=l.get_color(), **plot_kwargs)\n",
    "        #AX[1].plot([x0*1e6], [y0*1e6], 'o', color=l.get_color())\n",
    "        #AX[0].plot([x0*1e6], [z0*1e6], 'o', color=l.get_color())\n",
    "\n",
    "    AX[0].set_xlabel(\"x (um)\")\n",
    "    AX[0].set_ylabel(\"y (um)\")\n",
    "    AX[1].set_xlabel(\"z (um)\")\n",
    "\n",
    "    sig = np.array([0.2, 0.2, 0.6])*1e-6\n",
    "    ## Draw an outline of the PSF\n",
    "    a = np.arange(360)/360.*2*np.pi\n",
    "    rx, ry, rz = (sig)  # draw radius at 3 sigma\n",
    "    AX[0].plot((rx*np.cos(a))*1e6, (ry*np.sin(a))*1e6, lw=2, color='k')\n",
    "    AX[1].plot((rz*np.cos(a))*1e6, (ry*np.sin(a))*1e6, lw=2, color='k')\n",
    "    \n",
    "    AX[0].set_xlim(-4, 4)\n",
    "    AX[0].set_ylim(-4, 4)\n",
    "    AX[1].set_xlim(-4, 4)\n",
    "    \n",
    "    if len(particles) <= 20:\n",
    "        AX[1].legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_tracks(S) #particles=[2, 5, 7, 22, 30])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S.store.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Low-level timestamp simulation\n",
    "\n",
    "Here we simulate some timestamps array one by one. To generate smFRET data in one step (donor + acceptor, single or multiple populations and create [Photon-HDF5](http://photon-hdf5.org) files see the next notebook:\n",
    "\n",
    "- [PyBroMo - 2. Generate smFRET data, including mixtures](PyBroMo - 2. Generate smFRET data, including mixtures.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "code = S.compact_name()[:4]\n",
    "code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = pbm.ParticlesSimulation.from_datafile(code, mode='w')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulate timestamps treating all particles as a single population with same max \n",
    "brightness (particles may still have different diffusion coefficient though):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts_params1 = dict(\n",
    "    max_rates=(200e3,), \n",
    "    populations=(slice(0, 3),),\n",
    "    bg_rate=1e3,\n",
    "    )\n",
    "\n",
    "S.simulate_timestamps_mix(save_pos=True, overwrite=True, **ts_params1)\n",
    "S.simulate_timestamps_mix(save_pos=False, overwrite=True, **ts_params1)\n",
    "S.simulate_timestamps_mix(save_pos=False, overwrite=True, **ts_params1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulate timestamps giving different brightness to each populations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts_params2 = dict(\n",
    "    max_rates=(250e3, 180e3), \n",
    "    populations=(slice(0,3), slice(3, 3+1)),\n",
    "    bg_rate=1.2e3\n",
    ")\n",
    "\n",
    "S.simulate_timestamps_mix(save_pos=True, overwrite=True, **ts_params2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the timestamps and particles (pytables) arrays, we need to find the name first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S.timestamp_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Picking the last of the names above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts, part, pos = S.get_timestamp_data(S.timestamp_names[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts, part, pos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the data in memory as a numpy array use `[:]`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts[:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can find the name of a timestamps array with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S.timestamps_match_mix(**ts_params2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S.timestamps_match_mix(**ts_params2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts, part, pos = S.get_timestamp_data(S.timestamp_names[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos, ts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3.6 (py36)",
   "language": "python",
   "name": "py36"
  },
  "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.6.7"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {
    "height": "591.867px",
    "left": "0px",
    "right": "959px",
    "top": "111.133px",
    "width": "212px"
   },
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}