axon delay fit + presentation

sam_experiments.py

@@ -173,6 +173,11 @@ def sam_analysis(fit_path):
 
     cell_stds = []
     model_stds = []
+
+    approx_offset = approximate_axon_delay_in_idx(cell_data, model)
+    print("Approx offset idx:", approx_offset)
+    print("Approx offset ms:", (approx_offset * step_size) * 1000)
+
     for mod_freq in sorted(u_delta_freqs):
         # TODO problem of cutting the pdf as in some cases the pdf is shorter than 1 modulation frequency period!
         #  length info wrong ? always at least one period?
@@ -213,7 +218,11 @@ def sam_analysis(fit_path):
         final_model_mean = np.mean(model_means, axis=0)
         cell_stds.append(np.std(final_cell_mean))
         model_stds.append(np.std(final_model_mean))
-        final_model_mean_phase_corrected = correct_phase(final_cell_mean, final_model_mean, step_size)
+        # offset, final_model_mean_phase_corrected = correct_phase(final_cell_mean, final_model_mean, step_size)
+        # print("Offset:", offset)
+        # print("modfreq:", mod_freq)
+
+        final_model_mean_phase_corrected = np.roll(final_model_mean, approx_offset)
 
         # PLOT EVERY MOD FREQ
         fig, axes = plt.subplots(1, 5, figsize=(15, 5), sharex=True)
@@ -280,7 +289,73 @@ def correct_phase(cell_mean, model_mean, step_size):
             lowest_err = abs
             roll_idx = roll_by
 
-    return np.roll(model_mean, roll_idx)
+    return roll_idx, np.roll(model_mean, roll_idx)
+
+
+def approximate_axon_delay_in_idx(cell_data, model):
+
+    lowest_mod_freq = 80
+    highest_mod_freq = 150
+
+    eod_freq = cell_data.get_eod_frequency()
+    step_size = cell_data.get_sampling_interval()
+
+    durations = cell_data.get_sam_durations()
+    contrasts = cell_data.get_sam_contrasts()
+    u_contrasts = np.unique(contrasts)
+    spiketimes = cell_data.get_sam_spiketimes()
+    delta_freqs = cell_data.get_sam_delta_frequencies()
+    u_delta_freqs = np.unique(delta_freqs)
+
+    used_mod_freqs = []
+    axon_delays = []
+
+    for mod_freq in sorted(u_delta_freqs):
+        # Only use "stable" mod_freqs to approximate the axon delay
+        if not lowest_mod_freq <= mod_freq <= highest_mod_freq:
+            continue
+
+        if 1/mod_freq > durations[0] / 4:
+            print("skipped mod_freq: {}".format(mod_freq))
+            print("Duration: {} while mod_freq period: {:.2f}".format(durations[0], 1/mod_freq))
+            continue
+        mfreq_data = {}
+        cell_means = []
+        model_means = []
+        for c in u_contrasts:
+            mfreq_data[c] = []
+
+        for i in range(len(delta_freqs)):
+            if delta_freqs[i] != mod_freq:
+                continue
+
+            if len(spiketimes[i]) > 1:
+                print("There are more spiketimes in one 'point'! Only the first was used! ")
+            spikes = spiketimes[i][0]
+
+            cell_pdf = spiketimes_calculate_pdf(spikes, step_size)
+
+            cell_cuts = cut_pdf_into_periods(cell_pdf, 1/mod_freq, step_size, factor=1.0)
+            cell_mean = np.mean(cell_cuts, axis=0)
+            cell_means.append(cell_mean)
+
+            stimulus = SAM(eod_freq, contrasts[i] / 100, mod_freq)
+            v1, spikes_model = model.simulate(stimulus, durations[i] * 4)
+            model_pdf = spiketimes_calculate_pdf(spikes_model, step_size)
+            model_cuts = cut_pdf_into_periods(model_pdf, 1/mod_freq, step_size, factor=1.0)
+            model_mean = np.mean(model_cuts, axis=0)
+            model_means.append(model_mean)
+
+        final_cell_mean = np.mean(cell_means, axis=0)
+        final_model_mean = np.mean(model_means, axis=0)
+
+        offset, final_model_mean_phase_corrected = correct_phase(final_cell_mean, final_model_mean, step_size)
+
+        used_mod_freqs.append(mod_freq)
+        axon_delays.append(offset)
+
+    mean_delay = np.mean(axon_delays)
+    return int(round(mean_delay))
 
 
 def generate_pdf(model, stimulus, trials=4, sim_length=3, kernel_width=0.005):