446 lines
19 KiB
Python
446 lines
19 KiB
Python
|
|
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 warnings import warn
|
|
from scipy.optimize import minimize
|
|
import time
|
|
import os
|
|
|
|
|
|
SAVE_PATH_PREFIX = ""
|
|
|
|
|
|
def main():
|
|
run_with_real_data()
|
|
|
|
|
|
def iget_start_parameters():
|
|
# mem_tau, input_scaling, noise_strength, dend_tau,
|
|
# expand by tau_a, delta_a ?
|
|
|
|
mem_tau_list = [0.01]
|
|
input_scaling_list = [40, 60, 80]
|
|
noise_strength_list = [0.03] # [0.02, 0.06]
|
|
dend_tau_list = [0.001, 0.002]
|
|
|
|
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}
|
|
|
|
|
|
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
|
|
print("START PARAMETERS:", start_par_count)
|
|
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),
|
|
"\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)))
|
|
# 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, 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]
|
|
|
|
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)
|
|
c_cv = cell_data.get_coefficient_of_variation()
|
|
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("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)
|
|
|
|
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())
|
|
|
|
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.coefficient_of_variation = 0
|
|
|
|
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, start_parameters=None):
|
|
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)
|
|
self.coefficient_of_variation = data.get_coefficient_of_variation()
|
|
|
|
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()
|
|
|
|
# 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())
|
|
|
|
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))
|
|
|
|
return self.fit_routine_5(data, start_parameters)
|
|
# return self.fit_model(fit_adaption=False)
|
|
|
|
def fit_routine_5(self, cell_data=None, start_parameters=None):
|
|
global SAVE_PATH_PREFIX
|
|
SAVE_PATH_PREFIX = "fit_routine_5_"
|
|
# 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["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, 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})
|
|
|
|
if cell_data is not None:
|
|
print("##### After step 1: (Everything)")
|
|
# print_comparision_cell_model(cell_data, res_parameters_step2)
|
|
|
|
return fmin, self.base_model.get_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])
|
|
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
|
|
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])
|
|
self.base_model.set_variable("mem_tau", X[2])
|
|
|
|
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_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 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
|
|
model.set_variable("mem_tau", X[0])
|
|
model.set_variable("noise_strength", X[1])
|
|
model.set_variable("input_scaling", X[2])
|
|
model.set_variable("dend_tau", X[3])
|
|
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, 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)
|
|
f_baselines, f_zeros, f_infinities = self.base_model.calculate_fi_curve(self.fi_contrasts, self.eod_freq)
|
|
try:
|
|
f_infinities_fit = hF.fit_clipped_line(self.fi_contrasts, f_infinities)
|
|
except Exception as e:
|
|
print("EXCEPTION IN FIT LINE!")
|
|
print(e)
|
|
f_infinities_fit = [0, 0]
|
|
|
|
f_infinities_slope = f_infinities_fit[0]
|
|
try:
|
|
f_zeros_fit = hF.fit_boltzmann(self.fi_contrasts, f_zeros)
|
|
except Exception as e:
|
|
print("EXCEPTION IN FIT BOLTZMANN!")
|
|
print(e)
|
|
f_zeros_fit = [0, 0, 0, 0]
|
|
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_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
|
|
|
|
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_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]
|
|
|
|
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:
|
|
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),
|
|
"Coefficient of variation - expected: {:.2f}, current: {:.2f}, error: {:.3f}\n".format(
|
|
self.coefficient_of_variation, coefficient_of_variation, error_cv),
|
|
"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 = 10
|
|
error += abs((fit[i] - reference[i])+constant) / (abs(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()
|