from models.LIFACnoise import LifacNoiseModel
from stimuli.SinusoidalStepStimulus import SinusoidalStepStimulus
from CellData import CellData
from Baseline import get_baseline_class
from FiCurve import get_fi_curve_class
from AdaptionCurrent import Adaption
import numpy as np
from warnings import warn
from scipy.optimize import minimize
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 = 2
# values to be replicated:
self.baseline_freq = 0
self.vector_strength = -1
self.serial_correlation = []
self.coefficient_of_variation = 0
self.burstiness = -1
self.f_inf_values = []
self.f_inf_slope = 0
self.f_zero_values = []
# self.f_zero_slopes = []
self.f_zero_slope_at_straight = 0
self.f_zero_straight_contrast = 0
self.f_zero_fit = []
self.tau_a = 0
# counts how often the cost_function was called
self.counter = 0
def set_data_reference_values(self, cell_data: CellData):
self.eod_freq = cell_data.get_eod_frequency()
data_baseline = get_baseline_class(cell_data)
self.baseline_freq = data_baseline.get_baseline_frequency()
self.vector_strength = data_baseline.get_vector_strength()
self.serial_correlation = data_baseline.get_serial_correlation(self.sc_max_lag)
self.coefficient_of_variation = data_baseline.get_coefficient_of_variation()
self.burstiness = data_baseline.get_burstiness()
fi_curve = get_fi_curve_class(cell_data, cell_data.get_fi_contrasts())
self.fi_contrasts = fi_curve.stimulus_values
self.f_inf_values = fi_curve.f_inf_frequencies
self.f_inf_slope = fi_curve.get_f_inf_slope()
self.f_zero_values = fi_curve.f_zero_frequencies
self.f_zero_fit = fi_curve.f_zero_fit
# self.f_zero_slopes = [fi_curve.get_f_zero_fit_slope_at_stimulus_value(c) for c in self.fi_contrasts]
self.f_zero_slope_at_straight = fi_curve.get_f_zero_fit_slope_at_straight()
self.f_zero_straight_contrast = self.f_zero_fit[3]
# 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())
adaption = Adaption(fi_curve)
self.tau_a = adaption.get_tau_real()
def fit_model_to_data(self, data: CellData, start_parameters, fit_routine_func: callable):
self.set_data_reference_values(data)
return fit_routine_func(start_parameters)
def fit_routine_1(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)
# 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_all,
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()
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})
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])
self.base_model.set_variable("refractory_period", X[6])
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_without_ref_period(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("mem_tau", X[0])
self.base_model.set_variable("input_scaling", X[1])
self.base_model.set_variable("delta_a", X[2])
self.base_model.set_variable("dend_tau", X[3])
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_tau(self, X, tau_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("delta_a", X[3])
model.set_variable("dend_tau", X[4])
model.set_variable("tau_a", tau_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, model=None):
if model is None:
model = self.base_model
model_baseline = get_baseline_class(model, self.eod_freq)
baseline_freq = model_baseline.get_baseline_frequency()
vector_strength = model_baseline.get_vector_strength()
serial_correlation = model_baseline.get_serial_correlation(self.sc_max_lag)
coefficient_of_variation = model_baseline.get_coefficient_of_variation()
burstiness = model_baseline.get_burstiness()
fi_curve_model = get_fi_curve_class(model, self.fi_contrasts, self.eod_freq)
f_zeros = fi_curve_model.get_f_zero_frequencies()
f_infinities = fi_curve_model.get_f_inf_frequencies()
f_infinities_slope = fi_curve_model.get_f_inf_slope()
# f_zero_slopes = [fi_curve_model.get_f_zero_fit_slope_at_stimulus_value(x) for x in self.fi_contrasts]
f_zero_slope_at_straight = fi_curve_model.get_f_zero_fit_slope_at_stimulus_value(self.f_zero_straight_contrast)
# calculate errors with reference values
error_bf = abs((baseline_freq - self.baseline_freq) / self.baseline_freq)
error_vs = abs((vector_strength - self.vector_strength) / 0.1)
error_cv = abs((coefficient_of_variation - self.coefficient_of_variation) / 0.1)
error_bursty = (abs(burstiness - self.burstiness) / 0.02)
error_sc = 0
for i in range(self.sc_max_lag):
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) / (self.f_inf_slope/20))
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)
error_f_zero_slope_at_straight = abs(self.f_zero_slope_at_straight - f_zero_slope_at_straight) \
/ (self.f_zero_slope_at_straight / 10)
error_f_zero = calculate_list_error(f_zeros, self.f_zero_values)
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]
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 error_weights is not None:
warn("Error: weights had different length than errors and were ignored!")
return error_list
def calculate_list_error(fit, reference):
error = 0
for i in range(len(reference)):
error += abs_freq_error(fit[i] - reference[i])
norm_error = error / len(reference)
return norm_error
def abs_freq_error(diff, factor=10):
return abs(diff) / factor
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__':
print("use run_fitter.py to run the Fitter.")