continue playing around with inter-session displacement generator.

Author: JoeZiminski

src/spikeinterface/core/generate.py

@@ -1633,6 +1633,10 @@ def generate_templates(
         channel_factors = alpha * np.exp(-distances / spatial_decay)
         wfs = wf[:, np.newaxis] * channel_factors[np.newaxis, :]
 
+        #   if u in [3, 4]:
+        #      import matplotlib.pyplot as plt
+        #     breakpoint()
+
         # This mimic a propagation delay for distant channel
         propagation_speed = params["propagation_speed"][u]
         if propagation_speed is not None:
@@ -1901,6 +1905,13 @@ def get_traces(
                 upsample_ind = self.upsample_vector[i]
                 template = self.templates[unit_ind, :, :, upsample_ind]
 
+            #       if unit_ind == 4:
+            #          from spikeinterface.core import order_channels_by_depth
+            #         import matplotlib.pyplot as plt
+            #        off = np.atleast_2d(np.arange(template.shape[1])* 10)
+            #       plt.plot(template); plt.show()
+            #      breakpoint()
+
             if channel_indices is not None:
                 template = template[:, channel_indices]
 
@@ -1926,6 +1937,10 @@ def get_traces(
                 wf = wf * self.amplitude_vector[i]
             traces[start_traces:end_traces] += wf.astype(traces.dtype, copy=False)
 
+        # TODO: does templates on one channel not overwrite template on another channel?
+        # no, because the waveform is ADDED, but I guess the noise is added multiple times?
+        # but it is gaussian, so shouldn't matter.
+
         return traces.astype(self.dtype, copy=False)
 
     def get_num_samples(self) -> int:

src/spikeinterface/generation/debugging/debugging_session_displacement.py

@@ -0,0 +1,104 @@
+import spikeinterface.full as si
+from spikeinterface.generation.drifting_generator import generate_drifting_recording
+from spikeinterface.generation.session_displacement_generator import generate_inter_session_displacement_recordings
+import matplotlib.pyplot as plt
+
+if False:
+
+    _, raw_recording, _ = generate_drifting_recording(
+        num_units=10,
+        duration=50,
+        generate_sorting_kwargs=dict(firing_rates=(15, 25), refractory_period_ms=4.0),
+        seed=42,
+        generate_displacement_vector_kwargs=dict(
+            motion_list=[
+                dict(drift_mode="zigzag", non_rigid_gradient=0.01, t_start_drift=5, t_end_drift=45, period_s=10),
+            ],
+        ),
+    )
+
+# TODO: understand why sometimes the units are spraed all over the probe.
+# In my example, for example, they are clearly ok on the raw image.
+# but detect peaks is terrible.
+import numpy as np
+
+default_unit_params_range = dict(
+    alpha=(100.0, 500.0),
+    depolarization_ms=(0.09, 0.14),
+    repolarization_ms=(0.5, 0.8),
+    recovery_ms=(1.0, 1.5),
+    positive_amplitude=(0.1, 0.25),
+    smooth_ms=(0.03, 0.07),
+    spatial_decay=(20, 40),
+    propagation_speed=(250.0, 350.0),  # um  / ms
+    b=(0.5, 1),  # (0.5, 1)
+    c=(0.5, 1),  # (0.5, 1)
+    x_angle=(0, np.pi),
+    y_angle=(0, np.pi),
+    z_angle=(0, np.pi),  # (0, 2)
+)
+
+if False:
+    rec, _ = si.generate_ground_truth_recording(
+        num_units=5,
+        durations=[25],
+        num_channels=128,
+        #  generate_sorting_kwargs=dict(firing_rates=(200, 300), refractory_period_ms=2.0),
+        #  noise_kwargs=dict(noise_levels=5),
+        #  generate_templates_kwargs=dict(unit_params=dict(alpha=300)),  # TODO: increase alpha default?
+        seed=44,
+    )
+    rec_list = [rec]
+
+default_unit_params_range["alpha"] = (400, 500)  # do this or change the margin on generate_unit_locations_kwargs
+
+rec_list, _ = generate_inter_session_displacement_recordings(
+    non_rigid_gradient=None,  # 0.05,
+    num_units=5,
+    rec_durations=(25, 25, 25),  # TODO: checks on inputs
+    rec_shifts=(
+        0,
+        200,
+        400,
+    ),  # WTF happening at +100 um?? maybe needs to be discrete in terms of channels or it is ignored? shouldnt matter...
+    generate_sorting_kwargs=dict(firing_rates=(149, 150), refractory_period_ms=4.0),
+    generate_templates_kwargs=dict(unit_params=default_unit_params_range, ms_before=1.5, ms_after=3),  # mode=sphere
+    seed=44,
+    generate_unit_locations_kwargs=dict(
+        margin_um=0.0,  # if this is say 20, then units go off the edge of the probe and are such low amplitude they are not picked up.
+        minimum_z=5.0,
+        maximum_z=45.0,
+        minimum_distance=18.0,
+        max_iteration=100,
+        distance_strict=False,
+    ),
+)
+
+from spikeinterface.sortingcomponents.peak_detection import detect_peaks
+from spikeinterface.sortingcomponents.peak_localization import localize_peaks
+
+# < 150 might be a problem for detecting alpha
+for rec in rec_list:
+
+    si.plot_traces(rec, time_range=(0, 1))
+    # plt.show()
+
+    peaks = detect_peaks(rec, method="locally_exclusive")
+    peak_locs = localize_peaks(rec, peaks, method="grid_convolution")
+
+    print(peaks.shape)
+
+    si.plot_drift_raster_map(peaks=peaks, peak_locations=peak_locs, recording=rec, clim=(-300, 0))
+    plt.show()
+
+    # TODO: could try and create peaks, peak locs from
+    # unit locations and rec.spike_vector...
+    rec.spike_vector.size  # (sample, unit, segment)
+    # peaks (sample, channel, amplitude, segment)
+
+# Why are units dropping out (distance too far from probe)
+# Why are the amplitudes changing (clim)
+# why are the traces look like they have too many units (current settings)
+# TODO: I think it is OK when it goes over the edge, just need to check explicitly what happens.
+# It will be greater than margin, etc. Is it explicitly clipped or just allowed to go outside of range?
+# I guess this will ignore the margin. That's fine. I think here is it just going outside of range.