{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic FISSA usage\n",
    "\n",
    "This notebook contains a step-by-step example of how to use the function-based high-level interface to the [FISSA](https://github.com/rochefort-lab/fissa) toolbox, [fissa.run_fissa](https://fissa.readthedocs.io/en/1.0.0/source/packages/fissa.core.html#fissa.core.run_fissa).\n",
    "\n",
    "For more details about the methodology behind FISSA, please see our paper:\n",
    "\n",
    "Keemink, S. W., Lowe, S. C., Pakan, J. M. P., Dylda, E., van Rossum, M. C. W., and Rochefort, N. L. FISSA: A neuropil decontamination toolbox for calcium imaging signals, *Scientific Reports*, **8**(1):3493, 2018. doi: [10.1038/s41598-018-21640-2](https://www.doi.org/10.1038/s41598-018-21640-2).\n",
    "\n",
    "See [basic_usage_func.py](https://github.com/rochefort-lab/fissa/blob/1.0.0/examples/basic_usage_func.py) (or [basic_usage_func_windows.py](https://github.com/rochefort-lab/fissa/blob/1.0.0/examples/basic_usage_func_windows.py) for Windows users) for a short example script outside of a notebook interface."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import packages\n",
    "\n",
    "First, we need to import fissa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the FISSA toolbox\n",
    "import fissa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to import some plotting dependencies which we'll make use in this notebook to display the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For plotting our results, import numpy and matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fetch the colormap object for Cynthia Brewer's Paired color scheme\n",
    "colors = plt.get_cmap(\"Paired\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running FISSA\n",
    "\n",
    "With the function-based interface to FISSA, everything is handled in a single function call to [fissa.run_fissa](https://fissa.readthedocs.io/en/1.0.0/source/packages/fissa.core.html#fissa.core.run_fissa). The function takes as its input is the raw signals, and returns the decontaminated signals.\n",
    "\n",
    "The mandatory inputs to [fissa.run_fissa](https://fissa.readthedocs.io/en/1.0.0/source/packages/fissa.core.html#fissa.core.run_fissa) are:\n",
    "\n",
    "- the experiment images\n",
    "- the regions of interest (ROIs) to extract\n",
    "\n",
    "Images can be given as a path to a folder containing tiff stacks:\n",
    "```python\n",
    "images = \"folder\"\n",
    "```\n",
    "Each of these tiff-stacks in the folder (e.g. `\"folder/trial_001.tif\"`) is a trial with many frames.\n",
    "Although we refer to one trial as an `image`, it is actually a video recording.\n",
    "\n",
    "Alternatively, the image data can be given as a list of paths to tiffs:\n",
    "```python\n",
    "images = [\"folder/trial_001.tif\", \"folder/trial_002.tif\", \"folder/trial_003.tif\"]\n",
    "```\n",
    "or as a list of arrays which you have already loaded into memory:\n",
    "```python\n",
    "images = [array1, array2, array3, ...]\n",
    "```\n",
    "\n",
    "For the regions of interest (ROIs) input, you can either provide a single set of ROIs, or a set of ROIs for every image.\n",
    "\n",
    "If the ROIs were defined using ImageJ, use ImageJ's export function to save them in a zip.\n",
    "Then, provide the ROI filename.\n",
    "```python\n",
    "rois = \"rois.zip\"  # for a single set of ROIs used across all images\n",
    "```\n",
    "The same set of ROIs will be used for every image in `images`.\n",
    "\n",
    "Sometimes there is motion between trials causing the alignment of the ROIs to drift.\n",
    "In such a situation, you may need to use a slightly different location of the ROIs for each trial.\n",
    "This can be handled by providing FISSA with a list of ROI sets — one ROI set (i.e. one ImageJ zip file) per trial.\n",
    "```python\n",
    "rois = [\"rois1.zip\", \"rois2.zip\", ...]  # for a unique roiset for each image\n",
    "```\n",
    "Please note that the ROIs defined in each ROI set must correspond to the same physical reigons across all trials, and that the order must be consistent.\n",
    "That is to say, the 1st ROI listed in each ROI set must correspond to the same item appearing in each trial, etc.\n",
    "\n",
    "In this notebook, we will demonstrate how to use FISSA with ImageJ ROI sets, saved as zip files.\n",
    "However, you are not restricted to providing your ROIs to FISSA in this format.\n",
    "FISSA will also accept ROIs which are arbitrarily defined by providing them as arrays (`numpy.ndarray` objects).\n",
    "ROIs provided in this way can be defined either as boolean-valued masks indicating the presence of a ROI per-pixel in the image, or defined as a list of coordinates defining the boundary of the ROI.\n",
    "For examples of such usage, see our [Suite2p](https://fissa.readthedocs.io/en/1.0.0/examples/Suite2p%20example.html), [CNMF](https://fissa.readthedocs.io/en/1.0.0/examples/cNMF%20example.html), and [SIMA](https://fissa.readthedocs.io/en/1.0.0/examples/SIMA%20example.html) example notebooks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example, we will run FISSA on a small test dataset.\n",
    "\n",
    "The test dataset can be found and downloaded from [the examples folder of the fissa repository](https://github.com/rochefort-lab/fissa/tree/1.0.0/examples), along with the source for this example notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define path to imagery and to the ROI set\n",
    "images_location = \"exampleData/20150529\"\n",
    "rois_location = \"exampleData/20150429.zip\"\n",
    "\n",
    "# Call FISSA using the function-based interface\n",
    "result, raw = fissa.run_fissa(images_location, rois_location)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function-based interface is very straight forward, but note that you can only access the result which is returned by the function.\n",
    "\n",
    "If you need to access the raw traces, ROI masks, or demixing matrix, you need to use the more flexible object-oriented (class based) interface using [fissa.Experiment](https://fissa.readthedocs.io/en/1.0.0/source/packages/fissa.core.html#fissa.core.Experiment) instead.\n",
    "An example of this is given in our [object-oriented example usage notebook](https://rochefort-lab.github.io/fissa/examples/Basic%20usage.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with results\n",
    "\n",
    "The output of ``fissa.run_fissa`` is structured as a 2-d array of 2-d arrays (it can't be a 4-d array because of trials generally don't have the same number of frames).\n",
    "\n",
    "The results from the cell (ROI) numbered `c` and the trial (TIFF) numbered `t` are located at `result[c, t][0, :]`.\n",
    "\n",
    "The fourth and final dimension works through frames within the TIFF file (time).\n",
    "\n",
    "The third dimension iterates over output signals.\n",
    "The 0-th entry of this is the signal which most closely corresponds to the raw signal within the ROI, and is FISSA's best guess for the decontaminated cell source signal.\n",
    "The other signals are the isolated signals from contaminants such as neuropil and neighbouring cells.\n",
    "\n",
    "Let's compare the raw signal to the separated signal for a single trial from an example ROI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot sample trace\n",
    "\n",
    "# Select the ROI and trial to plot\n",
    "roi = 2\n",
    "trial = 1\n",
    "\n",
    "# Create the figure\n",
    "plt.figure(figsize=(12, 6))\n",
    "\n",
    "plt.plot(\n",
    "    raw[roi, trial][0, :],\n",
    "    lw=2,\n",
    "    label=\"Raw\",\n",
    "    color=colors((roi * 2) % colors.N),\n",
    ")\n",
    "plt.plot(\n",
    "    result[roi, trial][0, :],\n",
    "    lw=2,\n",
    "    label=\"Decontaminated\",\n",
    "    color=colors((roi * 2 + 1) % colors.N),\n",
    ")\n",
    "\n",
    "plt.title(\"ROI {},  Trial {}\".format(roi, trial), fontsize=15)\n",
    "plt.xlabel(\"Time (frame number)\", fontsize=15)\n",
    "plt.ylabel(\"Signal intensity (candela per unit area)\", fontsize=15)\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's next plot the traces across all ROIs and trials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot all ROIs and trials\n",
    "\n",
    "# Get the number of ROIs and trials\n",
    "n_roi = result.shape[0]\n",
    "n_trial = result.shape[1]\n",
    "\n",
    "# Find the maximum signal intensities for each ROI\n",
    "roi_max_raw = [\n",
    "    np.max([np.max(raw[i_roi, i_trial][0]) for i_trial in range(n_trial)])\n",
    "    for i_roi in range(n_roi)\n",
    "]\n",
    "roi_max_result = [\n",
    "    np.max([np.max(result[i_roi, i_trial][0]) for i_trial in range(n_trial)])\n",
    "    for i_roi in range(n_roi)\n",
    "]\n",
    "roi_max = np.maximum(roi_max_raw, roi_max_result)\n",
    "\n",
    "# Plot our figure using subplot panels\n",
    "plt.figure(figsize=(16, 10))\n",
    "for i_roi in range(n_roi):\n",
    "    for i_trial in range(n_trial):\n",
    "        # Make subplot axes\n",
    "        i_subplot = 1 + i_trial * n_roi + i_roi\n",
    "        plt.subplot(n_trial, n_roi, i_subplot)\n",
    "        # Plot the data\n",
    "        plt.plot(\n",
    "            raw[i_roi][i_trial][0, :],\n",
    "            label=\"Raw\",\n",
    "            color=colors((i_roi * 2) % colors.N),\n",
    "        )\n",
    "        plt.plot(\n",
    "            result[i_roi][i_trial][0, :],\n",
    "            label=\"Decontaminated\",\n",
    "            color=colors((i_roi * 2 + 1) % colors.N),\n",
    "        )\n",
    "        # Labels and boiler plate\n",
    "        plt.ylim([-0.05 * roi_max[i_roi], roi_max[i_roi] * 1.05])\n",
    "        if i_roi == 0:\n",
    "            plt.ylabel(\n",
    "                \"Trial {}\\n\\nSignal intensity\\n(candela per unit area)\".format(\n",
    "                    i_trial + 1\n",
    "                )\n",
    "            )\n",
    "        if i_trial == 0:\n",
    "            plt.title(\"ROI {}\".format(i_roi))\n",
    "            plt.legend()\n",
    "        if i_trial == n_trial - 1:\n",
    "            plt.xlabel(\"Time (frame number)\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing ROI signal to neuropil region signals\n",
    "It can be very instructive to compare the signal in the central ROI with the surrounding neuropil regions. These can be found for cell `c` and trial `t` in `raw[c, t][i, :]`, with `i=0` being the cell, and `i=1,2,3,...` indicating the surrounding regions.\n",
    "\n",
    "Below we compare directly the raw ROI trace, the decontaminated trace, and the surrounding neuropil region traces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the number of neuropil/surrounding regions.\n",
    "# The raw data has the raw ROI signal in raw[:, :][0] and raw surround\n",
    "# signals in the rest of the 3rd dimension.\n",
    "nRegions = raw[0, 0].shape[0] - 1\n",
    "\n",
    "# Select the ROI and trial to plot\n",
    "roi = 2\n",
    "trial = 1\n",
    "\n",
    "# Create the figure\n",
    "plt.figure(figsize=(12, 12))\n",
    "\n",
    "# Plot extracted traces for each neuropil subregion\n",
    "plt.subplot(2, 1, 1)\n",
    "# Plot trace of raw ROI signal\n",
    "plt.plot(\n",
    "    raw[roi, trial][0, :],\n",
    "    lw=2,\n",
    "    label=\"Raw ROI signal\",\n",
    "    color=colors((roi * 2) % colors.N),\n",
    ")\n",
    "# Plot traces from each neuropil region\n",
    "for i_neuropil in range(1, nRegions + 1):\n",
    "    alpha = i_neuropil / nRegions\n",
    "    plt.plot(\n",
    "        raw[roi, trial][i_neuropil, :],\n",
    "        lw=2,\n",
    "        label=\"Neuropil region {}\".format(i_neuropil),\n",
    "        color=\"k\",\n",
    "        alpha=alpha,\n",
    "    )\n",
    "plt.ylim([0, 125])\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.ylabel(\"Signal intensity (candela per unit area)\", fontsize=15)\n",
    "plt.title(\"ROI {},  Trial {}, neuropil region traces\".format(roi, trial), fontsize=15)\n",
    "\n",
    "# Plot the ROI signal\n",
    "plt.subplot(2, 1, 2)\n",
    "# Plot trace of raw ROI signal\n",
    "plt.plot(raw[roi, trial][0, :], lw=2, label=\"Raw\", color=colors((roi * 2) % colors.N))\n",
    "# Plot decontaminated signal matched to the ROI\n",
    "plt.plot(\n",
    "    result[roi, trial][0, :],\n",
    "    lw=2,\n",
    "    label=\"Decontaminated\",\n",
    "    color=colors((roi * 2 + 1) % colors.N),\n",
    ")\n",
    "\n",
    "plt.ylim([0, 125])\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.xlabel(\"Time (frame number)\", fontsize=15)\n",
    "plt.ylabel(\"Signal intensity (candela per unit area)\", fontsize=15)\n",
    "plt.title(\"ROI {},  Trial {}, raw and decontaminated\".format(roi, trial), fontsize=15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### df/f<sub>0</sub>\n",
    "\n",
    "The default output from `fissa.run_fissa` is in the same units as the raw input (candelas per unit area).\n",
    "\n",
    "It is often desirable to calculate the intensity of a signal relative to the baseline value, df/f<sub>0</sub>, for the traces.\n",
    "`fissa.run_fissa` will do this for you provide the argument `return_deltaf=True`, and the sampling frequency of your TIFF files with `freq=sample_frequency`.\n",
    "The sampling frequency must be provided because the data is smoothed in order to determine the baseline value f<sub>0</sub>.\n",
    "\n",
    "When `return_deltaf=True`, `run_fissa` will return the df/f<sub>0</sub> output *instead* of the source signal traces scaled at the recording intensity.\n",
    "If you need to access both the standard FISSA output *and* the df/f<sub>0</sub> output at the same time, you need to use the more flexible [fissa.Experiment](https://fissa.readthedocs.io/en/1.0.0/source/packages/fissa.core.html#fissa.core.Experiment) FISSA interface instead, as described in [this example](https://fissa.readthedocs.io/en/1.0.0/examples/Basic%20usage.html#df/f0)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_frequency = 10  # Hz\n",
    "\n",
    "deltaf_result, deltaf_raw = fissa.run_fissa(\n",
    "    images_location, rois_location, freq=sample_frequency, return_deltaf=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that by default, f<sub>0</sub> is determined as the minimum across all trials (all TIFFs) to ensure that results are directly comparable between trials, but you can normalise each trial individually instead if you prefer by providing the parameter `deltaf_across_trials=False`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot sample trace\n",
    "\n",
    "# Select the ROI and trial to plot\n",
    "roi = 2\n",
    "trial = 1\n",
    "\n",
    "# Create the figure\n",
    "plt.figure(figsize=(12, 6))\n",
    "\n",
    "plt.plot(\n",
    "    deltaf_raw[roi, trial][0, :],\n",
    "    lw=2,\n",
    "    label=\"Raw\",\n",
    "    color=colors((roi * 2) % colors.N),\n",
    ")\n",
    "plt.plot(\n",
    "    deltaf_result[roi, trial][0, :],\n",
    "    lw=2,\n",
    "    label=\"Decontaminated\",\n",
    "    color=colors((roi * 2 + 1) % colors.N),\n",
    ")\n",
    "\n",
    "plt.title(\"ROI {},  Trial {}\".format(roi, trial), fontsize=15)\n",
    "plt.xlabel(\"Time (frame number)\", fontsize=15)\n",
    "plt.ylabel(r\"$\\Delta f\\,/\\,f_0$\", fontsize=15)\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since FISSA is very good at removing contamination from the ROI signals, the minimum value on the decontaminated trace will typically be `0.`.\n",
    "Consequently, we use the minimum value of the (smoothed) raw signal to provide the f<sub>0</sub> from the raw trace for both the raw and decontaminated df/f<sub>0</sub>.\n",
    "\n",
    "We can plot the df/f<sub>0</sub> for every cell during every trial as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the number of ROIs and trials\n",
    "n_roi = result.shape[0]\n",
    "n_trial = result.shape[1]\n",
    "\n",
    "# Find the maximum signal intensities for each ROI,\n",
    "# so we can scale ylim the same across subplots\n",
    "roi_max = [\n",
    "    np.max([np.max(result[i_roi, i_trial][0]) for i_trial in range(n_trial)])\n",
    "    for i_roi in range(n_roi)\n",
    "]\n",
    "\n",
    "# Plot our figure using subplot panels\n",
    "plt.figure(figsize=(16, 10))\n",
    "for i_roi in range(n_roi):\n",
    "    for i_trial in range(n_trial):\n",
    "        # Make subplot axes\n",
    "        i_subplot = 1 + i_trial * n_roi + i_roi\n",
    "        plt.subplot(n_trial, n_roi, i_subplot)\n",
    "        # Plot the data\n",
    "        plt.plot(\n",
    "            result[i_roi][i_trial][0, :],\n",
    "            color=colors((i_roi * 2 + 1) % colors.N),\n",
    "        )\n",
    "        # Labels and boiler plate\n",
    "        plt.ylim([-0.05 * roi_max[i_roi], roi_max[i_roi] * 1.05])\n",
    "        if i_roi == 0:\n",
    "            plt.ylabel(\"Trial {}\\n\\n\".format(i_trial + 1) + r\"$\\Delta f\\,/\\,f_0$\")\n",
    "        if i_trial == 0:\n",
    "            plt.title(\"ROI {}\".format(i_roi))\n",
    "        if i_trial == n_trial - 1:\n",
    "            plt.xlabel(\"Time (frame number)\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For comparison purposes, we can also plot the df/f<sub>0</sub> for the raw data against the decontaminated signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot all ROIs and trials\n",
    "\n",
    "# Get the number of ROIs and trials\n",
    "n_roi = deltaf_result.shape[0]\n",
    "n_trial = deltaf_result.shape[1]\n",
    "\n",
    "# Find the maximum signal intensities for each ROI\n",
    "roi_max_raw = [\n",
    "    np.max([np.max(deltaf_raw[i_roi, i_trial][0]) for i_trial in range(n_trial)])\n",
    "    for i_roi in range(n_roi)\n",
    "]\n",
    "roi_max_result = [\n",
    "    np.max([np.max(deltaf_result[i_roi, i_trial][0]) for i_trial in range(n_trial)])\n",
    "    for i_roi in range(n_roi)\n",
    "]\n",
    "roi_max = np.maximum(roi_max_raw, roi_max_result)\n",
    "\n",
    "# Plot our figure using subplot panels\n",
    "plt.figure(figsize=(16, 10))\n",
    "for i_roi in range(n_roi):\n",
    "    for i_trial in range(n_trial):\n",
    "        # Make subplot axes\n",
    "        i_subplot = 1 + i_trial * n_roi + i_roi\n",
    "        plt.subplot(n_trial, n_roi, i_subplot)\n",
    "        # Plot the data\n",
    "        plt.plot(\n",
    "            deltaf_raw[i_roi][i_trial][0, :],\n",
    "            label=\"Raw\",\n",
    "            color=colors((i_roi * 2) % colors.N),\n",
    "        )\n",
    "        plt.plot(\n",
    "            deltaf_result[i_roi][i_trial][0, :],\n",
    "            label=\"Decontaminated\",\n",
    "            color=colors((i_roi * 2 + 1) % colors.N),\n",
    "        )\n",
    "        # Labels and boiler plate\n",
    "        plt.ylim([-0.05 * roi_max[i_roi], roi_max[i_roi] * 1.05])\n",
    "        if i_roi == 0:\n",
    "            plt.ylabel(\"Trial {}\\n\\n\".format(i_trial + 1) + r\"$\\Delta f\\,/\\,f_0$\")\n",
    "        if i_trial == 0:\n",
    "            plt.title(\"ROI {}\".format(i_roi))\n",
    "            plt.legend()\n",
    "        if i_trial == n_trial - 1:\n",
    "            plt.xlabel(\"Time (frame number)\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Caching\n",
    "\n",
    "After using FISSA to decontaminate the data collected in an experiment, you will probably want to save the output for later use, so you don't have to keep re-running FISSA on the data.\n",
    "\n",
    "To facilitate this, an option to cache the outputs is built into FISSA.\n",
    "If you provide ``fissa.run_fissa`` with an identifier to the experiment being analysed in the ``folder`` argument, it will cache results into the corresponding directory.\n",
    "Later, if you call ``fissa.run_fissa`` again with the same ``folder`` argument, it will load the saved results from that cache folder instead of recomputing them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the folder where FISSA's outputs will be cached, so they can be\n",
    "# quickly reloaded in the future without having to recompute them.\n",
    "#\n",
    "# This argument is optional; if it is not provided, FISSA will not save its\n",
    "# results for later use.\n",
    "#\n",
    "# If the output directory already exists, FISSA will load the contents of\n",
    "# the cache instead of recomputing it.\n",
    "#\n",
    "# Note: you *must* use a different folder for each experiment, otherwise\n",
    "# FISSA will load the existing data instead of computing results for the\n",
    "# new experiment.\n",
    "#\n",
    "# In this example, we will use the current datetime as the  name of the\n",
    "# experiment, but you can name your experiments however you want to.\n",
    "# If you want to take advantage of the caching of results, you should use\n",
    "# a more descriptive name than this so you can identify the actual\n",
    "# dataset that the FISSA results correspond to, and load them appropriately.\n",
    "\n",
    "import datetime\n",
    "\n",
    "output_folder = \"fissa-example_{}\".format(\n",
    "    datetime.datetime.now().strftime(\"%Y-%m-%d_%H-%M-%S\")\n",
    ")\n",
    "\n",
    "print(output_folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's run FISSA on this experiment again, but this time save the results to the experiment's output directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run FISSA, saving to results to output_folder\n",
    "result, raw = fissa.run_fissa(images_location, rois_location, folder=output_folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A subsequent call to ``fissa.run_fissa`` which uses the same experiment folder argument will load the cached data instead of re-running the FISSA signal separation routine from scratch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run FISSA, loading results from cache in output_folder\n",
    "result, raw = fissa.run_fissa(images_location, rois_location, folder=output_folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exporting to MATLAB\n",
    "\n",
    "The results can easily be exported to a MATLAB-compatible [MAT-file](https://mathworks.com/help/matlab/import_export/mat-file-versions.html) as follows.\n",
    "\n",
    "If we provide `export_to_matfile=True` to `fissa.run_fissa`, it will export the data a matfile named `\"separated.mat\"` within the cache directory (the cache directory as provided with the `folder` argument)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result, raw = fissa.run_fissa(\n",
    "    images_location, rois_location, folder=output_folder, export_to_matfile=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can export to a matfile with a custom file name by setting the `export_to_matfile` argument to the target path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result, raw = fissa.run_fissa(\n",
    "    images_location, rois_location, export_to_matfile=\"experiment_results.mat\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loading the generated file (e.g. `\"output_folder/separated.mat\"`) in MATLAB will provide you with all of FISSA's outputs.\n",
    "\n",
    "These are structured in the same way as the `raw` and `result` variables returned by `fissa.run_fissa`.\n",
    "With the python interface, the outputs are 2d numpy.ndarrays each element of which is itself a 2d numpy.ndarrays.\n",
    "Meanwhile, when the output is loaded into MATLAB the data is structured as a 2d cell-array each element of which is a 2d matrix.\n",
    "\n",
    "Additionally, note that whilst Python indexes from 0, MATLAB indexes from 1 instead.\n",
    "As a consequence of this, the results seen on Python for a given roi and trial `experiment.result[roi, trial]` correspond to the index `S.result{roi + 1, trial + 1}` on MATLAB.\n",
    "\n",
    "Our first plot in this notebook can be replicated in MATLAB as follows:\n",
    "```octave\n",
    "%% Plot example traces in MATLAB\n",
    "% Load FISSA output data in MATLAB\n",
    "% ... either from the automatic file name within a cache\n",
    "% S = load('fissa-example/separated.mat')\n",
    "% ... or from a custom output path\n",
    "S = load('experiment_results.mat')\n",
    "% Select the second trial\n",
    "% (On Python, this would be trial = 1)\n",
    "trial = 2;\n",
    "% Plot the result traces for each ROI\n",
    "figure;\n",
    "hold on;\n",
    "for i_roi = 1:size(S.result, 1);\n",
    "    plot(S.result{i_roi, trial}(1, :));\n",
    "end\n",
    "xlabel('Time (frame number)');\n",
    "ylabel('Signal intensity (candela per unit area)');\n",
    "grid on;\n",
    "box on;\n",
    "set(gca,'TickDir','out');\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Customisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Controlling verbosity\n",
    "\n",
    "The level of verbosity of FISSA can be controlled with the `verbosity` parameter.\n",
    "\n",
    "The default is `verbosity=1`.\n",
    "\n",
    "If the verbosity parameter is higher, FISSA will print out more information while it is processing.\n",
    "This can be helpful for debugging puproses.\n",
    "The verbosity reaches its maximum at `verbosity=6`.\n",
    "\n",
    "If `verbosity=0`, FISSA will run silently."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Call FISSA with elevated verbosity\n",
    "result = fissa.run_fissa(images_location, rois_location, verbosity=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis parameters\n",
    "\n",
    "FISSA has several user-definable settings, which can be set as optional arguments to `fissa.run_fissa`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# FISSA uses multiprocessing to speed up its processing.\n",
    "# By default, it will spawn one worker per CPU core on your machine.\n",
    "# However, if you have a lot of cores and not much memory, you many not\n",
    "# be able to suport so many workers simultaneously.\n",
    "# In particular, this can be problematic during the data preparation step\n",
    "# in which TIFFs are loaded into memory.\n",
    "# The default number of cores for the data preparation and separation steps\n",
    "# can be changed as follows.\n",
    "ncores_preparation = 4  # If None, uses all available cores\n",
    "ncores_separation = None  # if None, uses all available cores\n",
    "\n",
    "# By default, FISSA uses 4 subregions for the neuropil region.\n",
    "# If you have very dense data with a lot of different signals per unit area,\n",
    "# you may wish to increase the number of regions.\n",
    "n_regions = 8\n",
    "\n",
    "# By default, each surrounding region has the same area as the central ROI.\n",
    "# i.e. expansion = 1\n",
    "# However, you may wish to increase or decrease this value.\n",
    "expansion = 0.75\n",
    "\n",
    "# The degree of signal sparsity can be controlled with the alpha parameter.\n",
    "alpha = 0.02\n",
    "\n",
    "# If you change the experiment parameters, you need to change the cache directory too.\n",
    "# Otherwise FISSA will try to reload the results from the previous run instead of\n",
    "# computing the new results. FISSA will throw an error if you try to load data which\n",
    "# was generated with different analysis parameters to the current parameters.\n",
    "output_folder2 = output_folder + \"_alt\"\n",
    "\n",
    "# Run FISSA with these parameters\n",
    "result, raw = fissa.run_fissa(\n",
    "    images_location,\n",
    "    rois_location,\n",
    "    output_folder2,\n",
    "    nRegions=n_regions,\n",
    "    expansion=expansion,\n",
    "    alpha=alpha,\n",
    "    ncores_preparation=ncores_preparation,\n",
    "    ncores_separation=ncores_separation,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the new results for our example trace from before. Although we doubled the number of neuropil regions around the cell, very little has changed for this example because there were not many sources of contamination.\n",
    "\n",
    "However, there will be more of a difference if your data has more neuropil sources per unit area within the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_roi = result.shape[0]\n",
    "n_trial = result.shape[1]\n",
    "\n",
    "i_roi = 3\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "\n",
    "for i_trial in range(n_trial):\n",
    "    plt.plot(result[i_roi, i_trial][0, :], label=\"Trial {}\".format(i_trial + 1))\n",
    "\n",
    "plt.title(\"ROI {}\".format(i_roi), fontsize=15)\n",
    "plt.xlabel(\"Time (Frame number)\", fontsize=15)\n",
    "plt.ylabel(\"Signal intensity (candela per unit area)\", fontsize=15)\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with very large tiff files\n",
    "\n",
    "By default, FISSA loads entire TIFF files into memory at once and then manipulates all ROIs within the TIFF.\n",
    "This is more efficient, but can be problematic when working with very large TIFF files which are too big to be loaded into memory all at once.\n",
    "\n",
    "If you run out of memory when running FISSA, you can try reducing the number of workers during the memory-intensive preparation step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = fissa.run_fissa(images_location, rois_location, ncores_preparation=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, you can activate FISSA's low memory mode.\n",
    "In this mode, it will load and process frames one at a time.\n",
    "This will run a fair bit slower than the default mode, but has a much lower memory requirement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result, raw = fissa.run_fissa(images_location, rois_location, lowmemory_mode=True)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python",
   "language": "python",
   "name": "python"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}