demos/notebooks/demo_pipeline_cnmfE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CNMF-E demo pipeline: Intro \n",
    "This notebook demonstrates how to use Caiman for processing 1p microendoscopic data. It shows how to use Caiman for the following steps:\n",
    "\n",
    "![workflow diagram](images/full_cnmfe_workflow.jpg)\n",
    "\n",
    "1. Apply the nonrigid motion correction (NoRMCorre) algorithm for motion correction to the original movie.\n",
    "2. Apply the constrained nonnegative matrix factorization endoscopic (CNMF-E) source separation algorithm to extract an initial estimate of neuronal spatial footprint and calcium traces.\n",
    "3. Apply quality control metrics to evaluate the initial estimates to narrow them down to a final set of estimates.\n",
    "\n",
    "Tools for visualization of movies and results are also included. \n",
    "\n",
    "> This demo follows a similar pattern to the CNMF demo in `demo_pipeline.ipynb`. It includes less explanation except where there are important differences. If you want to get a more explanation-heavy picture of the fundamentals, we suggest starting with `demo_pipeline.ipynb`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imports and general setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bokeh.plotting as bpl\n",
    "import cv2\n",
    "import glob\n",
    "import holoviews as hv\n",
    "from IPython import get_ipython\n",
    "import logging\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import psutil\n",
    "\n",
    "import caiman as cm\n",
    "from caiman.source_extraction import cnmf\n",
    "from caiman.source_extraction.cnmf.cnmf import load_CNMF\n",
    "from caiman.utils.utils import download_demo\n",
    "from caiman.utils.visualization import inspect_correlation_pnr, nb_inspect_correlation_pnr\n",
    "from caiman.motion_correction import MotionCorrect\n",
    "from caiman.source_extraction.cnmf import params as params\n",
    "from caiman.utils.visualization import plot_contours, nb_view_patches, nb_plot_contour\n",
    "from caiman.utils.visualization import view_quilt\n",
    "\n",
    "try:\n",
    "    if __IPYTHON__:\n",
    "        get_ipython().run_line_magic('load_ext', 'autoreload')\n",
    "        get_ipython().run_line_magic('autoreload', '2')\n",
    "        # get_ipython().run_line_magic('matplotlib', 'qt')  #uncomment to run in qt mode\n",
    "except NameError:\n",
    "    pass\n",
    "\n",
    "try:\n",
    "    cv2.setNumThreads(0)\n",
    "except:\n",
    "    pass\n",
    "\n",
    "bpl.output_notebook()\n",
    "hv.notebook_extension('bokeh')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up logger and environment variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up logging\n",
    "logfile = None # Replace with a path if you want to log to a file\n",
    "logger = logging.getLogger('caiman')\n",
    "# Set to logging.INFO if you want much output, potentially much more output\n",
    "logger.setLevel(logging.WARNING)\n",
    "logfmt = logging.Formatter('%(relativeCreated)12d [%(filename)s:%(funcName)20s():%(lineno)s] [%(process)d] %(message)s')\n",
    "if logfile is not None:\n",
    "    handler = logging.FileHandler(logfile)\n",
    "else:\n",
    "    handler = logging.StreamHandler()\n",
    "handler.setFormatter(logfmt)\n",
    "logger.addHandler(handler)\n",
    "\n",
    "# set env variables in case they weren't already set\n",
    "os.environ[\"MKL_NUM_THREADS\"] = \"1\"\n",
    "os.environ[\"OPENBLAS_NUM_THREADS\"] = \"1\"\n",
    "os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"1\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set path to file to be processed\n",
    "In this demo, we analyze the data in `data_endoscope.tif`. This movie provides sample microendoscopic data from mouse dorsal striatum where the neurons are expressing GCaMP6f. The original movie was sampled at 30Hz but spatially cropped and downsampled to 10Hz to save space for the demo -- see [the original CNMFE paper](https://elifesciences.org/articles/28728#s3) for more details about the data.\n",
    "\n",
    "The `download_demo` function will download the file and return the filepath, which will be stored in your `caiman_data` directory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "movie_path = download_demo('data_endoscope.tif')  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> When adapting this demo for your data make sure to store the full path to your movie in the `movie_path` variable. Also, for reasons we will explain below, the memory requirement of the CNMF-E algorithm are much higher compared to the CNMF algorithm. If your datasets are too large, you might need to use our online algorithms. We discuss this more in `demo_pipeline.ipynb`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load and visualize raw data\n",
    "We visualize using the built-in movie object, which is described in detail in `demo_pipeline.ipynb`. In addition to neural activity, you can also see blood flow through blood vessels in the movie. Such background activity, standard for 1p data, brings new analysis challenges. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# press q to close\n",
    "movie_orig = cm.load(movie_path) \n",
    "downsampling_ratio = 0.2  # subsample 5x\n",
    "movie_orig.resize(fz=downsampling_ratio).play(gain=0.9,\n",
    "                                              q_max=99.5, \n",
    "                                              fr=10,\n",
    "                                              plot_text=True,\n",
    "                                              magnification=4,\n",
    "                                              do_loop=True,\n",
    "                                              backend='opencv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set up multicore processing\n",
    "To enable parallel computing we will enable multiple CPUs for processing. The resulting `cluster` variable contains the pool of processors (CPUs) that will be used in later steps. If you use `dview=cluster` in later steps, then parallel processing will be used. If you use `dview=None` then no parallel processing will be used. The `num_processors_to_use` variable determines how many CPU cores you will use (when set to `None` it goes to the default of one less than the number available):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"You have {psutil.cpu_count()} CPUs available in your current environment\")\n",
    "num_processors_to_use = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have already set up multiprocessing (the `cluster` variable is already in your namespace), then that cluster will be closed and a new one created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% start a cluster for parallel processing (if a cluster already exists it will be closed and a new session will be opened)\n",
    "if 'cluster' in locals():  # 'locals' contains list of current local variables\n",
    "    print('Closing previous cluster')\n",
    "    cm.stop_server(dview=cluster)\n",
    "print(\"Setting up new cluster\")\n",
    "_, cluster, n_processes = cm.cluster.setup_cluster(backend='multiprocessing', \n",
    "                                                 n_processes=num_processors_to_use, \n",
    "                                                 ignore_preexisting=False)\n",
    "print(f\"Successfully set up cluster with {n_processes} processes\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define parameters\n",
    "We first set some parameters related to the data and motion correction and create a `params` object. We'll modify this parameter object later on with settings for source extraction. You can also set all the parameters at once as demonstrated in the `demo_pipeline.ipynb` notebook.\n",
    "\n",
    "Note here we are setting `pw_rigid` to `False`, as our data seems to mainly contain large-scale translational motion. We can always redo this later if it turns out to be a mistake."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# dataset dependent parameters\n",
    "frate = 10                       # movie frame rate\n",
    "decay_time = 0.4                 # length of a typical transient in seconds\n",
    "\n",
    "# motion correction parameters\n",
    "motion_correct = True    # flag for performing motion correction\n",
    "pw_rigid = False         # flag for performing piecewise-rigid motion correction (otherwise just rigid)\n",
    "gSig_filt = (3, 3)       # sigma for high pass spatial filter applied before motion correction, used in 1p data\n",
    "max_shifts = (5, 5)      # maximum allowed rigid shift\n",
    "strides = (48, 48)       # start a new patch for pw-rigid motion correction every x pixels\n",
    "overlaps = (24, 24)      # overlap between patches (size of patch = strides + overlaps)\n",
    "max_deviation_rigid = 3  # maximum deviation allowed for patch with respect to rigid shifts\n",
    "border_nan = 'copy'      # replicate values along the boundaries\n",
    "\n",
    "mc_dict = {\n",
    "    'fnames': movie_path,\n",
    "    'fr': frate,\n",
    "    'decay_time': decay_time,\n",
    "    'pw_rigid': pw_rigid,\n",
    "    'max_shifts': max_shifts,\n",
    "    'gSig_filt': gSig_filt,\n",
    "    'strides': strides,\n",
    "    'overlaps': overlaps,\n",
    "    'max_deviation_rigid': max_deviation_rigid,\n",
    "    'border_nan': border_nan\n",
    "}\n",
    "\n",
    "parameters = params.CNMFParams(params_dict=mc_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Motion Correction\n",
    "The background signal in micro-endoscopic data is very strong and makes motion correction challenging. As a first step the algorithm performs a high pass spatial filtering with a Gaussian kernel to remove the bulk of the lower-frequency background activity and enhance spatial landmarks. The size of the kernel is given from the parameter `gSig_filt`. If this is left to the default value of `None` then no preprocessing is performed (default option, used in 2p data for CNMF). \n",
    "\n",
    "After spatial filtering, the NoRMCorre algorithm is used to determine the motion in each frame. The inferred motion is then applied to the *original* data, not the preprocessed data, so no information is lost before source separation. The motion corrected files are saved in memory mapped format. If no motion correction is performed (i.e., `motion_correct` was set to `False`), then the file gets directly memory mapped.\n",
    "\n",
    "> For a more detailed exploration of Caiman's motion correction pipeline, see `demo_motion_correction.ipynb`. \n",
    "\n",
    "The following also plots the discovered displacements in x- and y- directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "if motion_correct:\n",
    "    # do motion correction rigid\n",
    "    mot_correct = MotionCorrect(movie_path, dview=cluster, **parameters.get_group('motion'))\n",
    "    mot_correct.motion_correct(save_movie=True)\n",
    "    fname_mc = mot_correct.fname_tot_els if pw_rigid else mot_correct.fname_tot_rig\n",
    "    if pw_rigid:\n",
    "        bord_px = np.ceil(np.maximum(np.max(np.abs(mot_correct.x_shifts_els)),\n",
    "                                     np.max(np.abs(mot_correct.y_shifts_els)))).astype(int)\n",
    "    else:\n",
    "        bord_px = np.ceil(np.max(np.abs(mot_correct.shifts_rig))).astype(int)\n",
    "        # Plot shifts\n",
    "        plt.plot(mot_correct.shifts_rig)  # % plot rigid shifts\n",
    "        plt.legend(['x shifts', 'y shifts'])\n",
    "        plt.xlabel('frames')\n",
    "        plt.ylabel('pixels')\n",
    "        plt.gcf().set_size_inches(6,3)\n",
    "\n",
    "    bord_px = 0 if border_nan == 'copy' else bord_px\n",
    "    fname_new = cm.save_memmap(fname_mc, base_name='memmap_', order='C',\n",
    "                               border_to_0=bord_px)\n",
    "else:  # if no motion correction just memory map the file\n",
    "    fname_new = cm.save_memmap(movie_path, base_name='memmap_',\n",
    "                               order='C', border_to_0=0, dview=dview)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare original (left panel) and motion corrected movie (right panel). You will notice they look similar, as there wasn't much motion to begin with. You can see from the shift plot (plotted above) that the extracted shifts were all very small."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "movie_corrected = cm.load(mot_correct.mmap_file) # load motion corrected movie\n",
    "ds_ratio = 0.2\n",
    "cm.concatenate([movie_orig.resize(1, 1, ds_ratio) - mot_correct.min_mov*mot_correct.nonneg_movie,\n",
    "                movie_corrected.resize(1, 1, ds_ratio)], \n",
    "                axis=2).play(fr=10, \n",
    "                             gain=0.9, \n",
    "                             magnification=2) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load memory mapped file\n",
    "Memory mapping is discussed in more detail in `demo_pipeline.ipynb`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load memory mappable file\n",
    "Yr, dims, T = cm.load_memmap(fname_new)\n",
    "images = Yr.T.reshape((T,) + dims, order='F')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameter setting for CNMF-E\n",
    "Everything is now set up to run source extraction with CNMFE. We will construct a new parameter dictionary and use this to modify the *existing* `parameters` object, using the `change_params()` method.\n",
    "\n",
    "There are *two* main differences between the CNMF and CNMFE source separation algorithms. The first is the background model (this is discussed in the sidebar below on the Ring Model). The second difference is in how the models are initialized. This is addressed below when we go over setting corr/pnr thresholds for initialization, which we did not have to do for our 2p data.\n",
    "\n",
    "For now, note that we have set `gnb` to `0`: this is effectively the flag telling Caiman to use CNMFE instead of CNMF. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters for source extraction and deconvolution\n",
    "p = 1               # order of the autoregressive system\n",
    "K = None            # upper bound on number of components per patch, in general None for CNMFE\n",
    "gSig = np.array([3, 3])  # expected half-width of neurons in pixels \n",
    "gSiz = 2*gSig + 1     # half-width of bounding box created around neurons during initialization\n",
    "merge_thr = .7      # merging threshold, max correlation allowed\n",
    "rf = 40             # half-size of the patches in pixels. e.g., if rf=40, patches are 80x80\n",
    "stride_cnmf = 20    # amount of overlap between the patches in pixels \n",
    "tsub = 2            # downsampling factor in time for initialization, increase if you have memory problems\n",
    "ssub = 1            # downsampling factor in space for initialization, increase if you have memory problems\n",
    "gnb = 0             # number of background components (rank) if positive, set to 0 for CNMFE\n",
    "low_rank_background = None  # None leaves background of each patch intact (use True if gnb>0)\n",
    "nb_patch = 0        # number of background components (rank) per patch (0 for CNMFE)\n",
    "min_corr = .8       # min peak value from correlation image\n",
    "min_pnr = 10        # min peak to noise ration from PNR image\n",
    "ssub_B = 2          # additional downsampling factor in space for background (increase to 2 if slow)\n",
    "ring_size_factor = 1.4  # radius of ring is gSiz*ring_size_factor\n",
    "\n",
    "parameters.change_params(params_dict={'method_init': 'corr_pnr',  # use this for 1 photon\n",
    "                                'K': K,\n",
    "                                'gSig': gSig,\n",
    "                                'gSiz': gSiz,\n",
    "                                'merge_thr': merge_thr,\n",
    "                                'p': p,\n",
    "                                'tsub': tsub,\n",
    "                                'ssub': ssub,\n",
    "                                'rf': rf,\n",
    "                                'stride': stride_cnmf,\n",
    "                                'only_init': True,    # set it to True to run CNMF-E\n",
    "                                'nb': gnb,\n",
    "                                'nb_patch': nb_patch,\n",
    "                                'method_deconvolution': 'oasis',       # could use 'cvxpy' alternatively\n",
    "                                'low_rank_background': low_rank_background,\n",
    "                                'update_background_components': True,  # sometimes setting to False improve the results\n",
    "                                'min_corr': min_corr,\n",
    "                                'min_pnr': min_pnr,\n",
    "                                'normalize_init': False,               # just leave as is\n",
    "                                'center_psf': True,                    # True for 1p\n",
    "                                'ssub_B': ssub_B,\n",
    "                                'ring_size_factor': ring_size_factor,\n",
    "                                'del_duplicates': True,                # whether to remove duplicates from initialization\n",
    "                                'border_pix': bord_px});                # number of pixels to not consider in the borders)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize the model using the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnmfe_model = cnmf.CNMF(n_processes=n_processes, \n",
    "                        dview=cluster, \n",
    "                        params=parameters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    <h2 >CNMF-E: The Ring Model</h2>\n",
    "   Background activity is very ill-behaved with 1p recordings: it often fluctuates locally and is much larger in magnitude than the neural signals we want to extract. In other words, the large-scale background model used for CNMF is not sufficient for most 1p data. Hence, Pengcheng Zhou and others came up with a localized model of background activity for CNMFE: the background at each pixel is represented as the weighted sum of activity from a circle (or ring) of pixels a certain distance from that pixel. The distance of this ring from the reference pixel is set by the <em>ring_size_factor</em> parameter. This more complex pixel-by-pixel background model explains why CNMFE is computationally more expensive than CNMF, and also why it works better to mop up large-scale localized background noise in  1p data. \n",
    "    \n",
    "<p>When you set <em>gnb</em> in the CNMF model (usually to 1 or 2), you are setting the number of global background components. The fact that you can get away with so few is testament to how well-behaved the background activity is in 2p recordings. When we set <em>gnb</em> to 0 in Caiman, this is a flag telling Caiman's back end to switch to the more complicated ring model of the background activity.</p>\n",
    "\n",
    "For more details on CNMFE you can see the <a href=\"https://elifesciences.org/articles/28728\">original paper</a> and the <a href=\"https://elifesciences.org/articles/38173\">Caiman paper</a>. \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Key parameters for CNMFE\n",
    "The key parameters for CNMFE are slightly different than for CNMF, but with some overlap. As we'll see, because of the high levels of background activity, we can't initialize the same way as with CNMF. We have two new important parameters directly related to initialization that come into play: `min_corr` and `min_pnr`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`rf` (int): *patch half-width*\n",
    "> `rf`, which stands for 'receptive field', is the half width of patches in pixels. The patch width is `2*rf + 1`. `rf` should be *at least* 3-4 times larger than the observed neuron diameter. The larger the patch size, the less parallelization will be used by Caiman. If `rf` is set to `None`, then CNMFE will be run on the entire field of view.\n",
    "\n",
    "`stride` (int): *patch overlap*\n",
    "> `stride` is the overlap between patches in pixels (the actual overlap is `stride_cnmf + 1`). This should be at least the diameter of a neuron. The larger the overlap, the greater the computational load, but the results will be more accurate when stitching together results from different patches. This param should probably have been called 'overlap' instead of 'stride'.\n",
    "\n",
    "`gSig (int, int)`: *half-width of neurons*\n",
    "> `gSig` is roughly the half-width of neurons in your movie in pixels (height, width). It is the standard deviation of the mean-centered Gaussian used to filter the movie before initialization for CNMFE. It is related to the `gSiz` parameter, which is the width of the entire kernel filter.\n",
    "\n",
    "`merge_thr (float)`: *merge threshold* \n",
    "> If the correlation between two spatially overlapping components is above `merge_thr`, they will be merged into one component. \n",
    "\n",
    "`min_corr` (float): *minimum correlation*\n",
    "> Pixels from neurons tend to be correlated with their neighbors. During initialization, Caiman filters out those pixels below `min_corr` to help select seed pixels. We discuss this more below.\n",
    "\n",
    "`min_pnr` (float): *minimum peak to noise ratio*\n",
    "> Pixels from neurons tend to have a high signal-to-noise ratio. During initialization, Caiman filters out those pixels below `min_pnr` to help select seed pixels. We discuss this more below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspect summary images and set parameters\n",
    "### Correlation-pnr plot\n",
    "For CNMFE, Caiman uses the correlation and peak-to-noise ratio (PNR) for initialization, which will both tend to be high in regions that contain neurons. Hence, we set a threshold for both quantitites to remove the low correlation/low pnr regions and seed initialization with regions most likely to contain neuronal activity. \n",
    "\n",
    "First, we calculate the correlation and pnr maps of the raw motion corrected movie after filtering with a mean-centered Gaussian with standard deviation `gSig` (for more information on this preprocessing step, see the sidebar below on CNMFE initialization). These calculation can be computationally and memory demanding for large datasets, so we subsample if there are many thousands of frames:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(gSig)\n",
    "gsig_tmp = (3,3)\n",
    "correlation_image, peak_to_noise_ratio = cm.summary_images.correlation_pnr(images[::max(T//1000, 1)], # subsample if needed\n",
    "                                                                           gSig=gsig_tmp[0], # used for filter\n",
    "                                                                           swap_dim=False) # change swap dim if output looks weird, it is a problem with tiffile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/bokeh_menu.jpg\" align=\"right\" width=\"200\"></img>\n",
    "Using `nb_inspect_correlation_pnr()`, you can inspect the correlation and PNR images to find threshold values for these quantities: `min_corr` and `min_pnr`. You can adjust the range of values displayed in the two subplots below by choosing the Y-box select tool (third button from the left -- highlighted in yellow in the accompanying image to the right) and selecting the desired region in the histograms on the right hand side of each image. You can also use the pan button (first button on the left) to zoom/adjust the axis limits in the histogram to make it easier to see the limits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nb_inspect_correlation_pnr(correlation_image, peak_to_noise_ratio, cmap='inferno') # jet, fire are also good cmaps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are looking for a couple of things in the above plot:\n",
    "1) Did we filter with a `gSig` value small enough so that we aren't blending different neurons together? To see what it is like when `gSig` is too large, set `gsig_tmp` to `(6,6)` in the above cell and then inspect the resulting correlation-pnr plots. \n",
    "2) More importantly, we want to find the threshold correlation and pnr values so that the *lower* threshold eliminates most of the noise and blood vessels from the plots, leaving behind as many of the *neural* pixels as possible. For this data it will be at a correlation value lower bound between 0.8 and 0.9, and a pnr lower bound somewhere between 10 and 20. As with CNMF, there is no perfect value. Just keep in mind it is better to have false positives later than false negatives.\n",
    "\n",
    "You can tweak the parameters in the following cell, using the `change_params()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gsig_new = gSig # unchanged\n",
    "min_corr_new  = 0.85 \n",
    "min_pnr_new = 12     \n",
    "\n",
    "cnmfe_model.params.change_params(params_dict={'gSig': gsig_new,\n",
    "                                              'min_corr': min_corr_new, \n",
    "                                              'min_pnr': min_pnr_new});"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Caiman also includes a Qt-based corr-pnr viewer that some people prefer: `inspect_correlation_pnr()`. It provides what some say is a more intuitive interface than the notebook version above. It requires you to be in Qt mode (which you can enable using the cell magic `%matplotlib qt`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    <h2>CNMFE initialization: More on correlation and peak-to-noise-ratio</h2>\n",
    "     <img src=\"images/mn_centered_gaussian.jpg\" align=\"right\" width=\"200\"></img>\n",
    "How are correlation and peak-to-noise ratio actually calculated? First Caiman convolves the motion corrected movie with a <i>mean-centered Gaussian</i> (example to the right). The sigma of the Gaussian is <em>gSig</em>, and mean centering is turned on by setting <em>center_psf</em> to <em>True</em>. This mean centering creates a Gaussian with a positive peak in the middle of width <i>approximately</i> <em>gSig/2</em>, surrounded by a negative trench, and a ledge of zeros around the outer edges. This crucial preprocessing filter serves to highlight neuronal peaks and smooth away low-frequency background activity.\n",
    "\n",
    "<p>The function <em>correlation_pnr()</em> applies this mean-centered Gaussian to each frame of the motion corrected movie and returns the correlation image of that movie, as well as the peak-to-noise-ratio (PNR). The correlation image is the correlation of each pixel with its neighbors. The PNR is the ratio of the maximum magnitude at a pixel to the noise value at that pixel (it is a fast and rough measure of signal-to-noise). As mentioned above, both of these values tend to be higher in pixels that contain neurons. The CNMFE initialization procedure is to set a threshold for both quantities, take their <i>product</i>, and use the peaks in this product map to find <b>seed pixels</b> for initialization of the CNMFE source separation algorithm.</p>\n",
    "\n",
    "More details on the initialization procedure used here can be found in the <a href=\"https://elifesciences.org/articles/28728\">CNMFE paper</a>, or by exploring the code.         \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate our spatial parameters\n",
    "As discussed in `demo_pipeline.ipynb`, the other important parameters are those used for dividing the movie into patches for parallelization of the algorithm. Namely, we want to select `rf` and `stride` parameters so that at least 3-4 neuron diameters can fit into each patch, and at least one neuron fits in the overlap region between patches. You can visualize the patches using the `view_quilt()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate stride and overlap from parameters\n",
    "cnmfe_patch_width = cnmfe_model.params.patch['rf']*2 + 1\n",
    "cnmfe_patch_overlap = cnmfe_model.params.patch['stride'] + 1\n",
    "cnmfe_patch_stride = cnmfe_patch_width - cnmfe_patch_overlap\n",
    "print(f'Patch width: {cnmfe_patch_width} , Stride: {cnmfe_patch_stride}, Overlap: {cnmfe_patch_overlap}');\n",
    "\n",
    "# plot the patches\n",
    "patch_ax = view_quilt(correlation_image, \n",
    "                      cnmfe_patch_stride, \n",
    "                      cnmfe_patch_overlap, \n",
    "                      vmin=np.percentile(np.ravel(correlation_image), 50), \n",
    "                      vmax=np.percentile(np.ravel(correlation_image), 99.5),\n",
    "                      color='yellow',\n",
    "                      figsize=(4,4));\n",
    "patch_ax.set_title(f'CNMFE Patch Width {cnmfe_patch_width}, Overlap {cnmfe_patch_overlap}');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These patches and overlaps may seem a bit large, but that is ok: our main concern is that they not be too small. If you wanted to change them, you could use `change_params()`.\n",
    "\n",
    "Now that we are happy with our parameters, let's run the algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the CNMF-E algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "cnmfe_model.fit(images);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running the algorithm creates an `estimates` class, which we discuss in detail in `demo_pipeline.ipynb`. The CNMFE `estimates` class includes almost all the same attributes as with CNMF, such as the neural spatial and temporal components `A` and `C`. \n",
    "\n",
    "It also includes the discovered model of background activity, which is different from the CNMF model. For *CNMF* the background model is returned as low-rank matrices `b` and `f`. For CNMFE, the background model parameters are represented in the matrix `W` (the weights of the *ring model* -- discussed above -- for each pixel) as well as `b0` (a constant offset for each pixel). We will show how to reconstruct the background activity below. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Component Evaluation\n",
    "Source extraction typically produces many false positives. Our next step is quality control: separating the results into \"good\" and \"bad\" neurons using two different metrics (discussed in detail in `demo_notebook.ipynb`):\n",
    "\n",
    "- **Signal-to-noise ratio (SNR)**: a minimum SNR is set for the calcium transients (`min_SNR`).\n",
    "- **Spatial correlation**:  a minimum correlation is set between the shape of each component and the frames in the movie when that component is active (`rval_thr`). \n",
    "\n",
    "> Caiman does *not* use the CNN classifier to sort neurons based on shape for 1p data: the network was trained on 2p data. Hence, we set the `use_cnn` param to `False`. \n",
    "\n",
    "Here we set the two parameters and run `evaluate_components()` to see which pass muster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_SNR = 3            # SNR threshold\n",
    "rval_thr = 0.85    # spatial correlation threshold\n",
    "\n",
    "quality_params = {'min_SNR': min_SNR,\n",
    "                  'rval_thr': rval_thr,\n",
    "                  'use_cnn': False}\n",
    "cnmfe_model.params.change_params(params_dict=quality_params)\n",
    "\n",
    "cnmfe_model.estimates.evaluate_components(images, cnmfe_model.params, dview=cluster)\n",
    "\n",
    "print('*****')\n",
    "print(f\"Total number of components: {len(cnmfe_model.estimates.C)}\")\n",
    "print(f\"Number accepted: {len(cnmfe_model.estimates.idx_components)}\")\n",
    "print(f\"Number rejected: {len(cnmfe_model.estimates.idx_components_bad)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualize results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% plot contour plots of accepted and rejected components\n",
    "cnmfe_model.estimates.plot_contours_nb(img=correlation_image, \n",
    "                                       idx=cnmfe_model.estimates.idx_components);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "View traces of accepted and rejected components:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cnmfe_model.estimates.nb_view_components(img=correlation_image, \n",
    "                                        idx=cnmfe_model.estimates.idx_components,\n",
    "                                        cmap='viridis', #gray\n",
    "                                        thr=.9); #increase to see full footprint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rejected components\n",
    "cnmfe_model.estimates.nb_view_components(img=correlation_image, \n",
    "                                        idx=cnmfe_model.estimates.idx_components_bad,\n",
    "                                        cmap='viridis', #gray\n",
    "                                        denoised_color='red',\n",
    "                                        thr=0.9); #increase to see full footprint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save and load results (optional)\n",
    "If you want to save your results so you don't have to run CNMFE again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_results = True\n",
    "if save_results:\n",
    "    save_path =  r'demo_pipeline_cnmfe_results.hdf5'  # or add full/path/to/file.hdf5\n",
    "    cnmfe_model.estimates.Cn = correlation_image # squirrel away correlation image with cnmf object\n",
    "    cnmfe_model.save(save_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To load results and pick up where you left off (note this assumes you have done preliminaries like imports and started a cluster)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "load_results = True\n",
    "if load_results:\n",
    "    save_path =  r'demo_pipeline_cnmfe_results.hdf5'  # or add full/path/to/file.hdf5\n",
    "    cnmfe_model = load_CNMF(save_path, \n",
    "                                n_processes=num_processors_to_use, \n",
    "                                dview=cluster)\n",
    "    correlation_image = cnmfe_model.estimates.Cn\n",
    "    print(f\"Successfully loaded data.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A few loose ends\n",
    "We have extracted the calcium traces C, spatial footprints A, and estimated spike counts S, which is the main goal with CNMF. But there are a few important things remaining. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract $\\Delta F/F$ values\n",
    "Currently in Caiman, we don't return a true dff value for 1p data. This is because, as mentioned in `demo_pipeline.ipynb`,  Caiman normalizes to both the baseline fluorescence and background activity. The background activity in 1p is so ill-behaved (as discussed above in the sidebar on the ring model) that Caiman currently only *detrends* the data by subtracting away the baseline but not normalizing the data. This explains the warning you will see when you run the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if cnmfe_model.estimates.F_dff is None:\n",
    "    print('Calculating estimates.F_dff')\n",
    "    cnmfe_model.estimates.detrend_df_f(quantileMin=8, \n",
    "                                      frames_window=250,\n",
    "                                      flag_auto=False,\n",
    "                                      use_residuals=False);  # use denoised data\n",
    "else:\n",
    "    print(\"estimates.F_dff already defined\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Deconvolution for 1p?\n",
    "While we haven't discussed deconvolution (the estimation of the spikes that generated the calcium traces in `C`), we suggest treating the spike counts returned for 1p data (in `estimates.S`) with CNMFE with some caution. Currently (as of early 2024) we are aware of no simultaneous ground-truth recordings (electrophysiology + 1p imaging) that could be used to evaluate deconvolution models for 1p data. There are quite a few such recordings for 2p data (for instance see [the Spikefinder paper](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1006157)).\n",
    "\n",
    "Because of this, most researchers analyze the calcium traces directly for 1p recordings (the data in `estimates.C`) or a normalized version of them, rather than `estimates.S`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## View some different movie results\n",
    "As with CNMF, the CNMFE model of the original movie is:\n",
    "\n",
    "    original_movie = neural_activity + background + residual\n",
    "    \n",
    "The main between CNMF and CNMFE is the model of the background. We can reconstruct the neural movie as `AC` just as we did in `demo_pipeline.ipynb`. Unfortunately, reconstructing the background activity via the ring model is much more complicated for CNMFE, so we will just punt to a built-in function for that in what follows, `compute_background()`, which you can explore if you are interested:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in case you are working from loaded data, recover the raw movie\n",
    "Yr, dims, num_frames = cm.load_memmap(cnmfe_model.mmap_file)\n",
    "images = np.reshape(Yr.T, [num_frames] + list(dims), order='F')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get model of neural activity and background activity (note for the neural model we are just including the accepted components):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neural_activity = cnmfe_model.estimates.A[:, cnmfe_model.estimates.idx_components] @\\\n",
    "                  cnmfe_model.estimates.C[cnmfe_model.estimates.idx_components, :]  # AC\n",
    "neural_movie = cm.movie(neural_activity).reshape(dims + (-1,), order='F').transpose([2, 0, 1])\n",
    "background_model = cnmfe_model.estimates.compute_background(Yr);  # build in function -- explore source code for details\n",
    "bg_movie = cm.movie(background_model).reshape(dims + (-1,), order='F').transpose([2, 0, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### View neural movie"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "downsampling_ratio = 0.4 \n",
    "neural_movie.resize(fz=downsampling_ratio).play(gain=1.1,\n",
    "                                            q_max=99.5, \n",
    "                                            fr=20,\n",
    "                                            plot_text=True,\n",
    "                                            magnification=2,\n",
    "                                            do_loop=False,\n",
    "                                            backend='opencv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### View background model\n",
    "You will see many regions that are constant such as blood vessels, but also lots of large-scale background flourescence, in addition to some local activity which is is on spatial scales larger than `gSig`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "downsampling_ratio = 0.8 \n",
    "bg_movie.resize(fz=downsampling_ratio).play(gain=1.1,\n",
    "                                            q_max=99.5, \n",
    "                                            fr=10,\n",
    "                                            plot_text=True,\n",
    "                                            magnification=2,\n",
    "                                            do_loop=False,\n",
    "                                            backend='opencv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### View model and residual\n",
    "We can also use the built-in `play_movie()` method to view the original movie, predicted movie, and the residual simultaneously as discussed in more detail in `demo_pipeline.ipynb`: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# without background\n",
    "cnmfe_model.estimates.play_movie(images, \n",
    "                                 q_max=99.9, \n",
    "                                 magnification=2,\n",
    "                                 include_bck=False,\n",
    "                                 gain_res=5,\n",
    "                                 use_color=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A couple of things to note about that movie:\n",
    "- The middle panel of neural activity includes *all* components (accepted and rejected), so you will see some of the blood vessel \"activity\" that was discovered and later rejected.\n",
    "- The residual includes some activity that looks neural in origin. You can try playing with different params to get them (which do you think you would try?). The [Mesmerize](https://github.com/nel-lab/mesmerize-core) package is a great way to search parameter space and visualize results with more sophisticated visualization tools if you have a tricky data set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clean up open resourses\n",
    "Shut down server, close logger."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm.stop_server(dview=cluster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Shut down logger (otherwise will not be able to delete it)\n",
    "logging.shutdown()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delete_logs = True\n",
    "logging_dir = cm.paths.get_tempdir() \n",
    "if delete_logs:\n",
    "    log_files = glob.glob(logging_dir + '\\\\demo_pipeline' + '*' + '.log')\n",
    "    for log_file in log_files:\n",
    "        print(f\"Deleting {log_file}\")\n",
    "        os.remove(log_file)\n",
    "else:\n",
    "    print(f\"If you want to inspect your logs they are in {logging_dir}\")"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}