From 5af4117d4ca572eff4ee24f89ffedbb86dd428ac Mon Sep 17 00:00:00 2001 From: "a.ott" Date: Fri, 11 Sep 2020 17:21:47 +0200 Subject: [PATCH] add option to fit without dendTau or ref_period --- Fitter.py | 199 +++++++++++++++++++++++++++++++++++++++++-- models/LIFACnoise.py | 54 +++++++++++- run_Fitter.py | 31 +++---- 3 files changed, 256 insertions(+), 28 deletions(-) diff --git a/Fitter.py b/Fitter.py index 60f5119..4c914cb 100644 --- a/Fitter.py +++ b/Fitter.py @@ -136,9 +136,16 @@ 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") + # tau mins: + tau_min = 0.001 + for i in (0, 3, 5): + if X[i] < tau_min: + print("tried too small tau value") + return 1000 + abs(X[i] - tau_min) * 10000 + + for i in (1, 2, 4, 6): + if X[i] < 0: + print("tried negative parameter value") return 1000 + abs(X[i]) * 10000 if X[6] > 1.05/self.eod_freq: # refractory period shouldn't be larger than one eod period @@ -170,6 +177,184 @@ class Fitter: self.best_parameters_found = X return sum(error_list) + def fit_routine_no_dend_tau(self, start_parameters, error_weights=None): + self.counter = 0 + # fit only v_offset, mem_tau, input_scaling, dend_tau + + x0 = np.array([start_parameters["mem_tau"], start_parameters["noise_strength"], + start_parameters["input_scaling"], start_parameters["tau_a"], + start_parameters["delta_a"], start_parameters["refractory_period"]]) + initial_simplex = create_init_simples(x0, search_scale=3) + + # 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_f0_curve] + + 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": 400}) + + return fmin, self.base_model.get_parameters() + + def cost_function_no_dend_tau(self, X, error_weights=None): + # tau mins: + tau_min = 0.001 + for i in (0, 3): + if X[i] < tau_min: + print("tried too small tau value") + return 1000 + abs(X[i] - tau_min) * 10000 + + for i in (1, 2, 4, 5): + if X[i] < 0: + print("tried negative parameter value") + return 1000 + abs(X[i]) * 10000 + + if X[6] > 1.05/self.eod_freq: # refractory period shouldn't be larger than one eod period + print("tried too large ref period") + return 1000 + abs(X[6]) * 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]) + self.base_model.set_variable("tau_a", X[3]) + self.base_model.set_variable("delta_a", X[4]) + self.base_model.set_variable("refractory_period", 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) + + # time1 = time.time() + v_offset = test_model.find_v_offset(self.baseline_freq, base_stimulus) + self.base_model.set_variable("v_offset", v_offset) + # time2 = time.time() + # print("time taken for finding v_offset: {:.2f}s".format(time2-time1)) + + error_list = self.calculate_errors(error_weights) + # 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 fit_routine_no_ref_period(self, start_parameters, error_weights=None): + self.counter = 0 + # fit only v_offset, mem_tau, input_scaling, dend_tau + + self.base_model.set_variable("refractory_period", 0) + + x0 = np.array([start_parameters["mem_tau"], start_parameters["noise_strength"], + start_parameters["input_scaling"], start_parameters["tau_a"], start_parameters["delta_a"], + start_parameters["dend_tau"]]) + initial_simplex = create_init_simples(x0, search_scale=3) + + # 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_f0_curve] + + 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": 400}) + + return fmin, self.base_model.get_parameters() + + def cost_function_no_ref_period(self, X, error_weights=None): + # tau mins: + tau_min = 0.001 + for i in (0, 3, 5): + if X[i] < tau_min: + print("tried too small tau value") + return 1000 + abs(X[i] - tau_min) * 10000 + + for i in (1, 2, 4): + if X[i] < 0: + print("tried negative parameter value") + return 1000 + abs(X[i]) * 10000 + + if X[6] > 1.05/self.eod_freq: # refractory period shouldn't be larger than one eod period + print("tried too large ref period") + return 1000 + abs(X[6]) * 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]) + 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) + + # time1 = time.time() + v_offset = test_model.find_v_offset(self.baseline_freq, base_stimulus) + self.base_model.set_variable("v_offset", v_offset) + # time2 = time.time() + # print("time taken for finding v_offset: {:.2f}s".format(time2-time1)) + + error_list = self.calculate_errors(error_weights) + # 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 fit_routine_no_dend_tau_and_no_ref_period(self, start_parameters, error_weights=None): + self.counter = 0 + # fit only v_offset, mem_tau, input_scaling, dend_tau + self.base_model.set_variable("refractory_period", 0) + self.base_model.parameters["dend_tau"] = 0 + x0 = np.array([start_parameters["mem_tau"], start_parameters["noise_strength"], + start_parameters["input_scaling"], start_parameters["tau_a"], start_parameters["delta_a"]]) + initial_simplex = create_init_simples(x0, search_scale=3) + + # 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_f0_curve] + + 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": 400}) + + return fmin, self.base_model.get_parameters() + + def cost_function_no_dend_tau_and_no_ref_period(self, X, error_weights=None): + # tau mins: + tau_min = 0.001 + for i in (0, 3): + if X[i] < tau_min: + print("tried too small tau value") + return 1000 + abs(X[i] - tau_min) * 10000 + + for i in (1, 2, 4): + if X[i] < 0: + print("tried negative parameter value") + return 1000 + abs(X[i]) * 10000 + + if X[6] > 1.05/self.eod_freq: # refractory period shouldn't be larger than one eod period + print("tried too large ref period") + return 1000 + abs(X[6]) * 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]) + 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) + + # time1 = time.time() + v_offset = test_model.find_v_offset(self.baseline_freq, base_stimulus) + self.base_model.set_variable("v_offset", v_offset) + # time2 = time.time() + # print("time taken for finding v_offset: {:.2f}s".format(time2-time1)) + + error_list = self.calculate_errors(error_weights) + # 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 calculate_errors(self, error_weights=None, model=None): if model is None: model = self.base_model @@ -202,7 +387,7 @@ class Fitter: error_vs = abs((vector_strength - self.vector_strength) / 0.01) error_cv = abs((coefficient_of_variation - self.coefficient_of_variation) / 0.05) error_bursty = (abs(burstiness - self.burstiness) / 0.2) - error_hist = np.mean((isi_bins - self.isi_bins) ** 2) / 600 + error_hist = np.sqrt(np.mean((isi_bins - self.isi_bins) ** 2)) / 10 # print("error hist: {:.2f}".format(error_hist)) # print("Burstiness: cell {:.2f}, model: {:.2f}, error: {:.2f}".format(self.burstiness, burstiness, error_bursty)) @@ -211,7 +396,7 @@ class Fitter: error_sc += abs((serial_correlation[i] - self.serial_correlation[i]) / 0.1) # 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_slope = abs((f_infinities_slope - self.f_inf_slope) / abs(self.f_inf_slope+1)) * 25 error_f_inf = calculate_list_error(f_infinities, self.f_inf_values) # error_f_zero_slopes = calculate_list_error(f_zero_slopes, self.f_zero_slopes) @@ -219,7 +404,7 @@ class Fitter: / abs(self.f_zero_slope_at_straight+1) error_f_zero = calculate_list_error(f_zeros, self.f_zero_values) / 10 - error_f0_curve = self.calculate_f0_curve_error(model, fi_curve_model) / 1000 + error_f0_curve = self.calculate_f0_curve_error(model, fi_curve_model) / 20 error_list = [error_vs, error_sc, error_cv, error_hist, error_bursty, error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope_at_straight, error_f0_curve] @@ -290,7 +475,7 @@ class Fitter: elif len(model_curve) < len(cell_curve): cell_curve = cell_curve[:len(model_curve)] - error_f0_curve = np.mean((model_curve - cell_curve) ** 2) + error_f0_curve = np.sqrt(np.mean((model_curve - cell_curve) ** 2)) return error_f0_curve diff --git a/models/LIFACnoise.py b/models/LIFACnoise.py index 28a32aa..544fc2e 100644 --- a/models/LIFACnoise.py +++ b/models/LIFACnoise.py @@ -128,8 +128,10 @@ class LifacNoiseModel(AbstractModel): parameters = np.array( [v_zero, a_zero, step_size, threshold, v_base, delta_a, tau_a, v_offset, mem_tau, noise_strength, time_start, input_scaling, dend_tau, ref_period]) - - voltage_trace, adaption, spiketimes, input_voltage = simulate_fast(rectified_stimulus, total_time_s, parameters) + if dend_tau >= step_size: + voltage_trace, adaption, spiketimes, input_voltage = simulate_fast(rectified_stimulus, total_time_s, parameters) + else: + voltage_trace, adaption, spiketimes, input_voltage = simulate_fast_no_dend_tau(rectified_stimulus, total_time_s, parameters) self.stimulus = stimulus self.input_voltage = input_voltage @@ -304,3 +306,51 @@ def simulate_fast(rectified_stimulus_array, total_time_s, parameters: np.ndarray adaption[i] += delta_a / tau_a return output_voltage, adaption, spiketimes, input_voltage + + +@jit(nopython=True) +def simulate_fast_no_dend_tau(rectified_stimulus_array, total_time_s, parameters: np.ndarray): + v_zero = parameters[0] + a_zero = parameters[1] + step_size = parameters[2] + threshold = parameters[3] + v_base = parameters[4] + delta_a = parameters[5] + tau_a = parameters[6] + v_offset = parameters[7] + mem_tau = parameters[8] + noise_strength = parameters[9] + time_start = parameters[10] + input_scaling = parameters[11] + dend_tau = parameters[12] + ref_period = parameters[13] + + time = np.arange(time_start, total_time_s, step_size) + length = len(time) + output_voltage = np.zeros(length) + adaption = np.zeros(length) + input_voltage = rectified_stimulus_array + + spiketimes = [] + output_voltage[0] = v_zero + adaption[0] = a_zero + + for i in range(1, len(time), 1): + + noise_value = np.random.normal() + noise = noise_strength * noise_value / np.sqrt(step_size) + + output_voltage[i] = output_voltage[i - 1] + ((v_base - output_voltage[i - 1] + v_offset + ( + input_voltage[i] * input_scaling) - adaption[i - 1] + noise) / mem_tau) * step_size + + adaption[i] = adaption[i - 1] + ((-adaption[i - 1]) / tau_a) * step_size + + if len(spiketimes) > 0 and time[i] - spiketimes[-1] < ref_period + step_size/2: + output_voltage[i] = v_base + + if output_voltage[i] > threshold: + output_voltage[i] = v_base + spiketimes.append((i * step_size) + time_start) + adaption[i] += delta_a / tau_a + + return output_voltage, adaption, spiketimes, input_voltage \ No newline at end of file diff --git a/run_Fitter.py b/run_Fitter.py index 2afa107..8667052 100644 --- a/run_Fitter.py +++ b/run_Fitter.py @@ -4,7 +4,7 @@ from CellData import CellData from Baseline import get_baseline_class from FiCurve import get_fi_curve_class from Fitter import Fitter -from ModelFit import ModelFit +from ModelFit import get_best_fit, ModelFit import time import os @@ -14,8 +14,8 @@ from helperFunctions import plot_errors import multiprocessing as mp -SAVE_DIRECTORY = "./results/final_2/" -SAVE_DIRECTORY_BEST = "./results/final_2_best/" +SAVE_DIRECTORY = "./results/final_3/" +SAVE_DIRECTORY_BEST = "./results/final_3_best/" # [vs, sc, cv, isi_hist, bursty, f_inf, f_inf_slope, f_zero, f_zero_slope, f0_curve] ERROR_WEIGHTS = (1, 1, 1, 1, 1, 1, 1, 1, 0, 1) @@ -45,7 +45,7 @@ def test_single_cell(path): for i, p in enumerate(start_parameters): fitter = Fitter() fitter.set_data_reference_values(cell_data) - fmin, res_par = fitter.fit_routine(p, ERROR_WEIGHTS) + fmin, res_par = fitter.fit_routine_no_ref_period(p, ERROR_WEIGHTS) cell_path = os.path.split(cell_data.get_data_path())[-1] @@ -63,17 +63,19 @@ def fit_cell_base(parameters): fmin, res_par = fitter.fit_routine(parameters[2], ERROR_WEIGHTS) cell_data = parameters[0] - cell_path = os.path.split(cell_data.get_data_path())[-1] + cell_name = os.path.split(cell_data.get_data_path())[-1] error = fitter.calculate_errors(model=LifacNoiseModel(res_par)) - save_path = SAVE_DIRECTORY + "/" + cell_path + "/start_parameter_{:}_err_{:.2f}/".format(parameters[1], sum(error)) + save_path = SAVE_DIRECTORY + "/" + cell_name + "/start_parameter_{:}_err_{:.2f}/".format(parameters[1], sum(error)) save_fitting_run_info(parameters[0], res_par, parameters[2], plot=True, save_path=save_path) plot_errors(fitter.errors, save_path) + fit = ModelFit(save_path) + fit.generate_master_plot(save_path) time2 = time.time() del fitter - print("Time taken for " + cell_path + + print("Time taken for " + cell_name + "\n and start parameters ({:}): {:.2f}s thread time".format(parameters[1]+1, time2 - time1) + "\n error: {:.2f}".format(sum(error))) @@ -97,19 +99,10 @@ def fit_cell_parallel(cell_data, start_parameters): del pool del cell_data - save_master_plot(save_path_cell) - - -def save_master_plot(save_path_cell): - best_fit = None - min_err = np.inf - for fit in os.listdir(save_path_cell): - cur_fit = ModelFit(os.path.join(save_path_cell, fit)) - if cur_fit.comparable_error() < min_err: - min_err = cur_fit.comparable_error() - best_fit = cur_fit - + best_fit = get_best_fit(save_path_cell) best_fit.generate_master_plot(SAVE_DIRECTORY_BEST) + best_fit.generate_master_plot(SAVE_DIRECTORY) + def iget_start_parameters():