Core Experiment Classes

This module contains the core data structures for managing neural experiments.

Classes

class driada.experiment.exp_base.Experiment(signature, calcium, spikes, exp_identificators, static_features, dynamic_features, ground_truth=None, **kwargs)

Bases: object

Base class for calcium imaging and spike train experiments.

This class provides a unified interface for analyzing neural activity data (calcium imaging or spike trains) in relation to various experimental features (behavioral variables, stimuli, etc.). It handles data organization, feature extraction, mutual information analysis, and statistical significance testing.

Parameters:

• signature (str) – Unique identifier for the experiment.

• calcium (numpy.ndarray) – Calcium imaging data of shape (n_neurons, n_timepoints). Required parameter.

• spikes (numpy.ndarray or None) – Spike train data. Can be provided directly or reconstructed from calcium.

• exp_identificators (dict) – Experiment metadata and identifiers. Each key-value pair becomes an attribute of the object (e.g., exp_identificators={‘mouse_id’: ‘M1’} creates self.mouse_id = ‘M1’).

• static_features (dict) – Time-invariant features (e.g., cell types, anatomical properties). Should include ‘fps’, ‘t_rise_sec’, ‘t_off_sec’ or defaults will be used. The dict is stored as self.static_features, and each key also becomes an individual attribute (e.g., static_features={‘fps’: 20} creates both self.static_features[‘fps’] and self.fps = 20).

• dynamic_features (dict) – Time-varying features (e.g., behavior, stimuli). Keys are feature names, values are TimeSeries, MultiTimeSeries, or numpy arrays. All features must have the same length (number of timepoints).

• **kwargs (dict) –

  Additional parameters including:

  • optimize_kinetics : bool or str, optimize kinetics per neuron (default: False). If True, uses ‘lbfgs’ method. Can specify ‘lbfgs’ or ‘grid’ explicitly.

  • reconstruct_spikes : str or bool, spike reconstruction method (default: ‘wavelet’)

  • bad_frames_mask : array-like, boolean mask where True indicates bad frames

  • spike_kwargs : dict, parameters for spike reconstruction

  • verbose : bool, print progress messages (default: True)

signature

Experiment identifier.

Type:

str

neurons

List of Neuron objects, indexed by cell ID (0-based). Access with self.neurons[cell_id].

Type:

list

n_cells

Number of neurons in the experiment.

Type:

int

n_frames

Number of time points (frames) in the experiment.

Type:

int

calcium

Calcium imaging data as MultiTimeSeries object.

Type:

MultiTimeSeries

spikes

Spike train data as MultiTimeSeries object.

Type:

MultiTimeSeries

static_features

Time-invariant experimental features as originally provided.

Type:

dict

dynamic_features

Time-varying experimental features as TimeSeries/MultiTimeSeries objects.

Type:

dict

stats_tables

Nested dict storing mutual information statistics. Structure: stats_tables[mode][feat_id][cell_id] = stats_dict.

Type:

dict

significance_tables

Nested dict storing statistical significance data. Structure: significance_tables[mode][feat_id][cell_id] = sig_dict.

Type:

dict

embeddings

Stored dimensionality reduction results by data type and method. Structure: embeddings[data_type][method_name] = embedding_data.

Type:

dict

verbose

Whether to print progress messages.

Type:

bool

spike_reconstruction_method

Method used for spike reconstruction if applicable.

Type:

str or None

filtered_flag

Whether bad frames were filtered out.

Type:

bool

selectivity_tables_initialized

Whether selectivity tables have been initialized.

Type:

bool

exp_identificators

Original experiment identifiers dictionary.

Type:

dict

_data_hashes

Private attribute storing hash representations for caching.

Type:

dict

_rdm_cache

Private attribute caching representational dissimilarity matrices.

Type:

dict

check_ds(ds)

Validate downsampling rate for behavioral analysis.

get_neuron_feature_pair_stats(cell_id, feat_id, mode='calcium')

Get selectivity statistics for a neuron-feature pair.

get_neuron_feature_pair_significance(cell_id, feat_id, mode='calcium')

Get statistical significance data for a neuron-feature pair.

update_neuron_feature_pair_stats(stats, cell_id, feat_id, mode='calcium', ...)

Update statistics for a neuron-feature pair.

update_neuron_feature_pair_significance(sig, cell_id, feat_id, mode='calcium')

Update significance data for a neuron-feature pair.

get_multicell_shuffled_calcium(cbunch=None, method='roll_based', \*\*kwargs)

Get shuffled calcium data for specified neurons.

get_multicell_shuffled_spikes(cbunch=None, method='isi_based', \*\*kwargs)

Get shuffled spike data for specified neurons.

get_stats_slice(cell_ids, feat_ids, mode='calcium', vars=None)

Extract statistics for multiple neuron-feature pairs.

get_significance_slice(cell_ids, feat_ids, mode='calcium', vars=None)

Extract significance data for multiple neuron-feature pairs.

get_feature_entropy(feat_id, ds=1)

Calculate Shannon entropy of a feature.

get_significant_neurons(min_nspec=1, cbunch=None, fbunch=None, mode='calcium', ...)

Find neurons with significant selectivity to one or more features.

store_embedding(embedding, method_name, data_type='calcium', metadata=None)

Store dimensionality reduction results.

get_embedding(method_name, data_type='calcium')

Retrieve stored dimensionality reduction embedding.

compute_rdm(items, activity_type='calcium', metric='correlation', \*\*kwargs)

Compute representational dissimilarity matrix.

clear_rdm_cache()

Clear the RDM computation cache.

Notes

The class supports both calcium imaging and spike train analysis through the ‘mode’ parameter in various methods. Results are cached using hash-based lookups to avoid redundant computations. Statistical significance is determined using the INTENSE algorithm with two-stage hypothesis testing.

Spike reconstruction is performed automatically if spikes are not provided and reconstruct_spikes is not False or None. The ‘wavelet’ method is recommended for calcium imaging data. Set reconstruct_spikes to False or None to disable spike reconstruction entirely.

Individual static and dynamic features can be accessed as attributes. For example, if ‘position’ is a dynamic feature, access it via self.position. Protected attribute names will have an underscore prefix if conflicts occur.

__init__(signature, calcium, spikes, exp_identificators, static_features, dynamic_features, ground_truth=None, **kwargs)

Initialize experiment with neural data and behavioral features.

Creates an experiment object that integrates calcium imaging data, spike trains, and behavioral features for neural population analysis. Handles data validation, spike reconstruction, and sets up internal data structures for statistical analysis.

Parameters:

• signature (str) – Unique identifier for the experiment (e.g., ‘mouse123_session1’).

• calcium (array-like) – Calcium imaging data with shape (n_cells, n_frames). Required. Each row is a neuron’s calcium trace over time.

• spikes (array-like or None) – Spike train data with same shape as calcium. If None and reconstruct_spikes is specified in kwargs, spikes will be reconstructed from calcium data.

• exp_identificators (dict or None) – Metadata about the experiment (e.g., subject_id, session_date). Keys become attributes of the experiment object.

• static_features (dict or None) – Time-invariant features. Expected keys include: - ‘fps’: sampling rate (frames per second) - ‘t_rise_sec’: calcium rise time in seconds - ‘t_off_sec’: calcium decay time in seconds - Other experiment-wide parameters

• dynamic_features (dict or None) – Time-varying behavioral features (e.g., position, speed). Values should be array-like with time dimension matching n_frames. Keys become accessible via self.dynamic_features.

• ground_truth (dict or None, optional) – Ground truth for synthetic experiments. None for real data. For synthetic data, typically contains: - “expected_pairs”: list of (neuron_idx, feature_name) tuples - “tuning_parameters”: dict of per-neuron tuning params - “neuron_types”: dict mapping neuron_idx to group name - “population_config”: original population configuration

• **kwargs –

  Additional parameters:

  • optimize_kinetics (bool or str): Optimize kinetics per neuron. Default False. If True, uses ‘lbfgs’ method. Can specify ‘lbfgs’ or ‘grid’ explicitly.

  • reconstruct_spikes (str, False, or None): Method for spike reconstruction. Options: ‘wavelet’ (default), ‘threshold’, False, or None. If False or None, spike reconstruction is disabled. Only used if spikes is None.

  • bad_frames_mask (array-like): Boolean mask of frames to exclude.

  • spike_kwargs (dict): Parameters for spike reconstruction method.

  • verbose (bool): Print progress messages. Default True.

  • create_circular_2d (bool): If True, automatically create _2d versions of circular features (detected via type_info.is_circular) as (cos, sin) MultiTimeSeries. Original features are preserved. E.g., ‘headdirection’ -> also creates ‘headdirection_2d’. Default True.

Raises:

• ValueError – If calcium is None, if data shapes are inconsistent, if feature names conflict with protected attributes, or if data appears transposed (n_cells > n_frames).

• TypeError – If dynamic features have incompatible types.

Warning

UserWarning

If both spikes and reconstruct_spikes are provided (spikes will be overwritten), or if static feature names conflict with existing attributes (will be prefixed with underscore).

Notes

Protected attribute names that cannot be used as feature names: ‘spikes’, ‘calcium’, ‘neurons’, ‘n_cells’, ‘n_frames’, ‘static_features’, ‘dynamic_features’, ‘downsampling’, ‘significance_tables’, ‘stats_tables’, ‘_data_hashes’, ‘embeddings’, ‘_rdm_cache’, ‘intense_results’

The initialization process: 1. Validates and stores basic data (calcium required) 2. Handles spike data or reconstruction 3. Creates Neuron objects for each cell 4. Processes static and dynamic features 5. Builds internal data structures for caching 6. Validates data consistency

Examples

>>> import numpy as np
>>> from driada.information.info_base import TimeSeries
>>>
>>> # Basic initialization with calcium data only
>>> calcium_data = np.random.randn(10, 1000)  # 10 neurons, 1000 timepoints
>>> exp = Experiment('exp001', calcium_data, None, {}, {}, {},
...                  reconstruct_spikes=None, verbose=False)

>>> # With spikes and behavioral features
>>> spike_data = np.random.poisson(0.05, (10, 1000))  # spike trains
>>> speed_trace = TimeSeries(np.random.rand(1000))    # behavioral data
>>> exp = Experiment(
...     'exp002',
...     calcium_data,
...     spike_data,
...     {},
...     {'fps': 30.0},
...     {'speed': speed_trace},
...     verbose=False
... )

>>> # Without spike reconstruction (faster for doctests)
>>> exp = Experiment(
...     'exp003',
...     calcium_data,
...     None,
...     {},
...     {'fps': 30.0},
...     {},
...     reconstruct_spikes=None,
...     verbose=False
... )

check_ds(ds)

Check if downsampling rate is appropriate for behavior analysis.

Validates that the downsampling rate won’t cause time gaps larger than the minimum behavior time interval (0.25 seconds), which could lead to missed behavioral events.

Parameters:

ds (int) – Downsampling factor. The data will be sampled every ds frames. Must be a positive integer (>= 1).

Return type:

None

Raises:

• ValueError – If fps (frames per second) is not set for this experiment, or if ds is less than 1. Error messages include relevant context.

• TypeError – If ds is not an integer.

Warning

UserWarning

Issued if the time gap created by downsampling exceeds DEFAULT_MIN_BEHAVIOUR_TIME (0.25 seconds). The warning includes the current threshold, downsampling factor, and resulting time gap.

Notes

The time gap is calculated as (1/fps) * ds seconds. For example, with fps=20 and ds=10, the time gap would be 0.5 seconds, which exceeds the 0.25 second threshold and triggers a warning.

add_feature(name, data, ts_type=None, **ts_kwargs)

Add a dynamic feature to the experiment after initialization.

Mirrors the feature registration performed during __init__: wraps raw arrays, validates length, applies shuffle masks, stores in dynamic_features dict, and sets a direct attribute.

Parameters:

• name (str) – Feature name. Must not conflict with protected attributes.

• data (array_like or TimeSeries or MultiTimeSeries) – Feature data. Length must match n_frames. Raw arrays are wrapped in TimeSeries automatically.

• ts_type (str or None, optional) – Passed to TimeSeries(ts_type=...) when wrapping raw arrays. Common values: 'discrete', 'circular'. Ignored if data is already a TimeSeries/MultiTimeSeries.

• **ts_kwargs – Additional keyword arguments forwarded to the TimeSeries constructor (e.g. circular_period).

get_circular_2d_feature(name)

Get the _2d (cos, sin) version of a circular feature.

Parameters:

name (str) – Original circular feature name (without _2d suffix).

Returns:

The _2d version, or None if not available.

Return type:

MultiTimeSeries or None

Examples

>>> # Assuming exp has circular feature 'headdirection'
>>> mts = exp.get_circular_2d_feature('headdirection')  
>>> mts.n_dim  
2

has_circular_2d(name)

Check if a feature has a _2d (cos, sin) counterpart.

Parameters:

name (str) – Feature name to check.

Returns:

True if the _2d version exists.

Return type:

bool

get_circular_features()

Get all circular features (excluding _2d versions).

Returns:

Dictionary mapping feature names to TimeSeries objects for all circular features.

Return type:

dict

update_neuron_feature_pair_stats(stats, cell_id, feat_id, mode='calcium', force_update=False, stage2_only=False)

Updates calcium-feature pair statistics.

feat_id should be a string or an iterable of strings (in case of joint MI calculation). This function allows multifeatures.

Parameters:

• stats (dict) – Statistics dictionary containing the updates

• cell_id (int) – Cell identifier

• feat_id (str or iterable of str) – Feature identifier(s). Can be a string or iterable of strings for joint MI calculation

• mode (str, optional) – Data processing mode. Default is “calcium”

• force_update (bool, optional) – If True, force update even if data hashes match. Default is False

• stage2_only (bool, optional) – If True, only update stage 2 statistics. Default is False

update_neuron_feature_pair_significance(sig, cell_id, feat_id, mode='calcium')

Updates calcium-feature pair significance data.

feat_id should be a string or an iterable of strings (in case of joint MI calculation). This function allows multifeatures.

Parameters:

• sig (dict) – Significance data to update

• cell_id (int) – Cell identifier

• feat_id (str or iterable of str) – Feature identifier(s). Can be a string or iterable of strings for joint MI calculation

• mode (str, optional) – Data processing mode. Default is “calcium”

get_neuron_feature_pair_stats(cell_id, feat_id, mode='calcium')

Get selectivity statistics for a neuron-feature pair.

Retrieves pre-computed statistics measuring the relationship between neural activity and behavioral/experimental features. Supports both single features and multi-feature analysis.

Parameters:

• cell_id (int) – Neuron/cell identifier.

• feat_id (str or tuple of str) – Feature identifier(s). Can be a single feature name or tuple of feature names for joint analysis.

• mode ({'calcium', 'spikes'}, optional) – Type of neural activity. Default is ‘calcium’.

Returns:

Dictionary containing various statistical measures of the neuron-feature relationship. Returns None if statistics have not been computed or if data has changed since computation.

Return type:

dict or None

Notes

Statistics must be pre-computed using the selectivity analysis pipeline. The method checks data integrity using hashes to ensure statistics are up-to-date with the current data.

Note

Despite the ‘get’ name, this method can trigger side effects via _check_stats_relevance which may add new features to tables. This behavior will be removed after tuple multifeature deprecation.

get_neuron_feature_pair_significance(cell_id, feat_id, mode='calcium')

Get statistical significance data for a neuron-feature pair.

Retrieves significance testing results for the neuron-feature relationship, typically from shuffle-based permutation tests.

Parameters:

• cell_id (int) – Neuron/cell identifier.

• feat_id (str or tuple of str) – Feature identifier(s). Can be a single feature name or tuple of feature names for joint analysis.

• mode ({'calcium', 'spikes'}, optional) – Type of neural activity. Default is ‘calcium’.

Returns:

Dictionary containing significance test results including p-values, shuffle distributions, and multiple comparison corrections. Returns None if significance has not been computed or if underlying statistics are outdated.

Return type:

dict or None

Notes

Significance testing typically uses shuffle tests where temporal relationships are destroyed while preserving marginal distributions. Results include both single-stage and two-stage testing procedures.

Note

Despite the ‘get’ name, this method can trigger side effects via _check_stats_relevance which may add new features to tables. This behavior will be removed after tuple multifeature deprecation.

get_multicell_shuffled_calcium(cbunch=None, method='roll_based', return_array=True, **kwargs)

Get shuffled calcium data for multiple cells.

Parameters:

• cbunch (int, list, or None) – Cell indices. If None, all cells are used.

• method ({'roll_based', 'waveform_based', 'chunks_based'}, default='roll_based') – Shuffling method to use

• return_array (bool, default=True) – If True, return numpy array. If False, return MultiTimeSeries object.

• **kwargs – Additional parameters passed to the shuffling method

Returns:

If return_array=True: Shuffled calcium data with shape (n_cells, n_frames) If return_array=False: MultiTimeSeries object containing shuffled data

Return type:

np.ndarray or MultiTimeSeries

get_multicell_shuffled_spikes(cbunch=None, method='isi_based', return_array=True, **kwargs)

Get shuffled spike data for multiple cells.

Parameters:

• cbunch (int, list, or None) – Cell indices. If None, all cells are used.

• method ({'isi_based'}, default='isi_based') – Shuffling method. Currently only ‘isi_based’ is supported for spikes.

• return_array (bool, default=True) – If True, return numpy array. If False, return MultiTimeSeries object.

• **kwargs – Additional parameters passed to the shuffling method

Returns:

If return_array=True: Shuffled spike data with shape (n_cells, n_frames) If return_array=False: MultiTimeSeries object containing shuffled spike data

Return type:

np.ndarray or MultiTimeSeries

get_stats_slice(table_to_scan=None, cbunch=None, fbunch=None, sbunch=None, significance_mode=False, mode='calcium')

Returns slice of accumulated statistics data (or significance data if “significance_mode=True”).

Parameters:

• table_to_scan (dict, optional) – Specific table to scan. If None, uses default stats or significance table

• cbunch (int, list, or None, optional) – Cell identifiers to include. If None, includes all cells

• fbunch (str, list, or None, optional) – Feature identifiers to include. If None, includes all features

• sbunch (str, list, or None, optional) – Statistics keys to include. If None, includes all statistics

• significance_mode (bool, optional) – If True, returns significance data instead of statistics. Default is False

• mode (str, optional) – Data processing mode. Default is “calcium”

Returns:

Nested dictionary with structure table[feat_id][cell_id][stat_key]

Return type:

dict

get_significance_slice(cbunch=None, fbunch=None, sbunch=None, mode='calcium')

Extract significance test results for selected cells and features.

Convenience method that retrieves statistical significance data (p-values, test statistics, etc.) for specific cell-feature combinations. This is equivalent to calling get_stats_slice with significance_mode=True.

Parameters:

• cbunch (int, list of int, or None, optional) – Cell indices to include. None means all cells.

• fbunch (int, str, list, or None, optional) – Feature indices/names to include. None means all features.

• sbunch (str, list of str, or None, optional) – Significance measures to extract (e.g., ‘pval’, ‘qval’, ‘statistic’). None means all available measures.

• mode ({'calcium', 'spikes'}, default='calcium') – Which data type’s significance tables to use.

Returns:

Nested dictionary with structure: {feature: {cell: {measure: value}}}. Contains only the requested significance test results.

Return type:

dict

See also

get_stats_slice

More general method for extracting any statistics

Examples

>>> import numpy as np
>>> from driada.experiment.exp_base import Experiment
>>> from driada.information.info_base import TimeSeries
>>>
>>> # Create an example experiment with selectivity data
>>> calcium_data = np.random.randn(20, 1000)
>>> running_speed = TimeSeries(np.random.rand(1000))
>>> exp = Experiment(
...     'example',
...     calcium_data,
...     None,
...     {},
...     {'fps': 30.0},
...     {'running_speed': running_speed},
...     verbose=False
... )
>>>
>>> # Initialize stats tables
>>> exp._set_selectivity_tables('calcium')
>>>
>>> # Get p-values for cells 0-5 and feature 'running_speed'
>>> sig_data = exp.get_significance_slice(
...     cbunch=[0, 1, 2, 3, 4, 5],
...     fbunch=['running_speed'],
...     sbunch=['pval']
... )
>>> # Returns nested dict structure
>>> sorted(sig_data.keys())
['running_speed']
>>> sorted(sig_data['running_speed'].keys())
[0, 1, 2, 3, 4, 5]

get_feature_entropy(feat_id, ds=1)

Calculates entropy of a single dynamic feature or joint entropy of multiple features.

Parameters:

• feat_id (str or tuple) –

  • str: Name of a single feature (TimeSeries or MultiTimeSeries)

  • tuple: Names of exactly 2 features for joint entropy calculation

• ds (int) – Downsampling factor

Returns:

Entropy value in bits (or nats for continuous variables)

Return type:

float

Notes

• Single features use their native get_entropy() method

• Tuples calculate joint entropy for exactly 2 variables

• Joint entropy of 3+ variables is not supported (use MultiTimeSeries instead)

• Continuous variables may return negative entropy values

reconstruct_all_neurons(method='wavelet', n_iter=3, optimize_kinetics=True, hybrid_kinetics=True, wavelet=None, rel_wvt_times=None, show_progress=True, use_gpu=False, n_jobs=None, enable_parallelization=None, **kwargs)

Batch reconstruct spikes for all neurons with wavelet optimization.

Populates neuron.asp, neuron.sp, neuron.events for each neuron.

Parameters:

• method (str) – Reconstruction method (‘wavelet’ or ‘threshold’). Default ‘wavelet’.

• n_iter (int) – Number of iterations for iterative detection. Default 3.

• optimize_kinetics (bool) – Optimize calcium kinetics per neuron. Default True.

• hybrid_kinetics (bool) – Optimize kinetics BEFORE reconstruction (hybrid mode). Default True.

• wavelet (Wavelet, optional) – Pre-computed wavelet object. If None, creates once and reuses.

• rel_wvt_times (array-like, optional) – Pre-computed time resolutions. If None, computes once and reuses.

• show_progress (bool) – Show progress bar. Default True.

• use_gpu (bool) – Use GPU acceleration for wavelet transform. Default False. Requires PyTorch and CuPy. Ridge extraction remains CPU-only.

• n_jobs (int, optional) – Number of parallel jobs. Default uses experiment’s n_jobs setting. Use -1 for all available cores.

• enable_parallelization (bool, optional) – Enable parallel processing. Default uses experiment’s setting.

• **kwargs – Additional parameters passed to neuron.reconstruct_spikes()

Returns:

Updates neuron objects in-place and syncs to exp.spikes

Return type:

None

get_significant_neurons(min_nspec=1, cbunch=None, fbunch=None, mode='calcium', override_intense_significance=False, pval_thr=0.05, multicomp_correction=None, significance_update=False)

Returns a dict with neuron ids as keys and their significantly correlated features as values Only neurons with “min_nspec” or more significantly correlated features will be returned

Parameters:

• min_nspec (int) – Minimum number of significantly correlated features required

• cbunch (int, list or None) – Cell indices to analyze. By default (None), all neurons will be checked

• fbunch (str, list or None) – Feature names to check. By default (None), all features will be checked

• mode (str) – Data type: ‘calcium’ or ‘spikes’

• override_intense_significance (bool, optional) – If True, recompute significance using pval_thr and multicomp_correction instead of using pre-computed INTENSE significance. Default is False.

• pval_thr (float, optional) – P-value threshold for significance testing. Default is 0.05. Only used if override_intense_significance=True.

• multicomp_correction (str or None, optional) – Multiple comparison correction method. Default is None (no correction). Options: None, ‘bonferroni’, ‘holm’, ‘fdr_bh’ Only used if override_intense_significance=True.

• significance_update (bool, optional) – If True, update the significance tables with new thresholds. If False (default), only use new thresholds for current query. Only used if override_intense_significance=True.

Returns:

Dictionary with neuron IDs as keys and lists of significant features as values

Return type:

dict

store_embedding(embedding, method_name, data_type='calcium', metadata=None)

Store dimensionality reduction embedding in the experiment.

This method stores a computed embedding in the experiment’s internal embeddings dictionary. Previous embeddings with the same method_name and data_type will be overwritten without warning.

Parameters:

• embedding (np.ndarray) – The embedding array, shape (n_timepoints, n_components). The number of timepoints must match self.n_frames or self.n_frames//ds if downsampling was used.

• method_name (str) – Name of the DR method (e.g., ‘pca’, ‘umap’, ‘isomap’). This serves as the key for storing and retrieving the embedding.

• data_type (str, optional) – Type of data used (‘calcium’ or ‘spikes’). Default is ‘calcium’.

• metadata (dict, optional) – Additional metadata about the embedding. Common keys include: - ‘ds’: Downsampling factor used - ‘n_components’: Number of components - ‘neuron_indices’: Indices of neurons used - ‘method_params’: Parameters specific to the DR method

Raises:

• ValueError – If data_type is not ‘calcium’ or ‘spikes’.

• ValueError – If embedding timepoints don’t match expected frames. The expected number is self.n_frames//ds where ds is extracted from metadata (default 1).

Notes

The embedding is stored in self.embeddings[data_type][method_name] as a dictionary containing: - ‘data’: The embedding array - ‘metadata’: The provided metadata dict (or empty dict) - ‘timestamp’: Current time when stored (np.datetime64) - ‘shape’: Shape tuple of the embedding array

Previous embeddings with the same method_name are silently overwritten.

Examples

>>> import numpy as np
>>> from driada.experiment.exp_base import Experiment
>>>
>>> # Create a simple experiment
>>> calcium_data = np.random.randn(10, 1000)
>>> exp = Experiment('test', calcium_data, None, {},
...                  {'fps': 30.0}, {}, verbose=False)
>>>
>>> # Store a PCA embedding
>>> embedding = np.random.randn(1000, 3)  # 1000 timepoints, 3 components
>>> exp.store_embedding(embedding, 'pca', metadata={'n_components': 3})
>>>
>>> # Verify storage
>>> 'pca' in exp.embeddings['calcium']
True
>>>
>>> # Store a downsampled UMAP embedding
>>> # If experiment has 1000 frames and ds=5, embedding should have 200 rows
>>> downsampled_embedding = np.random.randn(200, 2)
>>> exp.store_embedding(downsampled_embedding, 'umap',
...                    metadata={'ds': 5, 'n_neighbors': 30})

See also

get_embedding

Retrieve stored embeddings

create_embedding

Create and store embeddings in one step

create_embedding(method, n_components=2, data_type='calcium', neuron_selection=None, **dr_kwargs)

Create dimensionality reduction embedding and store it.

Notes

This method modifies the experiment’s state by storing the computed embedding. Previous embeddings with the same method name will be overwritten. The method uses MultiTimeSeries internally for data handling and applies the dimensionality reduction through the MVData interface.

Parameters:

• method (str) – DR method name (‘pca’, ‘umap’, ‘isomap’, etc.).

• n_components (int, optional) – Number of embedding dimensions. Default is 2.

• data_type (str, optional) – Type of data to use (‘calcium’ or ‘spikes’). Default is ‘calcium’.

• neuron_selection (str, list or None, optional) – How to select neurons: - None or ‘all’: Use all neurons - ‘significant’: Use only significantly selective neurons - List of integers: Use specific neuron indices

• **dr_kwargs – Additional arguments for the DR method (e.g., n_neighbors, min_dist).

Returns:

embedding – The embedding array, shape (n_timepoints, n_components).

Return type:

np.ndarray

Raises:

• ValueError – If n_components is not positive, data_type is invalid, downsampling factor ‘ds’ is not an integer, neuron indices are out of bounds, significant neurons requested without selectivity analysis, or embedding method drops timepoints.

• AttributeError – If spike data requested but not available.

Examples

>>> import numpy as np
>>> from driada.experiment.exp_base import Experiment
>>> from driada.information.info_base import TimeSeries
>>>
>>> # Create experiment with some data
>>> calcium_data = np.random.randn(25, 1000)
>>> speed = TimeSeries(np.random.rand(1000))
>>> exp = Experiment('test', calcium_data, None, {},
...                  {'fps': 30.0}, {'speed': speed}, verbose=False)
>>>
>>> # Create PCA embedding using all neurons
>>> embedding = exp.create_embedding('pca', n_components=10)
Calculating PCA embedding...
>>> embedding.shape
(1000, 10)
>>>
>>> # Create downsampled PCA with specific neurons
>>> embedding = exp.create_embedding('pca', n_components=2,
...                                neuron_selection=[0, 1, 2, 3, 4], ds=10)
Calculating PCA embedding...

See also

store_embedding

Store computed embeddings

get_embedding

Retrieve stored embeddings

get_significant_neurons

Get neurons with significant selectivity

get_embedding(method_name, data_type='calcium')

Retrieve stored embedding.

This method retrieves a previously stored dimensionality reduction embedding from the experiment’s embeddings dictionary. The returned dictionary contains the embedding data along with metadata and timestamp.

Parameters:

• method_name (str) – Name of the DR method to retrieve (e.g., ‘pca’, ‘umap’).

• data_type (str, optional) – Type of data used (‘calcium’ or ‘spikes’). Default is ‘calcium’.

Returns:

Dictionary containing: - ‘data’: The embedding array (n_timepoints, n_components) - ‘metadata’: Dict with embedding parameters and settings - ‘timestamp’: np.datetime64 when the embedding was stored - ‘shape’: Tuple with shape of the embedding array

Return type:

dict

Raises:

• ValueError – If data_type is not ‘calcium’ or ‘spikes’.

• KeyError – If no embedding found for the specified method and data type.

Notes

To see available embeddings, check exp.embeddings[data_type].keys(). The returned dictionary is a reference to the stored data, so modifications will affect the stored embedding.

Examples

>>> import numpy as np
>>> from driada.experiment.exp_base import Experiment
>>>
>>> # Create experiment and store an embedding
>>> calcium_data = np.random.randn(10, 1000)
>>> exp = Experiment('test', calcium_data, None, {},
...                  {'fps': 30.0}, {}, verbose=False)
>>>
>>> # Store an embedding first
>>> embedding = np.random.randn(1000, 3)
>>> exp.store_embedding(embedding, 'pca', metadata={'n_components': 3})
>>>
>>> # Retrieve the stored PCA embedding
>>> embedding_dict = exp.get_embedding('pca')
>>> embedding_data = embedding_dict['data']
>>> print(f"Embedding shape: {embedding_dict['shape']}")
Embedding shape: (1000, 3)
>>>
>>> # Check available embeddings before retrieval
>>> available = list(exp.embeddings['calcium'].keys())
>>> print(f"Available embeddings: {available}")
Available embeddings: ['pca']

See also

store_embedding

Store embeddings in the experiment

create_embedding

Create and store embeddings in one step

compute_rdm(items, data_type='calcium', metric='correlation', average_method='mean', use_cache=True)

Compute RDM with caching support.

Parameters:

• items (str) – Name of dynamic feature to use as condition labels

• data_type (str, default 'calcium') – Type of data to use (‘calcium’ or ‘spikes’)

• metric (str, default 'correlation') – Distance metric for RDM computation

• average_method (str, default 'mean') – How to average within conditions (‘mean’ or ‘median’)

• use_cache (bool, default True) – Whether to use cached results

Returns:

• rdm (np.ndarray) – Representational dissimilarity matrix

• labels (np.ndarray) – The unique labels/conditions

clear_rdm_cache()

Clear the representational dissimilarity matrix (RDM) cache.

Removes all cached RDM computations to free memory or force recalculation with updated data. This is necessary after modifying the underlying neural data or when memory usage is a concern.

Notes

The RDM cache stores previously computed dissimilarity matrices to avoid expensive recomputation. Clear the cache when: - Neural data has been modified or reprocessed - Embeddings have been updated - Memory usage needs to be reduced - You want to force fresh computation with different parameters

After clearing, subsequent calls to compute_rdm() will recalculate the RDM from scratch, which may be computationally expensive for large datasets.

See also

compute_rdm

Method that uses and populates the RDM cache

Examples

>>> import numpy as np
>>> from driada.experiment.exp_base import Experiment
>>> from driada.information.info_base import TimeSeries
>>>
>>> # Create experiment with a categorical feature
>>> calcium_data = np.random.randn(10, 1000)
>>> conditions = TimeSeries(np.repeat([0, 1, 2], [333, 333, 334]), discrete=True)
>>> exp = Experiment('test', calcium_data, None, {},
...                  {'fps': 30.0}, {'conditions': conditions}, verbose=False)
>>>
>>> # Compute RDM (will be cached)
>>> rdm, labels = exp.compute_rdm('conditions')
>>>
>>> # Update embedding and clear cache
>>> new_embedding = np.random.randn(1000, 3)
>>> exp.store_embedding(new_embedding, 'pca')
>>> exp.clear_rdm_cache()
>>>
>>> # Verify cache is empty
>>> len(exp._rdm_cache)
0

class driada.experiment.neuron.Neuron(cell_id, ca, sp, default_t_rise=0.25, default_t_off=2.0, fps=20, fit_individual_t_off=False, optimize_kinetics=False, asp=None, wvt_ridges=None, seed=None, reconstructed=None, metrics=None)

Bases: object

Neural calcium and spike data processing.

This class handles calcium imaging time series data and spike trains, providing methods for preprocessing, spike-calcium deconvolution, and various shuffling techniques for statistical testing.

Parameters:

• cell_id (int or str) – Unique identifier for the neuron

• ca (array-like) – Calcium imaging time series data. Must be 1D array of finite values.

• sp (array-like, optional) – Spike train data (binary array where 1 indicates spike). If None, spike-related methods will not be available.

• default_t_rise (float, default=0.25) – Rise time constant in seconds for calcium transients. Must be positive.

• default_t_off (float, default=2.0) – Decay time constant in seconds for calcium transients. Must be positive.

• fps (float, default=20.0) – Sampling rate in frames per second. Must be positive.

• fit_individual_t_off (bool, default=False) – Whether to fit decay time for this specific neuron using optimization.

cell_id

The neuron identifier

Type:

int or str

t_rise

Rise time constant in frames (optimized, set by get_kinetics())

Type:

float

t_off

Decay time constant in frames (optimized, set by get_kinetics())

Type:

float

fps

Sampling rate in frames per second

Type:

float

ca_ts

Preprocessed calcium time series

Type:

TimeSeries

sp_ts

Spike train time series (if spikes provided)

Type:

TimeSeries or None

Notes

The class assumes spike data is binary (0 or 1 values). Non-binary spike data may produce incorrect results in spike counting.

static spike_form(t, t_rise, t_off)

Calculate normalized calcium response kernel shape.

Computes the double-exponential kernel used to model calcium indicator dynamics with separate rise and decay time constants.

Parameters:

• t (array-like) – Time points (in frames). Must be non-negative.

• t_rise (float) – Rise time constant (in frames). Must be positive.

• t_off (float) – Decay time constant (in frames). Must be positive.

Returns:

Normalized kernel values with peak = 1.

Return type:

ndarray

Raises:

ValueError – If t_rise or t_off are not positive.

Notes

The kernel has the form: (1 - exp(-t/τ_rise)) * exp(-t/τ_off) normalized to have maximum value of 1.

Examples

>>> t = np.linspace(0, 100, 100)
>>> kernel = spike_form(t, t_rise=5, t_off=20)
>>> kernel.max()  # Should be 1.0
1.0

static get_restored_calcium(sp, t_rise, t_off)

Reconstruct calcium signal from spike train.

Convolves spike train with double-exponential kernel to simulate calcium indicator dynamics. The output has the same length as the input spike train by truncating the convolution tail.

Parameters:

• sp (array-like) – Spike train. Can be binary (0/1) or amplitude-weighted (float). For best reconstruction fidelity, use neuron.asp.data with amplitude information. Must be 1D array.

• t_rise (float) – Rise time constant (in frames). Must be positive.

• t_off (float) – Decay time constant (in frames). Must be positive.

Returns:

Reconstructed calcium signal with same length as sp.

Return type:

ndarray

Raises:

ValueError – If t_rise or t_off are not positive, or if sp is empty.

Notes

Uses FFT-based convolution for optimal performance. Kernel length is adaptive: max(500, 5 × t_off frames) to ensure complete kernel capture for all decay time constants.

The convolution naturally handles amplitude-weighted spikes, where each spike value represents event strength in dF/F0 units.

Examples

>>> sp = np.zeros(1000)
>>> sp[100] = 1.0  # Single spike
>>> ca = get_restored_calcium(sp, t_rise=5, t_off=20)
>>> ca.shape
(1000,)

static ca_mse_error(t_off, ca, spk, t_rise)

Calculate RMSE between observed calcium and reconstructed from spikes.

This function is designed to be used with scipy.optimize.minimize, hence the parameter order with t_off first.

Parameters:

• t_off (float) – Decay time constant (in frames). Must be positive.

• ca (array-like) – Observed calcium signal. Must be 1D.

• spk (array-like) – Spike train. Must be 1D with same length as ca.

• t_rise (float) – Rise time constant (in frames). Must be positive.

Returns:

Root mean square error between observed and reconstructed calcium.

Return type:

float

Raises:

ValueError – If arrays have different lengths or time constants are invalid.

Notes

Parameter order (t_off first) is optimized for scipy.optimize.minimize where t_off is the parameter being optimized.

Examples

>>> ca = np.random.random(1000)
>>> spk = np.zeros(1000)
>>> spk[100:900:100] = 1
>>> error = ca_mse_error(t_off=20, ca=ca, spk=spk, t_rise=5)
>>> error > 0
True

static calcium_preprocessing(ca, seed=None)

Preprocess calcium signal for spike reconstruction.

Applies preprocessing steps: - Converts to float64 for numerical stability - Clips negative values to 0 (calcium cannot be negative) - Adds tiny noise to prevent numerical singularities

Parameters:

• ca (array-like) – Raw calcium signal. Must be 1D.

• seed (int, optional) – Random seed for reproducible noise. If None, uses current state.

Returns:

Preprocessed calcium signal as float64 array.

Return type:

ndarray

Raises:

ValueError – If ca is empty.

Notes

The small noise (1e-8 scale) prevents division by zero and other numerical issues in downstream spike reconstruction algorithms.

Examples

>>> ca = np.array([1.0, -0.5, 2.0, 0.5])
>>> processed = calcium_preprocessing(ca, seed=42)
>>> processed[1]  # Negative value clipped to ~0
0.0...
>>> (processed > 0).all()  # All values positive after noise
True

static extract_event_amplitudes(ca_signal, st_ev_inds, end_ev_inds, baseline_window=20, already_dff=False, baseline_offset=0, baseline_offset_sec=None, use_peak_refinement=False, t_rise_frames=None, t_off_frames=None, fps=None, peak_search_window_sec=0.2)

Extract amplitudes from calcium signal for detected events.

For raw fluorescence signals, applies dF/F0 normalization following Neugornet et al. 2021 (PMC8032960). For pre-normalized dF/F signals, directly extracts peak values.

Parameters:

• ca_signal (ndarray) – Calcium signal (1D). Can be raw fluorescence or pre-normalized dF/F.

• st_ev_inds (list of int) – Start indices for each event.

• end_ev_inds (list of int) – End indices for each event.

• baseline_window (int, optional) – Number of frames before event to use for F0 calculation. Default is 20 frames.

• already_dff (bool, optional) – If True, signal is already dF/F normalized and peak values are extracted directly. If False, applies dF/F0 normalization. Default is False for backward compatibility.

• baseline_offset (int, optional) – Number of frames to skip before baseline window (to avoid rise phase). Baseline is calculated from frames [start-baseline_window-baseline_offset, start-baseline_offset]. Default is 0 (baseline immediately before event). NOTE: Prefer baseline_offset_sec for FPS-independent specification.

• baseline_offset_sec (float, optional) – Alternative to baseline_offset: offset in seconds (converted to frames using fps). Useful for GCaMP indicators with non-zero rise time (~0.25s for GCaMP6f). If provided, overrides baseline_offset. Requires fps parameter.

• use_peak_refinement (bool, optional) – If True, refine peak location before calculating baseline (similar to ‘onset’ placement). Ensures baseline is calculated relative to TRUE peak, not wavelet detection start. Critical for overlapping events where wavelet detection timing is imprecise. Default is False for backward compatibility.

• t_rise_frames (float, optional) – Rise time in frames. Required if use_peak_refinement=True.

• t_off_frames (float, optional) – Decay time in frames. Required if use_peak_refinement=True.

• fps (float, optional) – Sampling rate in Hz. Required if use_peak_refinement=True.

• peak_search_window_sec (float, optional) – Time window in seconds to search for refined peak location. Only used if use_peak_refinement=True. Default is 0.2s.

Returns:

amplitudes – Amplitude for each event (dF/F0 if already_dff=False, peak dF/F if True).

Return type:

list of float

Notes

For events at the start where baseline_window extends before 0, uses available data. Returns 0 amplitude for events with invalid F0 (zero or negative baseline when already_dff=False).

When use_peak_refinement=True, searches for TRUE peak in expanded window (similar to ‘onset’ placement logic) and calculates baseline relative to refined peak instead of wavelet detection start. This is critical for overlapping events where baseline must avoid elevated signal from previous events.

References

Neugornet A, O’Donovan B, Ortinski PI (2021). Comparative Effects of Event Detection Methods on the Analysis and Interpretation of Ca2+ Imaging Data. Front Neurosci 15:620869.

static deconvolve_given_event_times(ca_signal, event_times, t_rise_frames, t_off_frames, event_mask=None)

Extract amplitudes via non-negative least squares deconvolution.

Given detected event times and known calcium kernel parameters, finds the optimal amplitudes that best reconstruct the observed signal. This solves: signal ≈ Σᵢ aᵢ · kernel(t - tᵢ)

Handles overlapping events naturally by jointly optimizing all amplitudes to minimize reconstruction error.

Parameters:

• ca_signal (ndarray) – Calcium signal (1D array, already dF/F normalized).

• event_times (array-like) – Event onset times as frame indices.

• t_rise_frames (float) – Rise time of calcium kernel in frames.

• t_off_frames (float) – Decay time of calcium kernel in frames.

• event_mask (ndarray of bool, optional) – Binary mask indicating frames to fit (typically expanded event regions). If provided, NNLS only fits these frames, reducing baseline noise contribution while covering full calcium transients. Default None fits all frames (legacy behavior).

Returns:

amplitudes – Optimal amplitudes for each event (non-negative).

Return type:

ndarray

Notes

Uses scipy.optimize.nnls for non-negative least squares solution.

When event_mask is provided, only masked frames contribute to the fit. The mask should cover full calcium transients (rise + decay), not just peaks. This reduces baseline noise while maintaining proper transient fitting.

References

Lawson CL, Hanson RJ (1995). Solving Least Squares Problems. SIAM, Philadelphia.

static compute_kernel_peak_offset(t_rise_frames, t_off_frames)

Compute the time offset from onset to peak for double-exponential kernel.

The kernel (1 - exp(-t/t_rise)) * exp(-t/t_off) starts at 0, rises to a peak, then decays. This function returns the time (in frames) from onset (t=0) to peak.

Parameters:

• t_rise_frames (float) – Rise time constant in frames.

• t_off_frames (float) – Decay time constant in frames.

Returns:

Time offset from onset to peak in frames. Returns 0 if kinetics are degenerate (t_rise ≈ t_off).

Return type:

float

static estimate_onset_times(ca_signal, st_inds, end_inds, t_rise_frames, t_off_frames)

Estimate spike onset times from detected event boundaries.

Event detectors (wavelet, threshold) find the rising/peak region of calcium transients, not the true spike onset. This function estimates onset times by finding the peak within each event and back-calculating using the kernel peak offset.

Parameters:

• ca_signal (ndarray) – Calcium signal (1D array).

• st_inds (list of int) – Event start frame indices (from detector).

• end_inds (list of int) – Event end frame indices (from detector).

• t_rise_frames (float) – Rise time constant in frames.

• t_off_frames (float) – Decay time constant in frames.

Returns:

Estimated onset times for each event.

Return type:

list of int

static amplitudes_to_point_events(length, ca_signal, st_ev_inds, end_ev_inds, amplitudes, placement='peak', t_rise_frames=None, t_off_frames=None, fps=None, peak_search_window_sec=0.2)

Convert event boundaries and amplitudes to point event array.

Stores amplitudes at specific positions as delta functions. Temporal spreading is handled by convolution with calcium kernel.

Parameters:

• length (int) – Length of output array (n_frames).

• ca_signal (ndarray) – Calcium signal (used for ‘peak’ placement).

• st_ev_inds (list of int) – Start indices for each event.

• end_ev_inds (list of int) – End indices for each event.

• amplitudes (list of float) – Amplitude (dF/F0) for each event.

• placement ({'start', 'peak', 'onset', 'onset_refined'}, optional) – Where to place amplitude: - ‘start’: at event start index (as detected by wavelet) - ‘peak’: at actual calcium peak within event window - ‘onset’: at estimated spike onset (peak - kernel_peak_offset), using peak within event boundaries - ‘onset_refined’: at estimated spike onset with peak refinement search outside event boundaries (legacy, causes lag on noisy data)

• t_rise_frames (float, optional) – Rise time in frames. Required for ‘onset’ placement.

• t_off_frames (float, optional) – Decay time in frames. Required for ‘onset’ placement.

• fps (float, optional) – Sampling rate in Hz. Required for ‘onset’ placement with peak refinement.

• peak_search_window_sec (float, optional) – Time window in seconds to search for refined peak location around wavelet-detected peak. Only used for ‘onset’ placement. Default is 0.2s.

Returns:

point_events – Array of zeros with amplitudes at event positions (1D, float). Sparse continuous array suitable for TimeSeries(discrete=False).

Return type:

ndarray

Notes

If multiple events have centers/peaks at same index, amplitudes are summed. This is rare but can occur with closely spaced events.

The ‘onset’ placement finds peak within event boundaries, then calculates the kernel peak offset to place spikes at estimated onset. This trusts wavelet detection and is robust to noise.

The ‘onset_refined’ placement (legacy) searches outside event boundaries for peak refinement. This can cause temporal lag on noisy data and is not recommended for real calcium imaging data.

property t_rise

Optimized rise time constant (frames).

Returns None if not yet optimized. Setter validates and clears cached metrics.

__init__(cell_id, ca, sp, default_t_rise=0.25, default_t_off=2.0, fps=20, fit_individual_t_off=False, optimize_kinetics=False, asp=None, wvt_ridges=None, seed=None, reconstructed=None, metrics=None)

Initialize Neuron object with calcium and spike data.

Parameters:

• cell_id (str or int) – Unique identifier for this neuron.

• ca (array-like) – Calcium fluorescence signal. Must be 1D.

• sp (array-like or None) – Binary spike train. If provided, must have same length as ca.

• default_t_rise (float, optional) – Default rise time constant in seconds. Default is DEFAULT_T_RISE.

• default_t_off (float, optional) – Default decay time constant in seconds. Default is DEFAULT_T_OFF.

• fps (float, optional) – Sampling rate in Hz. Must be positive. Default is DEFAULT_FPS.

• fit_individual_t_off (bool, optional) – DEPRECATED: Use optimize_kinetics instead. If True, fit individual decay time using old method. Default is False.

• optimize_kinetics (bool or str, optional) – If True or ‘direct’, optimize kinetics parameters (t_rise, t_off) using direct measurement from event shapes. If False, use default parameters. Default is False. Requires either asp parameter or prior call to reconstruct_spikes().

• asp (array-like or None, optional) – Pre-computed amplitude spikes (from prior reconstruction). If provided, enables kinetics optimization without re-running reconstruction.

• wvt_ridges (list or None, optional) – Pre-computed wavelet ridges (from prior reconstruction). Used for kinetics optimization to provide event boundaries.

• seed (int, optional) – Random seed for preprocessing reproducibility.

• reconstructed (array-like or None, optional) – Pre-computed reconstructed calcium signal. If provided, must have same length as ca. Skips reconstruction step and directly populates the internal reconstructed trace.

• metrics (dict or None, optional) – Pre-computed quality metrics dictionary. Recognized keys include ‘t_rise’ and ‘t_off’ (kinetics in seconds, converted to frames), ‘r2_score’ and ‘snr_recon’ (reconstruction quality), ‘event_r2_score’ and ‘event_snr’ (event-level quality), and ‘noise_level’ (noise amplitude). Unrecognized keys are stored but otherwise ignored.

Raises:

• ValueError – If ca is empty, fps is not positive, or array lengths don’t match. If optimize_kinetics=True but no spike data available.

• TypeError – If ca cannot be converted to numeric array.

Notes

The shuffle mask excludes CA_SHIFT_N_TOFF * t_off frames from each end to prevent near-zero circular shifts from contaminating the null distribution.

New workflow (recommended): >>> neuron = Neuron(cell_id=0, ca=calcium, sp=None, fps=20) >>> neuron.reconstruct_spikes(method=’wavelet’) >>> neuron.get_kinetics() # Optimize after reconstruction

Or pass pre-computed data: >>> neuron = Neuron(cell_id=0, ca=calcium, sp=None, fps=20, … asp=amp_spikes, wvt_ridges=ridges, … optimize_kinetics=True)

property t_off

Optimized decay time constant (frames).

Returns None if not yet optimized. Setter validates and clears cached metrics.

reconstruct_spikes(method='wavelet', iterative=True, n_iter=3, min_events_threshold=2, adaptive_thresholds=False, amplitude_method='deconvolution', show_progress=False, create_event_regions=False, event_mask_expansion_sec=5.0, wavelet=None, rel_wvt_times=None, use_gpu=False, **kwargs)

Reconstruct spikes from calcium signal.

Reconstructs discrete spike events from continuous calcium fluorescence traces using wavelet or threshold-based methods.

Parameters:

• method (str, optional) – Reconstruction method: ‘wavelet’ or ‘threshold’. Default is ‘wavelet’.

• iterative (bool, optional) – Use iterative wavelet detection with residual analysis (wavelet method only). Detects events in residuals across multiple iterations to handle overlapping events and improve detection of smaller events. Default is True.

• n_iter (int, optional) – Number of iterations (only if iterative=True). Default is 3.

• min_events_threshold (int, optional) – Stop iterating if fewer events detected (only if iterative=True). Default is 2.

• adaptive_thresholds (bool, optional) – Progressively relax detection thresholds across iterations (only if iterative=True). Default is False.

• amplitude_method (str, optional) – Method for extracting event amplitudes: ‘peak’ or ‘deconvolution’. - ‘peak’: Peak-based extraction with baseline subtraction (backward compatible) - ‘deconvolution’: Non-negative least squares deconvolution (default, optimal for overlapping events) Default is ‘deconvolution’.

• show_progress (bool, optional) – Whether to show progress bar during wavelet detection. Default is False (no progress bar for single neuron processing).

• create_event_regions (bool, optional) – If True, creates self.events (binary rectangular regions marking event durations). This is legacy behavior - modern code should use self.asp (amplitude spikes) instead. Only needed for backward compatibility or specific visualization needs. Default is False to avoid unnecessary computation and warnings.

• event_mask_expansion_sec (float, optional) – Time in seconds to expand event mask around detected events for NNLS deconvolution. The mask is expanded by ±event_mask_expansion_sec to cover the full calcium transient (rise + decay). Larger values include more of the decay but also more baseline noise. Default is 5.0 seconds (optimal balance for GCaMP6s with t_off ~2s).

• wavelet (Wavelet, optional) – Pre-computed wavelet object for batch processing optimization. If None, will be created by extract_wvt_events(). Default is None.

• rel_wvt_times (array-like, optional) – Pre-computed time resolutions for batch processing optimization. If None, will be computed by extract_wvt_events(). Default is None.

• use_gpu (bool, default=False) – Whether to use GPU acceleration for wavelet transform computation. Requires PyTorch and CuPy. Ridge extraction remains CPU-only.

• **kwargs –

  Additional parameters depend on method:

  For ‘wavelet’:

  • fpsfloat, optional

    Sampling rate in Hz. Default is DEFAULT_FPS (20.0 Hz).

  • min_event_durfloat, optional

    Minimum event duration in seconds. Default is 0.5.

  • max_event_durfloat, optional

    Maximum event duration in seconds. Default is 2.5.

  For ‘threshold’:

  • threshold_stdfloat, optional

    Number of standard deviations above mean. Default is 2.5.

  • smooth_sigmafloat, optional

    Gaussian smoothing sigma in frames. Default is 2.

  • min_spike_intervalfloat, optional

    Minimum interval between spikes in seconds. Default is 0.1.

Returns:

If create_event_regions=True: Binary event regions with shape (n_frames,). Values are 0 or 1. If create_event_regions=False: None (no binary regions created). Always populates self.asp (amplitude spikes) and self.sp (binary spikes) attributes. Optionally populates self.events (binary regions) if create_event_regions=True.

Return type:

ndarray or None

Raises:

• NotImplementedError – If method is not ‘wavelet’ or ‘threshold’.

• ValueError – If parameters are invalid for the chosen method.

Notes

The wavelet method uses continuous wavelet transform to detect calcium transient events. The threshold method uses derivative-based spike detection with Gaussian smoothing.

Iterative detection (default) performs multiple detection passes on residuals, removing detected events and searching for additional events in the remaining signal. This approach significantly improves detection of overlapping and smaller events compared to single-pass detection.

For single-pass detection (backward compatible), set iterative=False.

get_mad()

Get median absolute deviation of calcium signal.

Computes MAD as a robust measure of noise level in the calcium signal. Caches the result for efficiency.

Returns:

Median absolute deviation of the calcium signal, scaled to be consistent with standard deviation for normally distributed data.

Return type:

float

Notes

MAD is more robust to outliers than standard deviation, making it ideal for noise estimation in calcium imaging data which often contains spike-related transients.

get_snr(method='simple')

Get signal-to-noise ratio of calcium signal.

Unified interface for SNR calculation supporting multiple methods.

Parameters:

method ({'simple', 'wavelet'}, optional) –

Calculation method. Options:

• ’simple’: Peak-based SNR using spike locations (fast, default)

• ’wavelet’: Event-based SNR using wavelet regions (accurate)

Default is ‘simple’.

Returns:

Signal-to-noise ratio (dimensionless). Higher values indicate stronger signal relative to noise.

Return type:

float

Raises:

ValueError – If no spike/event data available, or if method is invalid.

Notes

Simple method:

• SNR = mean(calcium at spike peaks) / MAD(entire signal)

• Fast computation, uses asp (amplitude spikes) or sp (binary spikes)

• Caches result in self.snr

Wavelet method:

• SNR = median(event amplitudes) / std(baseline)

• More accurate, empirically validated against ground truth

• Requires prior wavelet reconstruction

• Caches result in self.wavelet_snr

Examples

>>> neuron.get_snr()                    # Simple SNR (default)
>>> neuron.get_snr(method='simple')     # Explicit simple
>>> neuron.get_snr(method='wavelet')    # Wavelet-based SNR

get_reconstruction_r2(event_only=False, n_mad=3.0, use_detected_events=True)

Get R² for calcium reconstruction quality.

Computes the coefficient of determination (R²) measuring how well the double-exponential model fits the observed calcium signal. Higher values indicate better model fit. Standard R² is cached.

Parameters:

• event_only (bool, default=False) – If True, compute R² only in detected event regions. If False, compute standard R² over entire signal (cached).

• n_mad (float, default=3.0) – Number of MAD (Median Absolute Deviation) units above median for event detection. Only used when event_only=True and use_detected_events=False.

• use_detected_events (bool, default=True) – If True, use self.events (wavelet-detected event regions) for event mask. If False, use MAD-based threshold on original signal (legacy behavior). Only used when event_only=True.

Returns:

R² value between -inf and 1. Values closer to 1 indicate better fit. R² > 0.9: Excellent fit (high-quality calcium transients) R² 0.7-0.9: Good fit (acceptable quality) R² 0.5-0.7: Moderate fit (check for artifacts) R² < 0.5: Poor fit (likely artifacts or model mismatch)

Return type:

float

Raises:

ValueError – If amplitude spike data is not available, or if use_detected_events=True but self.events is None.

Notes

Standard R² = 1 - (SS_residual / SS_total) computed over entire signal. Event R² = 1 - (SS_residual_events / SS_total_events) in event regions only.

When use_detected_events=True (recommended), event regions are defined by self.events (wavelet ridge detection), ensuring alignment with reconstruction.

When use_detected_events=False (legacy), event regions are defined by MAD threshold on the original signal.

Uses cached RMSE from get_noise_ampl() for efficiency in standard mode. Event-only mode is NOT cached as it depends on parameters.

get_snr_reconstruction()

Get reconstruction-based SNR.

Computes SNR as the ratio of calcium signal standard deviation to reconstruction error (RMSE). Provides a model-based quality metric complementary to the peak-based SNR. Caches the result.

Returns:

Reconstruction SNR value (positive). Higher is better. SNR_recon > 20: Excellent reconstruction fidelity SNR_recon 10-20: Good reconstruction fidelity SNR_recon 5-10: Fair reconstruction fidelity SNR_recon < 5: Poor reconstruction (check signal quality)

Return type:

float

Raises:

ValueError – If spike data is not available for reconstruction.

Notes

SNR_reconstruction = std(Ca) / RMSE

This metric differs from get_snr() which uses peak amplitudes at spike times. Reconstruction SNR measures overall model fit quality.

Uses cached RMSE from get_noise_ampl() for efficiency.

get_reconstruction_scaled(t_rise=None, t_off=None)

Get reconstruction transformed using ca.scdata’s scaler.

Applies the exact same MinMaxScaler transformation to the reconstruction that was used to create ca.scdata. This ensures consistency for information theory and dimensionality reduction analyses that use ca.scdata.

Caches result in self._reconstructed_scaled for efficiency. Cache is cleared when ASP or kinetics change (via _clear_cached_metrics).

Parameters:

• t_rise (float, optional) – Rise time in seconds. If None, uses optimized or default value.

• t_off (float, optional) – Decay time in seconds. If None, uses optimized or default value.

Returns:

Reconstruction scaled using ca_scaler.transform(). Shape matches ca.data. Note: Values may fall outside [0,1] if reconstruction differs from training data (ca.data) in amplitude or baseline.

Return type:

ndarray

Raises:

ValueError – If spike data (asp) is not available.

Notes

Scaling Behavior:

Uses scaler.transform() NOT scaler.fit_transform(): - Fitted on ca.data during __init__ - Applied to reconstruction (which may exceed ca.data range) - Result can fall outside [0,1] - this is EXPECTED and CORRECT

When values exceed [0,1]: - recon > data_max → scaled > 1.0 (e.g., missed event in original) - recon < data_min → scaled < 0.0 (e.g., baseline noise, wrong kinetics)

This is the mathematically correct behavior for scaler.transform() and preserves the relative scaling needed for information theory.

Use Cases: - Information theory analyses (mutual information with ca.scdata) - Dimensionality reduction (PCA/ICA on ca.scdata-like signals) - Comparing reconstruction vs ca.scdata in scaled space

DO NOT use for metrics: R², MAE, RMSE should use ca.data (original scale).

Examples

>>> neuron.reconstruct_spikes()
>>> recon_scaled = neuron.get_reconstruction_scaled()
>>>
>>> # Check range violations
>>> pct_below = 100 * np.sum(recon_scaled < 0) / len(recon_scaled)
>>> pct_above = 100 * np.sum(recon_scaled > 1.0) / len(recon_scaled)
>>> print(f'Below 0: {pct_below:.1f}%, Above 1: {pct_above:.1f}%')
>>>
>>> # Use with dimensionality reduction
>>> from sklearn.decomposition import PCA
>>> pca = PCA(n_components=5)
>>> pca.fit(ca.scdata_matrix)  # Fit on scaled data
>>> recon_proj = pca.transform(recon_scaled.reshape(1, -1))

See also

get_reconstruction_r2

R² using ca.data (for metrics)

get_mae

MAE using ca.data (for metrics)

get_mae()

Get Mean Absolute Error between observed and reconstructed calcium.

Computes MAE as the mean absolute deviation between observed calcium signal and reconstruction from detected spikes. Provides intuitive measure of reconstruction quality in original signal units (ΔF/F). Caches the result.

Returns:

MAE value (non-negative). Lower is better. MAE < 0.1: Excellent reconstruction (typical for high SNR) MAE 0.1-0.2: Good reconstruction (acceptable quality) MAE 0.2-0.3: Moderate reconstruction (check for issues) MAE > 0.3: Poor reconstruction (likely artifacts or failures)

Return type:

float

Raises:

ValueError – If amplitude spike data is not available.

Notes

MAE = mean(|Ca_observed - Ca_reconstructed|)

Unlike RMSE, MAE treats all errors equally (no squaring). Useful for understanding typical deviation and detecting outlier sensitivity:

• If RMSE >> MAE: Few large errors dominate (e.g., missed events)

• If RMSE ≈ MAE: Errors uniformly distributed (e.g., white noise)

The MAE/RMSE ratio can diagnose error distribution patterns.

get_event_rmse(n_mad=4.0, use_detected_events=True)

Get RMSE during event periods only.

Measures reconstruction error only during calcium transient events, ignoring baseline noise. Provides a more accurate measure of deconvolution quality than full-signal RMSE, which is dominated by baseline regions that reconstruction should NOT fit.

Parameters:

• n_mad (float, default=4.0) – Number of MAD (Median Absolute Deviation) units above median for event detection threshold. Only used if use_detected_events=False.

• use_detected_events (bool, default=True) – If True, use self.events (wavelet-detected) for event mask. If False, use MAD-based threshold (legacy behavior).

Returns:

Event RMSE (lower is better). Same units as calcium signal (ΔF/F). Event RMSE < 0.05: Excellent event reconstruction Event RMSE 0.05-0.10: Good event reconstruction Event RMSE 0.10-0.15: Moderate event reconstruction Event RMSE > 0.15: Poor event reconstruction

Return type:

float

Raises:

ValueError – If amplitude spike data is not available or no events detected.

Notes

Event_RMSE = sqrt(mean((Ca_events - Recon_events)^2))

Where events are defined as: Ca > median + n_mad * MAD

Complementary to Event R² - while Event R² shows proportion of variance explained, Event RMSE shows absolute reconstruction error magnitude.

get_event_mae(n_mad=4.0, use_detected_events=True)

Get MAE during event periods only.

Measures reconstruction error only during calcium transient events, ignoring baseline noise. Provides a more accurate measure of deconvolution quality than full-signal MAE, which is dominated by baseline regions that reconstruction should NOT fit.

Parameters:

• n_mad (float, default=4.0) – Number of MAD (Median Absolute Deviation) units above median for event detection threshold. Only used if use_detected_events=False.

• use_detected_events (bool, default=True) – If True, use self.events (wavelet-detected) for event mask. If False, use MAD-based threshold (legacy behavior).

Returns:

Event MAE (lower is better). Same units as calcium signal (ΔF/F). Event MAE < 0.05: Excellent event reconstruction Event MAE 0.05-0.10: Good event reconstruction Event MAE 0.10-0.15: Moderate event reconstruction Event MAE > 0.15: Poor event reconstruction

Return type:

float

Raises:

ValueError – If amplitude spike data is not available or no events detected.

Notes

Event_MAE = mean(|Ca_events - Recon_events|)

Where events are defined as: Ca > median + n_mad * MAD

Unlike Event RMSE, Event MAE does not square errors, making it less sensitive to outliers. Useful for understanding typical deviation:

• If Event_RMSE >> Event_MAE: Few large errors (missed/false events)

• If Event_RMSE ≈ Event_MAE: Errors uniformly distributed

get_event_count(n_mad=None)

Get count of detected spike events (ridges passing thresholds).

Counts the number of spike events detected by wavelet reconstruction that pass all thresholds. Returns count of non-zero entries in amplitude spike data. Result is cached for efficiency.

Parameters:

n_mad (float, default=3.0) – Not currently used - reserved for future threshold filtering. Kept for API consistency with other event-based methods.

Returns:

Number of detected spike events (ridges).

Return type:

int

Raises:

ValueError – If spike reconstruction has not been performed.

Notes

Counts non-zero entries in self.asp.data (amplitude spike array). Each non-zero entry represents a detected calcium event that passed wavelet ridge filtering and amplitude thresholds during reconstruction.

Useful for: - Assessing detection sensitivity - Computing event-based statistics - Validating reconstruction parameters

get_baseline_noise_std(n_mad=3.0)

Get baseline noise standard deviation from residual signal.

Estimates noise level as the standard deviation of reconstruction residuals (original - reconstructed). This measures the actual reconstruction error, including measurement noise, unmodeled dynamics, and reconstruction inaccuracies.

Parameters:

n_mad (float, default=3.0) – Not used in this implementation. Kept for API compatibility.

Returns:

Standard deviation of residual signal.

Return type:

float

Raises:

ValueError – If reconstruction has not been performed.

Notes

Residual = Ca_original - Ca_reconstructed

The residual contains: - Measurement noise - Unmodeled small events - Drift and baseline fluctuations - Reconstruction errors

This is the proper noise estimate for normalized metrics (NMAE, NRMSE).

get_wavelet_snr()

Get event-based signal-to-noise ratio.

Uses detected event regions to separate signal from baseline, providing accurate SNR measurement that accounts for event timing and shape. Works with both wavelet and threshold detection methods.

Returns:

Event SNR (signal_strength / baseline_noise). Higher values indicate better signal quality.

Return type:

float

Raises:

ValueError – If reconstruction not performed with create_event_regions=True, no events detected, insufficient data (< 3 events), or baseline noise is zero.

Notes

Requires prior call to:

neuron.reconstruct_spikes(method=’wavelet’, create_event_regions=True) OR neuron.reconstruct_spikes(method=’threshold’, create_event_regions=True)

SNR calculation: 1. Baseline: median and MAD from non-event frames 2. Signal: median of peak amplitudes across all events 3. SNR = (signal - baseline_median) / baseline_noise

Uses peak amplitudes (not event medians) to correctly handle sparse high-amplitude events.

See also

get_snr

Simple SNR based on spike times

get_event_snr

Alias for this method

reconstruct_spikes

Spike reconstruction (wavelet or threshold)

Examples

>>> neuron = Neuron(cell_id=0, ca=calcium_data, sp=None)
>>> neuron.reconstruct_spikes(method='threshold', create_event_regions=True)
>>> snr = neuron.get_event_snr()  # or get_wavelet_snr()
>>> print(f"Signal quality (SNR): {snr:.2f}")

get_event_snr()

Get event-based signal-to-noise ratio.

Uses detected event regions to separate signal from baseline, providing accurate SNR measurement that accounts for event timing and shape. Works with both wavelet and threshold detection methods.

Returns:

Event SNR (signal_strength / baseline_noise). Higher values indicate better signal quality.

Return type:

float

Raises:

ValueError – If reconstruction not performed with create_event_regions=True, no events detected, insufficient data (< 3 events), or baseline noise is zero.

Notes

Requires prior call to:

neuron.reconstruct_spikes(method=’wavelet’, create_event_regions=True) OR neuron.reconstruct_spikes(method=’threshold’, create_event_regions=True)

SNR calculation: 1. Baseline: median and MAD from non-event frames 2. Signal: median of peak amplitudes across all events 3. SNR = (signal - baseline_median) / baseline_noise

Uses peak amplitudes (not event medians) to correctly handle sparse high-amplitude events.

See also

get_snr

Simple SNR based on spike times

get_event_snr

Alias for this method

reconstruct_spikes

Spike reconstruction (wavelet or threshold)

Examples

>>> neuron = Neuron(cell_id=0, ca=calcium_data, sp=None)
>>> neuron.reconstruct_spikes(method='threshold', create_event_regions=True)
>>> snr = neuron.get_event_snr()  # or get_wavelet_snr()
>>> print(f"Signal quality (SNR): {snr:.2f}")

get_nmae(n_mad=3.0)

Get Normalized Mean Absolute Error (MAE / baseline_noise_std).

Computes MAE divided by baseline noise standard deviation. This normalization accounts for varying noise levels across neurons, enabling fair quality comparison. Values indicate error magnitude relative to baseline noise level.

Parameters:

n_mad (float, default=3.0) – Number of MAD units above median for baseline detection.

Returns:

Normalized MAE (dimensionless ratio). Lower is better. NMAE < 1.0: Error less than baseline noise (excellent) NMAE < 2.0: Error ~2× baseline noise (good) NMAE < 3.0: Error ~3× baseline noise (moderate) NMAE > 3.0: Error >> baseline noise (poor)

Return type:

float

Raises:

ValueError – If baseline noise cannot be estimated or is zero.

Notes

NMAE = MAE / std(baseline)

Unlike raw MAE, NMAE is comparable across neurons with different noise levels. Interpretation: error magnitude in “noise units”.

get_nrmse(n_mad=3.0)

Get Normalized RMSE (RMSE / baseline_noise_std).

Computes RMSE divided by baseline noise standard deviation. This normalization accounts for varying noise levels across neurons. More sensitive to large errors than NMAE due to squaring.

Parameters:

n_mad (float, default=3.0) – Number of MAD units above median for baseline detection.

Returns:

Normalized RMSE (dimensionless ratio). Lower is better. NRMSE < 1.0: Error less than baseline noise (excellent) NRMSE < 2.0: Error ~2× baseline noise (good) NRMSE < 3.0: Error ~3× baseline noise (moderate) NRMSE > 3.0: Error >> baseline noise (poor)

Return type:

float

Raises:

ValueError – If baseline noise cannot be estimated or is zero.

Notes

NRMSE = RMSE / std(baseline)

NRMSE > NMAE indicates few large errors dominate (e.g., missed events). NRMSE ≈ NMAE indicates uniformly distributed errors.

get_kinetics(method='direct', fps=20, use_cached=True, update_reconstruction=True, **kwargs)

Get optimized calcium kinetics parameters (t_rise, t_off).

Simple wrapper around optimize_kinetics() for easy access to optimized parameters. Caches results for efficiency - only runs optimization once per neuron.

Updates neuron attributes: - self.t_rise: optimized rise time (frames) - self.t_off: optimized decay time (frames) - self._kinetics_info: full optimization results dict

Parameters:

• method (str, optional) – Optimization method. Currently only ‘direct’ is supported. Default: ‘direct’.

• fps (float, optional) – Sampling rate in frames per second. Default: 20.0.

• use_cached (bool, optional) – If True and optimization already run, return cached results. If False, always re-run optimization. Default: True.

• update_reconstruction (bool, optional) – If True, recompute reconstruction with optimized kinetics. Default: True.

• **kwargs (dict, optional) – Additional arguments passed to optimize_kinetics() (e.g., t_rise_range, t_off_range, ftol, gtol, maxiter, etc.)

Returns:

Dictionary with keys: - ‘t_rise’: float, rise time in seconds - ‘t_off’: float, decay time in seconds - ‘t_rise_frames’: float, rise time in frames - ‘t_off_frames’: float, decay time in frames - ‘event_r2’: float, event R² with optimized parameters - ‘default_r2’: float, event R² with default parameters - ‘improvement’: float, R² improvement - ‘method’: str, optimization method used - (additional fields from optimize_kinetics())

Return type:

dict

Raises:

ValueError – If no events detected (call reconstruct_spikes() first).

Examples

>>> # Optimize kinetics and access via attributes
>>> neuron = Neuron(cell_id=0, ca=calcium_trace, sp=None, fps=20)
>>> neuron.reconstruct_spikes(method='wavelet')
>>> result = neuron.get_kinetics(fps=20)
>>>
>>> # Access optimized parameters via attributes (in frames)
>>> t_rise_frames = neuron.t_rise
>>> t_off_frames = neuron.t_off
>>>
>>> # Or via returned dict (in seconds)
>>> t_rise_sec = result['t_rise']
>>> t_off_sec = result['t_off']
>>> print(f"R² improvement: {result['improvement']:+.3f}")
>>>
>>> # Access full optimization info
>>> print(f"Converged: {neuron._kinetics_info['converged']}")
>>> print(f"Time: {neuron._kinetics_info['time']:.2f}s")

Notes

• Results are cached in self._kinetics_info after first call

• Optimized parameters stored in self.t_rise and self.t_off (frames)

• Uses fast direct derivative measurement method

• For backward compatibility: neuron.t_off now contains optimized value

get_t_off()

Get calcium decay time constant.

Deprecated since version 0.5.0: Use get_kinetics() instead for better optimization that jointly optimizes both t_rise and t_off using correct event R² metric. This method only fits t_off using simple MSE and is significantly slower.

Fits the decay time constant by optimizing the match between observed calcium and reconstructed calcium from spikes. Caches the result for efficiency.

Returns:

Decay time constant in frames.

Return type:

float

Raises:

ValueError – If spike data is not available for fitting.

Notes

Uses scipy.optimize.minimize to find the optimal t_off value that minimizes the RMSE between observed and reconstructed calcium.

DEPRECATED: Use get_kinetics() instead:

>>> # Old way (deprecated)
>>> t_off = neuron.get_t_off()

>>> # New way (recommended)
>>> kinetics = neuron.get_kinetics()
>>> t_off = kinetics['t_off_frames']  # in frames
>>> t_rise = kinetics['t_rise_frames']  # bonus: also get t_rise

get_noise_ampl()

Get noise amplitude estimate from calcium-spike reconstruction.

Returns the root mean square error (RMSE) between observed calcium and reconstructed calcium from spikes using current kinetics parameters. This provides an estimate of the noise level in the calcium signal after accounting for spike-related transients.

Returns:

RMSE between observed and reconstructed calcium signal.

Return type:

float

Raises:

ValueError – If spike data is not available for reconstruction.

Notes

Uses cached reconstruction from get_reconstructed() which respects optimized kinetics (self.t_rise, self.t_off). This is different from MAD, as it specifically measures the residual error after spike-to-calcium reconstruction. The value is cached after first computation.

get_shuffled_calcium(method='roll_based', return_array=True, seed=None, **kwargs)

Get shuffled calcium signal using various randomization methods.

Creates surrogate data that preserves certain statistical properties of the original calcium signal while destroying temporal relationships.

Parameters:

• method ({'roll_based', 'waveform_based', 'chunks_based'}, optional) – Shuffling method to use: - ‘roll_based’: Circular shift by random offset - ‘waveform_based’: Shuffle spikes then reconstruct calcium - ‘chunks_based’: Divide signal into chunks and reorder Default is ‘roll_based’.

• return_array (bool, optional) – If True, return numpy array. If False, return TimeSeries object. Default is True.

• seed (int, optional) – Random seed for reproducible shuffling.

• **kwargs – Additional arguments passed to shuffling method: - For chunks_based: n (int) - number of chunks

Returns:

Shuffled calcium signal with same length as original.

Return type:

ndarray or TimeSeries

Raises:

AttributeError – If the specified method does not exist.

Notes

Different methods preserve different signal properties: - roll_based: Preserves all autocorrelations - waveform_based: Preserves spike waveform shapes - chunks_based: Preserves local signal structure within chunks

get_shuffled_spikes(method='isi_based', return_array=True, seed=None, **kwargs)

Get shuffled spike train.

Creates surrogate spike data that preserves certain statistical properties while destroying temporal relationships.

Parameters:

• method ({'isi_based'}, optional) – Shuffling method to use. Currently only ‘isi_based’ is supported, which preserves inter-spike interval statistics. Default is ‘isi_based’.

• return_array (bool, optional) – If True, return numpy array. If False, return TimeSeries object. Default is True.

• seed (int, optional) – Random seed for reproducible shuffling.

• **kwargs – Additional arguments passed to shuffling method.

Returns:

Shuffled spike train with same number of spikes as original.

Return type:

ndarray or TimeSeries

Raises:

• AttributeError – If no spike data is available.

• ValueError – If method is not recognized.

Notes

The ISI-based method preserves the distribution of inter-spike intervals while randomizing spike positions.

property reconstructed

Get reconstructed calcium signal from amplitude spikes (cached).

Convenience property that returns cached reconstruction. For more control, use get_reconstructed() method instead.

Returns:

Reconstructed calcium signal as TimeSeries object. Returns None if amplitude spike data is not available.

Return type:

TimeSeries or None

See also

get_reconstructed

Method with force_reconstruction and custom parameters

get_reconstructed(force_reconstruction=False, **kwargs)

Get reconstructed calcium signal from amplitude spikes.

Lazily computes and caches the reconstructed calcium signal by convolving detected amplitude spikes with the calcium kernel.

Parameters:

• force_reconstruction (bool, optional) – If True, bypass cache and force recomputation. Default is False.

• **kwargs (dict, optional) –

  Custom reconstruction parameters:

  • t_rise_framesfloat, optional

    Custom rise time in frames. If not provided, uses optimized t_rise from get_kinetics().

  • t_off_framesfloat, optional

    Custom decay time in frames. If not provided, uses optimized t_off from get_kinetics().

  • spike_dataarray-like, optional

    Custom spike data to reconstruct from. If not provided, uses self.asp.data.

Returns:

Reconstructed calcium signal as TimeSeries object. Returns None if amplitude spike data is not available (call reconstruct_spikes() first) and no custom spike_data provided.

Return type:

TimeSeries or None

Notes

The default reconstruction uses: - self.asp.data: Amplitude spikes (dF/F units) - self.t_rise: Optimized rise time from get_kinetics() (lazy, cached) - self.t_off: Optimized decay time from get_kinetics() (lazy, cached)

The reconstructed signal can be compared with self.ca to assess reconstruction quality. Use get_reconstruction_r2(), get_mae(), or get_snr_reconstruction() for quantitative metrics.

Custom parameters are NOT cached. Only default reconstruction is cached.

Examples

>>> neuron = Neuron(cell_id=1, ca=calcium_data, sp=None, fps=20)
>>> neuron.reconstruct_spikes(method='wavelet')
>>>
>>> # Get cached default reconstruction
>>> recon = neuron.reconstructed
>>>
>>> # Force recomputation
>>> recon = neuron.get_reconstructed(force_reconstruction=True)
>>>
>>> # Use custom decay time
>>> recon_fast = neuron.get_reconstructed(t_off_frames=20)
>>> recon_slow = neuron.get_reconstructed(t_off_frames=60)

optimize_kinetics(method='direct', fps=20, update_reconstruction=True, max_event_dur_multiplier=4, detection_method='auto', **kwargs)

Universal kinetics optimization with automatic reconstruction update.

Single entry point for all optimization methods. Automatically updates instance kinetics (self.t_rise, self.t_off) and optionally re-runs spike detection with new parameters.

Parameters:

• method (str, optional) – Optimization method. Currently only ‘direct’ is supported. Default: ‘direct’

• fps (float, optional) – Sampling rate in frames per second. Default: 20.0

• update_reconstruction (bool, optional) – If True, automatically re-run spike detection with optimized kinetics. This ensures events match the new parameters. Default: True.

• max_event_dur_multiplier (float, optional) – Multiplier for calculating max_event_dur when update_reconstruction=True. Formula: max_event_dur = t_rise + multiplier * t_off. Higher values detect longer events but may merge overlapping events. Lower values improve precision but may miss event tails. Recommended range: 3.0-5.0. Default: 4.0 (optimal balance).

• detection_method ({'auto', 'wavelet', 'threshold'}, optional) –

  Method to use for event re-detection when update_reconstruction=True:

  • ’auto’: Use threshold if threshold_events exist, else wavelet (default)

  • ’wavelet’: Always use wavelet detection (slower, more sensitive)

  • ’threshold’: Always use threshold detection (faster, requires high SNR)

  Default: ‘auto’

• **kwargs (dict) –

  Method-specific parameters including:

  • min_r2 : float, minimum R² for t_off fit quality (default: 0.8). Events with poor exponential fit are rejected.

  See _optimize_kinetics_direct() for full list.

Returns:

Optimization results with keys:

• ’optimized’: bool, whether optimization succeeded

• ’t_rise’: float, optimized rise time (seconds)

• ’t_off’: float, optimized decay time (seconds)

• ’method’: str, method used

• Additional method-specific metrics

Return type:

dict

Examples

>>> # Fast threshold-based workflow (100-500x faster)
>>> neuron.detect_events_threshold(n_mad=4.0)
>>> result = neuron.optimize_kinetics(method='direct', fps=30)
>>> # Auto-detects threshold mode, re-runs threshold detection with optimized kinetics

>>> # Explicit threshold mode for iterative refinement
>>> neuron.detect_events_threshold(n_mad=4.0)
>>> result = neuron.optimize_kinetics(
...     method='direct', fps=30, detection_method='threshold'
... )

>>> # Traditional wavelet workflow (slower but more sensitive)
>>> neuron.reconstruct_spikes(method='wavelet')
>>> result = neuron.optimize_kinetics(method='direct', fps=30)
>>> # Uses wavelet detection for re-detection

>>> # Skip auto-reconstruction for speed (e.g., in batch processing)
>>> result = neuron.optimize_kinetics(
...     method='direct', fps=30, update_reconstruction=False
... )
>>> # Later: neuron.reconstruct_spikes(method='wavelet')

Notes

• Consistently updates self.t_rise/t_off and reconstructs events

• Setting update_reconstruction=False allows manual control of reconstruction timing

detect_events_threshold(threshold=None, n_mad=4.0, min_duration_frames=3, merge_gap_frames=2, use_scaled=True, signal=None)

Detect calcium events using threshold crossings (fast alternative to wavelet).

This method finds event boundaries by detecting when the calcium signal crosses above/below a threshold. Returns SimpleEvent objects compatible with optimize_kinetics(), providing ~100-500x speedup vs wavelet detection.

Parameters:

• threshold (float, optional) – Absolute threshold value. If None, computed as: median(signal) + n_mad * MAD(signal) where MAD = median absolute deviation (robust noise estimate).

• n_mad (float, optional) – Number of MAD units above median for auto-threshold. Only used if threshold=None. Default: 4.0 (robust detection).

• min_duration_frames (int, optional) – Minimum event duration in frames. Events shorter than this are discarded. Default: 3 frames.

• merge_gap_frames (int, optional) – Merge events separated by fewer than this many frames. Prevents event fragmentation. Default: 2 frames.

• use_scaled (bool, optional) – If True, use scaled calcium data (self.ca.scdata). If False, use raw calcium data (self.ca.data). Default: True (recommended for consistent thresholds).

• signal (ndarray, optional) – Custom signal to use for detection. If provided, use_scaled is ignored. Used for iterative detection on residual signals.

Returns:

Detected events with .start and .end attributes (frame indices). Empty list if no events detected.

Return type:

list of SimpleEvent

Raises:

• AttributeError – If calcium data not available or lacks scdata attribute.

• ValueError – If parameters are invalid (negative values, percentile out of range).

Examples

>>> # Automatic threshold (recommended)
>>> neuron = Neuron(cell_id=0, ca=calcium_data, sp=None)
>>> events = neuron.detect_events_threshold(n_mad=4.0)
>>> len(events)
42

>>> # Then use with optimize_kinetics for fast kinetics estimation
>>> neuron.threshold_events = events  # Store for optimize_kinetics
>>> result = neuron.optimize_kinetics(method='direct', fps=20)

>>> # Manual threshold
>>> events = neuron.detect_events_threshold(threshold=0.3, use_scaled=True)

>>> # More sensitive detection (lower threshold)
>>> events = neuron.detect_events_threshold(n_mad=3.0)

>>> # Access event properties
>>> for event in events[:3]:
...     print(f"Event: frames {event.start:.0f}-{event.end:.0f}, duration={event.duration:.0f}")

Notes

Performance: O(N) complexity vs O(N²) for wavelet detection - Typical speedup: 100-500x faster - Example: 1000 frames: ~0.01s (threshold) vs 1-5s (wavelet)

Algorithm: 1. Compute threshold (auto or manual) 2. Find upward crossings (signal goes above threshold) → event starts 3. Find downward crossings (signal goes below threshold) → event ends 4. Filter by minimum duration 5. Merge events with small gaps

Comparison with wavelet detection: - Threshold: Fast, simple, good for high SNR data - Wavelet: Slower, more sensitive, better for low SNR or overlapping events

The detected events are stored in self.threshold_events for later use.

See also

optimize_kinetics

Use detected events for kinetics estimation

reconstruct_spikes

Wavelet-based spike detection (slower but more sensitive)

Usage Examples

Creating an Experiment

from driada.experiment import Experiment, load_demo_experiment
from driada.information import MultiTimeSeries
import numpy as np

# Load sample experiment
exp = load_demo_experiment()

# Access and modify experiment properties
print(f"Sampling rate: {exp.fps} Hz")

# Add behavior data as dynamic features
n_timepoints = exp.n_frames  # Get number of timepoints from experiment
position = np.random.randn(2, n_timepoints)  # x, y coordinates
exp.dynamic_features['position'] = position

# Access metadata
# Note: exp_identificators is set during experiment creation
# and contains experiment parameters like mouse_id, session, etc.

Working with Neurons

from driada.experiment import Neuron, load_demo_experiment

# Load sample experiment
exp = load_demo_experiment()

# Create neuron analyzer for first cell
# Neuron takes calcium data and spike data directly
neuron = Neuron(
    cell_id=0,
    ca=exp.calcium.scdata[0],  # calcium trace for neuron 0
    sp=None,  # spike data (optional)
    fps=exp.fps
)

# Access neuron's data
ca_trace = neuron.ca  # Raw calcium data
sp_trace = neuron.sp  # Spike data (if provided)

Experiment Properties

from driada.experiment import load_demo_experiment

# Load sample experiment
exp = load_demo_experiment()

# Core data access
calcium_matrix = exp.calcium.data  # (n_neurons, n_time)
sampling_rate = exp.fps

# Duration and time
duration = exp.n_frames / exp.fps  # compute duration in seconds
time_vector = np.arange(exp.n_frames) / exp.fps  # compute time vector

# Neuron count
n_neurons = exp.n_cells

# Check for spike data
if exp.spikes is not None:
    spikes = exp.spikes.data

Data Organization

The Experiment class organizes different data modalities:

• calcium: MultiTimeSeries of calcium imaging data

• spikes: MultiTimeSeries of spike data (optional)

• behavior: Dictionary of behavioral variables

• info: Metadata dictionary

• fps: Sampling rate in Hz

Time alignment is automatic - all time series are assumed to be synchronously sampled at the specified frame rate.