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 scipy.optimize import minimize import time import os def main(): run_with_real_data() def run_with_real_data(): count = 0 for cell_data in icelldata_of_dir("./data/"): count += 1 if count <= 3: 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) 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.txt", "w") as file: file.writelines([str(parameters)]) results_path += "fit_routine_2_" print('Fitting of cell took function took {:.3f} s'.format((end_time - start_time))) print_comparision_cell_model(cell_data, parameters, plot=True, savepath=results_path) pass def print_comparision_cell_model(cell_data, parameters, plot=False, savepath=None): res_model = LifacNoiseModel(parameters) m_bf, m_vs, m_sc = res_model.calculate_baseline_markers(cell_data.get_eod_frequency()) m_f_values, m_f_slope = res_model.calculate_fi_markers(cell_data.get_fi_contrasts(), cell_data.get_eod_frequency()) c_bf = cell_data.get_base_frequency() c_vs = cell_data.get_vector_strength() c_sc = cell_data.get_serial_correlation(1) fi_curve = FICurve(cell_data) c_f_slope = fi_curve.get_f_infinity_slope() c_f_values = fi_curve.f_infinities 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("f_slope: cell - {:.2f} vs model {:.2f}".format(c_f_slope, m_f_slope)) print("f values:\n cell -", c_f_values, "\n model -", m_f_values) 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.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): 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) 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.delta_a = (self.f_zero_slope / self.f_inf_slope) / 1000 adaption = Adaption(data, fi_curve) self.tau_a = adaption.get_tau_real() return self.fit_routine_2(data) # return self.fit_model(fit_adaption=False) def fit_routine_1(self, cell_data=None): # 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): # 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_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]) 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_and_v_offset(self, X, error_weights=None): self.base_model.set_variable("tau_a", X[0]) self.base_model.set_variable("delta_a", X[1]) 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]) 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 calculate_errors(self, error_weights=None): baseline_freq, vector_strength, serial_correlation = 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) f_infinities_fit = hF.fit_clipped_line(self.fi_contrasts, f_infinities) f_infinities_slope = f_infinities_fit[0] f_zeros_fit = hF.fit_boltzmann(self.fi_contrasts, f_zeros) 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_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_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] error = sum(error_list) self.counter += 1 if self.counter % 200 == 0: 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 = 0 error += abs((fit[i] - reference[i]) / (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()