demo_caiman_pipeline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<html><head><meta content=\"text/html; charset=UTF-8\" http-equiv=\"content-type\"><style type=\"text/css\">ol</style></head><body class=\"c5\"><p class=\"c0 c4\"><span class=\"c3\"></span></p><p class=\"c2 title\" id=\"h.rrbabt268i6e\"><h1>CaImAn&rsquo;s Demo pipeline</h1></p><p class=\"c0\"><span class=\"c3\">This notebook will help to demonstrate the process of CaImAn and how it uses different functions to denoise, deconvolve and demix neurons from a Calcium Imaging Video. </span></p>\n",
    "<p><img src=\"docs/img/quickintro.png\" /></p>\n",
    "<p class=\"c0\"><span class=\"c3\">More information can be found in CaImAn&rsquo;s documentation. </span></p>\n",
    "</html>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import YouTubeVideo as yt\n",
    "from __future__ import division\n",
    "from __future__ import print_function\n",
    "from builtins import zip\n",
    "from builtins import str\n",
    "from builtins import map\n",
    "from builtins import range\n",
    "from past.utils import old_div\n",
    "import cv2\n",
    "try:\n",
    "    cv2.setNumThreads(1)\n",
    "except:\n",
    "    print('Open CV is naturally single threaded')\n",
    "\n",
    "try:\n",
    "    if __IPYTHON__:\n",
    "        print((1))\n",
    "        # this is used for debugging purposes only. allows to reload classes\n",
    "        # when changed\n",
    "        get_ipython().magic('load_ext autoreload')\n",
    "        get_ipython().magic('autoreload 2')\n",
    "except NameError:\n",
    "    print('Not IPYTHON')\n",
    "    pass\n",
    "\n",
    "import caiman as cm\n",
    "import numpy as np\n",
    "import os\n",
    "import glob\n",
    "import time\n",
    "import pylab as pl\n",
    "import psutil\n",
    "import sys\n",
    "from ipyparallel import Client\n",
    "from skimage.external.tifffile import TiffFile\n",
    "import scipy\n",
    "\n",
    "from caiman.motion_correction import tile_and_correct, motion_correction_piecewise\n",
    "from caiman.source_extraction.cnmf import cnmf as cnmf\n",
    "from caiman.motion_correction import MotionCorrect\n",
    "from caiman.components_evaluation import evaluate_components \n",
    "from caiman.utils.visualization import plot_contours, view_patches_bar\n",
    "from caiman.base.rois import extract_binary_masks_blob\n",
    "from caiman.utils.utils import download_demo\n",
    "from caiman.utils.visualization import plot_contours,view_patches_bar,nb_plot_contour,nb_view_patches\n",
    "import bokeh.plotting as bpl\n",
    "from IPython.lib.display import YouTubeVideo as yt\n",
    "bpl.output_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "params_movie = {'fname': ['Sue_2x_3000_40_-46.tif'],\n",
    "                'niter_rig': 1,\n",
    "                'max_shifts': (6, 6),  # maximum allow rigid shift\n",
    "                'splits_rig': 56,  # for parallelization split the movies in  num_splits chuncks across time\n",
    "                # if none all the splits are processed and the movie is saved\n",
    "                'num_splits_to_process_rig': None,\n",
    "                # intervals at which patches are laid out for motion correction\n",
    "                'strides': (48, 48),\n",
    "                # overlap between pathes (size of patch strides+overlaps)\n",
    "                'overlaps': (24, 24),\n",
    "                'splits_els': 56,  # for parallelization split the movies in  num_splits chuncks across time\n",
    "                # if none all the splits are processed and the movie is saved\n",
    "                'num_splits_to_process_els': [28, None],\n",
    "                'upsample_factor_grid': 4,  # upsample factor to avoid smearing when merging patches\n",
    "                # maximum deviation allowed for patch with respect to rigid\n",
    "                # shift\n",
    "                'max_deviation_rigid': 3,\n",
    "                'p': 1,  # order of the autoregressive system\n",
    "                'merge_thresh': 0.8,  # merging threshold, max correlation allowed\n",
    "                'rf': 15,  # half-size of the patches in pixels. rf=25, patches are 50x50\n",
    "                'stride_cnmf': 6,  # amounpl.it of overlap between the patches in pixels\n",
    "                'K': 4,  # number of components per patch\n",
    "                # if dendritic. In this case you need to set init_method to\n",
    "                # sparse_nmf\n",
    "                'is_dendrites': False,\n",
    "                'init_method': 'greedy_roi',\n",
    "                'gSig': [4, 4],  # expected half size of neurons\n",
    "                'alpha_snmf': None,  # this controls sparsity\n",
    "                'final_frate': 30\n",
    "                }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#%% parameters from dictionary\n",
    "fname = params_movie['fname']\n",
    "niter_rig = params_movie['niter_rig']\n",
    "# maximum allow rigid shift\n",
    "max_shifts = params_movie['max_shifts']  \n",
    "# for parallelization split the movies in  num_splits chuncks across time\n",
    "splits_rig = params_movie['splits_rig']  \n",
    "# if none all the splits are processed and the movie is saved\n",
    "num_splits_to_process_rig = params_movie['num_splits_to_process_rig']\n",
    "# intervals at which patches are laid out for motion correction\n",
    "strides = params_movie['strides']\n",
    "# overlap between pathes (size of patch strides+overlaps)\n",
    "overlaps = params_movie['overlaps']\n",
    "# for parallelization split the movies in  num_splits chuncks across time\n",
    "splits_els = params_movie['splits_els'] \n",
    "# if none all the splits are processed and the movie is saved\n",
    "num_splits_to_process_els = params_movie['num_splits_to_process_els']\n",
    "# upsample factor to avoid smearing when merging patches\n",
    "upsample_factor_grid = params_movie['upsample_factor_grid'] \n",
    "# maximum deviation allowed for patch with respect to rigid\n",
    "# shift\n",
    "max_deviation_rigid = params_movie['max_deviation_rigid']\n",
    "#%% some parameter settings\n",
    "# order of the autoregressive fit to calcium imaging in general one (slow gcamps) or two (fast gcamps fast scanning)\n",
    "p = params_movie['p']  \n",
    "# merging threshold, max correlation allowed\n",
    "merge_thresh= params_movie['merge_thresh'] \n",
    "# half-size of the patches in pixels. rf=25, patches are 50x50\n",
    "rf = params_movie['rf']  \n",
    "# amounpl.it of overlap between the patches in pixels\n",
    "stride_cnmf = params_movie['stride_cnmf'] \n",
    " # number of components per patch\n",
    "K =  params_movie['K'] \n",
    "# if dendritic. In this case you need to set init_method to sparse_nmf\n",
    "is_dendrites = params_movie['is_dendrites']\n",
    "# iinit method can be greedy_roi for round shapes or sparse_nmf for denritic data\n",
    "init_method = params_movie['init_method']\n",
    "# expected half size of neurons\n",
    "gSig = params_movie['gSig']  \n",
    "# this controls sparsity\n",
    "alpha_snmf = params_movie['alpha_snmf']  \n",
    "#frame rate of movie (even considering eventual downsampling)\n",
    "final_frate = params_movie['final_frate']\n",
    "#%% Extract spatial and temporal components on patches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# %% download movie if not there                                                                                                                                                                                \n",
    "if fname[0] in ['Sue_2x_3000_40_-46.tif','demoMovieJ.tif']:\n",
    "    download_demo(fname[0])\n",
    "    fname = [os.path.join('example_movies',fname[0])]\n",
    "m_orig = cm.load_movie_chain(fname[:1])\n",
    "yt('mLvNmPs9onE')\n",
    "print(\"Data from Sue Ann Koay, Princeton Neuroscience Institute, Tank Lab.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% load movie (in memory!)\n",
    "m_orig = cm.load_movie_chain(fname)\n",
    "#%% play movie\n",
    "downsample_ratio = .2\n",
    "offset_mov = -np.min(m_orig[:100])\n",
    "m_orig.resize(1, 1, downsample_ratio).play(\n",
    "    gain=10, offset = offset_mov, fr=30, magnification=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<em>some of the pixels in this movie are negative, we then need to make them positive </em>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Motion correction is performed in parallel on chunks taken across times. </h2>\n",
    "<p> Create temporal chunks of the movie for  parallel processing on all cores </p>\n",
    "<p> ipyparallel is used to create a cluster that handles the parallelization eithr on the PC ( see : https://ipyparallel.readthedocs.io/en/latest/intro.html) or on clusters interfacing with slurm ( see : https://slurm.schedmd.com/quickstart.html ) </p>\n",
    "<p><img src=\"docs/img/cordermmap.png\" /> </p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% start the cluster\n",
    "c, dview, n_processes = cm.cluster.setup_cluster(\n",
    "    backend='local', n_processes=None, single_thread=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> Rigid motion correction</h1>\n",
    "<p> We are now left with a really shacky video. In order to correct the image movement we are using a simple rigid motion correction algorithm. This algorithm first creates a correlation image over time using frames from the video. It then tries to match each frame to this template. In addition the template will get updated during the matching process. Making it more precise and so does the template matching.  </p>\n",
    "<img src=\"docs/img/rigidcorrection.png\" />\n",
    "more info : <em> http://opencv.org/about.html </em>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# movie must be mostly positive for this to work\n",
    "min_mov = cm.load(fname[0], subindices=range(400)).min()\n",
    "\n",
    "mc = MotionCorrect(fname[0], min_mov,\n",
    "                   dview=dview, max_shifts=max_shifts, niter_rig=niter_rig,\n",
    "                   splits_rig=splits_rig, \n",
    "                   num_splits_to_process_rig=num_splits_to_process_rig, \n",
    "                strides= strides, overlaps= overlaps, splits_els=splits_els,\n",
    "                num_splits_to_process_els=num_splits_to_process_els, \n",
    "                upsample_factor_grid=upsample_factor_grid,\n",
    "                   max_deviation_rigid=max_deviation_rigid, \n",
    "                shifts_opencv = True, nonneg_movie = True)\n",
    "\n",
    "mc.motion_correct_rigid(save_movie=True)\n",
    "# load motion corrected movie\n",
    "m_rig = cm.load(mc.fname_tot_rig)\n",
    "bord_px_rig = np.ceil(np.max(mc.shifts_rig)).astype(np.int)\n",
    "#%% visualize templates\n",
    "pl.figure(figsize = (20,10))\n",
    "pl.imshow(mc.total_template_rig, cmap = 'gray')\n",
    "#%% inspect movie\n",
    "m_rig.resize(1, 1, downsample_ratio).play(\n",
    "    gain=10, offset = offset_mov*.25, fr=30, magnification=2,bord_px = bord_px_rig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Inspect motion correction results</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% plot rigid shifts\n",
    "pl.close()\n",
    "pl.figure(figsize = (20,10))\n",
    "pl.plot(mc.shifts_rig)\n",
    "pl.legend(['x shifts','y shifts'])\n",
    "pl.xlabel('frames')\n",
    "pl.ylabel('pixels')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> Piecewise rigid motion correction </h1>\n",
    "<p> our video is now way less shaky but since the brain tissue is a compressible object, a non  rigid deformation still exist and in order to do a really nice source separation, needs to be corrected : </p>\n",
    "<p> Here we introduce a function for a fast Non-Rigid Motion Correction based on template matching.<br/>It is a similar approach than the one presented before but on overlapping spatial patches. The estimated alignments are subsequently up-sampled to create a smooth motion field for each frame that can efficiently approximate non-rigid motion in a piecewise-rigid manner.<br/>It can be run in an online mode resulting in comparable to or even\n",
    "faster than real time motion registration on streaming data.<p>\n",
    "<img src=\"docs\\img\\pwrigidcorrection.png\" />\n",
    "<p> more info : (Eftychios A. Pnevmatikakis, Andrea Giovannucci, NormCorre) </p>\n",
    "<em> http://biorxiv.org/content/biorxiv/early/2017/02/14/108514.full.pdf </em>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% motion correct piecewise rigid\n",
    "mc.motion_correct_pwrigid(save_movie=True, template=mc.total_template_rig,\n",
    "                          show_template = True)\n",
    "m_els = cm.load(mc.fname_tot_els)\n",
    "pl.figure(figsize = (20,10))\n",
    "pl.imshow(mc.total_template_els, cmap = 'gray')\n",
    "yt('JzierjDuBzg')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% visualize elastic shifts\n",
    "pl.close()\n",
    "pl.figure(figsize = (20,10))\n",
    "pl.subplot(2, 1, 1)\n",
    "pl.plot(mc.x_shifts_els)\n",
    "pl.ylabel('x shifts (pixels)')\n",
    "pl.subplot(2, 1, 2)\n",
    "pl.plot(mc.y_shifts_els)\n",
    "pl.ylabel('y_shifts (pixels)')\n",
    "pl.xlabel('frames')\n",
    "#%% compute borders to exclude\n",
    "bord_px_els = np.ceil(np.maximum(np.max(np.abs(mc.x_shifts_els)),\n",
    "                                 np.max(np.abs(mc.y_shifts_els)))).astype(np.int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check movie\n",
    "downsample_ratio = .2\n",
    "m_els.resize(1, 1, downsample_ratio).play(\n",
    "    gain=10, offset = 0, fr=30, magnification=2,bord_px = bord_px_els)\n",
    "# compare with original and rigid corrected movies\n",
    "downsample_factor = .2\n",
    "cm.concatenate([m_orig.resize(1, 1, downsample_factor)+offset_mov,\n",
    "                m_rig.resize(1, 1, downsample_factor), m_els.resize(\n",
    "    1, 1, downsample_factor)], axis=2).play(fr=60, gain=15, magnification=2, offset=0)\n",
    "#%% local correlation\n",
    "pl.figure(figsize = (20,10))\n",
    "pl.imshow(m_els.local_correlations(eight_neighbours=True, swap_dim=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> Assessing Motion correction Quality </h1>\n",
    "for the raw, rigid corrected and piecewise rigid corrected videos\n",
    "<p> using smoothness <em> (mean over time) </em> </p>\n",
    "--------------\n",
    "<p> using correlation to the template <em> (Pearson Correlation coefficient)</em> </p> \n",
    "--------------\n",
    "see : http://docs.opencv.org/2.4/doc/tutorials/imgproc/histograms/template_matching/template_matching.html ( like a normalized SSD )\n",
    "<p> optical flow : </p>\n",
    "--------------\n",
    "<img src=\"docs/img/opticalflow.png\" />\n",
    "more info :<em> http://docs.opencv.org/trunk/d7/d8b/tutorial_py_lucas_kanade.html </em>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#% compute metrics for the results (TAKES TIME!!)\n",
    "final_size = np.subtract(mc.total_template_els.shape, 2 * bord_px_els)\n",
    "winsize = 100\n",
    "swap_dim = False\n",
    "resize_fact_flow = .2\n",
    "tmpl, correlations, flows_orig, norms, smoothness = cm.motion_correction.compute_metrics_motion_correction(\n",
    "    mc.fname_tot_els, final_size[0], final_size[1],\n",
    "    swap_dim, winsize=winsize, play_flow=False, resize_fact_flow=resize_fact_flow)\n",
    "\n",
    "tmpl, correlations, flows_orig, norms, smoothness = cm.motion_correction.compute_metrics_motion_correction(\n",
    "    mc.fname_tot_rig, final_size[0], final_size[1],\n",
    "    swap_dim, winsize=winsize, play_flow=False, resize_fact_flow=resize_fact_flow)\n",
    "\n",
    "tmpl, correlations, flows_orig, norms, smoothness = cm.motion_correction.compute_metrics_motion_correction(\n",
    "    fname[0], final_size[0], final_size[1], swap_dim, winsize=winsize, play_flow=False, resize_fact_flow=resize_fact_flow)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% plot the results of metrics\n",
    "fls = [mc.fname_tot_els[:-4] + '_metrics.npz', mc.fname_tot_rig[:-4] +\n",
    "       '_metrics.npz', mc.fname[:-4] + '_metrics.npz']\n",
    "\n",
    "pl.figure(figsize = (20,10))\n",
    "for cnt, fl, metr in zip(range(len(fls)),fls,['pw_rigid','rigid','raw']):\n",
    "    with np.load(fl) as ld:\n",
    "        print(ld.keys())\n",
    "        print(fl)\n",
    "        print(str(np.mean(ld['norms'])) + '+/-' + str(np.std(ld['norms'])) +\n",
    "              ' ; ' + str(ld['smoothness']) + ' ; ' + str(ld['smoothness_corr']))\n",
    "        \n",
    "        pl.subplot(len(fls), 4, 1 + 4 * cnt)\n",
    "        pl.ylabel(metr)\n",
    "        try:\n",
    "            mean_img = np.mean(\n",
    "            cm.load(fl[:-12] + 'mmap'), 0)[12:-12, 12:-12]\n",
    "        except:\n",
    "            try:\n",
    "                mean_img = np.mean(\n",
    "                    cm.load(fl[:-12] + '.tif'), 0)[12:-12, 12:-12]\n",
    "            except:\n",
    "                mean_img = np.mean(\n",
    "                    cm.load(fl[:-12] + 'hdf5'), 0)[12:-12, 12:-12]\n",
    "                    \n",
    "        lq, hq = np.nanpercentile(mean_img, [.5, 99.5])\n",
    "        pl.imshow(mean_img, vmin=lq, vmax=hq)\n",
    "        pl.title('Mean')\n",
    "        pl.subplot(len(fls), 4, 4 * cnt + 2)\n",
    "        pl.imshow(ld['img_corr'], vmin=0, vmax=.35)\n",
    "        pl.title('Corr image')\n",
    "        pl.subplot(len(fls), 4, 4 * cnt + 3)\n",
    "        pl.plot(ld['norms'])\n",
    "        pl.xlabel('frame')\n",
    "        pl.ylabel('norm opt flow')\n",
    "        pl.subplot(len(fls), 4, 4 * cnt + 4)\n",
    "        flows = ld['flows']\n",
    "        pl.imshow(np.mean(\n",
    "        np.sqrt(flows[:, :, :, 0]**2 + flows[:, :, :, 1]**2), 0), vmin=0, vmax=0.3)\n",
    "        pl.colorbar()\n",
    "        pl.title('Mean optical flow')      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> Memory mapping </h1>\n",
    "\n",
    "<p> We want the parallel processes to access and our video matrix without having it in memory and duplicating it. we are using the memapping method to fullfill such a goal</p>\n",
    "<p> learn more : https://en.wikipedia.org/wiki/Memory-mapped_file </p>\n",
    "<p><img src=\"docs/img/fordermmap.png\" /></p>\n",
    "<h3> Data is saved on drive so that extracting blocks of movies is efficient.  </h3>\n",
    "<p>Data can be either quickly read per columns or per rows but not in both directions. So, we save the movie in such a way that it is easy to read across time.</p>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MEMORY MAPPING: save each chunk in F format on memory mapped files\n",
    "if 'max_shifts' not in params_movie:\n",
    "    fnames = [params_movie['fname']]\n",
    "    border_to_0 = 0\n",
    "elif not 'overlaps'in params_movie:\n",
    "    fnames = [mc.fname_tot_rig]\n",
    "    border_to_0 = bord_px_rig\n",
    "    m_els = m_rig\n",
    "else:\n",
    "    fnames = [mc.fname_tot_els]\n",
    "    border_to_0 = bord_px_els\n",
    "    \n",
    "idx_xy = None\n",
    "add_to_movie = -np.nanmin(m_els) + 1  # movie must be positive\n",
    "# if you need to remove frames from the beginning of each file\n",
    "remove_init = 0\n",
    "# downsample movie in time: use .2 or .1 if file is large and you want a quick answer             \n",
    "downsample_factor = 1 \n",
    "base_name = fname[0].split('/')[-1][:-4]\n",
    "name_new = cm.save_memmap_each(fnames, dview=dview, base_name=base_name, resize_fact=(\n",
    "    1, 1, downsample_factor), remove_init=remove_init, idx_xy=idx_xy,\n",
    "                               add_to_movie=add_to_movie, border_to_0=border_to_0)\n",
    "name_new.sort()\n",
    "print(name_new)\n",
    "\n",
    "#%% concatenate chunks if needed\n",
    "if len(name_new) > 1:\n",
    "    fname_new = cm.save_memmap_join(\n",
    "        name_new, base_name='Yr', n_chunks=12, dview=dview)\n",
    "else:\n",
    "    print('One file only, not saving!')\n",
    "    fname_new = name_new[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% LOAD MEMMAP FILE\n",
    "Yr, dims, T = cm.load_memmap(fname_new)\n",
    "d1, d2 = dims\n",
    "images = np.reshape(Yr.T, [T] + list(dims), order='F')\n",
    "Y = np.reshape(Yr, dims + (T,), order='F')\n",
    "m_images = cm.movie(images)\n",
    "#%%  checks on movies (might take time if large!)\n",
    "if np.min(images) < 0:\n",
    "    raise Exception('Movie too negative, add_to_movie should be larger')\n",
    "if np.sum(np.isnan(images)) > 0:\n",
    "    raise Exception('Movie contains nan! You did not remove enough borders')\n",
    "#%% correlation image\n",
    "Cn = cm.local_correlations(Y)\n",
    "Cn[np.isnan(Cn)] = 0\n",
    "pl.imshow(Cn, cmap='gray', vmax=.35)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> CNMF on patches </h1>\n",
    "<p> Our method builds upon and extend (Pnevmatikakis et al, Neuron, 2016)  <em> http://www.cell.com/neuron/fulltext/S0896-6273(15)01084-3 </em> </p>\n",
    "\n",
    "<p> <img src=\"docs/img/cnmf1.png\" /> </p>\n",
    "\n",
    "<p> We run portion of the CNMF algorithm  on patches </p>\n",
    "<p> We merge the patches with special attention to neurons on the border. </p>\n",
    "<p> We refine the result by rerunning CNMF on the whole movie </p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnm = cnmf.CNMF(n_processes, k=K, gSig=gSig, merge_thresh=0.8, p=0, \n",
    "                dview=dview, Ain=None, rf=rf, stride=stride_cnmf, memory_fact=1,\n",
    "                method_init=init_method, alpha_snmf=alpha_snmf,\n",
    "                only_init_patch=True, gnb=1, method_deconvolution='oasis')\n",
    "cnm = cnm.fit(images)\n",
    "\n",
    "A_tot = cnm.A\n",
    "C_tot = cnm.C\n",
    "YrA_tot = cnm.YrA\n",
    "b_tot = cnm.b\n",
    "f_tot = cnm.f\n",
    "sn_tot = cnm.sn\n",
    "print(('Number of components:' + str(A_tot.shape[-1])))\n",
    "pl.figure()\n",
    "crd = plot_contours(A_tot, Cn, thr=0.9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>DISCARD LOW QUALITY COMPONENT </h1>\n",
    "<p> The patch dubdivision creates several spurious components that are not neurons </p>\n",
    "<p>We select the components according to criteria examining spatial and temporal components</p>\n",
    "<img src=\"docs/img/evaluationcomponent.png\"/>\n",
    "\n",
    "<p> Temporal components, for each trace: </p>\n",
    "\n",
    "<li>  compute the robust mode, corresponding to the baseline value</li>\n",
    "<li> use the values under the mode to estimate noise variance</li>\n",
    "<li> compute the probability of having large transients  given the noise distribution estimated </li>\n",
    "<li> Threshold on this probability s.t. some of the component are discarded because lacking large enough positive transients </li>\n",
    "\n",
    "<p> Spatial components, for each components: </p>\n",
    "\n",
    "<li> average the frames in the moveie where the neurons is active (from temporal component), this provides a nice image of the neuron</li>\n",
    "<li> compare this image with the corresponding spatial component (Person's correlation coefficient)</li>\n",
    "<li> threshold the correlation coefficient  </li>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% DISCARD LOW QUALITY COMPONENT\n",
    "from caiman.components_evaluation import estimate_components_quality\n",
    "final_frate = params_movie['final_frate']\n",
    "r_values_min = .7  # threshold on space consistency\n",
    "fitness_min = -40  # threshold on time variability\n",
    "# threshold on time variability (if nonsparse activity)\n",
    "fitness_delta_min = -40\n",
    "Npeaks = 10\n",
    "traces = C_tot + YrA_tot\n",
    "idx_components, idx_components_bad = estimate_components_quality(\n",
    "    traces, Y, A_tot, C_tot, b_tot, f_tot, final_frate=final_frate, \n",
    "    Npeaks=Npeaks, r_values_min=r_values_min, fitness_min=fitness_min,\n",
    "    fitness_delta_min=fitness_delta_min)\n",
    "print(('Keeping ' + str(len(idx_components)) +\n",
    "       ' and discarding  ' + str(len(idx_components_bad))))\n",
    "#%%\n",
    "pl.figure(figsize=(20,10))\n",
    "crd = plot_contours(A_tot.tocsc()[:, idx_components], Cn, thr=0.9)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> CNMF full Field of View </h1> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%%\n",
    "#A_tot = A_tot.tocsc()[:, idx_components]\n",
    "#C_tot = C_tot[idx_components]\n",
    "#%% rerun updating the components to refine\n",
    "cnm = cnmf.CNMF(n_processes, k=A_tot.shape, gSig=gSig, \n",
    "                merge_thresh=merge_thresh, p=0, dview=dview, Ain=A_tot, Cin=C_tot,\n",
    "                f_in=f_tot, rf=None, stride=None, method_deconvolution='oasis')\n",
    "cnm = cnm.fit(images)\n",
    "#%%\n",
    "A, C, b, f, YrA, sn = cnm.A, cnm.C, cnm.b, cnm.f, cnm.YrA, cnm.sn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> Discard low quality components on full field of view</h1>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% again recheck quality of components, stricter criteria\n",
    "final_frate = params_movie['final_frate']\n",
    "r_values_min = .75\n",
    "fitness_min = - 50\n",
    "fitness_delta_min = - 50\n",
    "Npeaks = 10\n",
    "traces = C + YrA\n",
    "idx_components, idx_components_bad = estimate_components_quality(\n",
    "    traces, Y, A, C, b, f, final_frate=final_frate,\n",
    "    Npeaks=Npeaks, r_values_min=r_values_min, fitness_min=fitness_min,\n",
    "    fitness_delta_min=fitness_delta_min)\n",
    "print(' ***** ')\n",
    "print((len(traces)))\n",
    "print((len(idx_components)))\n",
    "print(len(idx_components_bad))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% visualize included and excluded components\n",
    "pl.figure(figsize=(20,10))\n",
    "pl.subplot(1, 2, 1)\n",
    "crd = plot_contours(A.tocsc()[:, idx_components], Cn, thr=0.9)\n",
    "pl.subplot(1, 2, 2)\n",
    "crd = plot_contours(A.tocsc()[:, idx_components_bad], Cn, thr=0.9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# kept components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "discard_traces_fluo=nb_view_patches(Yr,A.tocsc()[:,idx_components],\n",
    "                                    C[idx_components],b,f,dims[0],\n",
    "                                    dims[1],thr = 0.8,image_neurons=Cn,\n",
    "                                    denoised_color='red')\n",
    "print('you will may be need to change the data rate to generate this one: use jupyter notebook --NotebookApp.iopub_data_rate_limit=1.0e10')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# discarded components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "discard_traces_fluo=nb_view_patches(Yr,A.tocsc()[:,idx_components_bad],\n",
    "                                    C[idx_components_bad],b,f,dims[0],dims[1],\n",
    "                                    thr = 0.8,image_neurons=Cn, denoised_color='red')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we can see here a display of the components that have been kept throught the cnmf process and the one that have been discarded. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> closing, saving, and creating denoised version </h1>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.savez(os.path.join(os.path.split(fname_new)[0], \n",
    "                      os.path.split(fname_new)[1][:-4] + 'results_analysis.npz'),\n",
    "         Cn=Cn, A=A.todense(), C=C,\n",
    "         b=b, f=f, YrA=YrA, sn=sn, d1=d1, d2=d2,\n",
    "         idx_components=idx_components, idx_components_bad=idx_components_bad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% STOP CLUSTER and clean up log files\n",
    "cm.stop_server()\n",
    "log_files = glob.glob('*_LOG_*')\n",
    "for log_file in log_files:\n",
    "    os.remove(log_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#%% reconstruct denoised movie without background\n",
    "denoised = cm.movie(A.dot(C)).reshape(dims+(-1,),order = 'F').transpose([2,0,1])\n",
    "#%%\n",
    "denoised.play(gain = 30, offset = 0,fr =100, magnification = 2)\n"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}