29037613829

Committed 09 Jul 2026 05:36PM UTC coverage: 61.575% (-3.3%) from 64.888%

Build # 29037613829

Build Type

push

github

Committed by

web-flow

Commit Message

Merge 88653ad51 into 087557586

Coverage Stats

1822 of 3048 new or added lines in 81 files covered. (59.78%)

3091 existing lines in 75 files now uncovered.

12243 of 19883 relevant lines covered (61.58%)

0.62 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

44.79

/brainbox/processing.py

"""Process data from one form into another.

For example, taking spike times and binning them into non-overlapping bins and convolving spike
times with a gaussian kernel.
"""

import numpy as np
import pandas as pd
from scipy import interpolate, sparse
from brainbox import core
from iblutil.numerical import bincount2D
from iblutil.util import Bunch
import logging

_logger = logging.getLogger(__name__)


def sync(dt, times=None, values=None, timeseries=None, offsets=None, interp='zero', fillval=np.nan):
    """
    Function for resampling a single or multiple time series to a single, evenly-spaced, delta t
    between observations. Uses interpolation to find values.

    Can be used on raw numpy arrays of timestamps and values using the 'times' and 'values' kwargs
    and/or on brainbox.core.TimeSeries objects passed to the 'timeseries' kwarg. If passing both
    TimeSeries objects and numpy arrays, the offsets passed should be for the TS objects first and
    then the numpy arrays.

    Uses scipy's interpolation library to perform interpolation.
    See scipy.interp1d for more information regarding interp and fillval parameters.

    :param dt: Separation of points which the output timeseries will be sampled at
    :type dt: float
    :param timeseries: A group of time series to perform alignment or a single time series.
        Must have time stamps.
    :type timeseries: tuple of TimeSeries objects, or a single TimeSeries object.
    :param times: time stamps for the observations in 'values']
    :type times: np.ndarray or list of np.ndarrays
    :param values: observations corresponding to the timestamps in 'times'
    :type values: np.ndarray or list of np.ndarrays
    :param offsets: tuple of offsets for time stamps of each time series. Offsets for passed
        TimeSeries objects first, then offsets for passed numpy arrays. defaults to None
    :type offsets: tuple of floats, optional
    :param interp: Type of interpolation to use. Refer to scipy.interpolate.interp1d for possible
        values, defaults to np.nan
    :type interp: str
    :param fillval: Fill values to use when interpolating outside of range of data. See interp1d
        for possible values, defaults to np.nan
    :return: TimeSeries object with each row representing synchronized values of all
        input TimeSeries. Will carry column names from input time series if all of them have column
        names.
    """
    #########################################
    # Checks on inputs and input processing #
    #########################################

    # Initialize a list to contain times/values pairs if no TS objs are passed
    if timeseries is None:
        timeseries = []
    # If a single time series is passed for resampling, wrap it in an iterable
    elif isinstance(timeseries, core.TimeSeries):
        timeseries = [timeseries]
    # Yell at the user if they try to pass stuff to timeseries that isn't a TimeSeries object
    elif not all([isinstance(ts, core.TimeSeries) for ts in timeseries]):
        raise TypeError(
            "All elements of 'timeseries' argument must be brainbox.core.TimeSeries "
            "objects. Please uses 'times' and 'values' for np.ndarray args."
        )
    # Check that if something is passed to times or values, there is a corresponding equal-length
    # argument for the other element.
    if (times is not None) or (values is not None):
        if len(times) != len(values):
            raise ValueError("'times' and 'values' must have the same number of elements.")
        if type(times[0]) is np.ndarray:
            if not all([t.shape == v.shape for t, v in zip(times, values)]):
                raise ValueError("All arrays in 'times' must match the shape of the corresponding entry in 'values'.")
            # If all checks are passed, convert all times and values args into TimeSeries objects
            timeseries.extend([core.TimeSeries(t, v) for t, v in zip(times, values)])
        else:
            # If times and values are only numpy arrays and lists of arrays, pair them and add
            timeseries.append(core.TimeSeries(times, values))

    # Adjust each timeseries by the associated offset if necessary then load into a list
    if offsets is not None:
        tstamps = [ts.times + os for ts, os in zip(timeseries, offsets)]
    else:
        tstamps = [ts.times for ts in timeseries]
    # If all input timeseries have column names, put them together for the output TS
    if all([ts.columns is not None for ts in timeseries]):
        colnames = []
        for ts in timeseries:
            colnames.extend(ts.columns)
    else:
        colnames = None

    #################
    # Main function #
    #################

    # Get the min and max values for all timeseries combined after offsetting
    tbounds = np.array([(np.amin(ts), np.amax(ts)) for ts in tstamps])
    if not np.all(np.isfinite(tbounds)):
        # If there is a np.inf or np.nan in the time stamps for any of the timeseries this will
        # break any further code so we check for all finite values and throw an informative error.
        raise ValueError(
            'NaN or inf encountered in passed timeseries.\
                          Please either drop or fill these values.'
        )
    tmin, tmax = np.amin(tbounds[:, 0]), np.amax(tbounds[:, 1])
    if fillval == 'extrapolate':
        # If extrapolation is enabled we can ensure we have a full coverage of the data by
        # extending the t max to be an whole integer multiple of dt above tmin.
        # The 0.01% fudge factor is to account for floating point arithmetic errors.
        newt = np.arange(tmin, tmax + 1.0001 * (dt - (tmax - tmin) % dt), dt)
    else:
        newt = np.arange(tmin, tmax, dt)
    tsinterps = [interpolate.interp1d(ts.times, ts.values, kind=interp, fill_value=fillval, axis=0) for ts in timeseries]
    syncd = core.TimeSeries(newt, np.hstack([tsi(newt) for tsi in tsinterps]), columns=colnames)
    return syncd


def compute_cluster_average(spike_clusters, spike_var):
    """
    Quickish way to compute the average of some quantity across spikes in each cluster given
    quantity for each spike

    :param spike_clusters: cluster idx of each spike
    :param spike_var: variable of each spike (e.g spike amps or spike depths)
    :return: cluster id, average of quantity for each cluster, no. of spikes per cluster
    """
    clust, inverse, counts = np.unique(spike_clusters, return_inverse=True, return_counts=True)
    _spike_var = sparse.csr_matrix((spike_var, (inverse, np.zeros(inverse.size, dtype=int))))
    spike_var_avg = np.ravel(_spike_var.toarray()) / counts

    return clust, spike_var_avg, counts


def bin_spikes(spikes, binsize, interval_indices=False):
    """
    Wrapper for bincount2D which is intended to take in a TimeSeries object of spike times
    and cluster identities and spit out spike counts in bins of a specified width binsize, also in
    another TimeSeries object. Can either return a TS object with each row labeled with the
    corresponding interval or the value of the left edge of the bin.

    :param spikes: Spike times and cluster identities of sorted spikes
    :type spikes: TimeSeries object with \'clusters\' column and timestamps
    :param binsize: Width of the non-overlapping bins in which to bin spikes
    :type binsize: float
    :param interval_indices: Whether to use intervals as the time stamps for binned spikes, rather
        than the left edge value of the bins, defaults to False
    :type interval_indices: bool, optional
    :return: Object with 2D array of shape T x N, for T timesteps and N clusters, and the
        associated time stamps.
    :rtype: TimeSeries object
    """
    if type(spikes) is not core.TimeSeries:
        raise TypeError('Input spikes need to be in TimeSeries object format')

    if not hasattr(spikes, 'clusters'):
        raise AttributeError(
            "Input spikes need to have a clusters attribute. Make sure you set columns=('clusters',)) when constructing spikes."
        )

    rates, tbins, clusters = bincount2D(spikes.times, spikes.clusters, binsize)
    if interval_indices:
        intervals = pd.interval_range(tbins[0], tbins[-1], freq=binsize, closed='left')
        return core.TimeSeries(times=intervals, values=rates.T[:-1], columns=clusters)
    else:
        return core.TimeSeries(times=tbins, values=rates.T, columns=clusters)


def get_units_bunch(spks_b, *args):
    """
    Returns a bunch, where the bunch keys are keys from `spks` with labels of spike information
    (e.g. unit IDs, times, features, etc.), and the values for each key are arrays with values for
    each unit: these arrays are ordered and can be indexed by unit id.

    Parameters
    ----------
    spks_b : bunch
        A spikes bunch containing fields with spike information (e.g. unit IDs, times, features,
        etc.) for all spikes.
    features : list of strings (optional positional arg)
        A list of names of labels of spike information (which must be keys in `spks`) that specify
        which labels to return as keys in `units`. If not provided, all keys in `spks` are returned
        as keys in `units`.

    Returns
    -------
    units_b : bunch
        A bunch with keys of labels of spike information (e.g. cluster IDs, times, features, etc.)
        whose values are arrays that hold values for each unit. The arrays for each key are ordered
        by unit ID.

    Examples
    --------
    1) Create a units bunch given a spikes bunch, and get the amps for unit #4 from the units
    bunch.
        >>> from brainbox import processing
        >>> import one.alf.io as alfio
        >>> import ibllib.ephys.spikes as e_spks
        (*Note, if there is no 'alf' directory, make 'alf' directory from 'ks2' output directory):
        >>> e_spks.ks2_to_alf(path_to_ks_out, path_to_alf_out)
        >>> spks_b = alfio.load_object(path_to_alf_out, 'spikes')
        >>> units_b = processing.get_units_bunch(spks_b)
        # Get amplitudes for unit 4.
        >>> amps = units_b['amps']['4']

    TODO add computation time estimate?
    """

    # Initialize `units`
    units_b = Bunch()
    # Get the keys to return for `units`:
    if not args:
        feat_keys = list(spks_b.keys())
    else:
        feat_keys = args[0]
    # Get unit id for each spike and number of units. *Note: `n_units` might not equal `len(units)`
    # because some clusters may be empty (due to a "wontfix" bug in ks2).
    spks_unit_id = spks_b['clusters']
    n_units = np.max(spks_unit_id)
    units = np.unique(spks_b['clusters'])
    # For each key in `units`, iteratively get each unit's values and add as a key to a bunch,
    # `feat_bunch`. After iterating through all units, add `feat_bunch` as a key to `units`:
    for feat in feat_keys:
        # Initialize `feat_bunch` with a key for each unit.
        feat_bunch = Bunch((str(unit), np.array([])) for unit in np.arange(n_units))
        for unit in units:
            unit_idxs = np.where(spks_unit_id == unit)[0]
            feat_bunch[str(unit)] = spks_b[feat][unit_idxs]
        units_b[feat] = feat_bunch
    return units_b


def filter_units(units_b, t, **kwargs):
    """
    Filters units according to some parameters. **kwargs are the keyword parameters used to filter
    the units.

    Parameters
    ----------
    units_b : bunch
        A bunch with keys of labels of spike information (e.g. cluster IDs, times, features, etc.)
        whose values are arrays that hold values for each unit. The arrays for each key are ordered
        by unit ID.
    t : float
        Duration of time over which to calculate the firing rate and false positive rate.

    Keyword Parameters
    ------------------
    min_amp : float
        The minimum mean amplitude (in V) of the spikes in the unit. Default value is 50e-6.
    min_fr : float
        The minimum firing rate (in Hz) of the unit. Default value is 0.5.
    max_fpr : float
        The maximum false positive rate of the unit (using the fp formula in Hill et al. (2011)
        J Neurosci 31: 8699-8705). Default value is 0.2.
    rp : float
        The refractory period (in s) of the unit. Used to calculate `max_fp`. Default value is
        0.002.

    Returns
    -------
    filt_units : ndarray
        The ids of the filtered units.

    See Also
    --------
    get_units_bunch

    Examples
    --------
    1) Filter units according to the default parameters.
        >>> from brainbox import processing
        >>> import one.alf.io as alfio
        >>> import ibllib.ephys.spikes as e_spks
        (*Note, if there is no 'alf' directory, make 'alf' directory from 'ks2' output directory):
        >>> e_spks.ks2_to_alf(path_to_ks_out, path_to_alf_out)
        # Get a spikes bunch, units bunch, and filter the units.
        >>> spks_b = alfio.load_object(path_to_alf_out, 'spikes')
        >>> units_b = processing.get_units_bunch(spks_b, ['times', 'amps', 'clusters'])
        >>> T = spks_b['times'][-1] - spks_b['times'][0]
        >>> filtered_units = processing.filter_units(units_b, T)

    2) Filter units with no minimum amplitude, a minimum firing rate of 1 Hz, and a max false
    positive rate of 0.2, given a refractory period of 2 ms.
        >>> filtered_units  = processing.filter_units(units_b, T, min_amp=0, min_fr=1)

    TODO: `units_b` input arg could eventually be replaced by `clstrs_b` if the required metrics
          are in `clstrs_b['metrics']`
    """

    # Set params
    params = {'min_amp': 50e-6, 'min_fr': 0.5, 'max_fpr': 0.2, 'rp': 0.002}  # defaults
    params.update(kwargs)  # update from **kwargs

    # Iteratively filter the units for each filter param #
    # -------------------------------------------------- #
    units = np.asarray(list(units_b.amps.keys()))
    # Remove empty clusters
    empty_cl = np.where([len(units_b.amps[unit]) == 0 for unit in units])[0]
    filt_units = np.delete(units, empty_cl)
    for param in params.keys():
        if param == 'min_amp':  # return units above with amp > `'min_amp'`
            mean_amps = np.asarray([np.mean(units_b.amps[unit]) for unit in filt_units])
            filt_idxs = np.where(mean_amps > params['min_amp'])[0]
            filt_units = filt_units[filt_idxs]
        elif param == 'min_fr':  # return units with fr > `'min_fr'`
            fr = np.asarray([len(units_b.amps[unit]) / (units_b.times[unit][-1] - units_b.times[unit][0]) for unit in filt_units])
            filt_idxs = np.where(fr > params['min_fr'])[0]
            filt_units = filt_units[filt_idxs]
        elif param == 'max_fpr':  # return units with fpr < `'max_fpr'`
            fpr = np.zeros_like(filt_units, dtype='float')
            for i, unit in enumerate(filt_units):
                n_spks = len(units_b.amps[unit])
                n_isi_viol = len(np.where(np.diff(units_b.times[unit]) < params['rp'])[0])
                # fpr is min of roots of solved quadratic equation (Hill, et al. 2011).
                c = (t * n_isi_viol) / (2 * params['rp'] * n_spks**2)  # 3rd term in quadratic
                fpr[i] = np.min(np.abs(np.roots([-1, 1, c])))  # solve quadratic
            filt_idxs = np.where(fpr < params['max_fpr'])[0]
            filt_units = filt_units[filt_idxs]
    return filt_units.astype(int)

1	"""Process data from one form into another.
2
3	For example, taking spike times and binning them into non-overlapping bins and convolving spike
4	times with a gaussian kernel.
5	"""
6
7	import numpy as np	1✔
8	import pandas as pd	1✔
9	from scipy import interpolate, sparse	1✔
10	from brainbox import core	1✔
11	from iblutil.numerical import bincount2D	1✔
12	from iblutil.util import Bunch	1✔
13	import logging	1✔
14
15	_logger = logging.getLogger(__name__)	1✔
16
17
18	def sync(dt, times=None, values=None, timeseries=None, offsets=None, interp='zero', fillval=np.nan):	1✔
19	"""
20	Function for resampling a single or multiple time series to a single, evenly-spaced, delta t
21	between observations. Uses interpolation to find values.
22
23	Can be used on raw numpy arrays of timestamps and values using the 'times' and 'values' kwargs
24	and/or on brainbox.core.TimeSeries objects passed to the 'timeseries' kwarg. If passing both
25	TimeSeries objects and numpy arrays, the offsets passed should be for the TS objects first and
26	then the numpy arrays.
27
28	Uses scipy's interpolation library to perform interpolation.
29	See scipy.interp1d for more information regarding interp and fillval parameters.
30
31	:param dt: Separation of points which the output timeseries will be sampled at
32	:type dt: float
33	:param timeseries: A group of time series to perform alignment or a single time series.
34	Must have time stamps.
35	:type timeseries: tuple of TimeSeries objects, or a single TimeSeries object.
36	:param times: time stamps for the observations in 'values']
37	:type times: np.ndarray or list of np.ndarrays
38	:param values: observations corresponding to the timestamps in 'times'
39	:type values: np.ndarray or list of np.ndarrays
40	:param offsets: tuple of offsets for time stamps of each time series. Offsets for passed
41	TimeSeries objects first, then offsets for passed numpy arrays. defaults to None
42	:type offsets: tuple of floats, optional
43	:param interp: Type of interpolation to use. Refer to scipy.interpolate.interp1d for possible
44	values, defaults to np.nan
45	:type interp: str
46	:param fillval: Fill values to use when interpolating outside of range of data. See interp1d
47	for possible values, defaults to np.nan
48	:return: TimeSeries object with each row representing synchronized values of all
49	input TimeSeries. Will carry column names from input time series if all of them have column
50	names.
51	"""
52	#########################################
53	# Checks on inputs and input processing #
54	#########################################
55
56	# Initialize a list to contain times/values pairs if no TS objs are passed
57	if timeseries is None:	1✔
58	timeseries = []	1✔
59	# If a single time series is passed for resampling, wrap it in an iterable
60	elif isinstance(timeseries, core.TimeSeries):	1✔
61	timeseries = [timeseries]	1✔
62	# Yell at the user if they try to pass stuff to timeseries that isn't a TimeSeries object
63	elif not all([isinstance(ts, core.TimeSeries) for ts in timeseries]):	1✔
NEW 64	raise TypeError(	×
65	"All elements of 'timeseries' argument must be brainbox.core.TimeSeries "
66	"objects. Please uses 'times' and 'values' for np.ndarray args."
67	)
68	# Check that if something is passed to times or values, there is a corresponding equal-length
69	# argument for the other element.
70	if (times is not None) or (values is not None):	1✔
71	if len(times) != len(values):	1✔
NEW 72	raise ValueError("'times' and 'values' must have the same number of elements.")	×
73	if type(times[0]) is np.ndarray:	1✔
74	if not all([t.shape == v.shape for t, v in zip(times, values)]):	1✔
NEW 75	raise ValueError("All arrays in 'times' must match the shape of the corresponding entry in 'values'.")	×
76	# If all checks are passed, convert all times and values args into TimeSeries objects
77	timeseries.extend([core.TimeSeries(t, v) for t, v in zip(times, values)])	1✔
78	else:
79	# If times and values are only numpy arrays and lists of arrays, pair them and add
80	timeseries.append(core.TimeSeries(times, values))	1✔
81
82	# Adjust each timeseries by the associated offset if necessary then load into a list
83	if offsets is not None:	1✔
84	tstamps = [ts.times + os for ts, os in zip(timeseries, offsets)]	×
85	else:
86	tstamps = [ts.times for ts in timeseries]	1✔
87	# If all input timeseries have column names, put them together for the output TS
88	if all([ts.columns is not None for ts in timeseries]):	1✔
89	colnames = []	1✔
90	for ts in timeseries:	1✔
91	colnames.extend(ts.columns)	1✔
92	else:
93	colnames = None	1✔
94
95	#################
96	# Main function #
97	#################
98
99	# Get the min and max values for all timeseries combined after offsetting
100	tbounds = np.array([(np.amin(ts), np.amax(ts)) for ts in tstamps])	1✔
101	if not np.all(np.isfinite(tbounds)):	1✔
102	# If there is a np.inf or np.nan in the time stamps for any of the timeseries this will
103	# break any further code so we check for all finite values and throw an informative error.
NEW 104	raise ValueError(	×
105	'NaN or inf encountered in passed timeseries.\
106	Please either drop or fill these values.'
107	)
108	tmin, tmax = np.amin(tbounds[:, 0]), np.amax(tbounds[:, 1])	1✔
109	if fillval == 'extrapolate':	1✔
110	# If extrapolation is enabled we can ensure we have a full coverage of the data by
111	# extending the t max to be an whole integer multiple of dt above tmin.
112	# The 0.01% fudge factor is to account for floating point arithmetic errors.
113	newt = np.arange(tmin, tmax + 1.0001 * (dt - (tmax - tmin) % dt), dt)	1✔
114	else:
UNCOV 115	newt = np.arange(tmin, tmax, dt)	×
116	tsinterps = [interpolate.interp1d(ts.times, ts.values, kind=interp, fill_value=fillval, axis=0) for ts in timeseries]	1✔
117	syncd = core.TimeSeries(newt, np.hstack([tsi(newt) for tsi in tsinterps]), columns=colnames)	1✔
118	return syncd	1✔
119
120
121	def compute_cluster_average(spike_clusters, spike_var):	1✔
122	"""
123	Quickish way to compute the average of some quantity across spikes in each cluster given
124	quantity for each spike
125
126	:param spike_clusters: cluster idx of each spike
127	:param spike_var: variable of each spike (e.g spike amps or spike depths)
128	:return: cluster id, average of quantity for each cluster, no. of spikes per cluster
129	"""
130	clust, inverse, counts = np.unique(spike_clusters, return_inverse=True, return_counts=True)	1✔
131	_spike_var = sparse.csr_matrix((spike_var, (inverse, np.zeros(inverse.size, dtype=int))))	1✔
132	spike_var_avg = np.ravel(_spike_var.toarray()) / counts	1✔
133
134	return clust, spike_var_avg, counts	1✔
135
136
137	def bin_spikes(spikes, binsize, interval_indices=False):	1✔
138	"""
139	Wrapper for bincount2D which is intended to take in a TimeSeries object of spike times
140	and cluster identities and spit out spike counts in bins of a specified width binsize, also in
141	another TimeSeries object. Can either return a TS object with each row labeled with the
142	corresponding interval or the value of the left edge of the bin.
143
144	:param spikes: Spike times and cluster identities of sorted spikes
145	:type spikes: TimeSeries object with \'clusters\' column and timestamps
146	:param binsize: Width of the non-overlapping bins in which to bin spikes
147	:type binsize: float
148	:param interval_indices: Whether to use intervals as the time stamps for binned spikes, rather
149	than the left edge value of the bins, defaults to False
150	:type interval_indices: bool, optional
151	:return: Object with 2D array of shape T x N, for T timesteps and N clusters, and the
152	associated time stamps.
153	:rtype: TimeSeries object
154	"""
UNCOV 155	if type(spikes) is not core.TimeSeries:	×
156	raise TypeError('Input spikes need to be in TimeSeries object format')	×
157
UNCOV 158	if not hasattr(spikes, 'clusters'):	×
NEW 159	raise AttributeError(	×
160	"Input spikes need to have a clusters attribute. Make sure you set columns=('clusters',)) when constructing spikes."
161	)
162
UNCOV 163	rates, tbins, clusters = bincount2D(spikes.times, spikes.clusters, binsize)	×
UNCOV 164	if interval_indices:	×
165	intervals = pd.interval_range(tbins[0], tbins[-1], freq=binsize, closed='left')	×
166	return core.TimeSeries(times=intervals, values=rates.T[:-1], columns=clusters)	×
167	else:
168	return core.TimeSeries(times=tbins, values=rates.T, columns=clusters)	×
169
170
171	def get_units_bunch(spks_b, *args):	1✔
172	"""
173	Returns a bunch, where the bunch keys are keys from `spks` with labels of spike information
174	(e.g. unit IDs, times, features, etc.), and the values for each key are arrays with values for
175	each unit: these arrays are ordered and can be indexed by unit id.
176
177	Parameters
178	----------
179	spks_b : bunch
180	A spikes bunch containing fields with spike information (e.g. unit IDs, times, features,
181	etc.) for all spikes.
182	features : list of strings (optional positional arg)
183	A list of names of labels of spike information (which must be keys in `spks`) that specify
184	which labels to return as keys in `units`. If not provided, all keys in `spks` are returned
185	as keys in `units`.
186
187	Returns
188	-------
189	units_b : bunch
190	A bunch with keys of labels of spike information (e.g. cluster IDs, times, features, etc.)
191	whose values are arrays that hold values for each unit. The arrays for each key are ordered
192	by unit ID.
193
194	Examples
195	--------
196	1) Create a units bunch given a spikes bunch, and get the amps for unit #4 from the units
197	bunch.
198	>>> from brainbox import processing
199	>>> import one.alf.io as alfio
200	>>> import ibllib.ephys.spikes as e_spks
201	(*Note, if there is no 'alf' directory, make 'alf' directory from 'ks2' output directory):
202	>>> e_spks.ks2_to_alf(path_to_ks_out, path_to_alf_out)
203	>>> spks_b = alfio.load_object(path_to_alf_out, 'spikes')
204	>>> units_b = processing.get_units_bunch(spks_b)
205	# Get amplitudes for unit 4.
206	>>> amps = units_b['amps']['4']
207
208	TODO add computation time estimate?
209	"""
210
211	# Initialize `units`
UNCOV 212	units_b = Bunch()	×
213	# Get the keys to return for `units`:
214	if not args:	×
UNCOV 215	feat_keys = list(spks_b.keys())	×
216	else:
217	feat_keys = args[0]	×
218	# Get unit id for each spike and number of units. *Note: `n_units` might not equal `len(units)`
219	# because some clusters may be empty (due to a "wontfix" bug in ks2).
UNCOV 220	spks_unit_id = spks_b['clusters']	×
UNCOV 221	n_units = np.max(spks_unit_id)	×
222	units = np.unique(spks_b['clusters'])	×
223	# For each key in `units`, iteratively get each unit's values and add as a key to a bunch,
224	# `feat_bunch`. After iterating through all units, add `feat_bunch` as a key to `units`:
UNCOV 225	for feat in feat_keys:	×
226	# Initialize `feat_bunch` with a key for each unit.
227	feat_bunch = Bunch((str(unit), np.array([])) for unit in np.arange(n_units))	×
UNCOV 228	for unit in units:	×
229	unit_idxs = np.where(spks_unit_id == unit)[0]	×
230	feat_bunch[str(unit)] = spks_b[feat][unit_idxs]	×
231	units_b[feat] = feat_bunch	×
232	return units_b	×
233
234
235	def filter_units(units_b, t, **kwargs):	1✔
236	"""
237	Filters units according to some parameters. **kwargs are the keyword parameters used to filter
238	the units.
239
240	Parameters
241	----------
242	units_b : bunch
243	A bunch with keys of labels of spike information (e.g. cluster IDs, times, features, etc.)
244	whose values are arrays that hold values for each unit. The arrays for each key are ordered
245	by unit ID.
246	t : float
247	Duration of time over which to calculate the firing rate and false positive rate.
248
249	Keyword Parameters
250	------------------
251	min_amp : float
252	The minimum mean amplitude (in V) of the spikes in the unit. Default value is 50e-6.
253	min_fr : float
254	The minimum firing rate (in Hz) of the unit. Default value is 0.5.
255	max_fpr : float
256	The maximum false positive rate of the unit (using the fp formula in Hill et al. (2011)
257	J Neurosci 31: 8699-8705). Default value is 0.2.
258	rp : float
259	The refractory period (in s) of the unit. Used to calculate `max_fp`. Default value is
260	0.002.
261
262	Returns
263	-------
264	filt_units : ndarray
265	The ids of the filtered units.
266
267	See Also
268	--------
269	get_units_bunch
270
271	Examples
272	--------
273	1) Filter units according to the default parameters.
274	>>> from brainbox import processing
275	>>> import one.alf.io as alfio
276	>>> import ibllib.ephys.spikes as e_spks
277	(*Note, if there is no 'alf' directory, make 'alf' directory from 'ks2' output directory):
278	>>> e_spks.ks2_to_alf(path_to_ks_out, path_to_alf_out)
279	# Get a spikes bunch, units bunch, and filter the units.
280	>>> spks_b = alfio.load_object(path_to_alf_out, 'spikes')
281	>>> units_b = processing.get_units_bunch(spks_b, ['times', 'amps', 'clusters'])
282	>>> T = spks_b['times'][-1] - spks_b['times'][0]
283	>>> filtered_units = processing.filter_units(units_b, T)
284
285	2) Filter units with no minimum amplitude, a minimum firing rate of 1 Hz, and a max false
286	positive rate of 0.2, given a refractory period of 2 ms.
287	>>> filtered_units = processing.filter_units(units_b, T, min_amp=0, min_fr=1)
288
289	TODO: `units_b` input arg could eventually be replaced by `clstrs_b` if the required metrics
290	are in `clstrs_b['metrics']`
291	"""
292
293	# Set params
UNCOV 294	params = {'min_amp': 50e-6, 'min_fr': 0.5, 'max_fpr': 0.2, 'rp': 0.002} # defaults	×
UNCOV 295	params.update(kwargs) # update from **kwargs	×
296
297	# Iteratively filter the units for each filter param #
298	# -------------------------------------------------- #
UNCOV 299	units = np.asarray(list(units_b.amps.keys()))	×
300	# Remove empty clusters
301	empty_cl = np.where([len(units_b.amps[unit]) == 0 for unit in units])[0]	×
UNCOV 302	filt_units = np.delete(units, empty_cl)	×
303	for param in params.keys():	×
304	if param == 'min_amp': # return units above with amp > `'min_amp'`	×
305	mean_amps = np.asarray([np.mean(units_b.amps[unit]) for unit in filt_units])	×
306	filt_idxs = np.where(mean_amps > params['min_amp'])[0]	×
307	filt_units = filt_units[filt_idxs]	×
308	elif param == 'min_fr': # return units with fr > `'min_fr'`	×
NEW 309	fr = np.asarray([len(units_b.amps[unit]) / (units_b.times[unit][-1] - units_b.times[unit][0]) for unit in filt_units])	×
310	filt_idxs = np.where(fr > params['min_fr'])[0]	×
311	filt_units = filt_units[filt_idxs]	×
312	elif param == 'max_fpr': # return units with fpr < `'max_fpr'`	×
313	fpr = np.zeros_like(filt_units, dtype='float')	×
314	for i, unit in enumerate(filt_units):	×
315	n_spks = len(units_b.amps[unit])	×
316	n_isi_viol = len(np.where(np.diff(units_b.times[unit]) < params['rp'])[0])	×
317	# fpr is min of roots of solved quadratic equation (Hill, et al. 2011).
318	c = (t * n_isi_viol) / (2 * params['rp'] * n_spks**2) # 3rd term in quadratic	×
319	fpr[i] = np.min(np.abs(np.roots([-1, 1, c]))) # solve quadratic	×
320	filt_idxs = np.where(fpr < params['max_fpr'])[0]	×
321	filt_units = filt_units[filt_idxs]	×
322	return filt_units.astype(int)	×

int-brain-lab / ibllib / 29037613829

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous