From fa5b8595c17f8a2d8f4a949ecf02be34c8c3927f Mon Sep 17 00:00:00 2001 From: Alessio Buccino Date: Wed, 21 Oct 2026 19:11:00 +0200 Subject: [PATCH 01/13] Fix: do not include extra_metadata in zarr provenance dump (#4885) --- src/spikeinterface/core/base.py | 1 + src/spikeinterface/core/zarrextractors.py | 4 +++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/spikeinterface/core/base.py b/src/spikeinterface/core/base.py index 7fa7fb24c4..1e525885e3 100644 --- a/src/spikeinterface/core/base.py +++ b/src/spikeinterface/core/base.py @@ -556,6 +556,7 @@ def to_dict( to_dict_kwargs = dict( include_annotations=include_annotations, include_properties=include_properties, + include_extra_metadata=include_extra_metadata, # make_paths_relative() will make the recursivity later: relative_to=None, recursive=recursive, diff --git a/src/spikeinterface/core/zarrextractors.py b/src/spikeinterface/core/zarrextractors.py index d2a604e38b..128cb297da 100644 --- a/src/spikeinterface/core/zarrextractors.py +++ b/src/spikeinterface/core/zarrextractors.py @@ -898,7 +898,9 @@ def add_recording_to_zarr_group( zarr_kwargs, job_kwargs = split_job_kwargs(kwargs) if recording.check_serializability("json"): - zarr_group.attrs["provenance"] = check_json(recording.to_dict(recursive=True, relative_to=relative_to)) + zarr_group.attrs["provenance"] = check_json( + recording.to_dict(recursive=True, include_extra_metadata=False, relative_to=relative_to) + ) else: zarr_group.attrs["provenance"] = None From 3d287c71859315cec9df95efe59f0875ac9f1523 Mon Sep 17 00:00:00 2001 From: Charlie Windolf Date: Fri, 23 Oct 2026 10:55:22 -0400 Subject: [PATCH 02/13] dredge fixes (#4889) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- doc/modules/sortingcomponents.rst | 2 +- src/spikeinterface/preprocessing/motion.py | 2 +- .../sortingcomponents/motion/dredge.py | 348 +++++++----------- .../motion/motion_estimation.py | 4 +- .../sortingcomponents/motion/motion_utils.py | 9 +- 5 files changed, 136 insertions(+), 229 deletions(-) diff --git a/doc/modules/sortingcomponents.rst b/doc/modules/sortingcomponents.rst index 9cab4d8a08..6612d84948 100644 --- a/doc/modules/sortingcomponents.rst +++ b/doc/modules/sortingcomponents.rst @@ -298,7 +298,7 @@ Here is an example with non-rigid motion estimation: method="dredge_ap", rigid=False, win_shape="gaussian", - win_step_um=200.0, + win_step_um=400.0, win_scale_um=300.0, win_margin_um=None, bin_um=1.0, diff --git a/src/spikeinterface/preprocessing/motion.py b/src/spikeinterface/preprocessing/motion.py index bd45d3832a..388605c166 100644 --- a/src/spikeinterface/preprocessing/motion.py +++ b/src/spikeinterface/preprocessing/motion.py @@ -38,7 +38,7 @@ win_step_um=400.0, win_scale_um=400.0, win_margin_um=None, - batching_mode="offline", + batching_mode="full", ), "interpolate_motion_kwargs": dict( border_mode="force_extrapolate", spatial_interpolation_method="kriging", sigma_um=20.0, p=2 diff --git a/src/spikeinterface/sortingcomponents/motion/dredge.py b/src/spikeinterface/sortingcomponents/motion/dredge.py index 301deb752e..6d9cc1afdf 100644 --- a/src/spikeinterface/sortingcomponents/motion/dredge.py +++ b/src/spikeinterface/sortingcomponents/motion/dredge.py @@ -23,16 +23,16 @@ """ import gc -import warnings import importlib.util +import warnings import numpy as np from tqdm.auto import trange from spikeinterface.core import BaseRecording from spikeinterface.core.motion import Motion + from .motion_utils import ( - get_spatial_bin_edges, get_spatial_windows, get_window_domains, make_2d_motion_histogram, @@ -91,7 +91,7 @@ class DredgeApRegistration: What torch device to run on? E.g., "cpu" or "cuda" or "cuda:1". batching_mode : Literal["full", "online"], default: "full" Which batching_mode to run when computing displacement vector. - chunk_len_s=None, + chunk_len_s : float, default: 3000.0 If `batching_mode` is "online", how long should each batch be, in seconds. In general, longer is better but will use more RAM. """ @@ -113,10 +113,7 @@ def run( extra, **method_kwargs, ): - torch_spec = importlib.util.find_spec("torch") - if torch_spec is not None: - import torch - else: + if importlib.util.find_spec("torch") is None: raise ImportError("The DREDGE method requires torch: pip install torch") outs = dredge_ap( @@ -142,7 +139,6 @@ def run( return motion -# @TODO : Charlie I started very small refactoring, I let you continue def dredge_ap( recording, peaks, @@ -163,11 +159,8 @@ def dredge_ap( do_window_weights=True, weights_threshold_low=0.2, weights_threshold_high=0.2, - mincorr_percentile=None, - mincorr_percentile_nneighbs=None, # raster arguments - amp_scale_fn=None, ## @Charlie this one is not used anymore - post_transform=np.log1p, ###@this one is directly transimited to weight_correlation_matrix() and so get_wieiith() + post_transform=np.log1p, histogram_depth_smooth_um=1, histogram_time_smooth_s=1, avg_in_bin=False, @@ -179,7 +172,7 @@ def dredge_ap( progress_bar=True, extra_outputs=False, precomputed_D_C_maxdisp=None, - chunk_len_s=None, + chunk_len_s=3000.0, batching_mode="full", ): """Estimate motion from spikes @@ -241,39 +234,9 @@ def dredge_ap( if time_horizon_s: xcorr_kw["max_dt_bins"] = np.ceil(time_horizon_s / bin_s) - # TODO @charlie I think this is a bad to have the dict which is transported to every function - # this should be used only in histogram function but not in weight_correlation_matrix() - # only important kwargs should be explicitly reported - # raster_kw = dict( - # amp_scale_fn=amp_scale_fn, - # post_transform=post_transform, - # histogram_depth_smooth_um=histogram_depth_smooth_um, - # histogram_time_smooth_s=histogram_time_smooth_s, - # bin_s=bin_s, - # bin_um=bin_um, - # avg_in_bin=avg_in_bin, - # return_counts=count_masked_correlation, - # count_bins=count_bins, - # count_bin_min=count_bin_min, - # ) - # this will store return values other than the MotionEstimate extra = {} - # TODO charlie I switch this to make_2d_motion_histogram - # but we need to add all options from the original spike_raster() - # but I think this is OK - # raster_res = spike_raster( - # amps, - # depths_um, - # times_s, - # **raster_kw, - # ) - # if count_masked_correlation: - # raster, spatial_bin_edges_um, time_bin_edges_s, counts = raster_res - # else: - # raster, spatial_bin_edges_um, time_bin_edges_s = raster_res - motion_histogram, time_bin_edges_s, spatial_bin_edges_um = make_2d_motion_histogram( recording, peaks, @@ -286,25 +249,18 @@ def dredge_ap( hist_margin_um=0.0, # @charlie maybe we should expose this and set +20. for instance spatial_bin_edges=None, depth_smooth_um=histogram_depth_smooth_um, - time_smooth_s=histogram_time_smooth_s, + time_smooth_s=None, ) raster = motion_histogram.T - # TODO charlie : put the log for hitstogram + # as in dredge: spatial smoothing, then post_transform, then temporal smoothing + if post_transform is not None: + raster = post_transform(raster) + if histogram_time_smooth_s: + from scipy.ndimage import gaussian_filter1d - # TODO @charlie you should check that we are doing the same thing - # windows, window_centers = get_spatial_windows( - # np.c_[np.zeros_like(spatial_bin_edges_um), spatial_bin_edges_um], - # win_step_um, - # win_scale_um, - # spatial_bin_edges=spatial_bin_edges_um, - # margin_um=-win_scale_um / 2 if win_margin_um is None else win_margin_um, - # win_shape=win_shape, - # zero_threshold=1e-5, - # rigid=rigid, - # ) + raster = gaussian_filter1d(raster, histogram_time_smooth_s / bin_s, axis=1, mode="constant") - dim = ["x", "y", "z"].index(direction) contact_depths = recording.get_channel_locations()[:, dim] spatial_bin_centers = 0.5 * (spatial_bin_edges_um[1:] + spatial_bin_edges_um[:-1]) @@ -319,16 +275,13 @@ def dredge_ap( zero_threshold=1e-5, ) - # TODO charlie : the count has disapeared - # if extra_outputs and count_masked_correlation: - # extra["counts"] = counts - weights_kw = dict( mincorr=mincorr, time_horizon_s=time_horizon_s, do_window_weights=do_window_weights, weights_threshold_low=weights_threshold_low, weights_threshold_high=weights_threshold_high, + post_transform=post_transform, ) full_xcorr_kw = dict( @@ -340,16 +293,15 @@ def dredge_ap( **xcorr_kw, ) - if batching_mode == "full": + threshold_kw = dict(bin_s=bin_s) - threshold_kw = dict(bin_s=bin_s) + if batching_mode == "full": displacement, extra = compute_displacement_full( raster, windows, spatial_bin_edges_um, win_scale_um, - mincorr_percentile, extra, extra_outputs, thomas_kw, @@ -357,46 +309,36 @@ def dredge_ap( full_xcorr_kw, threshold_kw, precomputed_D_C_maxdisp, - post_transform, + progress_bar=progress_bar, ) elif batching_mode == "online": # T_total is number of bin_s in recording T_total = raster.shape[1] - if chunk_len_s is None: - chunk_len_s = max(time_horizon_s, int(np.floor(1000 / bin_s))) - - threshold_kw = dict( - mincorr_percentile_nneighbs=mincorr_percentile_nneighbs, - in_place=True, - soft=False, - time_horizon_s=time_horizon_s, - bin_s=bin_s, - ) + T_chunk = max(1, int(np.ceil(chunk_len_s / bin_s))) displacement, extra = compute_displacement_online( raster, windows, T_total, - chunk_len_s, + T_chunk, spatial_bin_edges_um, win_scale_um, - mincorr_percentile, extra, extra_outputs, thomas_kw, weights_kw, full_xcorr_kw, threshold_kw, + progress_bar=progress_bar, ) else: - raise ValueError(f"No `batching_mode` called {batching_mode}. Available modes are ['batch', 'offline'].") + raise ValueError(f"No `batching_mode` called {batching_mode}. Available modes are ['full', 'online'].") if extra_outputs: extra["windows"] = windows extra["window_centers"] = window_centers - extra["max_disp_um"] = max_disp_um time_bin_centers = 0.5 * (time_bin_edges_s[1:] + time_bin_edges_s[:-1]) @@ -461,10 +403,6 @@ class DredgeLfpRegistration: mincorr : float in [0,1] Minimum correlation between pairs of frames such that they will be included in the optimization of the displacement estimates. - mincorr_percentile, mincorr_percentile_nneighbs - If mincorr_percentile is set to a number in [0, 100], then mincorr will be replaced - by this percentile of the correlations of neighbors within mincorr_percentile_nneighbs - time bins of each other. device : string or torch.device Controls torch device """ @@ -523,8 +461,6 @@ def dredge_online_lfp( time_horizon_s=None, # weighting arguments mincorr=0.8, - mincorr_percentile=None, - mincorr_percentile_nneighbs=20, soft=False, # low-level arguments thomas_kw=None, @@ -593,6 +529,7 @@ def dredge_online_lfp( weights_kw = dict( mincorr=mincorr, time_horizon_s=time_horizon_s, + do_window_weights=False, ) full_xcorr_kw = dict( @@ -605,10 +542,7 @@ def dredge_online_lfp( ) threshold_kw = dict( - mincorr_percentile_nneighbs=mincorr_percentile_nneighbs, - in_place=True, soft=soft, - time_horizon_s=time_horizon_s, bin_s=1 / lfp_recording.sampling_frequency, # only relevant for time_horizon_s ) @@ -619,13 +553,13 @@ def dredge_online_lfp( T_chunk, contact_depths, win_scale_um, - mincorr_percentile, extra, extra_outputs, thomas_kw, weights_kw, full_xcorr_kw, threshold_kw, + progress_bar=progress_bar, ) motion = Motion([P_online.T], [lfp_recording.get_times(0)], window_centers, direction=direction) @@ -643,19 +577,34 @@ def compute_displacement_online( T_chunk, spatial_bin_edges_um, win_scale_um, - mincorr_percentile, extra, extra_outputs, thomas_kw, weights_kw, full_xcorr_kw, threshold_kw, + progress_bar=False, ): B = len(windows) + def get_Us(Ds, Cs, raster, raster_b=None, t_offset_bins=0): + Us, _ = weight_correlation_matrix( + Ds, + Cs, + windows, + raster, + spatial_bin_edges_um, + threshold_kw, + raster_b=raster_b, + t_offset_bins=t_offset_bins, + in_place=not extra_outputs, + **weights_kw, + ) + return Us + # -- allocate output and initialize first chunk - P_online = np.empty((B, T_total), dtype=np.float32) + P_online = np.empty((B, T_total)) # below, t0 is start of prev chunk, t1 start of cur chunk, t2 end of cur t0, t1 = 0, T_chunk @@ -671,13 +620,7 @@ def compute_displacement_online( win_scale_um, **full_xcorr_kw, ) - mincorr = weights_kw["mincorr"] - Ss0, mincorr0 = threshold_correlation_matrix( - Cs0, - mincorr=mincorr, - mincorr_percentile=mincorr_percentile, - **threshold_kw, - ) + Ss0 = get_Us(Ds0, Cs0, traces0.T) if extra_outputs: extra["D"] = [Ds0] extra["C"] = [Cs0] @@ -685,20 +628,18 @@ def compute_displacement_online( extra["D01"] = [] extra["C01"] = [] extra["S01"] = [] - extra["mincorrs"] = [mincorr0] extra["max_disp_um"] = max_disp_um P_online[:, t0:t1], _ = thomas_solve(Ds0, Ss0, **thomas_kw) # -- loop through chunks chunk_starts = range(T_chunk, T_total, T_chunk) - progress_bar = full_xcorr_kw["progress_bar"] if progress_bar: chunk_starts = trange( T_chunk, T_total, T_chunk, - desc=f"Online chunks [{T_chunk}s each]", + desc="Online chunks", ) for t1 in chunk_starts: t2 = min(T_total, t1 + T_chunk) @@ -716,21 +657,16 @@ def compute_displacement_online( spatial_bin_edges_um, win_scale_um, raster_b=traces0.T, + t_offset_bins=T_chunk, **full_xcorr_kw, ) # cross-correlation in current chunk Ds1, Cs1, _ = xcorr_windows(traces1.T, windows, spatial_bin_edges_um, win_scale_um, **full_xcorr_kw) - Ss1, mincorr1 = threshold_correlation_matrix( - Cs1, - mincorr_percentile=mincorr_percentile, - mincorr=mincorr, - **threshold_kw, - ) - Ss10, _ = threshold_correlation_matrix(Cs10, mincorr=mincorr1, t_offset_bins=T_chunk, **threshold_kw) + Ss1 = get_Us(Ds1, Cs1, traces1.T) + Ss10 = get_Us(Ds10, Cs10, traces1.T, raster_b=traces0.T, t_offset_bins=T_chunk) if extra_outputs: - extra["mincorrs"].append(mincorr1) extra["D"].append(Ds1) extra["C"].append(Cs1) extra["S"].append(Ss1) @@ -762,7 +698,6 @@ def compute_displacement_full( windows, spatial_bin_edges_um, win_scale_um, - mincorr_percentile, extra, extra_outputs, thomas_kw, @@ -770,7 +705,7 @@ def compute_displacement_full( full_xcorr_kw, threshold_kw, precomputed_D_C_maxdisp, - post_transform, + progress_bar=False, ): # cross-correlate to get D and C @@ -780,15 +715,12 @@ def compute_displacement_full( windows, spatial_bin_edges_um, win_scale_um, - # TODO charlie : put back the count for the mask - # masks=(counts > 0) if count_masked_correlation else None, - **full_xcorr_kw, + **dict(full_xcorr_kw, progress_bar=progress_bar), ) else: Ds, Cs, max_disp_um = precomputed_D_C_maxdisp # turn Cs into weights - progress_bar = full_xcorr_kw["progress_bar"] Us, wextra = weight_correlation_matrix( Ds, Cs, @@ -796,13 +728,7 @@ def compute_displacement_full( raster, spatial_bin_edges_um, threshold_kw, - # raster_kw, #@charlie this is removed - post_transform=post_transform, # @charlie this isnew - lambda_t=thomas_kw.get("lambda_t", DEFAULT_LAMBDA_T), - eps=thomas_kw.get("eps", DEFAULT_EPS), - progress_bar=progress_bar, in_place=not extra_outputs, - mincorr_percentile=mincorr_percentile, **weights_kw, ) extra.update({k: wextra[k] for k in wextra if k not in ("S", "U")}) @@ -812,9 +738,9 @@ def compute_displacement_full( if extra_outputs: extra["D"] = Ds extra["C"] = Cs + extra["max_disp_um"] = max_disp_um del Cs - # @charlie : is this needed ? gc.collect() # solve for P @@ -838,7 +764,7 @@ def compute_displacement_full( # -- linear algebra, Newton method solver, block tridiagonal (Thomas) solver -def laplacian(n, wink=True, eps=DEFAULT_EPS, lambd=1.0, ridge_mask=None): +def laplacian(n, eps=DEFAULT_EPS, lambd=1.0, ridge_mask=None, wink_start=True, wink_end=True): """Construct a discrete Laplacian operator (plus eps*identity).""" lap = np.zeros((n, n)) if ridge_mask is None: @@ -846,8 +772,9 @@ def laplacian(n, wink=True, eps=DEFAULT_EPS, lambd=1.0, ridge_mask=None): else: diag = lambd + eps * ridge_mask np.fill_diagonal(lap, diag) - if wink: + if wink_start: lap[0, 0] -= 0.5 * lambd + if wink_end: lap[-1, -1] -= 0.5 * lambd # fill diagonal using a for loop for space reasons when this is large for i in range(n - 1): @@ -871,6 +798,7 @@ def neg_hessian_likelihood_term(Ub, Ub_prevcur=None, Ub_curprev=None): if Ub_prevcur is None: np.fill_diagonal(negHUb, diagonal_terms) else: + assert Ub_curprev is not None diagonal_terms += Ub_prevcur.sum(0) + Ub_curprev.sum(1) np.fill_diagonal(negHUb, diagonal_terms) @@ -885,6 +813,7 @@ def newton_rhs( Ub_prevcur=None, Db_curprev=None, Ub_curprev=None, + prior_rhs=None, ): """Newton step right hand side @@ -898,8 +827,12 @@ def newton_rhs( return grad_at_0 # online case + assert Ub_prevcur is not None + assert Db_curprev is not None align_term = (Ub_prevcur.T + Ub_curprev) @ Pb_prev rhs = align_term + grad_at_0 + (Ub_curprev * Db_curprev).sum(1) - (Ub_prevcur * Db_prevcur).sum(0) + if prior_rhs is not None: + rhs += prior_rhs return rhs @@ -913,12 +846,13 @@ def newton_solve_rigid( Ub_prevcur=None, Db_curprev=None, Ub_curprev=None, + prior_rhs=None, ): """Solve the rigid Newton step D is TxT displacement, U is TxT subsampling or soft weights matrix. """ - from scipy.linalg import solve, lstsq + from scipy.linalg import lstsq, solve negHU = neg_hessian_likelihood_term( U, @@ -933,6 +867,7 @@ def newton_solve_rigid( Ub_prevcur=Ub_prevcur, Db_curprev=Db_curprev, Ub_curprev=Ub_curprev, + prior_rhs=prior_rhs, ) try: p = solve(Sigma0inv + negHU, targ, assume_a="pos") @@ -954,7 +889,6 @@ def thomas_solve( Ds_curprev=None, Us_curprev=None, progress_bar=False, - bandwidth=None, ): """Block tridiagonal algorithm, special cased to our setting @@ -977,27 +911,43 @@ def thomas_solve( Ds = np.asarray(Ds, dtype=np.float64) Us = np.asarray(Us, dtype=np.float64) + B, T, T_ = Ds.shape + assert T == T_ + assert Us.shape == Ds.shape + + # spatial prior's diagonal and off-diagonal block laplacian coefficients + lambd_s_diag = np.full(B, float(lambda_s)) + lambd_s_diag[[0, -1]] = lambda_s / 2 + lambd_s_offdiag = -lambda_s / 2 + online = P_prev is not None online_kw_rhs = online_kw_hess = lambda b: {} if online: assert Ds_prevcur is not None assert Us_prevcur is not None + assert Ds_curprev is not None + assert Us_curprev is not None + # priors linking cur first bin to prev last bin + p_last = P_prev[:, -1].astype(np.float64) + prior_rhs = np.zeros((B, T)) + prior_rhs[:, 0] = 0.5 * lambda_t * p_last + if B > 1 and lambda_s: + prior_rhs[:, 0] += 0.5 * lambd_s_diag * p_last + prior_rhs[1:, 0] += 0.5 * lambd_s_offdiag * p_last[:-1] + prior_rhs[:-1, 0] += 0.5 * lambd_s_offdiag * p_last[1:] online_kw_rhs = lambda b: dict( # noqa Pb_prev=P_prev[b].astype(np.float64, copy=False), Db_prevcur=Ds_prevcur[b].astype(np.float64, copy=False), Ub_prevcur=Us_prevcur[b].astype(np.float64, copy=False), Db_curprev=Ds_curprev[b].astype(np.float64, copy=False), Ub_curprev=Us_curprev[b].astype(np.float64, copy=False), + prior_rhs=prior_rhs[b], ) online_kw_hess = lambda b: dict( # noqa Ub_prevcur=Us_prevcur[b].astype(np.float64, copy=False), Ub_curprev=Us_curprev[b].astype(np.float64, copy=False), ) - B, T, T_ = Ds.shape - assert T == T_ - assert Us.shape == Ds.shape - # figure out which temporal bins are included in the problem # these are used to figure out where epsilon can be added # for numerical stability without changing the solution @@ -1005,8 +955,8 @@ def thomas_solve( had_weights[~had_weights.any(axis=1)] = 1 # temporal prior matrix - L_t = [laplacian(T, eps=eps, lambd=lambda_t, ridge_mask=w) for w in had_weights] - extra = dict(L_t=L_t) + L_t = [laplacian(T, eps=eps, lambd=lambda_t, ridge_mask=w, wink_start=not online) for w in had_weights] + extra: dict[str, list[np.ndarray] | np.ndarray] = dict(L_t=L_t) # just solve independent problems when there's no spatial regularization # not that there's much overhead to the backward pass etc but might as well @@ -1018,10 +968,8 @@ def thomas_solve( return P, extra # spatial prior is a sparse, block tridiagonal kronecker product - # the first and last diagonal blocks are - Lambda_s_diagb = laplacian(T, eps=eps, lambd=lambda_s / 2, ridge_mask=had_weights[0]) - # and the off-diagonal blocks are - Lambda_s_offdiag = laplacian(T, eps=0, lambd=-lambda_s / 2) + Lambda_s_diagb = laplacian(T, eps=eps, lambd=lambd_s_diag[0], ridge_mask=had_weights[0], wink_start=not online) + Lambda_s_offdiag = laplacian(T, eps=0, lambd=lambd_s_offdiag, wink_start=not online) # initialize block-LU stuff and forward variable alpha_hat_b = neg_hessian_likelihood_term(Us[0], **online_kw_hess(0)) @@ -1035,10 +983,7 @@ def thomas_solve( # forward pass for b in trange(1, B, desc="Solve") if progress_bar else range(1, B): - if b < B - 1: - Lambda_s_diagb = laplacian(T, eps=eps, lambd=lambda_s, ridge_mask=had_weights[b]) - else: - Lambda_s_diagb = laplacian(T, eps=eps, lambd=lambda_s / 2, ridge_mask=had_weights[b]) + Lambda_s_diagb = laplacian(T, eps=eps, lambd=lambd_s_diag[b], ridge_mask=had_weights[b], wink_start=not online) alpha_hat_b = neg_hessian_likelihood_term(Us[b], **online_kw_hess(b)) alpha_hat_b += L_t[b] @@ -1066,22 +1011,12 @@ def thomas_solve( def threshold_correlation_matrix( Cs, mincorr=0.0, - mincorr_percentile=None, - mincorr_percentile_nneighbs=20, time_horizon_s=0, in_place=False, bin_s=1, - t_offset_bins=None, - T=None, + t_offset_bins=0, soft=True, ): - if mincorr_percentile is not None: - diags = [np.diagonal(Cs, offset=j, axis1=1, axis2=2).ravel() for j in range(1, mincorr_percentile_nneighbs)] - mincorr = np.percentile( - np.concatenate(diags), - mincorr_percentile, - ) - # need abs to avoid -0.0s which cause numerical issues if in_place: Ss = Cs @@ -1095,15 +1030,11 @@ def threshold_correlation_matrix( Ss = np.square((Cs >= mincorr) * Cs) else: Ss = (Cs >= mincorr).astype(Cs.dtype) - if time_horizon_s is not None and time_horizon_s > 0 and T is not None and time_horizon_s < T: - tt0 = bin_s * np.arange(T) - tt1 = tt0 - if t_offset_bins: - tt1 = tt0 + t_offset_bins - dt = tt1[:, None] - tt0[None, :] - mask = (np.abs(dt) <= time_horizon_s).astype(Ss.dtype) - Ss *= mask[None] - return Ss, mincorr + if time_horizon_s: + Ta, Tb = Ss.shape[-2:] + dt_bins = np.arange(Ta)[:, None] + t_offset_bins - np.arange(Tb)[None, :] + Ss *= np.abs(dt_bins) * bin_s <= time_horizon_s + return Ss def xcorr_windows( @@ -1116,10 +1047,10 @@ def xcorr_windows( bin_um=1, max_disp_um=None, max_dt_bins=None, + t_offset_bins=0, progress_bar=True, centered=True, normalized=True, - masks=None, device=None, ): """Main computational function @@ -1155,8 +1086,6 @@ def xcorr_windows( else: T1 = T0 raster_b_ = raster_a_ - if masks is not None: - masks = torch.as_tensor(masks, dtype=torch.float, device=device) # estimate each window's displacement Ds = np.zeros((B, T0, T1), dtype=np.float32) @@ -1183,14 +1112,13 @@ def xcorr_windows( raster_a_[slices[b]], raster_b_[b_low:b_high], weights=window[slices[b]], - masks=None if masks is None else masks[slices[b]], - xmasks=None if masks is None else masks[b_low:b_high], disp=padding, possible_displacement=poss_disp, device=device, centered=centered, normalized=normalized, max_dt_bins=max_dt_bins, + t_offset_bins=t_offset_bins, ) return Ds, Cs, max_disp_um @@ -1200,14 +1128,13 @@ def calc_corr_decent_pair( raster_a, raster_b, weights=None, - masks=None, - xmasks=None, disp=None, batch_size=512, normalized=True, centered=True, possible_displacement=None, max_dt_bins=None, + t_offset_bins=0, device=None, ): """Weighted pairwise cross-correlation @@ -1240,7 +1167,7 @@ def calc_corr_decent_pair( # sensible default: at most half the domain. if disp is None: - disp == D // 2 + disp = D // 2 # range of displacements if D == D_: @@ -1264,17 +1191,15 @@ def calc_corr_decent_pair( C = np.zeros((Ta, Tb), dtype=np.float32) for i in range(0, Ta, batch_size): for j in range(0, Tb, batch_size): - dt_bins = min(abs(i - j), abs(i + batch_size - j), abs(i - j - batch_size)) + dt_low = i + t_offset_bins - (j + batch_size - 1) + dt_high = i + batch_size - 1 + t_offset_bins - j + dt_bins = 0 if dt_low <= 0 <= dt_high else min(abs(dt_low), abs(dt_high)) if max_dt_bins and dt_bins > max_dt_bins: continue - weights_ = weights - if masks is not None: - weights_ = masks.T[i : i + batch_size] * weights corr = normxcorr1d( raster_a[i : i + batch_size], raster_b[j : j + batch_size], - weights=weights_, - xmasks=None if xmasks is None else xmasks.T[j : j + batch_size], + weights=weights, padding=disp, normalized=normalized, centered=centered, @@ -1291,7 +1216,6 @@ def normxcorr1d( template, x, weights=None, - xmasks=None, centered=True, normalized=True, padding="same", @@ -1353,17 +1277,12 @@ def normxcorr1d( x = npx.atleast_2d(x) num_templates, lengtht = template.shape - num_inputs, lengthx = x.shape + _num_inputs, lengthx = x.shape # generalize over weighted / unweighted case device_kw = {} if conv_engine == "numpy" else dict(device=x.device) - if xmasks is None: - onesx = npx.ones((1, 1, lengthx), dtype=x.dtype, **device_kw) - wx = x[:, None, :] - else: - assert xmasks.shape == x.shape - onesx = xmasks[:, None, :] - wx = x[:, None, :] * onesx + onesx = npx.ones((1, 1, lengthx), dtype=x.dtype, **device_kw) + wx = x[:, None, :] no_weights = weights is None if no_weights: weights = npx.ones((1, 1, lengtht), dtype=x.dtype, **device_kw) @@ -1430,33 +1349,25 @@ def normxcorr1d( def get_weights( - Ds, - Ss, - Sigma0inv_t, windows, raster, dbe, - # @charlie raster_kw is removed in favor of post_transform only is this OK ??? - # raster_kw, post_transform=np.log1p, weights_threshold_low=0.0, weights_threshold_high=np.inf, - progress_bar=False, ): """Compute per-time-bin weighting for each nonrigid window""" # determine window-weighted raster "heat" in each nonrigid window # as a function of time assert windows.shape[1] == dbe.size - 1 weights = [] - p_inds = [] - for b in range((len(Ds))): + for b in range(len(windows)): ilow, ihigh = np.flatnonzero(windows[b])[[0, -1]] ihigh += 1 window_sliced = windows[b, ilow:ihigh] weights.append(window_sliced @ raster[ilow:ihigh]) weights_orig = np.array(weights) - # scale_fn = raster_kw["post_transform"] or raster_kw["amp_scale_fn"] scale_fn = post_transform if isinstance(weights_threshold_low, tuple): nspikes_threshold_low, amp_threshold_low = weights_threshold_low @@ -1472,7 +1383,7 @@ def get_weights( weights_thresh[weights_orig < weights_threshold_low] = 0 weights_thresh[weights_orig > weights_threshold_high] = np.inf - return weights, weights_thresh, p_inds + return weights_orig, weights_thresh def weight_correlation_matrix( @@ -1482,22 +1393,21 @@ def weight_correlation_matrix( raster, depth_bin_edges, threshold_kw, - # @charlie raster_kw is remove in favor of post_transform only - # raster_kw, + raster_b=None, + t_offset_bins=0, post_transform=np.log1p, mincorr=0.0, - mincorr_percentile=None, - mincorr_percentile_nneighbs=20, time_horizon_s=None, - lambda_t=DEFAULT_LAMBDA_T, - eps=DEFAULT_EPS, do_window_weights=True, weights_threshold_low=0.0, weights_threshold_high=np.inf, - progress_bar=True, in_place=False, ): - """Transform the correlation matrix into the weights used in optimization.""" + """Transform the correlation matrix into the weights used in optimization. + + If raster_b is supplied, Cs[:, i, j] correlates time bin i of raster with + time bin j of raster_b, as in xcorr_windows(). + """ extra = {} Ds = np.asarray(Ds) @@ -1505,54 +1415,50 @@ def weight_correlation_matrix( if Ds.ndim == 2: Ds = Ds[None] Cs = Cs[None] - B, T, T_ = Ds.shape - assert T == T_ assert Ds.shape == Cs.shape - extra = {} - Ss, mincorr = threshold_correlation_matrix( + Ss = threshold_correlation_matrix( Cs, mincorr=mincorr, - mincorr_percentile=mincorr_percentile, - mincorr_percentile_nneighbs=mincorr_percentile_nneighbs, time_horizon_s=time_horizon_s, - T=T, in_place=in_place, + t_offset_bins=t_offset_bins, **threshold_kw, ) extra["S"] = Ss - extra["mincorr"] = mincorr if not do_window_weights: return Ss, extra # get weights - L_t = lambda_t * laplacian(T, eps=max(1e-5, eps)) - weights_orig, weights_thresh, Pind = get_weights( - Ds, - Ss, - L_t, + weights_orig, weights_thresh = get_weights( windows, raster, depth_bin_edges, - # raster_kw, post_transform=post_transform, weights_threshold_low=weights_threshold_low, weights_threshold_high=weights_threshold_high, - progress_bar=progress_bar, ) + weights_thresh_b = weights_thresh + if raster_b is not None: + _, weights_thresh_b = get_weights( + windows, + raster_b, + depth_bin_edges, + post_transform=post_transform, + weights_threshold_low=weights_threshold_low, + weights_threshold_high=weights_threshold_high, + ) extra["weights_orig"] = weights_orig extra["weights_thresh"] = weights_thresh - extra["Pind"] = Pind # update noise model. we deliberately divide by zero and inf here. Us = Ss if in_place else np.zeros_like(Ss) with np.errstate(divide="ignore"): - # low mem impl of U = abs(1/(1/weights_thresh+1/weights_thresh'+1/S)) + # low mem impl of U = abs(1/(1/weights_thresh+1/weights_thresh_b'+1/S)) np.reciprocal(Ss, out=Us) - invW = 1.0 / weights_thresh - Us += invW[:, :, None] - Us += invW[:, None, :] + Us += 1.0 / weights_thresh[:, :, None] + Us += 1.0 / weights_thresh_b[:, None, :] np.reciprocal(Us, out=Us) # handles possible -0s that cause issues elsewhere np.abs(Us, out=Us) diff --git a/src/spikeinterface/sortingcomponents/motion/motion_estimation.py b/src/spikeinterface/sortingcomponents/motion/motion_estimation.py index b1faad758b..7e48ef146e 100644 --- a/src/spikeinterface/sortingcomponents/motion/motion_estimation.py +++ b/src/spikeinterface/sortingcomponents/motion/motion_estimation.py @@ -19,7 +19,7 @@ def estimate_motion( direction="y", rigid=False, win_shape="gaussian", - win_step_um=200.0, + win_step_um=400.0, win_scale_um=300.0, win_margin_um=None, method="dredge_ap", @@ -64,7 +64,7 @@ def estimate_motion( width controlled by win_scale_um. When win_margin_um is None the margin is automatically set to -win_scale_um/2. See get_spatial_windows. - win_step_um : float, default: 50 + win_step_um : float, default: 400 See win_shape win_scale_um : float, default: 150 See win_shape diff --git a/src/spikeinterface/sortingcomponents/motion/motion_utils.py b/src/spikeinterface/sortingcomponents/motion/motion_utils.py index 2f85a54cd0..c58322d2c6 100644 --- a/src/spikeinterface/sortingcomponents/motion/motion_utils.py +++ b/src/spikeinterface/sortingcomponents/motion/motion_utils.py @@ -75,17 +75,17 @@ def get_spatial_windows( min_ = np.min(contact_depths) - win_margin_um max_ = np.max(contact_depths) + win_margin_um - num_windows = int((max_ - min_) // win_step_um) + num_windows = int((max_ - min_) // win_step_um) + 1 - if num_windows < 1: - num_windows = 1 + if num_windows < 2: warnings.warn( f"You are trying to estimate motion with `non-rigid` on a probe that is too short!" f"get_spatial_windows(): {win_step_um=}/{win_scale_um=}/{win_margin_um=} are too large for the " f"probe size (depth range={np.ptp(contact_depths)}). Switching to rigid motion." ) + return get_spatial_windows(contact_depths, spatial_bin_centers, rigid=True, zero_threshold=zero_threshold) border = ((max_ - min_) % win_step_um) / 2 - window_centers = np.arange(num_windows + 1) * win_step_um + min_ + border + window_centers = np.arange(num_windows) * win_step_um + min_ + border windows = [] for win_center in window_centers: @@ -106,6 +106,7 @@ def get_spatial_windows( windows = np.array(windows) if zero_threshold is not None: + windows /= windows.sum(axis=1, keepdims=True) windows[windows < zero_threshold] = 0 windows /= windows.sum(axis=1, keepdims=True) From 2d26e2fc572926da75997153ebeb5456f8e32775 Mon Sep 17 00:00:00 2001 From: Heberto Mayorquin Date: Wed, 21 Oct 2026 04:23:28 -0600 Subject: [PATCH 03/13] Save `dtype` in `phase_shift` kwargs (#4886) Co-authored-by: Galen Lynch --- src/spikeinterface/preprocessing/phase_shift.py | 7 ++++++- .../preprocessing/tests/test_phase_shift.py | 14 +++++++++++++- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/src/spikeinterface/preprocessing/phase_shift.py b/src/spikeinterface/preprocessing/phase_shift.py index c01a9acf59..b7cdb456fd 100644 --- a/src/spikeinterface/preprocessing/phase_shift.py +++ b/src/spikeinterface/preprocessing/phase_shift.py @@ -65,7 +65,12 @@ def __init__(self, recording, margin_ms=40.0, inter_sample_shift=None, dtype=Non # for dumpability if inter_sample_shift is not None: inter_sample_shift = list(inter_sample_shift) - self._kwargs = dict(recording=recording, margin_ms=float(margin_ms), inter_sample_shift=inter_sample_shift) + self._kwargs = dict( + recording=recording, + margin_ms=float(margin_ms), + inter_sample_shift=inter_sample_shift, + dtype=np.dtype(dtype).str, + ) class PhaseShiftRecordingSegment(BasePreprocessorSegment): diff --git a/src/spikeinterface/preprocessing/tests/test_phase_shift.py b/src/spikeinterface/preprocessing/tests/test_phase_shift.py index 23718f6ac6..5d4b74a567 100644 --- a/src/spikeinterface/preprocessing/tests/test_phase_shift.py +++ b/src/spikeinterface/preprocessing/tests/test_phase_shift.py @@ -1,5 +1,6 @@ import numpy as np -from spikeinterface import NumpyRecording +from spikeinterface import NumpyRecording, load +from spikeinterface.core.testing import check_recordings_equal from spikeinterface.preprocessing import phase_shift @@ -112,5 +113,16 @@ def signal(time): np.testing.assert_allclose(shifted_traces[interior], expected[interior], atol=1e-3) +def test_phase_shift_dtype_survives_reload(): + traces, sampling_frequency, inter_sample_shift = create_shifted_channel() + recording = NumpyRecording([traces.astype("int16")], sampling_frequency) + recording.set_property("inter_sample_shift", inter_sample_shift) + + shifted = phase_shift(recording, dtype="float32") + reloaded = load(shifted.to_dict()) + + check_recordings_equal(shifted, reloaded, return_in_uV=False) + + if __name__ == "__main__": test_phase_shift() From 2c88b4c70885d0ff36beedb2d00865dd79be44d2 Mon Sep 17 00:00:00 2001 From: Alessio Buccino Date: Thu, 22 Oct 2026 10:39:43 +0200 Subject: [PATCH 04/13] Fix: openephys - check sync_times only if not None (#4890) --- src/spikeinterface/extractors/neoextractors/openephys.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/spikeinterface/extractors/neoextractors/openephys.py b/src/spikeinterface/extractors/neoextractors/openephys.py index 1896ed34d0..729e6f15a3 100644 --- a/src/spikeinterface/extractors/neoextractors/openephys.py +++ b/src/spikeinterface/extractors/neoextractors/openephys.py @@ -354,7 +354,7 @@ def __init__( sync_times = np.load(stream_folder / "synchronized_timestamps.npy", mmap_mode="r") else: sync_times = None - if len(sync_times) == self.segments[segment_index].get_num_samples(): + if sync_times is not None and len(sync_times) == self.segments[segment_index].get_num_samples(): self.segments[segment_index]._time_vector = sync_times else: warnings.warn(f"Could not load synchronized timestamps for {stream_name}") From 43c1015966b7b96002a2bdeffab34a279744d9d6 Mon Sep 17 00:00:00 2001 From: Pierre Yger Date: Mon, 26 Oct 2026 15:23:49 +0100 Subject: [PATCH 05/13] Add main_channel_id to Lupin (#4894) Co-authored-by: Chris Halcrow <57948917+chrishalcrow@users.noreply.github.com> Co-authored-by: Alessio Buccino --- src/spikeinterface/sorters/internal/lupin.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/spikeinterface/sorters/internal/lupin.py b/src/spikeinterface/sorters/internal/lupin.py index 6d79bb036d..b0e6d1dbbc 100644 --- a/src/spikeinterface/sorters/internal/lupin.py +++ b/src/spikeinterface/sorters/internal/lupin.py @@ -405,6 +405,9 @@ def _run_from_folder(cls, sorter_output_folder, params, verbose): final_spikes["segment_index"] = spikes["segment_index"] sorting = NumpySorting(final_spikes, sampling_frequency, templates.unit_ids) + main_channel_ids = templates.get_main_channels(peak_sign=params["peak_sign"], outputs="id") + sorting.set_property("main_channel_id", main_channel_ids) + auto_merge = True analyzer_final = None if auto_merge: @@ -425,7 +428,10 @@ def _run_from_folder(cls, sorter_output_folder, params, verbose): debug_folder=None, job_kwargs=job_kwargs, ) - sorting = NumpySorting.from_sorting(analyzer_final.sorting) + # This is an internal merge, we don't need to expose the "is_merged" property to the user + if "is_merged" in analyzer_final.sorting.get_property_keys(): + analyzer_final.sorting.delete_property("is_merged") + sorting = NumpySorting.from_sorting(analyzer_final.sorting, with_metadata=True) if params["save_array"]: sorting_pre_peeler = sorting_pre_peeler.save(folder=sorter_output_folder / "sorting_pre_peeler") From a9524c808517572b4487f23fefdf8b90e0f8b40e Mon Sep 17 00:00:00 2001 From: Alessio Buccino Date: Mon, 26 Oct 2026 15:24:12 +0100 Subject: [PATCH 06/13] Remove legacy location property when loading old recordings (#4896) --- src/spikeinterface/core/baserecordingsnippets.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/spikeinterface/core/baserecordingsnippets.py b/src/spikeinterface/core/baserecordingsnippets.py index 11eb909b59..648074d3a7 100644 --- a/src/spikeinterface/core/baserecordingsnippets.py +++ b/src/spikeinterface/core/baserecordingsnippets.py @@ -321,14 +321,13 @@ def _extra_metadata_copy(self, other): other._probegroup = self._probegroup.copy() def _extra_metadata_from_dict(self, dump_dict): - # load probe and handle backward-compatibility with legacy "contact_vector"/"location" property + # load probegroup if "probegroup" in dump_dict: - # this is for SI>=0.105.0 probegroup = dump_dict["probegroup"] self._probegroup = ProbeGroup.from_dict(probegroup) def _extra_metadata_to_dict(self, dump_dict): - # save probe + # save probegroup if self.has_probe(): probegroup = self.get_probegroup() dump_dict["probegroup"] = probegroup.to_dict() @@ -360,6 +359,7 @@ def _handle_extractor_backward_compatibility(self): locations = self.get_property("location") self.set_dummy_probe_from_locations(locations) check_for_probes_info = True + self.delete_property("location") if check_for_probes_info: for i, probe in enumerate(self._probegroup.probes): From 8765aa8f519b91461707f97e55d677d8a35aa025 Mon Sep 17 00:00:00 2001 From: Heberto Mayorquin Date: Mon, 19 Oct 2026 02:46:05 -0600 Subject: [PATCH 07/13] Make motion-corrected and silenced recordings JSON-serializable (#4877) --- src/spikeinterface/core/motion.py | 7 ++++--- src/spikeinterface/preprocessing/silence_periods.py | 6 +++--- .../motion/motion_interpolation.py | 13 ++++++++----- 3 files changed, 15 insertions(+), 11 deletions(-) diff --git a/src/spikeinterface/core/motion.py b/src/spikeinterface/core/motion.py index 6a99ca29c8..50622b2e7d 100644 --- a/src/spikeinterface/core/motion.py +++ b/src/spikeinterface/core/motion.py @@ -166,10 +166,11 @@ def to_dict(self): @staticmethod def from_dict(d): + # arrays are lists when the dict comes from json return Motion( - d["displacement"], - d["temporal_bins_s"], - d["spatial_bins_um"], + [np.asarray(displacement) for displacement in d["displacement"]], + [np.asarray(temporal_bins_s) for temporal_bins_s in d["temporal_bins_s"]], + np.asarray(d["spatial_bins_um"]), direction=d["direction"], interpolation_method=d["interpolation_method"], ) diff --git a/src/spikeinterface/preprocessing/silence_periods.py b/src/spikeinterface/preprocessing/silence_periods.py index 430348bbb3..912b606174 100644 --- a/src/spikeinterface/preprocessing/silence_periods.py +++ b/src/spikeinterface/preprocessing/silence_periods.py @@ -79,6 +79,9 @@ def __init__( "list_periods is deprecated. Please specify periods, which should be a np.array with " f"dtype {base_period_dtype}" ) + if isinstance(periods, list): + # periods reloaded from json are a list of [segment_index, start_sample_index, end_sample_index] + periods = np.array([tuple(period) for period in periods], dtype=base_period_dtype) if not isinstance(periods, np.ndarray): raise ValueError(f"periods must be a np.array with dtype {base_period_dtype}") @@ -132,9 +135,6 @@ def __init__( ) self.add_recording_segment(rec_segment) - # the base_period_dtype is a structured dtype, which is not json serializable - self._serializability["json"] = False - self._kwargs = dict( recording=recording, periods=periods, diff --git a/src/spikeinterface/sortingcomponents/motion/motion_interpolation.py b/src/spikeinterface/sortingcomponents/motion/motion_interpolation.py index a144940853..0eb63b8b6b 100644 --- a/src/spikeinterface/sortingcomponents/motion/motion_interpolation.py +++ b/src/spikeinterface/sortingcomponents/motion/motion_interpolation.py @@ -1,6 +1,6 @@ import numpy as np from spikeinterface.core.core_tools import define_function_handling_dict_from_class -from spikeinterface.core.motion import ensure_time_bins +from spikeinterface.core.motion import Motion, ensure_time_bins from spikeinterface.preprocessing import get_spatial_interpolation_kernel from spikeinterface.preprocessing.basepreprocessor import BasePreprocessor, BasePreprocessorSegment from spikeinterface.preprocessing.filter import fix_dtype @@ -353,6 +353,10 @@ def __init__( ): # assert recording.get_num_segments() == 1, "correct_motion() is only available for single-segment recordings" + if isinstance(motion, dict): + # a Motion reloaded from json is a dict + motion = Motion.from_dict(motion) + channel_locations = recording.get_channel_locations() assert channel_locations.ndim >= motion.dim, ( f"'direction' {motion.direction} not available. " @@ -419,6 +423,9 @@ def __init__( interpolation_time_bin_centers_s = motion.temporal_bins_s interpolation_time_bin_edges_s = motion.temporal_bin_edges_s else: + if isinstance(interpolation_time_bin_centers_s, list): + # time bins reloaded from json are lists + interpolation_time_bin_centers_s = [np.asarray(bins) for bins in interpolation_time_bin_centers_s] interpolation_time_bin_centers_s, interpolation_time_bin_edges_s = ensure_time_bins( interpolation_time_bin_centers_s, interpolation_time_bin_edges_s ) @@ -453,10 +460,6 @@ def __init__( ) self.add_recording_segment(rec_segment) - # this object is currently not JSON-serializable because the Motion obejct cannot be reloaded properly - # see issue #3313 - self._serializability["json"] = False - self._kwargs = dict( recording=recording, motion=motion, From 1851d54bd3d01db9c4241466d46160b6809daff8 Mon Sep 17 00:00:00 2001 From: Chris Halcrow <57948917+chrishalcrow@users.noreply.github.com> Date: Tue, 27 Oct 2026 11:35:28 +0000 Subject: [PATCH 08/13] Allow for `main_channel_indices` to use rec attributes (#4867) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Alessio Buccino --- src/spikeinterface/core/sortinganalyzer.py | 12 ++++++-- .../core/tests/test_sortinganalyzer.py | 28 +++++++++---------- 2 files changed, 23 insertions(+), 17 deletions(-) diff --git a/src/spikeinterface/core/sortinganalyzer.py b/src/spikeinterface/core/sortinganalyzer.py index b702237445..fe8431dc3b 100644 --- a/src/spikeinterface/core/sortinganalyzer.py +++ b/src/spikeinterface/core/sortinganalyzer.py @@ -817,8 +817,16 @@ def _handle_backward_compatibility_settings_pre_init(cls, settings: dict[str, An def main_channel_indices(self): if self._main_channel_indices is None: sorting_main_channel_ids = self.get_sorting_property("main_channel_id") - if sorting_main_channel_ids is not None and self.has_recording(): - main_channel_indices = self.recording.ids_to_indices(sorting_main_channel_ids) + if sorting_main_channel_ids is not None and (self.has_recording() or self.rec_attributes): + if self.has_recording(): + main_channel_indices = self.recording.ids_to_indices(sorting_main_channel_ids) + elif self.rec_attributes.get("channel_ids") is not None: + channel_ids = np.array(self.rec_attributes.get("channel_ids")) + main_channel_indices = np.array( + [list(channel_ids).index(id) for id in sorting_main_channel_ids], dtype=int + ) + else: + main_channel_indices = None self._main_channel_indices = main_channel_indices else: self._main_channel_indices = self._compute_main_channel_backwards_compatibility() diff --git a/src/spikeinterface/core/tests/test_sortinganalyzer.py b/src/spikeinterface/core/tests/test_sortinganalyzer.py index 527e053184..20293f6fd9 100644 --- a/src/spikeinterface/core/tests/test_sortinganalyzer.py +++ b/src/spikeinterface/core/tests/test_sortinganalyzer.py @@ -1160,10 +1160,10 @@ def test_select_channels_independent(dataset): assert np.all(sa_one.get_channel_locations() == sa_two.get_channel_locations()[shuffle_order]) -def test_main_channel_from_templates_dense_recordingless(tmp_path): - """When a dense analyzer has a `templates` extension but no attached recording, `get_main_channels` - recovers each unit's main channel from the templates (its peak channel) and reports it consistently - both as a positional index and as a channel id. +def test_main_channel_from_templates_dense_no_sorting_main_channel_ids(tmp_path): + """When a dense analyzer has a `templates` extension but the sorting has no main channel + ids, `get_main_channels, recovers each unit's main channel from the templates (its peak + channel) and reports it consistently both as a positional index and as a channel id. """ recording, sorting = generate_ground_truth_recording( num_channels=8, @@ -1189,18 +1189,16 @@ def test_main_channel_from_templates_dense_recordingless(tmp_path): expected_main_channel_ids = [unit0_main_channel_id, unit1_main_channel_id] folder = tmp_path / "analyzer_recordingless" - analyzer = create_sorting_analyzer(sorting, recording, format="binary_folder", folder=folder, sparse=False) + analyzer = create_sorting_analyzer(sorting, recording, format="memory", folder=folder, sparse=False) analyzer.compute(["random_spikes", "templates"]) - # simulate the recording being unavailable at load time by removing its serialized copy on disk - for file_ext in ("json", "pickle"): - recording_file = folder / f"recording.{file_ext}" - if recording_file.exists(): - recording_file.unlink() + analyzer._recording = None + analyzer.sorting._properties.pop("main_channel_id") + + analyzer.save_as(format="binary_folder", folder=folder) reloaded_analyzer = load_sorting_analyzer(folder) - assert not reloaded_analyzer.has_recording() - assert reloaded_analyzer.has_extension("templates") + assert "main_channel_id" not in reloaded_analyzer.get_sorting_property_keys() # craft a clear negative peak on each unit's main channel templates = reloaded_analyzer.get_extension("templates") @@ -1217,8 +1215,8 @@ def test_main_channel_from_templates_dense_recordingless(tmp_path): assert np.array_equal(recovered_main_channel_ids, expected_main_channel_ids) -def test_main_channel_from_templates_sparse_recordingless(): - """When a sparse analyzer has a `templates` extension but no attached recording, `get_main_channels` +def test_main_channel_from_templates_sparse_no_sorting_main_channel_ids(): + """When a sparse analyzer has a `templates` extension but sorting has no main channel ids, `get_main_channels` reports each unit's main channel as its template peak channel, not the geometric centroid of the unit's sparse channel set. """ @@ -1272,8 +1270,8 @@ def test_main_channel_from_templates_sparse_recordingless(): templates.data["average"] = average # recordingless, and force the backwards-compatibility recompute - analyzer._recording = None analyzer._main_channel_indices = None + analyzer.sorting._properties.pop("main_channel_id") # must be the template peak channels, not the geometric centroids ("c", "f" / indices [2, 5]) recovered_main_channel_indices = np.asarray(analyzer.get_main_channels(outputs="index")) From 1a7fba3991107f7b5e31e8a43b1678cc9d2745e4 Mon Sep 17 00:00:00 2001 From: Alessio Buccino Date: Tue, 27 Oct 2026 12:42:56 +0100 Subject: [PATCH 09/13] Prepare 0.105.4 release --- doc/releases/0.105.4.rst | 37 +++++++++++++++++++++++++++++++++++++ doc/whatisnew.rst | 3 ++- pyproject.toml | 2 +- 3 files changed, 40 insertions(+), 2 deletions(-) create mode 100644 doc/releases/0.105.4.rst diff --git a/doc/releases/0.105.4.rst b/doc/releases/0.105.4.rst new file mode 100644 index 0000000000..00d16a4a9f --- /dev/null +++ b/doc/releases/0.105.4.rst @@ -0,0 +1,37 @@ +.. _release0.105.4: + +SpikeInterface 0.105.4 release notes +------------------------------------ + +October 27th 2026 + +Minor release with bug fixes + +core: + +* Allow for main_channel_indices to use rec attributes (#4867) +* Remove legacy location property when loading old recordings (#4896) +* Fix: do not include extra_metadata in zarr provenance dump (#4885) + +preprocessing: + +* Feat: expose periods property in silence_periods-extractor (#4897) +* Dredge fixes (#4889) +* Save dtype in phase_shift kwargs (#4886) +* Make motion-corrected and silenced recordings JSON-serializable (#4877) + +extractors: + +* Fix: openephys - check sync_times only if not None (#4890) + +sorters: + +* Add main_channel_id to Lupin (#4894) + + +Contributors: + +* @alejoe91 +* @cwindolf +* @h-mayorquin +* @yger diff --git a/doc/whatisnew.rst b/doc/whatisnew.rst index 9b2f204339..cd841b3a6c 100644 --- a/doc/whatisnew.rst +++ b/doc/whatisnew.rst @@ -8,6 +8,7 @@ Release notes .. toctree:: :maxdepth: 1 + releases/0.105.4.rst releases/0.105.3.rst releases/0.105.2.rst releases/0.105.1.rst @@ -66,7 +67,7 @@ Release notes releases/0.9.1.rst -Versions 0.105.1/3 +Versions 0.105.1/4 ================== * Minor releases with bug fixes and performance improvements diff --git a/pyproject.toml b/pyproject.toml index 66bf4280b7..1741e6b873 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "spikeinterface" -version = "0.105.3" +version = "0.105.4" authors = [ { name="Alessio Buccino", email="alessiop.buccino@gmail.com" }, { name="Samuel Garcia", email="sam.garcia.die@gmail.com" }, From c75aa562a0d1833669498e5e32ac20e3357ffba5 Mon Sep 17 00:00:00 2001 From: Alessio Buccino Date: Tue, 27 Oct 2026 12:48:00 +0100 Subject: [PATCH 10/13] feat: expose periods property in silence_periods (#4897) --- src/spikeinterface/preprocessing/silence_periods.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/spikeinterface/preprocessing/silence_periods.py b/src/spikeinterface/preprocessing/silence_periods.py index 912b606174..60ac188b98 100644 --- a/src/spikeinterface/preprocessing/silence_periods.py +++ b/src/spikeinterface/preprocessing/silence_periods.py @@ -135,6 +135,8 @@ def __init__( ) self.add_recording_segment(rec_segment) + self._periods = periods + self._kwargs = dict( recording=recording, periods=periods, @@ -144,6 +146,10 @@ def __init__( apodization_samples=apodization_samples, ) + @property + def periods(self): + return self._periods + def _all_period_list_to_periods_vec(list_periods, num_seg): if num_seg == 1: From 26f84acec913956c65262b6eb1780fb639ccc9d6 Mon Sep 17 00:00:00 2001 From: Alessio Buccino Date: Tue, 27 Oct 2026 12:50:51 +0100 Subject: [PATCH 11/13] cosmetics --- doc/releases/0.105.4.rst | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/doc/releases/0.105.4.rst b/doc/releases/0.105.4.rst index 00d16a4a9f..96cdc41866 100644 --- a/doc/releases/0.105.4.rst +++ b/doc/releases/0.105.4.rst @@ -9,24 +9,24 @@ Minor release with bug fixes core: -* Allow for main_channel_indices to use rec attributes (#4867) -* Remove legacy location property when loading old recordings (#4896) -* Fix: do not include extra_metadata in zarr provenance dump (#4885) +* Allow for ``main_channel_indices`` to use ``rec_attributes`` (#4867) +* Remove legacy ``location`` property when loading old recordings (#4896) +* Fix: do not include ``extra_metadata`` in zarr provenance dump (#4885) preprocessing: -* Feat: expose periods property in silence_periods-extractor (#4897) +* Feat: expose ``periods`` property in silence_periods (#4897) * Dredge fixes (#4889) -* Save dtype in phase_shift kwargs (#4886) +* Save ``dtype`` in ``phase_shift`` kwargs (#4886) * Make motion-corrected and silenced recordings JSON-serializable (#4877) extractors: -* Fix: openephys - check sync_times only if not None (#4890) +* Fix: openephys - check ``sync_times`` only if not None (#4890) sorters: -* Add main_channel_id to Lupin (#4894) +* Add ``main_channel_id`` to Lupin (#4894) Contributors: From 4086a6957b6b152993cea832c7dd68c53862d131 Mon Sep 17 00:00:00 2001 From: Alessio Buccino Date: Tue, 27 Oct 2026 14:58:18 +0100 Subject: [PATCH 12/13] Apply suggestion from @alejoe91 --- doc/releases/0.105.4.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/releases/0.105.4.rst b/doc/releases/0.105.4.rst index 96cdc41866..d62d79970c 100644 --- a/doc/releases/0.105.4.rst +++ b/doc/releases/0.105.4.rst @@ -32,6 +32,7 @@ sorters: Contributors: * @alejoe91 +* @chrishalcrow * @cwindolf * @h-mayorquin * @yger From d9b30348f89dda66d000b9be715a071d09e9f488 Mon Sep 17 00:00:00 2001 From: Alessio Buccino Date: Tue, 27 Oct 2026 14:59:44 +0100 Subject: [PATCH 13/13] Apply suggestion from @alejoe91 --- doc/releases/0.105.4.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/releases/0.105.4.rst b/doc/releases/0.105.4.rst index d62d79970c..f19c3cf9cc 100644 --- a/doc/releases/0.105.4.rst +++ b/doc/releases/0.105.4.rst @@ -34,5 +34,6 @@ Contributors: * @alejoe91 * @chrishalcrow * @cwindolf +* @galenlynch * @h-mayorquin * @yger