diff --git a/.gitignore b/.gitignore index 0bb954b..97de299 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ /venv/ __pycache__/ .idea/ +/results/ diff --git a/Fitter.py b/Fitter.py new file mode 100644 index 0000000..82cca69 --- /dev/null +++ b/Fitter.py @@ -0,0 +1,344 @@ + +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()