diff --git a/Fitter.py b/Fitter.py index c097010..f13d19e 100644 --- a/Fitter.py +++ b/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 diff --git a/helperFunctions.py b/helperFunctions.py index ee28b6c..ba6c884 100644 --- a/helperFunctions.py +++ b/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])