SpikeInterface / spikeinterface

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

ValueError: could not broadcast input array - motion estimation #1417

Open bramn22 opened 1 year ago

bramn22 commented 1 year ago

Hi ,

When running the motion_estimation function, I get the following error:

"/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/spikeinterface/sortingcomponents/motion_estimation.py:1011: RuntimeWarning: invalid value encountered in true_divide corr /= npx.sqrt(var_x * var_template) Traceback (most recent call last): File "/home/farrowlab/conda-farrowlab/preprocess_npx.py", line 216, in main(cfg) File "/home/farrowlab/conda-farrowlab/preprocess_npx.py", line 208, in main rec = motion_correction_full(rec, out_path=out_path, shift_idx=shift_idx, visualize=True) File "/home/farrowlab/conda-farrowlab/preprocess_npx.py", line 118, in motion_correction_full motion, temporal_bins, spatial_bins, extra_check = spikeinterface.sortingcomponents.motion_estimation.estimate_motion( File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/spikeinterface/sortingcomponents/motion_estimation.py", line 127, in estimate_motion motion, temporal_bins = method_class.run(recording, peaks, peak_locations, direction, bin_duration_s, bin_um, File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/spikeinterface/sortingcomponents/motion_estimation.py", line 263, in run motion[:, i] = compute_global_displacement(pairwise_displacement, ValueError: could not broadcast input array from shape (1124,) into shape (1125,)"

function used: motion, temporal_bins, spatial_bins, extra_check = spikeinterface.sortingcomponents.motion_estimation.estimate_motion( rec, peaks, peak_locations=peak_locations, direction='y', bin_duration_s=5, bin_um=10., rigid=False, win_shape='gaussian', win_step_um=50., win_sigma_um=150, progress_bar=True, output_extra_check=True, verbose=True, )

Thank you, Bram

samuelgarcia commented 1 year ago

Hi bram, Thanks for the report. Which version do you have ? Is it installed from the github sources ?

bramn22 commented 1 year ago

Hi Samuel,

I am using 0.97.0. I just noticed that there is 0.97.1 so I am trying that now.

bramn22 commented 1 year ago

With version 0.97.1, an error now occurs during the writing of the binary recording:

"write_binary_recording: 20%|█▉ | 1124/5624 [06:04<24:20, 3.08it/s] concurrent.futures.process._RemoteTraceback: """ Traceback (most recent call last): File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/concurrent/futures/process.py", line 243, in _process_worker r = call_item.fn(*call_item.args, *call_item.kwargs) File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/concurrent/futures/process.py", line 202, in _process_chunk return [fn(args) for args in chunk] File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/concurrent/futures/process.py", line 202, in return [fn(*args) for args in chunk] File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/spikeinterface/core/job_tools.py", line 399, in function_wrapper return _func(segment_index, start_frame, end_frame, _worker_ctx) File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/spikeinterface/core/core_tools.py", line 209, in _write_binary_chunk traces = recording.get_traces(start_frame=start_frame, end_frame=end_frame, segment_index=segment_index, File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/spikeinterface/core/baserecording.py", line 160, in get_traces traces = rs.get_traces(start_frame=start_frame, end_frame=end_frame, channel_indices=channel_indices) File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/spikeinterface/sortingcomponents/motion_correction.py", line 335, in get_traces trace2 = correct_motion_on_traces(traces, times, self.channel_locations, self.motion, File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/spikeinterface/sortingcomponents/motion_correction.py", line 124, in correct_motion_on_traces f = scipy.interpolate.interp1d(spatial_bins, motion[bin_ind, :], kind='linear', IndexError: index 1125 is out of bounds for axis 0 with size 1125 """

The above exception was the direct cause of the following exception:

Traceback (most recent call last): File "/home/farrowlab/conda-farrowlab/preprocess_npx.py", line 216, in main(cfg) File "/home/farrowlab/conda-farrowlab/preprocess_npx.py", line 209, in main save_file(rec, os.path.join(out_path, 'data'), binary=False) File "/home/farrowlab/conda-farrowlab/preprocess_npx.py", line 44, in save_file rec.save(folder=out_path, n_jobs=10, chunk_duration='1s', progres_bar=True) File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/spikeinterface/core/base.py", line 621, in save loaded_extractor = self.save_to_folder(kwargs) File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/spikeinterface/core/base.py", line 700, in save_to_folder cached = self._save(folder=folder, verbose=verbose, save_kwargs) File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/spikeinterface/core/baserecording.py", line 297, in _save write_binary_recording(self, file_paths=file_paths, dtype=dtype, **job_kwargs) File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/spikeinterface/core/core_tools.py", line 280, in write_binary_recording executor.run() File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/spikeinterface/core/job_tools.py", line 364, in run for res in results: File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/site-packages/tqdm/std.py", line 1195, in iter for obj in iterable: File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/concurrent/futures/process.py", line 559, in _chain_from_iterable_of_lists for element in iterable: File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/concurrent/futures/_base.py", line 609, in result_iterator yield fs.pop().result() File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/concurrent/futures/_base.py", line 439, in result return self.__get_result() File "/home/farrowlab/conda-farrowlab/anaconda3/envs/si_bram/lib/python3.9/concurrent/futures/_base.py", line 391, in __get_result raise self._exception IndexError: index 1125 is out of bounds for axis 0 with size 1125"

alejoe91 commented 10 months ago

@bramn22 sorry for the super delay on this. We actually made quite some improvements to the motion/drift in more recent releases (>=0.98). Are you still experiencing issues?

borrepp commented 5 months ago

Hello @alejoe91 ,

I'm getting the same error, basically the peaks.py and peak_locations.py ended with different shapes. In my case is:

ValueError: could not broadcast input array from shape (588052,) into shape (588055,)

I'm using windows 11, spikeInterface version = 0.100.0

alejoe91 commented 5 months ago

ok, what is the output of print(len(sorting.to_spike_vector()))?

borrepp commented 5 months ago

I don't have a sorting object. The peaks.py and peak_locations.py were created while following the pipeline for motion correction.

Sorry in advance for the following comment, I'm not sure if this is a proper way to update this issue or if should I open a new one. I have encountered another error while trying to get the motion correction to work. Thanks in advance for your help.

I tried to run it step by step as suggested here: https://spikeinterface.readthedocs.io/en/0.100.1/modules/motion_correction.html#low-level-api

But, sometimes I get an error when running "localize_peaks".

I have the following filtered recording (originally extracted from NWB): #################################################################################################

print(rec_filter_file) BinaryFolderRecording: 32 channels - 30.0kHz - 1 segments - 44,771,040 samples 1,492.37s (24.87 minutes) - float32 dtype - 5.34 GiB #################################################################################################

One iteration ran without a problem with the following settings: ################################################################################################# peaks = detect_peaks(recording = rec_filter_file, method="by_channel_torch", peak_sign="both", detect_threshold=5, gather_mode="memory", **job_kwargs)

print(peaks.shape) (588055,) print(peaks) [( 115, 26, -74.86369324, 0) ( 116, 24, -44.0439415 , 0) ( 140, 17, 37.53662491, 0) ... (44770724, 26, -77.1676178 , 0) (44770724, 28, -40.46783829, 0) (44770859, 3, -64.48331451, 0)] #################################################################################################

After that, I ran localize_peaks with the following settings: ################################################################################################# peak_locations = localize_peaks(recording=rec_filter_file, peaks = peaks, method="monopolar_triangulation", radius_um= 75.0, max_distance_um= 150.0, optimizer="minimize_with_log_penality", enforce_decrease=True, feature="ptp", **job_kwargs)

print(peak_localization.shape) (588052,) print(peak_localization) [(-3.23900826e-08, -0.447295 , 1. , 43.91960107) ( 0.00000000e+00, 125.20039629, 11.46852403, 3703.97416739) (-2.10297796e-09, -1.14823871, 1.28532754, 35.90787696) ... ( 6.99432775e-09, -0.19747129, 1.20437011, 51.43084121) ( 0.00000000e+00, 1.30129915, 2.00380814, 95.6332904 ) ( 0.00000000e+00, -6.14027055, 2.34544175, 131.50203496)] #################################################################################################

As I mentioned, the "peaks" and "peak_locations" ended up with different shapes, which caused the original error while trying to run "estimate_motion"

However, when I ran the same steps but changed the "peak_sign" to "neg", I got fewer peaks (as expected?) But, I got an error when I tried to get peak_locations.

################################################################################################# peaks = detect_peaks(recording = rec_filter_file, method="by_channel_torch", peak_sign="neg", detect_threshold=5, gather_mode="memory", **job_kwargs)

print(peaks.shape) (434996,) print(peaks) [( 115, 26, -74.86369324, 0) ( 116, 24, -44.0439415 , 0) ( 272, 25, -42.79653168, 0) ... (44770724, 26, -77.1676178 , 0) (44770724, 28, -40.46783829, 0) (44770859, 3, -64.48331451, 0)] #################################################################################################

################################################################################################# peak_locations = localize_peaks(recording=rec_filter_file, peaks = peaks, method="monopolar_triangulation", radius_um= 75.0, max_distance_um= 150.0, optimizer="minimize_with_log_penality", enforce_decrease=True, feature="ptp", **job_kwargs)


_RemoteTraceback Traceback (most recent call last) _RemoteTraceback: """ Traceback (most recent call last): File "C:\Users\josev\AppData\Local\Programs\Python\Python38\lib\concurrent\futures\process.py", line 239, in _process_worker r = call_item.fn(*call_item.args, *call_item.kwargs) File "C:\Users\josev\AppData\Local\Programs\Python\Python38\lib\concurrent\futures\process.py", line 198, in _process_chunk return [fn(args) for args in chunk] File "C:\Users\josev\AppData\Local\Programs\Python\Python38\lib\concurrent\futures\process.py", line 198, in return [fn(args) for args in chunk] File "m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\core\job_tools.py", line 423, in function_wrapper return _func(segment_index, start_frame, end_frame, _worker_ctx) File "m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\core\node_pipeline.py", line 522, in _compute_peak_pipeline_chunk node_output = node.compute(traces_chunk, node_input_args) File "m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\core\node_pipeline.py", line 322, in compute waveforms = traces[peaks["sample_index"][:, None] + np.arange(-self.nbefore, self.nafter)] IndexError: index 30030 is out of bounds for axis 0 with size 30030 """

The above exception was the direct cause of the following exception:

IndexError Traceback (most recent call last) Cell In[12], line 2 1 # LOCALIZE Peaks ----> 2 peak_locations = localize_peaks(recording=rec_filter_file, peaks = peaks, method="monopolar_triangulation", 3 radius_um= 75.0, max_distance_um= 150.0, optimizer="minimize_with_log_penality", 4 enforce_decrease=True, feature="ptp", **job_kwargs) 5 """ 6 method="monopolar_triangulation", 7 radius_um=75.0, (...) 20 percentile=5.0, 21 """ 22 print('saving peak_locations')

File m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\sortingcomponents\peak_localization.py:113, in localize_peaks(recording, peaks, method, ms_before, ms_after, kwargs) 88 """Localize peak (spike) in 2D or 3D depending the method. 89 90 When a probe is 2D then: (...) 110 The dtype depends on the method. ("x", "y") or ("x", "y", "z", "alpha"). 111 """ 112 peak_retriever = PeakRetriever(recording, peaks) --> 113 peak_locations = _run_localization_from_peak_source( 114 recording, peak_retriever, method=method, ms_before=ms_before, ms_after=ms_after, kwargs 115 ) 116 return peak_locations

File m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\sortingcomponents\peak_localization.py:82, in _run_localization_from_peak_source(recording, peak_source, method, ms_before, ms_after, kwargs) 75 pipeline_nodes = [ 76 peak_source, 77 extract_dense_waveforms, 78 LocalizeGridConvolution(recording, parents=[peak_source, extract_dense_waveforms], method_kwargs), 79 ] 81 job_name = f"localize peaks using {method}" ---> 82 peak_locations = run_node_pipeline(recording, pipeline_nodes, job_kwargs, job_name=job_name, squeeze_output=True) 84 return peak_locations

File m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\core\node_pipeline.py:471, in run_node_pipeline(recording, nodes, job_kwargs, job_name, mp_context, gather_mode, gather_kwargs, squeeze_output, folder, names) 459 init_args = (recording, nodes) 461 processor = ChunkRecordingExecutor( 462 recording, 463 _compute_peak_pipeline_chunk, (...) 468 **job_kwargs, 469 ) --> 471 processor.run() 473 outs = gather_func.finalize_buffers(squeeze_output=squeeze_output) 474 return outs

File m:\Monkey_Python\SIenv\lib\site-packages\spikeinterface\core\job_tools.py:385, in ChunkRecordingExecutor.run(self) 382 if self.progress_bar: 383 results = tqdm(results, desc=self.job_name, total=len(all_chunks)) --> 385 for res in results: 386 if self.handle_returns: 387 returns.append(res)

File m:\Monkey_Python\SIenv\lib\site-packages\tqdm\notebook.py:250, in tqdm_notebook.iter(self) 248 try: 249 it = super(tqdm_notebook, self).iter() --> 250 for obj in it: 251 # return super(tqdm...) will not catch exception 252 yield obj 253 # NB: except ... [ as ...] breaks IPython async KeyboardInterrupt

File m:\Monkey_Python\SIenv\lib\site-packages\tqdm\std.py:1181, in tqdm.iter(self) 1178 time = self._time 1180 try: -> 1181 for obj in iterable: 1182 yield obj 1183 # Update and possibly print the progressbar. 1184 # Note: does not call self.update(1) for speed optimisation.

File ~\AppData\Local\Programs\Python\Python38\lib\concurrent\futures\process.py:484, in _chain_from_iterable_of_lists(iterable) 478 def _chain_from_iterable_of_lists(iterable): 479 """ 480 Specialized implementation of itertools.chain.from_iterable. 481 Each item in iterable should be a list. This function is 482 careful not to keep references to yielded objects. 483 """ --> 484 for element in iterable: 485 element.reverse() 486 while element:

File ~\AppData\Local\Programs\Python\Python38\lib\concurrent\futures_base.py:611, in Executor.map..result_iterator() 608 while fs: 609 # Careful not to keep a reference to the popped future 610 if timeout is None: --> 611 yield fs.pop().result() 612 else: 613 yield fs.pop().result(end_time - time.monotonic())

File ~\AppData\Local\Programs\Python\Python38\lib\concurrent\futures_base.py:432, in Future.result(self, timeout) 430 raise CancelledError() 431 elif self._state == FINISHED: --> 432 return self.__get_result() 434 self._condition.wait(timeout) 436 if self._state in [CANCELLED, CANCELLED_AND_NOTIFIED]:

File ~\AppData\Local\Programs\Python\Python38\lib\concurrent\futures_base.py:388, in Future.get_result(self) 386 def get_result(self): 387 if self._exception: --> 388 raise self._exception 389 else: 390 return self._result

IndexError: index 30030 is out of bounds for axis 0 with size 30030

alejoe91 commented 5 months ago

Ah I think we have a bug in the torch detection! Can you try with the locally_exclusive?

borrepp commented 5 months ago

With "locally_exclusive" it runs fine¡¡¡¡ This time, "peaks" and "peak_locations" have the same shape. I think I will just set the "radius_um" to be smaller than the spacing of the probe that I have, which should behave like "by_channel"

However, when I try to run the next step: "estimate_motion" with "method_motion" = "decentralized". I have the error that I posted here:

https://github.com/SpikeInterface/spikeinterface/issues/2603

I will wait on that issue to not mix channels.

On the other hand, if I run "estimate_motion" with "method_motion" = "iterative_template". I have no errors. But when I run "interpolate_motion":

####################################################################################### rec_corrected = interpolate_motion(recording=rec_filter_file, motion=motion, temporal_bins=temporal_bins, spatial_bins=spatial_bins, border_mode="remove_channels", spatial_interpolation_method="kriging", sigma_um=30.) #######################################################################################

I got zero channels:

print(rec_corrected) InterpolateMotionRecording: 0 channels - 30.0kHz - 1 segments - 44,771,040 samples 1,492.37s (24.87 minutes) - float32 dtype - 0.00 B

#######################################################################################

Did I do something wrong?

Thanks again for your help¡¡¡

borrepp commented 5 months ago

By the way, with locally_exclusive_torch. I also got the same number of peaks and peak_locations.