428 lines
18 KiB
Python
428 lines
18 KiB
Python
|
|
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
|
|
import time
|
|
from helperFunctions import plot_errors
|
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
|
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.best_parameters_found = []
|
|
self.smallest_error = np.inf
|
|
|
|
#
|
|
self.fi_contrasts = []
|
|
self.eod_freq = 0
|
|
self.data_sampling_interval = -1
|
|
|
|
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.f_zero_curve_contrast = 0
|
|
self.f_zero_curve_contrast_idx = -1
|
|
self.f_zero_curve_freq = np.array([])
|
|
self.f_zero_curve_time = np.array([])
|
|
|
|
self.errors = []
|
|
|
|
# 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()
|
|
self.data_sampling_interval = cell_data.get_sampling_interval()
|
|
|
|
data_baseline = get_baseline_class(cell_data)
|
|
data_baseline.load_values(cell_data.get_data_path())
|
|
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(), save_dir=cell_data.get_data_path())
|
|
self.f_inf_slope = fi_curve.get_f_inf_slope()
|
|
contrasts = np.array(cell_data.get_fi_contrasts())
|
|
if self.f_inf_slope < 0:
|
|
contrasts = contrasts * -1
|
|
# print("old contrasts:", cell_data.get_fi_contrasts())
|
|
# print("new contrasts:", contrasts)
|
|
contrasts = sorted(contrasts)
|
|
fi_curve = get_fi_curve_class(cell_data, 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]
|
|
|
|
max_contrast = max(contrasts)
|
|
test_contrast = 0.5 * max_contrast
|
|
diff_contrasts = np.abs(contrasts - test_contrast)
|
|
|
|
self.f_zero_curve_contrast_idx = np.argmin(diff_contrasts)
|
|
self.f_zero_curve_contrast = contrasts[self.f_zero_curve_contrast_idx]
|
|
times, freqs = fi_curve.get_mean_time_and_freq_traces()
|
|
self.f_zero_curve_freq = freqs[self.f_zero_curve_contrast_idx]
|
|
self.f_zero_curve_time = times[self.f_zero_curve_contrast_idx]
|
|
|
|
# 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"], start_parameters["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, 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": 400})
|
|
|
|
plot_errors(self.errors)
|
|
|
|
return fmin, self.base_model.get_parameters()
|
|
|
|
def cost_function_all(self, X, error_weights=None):
|
|
for i in range(len(X)):
|
|
if X[i] < 0:
|
|
print("tried impossible value")
|
|
return 1000 + abs(X[i]) * 10000
|
|
|
|
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)
|
|
time1 = time.time()
|
|
v_offset = test_model.find_v_offset(self.baseline_freq, base_stimulus)
|
|
self.base_model.set_variable("v_offset", v_offset)
|
|
time2 = time.time()
|
|
# print("time taken for finding v_offset: {:.2f}s".format(time2-time1))
|
|
|
|
# [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)
|
|
# print("sum: {:.2f}, ".format(sum(error_list)))
|
|
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):
|
|
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
|
|
|
|
time1 = time.time()
|
|
model_baseline = get_baseline_class(model, self.eod_freq, trials=5)
|
|
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()
|
|
time2 = time.time()
|
|
|
|
# print("Time taken for all baseline parameters: {:.2f}".format(time2-time1))
|
|
|
|
time1 = time.time()
|
|
fi_curve_model = get_fi_curve_class(model, self.fi_contrasts, self.eod_freq, trials=15)
|
|
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)
|
|
|
|
time2 = time.time()
|
|
|
|
# print("Time taken for all fi-curve parameters: {:.2f}".format(time2 - time1))
|
|
|
|
# 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.2)
|
|
|
|
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) / abs(self.f_inf_slope+1/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) \
|
|
/ abs(self.f_zero_slope_at_straight+1 / 10)
|
|
error_f_zero = calculate_list_error(f_zeros, self.f_zero_values)
|
|
|
|
times, freqs = fi_curve_model.get_mean_time_and_freq_traces()
|
|
freq_prediction = np.array(freqs[self.f_zero_curve_contrast_idx])
|
|
time_prediction = np.array(times[self.f_zero_curve_contrast_idx])
|
|
stimulus_start = fi_curve_model.get_stimulus_start() - time_prediction[0]
|
|
length = fi_curve_model.get_stimulus_duration() / 2
|
|
|
|
if model.get_sampling_interval() == self.data_sampling_interval:
|
|
start_idx = int(stimulus_start / fi_curve_model.get_sampling_interval())
|
|
end_idx = int((stimulus_start + length) / model.get_sampling_interval())
|
|
start_idx_cell = start_idx
|
|
start_idx_model = start_idx
|
|
end_idx_cell = end_idx
|
|
end_idx_model = end_idx
|
|
step_cell = 1
|
|
step_model = 1
|
|
else:
|
|
start_idx_cell = int(stimulus_start / self.data_sampling_interval)
|
|
start_idx_model = int(stimulus_start / fi_curve_model.get_sampling_interval())
|
|
end_idx_cell = int((stimulus_start + length) / self.data_sampling_interval)
|
|
end_idx_model = int((stimulus_start + length) / model.get_sampling_interval())
|
|
if round(model.get_sampling_interval() % self.data_sampling_interval, 4) == 0:
|
|
step_cell = int(model.get_sampling_interval() / self.data_sampling_interval)
|
|
step_model = 1
|
|
else:
|
|
raise ValueError("Model sampling interval is not a multiple of data sampling interval.")
|
|
|
|
if len(time_prediction) == 0 or len(time_prediction) < end_idx_model or time_prediction[0] > fi_curve_model.get_stimulus_start():
|
|
error_f0_curve = 200
|
|
else:
|
|
data_curve = self.f_zero_curve_freq[start_idx_cell:end_idx_cell:step_cell]
|
|
model_curve = freq_prediction[start_idx_model:end_idx_model:step_model]
|
|
if len(data_curve) < len(model_curve):
|
|
model_curve = model_curve[:len(data_curve)]
|
|
elif len(model_curve) < len(data_curve):
|
|
data_curve = data_curve[:len(model_curve)]
|
|
|
|
error_f0_curve = np.mean((model_curve - data_curve)**2) / 100
|
|
|
|
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, error_f0_curve]
|
|
|
|
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!")
|
|
if sum(error_list) < 0:
|
|
print("Error negative: ", error_list)
|
|
if np.isnan(sum(error_list)):
|
|
print("--------SOME ERROR VALUE(S) IS/ARE NaN:")
|
|
print(error_list)
|
|
return [50 for e in error_list]
|
|
# raise ValueError("Some error value(s) is/are NaN!")
|
|
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])
|
|
error += normed_quadratic_freq_error(fit[i], reference[i])
|
|
norm_error = error / len(reference)
|
|
|
|
return norm_error
|
|
|
|
def calculate_f0_curve_error(data_ficurve, model_ficurve):
|
|
return 0
|
|
|
|
def normed_quadratic_freq_error(fit, ref, factor=2):
|
|
return (abs(fit-ref)/factor)**2 / ref
|
|
|
|
|
|
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.")
|