usALEX-5samples-PR-raw-dir_ex_aa-fit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# usALEX-5samples - Template\n",
    "\n",
    "> *This notebook is executed through [8-spots paper analysis](8-spots paper analysis.ipynb).*\n",
    "> *For a direct execution, uncomment the cell below.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# ph_sel_name = \"all-ph\"\n",
    "# data_id = \"7d\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load software and filenames definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from fretbursts import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "init_notebook()\n",
    "from IPython.display import display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data folder:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_dir = './data/singlespot/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check that the folder exists:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "data_dir = os.path.abspath(data_dir) + '/'\n",
    "assert os.path.exists(data_dir), \"Path '%s' does not exist.\" % data_dir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "List of data files in `data_dir`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from glob import glob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "file_list = sorted(f for f in glob(data_dir + '*.hdf5') if '_BKG' not in f)\n",
    "file_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Selection for POLIMI 2012-12-6 dataset\n",
    "# file_list.pop(2)\n",
    "# file_list = file_list[1:-2]\n",
    "# display(file_list)\n",
    "# labels = ['22d', '27d', '17d', '12d', '7d']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Selection for P.E. 2012-12-6 dataset\n",
    "# file_list.pop(1)\n",
    "# file_list = file_list[:-1]\n",
    "# display(file_list)\n",
    "# labels = ['22d', '27d', '17d', '12d', '7d']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Selection for POLIMI 2012-11-26 datatset\n",
    "labels = ['17d', '27d', '7d', '12d', '22d']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "files_dict = {lab: fname for lab, fname in zip(labels, file_list)}\n",
    "files_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ph_sel_map = {'all-ph': Ph_sel('all'), 'AexAem': Ph_sel(Aex='Aem')}\n",
    "ph_sel = ph_sel_map[ph_sel_name]\n",
    "\n",
    "data_id, ph_sel_name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data load"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initial loading of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d = loader.photon_hdf5(filename=files_dict[data_id])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Laser alternation selection\n",
    "\n",
    "At this point we have only the timestamps and the detector numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d.ph_times_t, d.det_t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to define some parameters: donor and acceptor ch, excitation period and donor and acceptor excitiations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d.add(det_donor_accept=(0, 1), alex_period=4000, D_ON=(2850, 580), A_ON=(900, 2580), offset=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We should check if everithing is OK with an alternation histogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_alternation_hist(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the plot looks good we can apply the parameters with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "loader.alex_apply_period(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Measurements infos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the measurement data is in the `d` variable. We can print it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or check the **measurements duration**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d.time_max"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the background using automatic threshold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d.calc_bg(bg.exp_fit, time_s=60, tail_min_us='auto', F_bg=1.7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(d, timetrace_bg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d.rate_m, d.rate_dd, d.rate_ad, d.rate_aa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Burst search and selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from mpl_toolkits.axes_grid1 import AxesGrid\n",
    "import lmfit\n",
    "print('lmfit version:', lmfit.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "assert d.dir_ex == 0\n",
    "assert d.leakage == 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d.burst_search(m=10, F=6, ph_sel=ph_sel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(d.ph_sel, d.num_bursts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds_sa = d.select_bursts(select_bursts.naa, th1=30)\n",
    "ds_sa.num_bursts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preliminary selection and plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mask = (d.naa[0] - np.abs(d.na[0] + d.nd[0])) > 30\n",
    "ds_saw = d.select_bursts_mask_apply([mask])\n",
    "\n",
    "ds_sas0 = ds_sa.select_bursts(select_bursts.S, S2=0.10)\n",
    "ds_sas = ds_sa.select_bursts(select_bursts.S, S2=0.15)\n",
    "ds_sas2 = ds_sa.select_bursts(select_bursts.S, S2=0.20)\n",
    "ds_sas3 = ds_sa.select_bursts(select_bursts.S, S2=0.25)\n",
    "\n",
    "ds_st = d.select_bursts(select_bursts.size, add_naa=True, th1=30)\n",
    "ds_sas.num_bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dx = ds_sas0\n",
    "size = dx.na[0] + dx.nd[0]\n",
    "s_hist, s_bins = np.histogram(size, bins=np.r_[-15 : 25 : 1], density=True)\n",
    "s_ax = s_bins[:-1] + 0.5*(s_bins[1] - s_bins[0])\n",
    "plot(s_ax, s_hist, '-o', alpha=0.5)\n",
    "\n",
    "dx = ds_sas\n",
    "size = dx.na[0] + dx.nd[0]\n",
    "s_hist, s_bins = np.histogram(size, bins=np.r_[-15 : 25 : 1], density=True)\n",
    "s_ax = s_bins[:-1] + 0.5*(s_bins[1] - s_bins[0])\n",
    "plot(s_ax, s_hist, '-o', alpha=0.5)\n",
    "\n",
    "dx = ds_sas2\n",
    "size = dx.na[0] + dx.nd[0]\n",
    "s_hist, s_bins = np.histogram(size, bins=np.r_[-15 : 25 : 1], density=True)\n",
    "s_ax = s_bins[:-1] + 0.5*(s_bins[1] - s_bins[0])\n",
    "plot(s_ax, s_hist, '-o', alpha=0.5)\n",
    "\n",
    "dx = ds_sas3\n",
    "size = dx.na[0] + dx.nd[0]\n",
    "s_hist, s_bins = np.histogram(size, bins=np.r_[-15 : 25 : 1], density=True)\n",
    "s_ax = s_bins[:-1] + 0.5*(s_bins[1] - s_bins[0])\n",
    "plot(s_ax, s_hist, '-o', alpha=0.5)\n",
    "\n",
    "plt.title('(nd + na) for A-only population using different S cutoff');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dx = ds_sa\n",
    "\n",
    "alex_jointplot(dx);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(ds_sa, hist_S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A-direct excitation fitting\n",
    "\n",
    "To extract the A-direct excitation coefficient we need to fit the \n",
    "S values for the A-only population.\n",
    "\n",
    "The S value for the A-only population is fitted with different methods:\n",
    "- Histogram git with 2 Gaussians or with 2 asymmetric Gaussians \n",
    "(an asymmetric Gaussian has right- and left-side of the peak\n",
    "decreasing according to different sigmas).\n",
    "- KDE maximum\n",
    "\n",
    "In the following we apply these methods using different selection\n",
    "or weighting schemes to reduce amount of FRET population and make\n",
    "fitting of the A-only population easier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Even selection\n",
    "\n",
    "Here A-only and FRET population are evenly selected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dx = ds_sa\n",
    "\n",
    "bin_width = 0.03\n",
    "bandwidth = 0.03\n",
    "bins = np.r_[-0.2 : 1 : bin_width]\n",
    "x_kde = np.arange(bins.min(), bins.max(), 0.0002)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Weights\n",
    "weights = None\n",
    "\n",
    "## Histogram fit\n",
    "fitter_g = mfit.MultiFitter(dx.S)\n",
    "fitter_g.histogram(bins=np.r_[-0.2 : 1.2 : bandwidth])\n",
    "fitter_g.fit_histogram(model = mfit.factory_two_gaussians(p1_center=0.1, p2_center=0.4))\n",
    "S_hist_orig = fitter_g.hist_pdf\n",
    "\n",
    "S_2peaks = fitter_g.params.loc[0, 'p1_center']\n",
    "dir_ex_S2p = S_2peaks/(1 - S_2peaks)\n",
    "print('Fitted direct excitation (na/naa) [2-Gauss]:', dir_ex_S2p)\n",
    "\n",
    "## KDE\n",
    "fitter_g.calc_kde(bandwidth=bandwidth)\n",
    "fitter_g.find_kde_max(x_kde, xmin=0, xmax=0.15)\n",
    "\n",
    "S_peak = fitter_g.kde_max_pos[0]\n",
    "dir_ex_S_kde = S_peak/(1 - S_peak)\n",
    "print('Fitted direct excitation (na/naa) [KDE]:    ', dir_ex_S_kde)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(14, 4.5))\n",
    "\n",
    "mfit.plot_mfit(fitter_g, ax=ax[0])\n",
    "ax[0].set_title('2-Gaussians fit (S_fit = %.2f %%)' % (S_2peaks*100))\n",
    "\n",
    "mfit.plot_mfit(fitter_g, ax=ax[1], plot_model=False, plot_kde=True)\n",
    "ax[1].set_title('KDE fit (S_fit = %.2f %%)' % (S_peak*100));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## 2-Asym-Gaussian\n",
    "fitter_ag = mfit.MultiFitter(dx.S)\n",
    "fitter_ag.histogram(bins=np.r_[-0.2 : 1.2 : bandwidth])\n",
    "fitter_ag.fit_histogram(model = mfit.factory_two_asym_gaussians(p1_center=0.1, p2_center=0.4))\n",
    "#print(fitter_ag.fit_obj[0].model.fit_report())\n",
    "\n",
    "S_2peaks_a = fitter_ag.params.loc[0, 'p1_center']\n",
    "dir_ex_S2pa = S_2peaks_a/(1 - S_2peaks_a)\n",
    "print('Fitted direct excitation (na/naa) [2-Gauss]:', dir_ex_S2pa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(14, 4.5))\n",
    "\n",
    "mfit.plot_mfit(fitter_g, ax=ax[0])\n",
    "ax[0].set_title('2-Gaussians fit (S_fit = %.2f %%)' % (S_2peaks*100))\n",
    "\n",
    "mfit.plot_mfit(fitter_ag, ax=ax[1])\n",
    "ax[1].set_title('2-Asym-Gaussians fit (S_fit = %.2f %%)' % (S_2peaks_a*100));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zero threshold on nd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select bursts with:\n",
    "\n",
    "$$n_d < 0$$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dx = ds_sa.select_bursts(select_bursts.nd, th1=-100, th2=0)\n",
    "\n",
    "fitter = bext.bursts_fitter(dx, 'S')\n",
    "fitter.fit_histogram(model = mfit.factory_gaussian(center=0.1))\n",
    "S_1peaks_th = fitter.params.loc[0, 'center']\n",
    "dir_ex_S1p = S_1peaks_th/(1 - S_1peaks_th)\n",
    "print('Fitted direct excitation (na/naa) [2-Gauss]:', dir_ex_S1p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mfit.plot_mfit(fitter)\n",
    "plt.xlim(-0.1, 0.6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selection 1\n",
    "\n",
    "Bursts are weighted using $w = f(S)$, where the function $f(S)$ is a\n",
    "Gaussian fitted to the $S$ histogram of the FRET population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dx = ds_sa\n",
    "\n",
    "## Weights\n",
    "weights = 1 - mfit.gaussian(dx.S[0], fitter_g.params.loc[0, 'p2_center'], fitter_g.params.loc[0, 'p2_sigma'])\n",
    "weights[dx.S[0] >= fitter_g.params.loc[0, 'p2_center']] = 0\n",
    "\n",
    "## Histogram fit\n",
    "fitter_w1 = mfit.MultiFitter(dx.S)\n",
    "fitter_w1.weights = [weights]\n",
    "fitter_w1.histogram(bins=np.r_[-0.2 : 1.2 : bandwidth])\n",
    "fitter_w1.fit_histogram(model = mfit.factory_two_gaussians(p1_center=0.1, p2_center=0.4))\n",
    "S_2peaks_w1 = fitter_w1.params.loc[0, 'p1_center']\n",
    "dir_ex_S2p_w1 = S_2peaks_w1/(1 - S_2peaks_w1)\n",
    "print('Fitted direct excitation (na/naa) [2-Gauss]:', dir_ex_S2p_w1)\n",
    "\n",
    "## KDE\n",
    "fitter_w1.calc_kde(bandwidth=bandwidth)\n",
    "fitter_w1.find_kde_max(x_kde, xmin=0, xmax=0.15)\n",
    "S_peak_w1 = fitter_w1.kde_max_pos[0]\n",
    "dir_ex_S_kde_w1 = S_peak_w1/(1 - S_peak_w1)\n",
    "print('Fitted direct excitation (na/naa) [KDE]:    ', dir_ex_S_kde_w1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_weights(x, weights, ax):\n",
    "    ax2 = ax.twinx()\n",
    "    x_sort = x.argsort()\n",
    "    ax2.plot(x[x_sort], weights[x_sort], color='k', lw=4, alpha=0.4)\n",
    "    ax2.set_ylabel('Weights');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(14, 4.5))\n",
    "mfit.plot_mfit(fitter_w1, ax=ax[0])\n",
    "mfit.plot_mfit(fitter_g, ax=ax[0], plot_model=False, plot_kde=False)\n",
    "plot_weights(dx.S[0], weights, ax=ax[0])\n",
    "ax[0].set_title('2-Gaussians fit (S_fit = %.2f %%)' % (S_2peaks_w1*100))\n",
    "\n",
    "mfit.plot_mfit(fitter_w1, ax=ax[1], plot_model=False, plot_kde=True)\n",
    "mfit.plot_mfit(fitter_g, ax=ax[1], plot_model=False, plot_kde=False)\n",
    "plot_weights(dx.S[0], weights, ax=ax[1])\n",
    "ax[1].set_title('KDE fit (S_fit = %.2f %%)' % (S_peak_w1*100));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selection 2\n",
    "\n",
    "Bursts are here weighted using weights $w$:\n",
    "\n",
    "$$w = n_{aa} - |n_a + n_d|$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Weights\n",
    "sizes = dx.nd[0] + dx.na[0] #- dir_ex_S_kde_w3*dx.naa[0]\n",
    "weights = dx.naa[0] - abs(sizes)\n",
    "weights[weights < 0] = 0\n",
    "\n",
    "## Histogram\n",
    "fitter_w4 = mfit.MultiFitter(dx.S)\n",
    "fitter_w4.weights = [weights]\n",
    "fitter_w4.histogram(bins=np.r_[-0.2 : 1.2 : bandwidth])\n",
    "fitter_w4.fit_histogram(model = mfit.factory_two_gaussians(p1_center=0.1, p2_center=0.4))\n",
    "S_2peaks_w4 = fitter_w4.params.loc[0, 'p1_center']\n",
    "dir_ex_S2p_w4 = S_2peaks_w4/(1 - S_2peaks_w4)\n",
    "print('Fitted direct excitation (na/naa) [2-Gauss]:', dir_ex_S2p_w4)\n",
    "\n",
    "## KDE\n",
    "fitter_w4.calc_kde(bandwidth=bandwidth)\n",
    "fitter_w4.find_kde_max(x_kde, xmin=0, xmax=0.15)\n",
    "S_peak_w4 = fitter_w4.kde_max_pos[0]\n",
    "dir_ex_S_kde_w4 = S_peak_w4/(1 - S_peak_w4)\n",
    "print('Fitted direct excitation (na/naa) [KDE]:    ', dir_ex_S_kde_w4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(14, 4.5))\n",
    "\n",
    "mfit.plot_mfit(fitter_w4, ax=ax[0])\n",
    "mfit.plot_mfit(fitter_g, ax=ax[0], plot_model=False, plot_kde=False)\n",
    "#plot_weights(dx.S[0], weights, ax=ax[0])\n",
    "ax[0].set_title('2-Gaussians fit (S_fit = %.2f %%)' % (S_2peaks_w4*100))\n",
    "\n",
    "mfit.plot_mfit(fitter_w4, ax=ax[1], plot_model=False, plot_kde=True)\n",
    "mfit.plot_mfit(fitter_g, ax=ax[1], plot_model=False, plot_kde=False)\n",
    "#plot_weights(dx.S[0], weights, ax=ax[1])\n",
    "ax[1].set_title('KDE fit (S_fit = %.2f %%)' % (S_peak_w4*100));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selection 3\n",
    "\n",
    "Bursts are here selected according to:\n",
    "\n",
    "$$n_{aa} - |n_a + n_d| > 30$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mask = (d.naa[0] - np.abs(d.na[0] + d.nd[0])) > 30\n",
    "ds_saw = d.select_bursts_mask_apply([mask])\n",
    "print(ds_saw.num_bursts)\n",
    "\n",
    "dx = ds_saw\n",
    "\n",
    "## Weights\n",
    "weights = None\n",
    "\n",
    "## 2-Gaussians\n",
    "fitter_w5 = mfit.MultiFitter(dx.S)\n",
    "fitter_w5.histogram(bins=np.r_[-0.2 : 1.2 : bandwidth])\n",
    "fitter_w5.fit_histogram(model = mfit.factory_two_gaussians(p1_center=0.1, p2_center=0.4))\n",
    "S_2peaks_w5 = fitter_w5.params.loc[0, 'p1_center']\n",
    "dir_ex_S2p_w5 = S_2peaks_w5/(1 - S_2peaks_w5)\n",
    "print('Fitted direct excitation (na/naa) [2-Gauss]:', dir_ex_S2p_w5)\n",
    "\n",
    "## KDE\n",
    "fitter_w5.calc_kde(bandwidth=bandwidth)\n",
    "fitter_w5.find_kde_max(x_kde, xmin=0, xmax=0.15)\n",
    "S_peak_w5 = fitter_w5.kde_max_pos[0]\n",
    "S_2peaks_w5_fiterr = fitter_w5.fit_res[0].params['p1_center'].stderr\n",
    "dir_ex_S_kde_w5 = S_peak_w5/(1 - S_peak_w5)\n",
    "print('Fitted direct excitation (na/naa) [KDE]:    ', dir_ex_S_kde_w5)\n",
    "\n",
    "## 2-Asym-Gaussians\n",
    "fitter_w5a = mfit.MultiFitter(dx.S)\n",
    "fitter_w5a.histogram(bins=np.r_[-0.2 : 1.2 : bandwidth])\n",
    "fitter_w5a.fit_histogram(model = mfit.factory_two_asym_gaussians(p1_center=0.05, p2_center=0.3))\n",
    "S_2peaks_w5a = fitter_w5a.params.loc[0, 'p1_center']\n",
    "dir_ex_S2p_w5a = S_2peaks_w5a/(1 - S_2peaks_w5a)\n",
    "#print(fitter_w5a.fit_obj[0].model.fit_report(min_correl=0.5))\n",
    "print('Fitted direct excitation (na/naa) [2-Asym-Gauss]:', dir_ex_S2p_w5a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 3, figsize=(19, 4.5))\n",
    "\n",
    "mfit.plot_mfit(fitter_w5, ax=ax[0])\n",
    "mfit.plot_mfit(fitter_g, ax=ax[0], plot_model=False, plot_kde=False)\n",
    "ax[0].set_title('2-Gaussians fit (S_fit = %.2f %%)' % (S_2peaks_w5*100))\n",
    "\n",
    "mfit.plot_mfit(fitter_w5, ax=ax[1], plot_model=False, plot_kde=True)\n",
    "mfit.plot_mfit(fitter_g, ax=ax[1], plot_model=False, plot_kde=False)\n",
    "ax[1].set_title('KDE fit (S_fit = %.2f %%)' % (S_peak_w5*100));\n",
    "\n",
    "mfit.plot_mfit(fitter_w5a, ax=ax[2])\n",
    "mfit.plot_mfit(fitter_g, ax=ax[2], plot_model=False, plot_kde=False)\n",
    "ax[2].set_title('2-Asym-Gaussians fit (S_fit = %.2f %%)' % (S_2peaks_w5a*100));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save data to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sample = data_id\n",
    "n_bursts_aa = ds_sas.num_bursts[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following string contains the list of variables to be saved. When saving, the order of the variables is preserved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "variables = ('sample n_bursts_aa dir_ex_S1p dir_ex_S_kde dir_ex_S2p dir_ex_S2pa '\n",
    "             'dir_ex_S2p_w1 dir_ex_S_kde_w1 dir_ex_S_kde_w4 dir_ex_S_kde_w5 dir_ex_S2p_w5 dir_ex_S2p_w5a '\n",
    "             'S_2peaks_w5 S_2peaks_w5_fiterr\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is just a trick to format the different variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "variables_csv = variables.replace(' ', ',')\n",
    "fmt_float = '{%s:.6f}'\n",
    "fmt_int = '{%s:d}'\n",
    "fmt_str = '{%s}'\n",
    "fmt_dict = {**{'sample': fmt_str}, \n",
    "            **{k: fmt_int for k in variables.split() if k.startswith('n_bursts')}}\n",
    "var_dict = {name: eval(name) for name in variables.split()}\n",
    "var_fmt = ', '.join([fmt_dict.get(name, fmt_float) % name for name in variables.split()]) + '\\n'\n",
    "data_str = var_fmt.format(**var_dict)\n",
    "\n",
    "print(variables_csv)\n",
    "print(data_str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# NOTE: The file name should be the notebook name but with .csv extension\n",
    "with open('results/usALEX-5samples-PR-raw-dir_ex_aa-fit-%s.csv' % ph_sel_name, 'a') as f:\n",
    "    f.seek(0, 2)\n",
    "    if f.tell() == 0:\n",
    "        f.write(variables_csv)\n",
    "    f.write(data_str)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}