add f zero freq curve error

Fitter.py

@@ -10,6 +10,8 @@ from warnings import warn
 from scipy.optimize import minimize
 import time
 
+import matplotlib.pyplot as plt
+
 
 class Fitter:
 
@@ -21,9 +23,13 @@ class Fitter:
             if "step_size" not in params:
                 self.base_model.set_variable("step_size", 0.00005)
 
+        self.best_parameters_found = []
+        self.smallest_error = np.inf
+
         #
         self.fi_contrasts = []
         self.eod_freq = 0
+        self.data_sampling_interval = -1
 
         self.sc_max_lag = 2
 
@@ -42,6 +48,11 @@ class Fitter:
         self.f_zero_slope_at_straight = 0
         self.f_zero_straight_contrast = 0
         self.f_zero_fit = []
+        self.f_zero_curve_contrast = 0
+        self.f_zero_curve_contrast_idx = -1
+        self.f_zero_curve_freq = np.array([])
+        self.f_zero_curve_time = np.array([])
+
 
         # self.tau_a = 0
 
@@ -50,6 +61,7 @@ class Fitter:
 
     def set_data_reference_values(self, cell_data: CellData):
         self.eod_freq = cell_data.get_eod_frequency()
+        self.data_sampling_interval = cell_data.get_sampling_interval()
 
         data_baseline = get_baseline_class(cell_data)
         data_baseline.load_values(cell_data.get_data_path())
@@ -61,8 +73,9 @@ class Fitter:
 
         fi_curve = get_fi_curve_class(cell_data, cell_data.get_fi_contrasts(), save_dir=cell_data.get_data_path())
         self.f_inf_slope = fi_curve.get_f_inf_slope()
+        contrasts = np.array(cell_data.get_fi_contrasts())
         if self.f_inf_slope < 0:
-            contrasts = np.array(cell_data.get_fi_contrasts()) * -1
+            contrasts = contrasts * -1
             # print("old contrasts:", cell_data.get_fi_contrasts())
             # print("new contrasts:", contrasts)
             contrasts = sorted(contrasts)
@@ -72,13 +85,22 @@ class Fitter:
         self.f_inf_values = fi_curve.f_inf_frequencies
         self.f_inf_slope = fi_curve.get_f_inf_slope()
 
-
         self.f_zero_values = fi_curve.f_zero_frequencies
         self.f_zero_fit = fi_curve.f_zero_fit
         # self.f_zero_slopes = [fi_curve.get_f_zero_fit_slope_at_stimulus_value(c) for c in self.fi_contrasts]
         self.f_zero_slope_at_straight = fi_curve.get_f_zero_fit_slope_at_straight()
         self.f_zero_straight_contrast = self.f_zero_fit[3]
 
+        max_contrast = max(contrasts)
+        test_contrast = 0.5 * max_contrast
+        diff_contrasts = np.abs(contrasts - test_contrast)
+
+        self.f_zero_curve_contrast_idx = np.argmin(diff_contrasts)
+        self.f_zero_curve_contrast = contrasts[self.f_zero_curve_contrast_idx]
+        times, freqs = fi_curve.get_mean_time_and_freq_traces()
+        self.f_zero_curve_freq = freqs[self.f_zero_curve_contrast_idx]
+        self.f_zero_curve_time = times[self.f_zero_curve_contrast_idx]
+
         # around 1/3 of the value at straight
         # self.f_zero_slope = fi_curve.get_fi_curve_slope_at(fi_curve.get_f_zero_and_f_inf_intersection())
 
@@ -100,7 +122,7 @@ class Fitter:
 
         # error_list = [error_bf, error_vs, error_sc, error_cv,
         #                       error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope]
-        error_weights = (0, 2, 2, 2, 1, 1, 1, 1, 0)
+        error_weights = (0, 2, 2, 2, 1, 1, 1, 1, 0, 1)
         fmin = minimize(fun=self.cost_function_all,
                         args=(error_weights,), x0=x0, method="Nelder-Mead",
                         options={"initial_simplex": initial_simplex, "xatol": 0.001, "maxfev": 600, "maxiter": 1200})
@@ -152,6 +174,11 @@ class Fitter:
         return fmin, self.base_model.get_parameters()
 
     def cost_function_all(self, X, error_weights=None):
+        for i in range(len(X)):
+            if X[i] < 0:
+                print("tried impossible value")
+                return 1000 + abs(X[i]) * 10000
+
         self.base_model.set_variable("mem_tau", X[0])
         self.base_model.set_variable("noise_strength", X[1])
         self.base_model.set_variable("input_scaling", X[2])
@@ -160,6 +187,11 @@ class Fitter:
         self.base_model.set_variable("dend_tau", X[5])
         self.base_model.set_variable("refractory_period", X[6])
 
+
+
+
+        # TODO add tests for parameters punish impossible values (immediate high error) but also add a slope towards valid points
+
         base_stimulus = SinusoidalStepStimulus(self.eod_freq, 0)
         # find right v-offset
         test_model = self.base_model.get_model_copy()
@@ -172,7 +204,11 @@ class Fitter:
 
         # [error_bf, error_vs, error_sc, error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope]
         error_list = self.calculate_errors(error_weights)
-        # print(sum(error_list))
+        print("sum: {:.2f}, ".format(sum(error_list)))
+        if sum(error_list) < self.smallest_error:
+            self.smallest_error = sum(error_list)
+            self.best_parameters_found = X
+
         return sum(error_list)
 
     def cost_function_without_ref_period(self, X, error_weights=None):
@@ -325,7 +361,7 @@ class Fitter:
         error_bf = abs((baseline_freq - self.baseline_freq) / self.baseline_freq)
         error_vs = abs((vector_strength - self.vector_strength) / 0.1)
         error_cv = abs((coefficient_of_variation - self.coefficient_of_variation) / 0.1)
-        error_bursty = (abs(burstiness - self.burstiness) / 0.02)
+        error_bursty = (abs(burstiness - self.burstiness) / 0.2)
 
         error_sc = 0
         for i in range(self.sc_max_lag):
@@ -333,6 +369,8 @@ class Fitter:
         # error_sc = error_sc / self.sc_max_lag
 
 
+
+
         error_f_inf_slope = abs((f_infinities_slope - self.f_inf_slope) / abs(self.f_inf_slope+1/20))
         error_f_inf = calculate_list_error(f_infinities, self.f_inf_values)
 
@@ -341,8 +379,33 @@ class Fitter:
                                             / abs(self.f_zero_slope_at_straight+1 / 10)
         error_f_zero = calculate_list_error(f_zeros, self.f_zero_values)
 
+        # TODO
+
+        if model.get_sampling_interval() != self.data_sampling_interval:
+            raise ValueError("Sampling intervals not the same!")
+
+        times, freqs = fi_curve_model.get_mean_time_and_freq_traces()
+        freq_prediction = np.array(freqs[self.f_zero_curve_contrast_idx])
+        time_prediction = np.array(times[self.f_zero_curve_contrast_idx])
+        stimulus_start = fi_curve_model.get_stimulus_start() - time_prediction[0]
+        length = fi_curve_model.get_stimulus_duration() / 2
+
+        start_idx = int(stimulus_start / fi_curve_model.get_sampling_interval())
+        end_idx = int((stimulus_start + length) / model.get_sampling_interval())
+
+        if len(time_prediction) == 0 or len(time_prediction) < end_idx or time_prediction[0] > fi_curve_model.get_stimulus_start():
+            error_f0_curve = 200
+        else:
+            error_f0_curve = np.mean((self.f_zero_curve_freq[start_idx:end_idx] - freq_prediction[start_idx:end_idx])**2) / 100
+
+        # plt.plot(self.f_zero_curve_freq[start_idx:end_idx])
+        # plt.plot(freq_prediction[start_idx:end_idx])
+        # plt.plot((self.f_zero_curve_freq[start_idx:end_idx] - freq_prediction[start_idx:end_idx])**2)
+        # plt.show()
+        # plt.close()
+
         error_list = [error_bf, error_vs, error_sc, error_cv, error_bursty,
-                      error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope_at_straight]
+                      error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope_at_straight, error_f0_curve]
 
         if error_weights is not None and len(error_weights) == len(error_list):
             for i in range(len(error_weights)):
@@ -368,6 +431,8 @@ def calculate_list_error(fit, reference):
 
     return norm_error
 
+def calculate_f0_curve_error(data_ficurve, model_ficurve):
+    return 0
 
 def normed_quadratic_freq_error(fit, ref, factor=2):
     return (abs(fit-ref)/factor)**2 / ref

helperFunctions.py

@@ -481,7 +481,6 @@ def detect_f_zero_in_frequency_trace(time, frequency, stimulus_start, sampling_i
     if start_idx < 0:
         raise ValueError("Time window to detect f_zero starts in an negative index!")
 
-
     min_during_start_of_stim = min(frequency[start_idx:end_idx])
     max_during_start_of_stim = max(frequency[start_idx:end_idx])