From b37db0ea36b53940113a28d5d486306c29b3c1fd Mon Sep 17 00:00:00 2001 From: "a.ott" Date: Mon, 11 May 2020 10:21:17 +0200 Subject: [PATCH] add coefficient of variation calculation + in fitter --- CellData.py | 10 ++++++++++ Fitter.py | 36 +++++++++++++++++++++++------------- fit_lifacnoise.py | 14 ++++++++++---- helperFunctions.py | 9 --------- models/LIFACnoise.py | 3 ++- 5 files changed, 45 insertions(+), 27 deletions(-) diff --git a/CellData.py b/CellData.py index 82bd890..100fc7f 100644 --- a/CellData.py +++ b/CellData.py @@ -156,6 +156,16 @@ class CellData: return mean_sc + def get_coefficient_of_variation(self): + spiketimes = self.get_base_spikes() + # CV (stddev of ISI divided by mean ISI (np.diff(spiketimes)) + cvs = [] + for st in spiketimes: + st = np.array(st) + cvs.append(hf.calculate_coefficient_of_variation(st)) + + return np.mean(cvs) + def get_eod_frequency(self): eods = self.get_base_traces(self.EOD) sampling_interval = self.get_sampling_interval() diff --git a/Fitter.py b/Fitter.py index ea1b17f..7ad753d 100644 --- a/Fitter.py +++ b/Fitter.py @@ -7,6 +7,7 @@ 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 @@ -73,8 +74,8 @@ def run_with_real_data(): with open(results_path + "fit_parameters_start_{}.txt".format(start_par_count), "w") as file: file.writelines(["start_parameters:\t" + str(start_parameters), - "final_parameters:\t" + str(parameters), - "final_fmin:\t" + str(fmin)]) + "\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))) @@ -82,8 +83,8 @@ def run_with_real_data(): print_comparision_cell_model(cell_data, parameters, plot=True, savepath=results_path) break - from Sounds import play_finished_sound - play_finished_sound() + # from Sounds import play_finished_sound + # play_finished_sound() pass @@ -91,7 +92,7 @@ def print_comparision_cell_model(cell_data, parameters, plot=False, savepath=Non res_model = LifacNoiseModel(parameters) fi_curve = FICurve(cell_data) - m_bf, m_vs, m_sc = res_model.calculate_baseline_markers(cell_data.get_eod_frequency()) + 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] @@ -102,6 +103,7 @@ def print_comparision_cell_model(cell_data, parameters, plot=False, savepath=Non 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 @@ -110,6 +112,7 @@ def print_comparision_cell_model(cell_data, parameters, plot=False, savepath=Non 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) @@ -142,6 +145,7 @@ class Fitter: 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 @@ -162,6 +166,7 @@ class Fitter: 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 @@ -315,11 +320,11 @@ class Fitter: if start_parameters is None: x0 = np.array([0.02, 70, 0.001]) else: - x0 = np.array([start_parameters["mem_tau"], start_parameters["input_scaling"], + 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, 2, 1) - fmin = minimize(fun=self.cost_function_all_without_noise, + 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 @@ -465,8 +470,8 @@ class Fitter: 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) + 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) @@ -491,6 +496,7 @@ class Fitter: 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 @@ -498,13 +504,17 @@ class Fitter: 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] + 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] - error = sum(error_list) + 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), @@ -531,7 +541,7 @@ def calculate_f_values_error(fit, reference): 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 = 50 + constant = 10 error += abs((fit[i] - reference[i])+constant) / (abs(reference[i]) + constant) norm_error = error / len(reference) diff --git a/fit_lifacnoise.py b/fit_lifacnoise.py index 12c4882..339b6c1 100644 --- a/fit_lifacnoise.py +++ b/fit_lifacnoise.py @@ -34,18 +34,22 @@ def run_with_real_data(): print('Fitting of cell took function took {:.3f} s'.format((end_time - start_time))) res_model = LifacNoiseModel(parameters) - m_bf, m_vs, m_sc = res_model.calculate_baseline_markers(celldata.get_eod_frequency()) + m_bf, m_vs, m_sc, m_cv = res_model.calculate_baseline_markers(celldata.get_eod_frequency()) m_f_values, m_f_slope = res_model.calculate_fi_markers(celldata.get_fi_contrasts(), celldata.get_eod_frequency()) c_bf = celldata.get_base_frequency() c_vs = celldata.get_vector_strength() c_sc = celldata.get_serial_correlation(1) + c_cv = celldata.get_coefficient_of_variation() + fi_curve = FICurve(celldata) 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("cv: cell - {:.2f} vs model {:.2f}".format(c_cv[0], m_cv[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) @@ -76,9 +80,9 @@ def run_test_with_fixed_model(): model = LifacNoiseModel(parameters) eod_freq = 750 contrasts = np.arange(0.5, 1.51, 0.1) - baseline_freq, vector_strength, serial_correlation = model.calculate_baseline_markers(eod_freq) + baseline_freq, vector_strength, serial_correlation, coefficient_of_variation = model.calculate_baseline_markers(eod_freq) f_infinities, f_infinities_slope = model.calculate_fi_markers(contrasts, eod_freq) - print("Baseline freq:{:.2f}\nVector strength: {:.3f}\nSerial cor:".format(baseline_freq, vector_strength), serial_correlation) + print("Baseline freq:{:.2f}\nVector strength: {:.3f}\nCV: {:.3f}\nSerial cor:".format(baseline_freq, vector_strength, coefficient_of_variation), serial_correlation) print("FI-Curve\nSlope: {:.2f}\nValues:".format(f_infinities_slope), f_infinities) fitter = Fitter() starting_conditions = [[0.05, 0.02, 50], [0.01, 0.08, 75], [0.02, 0.03, 70], @@ -126,6 +130,7 @@ class Fitter: self.baseline_freq = 0 self.vector_strength = -1 self.serial_correlation = [] + self.coefficient_of_variation = 0 self.f_infinities = [] self.f_infinities_slope = 0 @@ -157,7 +162,7 @@ class Fitter: v_offset = test_model.find_v_offset(self.baseline_freq, base_stimulus) self.model.set_variable("v_offset", v_offset) - baseline_freq, vector_strength, serial_correlation = self.model.calculate_baseline_markers(self.eod_freq, self.sc_max_lag) + baseline_freq, vector_strength, serial_correlation, coefficient_of_variation = self.model.calculate_baseline_markers(self.eod_freq, self.sc_max_lag) f_infinities, f_infinities_slope = self.model.calculate_fi_markers(self.fi_contrasts, self.eod_freq) error_bf = abs((baseline_freq - self.baseline_freq) / self.baseline_freq) @@ -200,6 +205,7 @@ class Fitter: 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 diff --git a/helperFunctions.py b/helperFunctions.py index da217b2..406c195 100644 --- a/helperFunctions.py +++ b/helperFunctions.py @@ -220,13 +220,6 @@ def mean_freq_of_spiketimes_after_time_x(spiketimes, time_x, time_in_ms=False): return mean_freq - # freq = calculate_isi_frequency_trace(spiketimes, sampling_interval, time_in_ms) - # # returned frequency starts at the - # idx = int((time_x-spiketimes[0]) / sampling_interval) - # rest_array = freq[idx:] - # mean_freq = np.mean(rest_array) - # return mean_freq - def calculate_mean_isi_freq(spiketimes, time_in_ms=False): if len(spiketimes) < 2: @@ -241,8 +234,6 @@ def calculate_mean_isi_freq(spiketimes, time_in_ms=False): return sum(freqs * weights) / sum(weights) - - # @jit(nopython=True) # only faster at around 30 000 calls def calculate_coefficient_of_variation(spiketimes: np.ndarray) -> float: # CV (stddev of ISI divided by mean ISI (np.diff(spiketimes)) diff --git a/models/LIFACnoise.py b/models/LIFACnoise.py index 224de90..7757c13 100644 --- a/models/LIFACnoise.py +++ b/models/LIFACnoise.py @@ -189,8 +189,9 @@ class LifacNoiseModel(AbstractModel): vector_strength = hF.calculate_vector_strength_from_spiketimes(time_trace, stimulus_array, spiketimes, self.get_sampling_interval()) serial_correlation = hF.calculate_serial_correlation(np.array(spiketimes), max_lag) + coeffient_of_variation = hF.calculate_coefficient_of_variation(spiketimes) - return baseline_freq, vector_strength, serial_correlation + return baseline_freq, vector_strength, serial_correlation, coeffient_of_variation def calculate_fi_markers(self, contrasts, stimulus_freq): """