SpikeInterface / spikeinterface

A Python-based module for creating flexible and robust spike sorting pipelines.
https://spikeinterface.readthedocs.io
MIT License
501 stars 186 forks source link

estimate_motion: Only zero batch or zero channel inputs are supported, but got input shape: [1, 1, 1, 0] #2603

Closed borrepp closed 2 months ago

borrepp commented 6 months ago

Hello, I have a 32-channel recording with the following specs:

NwbRecordingExtractor: 32 channels - 30.0kHz - 1 segments - 44,771,040 samples 
                       1,492.37s (24.87 minutes) - float64 dtype - 10.67 GiB

job_kwargs = dict(chunk_duration="1s", n_jobs=28, progress_bar=True)

I have saved the filtered recording and now I'm trying to run "correc_motion" with the preset="nonrigid_accurate", but I got the following error:

I'm running on Windows 11, with spikeInterface version = 0.100.0

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[7], [line 6](vscode-notebook-cell:?execution_count=7&line=6)
      [4](vscode-notebook-cell:?execution_count=7&line=4)     os.mkdir(motion_folder)
      [5](vscode-notebook-cell:?execution_count=7&line=5) job_kwargs = dict(chunk_duration="1s", n_jobs=n_jobs, progress_bar=True)
----> [6](vscode-notebook-cell:?execution_count=7&line=6) rec_corrected = si.correct_motion(recording=rec_filter_file, preset="nonrigid_accurate", folder=motion_folder, 
      [7](vscode-notebook-cell:?execution_count=7&line=7)                                   **job_kwargs)
      [9](vscode-notebook-cell:?execution_count=7&line=9) #si.plot_traces({'filtered':rec_filter,  "corrected": rec_corrected}, time_range=[0, 100], backend="ipywidgets")

File [m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\preprocessing\motion.py:394](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/preprocessing/motion.py:394), in correct_motion(recording, preset, folder, output_motion_info, detect_kwargs, select_kwargs, localize_peaks_kwargs, estimate_motion_kwargs, interpolate_motion_kwargs, **job_kwargs)
    [391](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/preprocessing/motion.py:391)     np.save(folder / "peak_locations.npy", peak_locations)
    [393](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/preprocessing/motion.py:393) t0 = time.perf_counter()
--> [394](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/preprocessing/motion.py:394) motion, temporal_bins, spatial_bins = estimate_motion(recording, peaks, peak_locations, **estimate_motion_kwargs)
    [395](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/preprocessing/motion.py:395) t1 = time.perf_counter()
    [396](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/preprocessing/motion.py:396) run_times["estimate_motion"] = t1 - t0

File [m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\sortingcomponents\motion_estimation.py:152](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:152), in estimate_motion(recording, peaks, peak_locations, direction, bin_duration_s, bin_um, margin_um, rigid, win_shape, win_step_um, win_sigma_um, post_clean, speed_threshold, sigma_smooth_s, method, output_extra_check, progress_bar, upsample_to_histogram_bin, verbose, **method_kwargs)
    [150](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:150) # run method
    [151](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:151) method_class = estimate_motion_methods[method]
--> [152](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:152) motion, temporal_bins = method_class.run(
    [153](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:153)     recording,
    [154](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:154)     peaks,
    [155](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:155)     peak_locations,
    [156](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:156)     direction,
    [157](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:157)     bin_duration_s,
    [158](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:158)     bin_um,
    [159](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:159)     spatial_bin_edges,
    [160](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:160)     non_rigid_windows,
    [161](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:161)     verbose,
    [162](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:162)     progress_bar,
    [163](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:163)     extra_check,
    [164](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:164)     **method_kwargs,
    [165](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:165) )
    [167](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:167) # replace nan by zeros
    [168](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:168) motion[np.isnan(motion)] = 0

File [m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\sortingcomponents\motion_estimation.py:363](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:363), in DecentralizedRegistration.run(cls, recording, peaks, peak_locations, direction, bin_duration_s, bin_um, spatial_bin_edges, non_rigid_windows, verbose, progress_bar, extra_check, histogram_depth_smooth_um, histogram_time_smooth_s, pairwise_displacement_method, max_displacement_um, weight_scale, error_sigma, conv_engine, torch_device, batch_size, corr_threshold, time_horizon_s, convergence_method, soft_weights, normalized_xcorr, centered_xcorr, temporal_prior, spatial_prior, force_spatial_median_continuity, reference_displacement, reference_displacement_time_s, robust_regression_sigma, lsqr_robust_n_iter, weight_with_amplitude)
    [360](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:360) if verbose:
    [361](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:361)     print(f"Computing pairwise displacement: {i + 1} / {len(non_rigid_windows)}")
--> [363](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:363) pairwise_displacement, pairwise_displacement_weight = compute_pairwise_displacement(
    [364](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:364)     motion_histogram[:, window_slice],
    [365](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:365)     bin_um,
    [366](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:366)     window=win[window_slice],
    [367](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:367)     method=pairwise_displacement_method,
    [368](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:368)     weight_scale=weight_scale,
    [369](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:369)     error_sigma=error_sigma,
    [370](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:370)     conv_engine=conv_engine,
    [371](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:371)     torch_device=torch_device,
    [372](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:372)     batch_size=batch_size,
    [373](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:373)     max_displacement_um=max_displacement_um,
    [374](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:374)     normalized_xcorr=normalized_xcorr,
    [375](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:375)     centered_xcorr=centered_xcorr,
    [376](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:376)     corr_threshold=corr_threshold,
    [377](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:377)     time_horizon_s=time_horizon_s,
    [378](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:378)     bin_duration_s=bin_duration_s,
    [379](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:379)     progress_bar=False,
    [380](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:380) )
    [382](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:382) if spatial_prior:
    [383](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:383)     all_pairwise_displacements[i] = pairwise_displacement

File [m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\sortingcomponents\motion_estimation.py:865](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:865), in compute_pairwise_displacement(motion_hist, bin_um, method, weight_scale, error_sigma, conv_engine, torch_device, batch_size, max_displacement_um, corr_threshold, time_horizon_s, normalized_xcorr, centered_xcorr, bin_duration_s, progress_bar, window)
    [862](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:862) correlation = np.empty((size, size), dtype=motion_hist.dtype)
    [864](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:864) for i in xrange(0, size, batch_size):
--> [865](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:865)     corr = normxcorr1d(
    [866](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:866)         motion_hist_engine,
    [867](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:867)         motion_hist_engine[i : i + batch_size],
    [868](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:868)         weights=window_engine,
    [869](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:869)         padding=possible_displacement.size // 2,
    [870](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:870)         conv_engine=conv_engine,
    [871](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:871)         normalized=normalized_xcorr,
    [872](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:872)         centered=centered_xcorr,
    [873](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:873)     )
    [874](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:874)     if conv_engine == "torch":
    [875](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:875)         max_corr, best_disp_inds = torch.max(corr, dim=2)

File [m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\sortingcomponents\motion_estimation.py:1415](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:1415), in normxcorr1d(template, x, weights, centered, normalized, padding, conv_engine)
   [1407](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:1407)     wt = template[:, None, :] * weights
   [1409](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:1409) # conv1d valid rule:
   [1410](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:1410) # (B,1,L),(O,1,L)->(B,O,L)
   [1411](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:1411) 
   [1412](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:1412) # compute expectations
   [1413](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:1413) # how many points in each window? seems necessary to normalize
   [1414](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:1414) # for numerical stability.
-> [1415](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:1415) N = conv1d(ones, weights, padding=padding)
   [1416](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:1416) if centered:
   [1417](file:///M:/Monkey_Python/SIenv/lib/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:1417)     Et = conv1d(ones, wt, padding=padding)

RuntimeError: Only zero batch or zero channel inputs are supported, but got input shape: [1, 1, 1, 0]
borrepp commented 6 months ago

If I changed the conv_engine to "numpy", I go the following error:


ValueError Traceback (most recent call last) Cell In[47], line 18 13 peaks = np.load(os.path.join(motion_folder, 'peaks.npy')) 15 peak_localization = np.load(os.path.join(motion_folder, 'peak_locations.npy')) ---> 18 motion, temporal_bins, spatial_bins = estimate_motion(recording=rec_filter_file, 19 peaks=peaks, 20 peak_locations=peak_localization, 21 method="decentralized", 22 pairwise_displacement_method = 'conv', 23 conv_engine="numpy", 24 batch_size = 1, 25 progress_bar = True, 26 verbose = True) 28 print(motion, temporal_bins, spatial_bins)

File m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\sortingcomponents\motion_estimation.py:152, in estimate_motion(recording, peaks, peak_locations, direction, bin_duration_s, bin_um, margin_um, rigid, win_shape, win_step_um, win_sigma_um, post_clean, speed_threshold, sigma_smooth_s, method, output_extra_check, progress_bar, upsample_to_histogram_bin, verbose, method_kwargs) 150 # run method 151 method_class = estimate_motion_methods[method] --> 152 motion, temporal_bins = method_class.run( 153 recording, 154 peaks, 155 peak_locations, 156 direction, 157 bin_duration_s, 158 bin_um, 159 spatial_bin_edges, 160 non_rigid_windows, 161 verbose, 162 progress_bar, 163 extra_check, 164 method_kwargs, 165 ) 167 # replace nan by zeros 168 motion[np.isnan(motion)] = 0

File m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\sortingcomponents\motion_estimation.py:363, in DecentralizedRegistration.run(cls, recording, peaks, peak_locations, direction, bin_duration_s, bin_um, spatial_bin_edges, non_rigid_windows, verbose, progress_bar, extra_check, histogram_depth_smooth_um, histogram_time_smooth_s, pairwise_displacement_method, max_displacement_um, weight_scale, error_sigma, conv_engine, torch_device, batch_size, corr_threshold, time_horizon_s, convergence_method, soft_weights, normalized_xcorr, centered_xcorr, temporal_prior, spatial_prior, force_spatial_median_continuity, reference_displacement, reference_displacement_time_s, robust_regression_sigma, lsqr_robust_n_iter, weight_with_amplitude) 360 if verbose: 361 print(f"Computing pairwise displacement: {i + 1} / {len(non_rigid_windows)}") --> 363 pairwise_displacement, pairwise_displacement_weight = compute_pairwise_displacement( 364 motion_histogram[:, window_slice], 365 bin_um, 366 window=win[window_slice], 367 method=pairwise_displacement_method, 368 weight_scale=weight_scale, 369 error_sigma=error_sigma, 370 conv_engine=conv_engine, 371 torch_device=torch_device, 372 batch_size=batch_size, 373 max_displacement_um=max_displacement_um, 374 normalized_xcorr=normalized_xcorr, 375 centered_xcorr=centered_xcorr, 376 corr_threshold=corr_threshold, 377 time_horizon_s=time_horizon_s, 378 bin_duration_s=bin_duration_s, 379 progress_bar=False, 380 ) 382 if spatial_prior: 383 all_pairwise_displacements[i] = pairwise_displacement

File m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\sortingcomponents\motion_estimation.py:865, in compute_pairwise_displacement(motion_hist, bin_um, method, weight_scale, error_sigma, conv_engine, torch_device, batch_size, max_displacement_um, corr_threshold, time_horizon_s, normalized_xcorr, centered_xcorr, bin_duration_s, progress_bar, window) 862 correlation = np.empty((size, size), dtype=motion_hist.dtype) 864 for i in xrange(0, size, batch_size): --> 865 corr = normxcorr1d( 866 motion_hist_engine, 867 motion_hist_engine[i : i + batch_size], 868 weights=window_engine, 869 padding=possible_displacement.size // 2, 870 conv_engine=conv_engine, 871 normalized=normalized_xcorr, 872 centered=centered_xcorr, 873 ) 874 if conv_engine == "torch": 875 max_corr, best_disp_inds = torch.max(corr, dim=2)

File m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\sortingcomponents\motion_estimation.py:1415, in normxcorr1d(template, x, weights, centered, normalized, padding, conv_engine) 1407 wt = template[:, None, :] * weights 1409 # conv1d valid rule: 1410 # (B,1,L),(O,1,L)->(B,O,L) 1411 1412 # compute expectations 1413 # how many points in each window? seems necessary to normalize 1414 # for numerical stability. -> 1415 N = conv1d(ones, weights, padding=padding) 1416 if centered: 1417 Et = conv1d(ones, wt, padding=padding)

File m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\sortingcomponents\motion_estimation.py:1476, in scipy_conv1d(input, weights, padding) 1474 for m in range(n): 1475 for c in range(c_out): -> 1476 output[m, c] = correlate(input[m, 0], weights[c, 0], mode=mode) 1478 return output

File m:\Monkey_Python\SIenv\lib\site-packages\scipy\signal_signaltools.py:241, in correlate(in1, in2, mode, method) 239 # this either calls fftconvolve or this function with method=='direct' 240 if method in ('fft', 'auto'): --> 241 return convolve(in1, _reverse_and_conj(in2), mode, method) 243 elif method == 'direct': 244 # fastpath to faster numpy.correlate for 1d inputs when possible 245 if _np_conv_ok(in1, in2, mode):

File m:\Monkey_Python\SIenv\lib\site-packages\scipy\signal_signaltools.py:1428, in convolve(in1, in2, mode, method) 1425 elif method == 'direct': 1426 # fastpath to faster numpy.convolve for 1d inputs when possible 1427 if _np_conv_ok(volume, kernel, mode): -> 1428 return np.convolve(volume, kernel, mode) 1430 return correlate(volume, _reverse_and_conj(kernel), mode, 'direct') 1431 else:

File <__array_function__ internals>:200, in convolve(*args, **kwargs)

File m:\Monkey_Python\SIenv\lib\site-packages\numpy\core\numeric.py:848, in convolve(a, v, mode) 846 a, v = v, a 847 if len(a) == 0: --> 848 raise ValueError('a cannot be empty') 849 if len(v) == 0: 850 raise ValueError('v cannot be empty')

ValueError: a cannot be empty

borrepp commented 2 months ago

Hi all,

I think the error I had came from the wrong units of the channel location property (millimeters instead of micrometers). I'm using an NWB file with location units in mm. After extracting the NWB recording object, I updated the location property to microns and I was able to run motion correction (step by step).

Thanks¡¡