pipeline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "# Preface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "## about this document"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Welcome to the main pipeline of minian!\n",
    "The purpose of this annotated version of the minian pipeline is to guide the user through each step of the code, working with a short demo movie.\n",
    "The intention is to enable the user to understand the code as much as possible so that they are equipped with the knowledge necessary to customize the code for their own needs, regardless of prior programming skills.\n",
    "\n",
    "Before we start, it s highly recommended that you get familiar with basic python concepts and operations like [string manipulation](https://docs.python.org/3/library/string.html), [tuples, lists and dictionaries](https://docs.python.org/3/tutorial/datastructures.html).\n",
    "\n",
    "The notes in this pipeline are supposed to give you minimal amount of knowledge to walkthrough the pipeline and touch on parameters that are most commonly tweaked.\n",
    "The [Minian readthedocs](https://minian.readthedocs.io) site contains more comprehensive documentations.\n",
    "In particular, the [API reference](https://minian.readthedocs.io/page/api/index.html) contains detailed documentation of every minian function.\n",
    "Be sure to check them out whenever you are puzzled by how to specify parameters!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "## text styling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Note on the styling of this document: most of the sentences should hopefully make sense if taken literally.\n",
    "However, some special formatting of the text is used to demonstrate the close relationship between the concepts discussed and the code, as well as encouraging the reader to understand the Python syntax.\n",
    "Specifically:\n",
    "\n",
    "-  a [hyperlink](https://en.wikipedia.org/wiki/Hyperlink) usually points to a well-defined python module, class or methods, especially when that concept is first encountered in this document.\n",
    "    The link usually points to the official documentation of that concept, which in some cases might not be the best place to start for a beginner.\n",
    "    If you find the documentation puzzling, try to google the concept in question and find a tutorial that best suits you.\n",
    "-  an inline `code` usually refers to a name that already exsists in the [namespace](https://docs.python.org/3/tutorial/classes.html#python-scopes-and-namespaces) (i.e. the context where we run the codes in this document).\n",
    "    It can be a previously encountered concept, but more often it referes to variable names or method names that we [imported](https://docs.python.org/3/reference/import.html) or have defined along the way.\n",
    "-  **bold** texts are used more loosely to highlight anything that requires more attention.\n",
    "    Though they are not used as carefully as previous formats, they often refer to specific values that a variable or method arguments can assume.\n",
    "-  <div class=\"alert alert-info\">\n",
    "    \n",
    "    Info boxes are used to provide hints and tips to help users run through this pipeline smoothly.\n",
    "    </div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "## workflow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "There are 5 main sections in this pipeline: Setting up, Pre-processing, Motion Correction, Initialization, and CNMF, which are composed of interative spatial update and temporal update.\n",
    "\n",
    "![workflow](img/workflow.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setting Up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The cells under this section should be executed every time the kernel is restarted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## load modules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Loads the minian modules, usually this cell should not be modified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import itertools as itt\n",
    "import os\n",
    "import sys\n",
    "\n",
    "import holoviews as hv\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "from dask.distributed import Client, LocalCluster\n",
    "from holoviews.operation.datashader import datashade, regrid\n",
    "from holoviews.util import Dynamic\n",
    "from IPython.core.display import display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## set path and parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Set all of the parameters that control the notebook’s behavior.\n",
    "Ideally, the following cell is the only part of the code the user will have to change when analyzing different datasets.\n",
    "Here we briefly introduce only some of the initial parameters that are necessary to start the pipeline, and leave the discussion of specific parameters for later.\n",
    "\n",
    "* `minian_path` is the path that contains the **minian** folder , where the minian codebase (.py files) reside.\n",
    "    The default value `\".\"` means “current folder”, which should work in most cases, unless you want to try out another version of minian that is not in the same folder as this notebook.\n",
    "\n",
    "* `dpath` is the folder that contains the videos to be processed.\n",
    "\n",
    "* `interactive` controls whether interactive plots will be shown for parameters exploration.\n",
    "    Interactive plotting requires CPU/memory usage, and thus could require some time (in particular, those steps where video is played).\n",
    "    In principle, the user might want to visualize interactive plots during the initial parameters exploration, once the parameters are set and ready for batch processing, the user will set interactive as False to reduce processing time.\n",
    "\n",
    "* `output_size` controls the relative size of all the plots on a scale of 0-100 percent, though it can be set to values >100 without any problem. \n",
    "\n",
    "* `param_save_minian` specifies the destination folder and format of the saved data.\n",
    "    `dpath` is the folder path  where  the data will be saved.\n",
    "    `meta_dict` is a `dictionary` that is used to construct meta data for the final labeled data structure.\n",
    "    `overwrite` is a boolean value controlling whether the data is overwritten if a file already exists.\n",
    "    We set it to `True` here so you can easily play with the demo multiple times, but **use caution** with this option during actual analysis.\n",
    "    In addition to erasing prior data that may be important to you, overwritting data may cause compatibility issues with existing data from the same minian dataset folder.\n",
    "    If you want to re-analyze a video from scratch using different parameters, it is recommended that you delete existing data first.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "<strong>folder structure</strong>\n",
    "\n",
    "The defult `meta_dict` in `param_save_minian` assumes output minian datasets are stored in heirarchiically arranged folders, as shown below:\n",
    "\n",
    "```\n",
    "mice1  \n",
    "│\n",
    "└───session1\n",
    "│   │\n",
    "│   └───minian\n",
    "│       │   Y.zarr\n",
    "│       │   A.zarr\n",
    "│       │   ...\n",
    "│   \n",
    "└───session2\n",
    "    │\n",
    "    └───minian\n",
    "```\n",
    "\n",
    "The default value can be read as follows:\n",
    "The name of the last folder (`-1`) in `dpath` (the folder that directly contains the videos) will be used to designate the value of a metadata dimension named `\"session\"`.\n",
    "The name of the second-to-last folder (`-2`) in `dpath` will be used to designate the value for `\"animal\"` and so on.\n",
    "Both the keys (name of metadata dimension) and values (numbers indicating which level of folder name should be used) of `meta_dict` can be modified to represent your preferred way of data storage. \n",
    "Note that the metadata are determined by the folder structure of saved minian datasets, not by those of input movie data.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Set up Initial Basic Parameters#\n",
    "minian_path = \".\"\n",
    "dpath = \"./demo_movies/\"\n",
    "minian_ds_path = os.path.join(dpath, \"minian\")\n",
    "intpath = \"./minian_intermediate\"\n",
    "subset = dict(frame=slice(0, None))\n",
    "subset_mc = None\n",
    "interactive = True\n",
    "output_size = 100\n",
    "n_workers = int(os.getenv(\"MINIAN_NWORKERS\", 4))\n",
    "param_save_minian = {\n",
    "    \"dpath\": minian_ds_path,\n",
    "    \"meta_dict\": dict(session=-1, animal=-2),\n",
    "    \"overwrite\": True,\n",
    "}\n",
    "\n",
    "# Pre-processing Parameters#\n",
    "param_load_videos = {\n",
    "    \"pattern\": \"msCam[0-9]+\\.avi$\",\n",
    "    \"dtype\": np.uint8,\n",
    "    \"downsample\": dict(frame=1, height=1, width=1),\n",
    "    \"downsample_strategy\": \"subset\",\n",
    "}\n",
    "param_denoise = {\"method\": \"median\", \"ksize\": 7}\n",
    "param_background_removal = {\"method\": \"tophat\", \"wnd\": 15}\n",
    "\n",
    "# Motion Correction Parameters#\n",
    "subset_mc = None\n",
    "param_estimate_motion = {\"dim\": \"frame\"}\n",
    "\n",
    "# Initialization Parameters#\n",
    "param_seeds_init = {\n",
    "    \"wnd_size\": 1000,\n",
    "    \"method\": \"rolling\",\n",
    "    \"stp_size\": 500,\n",
    "    \"max_wnd\": 15,\n",
    "    \"diff_thres\": 2,\n",
    "}\n",
    "param_pnr_refine = {\"noise_freq\": 0.06, \"thres\": 1}\n",
    "param_ks_refine = {\"sig\": 0.05}\n",
    "param_seeds_merge = {\"thres_dist\": 5, \"thres_corr\": 0.8, \"noise_freq\": 0.06}\n",
    "param_initialize = {\"thres_corr\": 0.8, \"wnd\": 10, \"noise_freq\": 0.06}\n",
    "param_init_merge = {\"thres_corr\": 0.8}\n",
    "\n",
    "# CNMF Parameters#\n",
    "param_get_noise = {\"noise_range\": (0.06, 0.5)}\n",
    "param_first_spatial = {\n",
    "    \"dl_wnd\": 5,\n",
    "    \"sparse_penal\": 0.01,\n",
    "    \"update_background\": True,\n",
    "    \"size_thres\": (25, None),\n",
    "}\n",
    "param_first_temporal = {\n",
    "    \"noise_freq\": 0.06,\n",
    "    \"sparse_penal\": 1,\n",
    "    \"p\": 1,\n",
    "    \"add_lag\": 20,\n",
    "    \"jac_thres\": 0.2,\n",
    "}\n",
    "param_first_merge = {\"thres_corr\": 0.8}\n",
    "param_second_spatial = {\n",
    "    \"dl_wnd\": 5,\n",
    "    \"sparse_penal\": 0.01,\n",
    "    \"update_background\": True,\n",
    "    \"size_thres\": (25, None),\n",
    "}\n",
    "param_second_temporal = {\n",
    "    \"noise_freq\": 0.06,\n",
    "    \"sparse_penal\": 1,\n",
    "    \"p\": 1,\n",
    "    \"add_lag\": 20,\n",
    "    \"jac_thres\": 0.4,\n",
    "}\n",
    "\n",
    "os.environ[\"OMP_NUM_THREADS\"] = \"1\"\n",
    "os.environ[\"MKL_NUM_THREADS\"] = \"1\"\n",
    "os.environ[\"OPENBLAS_NUM_THREADS\"] = \"1\"\n",
    "os.environ[\"MINIAN_INTERMEDIATE\"] = intpath"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## import minian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The following cell loads **minian** and usually should not be modified. If you encounter an `ImportError`, check that you followed the installation instructions and that `minian_path` is pointing to the right place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "%%capture\n",
    "sys.path.append(minian_path)\n",
    "from minian.cnmf import (\n",
    "    compute_AtC,\n",
    "    compute_trace,\n",
    "    get_noise_fft,\n",
    "    smooth_sig,\n",
    "    unit_merge,\n",
    "    update_spatial,\n",
    "    update_temporal,\n",
    ")\n",
    "from minian.initialization import (\n",
    "    gmm_refine,\n",
    "    initA,\n",
    "    initbf,\n",
    "    initC,\n",
    "    intensity_refine,\n",
    "    ks_refine,\n",
    "    pnr_refine,\n",
    "    seeds_init,\n",
    "    seeds_merge,\n",
    ")\n",
    "from minian.motion_correction import apply_transform, estimate_motion\n",
    "from minian.preprocessing import denoise, remove_background\n",
    "from minian.utilities import (\n",
    "    TaskAnnotation,\n",
    "    get_optimal_chk,\n",
    "    load_videos,\n",
    "    open_minian,\n",
    "    save_minian,\n",
    ")\n",
    "from minian.visualization import (\n",
    "    CNMFViewer,\n",
    "    VArrayViewer,\n",
    "    generate_videos,\n",
    "    visualize_gmm_fit,\n",
    "    visualize_motion,\n",
    "    visualize_preprocess,\n",
    "    visualize_seeds,\n",
    "    visualize_spatial_update,\n",
    "    visualize_temporal_update,\n",
    "    write_video,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## module initialization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The following cell handles initialization of modules and parameters necessary for minian to be run and usually should not be modified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "dpath = os.path.abspath(dpath)\n",
    "hv.notebook_extension(\"bokeh\", width=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## start cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster = LocalCluster(\n",
    "    n_workers=n_workers,\n",
    "    memory_limit=\"2GB\",\n",
    "    resources={\"MEM\": 1},\n",
    "    threads_per_worker=2,\n",
    "    dashboard_address=\":8787\",\n",
    ")\n",
    "annt_plugin = TaskAnnotation()\n",
    "cluster.scheduler.add_plugin(annt_plugin)\n",
    "client = Client(cluster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pre-processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "In the pre-processing steps that follow, the pipeline will load and process the videos  (downsampling, subsetting, denoising).\n",
    "\n",
    "All functions are evaluated lazily, which means that initially only a “plan” for the actual computation will be created, without its execution. Actual computations are carried out only when results are being saved"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## loading videos and visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Recall the values of `param_load_videos`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_load_videos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The first argument of load_videos should be the path that contains the videos(`dpath`).\n",
    "We then pass the dictionary, `param_load_videos`, defined earlier, which specifies several relevant arguments.\n",
    "The argument `pattern` is optional and is the [regular expression](https://docs.python.org/3/library/re.html) used to filter files under the specified folder.\n",
    "The default value `r\"msCam[0-9]+\\.avi$\"` means that a file can only be loaded if its filename contains **'msCam'**, followed by at least one number, then **'.avi'** as the end of the filename.\n",
    "This can be changed to suit the naming convention of your videos.\n",
    "The resulting \"video array\" `varr` contains three dimensions: `height`, `width`, and `frame`.\n",
    "If you wish to downsample the video, pass in a dictionary to `downsample`, with the name of dimensions as keys and  the downsampling folds as integer value.\n",
    "For example, `downsample=dict(\"frame\"=2)` will temporally downsample the video with a factor of 2.\n",
    "`downsample_strategy` will assume two values: either `\"subset\"`, meaning downsampling are carried out simply by subsetting the data, or `\"mean\"`, meaning a mean will be calculated on the window of downsampling (the latter being slower).\n",
    "\n",
    "In addition to the video array `varr`, the following cell also try to estimate best chunk size `chk` to use for computations.\n",
    "This variable is needed for later steps since it's important to keep chunk size consistent within the pipeline.\n",
    "If for some reason you have to restart the kernel at some point, remember to either note down the content of `chk` or rerun the following cell.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "<strong>changing parameters</strong>\n",
    "\n",
    "All minian parameters are `dict` and you can freely change them in various ways.\n",
    "You can go back to the initial parameter setting cell and change things there.\n",
    "Alternatively, you can add a code cell before running the relevant step.\n",
    "For example, the following line will tell the function to load from `\"/my/data_path\"`:\n",
    "```python\n",
    "param_load_videos[\"vapth\"] = \"/my/data_path\"\n",
    "```\n",
    "While the following line will change the downsample setting (which is specified as a `dict` on its own) when loading the video:\n",
    "```python\n",
    "param_load_videos[\"downsample\"] = {\"frame\": 2}\n",
    "```\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "varr = load_videos(dpath, **param_load_videos)\n",
    "chk, _ = get_optimal_chk(varr, dtype=float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "We then immediately save the array representation to the intermediate folder to avoid repeatedly loading the video in later steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "varr = save_minian(\n",
    "    varr.chunk({\"frame\": chk[\"frame\"], \"height\": -1, \"width\": -1}).rename(\"varr\"),\n",
    "    intpath,\n",
    "    overwrite=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The variable `varr` is a [xarray.DataArray](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.html#xarray.DataArray).\n",
    "Now is a perfect time to familiarize yourself with this data structure and the [xarray](https://xarray.pydata.org/en/stable/) module in general, since we will be using these data structures throughout the analysis.\n",
    "Basically, a `xarray.DataArray` is N-dimensional array labeled with additional metadata, with many useful properties that make them easy to manipulate.\n",
    "We can ask the computer to print out some information of `varr` by calling it (as with any other variable):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "varr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "We can see now that `varr` is a `xarray.DataArray` with a [name](https://xarray.pydata.org/en/stable/generated/xarray.DataArray.name.html#xarray.DataArray.name) `\"demo_movies\"` and three dimensions: `frame`, `height` and `width`.\n",
    "Each dimension is labeled with ascending natural numbers.\n",
    "The [dtype](https://xarray.pydata.org/en/stable/generated/xarray.DataArray.dtype.html#xarray.DataArray.dtype) ([data type](https://docs.scipy.org/doc/numpy-1.14.0/user/basics.types.html)) of `varr` is `numpy.uint8`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## visualize raw data and optionally set roi for motion correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Once the data is loaded we can visualize the content of `varr` with the help of `VArrayViewer`, which shows the array as a movie.\n",
    "You can also plot summary traces like mean fluorescence across `frames` by passing a `list` of names of traces as inputs.\n",
    "Currently `\"mean\"`, `\"min\"`, `\"max\"` and `\"diff\"` are supported, where `\"diff\"` is mean fluorescent value difference across all pixels in a frame.\n",
    "\n",
    "`VArrayViewer` also supports a box drawing tool where you can draw a box in the field of view (FOV) and save it as a mask using the `“save mask”` button.\n",
    "The mask is saved as `vaviewer.mask`, and can be retrieved and used at later stages, for example, when you want to run motion correction on a sub-region of the FOV.\n",
    "See the [API reference](https://minian.readthedocs.io/page/api/minian.visualization.html#minian-visualization-VArrayViewer) for more detail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=output_size)\n",
    "if interactive:\n",
    "    vaviewer = VArrayViewer(varr, framerate=5, summary=[\"mean\", \"max\"])\n",
    "    display(vaviewer.show())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "If you decide to set a mask for motion correction, the following cell is an example of how to convert the mask into a `subset_mc` parameter that can be later passed into motion correction functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if interactive:\n",
    "    try:\n",
    "        subset_mc = list(vaviewer.mask.values())[0]\n",
    "    except IndexError:\n",
    "        pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## subset part of video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Before proceeding to pre-processing, it’s good practice to check if there is anything obviously wrong with the video (e.g. the camera suddenly dropped, resulting in dark frames).\n",
    "This can usually be observed by visualizing the video and plotting the timecourse of the mean fluorescence.\n",
    "We can utilize the [xarray.DataArray.sel](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.sel.html) method and [slice](https://docs.python.org/3/library/functions.html#slice) to subset any part of the data we want.\n",
    "By default `subset = None` will result in no subsetting.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "<strong>subsetting data</strong>\n",
    "\n",
    "The [xarray.DataArray.sel](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.sel.html) method takes in either a `dict` or keyword arguments.\n",
    "In both cases you want to specify the dimension names and the coordinates of the subset as key-value pairs.\n",
    "For example, say you want only the first 800 frames of the video, the following two lines will both work and they are equivalent:\n",
    "```python\n",
    "varr.sel(frame=slice(0, 799)) # slice object is inclusive on both ends\n",
    "varr.sel({\"frame\": slice(0, 799)})\n",
    "```\n",
    "This also works on arbitrary dimensions.\n",
    "For example, the following will give you a 100px x 100px chunk of your movie at a corner:\n",
    "```python\n",
    "varr.sel(height=slice(0, 99), width=slice(0, 99))\n",
    "```\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varr_ref = varr.sel(subset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## glow removal and visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Here we remove the general glow background caused by vignetting effect.\n",
    "We simply calculate a minimum projection across all `frame`s and subtract that projection from all `frame`s.\n",
    "A benefit of doing this is you could interpret the result as \"change of fluorescence from baseline\", while preserving the linear scale of the raw data, which is usually the range of a 8-bit integer -- 0-255.\n",
    "The result can be visualized again with `VArrayViewer`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "varr_min = varr_ref.min(\"frame\").compute()\n",
    "varr_ref = varr_ref - varr_min"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.7))\n",
    "if interactive:\n",
    "    vaviewer = VArrayViewer(\n",
    "        [varr.rename(\"original\"), varr_ref.rename(\"glow_removed\")],\n",
    "        framerate=5,\n",
    "        summary=None,\n",
    "        layout=True,\n",
    "    )\n",
    "    display(vaviewer.show())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## denoise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Recall that by default we use a median filter for denoising:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_denoise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "There is only one parameter controlling how the filtering is done: the kernel size (`ksize`) of the filter.\n",
    "The effect of this parameter can be visualized with the tool below.\n",
    "Alternatively other methods like gaussian filter, anisotropic filter etc. are implmented as well.\n",
    "See the [API reference](https://minian.readthedocs.io/page/api/minian.preprocessing.html#minian-preprocessing-denoise) for `denoise` for more detail.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "Generally `ksize=5` is good (approximately half the diamater of the largest cell).\n",
    "Note that if you do want to play with the ksize, it has to be odd number.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "if interactive:\n",
    "    display(\n",
    "        visualize_preprocess(\n",
    "            varr_ref.isel(frame=0).compute(),\n",
    "            denoise,\n",
    "            method=[\"median\"],\n",
    "            ksize=[5, 7, 9],\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell would carry out denoise step.\n",
    "Be sure to [change the parameters](https://minian.readthedocs.io/page/start_guide/faq.html#i-don-t-know-python-can-i-still-use-the-pipeline) based on visualization results before running the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varr_ref = denoise(varr_ref, **param_denoise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## background removal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Recall the parameters for background removal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_background_removal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "This step attempts to estimate background (everything except the fluorescent signal of in-focus cells) and subtracts it from the frame.\n",
    "By default we use a morphological tophat operation to estimate the background from each frame:\n",
    "First, a [disk element](http://scikit-image.org/docs/dev/api/skimage.morphology.html#disk) with a radius of `wnd` is created.\n",
    "Then, a [morphological erosion](https://homepages.inf.ed.ac.uk/rbf/HIPR2/erode.htm) using the disk element is applied to each frame, which eats away any bright \"features\" that are smaller than the disk element.\n",
    "Subsequently, a [morphological dilation](https://homepages.inf.ed.ac.uk/rbf/HIPR2/dilate.htm) is applied to the \"eroded\" image, which in theory undoes the erosion except the bright \"features\" that were completely eaten away.\n",
    "The overall effect of this process is to remove any bright feature that is smaller than a disk with radius `wnd`.\n",
    "Thus, when setting `wnd` to the **largest** expected radius of cell, this process can give us a good estimation of the background.\n",
    "Then finally the estimated background is subtracted from each frame.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "Pragmatically `wnd=15` works well.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "if interactive:\n",
    "    display(\n",
    "        visualize_preprocess(\n",
    "            varr_ref.isel(frame=0).compute(),\n",
    "            remove_background,\n",
    "            method=[\"tophat\"],\n",
    "            wnd=[10, 15, 20],\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell would carry out background removal step.\n",
    "Be sure to [change the parameters](https://minian.readthedocs.io/page/start_guide/faq.html#i-don-t-know-python-can-i-still-use-the-pipeline) based on visualization results before running the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varr_ref = remove_background(varr_ref, **param_background_removal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## save result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Here we are saving our pre-processed video (`varr_ref`) into the intermediate folder (`intpath`).\n",
    "Note that for every saved variable a separate folder will be created based on the `.name` attribute of that variable.\n",
    "And variables with the same `.name` attribute will be saved to same folder regardless the variable name, potentially overwritting each other!\n",
    "Here we [rename](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.rename.html) it to `\"varr_ref\"` so that the saved folder will be named \"varr_ref.zarr\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "varr_ref = save_minian(varr_ref.rename(\"varr_ref\"), dpath=intpath, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Motion Correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## estimate motion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Recall the parameters for `estimate_motion`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_estimate_motion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "By default the motion estimation process is simple: for each pair of frames it calculates a [phase correlation](https://en.wikipedia.org/wiki/Phase_correlation) between the two frames using [fft](https://en.wikipedia.org/wiki/Fast_Fourier_transform).\n",
    "Then the peak of the phase correlation will correspond to the translational shift between the two frames.\n",
    "The argument `dim` specifies along which dimension to run the motion estimation, and should always be set to `\"frame\"` here.\n",
    "Usually this step is parameter-free, but see [API reference](https://minian.readthedocs.io/page/api/minian.motion_correction.html#minian-motion_correction-estimate_motion) for `estimate_motion` for more advanced tweaking of the parameters.\n",
    "By default, the results from `estimate_motion` are saved in a two dimensional `DataArray` called `motion`, with two labels on the `shift_dim` dimension, representing the shifts along `\"height\"` and `\"width\"` directions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "motion = estimate_motion(varr_ref.sel(subset_mc), **param_estimate_motion)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## save motion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Here is the first time we save a variable to the final output folder using `param_save_minian`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_save_minian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "As mentioned before `param_save_minian` decides how your data will be saved and what metadata will be stored.\n",
    "Additionally we use the `chk` variable earlier to make sure all our data have same chunk size along same dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "motion = save_minian(\n",
    "    motion.rename(\"motion\").chunk({\"frame\": chk[\"frame\"]}), **param_save_minian\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## visualization of motion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Here we visualize `motion` as a fluctuating curve across `frame`s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=output_size)\n",
    "visualize_motion(motion)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## apply transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "After determining each frame's motion, we use the function `apply_transform` to correct for the motion.\n",
    "Notably, we have to decide what to do with pixels that are shifted from outside of the FOV.\n",
    "The default is to fill them with 0. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Y = apply_transform(varr_ref, motion, fill=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## save result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Here we save two versions of the motion-corrected movie `Y`.\n",
    "Their contents are identical.\n",
    "The only difference is how they are chunked.\n",
    "Also note that we convert the data to `float` type for better downstream processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "Y_fm_chk = save_minian(Y.astype(float).rename(\"Y_fm_chk\"), intpath, overwrite=True)\n",
    "Y_hw_chk = save_minian(\n",
    "    Y_fm_chk.rename(\"Y_hw_chk\"),\n",
    "    intpath,\n",
    "    overwrite=True,\n",
    "    chunks={\"frame\": -1, \"height\": chk[\"height\"], \"width\": chk[\"width\"]},\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## visualization of motion-correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Here we visualize the result of motion correction `Y_fm_chk` side by side with input `varr_ref`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.7))\n",
    "if interactive:\n",
    "    vaviewer = VArrayViewer(\n",
    "        [varr_ref.rename(\"before_mc\"), Y_fm_chk.rename(\"after_mc\")],\n",
    "        framerate=5,\n",
    "        summary=None,\n",
    "        layout=True,\n",
    "    )\n",
    "    display(vaviewer.show())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "A potentially better visualization is to look at max projcetion across all frames before and after the motion correction.\n",
    "You should see that the cells after motion correction have much more defined borders."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im_opts = dict(\n",
    "    frame_width=500,\n",
    "    aspect=varr_ref.sizes[\"width\"] / varr_ref.sizes[\"height\"],\n",
    "    cmap=\"Viridis\",\n",
    "    colorbar=True,\n",
    ")\n",
    "(\n",
    "    regrid(\n",
    "        hv.Image(\n",
    "            varr_ref.max(\"frame\").compute().astype(np.float32),\n",
    "            [\"width\", \"height\"],\n",
    "            label=\"before_mc\",\n",
    "        ).opts(**im_opts)\n",
    "    )\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            Y_hw_chk.max(\"frame\").compute().astype(np.float32),\n",
    "            [\"width\", \"height\"],\n",
    "            label=\"after_mc\",\n",
    "        ).opts(**im_opts)\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## generate video for motion-correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Finally, we can generate a mp4 video for the original movie and motion-corrected movie side-by-side so that we can play it back at original speed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "vid_arr = xr.concat([varr_ref, Y_fm_chk], \"width\").chunk({\"width\": -1})\n",
    "write_video(vid_arr, \"minian_mc.mp4\", dpath)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Initialization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "In order to run CNMF, we first need to generate an initial estimate of where our cells are likely to be and what their temporal activity is likely to look like.\n",
    "The whole initialization section is adapted from the [MIN1PIPE](https://github.com/JinghaoLu/MIN1PIPE) package.\n",
    "See the [paper](https://www.cell.com/cell-reports/fulltext/S2211-1247(18)30826-X) for full details about the theory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## compute max projection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Here we calculate a max projection that will be used later.\n",
    "We also save this max projection in the output data folder, since it will be useful when carrying out cross-session registrations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_proj = save_minian(\n",
    "    Y_fm_chk.max(\"frame\").rename(\"max_proj\"), **param_save_minian\n",
    ").compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## generating over-complete set of seeds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "We first generate an over-complete set of **seeds**.\n",
    "Recall the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_seeds_init"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The idea is that we select some subset of frames, compute a max projection of those frames, and find the local maxima of that max projection.\n",
    "We keep repeating this process and we collect all the local maxima until we obtain an overly-complete set of local maximas, which are the potential locations of cells, which we call **seeds**.\n",
    "The assumption here is that the center of cells are brighter than their surroundings on some, but not necessarily all, frames.\n",
    "There are several parameters controlling how we subset the frames:\n",
    "By default we use `method=\"rolling\"`, which use a rolling window across time to chunk and compute max projections.\n",
    "`wnd_size` controls the number of frames in each chunk.\n",
    "`stp_size` is the distance between the center of each chunk.\n",
    "For example, if `wnd_size=100` and `stp_size=50`, the windows will be as follows: (0, 100), (50, 150), (100, 200)...\n",
    "Alternatively you can use `method=\"random\"` to use random sampling of frames instead of rolling window.\n",
    "See the [API reference](https://minian.readthedocs.io/page/api/minian.initialization.html#minian-initialization-seeds_init) of `seeds_init` for details.\n",
    "Additionally we have two parameters controlling how the local maxima are found.\n",
    "`max_wnd` controls the window size within which a single pixel will be choosen as local maxima.\n",
    "In order to capture cells with all sizes, we actually find local maximas with different window size and merge all of them, starting from 2 all the way up to `max_wnd`.\n",
    "Hence `max_wnd` should be the radius of the **largest** cell you want to detect.\n",
    "Finally in order to get rid of local maxima with very little fluctuation, we set a `diff_thres` which is the minimal fluorescent diffrence of a seed across `frame`s.\n",
    "Since the linear scale of the raw data is preserved, we can set this threshold empirically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "seeds = seeds_init(Y_fm_chk, **param_seeds_init)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The `seeds` variable is a [pandas.DataFrame](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.html), with each row representing a seed.\n",
    "The column \"height\" and \"width\" defines the location of the seed.\n",
    "The column \"seeds\" is the number of chunks where the particular seed/pixel is considered a local maxima."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "seeds.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "We can visualize the seeds as points overlaid on top of the `max_proj` image. Each white dot is a seed and could potentially be the location of a cell. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=output_size)\n",
    "visualize_seeds(max_proj, seeds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## peak-noise-ratio refine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Next, we refine seeds based upon their temporal activity.\n",
    "This requires that we separate our signal from noise based upon frequency.\n",
    "To find the best cut-off frequency, we are going to take a few example seeds and separate their activity with different cut-off frequencies.\n",
    "We will then view the results and select a frequency which we believe best separates signal from noise.\n",
    "\n",
    "This is a complicated part in the pipeline, but it is important to understand since the idea is used later and it will allow you to do parameter exploration based on your need.\n",
    "The basic idea is to run some function on a small subset of the data using different parameters within a for loop, and then visualizing the results using `holoviews`.\n",
    "\n",
    "We will go line by line:\n",
    "\n",
    "1. First we create a `list` of frequencies we want to test -- `noise_freq_list`.\n",
    "    The \"frequency\" values here are a fraction of your **sampling rate**.\n",
    "    Note that if you have temporally downsampled, the fraction here is relative to the downsampled rate.\n",
    "1. Then we randomly select 6 rows from `seeds` and call them `example_seeds`.\n",
    "1. Then we pull out the corresponding temporal activities of the 6 selected seeds and assign it to `example_trace`.\n",
    "1. We then create an empty dictionary `smooth_dict` to store the resulting visualizations.\n",
    "1. After initializing these variables, we use a `for` loop to iterate through `noise_freq_list`, with one of the values from the list as `freq` during each iteration.\n",
    "1. Within the loop, we run `smooth_sig` twice on `example_trace` with the current `freq` we are testing out.\n",
    "    The low-passed result is assigned to `trace_smth_low`, while the high-pass result is assigned to `trace_smth_high`.\n",
    "1. Then we make sure to actually carry-out the computation by calling the `compute` method on the resulting `DataArray`s.\n",
    "1. Next, we turn the two traces into visualizations:\n",
    "    we construct interactive line plots ([hv.Curve](http://holoviews.org/reference/elements/bokeh/Curve.html)s) from them and put them in a container called a [hv.HoloMap](http://holoviews.org/reference/containers/bokeh/HoloMap.html).\n",
    "1. Finally, we store the whole visualization in `smooth_dict`, with the keys being the `freq` and values corresponding to the result of this iteration.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "So if you want to customize the parameter values to explore, in this case the cut-off frequency to test, you can edit `noise_freq_list`.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "if interactive:\n",
    "    noise_freq_list = [0.005, 0.01, 0.02, 0.06, 0.1, 0.2, 0.3, 0.45, 0.6, 0.8]\n",
    "    example_seeds = seeds.sample(6, axis=\"rows\")\n",
    "    example_trace = Y_hw_chk.sel(\n",
    "        height=example_seeds[\"height\"].to_xarray(),\n",
    "        width=example_seeds[\"width\"].to_xarray(),\n",
    "    ).rename(**{\"index\": \"seed\"})\n",
    "    smooth_dict = dict()\n",
    "    for freq in noise_freq_list:\n",
    "        trace_smth_low = smooth_sig(example_trace, freq)\n",
    "        trace_smth_high = smooth_sig(example_trace, freq, btype=\"high\")\n",
    "        trace_smth_low = trace_smth_low.compute()\n",
    "        trace_smth_high = trace_smth_high.compute()\n",
    "        hv_trace = hv.HoloMap(\n",
    "            {\n",
    "                \"signal\": (\n",
    "                    hv.Dataset(trace_smth_low)\n",
    "                    .to(hv.Curve, kdims=[\"frame\"])\n",
    "                    .opts(frame_width=300, aspect=2, ylabel=\"Signal (A.U.)\")\n",
    "                ),\n",
    "                \"noise\": (\n",
    "                    hv.Dataset(trace_smth_high)\n",
    "                    .to(hv.Curve, kdims=[\"frame\"])\n",
    "                    .opts(frame_width=300, aspect=2, ylabel=\"Signal (A.U.)\")\n",
    "                ),\n",
    "            },\n",
    "            kdims=\"trace\",\n",
    "        ).collate()\n",
    "        smooth_dict[freq] = hv_trace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "At the end of the process, we put together a `smooth_dict`.\n",
    "Here we convert that into an interactive plot, from which we can determine the frequency that best separates noise and signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.7))\n",
    "if interactive:\n",
    "    hv_res = (\n",
    "        hv.HoloMap(smooth_dict, kdims=[\"noise_freq\"])\n",
    "        .collate()\n",
    "        .opts(aspect=2)\n",
    "        .overlay(\"trace\")\n",
    "        .layout(\"seed\")\n",
    "        .cols(3)\n",
    "    )\n",
    "    display(hv_res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<strong>picking noise frequency</strong>\n",
    "\n",
    "We can now use the interactive visualization to pick the best cut-off frequency.\n",
    "Here is an example of what you might see:\n",
    "\n",
    "<div stype=\"clear:both\"><img src=\"img/param_pnr.png\" style=\"width: 50%\"/></div>\n",
    "\n",
    "We are looking for the frequency that can best seperate real signal from noise.\n",
    "Hence, `noise_freq=0.005` in the example is not ideal, since real calcium activities are overly smoothed as well.\n",
    "At the same time, `noise_freq=0.45` is not ideal either, since a lot of high-frequency noise are visible in the signal.\n",
    "Hence, `noise_freq=0.05` in the middle is a good choice in this example.\n",
    "Now, say you already found your parameters, it's time now to pass them in! Either go back to initial parameters setting step and modify them there, or call the parameter here and change its value/s accordingly.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Having determined the frequency that best separates signal from noise, we move on to peak-noise-ratio refining.\n",
    "Recall the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_pnr_refine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "First we filter the temporal activities for each seed using the `noise_freq` we choose.\n",
    "**signal** is defined as the low-pass filtered temporal signal, while **noise** is high-pass filtered signal.\n",
    "Then we compute the peak-to-peak value (max minus min) for both the **real** signal and **noise** signal.\n",
    "The peak-noise-ratio is defined as the ratio between the peak-to-peak value of **signal** and that of **noise**.\n",
    "We then threshold the seeds based on this peak-noise-ratio, with the assumption that temporal activities from real cells should have higher fluctuation in the low-frequency range and lower fluctuation in the high-frequency range.\n",
    "`thres` is the threshold for peak-noise-ratios.\n",
    "Pragmatically `thres=1` works fine and makes sense.\n",
    "You can also use `thres=\"auto\"`, where a gaussian mixture model with 2 components will be run on the peak-noise-ratios and seeds will be selected if they belong to the \"higher\" gaussian."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell would carry out pnr refinement step.\n",
    "Be sure to [change](https://minian.readthedocs.io/page/start_guide/faq.html#i-don-t-know-python-can-i-still-use-the-pipeline) the `noise_freq` based on visualization results before running the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "seeds, pnr, gmm = pnr_refine(Y_hw_chk, seeds, **param_pnr_refine)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The resulting `seeds` dataframe will have an additional column \"mask_pnr\" showing you whether a seed passed the peak-noise-ratio threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "seeds.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "If you chose `thres=\"auto\"` you can use the following cell to visualize the gaussian mixture model fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if gmm:\n",
    "    display(visualize_gmm_fit(pnr, gmm, 100))\n",
    "else:\n",
    "    print(\"nothing to show\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "We can visualize seeds on top of the max projection again.\n",
    "This time we can also pass in \"mask_pnr\" so that the visualization can show us which seeds are filtered out in `pnr_refine`.\n",
    "White dots are accepted seeds and red dots are taken out.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "If you see seeds being filtered out that you believe should be cells, either skip this step or try lower the threshold a bit.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=output_size)\n",
    "visualize_seeds(max_proj, seeds, \"mask_pnr\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ks refine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "`ks_refine` refines the seeds using [Kolmogorov-Smirnov test](https://en.wikipedia.org/wiki/Kolmogorov–Smirnov_test).\n",
    "Recall the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_ks_refine"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The idea is simple: if a seed corresponds to a cell, its fluorescence intensity across frames should be distributed following a [bimodal](https://en.wikipedia.org/wiki/Multimodal_distribution) distribution, with a large peak representing noise/little activity, and another peak representing when the seed/cell is active.\n",
    "Thus, we can carry out KS test on the distribution of fluorescence intensity of each seed, and keep only the seeds where the null hypothesis (that the fluorescence is just a normal distribution) is rejected.\n",
    "`sig` is the p value at which the null hypothesis is rejected. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "seeds = ks_refine(Y_hw_chk, seeds, **param_ks_refine)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "We can visualize seeds on top of the max projection.\n",
    "White dots are accepted seeds and red dots are taken out.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "The ks refine step tend to have minimal impact.\n",
    "It's quite normal to see no seeds being filtered out and you can proceed in such cases.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=output_size)\n",
    "visualize_seeds(max_proj, seeds, \"mask_ks\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## merge seeds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "At this point, much of our refined seeds likely reflect the position of an actual cell.\n",
    "However, we are likely to still have multiple seeds per cell, which we want to avoid.\n",
    "Here we discard redundant seeds through a process of merging.\n",
    "Recall the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_seeds_merge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The function `seeds_merge` attempts to merge seeds together which potentially come from the same cell, based upon their spatial distance and temporal correlation.\n",
    "Specifically, `thres_dist` is the threshold for euclidean distance between pairs of seeds, in pixels, and `thres_corr` is the threshold for pearson correlation between pairs of seeds.\n",
    "In addition, it's very beneficial to smooth the signals before computing the correlation, and `noise_freq` determines how smoothing should be done.\n",
    "Any pair of seeds that are within `thres_dist` **and** has a correlation higher than `thres_corr` will be merged together, such that only the seed with maximum intensity in the max projection of the video will be kept.\n",
    "\n",
    "Note that we apply the masking from previous refining steps here and use only those seeds as input to `seeds_merge`.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "Since the merge is carried out transitively, in practice it is beneficial to keep `thres_dist` small to avoid over-merging, especially when the density of cells is high.\n",
    "The goal of this step is to reduce the number of duplicated seeds as much as possible, but not over-merge the seeds, since merging is irreversible.\n",
    "Hence it is OK to have multiple seeds per cell at this stage.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "seeds_final = seeds[seeds[\"mask_ks\"] & seeds[\"mask_pnr\"]].reset_index(drop=True)\n",
    "seeds_final = seeds_merge(Y_hw_chk, max_proj, seeds_final, **param_seeds_merge)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "We visualize the result on top of the max projection.\n",
    "The red dots here indicate seeds that has been merged to nearby seeds (those shown in white)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=output_size)\n",
    "visualize_seeds(max_proj, seeds_final, \"mask_mrg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## initialize spatial matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Up till now, the seeds we have are only one-pixel dots.\n",
    "In order to kick start CNMF we need something more like the spatial footprint (`A`) and temporal activities (`C`) of real cells.\n",
    "Thus we need to `initilalize` `A` and `C` from the seeds we have.\n",
    "Recall the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_initialize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "To obtain the initial spatial matrix `A`, for each seed, we calculate Pearson correlation between the seed and surrounding pixels.\n",
    "Calculating correlation with all other pixels for every seed is time-consuming and unnecessary.\n",
    "Hence we use `wnd` to control the window size for calculating the correlation, and thus is the maximum possible size of any spatial footprint in the initial spatial matrix.\n",
    "At the same time we do not want pixels with low correlation value to influence our estimation of temporal signals, thus a `thres_corr` is also implemented where only pixels with correlation above this threshold are kept."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "A_init = initA(Y_hw_chk, seeds_final[seeds_final[\"mask_mrg\"]], **param_initialize)\n",
    "A_init = save_minian(A_init.rename(\"A_init\"), intpath, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## initialize temporal matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "After generating `A`, for each seed, we calculate a weighted average of pixels around the seed, where the weight are the initial spatial footprints in `A` we just generated.\n",
    "We use this weighted average as the initial estimation of temporal activities for each units in `C`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "C_init = initC(Y_fm_chk, A_init)\n",
    "C_init = save_minian(\n",
    "    C_init.rename(\"C_init\"), intpath, overwrite=True, chunks={\"unit_id\": 1, \"frame\": -1}\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## merge units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "We want to do another pass of merging before proceeding, and this time we want to merge based on initialized `A` and `C`.\n",
    "Recall the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_init_merge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The `unit_merge` consider all cells whose spatial footprints have at least one pixel in common to be potential target of merging.\n",
    "Hence, the only parameter you need to specify is a threshold for the correlation of the temporal activities of cells to be merged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "A, C = unit_merge(A_init, C_init, **param_init_merge)\n",
    "A = save_minian(A.rename(\"A\"), intpath, overwrite=True)\n",
    "C = save_minian(C.rename(\"C\"), intpath, overwrite=True)\n",
    "C_chk = save_minian(\n",
    "    C.rename(\"C_chk\"),\n",
    "    intpath,\n",
    "    overwrite=True,\n",
    "    chunks={\"unit_id\": -1, \"frame\": chk[\"frame\"]},\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## initialize background terms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Finally, we need two more terms: `b` and `f`, representing the spatial footprint and temporal dynamics of the background, respectively.\n",
    "We first compute an estimation of cellular activities by taking the outer product of `A` and `C`, resulting in an video array with dimesion `height`, `width` and `frame`.\n",
    "We then subtract this array from `Y_fm_chk`, resulting in a \"residule\" movie.\n",
    "Then `b` is estimated as the mean projection of the \"residule\" movie, while `f` is estimated as the mean fluorescence fluctuation across time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "b, f = initbf(Y_fm_chk, A, C_chk)\n",
    "f = save_minian(f.rename(\"f\"), intpath, overwrite=True)\n",
    "b = save_minian(b.rename(\"b\"), intpath, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## visualization of initialization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Finally we visualize the result of our initialization by plotting a projection of the spatial matrix `A`, a raster of the temporal matrix `C`, as well as background terms `b` and `f`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.55))\n",
    "im_opts = dict(\n",
    "    frame_width=500,\n",
    "    aspect=A.sizes[\"width\"] / A.sizes[\"height\"],\n",
    "    cmap=\"Viridis\",\n",
    "    colorbar=True,\n",
    ")\n",
    "cr_opts = dict(frame_width=750, aspect=1.5 * A.sizes[\"width\"] / A.sizes[\"height\"])\n",
    "(\n",
    "    regrid(\n",
    "        hv.Image(\n",
    "            A.sum(\"unit_id\").rename(\"A\").compute().astype(np.float32),\n",
    "            kdims=[\"width\", \"height\"],\n",
    "        ).opts(**im_opts)\n",
    "    ).relabel(\"Initial Spatial Footprints\")\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            C.rename(\"C\").compute().astype(np.float32), kdims=[\"frame\", \"unit_id\"]\n",
    "        ).opts(cmap=\"viridis\", colorbar=True, **cr_opts)\n",
    "    ).relabel(\"Initial Temporal Components\")\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            b.rename(\"b\").compute().astype(np.float32), kdims=[\"width\", \"height\"]\n",
    "        ).opts(**im_opts)\n",
    "    ).relabel(\"Initial Background Sptial\")\n",
    "    + datashade(hv.Curve(f.rename(\"f\").compute(), kdims=[\"frame\"]), min_alpha=200)\n",
    "    .opts(**cr_opts)\n",
    "    .relabel(\"Initial Background Temporal\")\n",
    ").cols(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CNMF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "This section assumes you already have some background knowledge about the CNMF algorithm.\n",
    "Please refer to the [the paper](https://www.sciencedirect.com/science/article/pii/S0896627315010843) for detailed information of the algorithm. \n",
    "\n",
    "As a quick recap, here is the essential idea of CNMF:\n",
    "We believe our movie, `Y`, with dimensions `height`, `width` and `frame`, can be written in (and thus broken down as) the following equation:\n",
    "\n",
    "$$\\mathbf{Y} = \\mathbf{A}^T \\mathbf{C} + \\mathbf{b}^T \\mathbf{f} + \\epsilon$$\n",
    "\n",
    "where:\n",
    "\n",
    "* `A` is the spatial footprint of each unit, with dimension `height`, `width` and `unit_id`.\n",
    "* `C` is the temporal activities of each unit, with dimension `unit_id` and `frame`.\n",
    "* `b` and `f` are the spatial footprint and temporal activities of some background, respectively.\n",
    "* $\\epsilon$ is the noise.\n",
    "\n",
    "Note that strictly speaking, matrix multiplication is usually only defined for two dimensional matrices, but our `A` here has three dimensions, so in fact we are taking the [tensor product](https://en.wikipedia.org/wiki/Tensor_product) of `A` and `C`, reducing the dimension `unit_id`.\n",
    "This might seem to complicate things (compared to just treating `height` and `width` as one flattened `spatial` dimension), but it ends up making some sense.\n",
    "When you take a dot product of any two \"matrices\" on a certain **dimension**, all that is happening is a **product** followed by a **sum** -- you take the product for all pairs of numbers with the the same indexes from the two \"matrices\", and then you take the sum of all those products along the dimension.\n",
    "Thus when we take the tensor product of `A` and `C`, we are actually multiplying all those numbers in dimension `height`, `width` and `frame`, matched by `unit_id`, and then take the sum.\n",
    "Conceptually, for each unit, we are weighting the spatial footprint (`height` and `width`) by the fluorecense of that unit on given `frame`, which is the **product**, and then we are collapsing all units together, which is the **sum**.\n",
    "With that, the equation above is trying to say that our movie is made up of a weighted sum of the spatial footprint and temporal activities of all units, plus some background and noise.\n",
    "\n",
    "Now, there is another rule about `C` that separates it from background and noise:\n",
    "Each \"row\" of `C`, which is the temporal trace for each unit, should be described as an [autoregressive process](https://en.wikipedia.org/wiki/Autoregressive_model) (AR process), with a parameter `p` defining the **order** of the AR process:\n",
    "\n",
    "$$ c(t) = \\sum_{i=0}^{p}\\gamma_i c(t-i) + s(t) + \\epsilon$$\n",
    "\n",
    "where:\n",
    "\n",
    "* $c(t)$ is the calcium concentration at time (`frame`) $t$.\n",
    "* $s(t)$ is spike/firing rate at time $t$.\n",
    "* $\\epsilon$ is noise.\n",
    "\n",
    "Basically, this equation is trying to say that at any given time $t$, the calcium concentration at that moment $c(t)$ depends on the spike at that moment $s(t)$, as well as its own history up to `p` time-steps back $c(t-i)$, scaled by some parameters $\\gamma_i$s, plus some noise $\\epsilon$.\n",
    "Another intuition of this equation comes from looking at different `p`s: when `p=0`, the calcium concentration is an exact copy of the spiking activities, which is probably not true.\n",
    "When `p=1`, the calcium concentration has an instant rise in response to a spike followed by an exponential decay.\n",
    "When `p=2`, calcium concentration has some rise time following a spike and an exponential decay.\n",
    "\n",
    "With all this in mind, CNMF tries to find the spatial matrix (`A`) and temporal activity (`C`) (along with `b` and `f`) that best describe `Y`.\n",
    "There are a few more important practical concerns: Firstly we cannot solve this problem in one shot -- we need to iteratively and separately update `A` and `C` to approach the true solution.\n",
    "Often enough,  two iterations  after the initialization seem to give good enough results, but you can always add more iterations.\n",
    "Secondly, by intuition you may define \"best describe `Y`\" as the results that minimize the noise/error $\\epsilon$.\n",
    "However we have to control for the [sparsity](https://en.wikipedia.org/wiki/Sparse_matrix) of our model as well, since we do not want every little random pixel that happens to correlate with a cell to be counted as part of the spatial footprint of the cell (non-sparse `A`), nor do we want a tiny spike at every frame trying to explain every noisy peak we observe (non-sparse `C`).\n",
    "Thus, the balance between fidelity (minimizing error) and sparsity (minimizing non-zero entries) is an important concern for both the spatial and temporal update."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## estimate spatial noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Prior to performing first spatial update, we need to get a sense of how much noise is expected, which we will then feed into CNMF.\n",
    "We compute a Fast Fourier transform (fft) for every pixel independently, and estimate noise from its power spectral density [power spectral density](https://en.wikipedia.org/wiki/Spectral_density).\n",
    "Recall the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_get_noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The `noise_range` is specified as a fraction sampling frequency.\n",
    "`sn_spatial` is the resulting estimation of noise power for each pixel.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "It is recommended to use the cut-off frequency you find during [peak-noise-ratio refine](#peak-noise-ratio-refine) as the lower bound of `noise_range` to be consistent.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "sn_spatial = get_noise_fft(Y_hw_chk, **param_get_noise)\n",
    "sn_spatial = save_minian(sn_spatial.rename(\"sn_spatial\"), intpath, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## first spatial update"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### parameter exploration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "We will now do some parameter exploration before actually performing the first spatial update.\n",
    "For this purpose, we randomly select 10 units to pass to `update_spatial` with different parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if interactive:\n",
    "    units = np.random.choice(A_init.coords[\"unit_id\"], 10, replace=False)\n",
    "    units.sort()\n",
    "    A_sub = A_init.sel(unit_id=units).persist()\n",
    "    C_sub = C_init.sel(unit_id=units).persist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "We again use a `for` loop and collect results from different runs into `A_dict` and `C_dict`, whose keys are different parameters and values are results corresponding to running `update_spatial` with those parameters.\n",
    "The function `visualize_spatial_update` can directly take dictionaries like these and handle all the visualization details.\n",
    "As a demonstration, the `sparse_penal` parameter is the only one parameter in `update_spatial` that we are interested in playing with.\n",
    "But you can add more to suit your need.\n",
    "See [the next section](#spatial-update) for discussion of parameters for `update_spatial`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "if interactive:\n",
    "    sprs_ls = [0.005, 0.01, 0.05]\n",
    "    A_dict = dict()\n",
    "    C_dict = dict()\n",
    "    for cur_sprs in sprs_ls:\n",
    "        cur_A, cur_b, cur_f, cur_mask = update_spatial(\n",
    "            Y_hw_chk,\n",
    "            A_sub,\n",
    "            b,\n",
    "            C_sub,\n",
    "            f,\n",
    "            sn_spatial,\n",
    "            in_memory=False,\n",
    "            update_background=False,\n",
    "            dl_wnd=param_first_spatial[\"dl_wnd\"],\n",
    "            sparse_penal=cur_sprs,\n",
    "        )\n",
    "        if cur_A.sizes[\"unit_id\"]:\n",
    "            A_dict[cur_sprs] = cur_A.compute()\n",
    "            C_dict[cur_sprs] = C_sub.sel(unit_id=cur_mask).compute()\n",
    "    hv_res = visualize_spatial_update(A_dict, C_dict, kdims=[\"sparse penalty\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The visualization for spatial update composed of the resulting spatial footprints and binarized spatil footprints on the left, as well as temporal components of units involved in the update.\n",
    "Note that if units are dropped during the update their temporal components won't be shown.\n",
    "Here we display the visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "if interactive:\n",
    "    display(hv_res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<strong>choosing sparse penalty for spatial update</strong>\n",
    "\n",
    "Here is an example of what you might see during the parameter exploration:\n",
    "    \n",
    "<div stype=\"clear:both\"><img src=\"img/param_spatial_update.png\" style=\"width: 70%\"/></div>\n",
    "\n",
    "As you can see the `sparse_penal` parameter directly controls the overal sparsity of the resulting spatial footprints.\n",
    "When `sparse_penal=0.01`, the spatial footprints extend far away from the centorids of cells, resulting in high overlap between cells and unnatrual shapes.\n",
    "This is more evident from the binary spatial footprints.\n",
    "At the same time, when `sparse_penal=1`, the algorithm become too strict, and only one cell is left as acceptable.\n",
    "This is also not desirable and we usually want to avoid dropping cells during spatial update.\n",
    "Hence in this example, `sparse_penal=0.3` is considered a good choice among the three cases.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### spatial update"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Given the data `Y` and our units' activity `C` from previous the update, we want to find the spatial footprints `A` such that:\n",
    "\n",
    "1. the error `Y - A.dot(C, 'unit_id')` is as small as possible.\n",
    "1. the [l1-norm](http://mathworld.wolfram.com/L1-Norm.html) of `A` is as small as possible.\n",
    "    Here the l1-norm is a proxy to control for the sparsity of `A`.\n",
    "    \n",
    "Now, updating the `A` matrix altogether is very computationally demanding, and it is much better if we can breakdown our big problem into smaller chunks that can be parallelized.\n",
    "One of the benefits of CNMF is its ability to demix cells, so it is best to keep updating all the cells together.\n",
    "However, it is fine to treat each pixel as independent and update different pixels independently.\n",
    "Please refer to the [API reference](https://minian.readthedocs.io/page/api/minian.cnmf.html#minian-cnmf-update_spatial) of `update_spatial` for details on all the parameters and the mathematical formulation of per-pixel optimization problem.\n",
    "Here we focus on a few parameters that requires more attention.\n",
    "Recall the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_first_spatial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "When carrying out spatial update, it is very inefficient to estimate a weight for all the cells for each pixel.\n",
    "Instead, when updating for each pixel, we only want to consider cells that are close by and ignore cells that are very far from the pixel being updated.\n",
    "For this purpose we carry out a [morphological dilation](https://homepages.inf.ed.ac.uk/rbf/HIPR2/dilate.htm) on the spatial footprints of each cell using the previous estimation of `A`.\n",
    "We then binarize this dilated spatial footprints matrix and use it as a mask.\n",
    "Then when updating for each pixel, only cells that have non-zero values in the mask on this pixel will be considered for update.\n",
    "The parameter `dl_wnd` controls the window size of the morphological dilation operation.\n",
    "\n",
    "The scalar `sparse_penal` controls the balance between error objective and the l1-norm term.\n",
    "The higher the `sparse_penal`, the sparser the result will become.\n",
    "It is hard to estimate theoretically, and the best way to set this is through parameter exploration.\n",
    "\n",
    "The boolean parameter `update_background` controls whether we want to update the background in this step.\n",
    "This is the only place in the pipeline that the background will be updated.\n",
    "The way it is updated is by essentially treating `b` as another `unit` and updating it according to the temporal activity `f`.\n",
    "Since the morphology-based [background removal](#backgroun-removal) works well at cleaning the backgrounds, this update usually result in an empty estimation of the background.\n",
    "\n",
    "Lastly, it is often convenient to filter out cells that has either too large or too small spatial footprints at this step.\n",
    "The `size_thres` controls the range of area (number of non-zero pixels) of the spatial footprints that will be accepted."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell would carry out spatial update.\n",
    "Be sure to [change the parameters](https://minian.readthedocs.io/page/start_guide/faq.html#i-don-t-know-python-can-i-still-use-the-pipeline) based on visualization results before running the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "A_new, b_new, f_new, mask = update_spatial(\n",
    "    Y_hw_chk, A, b, C, f, sn_spatial, **param_first_spatial\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### visualization of spatial footprints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "We visualize the result of spatial update by plotting both the spatial footprints and the binarized spatial footprints before and after the update."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "opts = dict(\n",
    "    plot=dict(height=A.sizes[\"height\"], width=A.sizes[\"width\"], colorbar=True),\n",
    "    style=dict(cmap=\"Viridis\"),\n",
    ")\n",
    "(\n",
    "    regrid(\n",
    "        hv.Image(\n",
    "            A.sum(\"unit_id\").compute().astype(np.float32).rename(\"A\"),\n",
    "            kdims=[\"width\", \"height\"],\n",
    "        ).opts(**opts)\n",
    "    ).relabel(\"Spatial Footprints Initial\")\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            (A.fillna(0) > 0).sum(\"unit_id\").compute().astype(np.uint8).rename(\"A\"),\n",
    "            kdims=[\"width\", \"height\"],\n",
    "        ).opts(**opts)\n",
    "    ).relabel(\"Binary Spatial Footprints Initial\")\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            A_new.sum(\"unit_id\").compute().astype(np.float32).rename(\"A\"),\n",
    "            kdims=[\"width\", \"height\"],\n",
    "        ).opts(**opts)\n",
    "    ).relabel(\"Spatial Footprints First Update\")\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            (A_new > 0).sum(\"unit_id\").compute().astype(np.uint8).rename(\"A\"),\n",
    "            kdims=[\"width\", \"height\"],\n",
    "        ).opts(**opts)\n",
    "    ).relabel(\"Binary Spatial Footprints First Update\")\n",
    ").cols(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### visualization of background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "We also visualize the spatial and temporal component of the background before and after the update."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.55))\n",
    "opts_im = dict(\n",
    "    plot=dict(height=b.sizes[\"height\"], width=b.sizes[\"width\"], colorbar=True),\n",
    "    style=dict(cmap=\"Viridis\"),\n",
    ")\n",
    "opts_cr = dict(plot=dict(height=b.sizes[\"height\"], width=b.sizes[\"height\"] * 2))\n",
    "(\n",
    "    regrid(\n",
    "        hv.Image(b.compute().astype(np.float32), kdims=[\"width\", \"height\"]).opts(\n",
    "            **opts_im\n",
    "        )\n",
    "    ).relabel(\"Background Spatial Initial\")\n",
    "    + hv.Curve(f.compute().rename(\"f\").astype(np.float16), kdims=[\"frame\"])\n",
    "    .opts(**opts_cr)\n",
    "    .relabel(\"Background Temporal Initial\")\n",
    "    + regrid(\n",
    "        hv.Image(b_new.compute().astype(np.float32), kdims=[\"width\", \"height\"]).opts(\n",
    "            **opts_im\n",
    "        )\n",
    "    ).relabel(\"Background Spatial First Update\")\n",
    "    + hv.Curve(f_new.compute().rename(\"f\").astype(np.float16), kdims=[\"frame\"])\n",
    "    .opts(**opts_cr)\n",
    "    .relabel(\"Background Temporal First Update\")\n",
    ").cols(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### save results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Once we are satisfied with the result, we can commit by saving them to the intermediate folder.\n",
    "Note that this overwrites the result of previous updates.\n",
    "If you want to keep previous results around, consider `rename` the variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "A = save_minian(\n",
    "    A_new.rename(\"A\"),\n",
    "    intpath,\n",
    "    overwrite=True,\n",
    "    chunks={\"unit_id\": 1, \"height\": -1, \"width\": -1},\n",
    ")\n",
    "b = save_minian(b_new.rename(\"b\"), intpath, overwrite=True)\n",
    "f = save_minian(\n",
    "    f_new.chunk({\"frame\": chk[\"frame\"]}).rename(\"f\"), intpath, overwrite=True\n",
    ")\n",
    "C = C.sel(unit_id=A.coords[\"unit_id\"].values)\n",
    "C_chk = C_chk.sel(unit_id=A.coords[\"unit_id\"].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## first temporal update"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### parameter exploration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "First off we randomly select 10 cells to do parameter exploring."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if interactive:\n",
    "    units = np.random.choice(A.coords[\"unit_id\"], 10, replace=False)\n",
    "    units.sort()\n",
    "    A_sub = A.sel(unit_id=units).persist()\n",
    "    C_sub = C_chk.sel(unit_id=units).persist()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Now we move on to the parameter exploring of temporal update.\n",
    "This time there is more than one parameter to play with for temporal update, and we now have four `list`s of parameters: `p_ls`, `sprs_ls`, `add_ls`, and `noise_ls`.\n",
    "We use [itertools.product](https://docs.python.org/3.7/library/itertools.html#itertools.product) to iterate through all possible combinations of parameter values, which save us from nested `for` loops.\n",
    "Please see [next section](#temporal-update) for details of the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "if interactive:\n",
    "    p_ls = [1]\n",
    "    sprs_ls = [0.1, 0.5, 1, 2]\n",
    "    add_ls = [20]\n",
    "    noise_ls = [0.06]\n",
    "    YA_dict, C_dict, S_dict, g_dict, sig_dict, A_dict = [dict() for _ in range(6)]\n",
    "    YrA = (\n",
    "        compute_trace(Y_fm_chk, A_sub, b, C_sub, f)\n",
    "        .persist()\n",
    "        .chunk({\"unit_id\": 1, \"frame\": -1})\n",
    "    )\n",
    "    for cur_p, cur_sprs, cur_add, cur_noise in itt.product(\n",
    "        p_ls, sprs_ls, add_ls, noise_ls\n",
    "    ):\n",
    "        ks = (cur_p, cur_sprs, cur_add, cur_noise)\n",
    "        print(\n",
    "            \"p:{}, sparse penalty:{}, additional lag:{}, noise frequency:{}\".format(\n",
    "                cur_p, cur_sprs, cur_add, cur_noise\n",
    "            )\n",
    "        )\n",
    "        cur_C, cur_S, cur_b0, cur_c0, cur_g, cur_mask = update_temporal(\n",
    "            A_sub,\n",
    "            C_sub,\n",
    "            YrA=YrA,\n",
    "            sparse_penal=cur_sprs,\n",
    "            p=cur_p,\n",
    "            use_smooth=True,\n",
    "            add_lag=cur_add,\n",
    "            noise_freq=cur_noise,\n",
    "        )\n",
    "        YA_dict[ks], C_dict[ks], S_dict[ks], g_dict[ks], sig_dict[ks], A_dict[ks] = (\n",
    "            YrA.compute(),\n",
    "            cur_C.compute(),\n",
    "            cur_S.compute(),\n",
    "            cur_g.compute(),\n",
    "            (cur_C + cur_b0 + cur_c0).compute(),\n",
    "            A_sub.compute(),\n",
    "        )\n",
    "    hv_res = visualize_temporal_update(\n",
    "        YA_dict,\n",
    "        C_dict,\n",
    "        S_dict,\n",
    "        g_dict,\n",
    "        sig_dict,\n",
    "        A_dict,\n",
    "        kdims=[\"p\", \"sparse penalty\", \"additional lag\", \"noise frequency\"],\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The visualization for temporal update composed of a main plot showing an overlay of multiple temporal signals.\n",
    "The \"Raw Signal\" represent the temporal fluctuation of the unit in the raw movie `YrA`, which is the target of the fitting process during temporal update.\n",
    "The \"Fitted Calcium Calcium Trace\" and \"Fitted Spikes\" represent the model-estimated calciums dynamic `C` and model-estimated deconvolved spikes `S` respectively.\n",
    "The \"Fitted Signal\" is mostly identical to `C` except for some baseline and initial calcium decay.\n",
    "Additionally, a modeled calcium response to a single spike is shown at bottom left, helping you visualize the modeled AR process.\n",
    "The spatial footprint of current unit is shown at bottom right."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "if interactive:\n",
    "    display(hv_res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "<strong>choosing sparse penalty for temporal update</strong>\n",
    "\n",
    "Here is an example of what you might see during the parameter exploration:\n",
    "    \n",
    "<div stype=\"clear:both\"><img src=\"img/param_temporal_update.png\" style=\"width: 50%\"/></div>\n",
    "\n",
    "As you can see the `sparse_penal` parameter directly controls the overal sparsity of the resulting calcium dynamic and deconvolved spikes.\n",
    "When `sparse_penal=1`, the fitted spikes contains lots of small valued fluctuations, that mostly correspond to high-frequency noise instead of real calcium dynamic.\n",
    "At the same time, when `sparse_penal=10`, the sparse penalty is too large, and a lot of real calcium dynamics, as evident in the raw signal, are left out in the fitted calcium and spike traces.\n",
    "Hence in this example, `sparse_penal=3` is considered a good choice among the three cases.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### temporal update"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Given the spatial footprints of units `A`, our goal is to find the calcium dynamic `C` as well as the deconvolved spike signal `S` for each cell such that:\n",
    "\n",
    "1. the error `Y - A.dot(C, 'unit_id')` is as small as possible.\n",
    "1. the relationship between `C` and `S` can be described by an AR process.\n",
    "1. the l1-norm of `C` is as small as possible, which serve as a proxy to control the sparsity of both `C` and `S`.\n",
    "\n",
    "Now, doing the matrix multiplication between `A` and `C` at each iteration is very expensive.\n",
    "Since presumably we already have a good estimation of the spatial footprints of units, it would be desirable if we can take advantage of that, and project our movie data `Y` using `A` to get a one dimensional activities across time for each unit.\n",
    "The following cell carry out that computation and return a variable `YrA` which only have dimensions `frame` and `unit_id`.\n",
    "The `YrA` variable can be interpreted as the demixed temporal activity for each unit.\n",
    "Then during temporal update, we can use `YrA` as the target of the fitting process (instead of `Y`) and we eliminate the need to carry out demanding multiplications with `A`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "YrA = save_minian(\n",
    "    compute_trace(Y_fm_chk, A, b, C_chk, f).rename(\"YrA\"),\n",
    "    intpath,\n",
    "    overwrite=True,\n",
    "    chunks={\"unit_id\": 1, \"frame\": -1},\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Once we have `YrA` we can proceed to temporal update for each unit.\n",
    "Please refer to the [API reference](https://minian.readthedocs.io/page/api/minian.cnmf.html#minian-cnmf-update_temporal) of `update_temporal` for details on all parameters and mathematical formulatioin of the per-unit optimization problem.\n",
    "Here we focus on a few parameters that requires more attention:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_first_temporal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Although we obtained `YrA` for each unit, it is still beneficial to group together units that heavily overlap, and carry out the update process group-wise instead of independently for each unit.\n",
    "This way the relative numerical relationship between heavily-overlapping units are better preserved.\n",
    "For this purpose, we compute pairwise [Jaccard index](https://en.wikipedia.org/wiki/Jaccard_index) for all units, which serve as a metric of overlap between spatial footprints of units.\n",
    "The parameter `jac_thres` is the threshold of Jaccard index above which units will be grouped together transitively.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "Since the grouping of units is transitive and uses binarized spatial footprints, in practive you may find large amount of units being grouped together if the spatial footprints is not sparse enough.\n",
    "This lead to significantly increased memory demand, and potential \"masking\" effect between units (activities of highly active units suppress the activities of less active units when grouped together).\n",
    "If this happens, consider either refine the spatial footprints or increase the `jac_thres`.\n",
    "</div>\n",
    "\n",
    "Regarding the actual optimization process, the first thing we want to determine is order of the AR process `p`.\n",
    "Usually, `p=1` is good enough and tend to result in multiple spikes in the deconvolved signal `S` accounting for a single rise of calcium concentration in `C`.\n",
    "However, if you believe the rise time of your calcium signal is not neglectable, and a single rise of calcium concentration should be modeled as a single calcium event in the deconvolved signal `S`, then `p=2` is a better choice since it allows for modeling of non-zero rise time.\n",
    "\n",
    "Next, we estimate the AR coefficients from the auto-covariance of the `YrA` signal for each cell.\n",
    "Two additional steps can improve the reliability of this estimation.\n",
    "First, the `YrA` can be smoothed for estimation of the AR coefficients.\n",
    "This can help prevent the high frequency noise in `YrA` biasing the estimation of AR proces in to very fast dynamics.\n",
    "The `noise_freq` parameter is the cut-off frequency of this low-pass filtering.\n",
    "Secondly, although in theory only `p` auto-covariances are needed to solve `p` number of AR coefficients, we can use auto-covariance more than `p` time lags and solve for AR coefficients using least square.\n",
    "This help make the estimation numerically more stable.\n",
    "The `add_lag` parameter is the number of additional auto-covariance to use for estimating AR coefficients.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "The `noise_freq` should be set to those determined in [peak-noise-ratio refine](#peak-noise-ratio-refine) to be consistent.\n",
    "The `add_lag` usually have limited impact on the result as long as it is large enough.\n",
    "Pragmatically we have found `add_lag=20` works for most cases.\n",
    "However calcium dynamic are very slow in your data, you might have to increase `add_lag` so that auto-covariance with further time lag may contribute to the estimation of AR coefficients.\n",
    "</div>\n",
    "\n",
    "Finally, the scalar `sparse_penal` controls the balance between the error and l1-norm of `C` in the optimization objective.\n",
    "The higher the value, the sparser both `C` and `S` will become.\n",
    "It is hard to estimate theoretically, and the best way to set this is through parameter exploration.\n",
    "Note that despite the name, this is a completely different parameter than the one in spatial updates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell would carry out temporal update.\n",
    "Be sure to [change the parameters](https://minian.readthedocs.io/page/start_guide/faq.html#i-don-t-know-python-can-i-still-use-the-pipeline) based on visualization results before running the following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "C_new, S_new, b0_new, c0_new, g, mask = update_temporal(\n",
    "    A, C, YrA=YrA, **param_first_temporal\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### visualization of temporal components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The cell below plot both the calcium dynamic and the deconvolved spikes before and after the temporal update."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "opts_im = dict(frame_width=500, aspect=2, colorbar=True, cmap=\"Viridis\")\n",
    "(\n",
    "    regrid(\n",
    "        hv.Image(\n",
    "            C.compute().astype(np.float32).rename(\"ci\"), kdims=[\"frame\", \"unit_id\"]\n",
    "        ).opts(**opts_im)\n",
    "    ).relabel(\"Temporal Trace Initial\")\n",
    "    + hv.Div(\"\")\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            C_new.compute().astype(np.float32).rename(\"c1\"), kdims=[\"frame\", \"unit_id\"]\n",
    "        ).opts(**opts_im)\n",
    "    ).relabel(\"Temporal Trace First Update\")\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            S_new.compute().astype(np.float32).rename(\"s1\"), kdims=[\"frame\", \"unit_id\"]\n",
    "        ).opts(**opts_im)\n",
    "    ).relabel(\"Spikes First Update\")\n",
    ").cols(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### visualization of dropped units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Since it is common to drop units during temporal update, it is helpful to check whether the dropped units contain any real temporal activities.\n",
    "The cell below plot both the `YrA` signal and spatial footprints for the dropped units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "if interactive:\n",
    "    h, w = A.sizes[\"height\"], A.sizes[\"width\"]\n",
    "    im_opts = dict(aspect=w / h, frame_width=500, cmap=\"Viridis\")\n",
    "    cr_opts = dict(aspect=3, frame_width=1000)\n",
    "    bad_units = mask.where(mask == False, drop=True).coords[\"unit_id\"].values\n",
    "    if len(bad_units) > 0:\n",
    "        hv_res = (\n",
    "            hv.NdLayout(\n",
    "                {\n",
    "                    \"Spatial Footprint\": Dynamic(\n",
    "                        hv.Dataset(A.sel(unit_id=bad_units).compute().rename(\"A\"))\n",
    "                        .to(hv.Image, kdims=[\"width\", \"height\"])\n",
    "                        .opts(**im_opts)\n",
    "                    ),\n",
    "                    \"Spatial Footprints of Accepted Units\": Dynamic(\n",
    "                        hv.Image(\n",
    "                            A.sel(unit_id=mask).sum(\"unit_id\").compute().rename(\"A\"),\n",
    "                            kdims=[\"width\", \"height\"],\n",
    "                        ).opts(**im_opts)\n",
    "                    ),\n",
    "                }\n",
    "            )\n",
    "            + datashade(\n",
    "                hv.Dataset(YrA.sel(unit_id=bad_units).rename(\"raw\")).to(\n",
    "                    hv.Curve, kdims=[\"frame\"]\n",
    "                )\n",
    "            )\n",
    "            .opts(**cr_opts)\n",
    "            .relabel(\"Temporal Trace\")\n",
    "        ).cols(1)\n",
    "        display(hv_res)\n",
    "    else:\n",
    "        print(\"No rejected units to display\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### visualization of accepted units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Lastly, we can go through all the accepted units and visualize the result of temporal update.\n",
    "The following cell visualize all the units in the same way as during parameter exploration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "if interactive:\n",
    "    sig = C_new + b0_new + c0_new\n",
    "    display(\n",
    "        visualize_temporal_update(\n",
    "            YrA.sel(unit_id=mask),\n",
    "            C_new,\n",
    "            S_new,\n",
    "            g,\n",
    "            sig,\n",
    "            A.sel(unit_id=mask),\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### save results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "C = save_minian(\n",
    "    C_new.rename(\"C\").chunk({\"unit_id\": 1, \"frame\": -1}), intpath, overwrite=True\n",
    ")\n",
    "C_chk = save_minian(\n",
    "    C.rename(\"C_chk\"),\n",
    "    intpath,\n",
    "    overwrite=True,\n",
    "    chunks={\"unit_id\": -1, \"frame\": chk[\"frame\"]},\n",
    ")\n",
    "S = save_minian(\n",
    "    S_new.rename(\"S\").chunk({\"unit_id\": 1, \"frame\": -1}), intpath, overwrite=True\n",
    ")\n",
    "b0 = save_minian(\n",
    "    b0_new.rename(\"b0\").chunk({\"unit_id\": 1, \"frame\": -1}), intpath, overwrite=True\n",
    ")\n",
    "c0 = save_minian(\n",
    "    c0_new.rename(\"c0\").chunk({\"unit_id\": 1, \"frame\": -1}), intpath, overwrite=True\n",
    ")\n",
    "A = A.sel(unit_id=C.coords[\"unit_id\"].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## merge units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "One thing CNMF cannot do is merge together units that belong to the same cell, and we have to do it manually between different iteration of updates.\n",
    "Recall the parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "outputs": [],
   "source": [
    "param_first_merge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The idea is straight-forward and based purely on correlation of temporal activities.\n",
    "Any units whose spatial footprints share at least one pixel are considered potential targets for merging, and any of these units that have a correlation of temporal activities higher than `thres_corr` will be merged. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "A_mrg, C_mrg, [sig_mrg] = unit_merge(A, C, [C + b0 + c0], **param_first_merge)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Now you can visualize the results of unit merging.\n",
    "The left panel shows the original temporal signal, while the right panel shows the temporal signal after merging."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "opts_im = dict(frame_width=500, aspect=2, colorbar=True, cmap=\"Viridis\")\n",
    "(\n",
    "    regrid(\n",
    "        hv.Image(\n",
    "            C.compute().astype(np.float32).rename(\"c1\"), kdims=[\"frame\", \"unit_id\"]\n",
    "        )\n",
    "        .relabel(\"Temporal Signals Before Merge\")\n",
    "        .opts(**opts_im)\n",
    "    )\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            C_mrg.compute().astype(np.float32).rename(\"c2\"), kdims=[\"frame\", \"unit_id\"]\n",
    "        )\n",
    "        .relabel(\"Temporal Signals After Merge\")\n",
    "        .opts(**opts_im)\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Once we are satisfied with the result of merging we can commit to saving them to intermediate folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "A = save_minian(A_mrg.rename(\"A_mrg\"), intpath, overwrite=True)\n",
    "C = save_minian(C_mrg.rename(\"C_mrg\"), intpath, overwrite=True)\n",
    "C_chk = save_minian(\n",
    "    C.rename(\"C_mrg_chk\"),\n",
    "    intpath,\n",
    "    overwrite=True,\n",
    "    chunks={\"unit_id\": -1, \"frame\": chk[\"frame\"]},\n",
    ")\n",
    "sig = save_minian(sig_mrg.rename(\"sig_mrg\"), intpath, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## second spatial update"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "This section analogous to the [first time](#first-spatial-update) we so spatial update except for changes in variable names."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### parameter exploration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if interactive:\n",
    "    units = np.random.choice(A.coords[\"unit_id\"], 10, replace=False)\n",
    "    units.sort()\n",
    "    A_sub = A.sel(unit_id=units).persist()\n",
    "    C_sub = sig.sel(unit_id=units).persist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "if interactive:\n",
    "    sprs_ls = [5e-3, 1e-2, 5e-2]\n",
    "    A_dict = dict()\n",
    "    C_dict = dict()\n",
    "    for cur_sprs in sprs_ls:\n",
    "        cur_A, cur_b, cur_f, cur_mask = update_spatial(\n",
    "            Y_hw_chk,\n",
    "            A_sub,\n",
    "            b,\n",
    "            C_sub,\n",
    "            f,\n",
    "            sn_spatial,\n",
    "            in_memory=False,\n",
    "            update_background=False,\n",
    "            dl_wnd=param_second_spatial[\"dl_wnd\"],\n",
    "            sparse_penal=cur_sprs,\n",
    "        )\n",
    "        if cur_A.sizes[\"unit_id\"]:\n",
    "            A_dict[cur_sprs] = cur_A.compute()\n",
    "            C_dict[cur_sprs] = C_sub.sel(unit_id=cur_mask).compute()\n",
    "    hv_res = visualize_spatial_update(A_dict, C_dict, kdims=[\"sparse penalty\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "if interactive:\n",
    "    display(hv_res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### spatial update"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "A_new, b_new, f_new, mask = update_spatial(\n",
    "    Y_hw_chk, A, b, sig, f, sn_spatial, **param_second_spatial\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### visualization of spatial footprints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "opts = dict(\n",
    "    plot=dict(height=A.sizes[\"height\"], width=A.sizes[\"width\"], colorbar=True),\n",
    "    style=dict(cmap=\"Viridis\"),\n",
    ")\n",
    "(\n",
    "    regrid(\n",
    "        hv.Image(\n",
    "            A.sum(\"unit_id\").compute().astype(np.float32).rename(\"A\"),\n",
    "            kdims=[\"width\", \"height\"],\n",
    "        ).opts(**opts)\n",
    "    ).relabel(\"Spatial Footprints Last\")\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            (A.fillna(0) > 0).sum(\"unit_id\").compute().astype(np.uint8).rename(\"A\"),\n",
    "            kdims=[\"width\", \"height\"],\n",
    "        ).opts(**opts)\n",
    "    ).relabel(\"Binary Spatial Footprints Last\")\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            A_new.sum(\"unit_id\").compute().astype(np.float32).rename(\"A\"),\n",
    "            kdims=[\"width\", \"height\"],\n",
    "        ).opts(**opts)\n",
    "    ).relabel(\"Spatial Footprints New\")\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            (A_new > 0).sum(\"unit_id\").compute().astype(np.uint8).rename(\"A\"),\n",
    "            kdims=[\"width\", \"height\"],\n",
    "        ).opts(**opts)\n",
    "    ).relabel(\"Binary Spatial Footprints New\")\n",
    ").cols(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### visualization of background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.55))\n",
    "opts_im = dict(\n",
    "    plot=dict(height=b.sizes[\"height\"], width=b.sizes[\"width\"], colorbar=True),\n",
    "    style=dict(cmap=\"Viridis\"),\n",
    ")\n",
    "opts_cr = dict(plot=dict(height=b.sizes[\"height\"], width=b.sizes[\"height\"] * 2))\n",
    "(\n",
    "    regrid(\n",
    "        hv.Image(b.compute().astype(np.float32), kdims=[\"width\", \"height\"]).opts(\n",
    "            **opts_im\n",
    "        )\n",
    "    ).relabel(\"Background Spatial Last\")\n",
    "    + hv.Curve(f.compute().rename(\"f\").astype(np.float16), kdims=[\"frame\"])\n",
    "    .opts(**opts_cr)\n",
    "    .relabel(\"Background Temporal Last\")\n",
    "    + regrid(\n",
    "        hv.Image(b_new.compute().astype(np.float32), kdims=[\"width\", \"height\"]).opts(\n",
    "            **opts_im\n",
    "        )\n",
    "    ).relabel(\"Background Spatial New\")\n",
    "    + hv.Curve(f_new.compute().rename(\"f\").astype(np.float16), kdims=[\"frame\"])\n",
    "    .opts(**opts_cr)\n",
    "    .relabel(\"Background Temporal New\")\n",
    ").cols(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### save results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "A = save_minian(\n",
    "    A_new.rename(\"A\"),\n",
    "    intpath,\n",
    "    overwrite=True,\n",
    "    chunks={\"unit_id\": 1, \"height\": -1, \"width\": -1},\n",
    ")\n",
    "b = save_minian(b_new.rename(\"b\"), intpath, overwrite=True)\n",
    "f = save_minian(\n",
    "    f_new.chunk({\"frame\": chk[\"frame\"]}).rename(\"f\"), intpath, overwrite=True\n",
    ")\n",
    "C = C.sel(unit_id=A.coords[\"unit_id\"].values)\n",
    "C_chk = C_chk.sel(unit_id=A.coords[\"unit_id\"].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## second temporal update"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### parameter exploration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "This section is almost analogous to the [first time](#first-temporal-update) except for changes in variable names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if interactive:\n",
    "    units = np.random.choice(A.coords[\"unit_id\"], 10, replace=False)\n",
    "    units.sort()\n",
    "    A_sub = A.sel(unit_id=units).persist()\n",
    "    C_sub = C_chk.sel(unit_id=units).persist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "if interactive:\n",
    "    p_ls = [1]\n",
    "    sprs_ls = [0.1, 0.5, 1, 2]\n",
    "    add_ls = [20]\n",
    "    noise_ls = [0.06]\n",
    "    YA_dict, C_dict, S_dict, g_dict, sig_dict, A_dict = [dict() for _ in range(6)]\n",
    "    YrA = (\n",
    "        compute_trace(Y_fm_chk, A_sub, b, C_sub, f)\n",
    "        .persist()\n",
    "        .chunk({\"unit_id\": 1, \"frame\": -1})\n",
    "    )\n",
    "    for cur_p, cur_sprs, cur_add, cur_noise in itt.product(\n",
    "        p_ls, sprs_ls, add_ls, noise_ls\n",
    "    ):\n",
    "        ks = (cur_p, cur_sprs, cur_add, cur_noise)\n",
    "        print(\n",
    "            \"p:{}, sparse penalty:{}, additional lag:{}, noise frequency:{}\".format(\n",
    "                cur_p, cur_sprs, cur_add, cur_noise\n",
    "            )\n",
    "        )\n",
    "        cur_C, cur_S, cur_b0, cur_c0, cur_g, cur_mask = update_temporal(\n",
    "            A_sub,\n",
    "            C_sub,\n",
    "            YrA=YrA,\n",
    "            sparse_penal=cur_sprs,\n",
    "            p=cur_p,\n",
    "            use_smooth=True,\n",
    "            add_lag=cur_add,\n",
    "            noise_freq=cur_noise,\n",
    "        )\n",
    "        YA_dict[ks], C_dict[ks], S_dict[ks], g_dict[ks], sig_dict[ks], A_dict[ks] = (\n",
    "            YrA.compute(),\n",
    "            cur_C.compute(),\n",
    "            cur_S.compute(),\n",
    "            cur_g.compute(),\n",
    "            (cur_C + cur_b0 + cur_c0).compute(),\n",
    "            A_sub.compute(),\n",
    "        )\n",
    "    hv_res = visualize_temporal_update(\n",
    "        YA_dict,\n",
    "        C_dict,\n",
    "        S_dict,\n",
    "        g_dict,\n",
    "        sig_dict,\n",
    "        A_dict,\n",
    "        kdims=[\"p\", \"sparse penalty\", \"additional lag\", \"noise frequency\"],\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "if interactive:\n",
    "    display(hv_res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### temporal update"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "YrA = save_minian(\n",
    "    compute_trace(Y_fm_chk, A, b, C_chk, f).rename(\"YrA\"),\n",
    "    intpath,\n",
    "    overwrite=True,\n",
    "    chunks={\"unit_id\": 1, \"frame\": -1},\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "C_new, S_new, b0_new, c0_new, g, mask = update_temporal(\n",
    "    A, C, YrA=YrA, **param_second_temporal\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### visualization of temporal components"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "\n",
    "You might notice \"Temporal Trace Last\" has different number of units than \"Spikes Last\".\n",
    "This is because `C` was merged in previous steps while `S` was not.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "opts_im = dict(frame_width=500, aspect=2, colorbar=True, cmap=\"Viridis\")\n",
    "(\n",
    "    regrid(\n",
    "        hv.Image(\n",
    "            C.compute().astype(np.float32).rename(\"c1\"), kdims=[\"frame\", \"unit_id\"]\n",
    "        ).opts(**opts_im)\n",
    "    ).relabel(\"Temporal Trace Last\")\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            S.compute().astype(np.float32).rename(\"s1\"), kdims=[\"frame\", \"unit_id\"]\n",
    "        ).opts(**opts_im)\n",
    "    ).relabel(\"Spikes Last\")\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            C_new.compute().astype(np.float32).rename(\"c2\"), kdims=[\"frame\", \"unit_id\"]\n",
    "        ).opts(**opts_im)\n",
    "    ).relabel(\"Temporal Trace New\")\n",
    "    + regrid(\n",
    "        hv.Image(\n",
    "            S_new.compute().astype(np.float32).rename(\"s2\"), kdims=[\"frame\", \"unit_id\"]\n",
    "        ).opts(**opts_im)\n",
    "    ).relabel(\"Spikes New\")\n",
    ").cols(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### visualization of dropped units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "if interactive:\n",
    "    h, w = A.sizes[\"height\"], A.sizes[\"width\"]\n",
    "    im_opts = dict(aspect=w / h, frame_width=500, cmap=\"Viridis\")\n",
    "    cr_opts = dict(aspect=3, frame_width=1000)\n",
    "    bad_units = mask.where(mask == False, drop=True).coords[\"unit_id\"].values\n",
    "    if len(bad_units) > 0:\n",
    "        hv_res = (\n",
    "            hv.NdLayout(\n",
    "                {\n",
    "                    \"Spatial Footprint\": Dynamic(\n",
    "                        hv.Dataset(A.sel(unit_id=bad_units).compute().rename(\"A\"))\n",
    "                        .to(hv.Image, kdims=[\"width\", \"height\"])\n",
    "                        .opts(**im_opts)\n",
    "                    ),\n",
    "                    \"Spatial Footprints of Accepted Units\": Dynamic(\n",
    "                        hv.Image(\n",
    "                            A.sel(unit_id=mask).sum(\"unit_id\").compute().rename(\"A\"),\n",
    "                            kdims=[\"width\", \"height\"],\n",
    "                        ).opts(**im_opts)\n",
    "                    ),\n",
    "                }\n",
    "            )\n",
    "            + datashade(\n",
    "                hv.Dataset(YrA.sel(unit_id=bad_units).rename(\"raw\")).to(\n",
    "                    hv.Curve, kdims=[\"frame\"]\n",
    "                )\n",
    "            )\n",
    "            .opts(**cr_opts)\n",
    "            .relabel(\"Temporal Trace\")\n",
    "        ).cols(1)\n",
    "        display(hv_res)\n",
    "    else:\n",
    "        print(\"No rejected units to display\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### visualization of accepted units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.6))\n",
    "if interactive:\n",
    "    sig = C_new + b0_new + c0_new\n",
    "    display(\n",
    "        visualize_temporal_update(\n",
    "            YrA.sel(unit_id=mask),\n",
    "            C_new,\n",
    "            S_new,\n",
    "            g,\n",
    "            sig,\n",
    "            A.sel(unit_id=mask),\n",
    "        )\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### save results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "C = save_minian(\n",
    "    C_new.rename(\"C\").chunk({\"unit_id\": 1, \"frame\": -1}), intpath, overwrite=True\n",
    ")\n",
    "C_chk = save_minian(\n",
    "    C.rename(\"C_chk\"),\n",
    "    intpath,\n",
    "    overwrite=True,\n",
    "    chunks={\"unit_id\": -1, \"frame\": chk[\"frame\"]},\n",
    ")\n",
    "S = save_minian(\n",
    "    S_new.rename(\"S\").chunk({\"unit_id\": 1, \"frame\": -1}), intpath, overwrite=True\n",
    ")\n",
    "b0 = save_minian(\n",
    "    b0_new.rename(\"b0\").chunk({\"unit_id\": 1, \"frame\": -1}), intpath, overwrite=True\n",
    ")\n",
    "c0 = save_minian(\n",
    "    c0_new.rename(\"c0\").chunk({\"unit_id\": 1, \"frame\": -1}), intpath, overwrite=True\n",
    ")\n",
    "A = A.sel(unit_id=C.coords[\"unit_id\"].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "The following cell calls `generate_videos` to create a video that can help us quickly visualize the overall quality of the whole pipeline.\n",
    "By default, this video is saved under `dpath`.\n",
    "Please see [API reference](https://minian.readthedocs.io/page/api/minian.visualization.html#minian-visualization-generate_videos) of `generate_videos` for the interpretation of resulting video."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "generate_videos(varr.sel(subset), Y_fm_chk, A=A, C=C_chk, vpath=dpath)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Here we use `CNMFViewer` to visualize the final results.\n",
    "Please see [API reference](https://minian.readthedocs.io/page/api/minian.visualization.html#minian-visualization-CNMFViewer) for the function of each panels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "if interactive:\n",
    "    cnmfviewer = CNMFViewer(A=A, C=C, S=S, org=Y_fm_chk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.output(size=int(output_size * 0.35))\n",
    "if interactive:\n",
    "    display(cnmfviewer.show())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## save unit labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "If manual manipulation of `unit_labels` are done during visualization, we should assign them as coordinates to our final results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if interactive:\n",
    "    A = A.assign_coords(unit_labels=(\"unit_id\", cnmfviewer.unit_labels))\n",
    "    C = C.assign_coords(unit_labels=(\"unit_id\", cnmfviewer.unit_labels))\n",
    "    S = S.assign_coords(unit_labels=(\"unit_id\", cnmfviewer.unit_labels))\n",
    "    c0 = c0.assign_coords(unit_labels=(\"unit_id\", cnmfviewer.unit_labels))\n",
    "    b0 = b0.assign_coords(unit_labels=(\"unit_id\", cnmfviewer.unit_labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## save final results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "tutorial"
    ]
   },
   "source": [
    "Finally we commit to saving final results using `param_save_minian`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "A = save_minian(A.rename(\"A\"), **param_save_minian)\n",
    "C = save_minian(C.rename(\"C\"), **param_save_minian)\n",
    "S = save_minian(S.rename(\"S\"), **param_save_minian)\n",
    "c0 = save_minian(c0.rename(\"c0\"), **param_save_minian)\n",
    "b0 = save_minian(b0.rename(\"b0\"), **param_save_minian)\n",
    "b = save_minian(b.rename(\"b\"), **param_save_minian)\n",
    "f = save_minian(f.rename(\"f\"), **param_save_minian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## close cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "client.close()\n",
    "cluster.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  },
  "name": "pipeline.ipynb",
  "toc-autonumbering": false,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": true
 },
 "nbformat": 4,
 "nbformat_minor": 4
}