from models.LIFACnoise import LifacNoiseModel from stimuli.SinusoidalStepStimulus import SinusoidalStepStimulus from CellData import CellData, icelldata_of_dir from FiCurve import FICurve from AdaptionCurrent import Adaption import helperFunctions as hF import functions as fu import numpy as np from warnings import warn from scipy.optimize import minimize import time import os SAVE_PATH_PREFIX = "" def main(): run_with_real_data() def iget_start_parameters(mem_tau_list=None, input_scaling_list=None, noise_strength_list=None, dend_tau_list=None, tau_a_list=None, delta_a_list=None): # mem_tau, input_scaling, noise_strength, dend_tau, # expand by tau_a, delta_a ? if mem_tau_list is None: mem_tau_list = [0.01] if input_scaling_list is None: input_scaling_list = [40, 60, 80] if noise_strength_list is None: noise_strength_list = [0.03] # [0.02, 0.06] if dend_tau_list is None: dend_tau_list = [0.001, 0.002] # if tau_a_list is None: # tau_a_list = # if delta_a_list is None: # delta_a_list = for mem_tau in mem_tau_list: for input_scaling in input_scaling_list: for noise_strength in noise_strength_list: for dend_tau in dend_tau_list: yield {"mem_tau": mem_tau, "input_scaling": input_scaling, "noise_strength": noise_strength, "dend_tau": dend_tau} def run_with_real_data(): for cell_data in icelldata_of_dir("./data/"): start_par_count = 0 for start_parameters in iget_start_parameters(): start_par_count += 1 print("START PARAMETERS:", start_par_count) if start_par_count <= 0: continue print("cell:", cell_data.get_data_path()) trace = cell_data.get_base_traces(trace_type=cell_data.V1) if len(trace) == 0: print("NO V1 TRACE FOUND") continue results_path = "results/" + os.path.split(cell_data.get_data_path())[-1] + "/" print("results at:", results_path) start_time = time.time() fitter = Fitter() fmin, parameters = fitter.fit_model_to_data(cell_data, start_parameters) print(fmin) print(parameters) end_time = time.time() if not os.path.exists(results_path): os.makedirs(results_path) with open(results_path + "fit_parameters_start_{}.txt".format(start_par_count), "w") as file: file.writelines(["start_parameters:\t" + str(start_parameters), "\nfinal_parameters:\t" + str(parameters), "\nfinal_fmin:\t" + str(fmin)]) results_path += SAVE_PATH_PREFIX + "par_set_" + str(start_par_count) + "_" print('Fitting of cell took function took {:.3f} s'.format((end_time - start_time))) #print(results_path) print_comparision_cell_model(cell_data, parameters, plot=True, savepath=results_path) break # from Sounds import play_finished_sound # play_finished_sound() pass def print_comparision_cell_model(cell_data, parameters, plot=False, savepath=None): res_model = LifacNoiseModel(parameters) fi_curve = FICurve(cell_data) m_bf, m_vs, m_sc, m_cv = res_model.calculate_baseline_markers(cell_data.get_eod_frequency()) f_baselines, f_zeros, m_f_infinities = res_model.calculate_fi_curve(fi_curve.stimulus_value, cell_data.get_eod_frequency()) f_infinities_fit = hF.fit_clipped_line(fi_curve.stimulus_value, m_f_infinities) m_f_infinities_slope = f_infinities_fit[0] f_zeros_fit = hF.fit_boltzmann(fi_curve.stimulus_value, f_zeros) m_f_zero_slope = fu.full_boltzmann_straight_slope(f_zeros_fit[0], f_zeros_fit[1], f_zeros_fit[2], f_zeros_fit[3]) c_bf = cell_data.get_base_frequency() c_vs = cell_data.get_vector_strength() c_sc = cell_data.get_serial_correlation(1) c_cv = cell_data.get_coefficient_of_variation() c_f_slope = fi_curve.get_f_infinity_slope() c_f_values = fi_curve.f_infinities c_f_zero_slope = fi_curve.get_fi_curve_slope_of_straight() c_f_zero_values = fi_curve.f_zeros print("bf: cell - {:.2f} vs model {:.2f}".format(c_bf, m_bf)) print("vs: cell - {:.2f} vs model {:.2f}".format(c_vs, m_vs)) print("sc: cell - {:.2f} vs model {:.2f}".format(c_sc[0], m_sc[0])) print("cv: cell - {:.2f} vs model {:.2f}".format(c_cv, m_cv)) print("f_inf_slope: cell - {:.2f} vs model {:.2f}".format(c_f_slope, m_f_infinities_slope)) print("f infinity values:\n cell -", c_f_values, "\n model -", m_f_infinities) print("f_zero_slope: cell - {:.2f} vs model {:.2f}".format(c_f_zero_slope, m_f_zero_slope)) print("f zero values:\n cell -", c_f_zero_values, "\n model -", f_zeros) if plot: f_b, f_zero, f_inf = res_model.calculate_fi_curve(cell_data.get_fi_contrasts(), cell_data.get_eod_frequency()) fi_curve.plot_fi_curve(savepath=savepath, comp_f_baselines=f_b, comp_f_zeros=f_zero, comp_f_infs=f_inf) class Fitter: def __init__(self, params=None): if params is None: self.base_model = LifacNoiseModel({"step_size": 0.00005}) else: self.base_model = LifacNoiseModel(params) if "step_size" not in params: self.base_model.set_variable("step_size", 0.00005) # self.fi_contrasts = [] self.eod_freq = 0 self.sc_max_lag = 1 # values to be replicated: self.baseline_freq = 0 self.vector_strength = -1 self.serial_correlation = [] self.coefficient_of_variation = 0 self.f_inf_values = [] self.f_inf_slope = 0 self.f_zero_values = [] self.f_zero_slope = 0 self.f_zero_fit = [] self.tau_a = 0 self.delta_a = 0 # counts how often the cost_function was called self.counter = 0 def fit_model_to_data(self, data: CellData, start_parameters=None): self.eod_freq = data.get_eod_frequency() self.baseline_freq = data.get_base_frequency() self.vector_strength = data.get_vector_strength() self.serial_correlation = data.get_serial_correlation(self.sc_max_lag) self.coefficient_of_variation = data.get_coefficient_of_variation() fi_curve = FICurve(data, contrast=True) self.fi_contrasts = fi_curve.stimulus_value self.f_inf_values = fi_curve.f_infinities self.f_inf_slope = fi_curve.get_f_infinity_slope() self.f_zero_values = fi_curve.f_zeros self.f_zero_fit = fi_curve.boltzmann_fit_vars self.f_zero_slope = fi_curve.get_fi_curve_slope_of_straight() # self.f_zero_slope = fi_curve.get_fi_curve_slope_at(fi_curve.get_f_zero_and_f_inf_intersection()) # around 1/3 of the value at straight self.delta_a = (self.f_zero_slope / self.f_inf_slope) / 1000 # seems to work if divided by 1000... adaption = Adaption(data, fi_curve) self.tau_a = adaption.get_tau_real() # print("delta_a: {:.3f}".format(self.delta_a), "tau_a: {:.3f}".format(self.tau_a)) return self.fit_routine_5(data, start_parameters) # return self.fit_model(fit_adaption=False) def fit_routine_1(self, cell_data=None): global SAVE_PATH_PREFIX SAVE_PATH_PREFIX = "fit_routine_1_" # errors: [error_bf, error_vs, error_sc, error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope] self.counter = 0 # fit only v_offset, mem_tau, noise_strength, input_scaling x0 = np.array([0.02, 0.03, 70]) initial_simplex = create_init_simples(x0, search_scale=2) error_weights = (1, 1, 1, 1, 1, 0, 0) fmin_step1 = minimize(fun=self.cost_function_with_fixed_adaption, args=(self.tau_a, self.delta_a, error_weights), x0=x0, method="Nelder-Mead", options={"initial_simplex": initial_simplex}) res_parameters_step1 = self.base_model.get_parameters() if cell_data is not None: print("##### After step 1: (fixed adaption)") print_comparision_cell_model(cell_data, res_parameters_step1) self.counter = 0 x0 = np.array([res_parameters_step1["mem_tau"], res_parameters_step1["noise_strength"], res_parameters_step1["input_scaling"], res_parameters_step1["tau_a"], res_parameters_step1["delta_a"]]) initial_simplex = create_init_simples(x0, search_scale=2) error_weights = (1, 1, 1, 1, 1, 2, 4) fmin_step2 = minimize(fun=self.cost_function_all, args=(error_weights), x0=x0, method="Nelder-Mead", options={"initial_simplex": initial_simplex}) res_parameters_step2 = self.base_model.get_parameters() if cell_data is not None: print("##### After step 2: (Everything)") # print_comparision_cell_model(cell_data, res_parameters_step2) return fmin_step2, res_parameters_step2 def fit_routine_2(self, cell_data=None): global SAVE_PATH_PREFIX SAVE_PATH_PREFIX = "fit_routine_2_" # errors: [error_bf, error_vs, error_sc, error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope] self.counter = 0 # fit only v_offset, mem_tau, noise_strength, input_scaling x0 = np.array([0.02, 0.03, 70]) initial_simplex = create_init_simples(x0, search_scale=2) error_weights = (1, 1, 5, 1, 2, 0, 0) fmin = minimize(fun=self.cost_function_with_fixed_adaption, args=(self.tau_a, self.delta_a, error_weights), x0=x0, method="Nelder-Mead", options={"initial_simplex": initial_simplex}) res_parameters = self.base_model.get_parameters() return fmin, res_parameters def fit_routine_3(self, cell_data=None): global SAVE_PATH_PREFIX SAVE_PATH_PREFIX = "fit_routine_3_" # errors: [error_bf, error_vs, error_sc, error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope] self.counter = 0 # fit only v_offset, mem_tau, noise_strength, input_scaling, dend_tau x0 = np.array([0.02, 0.03, 70, 0.001]) initial_simplex = create_init_simples(x0, search_scale=2) error_weights = (1, 1, 5, 1, 2, 0, 0) fmin = minimize(fun=self.cost_function_with_fixed_adaption_with_dend_tau, args=(self.tau_a, self.delta_a, error_weights), x0=x0, method="Nelder-Mead", options={"initial_simplex": initial_simplex}) res_parameters = self.base_model.get_parameters() return fmin, res_parameters def fit_routine_4(self, cell_data=None, start_parameters=None): global SAVE_PATH_PREFIX SAVE_PATH_PREFIX = "fit_routine_4_" # errors: [error_bf, error_vs, error_sc, error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope] self.counter = 0 # fit only v_offset, mem_tau, input_scaling, dend_tau if start_parameters is None: x0 = np.array([0.02, 70, 0.001]) else: x0 = np.array([start_parameters["mem_tau"], start_parameters["noise_strength"], start_parameters["input_scaling"], start_parameters["dend_tau"]]) initial_simplex = create_init_simples(x0, search_scale=2) error_weights = (0, 5, 15, 1, 2, 1, 0) fmin = minimize(fun=self.cost_function_with_fixed_adaption_with_dend_tau, args=(self.tau_a, self.delta_a, error_weights), x0=x0, method="Nelder-Mead", options={"initial_simplex": initial_simplex, "xatol": 0.001, "maxfev": 400, "maxiter": 400}) res_parameters = fmin.x # print_comparision_cell_model(cell_data, self.base_model.get_parameters()) self.counter = 0 x0 = np.array([self.tau_a, self.delta_a, res_parameters[0]]) initial_simplex = create_init_simples(x0, search_scale=2) error_weights = (0, 1, 1, 2, 2, 4, 2) fmin = minimize(fun=self.cost_function_only_adaption, args=(error_weights,), x0=x0, method="Nelder-Mead", options={"initial_simplex": initial_simplex, "xatol": 0.001}) res_parameters = fmin.x print(fmin) print_comparision_cell_model(cell_data, self.base_model.get_parameters()) # # # self.counter = 0 # # x0 = np.array([res_parameters[0], # # res_parameters[1], self.tau_a, # # self.delta_a, res_parameters[2]]) # # initial_simplex = create_init_simples(x0, search_scale=2) # # error_weights = (1, 3, 1, 2, 1, 3, 2) # # fmin = minimize(fun=self.cost_function_all_without_noise, # # args=(error_weights,), x0=x0, method="Nelder-Mead", # # options={"initial_simplex": initial_simplex, "xatol": 0.001}) # # res_parameters = self.base_model.get_parameters() # # # # print_comparision_cell_model(cell_data, self.base_model.get_parameters()) # # self.counter = 0 # x0 = np.array([res_parameters[0], start_parameters["noise_strength"], # res_parameters[1], res_parameters[2], # res_parameters[3], res_parameters[4]]) # initial_simplex = create_init_simples(x0, search_scale=2) # error_weights = (0, 1, 2, 1, 1, 3, 2) # fmin = minimize(fun=self.cost_function_all, # args=(error_weights,), x0=x0, method="Nelder-Mead", # options={"initial_simplex": initial_simplex, "xatol": 0.001, "maxiter": 599}) # res_parameters = self.base_model.get_parameters() return fmin, self.base_model.get_parameters() def fit_routine_5(self, cell_data=None, start_parameters=None): global SAVE_PATH_PREFIX SAVE_PATH_PREFIX = "fit_routine_5_" # errors: [error_bf, error_vs, error_sc, error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope] self.counter = 0 # fit only v_offset, mem_tau, input_scaling, dend_tau if start_parameters is None: x0 = np.array([0.02, 70, 0.001]) else: x0 = np.array([start_parameters["mem_tau"], start_parameters["noise_strength"], start_parameters["input_scaling"], self.tau_a, self.delta_a, start_parameters["dend_tau"]]) initial_simplex = create_init_simples(x0, search_scale=2) error_weights = (0, 1, 1, 1, 1, 1, 2, 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": 400, "maxiter": 400}) res_parameters = fmin.x return fmin, self.base_model.get_parameters() def fit_model(self, x0=None, initial_simplex=None, fit_adaption=False): self.counter = 0 if fit_adaption: if x0 is None: x0 = np.array([0.02, 0.03, 70, self.tau_a, self.delta_a]) if initial_simplex is None: initial_simplex = create_init_simples(x0) fmin = minimize(fun=self.cost_function_all, x0=x0, method="Nelder-Mead", options={"initial_simplex": initial_simplex}) else: if x0 is None: x0 = np.array([0.02, 0.03, 70]) if initial_simplex is None: initial_simplex = create_init_simples(x0) fmin = minimize(fun=self.cost_function_with_fixed_adaption, x0=x0, args=(self.tau_a, self.delta_a), method="Nelder-Mead", options={"initial_simplex": initial_simplex}) return fmin, self.base_model.get_parameters() def cost_function_all(self, X, error_weights=None): 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]) self.base_model.set_variable("tau_a", X[3]) self.base_model.set_variable("delta_a", X[4]) self.base_model.set_variable("dend_tau", X[5]) base_stimulus = SinusoidalStepStimulus(self.eod_freq, 0) # find right v-offset test_model = self.base_model.get_model_copy() test_model.set_variable("noise_strength", 0) v_offset = test_model.find_v_offset(self.baseline_freq, base_stimulus) self.base_model.set_variable("v_offset", v_offset) # [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) return sum(error_list) def cost_function_all_without_noise(self, X, error_weights=None): self.base_model.set_variable("mem_tau", X[0]) self.base_model.set_variable("input_scaling", X[1]) self.base_model.set_variable("tau_a", X[2]) self.base_model.set_variable("delta_a", X[3]) self.base_model.set_variable("dend_tau", X[4]) self.base_model.set_variable("noise_strength", 0) base_stimulus = SinusoidalStepStimulus(self.eod_freq, 0) # find right v-offset test_model = self.base_model.get_model_copy() test_model.set_variable("noise_strength", 0) v_offset = test_model.find_v_offset(self.baseline_freq, base_stimulus) self.base_model.set_variable("v_offset", v_offset) # [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) return sum(error_list) def cost_function_only_adaption(self, X, error_weights=None): self.base_model.set_variable("tau_a", X[0]) self.base_model.set_variable("delta_a", X[1]) self.base_model.set_variable("mem_tau", X[2]) base_stimulus = SinusoidalStepStimulus(self.eod_freq, 0) # find right v-offset test_model = self.base_model.get_model_copy() test_model.set_variable("noise_strength", 0) v_offset = test_model.find_v_offset(self.baseline_freq, base_stimulus) self.base_model.set_variable("v_offset", v_offset) # [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) return sum(error_list) def cost_function_with_fixed_adaption(self, X, tau_a, delta_a, error_weights=None): # set model parameters: model = self.base_model model.set_variable("mem_tau", X[0]) model.set_variable("noise_strength", X[1]) model.set_variable("input_scaling", X[2]) model.set_variable("tau_a", tau_a) model.set_variable("delta_a", delta_a) base_stimulus = SinusoidalStepStimulus(self.eod_freq, 0) # find right v-offset test_model = model.get_model_copy() test_model.set_variable("noise_strength", 0) v_offset = test_model.find_v_offset(self.baseline_freq, base_stimulus) model.set_variable("v_offset", v_offset) error_list = self.calculate_errors(error_weights) return sum(error_list) def cost_function_with_fixed_adaption_with_dend_tau_no_noise(self, X, tau_a, delta_a, error_weights=None): # set model parameters: model = self.base_model model.set_variable("mem_tau", X[0]) model.set_variable("input_scaling", X[1]) model.set_variable("dend_tau", X[2]) model.set_variable("tau_a", tau_a) model.set_variable("delta_a", delta_a) model.set_variable("noise_strength", 0) base_stimulus = SinusoidalStepStimulus(self.eod_freq, 0) # find right v-offset test_model = model.get_model_copy() test_model.set_variable("noise_strength", 0) v_offset = test_model.find_v_offset(self.baseline_freq, base_stimulus) model.set_variable("v_offset", v_offset) error_list = self.calculate_errors(error_weights) return sum(error_list) def cost_function_with_fixed_adaption_with_dend_tau(self, X, tau_a, delta_a, error_weights=None): # set model parameters: model = self.base_model model.set_variable("mem_tau", X[0]) model.set_variable("noise_strength", X[1]) model.set_variable("input_scaling", X[2]) model.set_variable("dend_tau", X[3]) model.set_variable("tau_a", tau_a) model.set_variable("delta_a", delta_a) base_stimulus = SinusoidalStepStimulus(self.eod_freq, 0) # find right v-offset test_model = model.get_model_copy() test_model.set_variable("noise_strength", 0) v_offset = test_model.find_v_offset(self.baseline_freq, base_stimulus) model.set_variable("v_offset", v_offset) error_list = self.calculate_errors(error_weights) return sum(error_list) def calculate_errors(self, error_weights=None): baseline_freq, vector_strength, serial_correlation, coefficient_of_variation =\ self.base_model.calculate_baseline_markers(self.eod_freq, self.sc_max_lag) # print("baseline features calculated!") # f_infinities, f_infinities_slope = self.base_model.calculate_fi_markers(self.fi_contrasts, self.eod_freq) f_baselines, f_zeros, f_infinities = self.base_model.calculate_fi_curve(self.fi_contrasts, self.eod_freq) try: f_infinities_fit = hF.fit_clipped_line(self.fi_contrasts, f_infinities) except Exception as e: print("EXCEPTION IN FIT LINE!") print(e) f_infinities_fit = [0, 0] f_infinities_slope = f_infinities_fit[0] try: f_zeros_fit = hF.fit_boltzmann(self.fi_contrasts, f_zeros) except Exception as e: print("EXCEPTION IN FIT BOLTZMANN!") print(e) f_zeros_fit = [0, 0, 0, 0] f_zero_slope = fu.full_boltzmann_straight_slope(f_zeros_fit[0], f_zeros_fit[1], f_zeros_fit[2], f_zeros_fit[3]) # print("fi-curve features calculated!") # calculate errors with reference values error_bf = abs((baseline_freq - self.baseline_freq) / self.baseline_freq) error_vs = abs((vector_strength - self.vector_strength) / self.vector_strength) error_sc = abs((serial_correlation[0] - self.serial_correlation[0]) / self.serial_correlation[0]) error_cv = abs((coefficient_of_variation - self.coefficient_of_variation) / self.coefficient_of_variation) error_f_inf_slope = abs((f_infinities_slope - self.f_inf_slope) / self.f_inf_slope) * 4 error_f_inf = calculate_f_values_error(f_infinities, self.f_inf_values) * .5 error_f_zero_slope = abs((f_zero_slope - self.f_zero_slope) / self.f_zero_slope) error_f_zero = calculate_f_values_error(f_zeros, self.f_zero_values) error_list = [error_bf, error_vs, error_sc, error_cv, error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope] if error_weights is not None and len(error_weights) == len(error_list): for i in range(len(error_weights)): error_list[i] = error_list[i] * error_weights[i] elif len(error_weights) != len(error_list): warn("Error weights had different length than errors and were ignored!") error = sum(error_list) self.counter += 1 if self.counter % 200 == 0: # and False: # TODO currently shut off! print("\nCost function run times: {:}\n".format(self.counter), "Total weighted error: {:.4f}\n".format(error), "Baseline frequency - expected: {:.0f}, current: {:.0f}, error: {:.3f}\n".format( self.baseline_freq, baseline_freq, error_bf), "Vector strength - expected: {:.2f}, current: {:.2f}, error: {:.3f}\n".format( self.vector_strength, vector_strength, error_vs), "Serial correlation - expected: {:.2f}, current: {:.2f}, error: {:.3f}\n".format( self.serial_correlation[0], serial_correlation[0], error_sc), "f-infinity slope - expected: {:.0f}, current: {:.0f}, error: {:.3f}\n".format( self.f_inf_slope, f_infinities_slope, error_f_inf_slope), "f-infinity values:\nexpected:", np.around(self.f_inf_values), "\ncurrent: ", np.around(f_infinities), "\nerror: {:.3f}\n".format(error_f_inf), "f-zero slope - expected: {:.0f}, current: {:.0f}, error: {:.3f}\n".format( self.f_zero_slope, f_zero_slope, error_f_zero_slope), "f-zero values:\nexpected:", np.around(self.f_zero_values), "\ncurrent: ", np.around(f_zeros), "\nerror: {:.3f}".format(error_f_zero)) return error_list def calculate_f_values_error(fit, reference): error = 0 for i in range(len(reference)): # TODO ??? add a constant to f_inf to allow for small differences in small values # example: 1 vs 3 would result in way smaller error. constant = 10 error += abs((fit[i] - reference[i])+constant) / (abs(reference[i]) + constant) norm_error = error / len(reference) return norm_error def create_init_simples(x0, search_scale=3.): dim = len(x0) simplex = [[x0[0]/search_scale], [x0[0]*search_scale]] for i in range(1, dim, 1): for vertex in simplex: vertex.append(x0[i]*search_scale) new_vertex = list(x0[:i]) new_vertex.append(x0[i]/search_scale) simplex.append(new_vertex) return simplex if __name__ == '__main__': main()