diff --git a/Fitter.py b/Fitter.py index 8b0c73b..82b6b6d 100644 --- a/Fitter.py +++ b/Fitter.py @@ -12,60 +12,104 @@ import time import os -def main(): - run_with_real_data() - +SAVE_PATH_PREFIX = "" -def run_with_real_data(): - for cell_data in icelldata_of_dir("./data/"): - 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) +def main(): + run_with_real_data() - start_time = time.time() - fitter = Fitter() - fmin, parameters = fitter.fit_model_to_data(cell_data) - print(fmin) - print(parameters) - end_time = time.time() +def iget_start_parameters(mem_tau_list=None, input_scaling_list=None, noise_strength_list=None, dend_tau_list=None): + # mem_tau, input_scaling, noise_strength, dend_tau, + # expand by tau_a, delta_a ? + if mem_tau_list is None: + mem_tau_list = [0.01] + if input_scaling_list is None: + input_scaling_list = [40, 60, 80] + if noise_strength_list is None: + noise_strength_list = [0.02, 0.06] + if dend_tau_list is None: + dend_tau_list = [0.001, 0.002] - if not os.path.exists(results_path): - os.makedirs(results_path) + for mem_tau in mem_tau_list: + for input_scaling in input_scaling_list: + for noise_strength in noise_strength_list: + for dend_tau in dend_tau_list: + yield {"mem_tau": mem_tau, "input_scaling": input_scaling, + "noise_strength": noise_strength, "dend_tau": dend_tau} - with open(results_path + "fit_parameters.txt", "w") as file: - file.writelines([str(parameters)]) - results_path += "fit_routine_3_" - 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) +def run_with_real_data(): + for cell_data in icelldata_of_dir("./data/"): + start_par_count = 0 + for start_parameters in iget_start_parameters(): + start_par_count += 1 + if start_par_count <= 4: + 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, start_parameters) + + 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_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)]) + + 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))) + #print(results_path) + print_comparision_cell_model(cell_data, parameters, plot=True, savepath=results_path) break + from Sounds import play_finished_sound + play_finished_sound() pass def print_comparision_cell_model(cell_data, parameters, plot=False, savepath=None): 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_f_values, m_f_slope = res_model.calculate_fi_markers(cell_data.get_fi_contrasts(), 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] + + f_zeros_fit = hF.fit_boltzmann(fi_curve.stimulus_value, f_zeros) + m_f_zero_slope = fu.full_boltzmann_straight_slope(f_zeros_fit[0], f_zeros_fit[1], f_zeros_fit[2], f_zeros_fit[3]) 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 + + c_f_zero_slope = fi_curve.get_fi_curve_slope_of_straight() + c_f_zero_values = fi_curve.f_zeros 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) + 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) + + print("f_zero_slope: cell - {:.2f} vs model {:.2f}".format(c_f_zero_slope, m_f_zero_slope)) + print("f zero values:\n cell -", c_f_zero_values, "\n model -", f_zeros) if plot: f_b, f_zero, f_inf = res_model.calculate_fi_curve(cell_data.get_fi_contrasts(), cell_data.get_eod_frequency()) @@ -108,7 +152,7 @@ class Fitter: # counts how often the cost_function was called self.counter = 0 - def fit_model_to_data(self, data: CellData): + def fit_model_to_data(self, data: CellData, start_parameters=None): self.eod_freq = data.get_eod_frequency() self.baseline_freq = data.get_base_frequency() @@ -123,18 +167,20 @@ class Fitter: 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.f_zero_slope = fi_curve.get_fi_curve_slope_at(fi_curve.get_f_zero_and_f_inf_intersection()) # around 1/3 of the value at straight + # self.f_zero_slope = fi_curve.get_fi_curve_slope_at(fi_curve.get_f_zero_and_f_inf_intersection()) # around 1/3 of the value at straight self.delta_a = (self.f_zero_slope / self.f_inf_slope) / 1000 # seems to work if divided by 1000... adaption = Adaption(data, fi_curve) self.tau_a = adaption.get_tau_real() - print("delta_a: {:.3f}".format(self.delta_a), "tau_a: {:.3f}".format(self.tau_a)) + # print("delta_a: {:.3f}".format(self.delta_a), "tau_a: {:.3f}".format(self.tau_a)) - return self.fit_routine_3(data) + return self.fit_routine_4(data, start_parameters) # return self.fit_model(fit_adaption=False) def fit_routine_1(self, cell_data=None): + global SAVE_PATH_PREFIX + SAVE_PATH_PREFIX = "fit_routine_1_" # 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 @@ -149,7 +195,6 @@ class Fitter: 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"], @@ -167,6 +212,8 @@ class Fitter: return fmin_step2, res_parameters_step2 def fit_routine_2(self, cell_data=None): + global SAVE_PATH_PREFIX + SAVE_PATH_PREFIX = "fit_routine_2_" # 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 @@ -181,6 +228,8 @@ class Fitter: return fmin, res_parameters def fit_routine_3(self, cell_data=None): + global SAVE_PATH_PREFIX + SAVE_PATH_PREFIX = "fit_routine_3_" # 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, dend_tau @@ -194,6 +243,64 @@ class Fitter: return fmin, res_parameters + def fit_routine_4(self, cell_data=None, start_parameters=None): + global SAVE_PATH_PREFIX + SAVE_PATH_PREFIX = "fit_routine_4_" + # 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, input_scaling, dend_tau + 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"], start_parameters["dend_tau"]]) + initial_simplex = create_init_simples(x0, search_scale=2) + error_weights = (0, 1, 5, 1, 2, 0, 0) + fmin = minimize(fun=self.cost_function_with_fixed_adaption_with_dend_tau_no_noise, + args=(self.tau_a, self.delta_a, error_weights), x0=x0, method="Nelder-Mead", + options={"initial_simplex": initial_simplex, "xatol": 0.001}) + res_parameters = fmin.x + + print_comparision_cell_model(cell_data, self.base_model.get_parameters()) + + self.counter = 0 + x0 = np.array([res_parameters[0], res_parameters[1], self.tau_a, + self.delta_a, res_parameters[2]]) + initial_simplex = create_init_simples(x0, search_scale=2) + error_weights = (0, 0, 0, 0, 0, 4, 2) + fmin = minimize(fun=self.cost_function_all_without_noise, + args=(error_weights,), x0=x0, method="Nelder-Mead", + options={"initial_simplex": initial_simplex, "xatol": 0.001}) + res_parameters = fmin.x + print(fmin) + print_comparision_cell_model(cell_data, self.base_model.get_parameters()) + + self.counter = 0 + x0 = np.array([res_parameters[0], + res_parameters[1], self.tau_a, + self.delta_a, res_parameters[2]]) + initial_simplex = create_init_simples(x0, search_scale=2) + error_weights = (0, 0, 1, 0, 0, 5, 2) + fmin = minimize(fun=self.cost_function_all_without_noise, + args=(error_weights,), x0=x0, method="Nelder-Mead", + options={"initial_simplex": initial_simplex, "xatol": 0.001}) + res_parameters = self.base_model.get_parameters() + + print_comparision_cell_model(cell_data, self.base_model.get_parameters()) + + # noise_strength = 0.03 + # self.counter = 0 + # x0 = np.array([res_parameters["mem_tau"], noise_strength, + # res_parameters["input_scaling"], res_parameters["tau_a"], + # res_parameters["delta_a"], res_parameters["dend_tau"]]) + # initial_simplex = create_init_simples(x0, search_scale=2) + # error_weights = (0, 2, 2, 1, 1, 3, 2) + # fmin = minimize(fun=self.cost_function_all, + # args=(error_weights,), x0=x0, method="Nelder-Mead", + # options={"initial_simplex": initial_simplex, "xatol": 0.001}) + # res_parameters = self.base_model.get_parameters() + + return fmin, self.base_model.get_parameters() + def fit_model(self, x0=None, initial_simplex=None, fit_adaption=False): self.counter = 0 @@ -220,6 +327,27 @@ class Fitter: 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) + 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_all_without_noise(self, X, error_weights=None): + self.base_model.set_variable("mem_tau", X[0]) + self.base_model.set_variable("input_scaling", X[1]) + self.base_model.set_variable("tau_a", X[2]) + self.base_model.set_variable("delta_a", X[3]) + self.base_model.set_variable("dend_tau", X[4]) + self.base_model.set_variable("noise_strength", 0) base_stimulus = SinusoidalStepStimulus(self.eod_freq, 0) # find right v-offset @@ -276,6 +404,27 @@ class Fitter: return sum(error_list) + def cost_function_with_fixed_adaption_with_dend_tau_no_noise(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("input_scaling", X[1]) + model.set_variable("dend_tau", X[2]) + model.set_variable("tau_a", tau_a) + model.set_variable("delta_a", delta_a) + model.set_variable("noise_strength", 0) + + 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 cost_function_with_fixed_adaption_with_dend_tau(self, X, tau_a, delta_a, error_weights=None): # set model parameters: model = self.base_model @@ -329,7 +478,7 @@ class Fitter: error = sum(error_list) self.counter += 1 - if self.counter % 200 == 0: + if self.counter % 200 == 0 and False: # TODO currently shut off! 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( @@ -354,8 +503,8 @@ 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 = 0 - error += abs((fit[i] - reference[i]) / (reference[i] + constant)) + constant = 50 + error += abs((fit[i] - reference[i])+constant) / (abs(reference[i]) + constant) norm_error = error / len(reference) diff --git a/Sounds.py b/Sounds.py new file mode 100644 index 0000000..afa1faa --- /dev/null +++ b/Sounds.py @@ -0,0 +1,40 @@ +from pyaudio import PyAudio +import numpy as np + +BITRATE = 20000 +LENGTH = 0.5 + + +def play_finished_sound(): + global BITRATE + global LENGTH + + frequency = 261 + num_of_frames = int(BITRATE*LENGTH) + + frames = np.arange(0, num_of_frames, 1) + + wave_data_numeric = np.sin(frames / ((BITRATE / frequency) / np.pi)) * 127 + 128 + wave_data_numeric = wave_data_numeric.astype(int) + wave_data_chr = "".join([chr(x) for x in wave_data_numeric]) + + rest_frames = num_of_frames % BITRATE + rest = [chr(128)]*rest_frames + + wave_data_chr.join(rest) + + p = PyAudio() + stream = p.open( + format=p.get_format_from_width(1), + channels=1, + rate=BITRATE, + output=True, + ) + stream.write(wave_data_chr) + stream.stop_stream() + stream.close() + p.terminate() + + +if __name__ == '__main__': + play_finished_sound() \ No newline at end of file diff --git a/helperFunctions.py b/helperFunctions.py index 1d7814b..15f6a12 100644 --- a/helperFunctions.py +++ b/helperFunctions.py @@ -427,6 +427,10 @@ def detect_f_zero_in_frequency_trace(time, frequency, stimulus_start, sampling_i stimulus_start = stimulus_start - time[0] # time start is generally != 0 and != delay freq_before = frequency[int(buffer/sampling_interval):int((stimulus_start - buffer) / sampling_interval)] + + if len(freq_before) < 3: + print("mäh") + min_before = min(freq_before) max_before = max(freq_before) mean_before = np.mean(freq_before) diff --git a/models/LIFACnoise.py b/models/LIFACnoise.py index 3241aa5..f77716e 100644 --- a/models/LIFACnoise.py +++ b/models/LIFACnoise.py @@ -227,6 +227,12 @@ class LifacNoiseModel(AbstractModel): # if c == contrasts[0] or c == contrasts[-1]: # plt.plot(frequency) # plt.show() + + if len(time) == 0 or time[0] >= stim_start or len(spiketimes) < 5: + f_infinities.append(0) + f_zeros.append(0) + f_baselines.append(0) + continue f_inf = hF.detect_f_infinity_in_freq_trace(time, frequency, stim_start, stim_duration, sampling_interval) f_infinities.append(f_inf)