# -*- coding: utf-8 -*-
"""
Main FISSA user interface.
Authors:
- Sander W Keemink <swkeemink@scimail.eu>
- Scott C Lowe <scott.code.lowe@gmail.com>
"""
from __future__ import print_function
import collections
import datetime
import functools
import glob
import itertools
import os.path
import sys
import time
import warnings
try:
import collections.abc as cabc
except ImportError:
import collections as cabc
import numpy as np
from joblib import Parallel, delayed
from past.builtins import basestring
from scipy.io import savemat
from tqdm.auto import tqdm
from . import deltaf, extraction
from . import neuropil as npil
from . import roitools
def _pretty_timedelta(td=None, **kwargs):
"""
Represent a difference in time as a human-readable string.
Parameters
----------
td : datetime.timedelta, optional
The amount of time elapsed.
**kwargs
Additional arguments as per :class:`datetime.timedelta` constructor.
Returns
-------
str
Representation of the amount of time elapsed.
"""
if td is None:
td = datetime.timedelta(**kwargs)
elif not isinstance(td, datetime.timedelta):
raise ValueError(
"First argument should be a datetime.timedelta instance,"
" but {} was given.".format(type(td))
)
elif kwargs:
raise ValueError(
"Either a timedelta object or its arguments should be given, not both."
)
if td.total_seconds() < 2:
return "{:.3f} seconds".format(td.total_seconds())
if td.total_seconds() < 10:
return "{:.2f} seconds".format(td.total_seconds())
if td.total_seconds() < 60:
return "{:.1f} seconds".format(td.total_seconds())
if td.total_seconds() < 3600:
s = td.total_seconds()
m = int(s // 60)
s -= m * 60
return "{:d} min, {:.0f} sec".format(m, s)
# For durations longer than one hour, we use the default string
# representation for a datetime.timedelta, H:MM:SS.microseconds
return str(td)
[docs]def separate_trials(
raw,
alpha=0.1,
max_iter=20000,
tol=1e-4,
max_tries=1,
method="nmf",
verbosity=1,
label=None,
total=None,
):
r"""
Separate signals within a set of 2d arrays.
.. versionadded:: 1.0.0
Parameters
----------
raw : list of n_trials :term:`array_like`, each shaped ``(nRegions + 1, observations)``
Raw signals.
A list of 2-d arrays, each of which contains observations of mixed
signals, mixed in the same way across all trials.
The `nRegions` signals must be the same for each trial, and the 0-th
region, ``raw[trial][0]``, should be from the region of interest for
which a matching source signal should be identified.
alpha : float, default=0.1
Sparsity regularizaton weight for NMF algorithm. Set to zero to
remove regularization. Default is ``0.1``.
(Only used for ``method="nmf"``.)
max_iter : int, default=20000
Maximum number of iterations before timing out on an attempt.
tol : float, default=1e-4
Tolerance of the stopping condition.
max_tries : int, default=1
Maximum number of random initial states to try. Each random state will
be optimized for `max_iter` iterations before timing out.
method : {"nmf", "ica"}, default="nmf"
Which blind source-separation method to use. Either ``"nmf"``
for non-negative matrix factorization, or ``"ica"`` for
independent component analysis. Default is ``"nmf"``.
verbosity : int, default=1
Level of verbosity. The options are:
- ``0``: No outputs.
- ``1``: Print separation start.
- ``2``: Print separation end.
- ``3``: Print progress details during separation.
label : str or int, optional
Label/name or index of the ROI currently being processed.
Only used for progress messages.
total : int, optional
Total number of ROIs. Only used for reporting progress.
Returns
-------
Xsep : list of n_trials :class:`numpy.ndarray`, each shaped ``(nRegions + 1, observations)``
The separated signals, unordered.
Xmatch : list of n_trials :class:`numpy.ndarray`, each shaped ``(nRegions + 1, observations)``
The separated traces, ordered by matching score against the raw ROI
signal.
Xmixmat : :class:`numpy.ndarray`, shaped ``(nRegions + 1, nRegions + 1)``
Mixing matrix.
convergence : dict
Metadata for the convergence result, with the following keys and
values:
converged : bool
Whether the separation model converged, or if it ended due to
reaching the maximum number of iterations.
iterations : int
The number of iterations which were needed for the separation model
to converge.
max_iterations : int
Maximum number of iterations to use when fitting the
separation model.
random_state : int or None
Random seed used to initialise the separation model.
"""
# Get the timestamp for program start
t0 = time.time()
header = ""
if verbosity >= 1:
# Set up message header
# Use the label, if this was provided
if label is None:
header = ""
elif isinstance(label, int) and isinstance(total, int):
# Pad left based on the total number of jobs, so it is [ 1/10] etc
fmtstr = "{:" + str(int(np.maximum(1, np.ceil(np.log10(total))))) + "d}"
header = fmtstr.format(label + 1)
else:
header = str(label)
# Try to label with [1/5] to indicate progess, if possible
if header and total is not None:
header += "/{}".format(total)
if header:
header = "[Separation " + header + "] "
# Include the ROI label as a footer
footer = ""
if isinstance(label, int) and isinstance(total, int):
# Include the ROI label as a footer
footer = " (ROI {})".format(label)
# Done with header and footer
# Build intro message
message = header + "Signal separation starting" + footer
# Wait briefly to prevent messages colliding when using multiprocessing
if isinstance(label, int) and label < 12:
time.sleep(label / 50.0)
print(message)
sys.stdout.flush()
# Join together the raw data across trials, collapsing down the trials
X = np.concatenate(raw, axis=1)
# Check for values below 0
if X.min() < 0:
message_extra = ""
if label is not None:
message_extra = " for ROI {}".format(label)
warnings.warn(
"{}Found values below zero in raw signal{}. Offsetting so minimum is 0."
"".format(header, message_extra)
)
X -= X.min()
# Separate the signals
Xsep, Xmatch, Xmixmat, convergence = npil.separate(
X,
method,
max_iter=max_iter,
tol=tol,
max_tries=max_tries,
alpha=alpha,
verbosity=verbosity - 2,
prefix=" " + header,
)
# Unravel observations from multiple trials into a list of arrays
trial_lengths = [r.shape[1] for r in raw]
indices = np.cumsum(trial_lengths[:-1])
Xsep = np.split(Xsep, indices, axis=1)
Xmatch = np.split(Xmatch, indices, axis=1)
# Report status
if verbosity >= 2:
# Build end message
message = header + "Signal separation finished" + footer
message += " in {}".format(_pretty_timedelta(seconds=time.time() - t0))
print(message)
sys.stdout.flush()
return Xsep, Xmatch, Xmixmat, convergence
[docs]class Experiment:
r"""
Class-based interface for running FISSA on experiments.
Uses the methodology described in
`FISSA: A neuropil decontamination toolbox for calcium imaging signals <doi_>`_.
.. _doi: https://www.doi.org/10.1038/s41598-018-21640-2
Parameters
----------
images : str or list
The raw imaging data.
Should be one of:
- the path to a directory containing TIFF files (string),
- a list of paths to TIFF files (list of strings),
- a list of :term:`array_like` data already loaded into memory,
each shaped ``(n_frames, height, width)``.
Note that each TIFF or array is considered a single trial.
rois : str or list
The region of interest (ROI) definitions.
Should be one of:
- the path to a directory containing ImageJ ZIP files (string),
- the path of a single ImageJ ZIP file (string),
- a list of ImageJ ZIP files (list of strings),
- a list of arrays, each encoding a ROI polygons,
- a list of lists of binary arrays, each representing a ROI mask.
This can either be a single roiset for all trials, or a different
roiset for each trial.
folder : str, optional
Path to a cache directory from which pre-extracted data will
be loaded if present, and saved to otherwise. If `folder` is
unset, the experiment data will not be saved.
nRegions : int, default=4
Number of neuropil regions and signals to use. Default is ``4``.
Use a higher number for densely labelled tissue.
expansion : float, default=1
Expansion factor for each neuropil region, relative to the
ROI area. Default is ``1``. The total neuropil area will be
``nRegions * expansion * area(ROI)``.
method : "nmf" or "ica", default="nmf"
Which blind source-separation method to use. Either ``"nmf"``
for non-negative matrix factorization, or ``"ica"`` for
independent component analysis. Default is ``"nmf"`` (recommended).
alpha : float, default=0.1
Sparsity regularizaton weight for NMF algorithm. Set to zero to
remove regularization. Default is ``0.1``.
max_iter : int, default=20000
Maximum number of iterations before timing out on an attempt.
.. versionadded:: 1.0.0
tol : float, default=1e-4
Tolerance of the stopping condition.
.. versionadded:: 1.0.0
max_tries : int, default=1
Maximum number of random initial states to try. Each random state will
be optimized for `max_iter` iterations before timing out.
.. versionadded:: 1.0.0
ncores_preparation : int or None, default=-1
The number of parallel subprocesses to use during the data
preparation steps of :meth:`separation_prep`.
These steps are ROI and neuropil subregion definitions, and extracting
raw signals from TIFFs.
If set to ``None`` or ``-1`` (default), the number of processes used
will equal the number of threads on the machine.
If this is set to ``-2``, the number of processes used will be one less
than the number of threads on the machine; etc.
Note that the preparation process can be quite memory-intensive and it
may be necessary to reduce the number of processes from the default.
ncores_separation : int or None, default=-1
The number of parallel subprocesses to use during the signal
separation steps of :meth:`separate`.
If set to ``None`` or ``-1`` (default), the number of processes used
will equal the number of threads on the machine.
If this is set to ``-2``, the number of processes used will be one less
than the number of threads on the machine; etc.
The separation routine requires less memory per process than
the preparation routine, and so `ncores_separation` be often be set
higher than `ncores_preparation`.
lowmemory_mode : bool, optional
If ``True``, FISSA will load TIFF files into memory frame-by-frame
instead of holding the entire TIFF in memory at once. This
option reduces the memory load, and may be necessary for very
large inputs. Default is ``False``.
datahandler : :class:`fissa.extraction.DataHandlerAbstract`, optional
A custom datahandler object for handling ROIs and calcium data can
be given here. See :mod:`fissa.extraction` for example datahandler
classes. The default datahandler is
:class:`~fissa.extraction.DataHandlerTifffile`.
If `datahandler` is set, the `lowmemory_mode` parameter is
ignored.
verbosity : int, default=1
How verbose the processing will be. Increase for more output messages.
Processing is silent if ``verbosity=0``.
.. versionadded:: 1.0.0
Attributes
----------
result : :class:`numpy.ndarray`
A :class:`numpy.ndarray` of shape ``(n_rois, n_trials)``, each element
of which is itself a :class:`numpy.ndarray` shaped
``(n_signals, n_timepoints)``.
The final output of FISSA, with separated signals ranked in order of
their weighting toward the raw cell ROI signal relative to their
weighting toward other mixed raw signals.
The ordering is such that ``experiment.result[roi, trial][0, :]``
is the signal with highest score in its contribution to the raw
neuronal signal.
Subsequent signals are sorted in order of diminishing score.
The units are same as `raw` (candelas per unit area).
This field is only populated after :meth:`separate` has been run; until
then, it is set to ``None``.
roi_polys : :class:`numpy.ndarray`
A :class:`numpy.ndarray` of shape ``(n_rois, n_trials)``, each element
of which is itself a list of length ``nRegions + 1``, each element of
which is a list of length ``n_contour`` containing a :class:`numpy.ndarray`
of shape ``(n_nodes, 2)``.
Polygon contours describing the outline of each region.
For contiguous ROIs, the outline of the ``i_roi``-th ROI used in the
``i_trial``-th trial is described by the array at
``experiment.roi_polys[i_roi, i_trial][0][0]``.
This array consists of ``n_nodes`` rows, each representing the
coordinate of a node in ``(y, x)`` format.
For non-contiguous ROIs, a contour is needed for each disconnected
polygon making up the total aggregate ROI. These contours are found at
``experiment.roi_polys[i_roi, i_trial][0][i_contour]``.
Similarly, the `nRegions` neuropil regions are each described by the
polygons
``experiment.roi_polys[i_roi, i_trial][i_neurpil + 1][i_contour]``,
respectively.
means : list of `n_trials` :class:`numpy.ndarray`, each shaped ``(height, width)``
The temporal-mean image for each trial (i.e. for each TIFF file,
the average image over all of its frames).
raw : :class:`numpy.ndarray`
A :class:`numpy.ndarray` of shape ``(n_rois, n_trials)``, each element
of which is itself a :class:`numpy.ndarray` shaped
``(n_signals, n_timepoints)``.
For each ROI and trial (``raw[i_roi, i_trial]``) we extract a temporal
trace of the average value within the spatial area of each of the
``nRegions + 1`` regions.
The 0-th region is the ``i_roi``-th ROI (``raw[i_roi, i_trial][0]``).
The subsequent ``nRegions`` vectors are the traces for each of the
neuropil regions.
The units are the same as the supplied imagery (candelas per unit
area).
sep : :class:`numpy.ndarray`
A :class:`numpy.ndarray` of shape ``(n_rois, n_trials)``, each element
of which is itself a :class:`numpy.ndarray` shaped
``(n_signals, n_timepoints)``.
The separated signals, before output signals are ranked according to
their matching against the raw signal from within the ROI.
Separated signal ``i`` for a specific ROI and trial can be found at
``experiment.sep[roi, trial][i, :]``.
This field is only populated after :meth:`separate` has been run; until
then, it is set to ``None``.
mixmat : :class:`numpy.ndarray`
A :class:`numpy.ndarray` of shape ``(n_rois, n_trials)``, each element
of which is itself a :class:`numpy.ndarray` shaped
``(n_rois, n_signals)``.
The mixing matrix, which maps from ``experiment.raw`` to
``experiment.sep``.
Because we use the collate the traces from all trials to determine
separate the signals, the mixing matrices for a given ROI are the
same across all trials.
This means all ``n_trials`` elements in ``mixmat[i_roi, :]`` are
identical.
This field is only populated after :meth:`separate` has been run; until
then, it is set to ``None``.
info : :class:`numpy.ndarray` shaped ``(n_rois, n_trials)`` of dicts
Information about the separation routine.
Each dictionary in the array has the following fields:
converged : bool
Whether the separation model converged, or if it ended due to
reaching the maximum number of iterations.
iterations : int
The number of iterations which were needed for the separation model
to converge.
max_iterations : int
Maximum number of iterations to use when fitting the
separation model.
random_state : int or None
Random seed used to initialise the separation model.
This field is only populated after :meth:`separate` has been run; until
then, it is set to ``None``.
deltaf_raw : :class:`numpy.ndarray`
A :class:`numpy.ndarray` of shape ``(n_rois, n_trials)``, each element
of which is itself a :class:`numpy.ndarray` shaped ``(1, n_timepoint)``.
The amount of change in fluorence relative to the baseline fluorence
(Δf/f\ :sub:`0`).
This field is only populated after :meth:`calc_deltaf` has been run;
until then, it is set to ``None``.
.. versionchanged:: 1.0.0
The shape of the interior arrays changed from ``(n_timepoint, )``
to ``(1, n_timepoint)``.
deltaf_result : :class:`numpy.ndarray`
A :class:`numpy.ndarray` of shape ``(n_rois, n_trials)``, each element
of which is itself a :class:`numpy.ndarray` shaped
``(n_signals, n_timepoints)``.
The amount of change in fluorence relative to the baseline fluorence
(Δf/f\ :sub:`0`).
By default, the baseline is taken from :attr:`raw` because the
minimum values in :attr:`result` are typically zero.
See :meth:`calc_deltaf` for details.
This field is only populated after :meth:`calc_deltaf` has been run;
until then, it is set to ``None``.
"""
_defaults = {
"nRegions": 4,
"expansion": 1,
"method": "nmf",
"alpha": 0.1,
"max_iter": 20000,
"tol": 1e-4,
"max_tries": 1,
}
_preparation_params = ["nRegions", "expansion"]
_separation_params = ["method", "alpha", "max_iter", "tol", "max_tries"]
_preparation_outputs = ["means", "raw", "roi_polys"]
_separation_outputs = ["info", "mixmat", "result", "sep"]
_deltaf_outputs = ["deltaf_raw", "deltaf_result"]
def __init__(
self,
images,
rois,
folder=None,
nRegions=None,
expansion=None,
method=None,
alpha=None,
max_iter=None,
tol=None,
max_tries=None,
ncores_preparation=-1,
ncores_separation=-1,
lowmemory_mode=False,
datahandler=None,
verbosity=1,
):
# Initialise internal variables
self.clear(verbosity=0)
if isinstance(images, basestring):
self.images = sorted(glob.glob(os.path.join(images, "*.tif*")))
elif isinstance(images, cabc.Sequence):
self.images = images
else:
raise ValueError("images should either be string or list")
if isinstance(rois, basestring):
if rois[-3:] == "zip":
self.rois = [rois] * len(self.images)
else:
self.rois = sorted(glob.glob(os.path.join(rois, "*.zip")))
elif isinstance(rois, cabc.Sequence):
self.rois = rois
if len(rois) == 1: # if only one roiset is specified
self.rois *= len(self.images)
else:
raise ValueError("rois should either be string or list")
if datahandler is not None and lowmemory_mode:
raise ValueError(
"Only one of lowmemory_mode and datahandler should be set."
)
elif lowmemory_mode:
self.datahandler = extraction.DataHandlerTifffileLazy()
else:
self.datahandler = datahandler
# define class variables
self.folder = folder
self.nRegions = nRegions
self.expansion = expansion
self.method = method
self.alpha = alpha
self.max_iter = max_iter
self.tol = tol
self.max_tries = max_tries
self.ncores_preparation = ncores_preparation
self.ncores_separation = ncores_separation
self.verbosity = verbosity
# check if any data already exists
if folder is None:
pass
elif folder and not os.path.exists(folder):
os.makedirs(folder)
else:
self.load()
@property
def nCell(self):
if getattr(self, "result", None) is not None:
return self.result.shape[0]
if getattr(self, "raw", None) is not None:
return self.raw.shape[0]
return None
@property
def nTrials(self):
return len(self.images)
def __setattr__(self, name, value):
def check_same_value():
if not hasattr(self, name):
return False
current = getattr(self, name)
if type(current) is not type(value):
return False
if isinstance(current, np.ndarray):
return np.array_equal(current, value)
return current == value
if getattr(self, name, None) is None:
# No need to clear if the current value is None
pass
elif name in ["images", "rois"]:
if not check_same_value():
self.clear()
elif hasattr(self, "_preparation_params") and name in self._preparation_params:
if not check_same_value():
self.clear()
elif hasattr(self, "_separation_params") and name in self._separation_params:
if not check_same_value():
self.clear_separated()
self.__dict__[name] = value
def __str__(self):
if isinstance(self.images, basestring):
str_images = repr(self.images)
elif isinstance(self.images, cabc.Sequence):
str_images = "<{} of length {}>".format(
self.images.__class__.__name__, len(self.images)
)
else:
str_images = repr(self.images)
if isinstance(self.rois, basestring):
str_rois = repr(self.rois)
elif isinstance(self.rois, cabc.Sequence):
str_rois = "<{} of length {}>".format(
self.rois.__class__.__name__, len(self.rois)
)
else:
str_images = repr(self.rois)
fields = (
["folder"]
+ self._preparation_params
+ self._separation_params
+ ["ncores_preparation", "ncores_separation", "datahandler", "verbosity"]
)
str_parts = [
"{}={}".format(field, repr(getattr(self, field))) for field in fields
]
return "{}.{}(images={}, rois={}, {})".format(
__name__,
self.__class__.__name__,
str_images,
str_rois,
", ".join(str_parts),
)
def __repr__(self):
fields = (
["images", "rois", "folder"]
+ self._preparation_params
+ self._separation_params
+ ["ncores_preparation", "ncores_separation", "datahandler", "verbosity"]
)
repr_parts = [
"{}={}".format(field, repr(getattr(self, field))) for field in fields
]
return "{}.{}({})".format(
__name__, self.__class__.__name__, ", ".join(repr_parts)
)
[docs] def clear(self, verbosity=None):
r"""
Clear prepared data, and all data downstream of prepared data.
.. versionadded:: 1.0.0
Parameters
----------
verbosity : int, optional
Whether to show the data fields which were cleared.
By default, the object's :attr:`verbosity` attribute is used.
"""
if verbosity is None:
verbosity = getattr(self, "verbosity", 1) - 1
keys = self._preparation_outputs + ["deltaf_raw"]
# Wipe outputs
keys_cleared = []
for key in keys:
if getattr(self, key, None) is not None:
keys_cleared.append(key)
setattr(self, key, None)
if verbosity >= 1 and keys_cleared:
print("Cleared {}".format(", ".join(repr(k) for k in keys_cleared)))
# Wipe outputs of separate(), as they no longer match self.raw
self.clear_separated(verbosity=verbosity)
[docs] def clear_separated(self, verbosity=None):
r"""
Clear separated data, and all data downstream of separated data.
.. versionadded:: 1.0.0
Parameters
----------
verbosity : int, optional
Whether to show the data fields which were cleared.
By default, the object's :attr:`verbosity` attribute is used.
"""
if verbosity is None:
verbosity = getattr(self, "verbosity", 1) - 1
keys = self._separation_outputs + ["deltaf_result"]
# Wipe outputs
keys_cleared = []
for key in keys:
if getattr(self, key, None) is not None:
keys_cleared.append(key)
setattr(self, key, None)
if verbosity >= 1 and keys_cleared:
print("Cleared {}".format(", ".join(repr(k) for k in keys_cleared)))
def _adopt_default_parameters(self, only_preparation=False, force=False):
r"""
Adopt default values for unset analysis parameters.
.. versionadded:: 1.0.0
Parameters
----------
only_preparation : bool, optional
Whether to restrict the parameters to only those used for data
extraction during the preparation step. Default is ``False``.
force : bool, optional
If `True`, all parameters will be overridden with default values
even if they had already been set. Default is ``False``.
"""
defaults = self._defaults
if only_preparation:
# Prune down to only the preparation parameters
defaults = {
k: v for k, v in defaults.items() if k in self._preparation_params
}
# Check through each parameter and set unset values from defaults
keys_adopted = []
for key, value in defaults.items():
if getattr(self, key, None) is not None and not force:
continue
setattr(self, key, value)
keys_adopted.append(key)
if self.verbosity >= 5 and keys_adopted:
print(
"Adopted default values for {}".format(
", ".join(repr(k) for k in keys_adopted)
)
)
[docs] def load(self, path=None, force=False, skip_clear=False):
r"""
Load data from cache file in npz format.
.. versionadded:: 1.0.0
Parameters
----------
path : str, optional
Path to cache file (.npz format) or a directory containing
``"prepared.npz"`` and/or ``"separated.npz"`` files.
Default behaviour is to use the :attr:`folder` parameter which was
provided when the object was initialised is used
(``experiment.folder``).
force : bool, optional
Whether to load the cache even if its experiment parameters differ
from the properties of this experiment. Default is ``False``.
skip_clear : bool, optional
Whether to skip clearing values before loading. Default is ``False``.
"""
dynamic_properties = ["nCell", "nTrials"]
ValGroup = collections.namedtuple(
"ValGroup",
["category", "validators", "fields", "clearif", "clearfn"],
)
validation_groups = [
ValGroup(
"prepared",
self._preparation_params,
self._preparation_outputs + ["deltaf_raw"],
["raw"],
self.clear,
),
ValGroup(
"separated",
self._preparation_params + self._separation_params,
self._separation_outputs + ["deltaf_result"],
["result"],
self.clear_separated,
),
]
if path is None:
if self.folder is None:
raise ValueError(
"path must be provided if experiment folder is not defined"
)
path = self.folder
if os.path.isdir(path) or path == "":
for fname in ("prepared.npz", "separated.npz"):
fullfname = os.path.join(path, fname)
if not os.path.exists(fullfname):
continue
self.load(fullfname)
return
if self.verbosity >= 1:
print("Loading data from cache {}".format(path))
cache = np.load(path, allow_pickle=True)
def _unpack_scalar(x):
if np.array_equal(x, None):
return None
if x.ndim == 0:
# Handle loading scalars
return x.item()
return x
if force:
for field in cache.files:
if field in dynamic_properties:
continue
setattr(self, field, _unpack_scalar(cache[field]))
return
set_fields = set()
for category, validators, fields, clearif, clearfn in validation_groups:
valid = True
validation_errors = []
for validator in validators:
if getattr(self, validator, None) is None:
# If the validator is not yet set locally, it is fine to
# overwrite it.
continue
if validator not in cache:
# If the validator is not set in the cache and is set
# locally, we can't verify that the cached data is
# compatible. We don't raise an error for this because the
# contents are probably not this category.
valid = False
break
value = _unpack_scalar(cache[validator])
if value is None:
valid = False
break
if not np.array_equal(getattr(self, validator), value):
# If the validator is set and doesn't match the value in
# the cache, we will raise an error.
validation_errors.append(
" {}: Experiment (ours) {}, Cache (theirs) {}".format(
validator,
getattr(self, validator),
value,
)
)
if len(validation_errors) > 0:
raise ValueError(
"Experiment parameter value(s) in {} do not match the"
" current experiment values:\n{}".format(
path, "\n".join(validation_errors)
)
)
if not valid:
continue
# Check the image and roi size is appropriate
for k in ["raw", "result"]:
if k not in cache.files or np.array_equal(cache[k], None):
continue
if cache[k].shape[1] != self.nTrials:
raise ValueError(
"Data mismatch between {} and our images."
" Cached {} has {} trials, but our Experiment has {}"
" trials.".format(path, k, cache[k].shape[1], self.nTrials)
)
if self.nCell is not None and cache[k].shape[0] != self.nCell:
raise ValueError(
"Data mismatch between {} and our roisets."
" Cached {} has {} ROIs, but our Experiment has {}"
" ROIs.".format(path, k, cache[k].shape[1], self.nCell)
)
# Wipe the values currently held before setting new values
if not skip_clear:
for field in clearif:
if field in cache.files:
clearfn()
break
# All the validators were valid, so we are okay to load the fields
#
# Check to see if there are any fields to load. If not, we won't
# load the validators.
any_field_to_load = False
for field in cache.files:
if field in dynamic_properties:
continue
if field in fields:
any_field_to_load = True
# If we don't have any data to load, no need to set the validators
# or print that we loaded something.
if not any_field_to_load:
continue
# Load all the validators, overwriting our local values if None.
# We do this before loading in the data fields because of automatic
# clear when parameter attributes change.
for validator in validators:
if validator not in cache.files:
continue
value = _unpack_scalar(cache[validator])
if getattr(self, validator, None) is None:
if self.verbosity >= 2:
print(
" Adopting value {}={} from {}".format(
validator, repr(value), path
)
)
setattr(self, validator, value)
set_fields.add(validator)
# Load all the fields
for field in fields:
if field not in cache or field in dynamic_properties:
continue
setattr(self, field, _unpack_scalar(cache[field]))
set_fields.add(field)
if self.verbosity >= 2:
print("Loaded {} data from {}".format(category, path))
# Check there weren't any left over fields in the cache which
# were left unloaded
unset_fields = []
for field in cache.files:
if field in dynamic_properties:
continue
if field not in set_fields:
unset_fields.append(field)
if len(unset_fields) > 0 and self.verbosity >= 1:
print(
"Warning: field(s) {} in {} were not loaded.".format(unset_fields, path)
)
[docs] def separation_prep(self, redo=False):
r"""
Prepare and extract the data to be separated.
For each trial, performs the following steps:
- Load in data as arrays.
- Load in ROIs as masks.
- Grow and seaparate ROIs to define neuropil regions.
- Using neuropil and original ROI regions, extract traces from data.
After running this you can access the raw data (i.e. pre-separation)
as ``experiment.raw`` and ``experiment.rois``.
``experiment.raw`` is a list of arrays.
``experiment.raw[roi, trial]`` gives you the traces of a specific ROI
and trial, across the ROI and neuropil regions.
``experiment.roi_polys`` is a list of lists of arrays.
``experiment.roi_polys[roi, trial][region][0]`` gives you the
polygon for the region for a specific ROI, trial and region.
``region=0`` is the ROI itself (i.e. the outline of the neuron cell),
and ``region>0`` gives the different neuropil regions.
For separable masks, it is possible multiple outlines are
found, which can be accessed as
``experiment.roi_polys[roi, trial][region][i]``,
where ``i`` is the outline index.
Parameters
----------
redo : bool, optional
If ``False``, we load previously prepared data when possible.
If ``True``, we re-run the preparation, even if it has previously
been run. Default is ``False``.
"""
# Get the timestamp for program start
t0 = time.time()
# define filename where data will be present
if self.folder is None:
fname = None
redo = True
else:
fname = os.path.join(self.folder, "prepared.npz")
# try to load data from filename
if fname is None or not os.path.isfile(fname):
redo = True
if not redo:
try:
self.clear()
self.load(fname)
if self.raw is not None:
return
except BaseException as err:
print("An error occurred while loading {}".format(fname))
print(err)
print("Extraction will be redone and {} overwritten".format(fname))
# Wipe outputs
self.clear()
# Adopt default values
self._adopt_default_parameters(only_preparation=True)
# Extract signals
n_trial = len(self.images)
if self.verbosity >= 2:
msg = "Doing region growing and data extraction for {} trials...".format(
n_trial
)
msg += "\n Images:"
for image in self.images:
if self.verbosity >= 4 or isinstance(image, basestring):
msg += "\n {}".format(image)
else:
msg += "\n {}".format(image.__class__)
msg += "\n ROI sets:"
for roiset in self.rois:
if self.verbosity >= 4 or isinstance(roiset, basestring):
msg += "\n {}".format(roiset)
else:
msg += "\n {}".format(roiset.__class__)
for key in self._preparation_params:
msg += "\n {}: {}".format(key, repr(getattr(self, key)))
print(msg)
sys.stdout.flush()
# Make a handle to the extraction function with parameters configured
_extract_cfg = functools.partial(
extract,
nRegions=self.nRegions,
expansion=self.expansion,
datahandler=self.datahandler,
verbosity=self.verbosity - 2,
total=n_trial,
)
# check whether we should show progress bars
disable_progressbars = self.verbosity <= 0 or 3 <= self.verbosity
# Check how many workers to spawn.
# Map the behaviour of ncores=None to one job per CPU core, like for
# multiprocessing.Pool(processes=None). With joblib, this is
# joblib.Parallel(n_jobs=-1) instead.
n_jobs = -1 if self.ncores_preparation is None else self.ncores_preparation
if 0 <= n_jobs <= 1:
# Don't use multiprocessing
outputs = [
_extract_cfg(image, rois, label=i)
for i, (image, rois) in tqdm(
enumerate(zip(self.images, self.rois)),
total=n_trial,
desc="Extracting traces",
disable=disable_progressbars,
)
]
else:
# Use multiprocessing
outputs = Parallel(
n_jobs=n_jobs, backend="threading", verbose=max(0, self.verbosity - 5)
)(
delayed(_extract_cfg)(image, rois, label=i)
for i, (image, rois) in tqdm(
enumerate(zip(self.images, self.rois)),
total=n_trial,
desc="Extracting traces",
disable=disable_progressbars,
)
)
# get number of cells
n_roi = len(outputs[0][1])
# predefine data structures
raw = np.empty((n_roi, n_trial), dtype=object)
roi_polys = np.empty_like(raw)
# Set outputs
means = []
for i_trial, (raw_i, polys_i, mean_i) in enumerate(outputs):
means.append(mean_i)
for i_roi in range(n_roi):
raw[i_roi][i_trial] = raw_i[i_roi]
roi_polys[i_roi][i_trial] = polys_i[i_roi]
self.raw = raw
self.roi_polys = roi_polys
self.means = means
if self.verbosity >= 1:
print(
"Finished extracting raw signals from {} ROIs across {} trials in {}.".format(
n_roi,
n_trial,
_pretty_timedelta(seconds=time.time() - t0),
)
)
sys.stdout.flush()
# Maybe save to cache file
if self.folder is not None:
self.save_prep()
[docs] def save_prep(self, destination=None):
r"""
Save prepared raw signals, extracted from images, to an npz file.
.. versionadded:: 1.0.0
Parameters
----------
destination : str, optional
Path to output file. The default destination is
``"prepared.npz"`` within the cache directory
``experiment.folder``.
"""
fields = set(self._preparation_params + self._preparation_outputs)
if destination is None:
if self.folder is None:
raise ValueError(
"The folder attribute must be declared in order to save"
" preparation outputs the cache."
)
destination = os.path.join(self.folder, "prepared.npz")
if self.verbosity >= 1:
print("Saving extracted traces to {}".format(destination))
sys.stdout.flush()
destdir = os.path.dirname(destination)
if destdir and not os.path.isdir(destdir):
os.makedirs(destdir)
np.savez_compressed(
destination,
**{
field: getattr(self, field)
for field in fields
if getattr(self, field) is not None
}
)
[docs] def separate(self, redo_prep=False, redo_sep=False):
r"""
Separate all the trials with FISSA algorithm.
After running ``separate``, data can be found as follows:
experiment.sep
Raw separation output, without being matched. Signal ``i`` for
a specific ROI and trial can be found as
``experiment.sep[roi, trial][i, :]``.
experiment.result
Final output, in order of presence in the ROI.
Signal ``i`` for a specific ROI and trial can be found at
``experiment.result[roi, trial][i, :]``.
Note that the ordering is such that ``i = 0`` is the signal
most strongly present in the ROI, and subsequent entries
are in diminishing order.
experiment.mixmat
The mixing matrix, which maps from ``experiment.raw`` to
``experiment.sep``.
experiment.info
Information about separation routine, iterations needed, etc.
Parameters
----------
redo_prep : bool, optional
Whether to redo the preparation. Default is ``False.`` Note that
if this is true, we set ``redo_sep = True`` as well.
redo_sep : bool, optional
Whether to redo the separation. Default is ``False``. Note that
this parameter is ignored if `redo_prep` is set to ``True``.
"""
# Get the timestamp for program start
t0 = time.time()
# Do data preparation
if redo_prep or self.raw is None:
self.separation_prep(redo_prep)
if redo_prep:
redo_sep = True
# Define filename to store data in
if self.folder is None:
fname = None
redo_sep = True
else:
fname = os.path.join(self.folder, "separated.npz")
if fname is None or not os.path.isfile(fname):
redo_sep = True
if not redo_sep:
try:
self.clear_separated()
self.load(fname)
if self.result is not None:
return
except BaseException as err:
print("An error occurred while loading {}".format(fname))
print(err)
print(
"Signal separation will be redone and {} overwritten"
"".format(fname)
)
# Wipe outputs
self.clear_separated()
# Adopt default values
self._adopt_default_parameters()
# Check size of the input arrays
n_roi = len(self.raw)
n_trial = len(self.raw[0])
# Print what data will be analysed
if self.verbosity >= 2:
msg = "Doing signal separation for {} ROIs over {} trials...".format(
n_roi, n_trial
)
msg += "\n method: {}".format(repr(self.method))
if "ica" not in self.method.lower():
msg += "\n alpha: {}".format(repr(self.alpha))
for key in ["max_iter", "max_tries", "tol"]:
msg += "\n {}: {}".format(key, repr(getattr(self, key)))
print(msg)
sys.stdout.flush()
# Make a handle to the separation function with parameters configured
_separate_cfg = functools.partial(
separate_trials,
alpha=self.alpha,
max_iter=self.max_iter,
tol=self.tol,
max_tries=self.max_tries,
method=self.method,
verbosity=self.verbosity - 2,
total=n_roi,
)
# check whether we should show progress bars
disable_progressbars = self.verbosity <= 0 or 3 <= self.verbosity
# Check how many workers to spawn.
# Map the behaviour of ncores=None to one job per CPU core, like for
# multiprocessing.Pool(processes=None). With joblib, this is
# joblib.Parallel(n_jobs=-1) instead.
n_jobs = -1 if self.ncores_separation is None else self.ncores_separation
# Do the extraction
if 0 <= n_jobs <= 1:
# Don't use multiprocessing
outputs = [
_separate_cfg(X, label=i)
for i, X in tqdm(
enumerate(self.raw),
total=len(self.raw),
desc="Separating data",
disable=disable_progressbars,
)
]
else:
# Use multiprocessing
outputs = Parallel(
n_jobs=n_jobs, backend="threading", verbose=max(0, self.verbosity - 5)
)(
delayed(_separate_cfg)(X, label=i)
for i, X in tqdm(
enumerate(self.raw),
total=len(self.raw),
desc="Separating data",
disable=disable_progressbars,
)
)
# Define output shape as an array of objects shaped (n_roi, n_trial)
sep = np.empty((n_roi, n_trial), dtype=object)
result = np.empty_like(sep)
mixmat = np.empty_like(sep)
info = np.empty_like(sep)
# Place our outputs into the initialised arrays
for i_roi, (sep_i, match_i, mixmat_i, conv_i) in enumerate(outputs):
sep[i_roi, :] = sep_i
result[i_roi, :] = match_i
mixmat[i_roi, :] = [mixmat_i] * n_trial
info[i_roi, :] = conv_i
# list non-converged cells
non_converged_rois = [
i_roi for i_roi, info_i in enumerate(info) if not info_i[0]["converged"]
]
if self.verbosity >= 1:
message = "Finished separating signals from {} ROIs across {} trials in {}".format(
n_roi,
n_trial,
_pretty_timedelta(seconds=time.time() - t0),
)
if len(non_converged_rois) > 0:
message += (
"\n"
"Separation did not converge for the following {} ROIs: {}."
"\nConsider increasing max_iter (currently set to {})"
" or other FISSA parameters if this happens often and/or"
" to a lot of cells.".format(
len(non_converged_rois), non_converged_rois, self.max_iter
)
)
print(message)
sys.stdout.flush()
# Set outputs
self.info = info
self.mixmat = mixmat
self.sep = sep
self.result = result
# Maybe save to cache file
if self.folder is not None:
self.save_separated()
[docs] def save_separated(self, destination=None):
r"""
Save separated signals to an npz file.
.. versionadded:: 1.0.0
Parameters
----------
destination : str, optional
Path to output file. The default destination is ``"separated.npz"``
within the cache directory ``experiment.folder``.
"""
fields = set(
self._preparation_params
+ self._separation_params
+ self._separation_outputs
+ self._deltaf_outputs
)
if destination is None:
if self.folder is None:
raise ValueError(
"The folder attribute must be declared in order to save"
" separation outputs to the cache."
)
destination = os.path.join(self.folder, "separated.npz")
if self.verbosity >= 1:
print("Saving results to {}".format(destination))
sys.stdout.flush()
destdir = os.path.dirname(destination)
if destdir and not os.path.isdir(destdir):
os.makedirs(destdir)
np.savez_compressed(
destination,
**{
field: getattr(self, field)
for field in fields
if getattr(self, field) is not None
}
)
[docs] def calc_deltaf(self, freq, use_raw_f0=True, across_trials=True):
r"""
Calculate Δf/f0 for raw and result traces.
The outputs are found in the :attr:`deltaf_raw` and
:attr:`deltaf_result` attributes, which can be accessed at
``experiment.deltaf_raw`` and ``experiment.deltaf_result``.
Parameters
----------
freq : float
Imaging frequency, in Hz.
use_raw_f0 : bool, optional
If ``True`` (default), use an f0 estimate from the raw ROI trace
for both raw and result traces. If ``False``, use individual f0
estimates for each of the traces.
across_trials : bool, optional
If ``True``, we estimate a single baseline f0 value across all
trials. If ``False``, each trial will have their own baseline f0,
and Δf/f\ :sub:`0` value will be relative to the trial-specific f0.
Default is ``True``.
"""
# Get the timestamp for program start
t0 = time.time()
if self.verbosity >= 2:
msg = "Calculating Δf/f0 for raw and result signals"
if self.verbosity < 3:
pass
elif across_trials:
msg += " (same f0 across all trials"
else:
msg += " (different f0 baseline for each trial"
if self.verbosity < 3:
pass
elif use_raw_f0:
msg += ", using f0 in raw data for result)"
else:
msg += ")"
msg += (
"\nCaution: Measuring baseline f0 from result may result"
" in division by zero."
)
print(msg)
sys.stdout.flush()
# Initialise output arrays
deltaf_raw = np.empty_like(self.raw)
deltaf_result = np.empty_like(self.result)
# Can't include Δ in the tqdm description on Python2
desc = "Calculating {}f/f0".format("d" if sys.version_info < (3, 0) else "Δ")
# Check size of the input arrays
n_roi = len(self.result)
n_trial = len(self.result[0])
# Loop over cells
for i_roi in tqdm(
range(n_roi),
total=n_roi,
desc=desc,
disable=self.verbosity < 1,
):
# if deltaf should be calculated across all trials
if across_trials:
# get concatenated traces
raw_conc = np.concatenate(self.raw[i_roi], axis=1)[0, :]
result_conc = np.concatenate(self.result[i_roi], axis=1)
# calculate Δf/f0
raw_f0 = deltaf.findBaselineF0(raw_conc, freq)
raw_conc = (raw_conc - raw_f0) / raw_f0
result_f0 = deltaf.findBaselineF0(result_conc, freq, 1).T[:, None]
if use_raw_f0:
result_conc = (result_conc - result_f0) / raw_f0
else:
result_conc = (result_conc - result_f0) / result_f0
# store Δf/f0
curTrial = 0
for i_trial in range(n_trial):
nextTrial = curTrial + self.raw[i_roi][i_trial].shape[1]
signal = raw_conc[curTrial:nextTrial]
deltaf_raw[i_roi][i_trial] = np.expand_dims(signal, axis=0)
signal = result_conc[:, curTrial:nextTrial]
deltaf_result[i_roi][i_trial] = signal
curTrial = nextTrial
else:
# loop across trials
for i_trial in range(n_trial):
# get current signals
raw_sig = self.raw[i_roi][i_trial][0, :]
result_sig = self.result[i_roi][i_trial]
# calculate Δf/fo
raw_f0 = deltaf.findBaselineF0(raw_sig, freq)
result_f0 = deltaf.findBaselineF0(result_sig, freq, 1).T[:, None]
result_f0[result_f0 < 0] = 0
raw_sig = (raw_sig - raw_f0) / raw_f0
if use_raw_f0:
result_sig = (result_sig - result_f0) / raw_f0
else:
result_sig = (result_sig - result_f0) / result_f0
# store Δf/f0
deltaf_raw[i_roi][i_trial] = np.expand_dims(raw_sig, axis=0)
deltaf_result[i_roi][i_trial] = result_sig
self.deltaf_raw = deltaf_raw
self.deltaf_result = deltaf_result
if self.verbosity >= 1:
print(
"Finished calculating Δf/f0 for raw and result signals in {}".format(
_pretty_timedelta(seconds=time.time() - t0)
)
)
sys.stdout.flush()
# Maybe save to cache file
if self.folder is not None:
self.save_separated()
[docs] def to_matfile(self, fname=None, legacy=False):
r"""Save the results to a MATLAB file.
.. versionadded:: 1.0.0
This will generate a MAT-file (.mat) which can be loaded into MATLAB.
The MAT-file contains structs for all the experiment output attributes
(:attr:`roi_polys`, :attr:`result`, :attr:`raw`, etc.)
and analysis parameters (:attr:`expansion`, :attr:`nRegions`,
:attr:`alpha`, etc.).
If Δf/f\ :sub:`0` was calculated with :meth:`calc_deltaf`,
:attr:`deltaf_result` and :attr:`deltaf_raw` are also included.
These can be interfaced with as illustrated below.
``result{1, 1}(1, :)``
The separated signal for the first ROI and first trial.
This is equivalent to ``experiment.result[0, 0][0, :]`` when
interacting with the :class:`Experiment` object in Python.
``result{roi, trial}(1, :)``
The separated signal for the ``roi``-th ROI and ``trial``-th trial.
This is equivalent to
``experiment.result[roi - 1, trial - 1][0, :]`` when
interacting with the :class:`Experiment` object in Python.
``result{roi, trial}(2, :)``
A contaminating signal.
``raw{roi, trial}(1, :)``
Raw measured neuronal signal, averaged over the ROI.
This is equivalent to ``experiment.raw[roi - 1, trial - 1][0, :]``
when interacting with the :class:`Experiment` object in Python.
``raw{roi, trial}(2, :)``
Raw signal from first neuropil region (of ``nRegions``).
``roi_polys{roi, trial}{1}``
Polygon outlining the ROI, as an n-by-2 array of coordinates.
``roi_polys{roi, trial}{2}``
Polygon outlining the first neuropil region (of ``nRegions``),
as an n-by-2 array of coordinates.
Examples
--------
Here are some example MATLAB plots.
Plotting raw and decontaminated traces:
.. code:: octave
% Load the FISSA output data
S = load('separated.mat')
% Separated signal for the third ROI, second trial
roi = 3; trial = 2;
% Plot the raw and result traces for the ROI signal
figure; hold on;
plot(S.raw{roi, trial}(1, :));
plot(S.result{roi, trial}(1, :));
title(sprintf('ROI %d, Trial %d', roi, trial));
xlabel('Time (frame number)');
ylabel('Signal intensity (candela per unit area)');
legend({'Raw', 'Result'});
If all ROIs are contiguous and described by a single contour,
the the mean image and ROI locations for one trial can be plotted as
follows:
.. code:: octave
% Load the FISSA output data
S = load('separated.mat')
trial = 1;
figure; hold on;
% Plot the mean image
imagesc(squeeze(S.means(trial, :, :)));
colormap('gray');
% Plot ROI locations
for i_roi = 1:size(S.result, 1);
contour = S.roi_polys{i_roi, trial}{1};
plot(contour(:, 2), contour(:, 1));
end
set(gca, 'YDir', 'reverse');
Parameters
----------
fname : str, optional
Destination for output file. The default is a file named
``"separated.mat"`` within the cache save directory for the
experiment (the :attr:`folder` argument when the
:class:`Experiment` instance was created).
legacy : bool, default=False
Whether to use the legacy format of :meth:`save_to_matlab`.
This also changes the default output name to ``"matlab.mat"``.
"""
default_name = "separated.mat"
if legacy:
default_name = "matlab.mat"
# define filename
if fname is None:
if self.folder is None:
raise ValueError(
"fname must be provided if experiment folder is undefined"
)
fname = os.path.join(self.folder, default_name)
# initialize dictionary to save
M = collections.OrderedDict()
def reformat_dict_for_legacy(orig_dict):
new_dict = collections.OrderedDict()
# loop over cells and trial
for i_roi in range(len(self.result)):
# get current cell label
c_lab = "cell" + str(i_roi)
# update dictionary
new_dict[c_lab] = collections.OrderedDict()
for i_trial in range(len(self.result[0])):
# get current trial label
t_lab = "trial" + str(i_trial)
# update dictionary
new_dict[c_lab][t_lab] = orig_dict[i_roi][i_trial]
return new_dict
if legacy:
M["ROIs"] = reformat_dict_for_legacy(self.roi_polys)
M["raw"] = reformat_dict_for_legacy(self.raw)
M["result"] = reformat_dict_for_legacy(self.result)
if getattr(self, "deltaf_raw", None) is not None:
M["df_raw"] = reformat_dict_for_legacy(self.deltaf_raw)
if getattr(self, "deltaf_result", None) is not None:
M["df_result"] = reformat_dict_for_legacy(self.deltaf_result)
else:
fields = (
self._preparation_params
+ self._separation_params
+ self._preparation_outputs
+ self._separation_outputs
+ self._deltaf_outputs
)
for field in fields:
x = getattr(self, field)
if x is None:
continue
M[field] = x
with warnings.catch_warnings():
warnings.filterwarnings(
"ignore",
message="Creating an ndarray from ragged nested sequences",
)
savemat(fname, M)
[docs] def save_to_matlab(self, fname=None):
r"""
Save the results to a MATLAB file.
.. deprecated:: 1.0.0
Use ``experiment.to_matfile(legacy=True)`` instead.
This will generate a .mat file which can be loaded into MATLAB to
provide structs: ROIs, result, raw.
If Δf/f\ :sub:`0` was calculated, these will also be stored as ``df_result``
and ``df_raw``, which will have the same format as ``result`` and
``raw``.
These can be interfaced with as follows, for ROI 0, trial 0:
``ROIs.cell0.trial0{1}``
Polygon outlining the ROI.
``ROIs.cell0.trial0{2}``
Polygon outlining the first (of ``nRegions``) neuropil region.
``result.cell0.trial0(1, :)``
Final extracted neuronal signal.
``result.cell0.trial0(2, :)``
Contaminating signal.
``raw.cell0.trial0(1, :)``
Raw measured cell signal, average over the ROI.
``raw.cell0.trial0(2, :)``
Raw signal from first (of ``nRegions``) neuropil region.
Parameters
----------
fname : str, optional
Destination for output file. Default is a file named
``"matlab.mat"`` within the cache save directory for the experiment
(the `folder` argument when the ``Experiment`` instance was created).
See Also
--------
Experiment.to_matfile
"""
warnings.warn(
"The experiment.save_to_matlab() method is deprecated."
" Please use experiment.to_matfile() instead.",
DeprecationWarning,
)
return self.to_matfile(fname=fname, legacy=True)
[docs]def run_fissa(
images,
rois,
folder=None,
freq=None,
return_deltaf=False,
deltaf_across_trials=True,
export_to_matfile=False,
**kwargs
):
r"""
Function-based interface to run FISSA on an experiment.
.. versionadded:: 1.0.0
Uses the methodology described in
`FISSA: A neuropil decontamination toolbox for calcium imaging signals <doi_>`_.
.. _doi: https://www.doi.org/10.1038/s41598-018-21640-2
Parameters
----------
images : str or list
The raw recording data.
Should be one of:
- the path to a directory containing TIFF files (string),
- a list of paths to TIFF files (list of strings),
- a list of :term:`array_like` data already loaded into memory, each
shaped ``(n_frames, height, width)``.
Note that each TIFF/array is considered a single trial.
rois : str or list
The roi definitions.
Should be one of:
- the path to a directory containing ImageJ ZIP files (string),
- the path of a single ImageJ ZIP file (string),
- a list of ImageJ ZIP files (list of strings),
- a list of arrays, each encoding a ROI polygons,
- a list of lists of binary arrays, each representing a ROI mask.
This can either be a single roiset for all trials, or a different
roiset for each trial.
folder : str, optional
Path to a cache directory from which pre-extracted data will
be loaded if present, and saved to otherwise. If `folder` is
unset, the experiment data will not be saved.
freq : float, optional
Imaging frequency, in Hz. Required if ``return_deltaf=True``.
return_deltaf : bool, optional
Whether to return Δf/f\ :sub:`0`. Otherwise, the decontaminated signal
is returned scaled against the raw recording. Default is ``False``.
deltaf_across_trials : bool, default=True
If ``True``, we estimate a single baseline f0 value across all
trials when computing Δf/f\ :sub:`0`.
If ``False``, each trial will have their own baseline f0, and
Δf/f\ :sub:`0` value will be relative to the trial-specific f0.
Default is ``True``.
export_to_matfile : bool or str or None, default=False
Whether to export the data to a MATLAB-compatible .mat file.
If `export_to_matfile` is a string, it is used as the path to the output
file. If ``export_to_matfile=True``, the matfile is saved to the
default path of ``"separated.mat"`` within the `folder` directory, and
`folder` must be set. If this is ``None``, the matfile is exported to
the default path if `folder` is set, and otherwise is not exported.
Default is ``False``.
**kwargs
Additional keyword arguments as per :class:`Experiment`.
Returns
-------
result : 2d numpy.ndarray of 2d numpy.ndarrays of np.float64
The vector ``result[roi, trial][0, :]`` is the trace from ROI ``roi``
in trial ``trial``.
If ``return_deltaf=True``, this is Δf/f\ :sub:`0`;
otherwise, it is the decontaminated signal scaled as per the raw
signal. f\ :sub:`0` is the baseline as calculated from the raw
signal.
raw : 2d numpy.ndarray of 2d numpy.ndarrays of np.float64
The raw traces without separation.
The vector ``raw[c, t][0, :]`` is the ROI trace from cell ``c`` in
trial ``t``. The vector ``raw[c, t][i, :]`` for i>=1 the trace from
cell ``c`` in trial ``t``, from neuropil region ``i-1``.
If ``return_deltaf=True``, this is Δf/f\ :sub:`0`; otherwise it's
the raw extracted signal.
See Also
--------
fissa.core.Experiment
"""
# Parse arguments
if export_to_matfile is None:
export_to_matfile = folder is not None
if return_deltaf and freq is None:
raise ValueError("The argument `freq` must be set to determine df/f0.")
# Make a new Experiment object
experiment = Experiment(images, rois, folder=folder, **kwargs)
# Run separation
experiment.separate()
# Calculate Δf/f0
if return_deltaf or (export_to_matfile and freq is not None):
experiment.calc_deltaf(freq=freq, across_trials=deltaf_across_trials)
# Save to matfile
if export_to_matfile:
matlab_fname = (
None if isinstance(export_to_matfile, bool) else export_to_matfile
)
experiment.to_matfile(matlab_fname)
# Return appropriate data
if return_deltaf:
return experiment.deltaf_result, experiment.deltaf_raw
return experiment.result, experiment.raw