{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# [PyBroMo](http://tritemio.github.io/PyBroMo/) - B.1 Disk-single-core - Generate photon timestamps" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "This notebook is part of PyBroMo a \n", "python-based single-molecule Brownian motion diffusion simulator \n", "that simulates confocal smFRET\n", "experiments.\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## *Overview*\n", "\n", "*In this notebook we show how to generated timestamps of emitted-photons from saved diffusion traces*." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "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": [ "# Timestamps simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a memo, let's write some formulas related to the FRET efficiency:\n", "\n", "$$ k = \\frac{F_a}{F_d} \\quad,\\qquad E = \\frac{k}{k+1} \\qquad\\Rightarrow\\qquad k = \\frac{E}{1-E}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "S = pbm.ParticlesSimulation.from_datafile('016', mode='w')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def em_rates_from_E(em_rate_tot, E_values):\n", " E_values = np.asarray(E_values)\n", " em_rates_a = E_values * em_rate_tot\n", " em_rates_d = em_rate_tot - em_rates_a\n", "\n", " k_values = E_values/(1 - E_values)\n", " assert np.allclose((em_rates_a/em_rates_d), k_values)\n", "\n", " em_rates = np.hstack([em_rates_a, em_rates_d])\n", " return em_rates" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "em_rate_tot = 200e3\n", "E_list = np.array([0.01, 0.02, 0.05, 0.1, 0.2, 0.4])\n", "\n", "em_rate_list = em_rates_from_E(em_rate_tot, E_list)\n", "em_rate_list" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Get the random state at the end of the diffusion simulation\n", "saved_rs_state = S.traj_group._v_attrs['last_random_state']\n", "pbm.hash_(saved_rs_state)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simulation of the series of emission rates\n", "\n", "Here we perform the timetsamps simulation for the list of emission rates previously computed.\n", "\n", "At this point we also choose the Poisson background level to add to the simulation." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "em_rate_list" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simulate timestamps for background = 1kcps:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "rs = np.random.RandomState()\n", "rs.set_state(saved_rs_state)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%timeit -n1 -r1\n", "for em_rate in em_rate_list:\n", " print(' Emission rate: ', em_rate, flush=True)\n", " S.simulate_timestamps_mix(max_rates=(em_rate,), populations=(slice(0, 20),), \n", " bg_rate=1e3, rs=rs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simulate timestamps for background = 4kcps:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%%timeit -n1 -r1\n", "for em_rate in em_rate_list:\n", " print(' Emission rate: ', em_rate, flush=True)\n", " S.simulate_timestamps_mix(max_rates=(em_rate,), populations=(slice(0, 20),), \n", " bg_rate=4e3, rs=rs)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "for k in S.ts_store.h5file.root.timestamps._v_children.keys():\n", " if not k.endswith('_par'):\n", " print(k)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ts, ts_par = S.get_timestamps_part('Pop1_P20_Pstart0_max_rate198000cps_BG4000cps_t_1s_rs_8798a6')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ts[:]" ] }, { "cell_type": "raw", "metadata": {}, "source": [ "array([ 8650, 19320, 22290, ..., 19990800, 19994140, 19996650])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "bins = np.arange(0, 1, 1e-3)\n", "plt.hist(ts*ts.attrs['clk_p'], bins=bins, histtype='step');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Verify the simulation\n", "\n", "Check that the new arrays show up in the data file:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "group = '/timestamps'\n", "\n", "print('Nodes in in %s:\\n' % group)\n", "\n", "print(S.ts_store.h5file.get_node(group))\n", "for node in S.ts_store.h5file.get_node(group)._f_list_nodes():\n", " print('\\t%s' % node.name)\n", " #print('\\t %s' % node.title)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "[t for t in S.timestamp_names if 'BG4000cps' in t]" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "S.ts_store.close()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "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.5.0" } }, "nbformat": 4, "nbformat_minor": 0 }