diff --git a/Fitter.py b/Fitter.py index 1715da0..c0b96c1 100644 --- a/Fitter.py +++ b/Fitter.py @@ -9,6 +9,7 @@ import numpy as np from warnings import warn from scipy.optimize import minimize import time +from helperFunctions import plot_errors import matplotlib.pyplot as plt @@ -53,6 +54,7 @@ class Fitter: self.f_zero_curve_freq = np.array([]) self.f_zero_curve_time = np.array([]) + self.errors = [] # self.tau_a = 0 @@ -125,54 +127,9 @@ class Fitter: 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": 800}) + options={"initial_simplex": initial_simplex, "xatol": 0.001, "maxfev": 600, "maxiter": 400}) - print("best that was returned: {}".format(fmin.fun)) - print("best that was visited: {}".format(self.smallest_error)) - - return fmin, self.base_model.get_parameters() - - def fit_routine_const_ref_period(self, start_parameters): - 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"], self.tau_a, start_parameters["delta_a"], - start_parameters["dend_tau"], start_parameters["refractory_period"]]) - initial_simplex = create_init_simples(x0, search_scale=2) - self.base_model.set_variable("refractory_period", start_parameters["refractory_period"]) - # 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, 1, 1, 1, 1, 1, 1, 1, 1) - fmin = minimize(fun=self.cost_function_without_ref_period, - args=(error_weights,), x0=x0, method="Nelder-Mead", - options={"initial_simplex": initial_simplex, "xatol": 0.001, "maxfev": 200, "maxiter": 400}) - - return fmin, self.base_model.get_parameters() - - # similar results to fit routine 1 - def fit_routine_2(self, start_parameters): - self.counter = 0 - - x0 = np.array([start_parameters["mem_tau"], start_parameters["input_scaling"], # mem_tau, input_scaling - start_parameters["delta_a"], start_parameters["dend_tau"]]) # delta_a, dend_tau - initial_simplex = create_init_simples(x0, search_scale=2) - - error_weights = (0, 1, 1, 1, 1, 3, 2, 3, 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, "maxfev": 100, "maxiter": 400}) - best_pars = fmin.x - x0 = np.array([best_pars[0], start_parameters["noise_strength"], - best_pars[1], self.tau_a, best_pars[2], - best_pars[3]]) - initial_simplex = create_init_simples(x0, search_scale=2) - # 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, 2, 1, 1, 1, 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": 100, "maxiter": 400}) + plot_errors(self.errors) return fmin, self.base_model.get_parameters() @@ -206,7 +163,7 @@ class Fitter: if sum(error_list) < self.smallest_error: self.smallest_error = sum(error_list) self.best_parameters_found = X - + self.errors.append(error_list) return sum(error_list) def cost_function_without_ref_period(self, X, error_weights=None): diff --git a/analysis.py b/analysis.py index 6356167..b3a7076 100644 --- a/analysis.py +++ b/analysis.py @@ -17,28 +17,28 @@ def main(): # parser.add_argument("dir", help="folder containing the cell folders with the fit results") # args = parser.parse_args() - # dir_path = "results/invivo_results/" # args.dir - dir_path = "results/results_add__trial_more_iter_NM/invivo_results" # args.dir + dir_path = "results/invivo_results/" # args.dir + # dir_path = "results/results_add__trial_more_iter_NM/invivo_results" # args.dir # if not os.path.isdir(dir_path): # print("Argument dir is not a directory.") # parser.print_usage() # exit(0) - sensitivity_analysis(dir_path, max_models=3) - - # fits_info = get_fit_info(dir_path) - # - # errors = calculate_percent_errors(fits_info) - # create_boxplots(errors) - # labels, corr_values, corrected_p_values = behaviour_correlations(fits_info, model_values=False) - # create_correlation_plot(labels, corr_values, corrected_p_values) - # - # labels, corr_values, corrected_p_values = parameter_correlations(fits_info) - # create_correlation_plot(labels, corr_values, corrected_p_values) - # - # create_parameter_distributions(get_parameter_values(fits_info)) - # cell_b, model_b = get_behaviour_values(fits_info) - # create_behaviour_distributions(cell_b, model_b) + # sensitivity_analysis(dir_path, max_models=3) + + fits_info = get_fit_info(dir_path) + + errors = calculate_percent_errors(fits_info) + create_boxplots(errors) + labels, corr_values, corrected_p_values = behaviour_correlations(fits_info, model_values=False) + create_correlation_plot(labels, corr_values, corrected_p_values) + + labels, corr_values, corrected_p_values = parameter_correlations(fits_info) + create_correlation_plot(labels, corr_values, corrected_p_values) + + create_parameter_distributions(get_parameter_values(fits_info)) + cell_b, model_b = get_behaviour_values(fits_info) + create_behaviour_distributions(cell_b, model_b) pass diff --git a/helperFunctions.py b/helperFunctions.py index ba6c884..c5a3fe2 100644 --- a/helperFunctions.py +++ b/helperFunctions.py @@ -5,6 +5,25 @@ from scipy.optimize import curve_fit import functions as fu from numba import jit import matplotlib.pyplot as plt +import time + + +def plot_errors(list_errors): + names = ["error_bf", "error_vs", "error_sc", "error_f_inf", + "error_f_inf_slope", "error_f_zero", "error_f_zero_s", "f_zero_curve"] + data = np.array(list_errors) + + fig, axes = plt.subplots(2, 4, figsize=(10, 8)) + + for i in range(8): + col = i % 4 + row = int(i/4.0) + + axes[row, col].hist(data[:, i]) + axes[row, col].set_title(names[i]) + + plt.savefig("figures/error_distributions/error_distribution_{}.png".format(time.strftime("%H:%M:%S"))) + plt.close() def fit_clipped_line(x, y):