add noise to the dendrite not the outputcd PycharmProjects/neuronModel/

models/LIFACnoise.py

@@ -300,10 +300,10 @@ def simulate_fast(rectified_stimulus_array, total_time_s, parameters: np.ndarray
         noise = noise_strength * noise_value / np.sqrt(step_size)
 
         input_voltage[i] = input_voltage[i - 1] + (
-                    (-input_voltage[i - 1] + rectified_stimulus_array[i]) / dend_tau) * step_size
+                    (-input_voltage[i - 1] + rectified_stimulus_array[i] + noise) / dend_tau) * step_size
 
         output_voltage[i] = output_voltage[i - 1] + ((v_base - output_voltage[i - 1] + v_offset + (
-                    input_voltage[i] * input_scaling) - adaption[i - 1] + noise) / mem_tau) * step_size
+                    input_voltage[i] * input_scaling) - adaption[i - 1]) / mem_tau) * step_size
 
         adaption[i] = adaption[i - 1] + ((-adaption[i - 1]) / tau_a) * step_size
 

sam_experiments.py

@@ -31,6 +31,8 @@ def run_sam_analysis_for_all_cells(folder):
     count = 0
     for item in sorted(os.listdir(folder)):
         cell_folder = os.path.join(folder, item)
+        if not os.path.isdir(cell_folder):
+            continue
         # fit = get_best_fit(cell_folder, use_comparable_error=False)
         # cell_data = fit.get_cell_data()
         #
@@ -42,7 +44,6 @@ def run_sam_analysis_for_all_cells(folder):
     print(count)
 
 
-
 def test_model_response(model: LifacNoiseModel, eod_freq, contrast, modulation_frequencies):
 
     stds = []

test.py

@@ -1,37 +1,81 @@
 import os
 from parser.CellData import CellData
+import numpy as np
+from fitting.ModelFit import ModelFit, get_best_fit
 # from plottools.axes import labelaxes_params
 import matplotlib.pyplot as plt
 
+colors = ["black", "red", "blue", "orange", "green"]
 
-data_folder = "./data/final/"
-for cell in sorted(os.listdir(data_folder)):
-    print(cell)
-    cell_folder = os.path.join(data_folder, cell)
-    if not os.path.exists(os.path.join(cell_folder, "samspikes1.dat")):
-        continue
 
-    cell_data = CellData(cell_folder)
-    sam_spikes = cell_data.get_sam_spiketimes()
-    delta_freqs = cell_data.get_sam_delta_frequencies()
-    # [time_traces, v1_traces, eod_traces, local_eod_traces, stimulus_traces]
-    [time_traces, v1_traces, eod_traces, local_eod_traces, stimulus_traces] = cell_data.get_sam_traces()
-    print(len(time_traces))
-    # for i in range(len(delta_freqs)):
-    #
-    #     fig, axes = plt.subplots(2, 1, sharex="all")
-    #
-    #     axes[0].plot(time_traces[i], local_eod_traces[i])
-    #     axes[0].set_title("Local EOD - dF {}".format(delta_freqs[i]))
-    #     axes[1].plot(time_traces[i], v1_traces[i])
-    #     axes[1].set_title("v1 trace")
-    #     axes[1].eventplot(sam_spikes[i], lineoffsets=max(v1_traces[i]))
-    #     plt.show()
-    #     plt.close()
-    #     break
+def main():
+    # sam_tests()
 
+    fit = get_best_fit("results/final_sam2/2012-12-20-ae-invivo-1/")
+    fit.generate_master_plot()
 
 
+def sam_tests():
+    data_folder = "./data/final/"
+    for cell in sorted(os.listdir(data_folder)):
+        print(cell)
+        cell_folder = os.path.join(data_folder, cell)
+        if not os.path.exists(os.path.join(cell_folder, "samspikes1.dat")):
+            continue
 
+        cell_data = CellData(cell_folder)
+        sampling_rate = int(round(1 / cell_data.get_sampling_interval()))
+        sam_spikes = cell_data.get_sam_spiketimes()
+        delta_freqs = cell_data.get_sam_delta_frequencies()
 
+        [time_traces, v1_traces, eod_traces, local_eod_traces, stimulus_traces] = cell_data.get_sam_traces()
+        print(len(time_traces))
+        for i in range(len(delta_freqs)):
 
+            fig, axes = plt.subplots(2, 1, sharex="all")
+
+            axes[0].plot(time_traces[i], local_eod_traces[i])
+            axes[0].set_title("Local EOD - dF {}".format(delta_freqs[i]))
+            axes[1].plot(time_traces[i], v1_traces[i])
+            axes[1].set_title("v1 trace")
+            ah_spike = average_spike_height(sam_spikes, v1_traces[i], sampling_rate)
+            for j, idx in enumerate(get_x_best(ah_spike)):
+                axes[1].eventplot(sam_spikes[idx], lineoffsets=max(v1_traces[i] + 1.5 * (j + 1)),
+                                  colors=colors[j % len(colors)])
+            plt.show()
+            plt.close()
+        break
+
+
+def average_spike_height(spike_trains, local_eod, sampling_rate):
+    average_height = []
+    for spikes_train in spike_trains:
+        indices = np.array([s * sampling_rate for s in spikes_train[0]], dtype=np.int)
+        local_eod = np.array(local_eod)
+        spike_values = [local_eod[i] for i in indices if i < len(local_eod)]
+        average_height.append(np.mean(spike_values))
+
+    return average_height
+
+
+def get_x_best(average_heights, x=5):
+    biggest_idx = []
+    biggest_heights = []
+
+    for i, height in enumerate(average_heights):
+
+        if len(biggest_idx) < x:
+            biggest_idx.append(i)
+            biggest_heights.append(height)
+        elif height > min(biggest_heights):
+            mini = np.argmin(biggest_heights)
+            biggest_heights[mini] = height
+            biggest_idx[mini] = i
+
+    biggest_heights, biggest_idx = (list(t) for t in zip(*sorted(zip(biggest_heights, biggest_idx), reverse=True)))
+    print(biggest_heights)
+    return biggest_idx
+
+
+if __name__ == '__main__':
+    main()