271 lines
11 KiB
Python
271 lines
11 KiB
Python
from models.LIFACnoise import LifacNoiseModel
|
|
from CellData import CellData, icelldata_of_dir
|
|
from FiCurve import FICurve
|
|
from AdaptionCurrent import Adaption
|
|
from stimuli.SinusoidalStepStimulus import SinusoidalStepStimulus
|
|
import helperFunctions as hF
|
|
import numpy as np
|
|
from scipy.optimize import curve_fit, minimize
|
|
import functions as fu
|
|
import time
|
|
import matplotlib.pyplot as plt
|
|
|
|
|
|
def main():
|
|
#run_test_with_fixed_model()
|
|
#quit()
|
|
|
|
# fitter = Fitter()
|
|
# fmin, params = fitter.fit_model_to_values(700, 1400, [-0.3], 1, [0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3], [1370, 1380, 1390, 1400, 1410, 1420, 1430], 100, 0.02, 0.01)
|
|
# print("calculated parameters:")
|
|
# print(params)
|
|
run_with_real_data()
|
|
|
|
|
|
def run_with_real_data():
|
|
for celldata in icelldata_of_dir("./data/"):
|
|
start_time = time.time()
|
|
fitter = Fitter()
|
|
fmin, parameters = fitter.fit_model_to_data(celldata)
|
|
|
|
print(fmin)
|
|
print(parameters)
|
|
end_time = time.time()
|
|
|
|
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_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)
|
|
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("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)
|
|
|
|
f_b, f_zero, f_inf = res_model.calculate_fi_curve(celldata.get_fi_contrasts(), celldata.get_eod_frequency())
|
|
|
|
fi_curve.plot_fi_curve(comp_f_baselines=f_b, comp_f_zeros=f_zero, comp_f_infs=f_inf)
|
|
break
|
|
|
|
pass
|
|
|
|
|
|
def run_test_with_fixed_model():
|
|
a_tau = 0.1
|
|
a_delta = 0.08
|
|
|
|
parameters = {"mem_tau": 0.015,
|
|
"v_base": 0,
|
|
"v_zero": 0,
|
|
"threshold": 1,
|
|
"v_offset": -10,
|
|
"input_scaling": 60,
|
|
"delta_a": a_delta,
|
|
"tau_a": a_tau,
|
|
"a_zero": 10,
|
|
"noise_strength": 0.05,
|
|
"step_size": 0.00005}
|
|
|
|
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)
|
|
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("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],
|
|
[0.02, 0.03, 70]]
|
|
for x0 in starting_conditions:
|
|
init_simplex = create_init_simples(x0)
|
|
|
|
fmin, fit_parameters = fitter.fit_model_to_values(eod_freq, baseline_freq, serial_correlation, vector_strength, contrasts, f_infinities, f_infinities_slope, a_delta, a_tau, x0, init_simplex)
|
|
print(x0)
|
|
print(fmin)
|
|
print("calculated parameters:")
|
|
print(fit_parameters)
|
|
|
|
print("ref parameters:")
|
|
print(parameters)
|
|
|
|
|
|
def create_init_simples(x0, search_scale=2):
|
|
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
|
|
|
|
|
|
class Fitter:
|
|
|
|
def __init__(self, step_size=None):
|
|
if step_size is not None:
|
|
self.model = LifacNoiseModel({"step_size": step_size})
|
|
else:
|
|
self.model = LifacNoiseModel({"step_size": 0.00005})
|
|
# self.data = data
|
|
self.fi_contrasts = []
|
|
self.eod_freq = 0
|
|
|
|
self.modulation_frequency = 10
|
|
self.sc_max_lag = 1
|
|
|
|
# expected values the model has to replicate
|
|
self.baseline_freq = 0
|
|
self.vector_strength = -1
|
|
self.serial_correlation = []
|
|
|
|
self.f_infinities = []
|
|
self.f_infinities_slope = 0
|
|
|
|
# fixed values needed to fit model
|
|
self.a_tau = 0
|
|
self.a_delta = 0
|
|
|
|
self.counter = 0
|
|
|
|
def calculate_needed_values_from_data(self, data: CellData):
|
|
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)
|
|
|
|
fi_curve = FICurve(data, contrast=True)
|
|
self.fi_contrasts = fi_curve.stimulus_value
|
|
print("Fitter: fi-contrasts", self.fi_contrasts)
|
|
self.f_infinities = fi_curve.f_infinities
|
|
self.f_infinities_slope = fi_curve.get_f_infinity_slope()
|
|
|
|
f_zero_slope = fi_curve.get_fi_curve_slope_of_straight()
|
|
self.a_delta = (f_zero_slope / self.f_infinities_slope) / 1000
|
|
|
|
adaption = Adaption(data, fi_curve)
|
|
self.a_tau = adaption.get_tau_real() *2
|
|
print()
|
|
|
|
# mem_tau, (threshold?), (v_offset), noise_strength, input_scaling
|
|
def cost_function(self, X, tau_a=0.1, delta_a=0.08, error_scaling=()):
|
|
# set model parameters to the given ones:
|
|
self.model.set_variable("mem_tau", X[0])
|
|
self.model.set_variable("noise_strength", X[1])
|
|
self.model.set_variable("input_scaling", X[2])
|
|
#self.model.set_variable("tau_a", X[3])
|
|
#self.model.set_variable("delta_a", X[4])
|
|
self.model.set_variable("tau_a", tau_a)
|
|
self.model.set_variable("delta_a", delta_a)
|
|
|
|
# minimize the difference in baseline_freq first by fitting v_offset
|
|
# v_offset = self.__fit_v_offset_to_baseline_frequency__()
|
|
base_stimulus = SinusoidalStepStimulus(self.eod_freq, 0)
|
|
test_model = self.model.get_model_copy()
|
|
test_model.set_variable("noise_strength", 0)
|
|
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)
|
|
# only eod with amplitude 1 and no modulation
|
|
# _, spiketimes = self.model.simulate_fast(base_stimulus, 30)
|
|
|
|
# baseline_freq = hF.mean_freq_of_spiketimes_after_time_x(spiketimes, 5)
|
|
# # print("model:", baseline_freq, "data:", self.baseline_freq)
|
|
#
|
|
# if len(spiketimes) > 10:
|
|
# relative_spiketimes = np.array([s % (1/self.eod_freq) for s in spiketimes if s > 0])
|
|
# eod_durations = np.full((len(relative_spiketimes)), 1/self.eod_freq)
|
|
# vector_strength = hF.__vector_strength__(relative_spiketimes, eod_durations)
|
|
# serial_correlation = hF.calculate_serial_correlation(np.array(spiketimes), self.sc_max_lag)
|
|
# else:
|
|
# vector_strength = 0
|
|
# serial_correlation = [0]*self.sc_max_lag
|
|
|
|
f_infinities = []
|
|
for contrast in self.fi_contrasts:
|
|
stimulus = SinusoidalStepStimulus(self.eod_freq, contrast)
|
|
_, spiketimes = self.model.simulate_fast(stimulus, 1)
|
|
|
|
if len(spiketimes) < 2:
|
|
f_infinities.append(0)
|
|
else:
|
|
f_infinity = hF.mean_freq_of_spiketimes_after_time_x(spiketimes, 0.5)
|
|
f_infinities.append(f_infinity)
|
|
|
|
popt, pcov = curve_fit(fu.clipped_line, self.fi_contrasts, f_infinities, maxfev=10000)
|
|
|
|
f_infinities_slope = popt[0]
|
|
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_f_inf_slope = abs((f_infinities_slope - self.f_infinities_slope) / self.f_infinities_slope) * 4
|
|
#print("vs:", vector_strength, self.vector_strength)
|
|
#print("sc", serial_correlation[0], self.serial_correlation[0])
|
|
#print("f slope:", f_infinities_slope, self.f_infinities_slope)
|
|
error_f_inf = 0
|
|
for i in range(len(f_infinities)):
|
|
f_inf = f_infinities[i]
|
|
if f_inf <= 0:
|
|
f_inf = 1
|
|
error_f_inf += abs((f_infinities[i] - self.f_infinities[i]) / f_inf)
|
|
|
|
error_f_inf = error_f_inf / len(f_infinities)
|
|
|
|
# print("mem_tau:", X[0], "noise:", X[0], "input_scaling:", X[2])
|
|
errors = [error_bf, error_vs, error_sc, error_f_inf_slope, error_f_inf]
|
|
|
|
self.counter += 1
|
|
if self.counter % 100 == 0:
|
|
print("\nCost function run times: {:}\n".format(self.counter),
|
|
"Total error: {:.4f}\n".format(sum(errors)),
|
|
"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),
|
|
"f-infinity slope - expected: {:.0f}, current: {:.0f}, error: {:.3f}\n".format(self.f_infinities_slope, f_infinities_slope, error_f_inf_slope),
|
|
"f-infinity values:\nexpected:", np.around(self.f_infinities), "\ncurrent: ", np.around(f_infinities), "\nerror: {:.3f}".format(error_f_inf))
|
|
return error_bf + error_vs + error_sc + error_f_inf_slope + error_f_inf
|
|
|
|
def fit_model_to_data(self, data: CellData):
|
|
self.calculate_needed_values_from_data(data)
|
|
return self.fit_model()
|
|
|
|
def fit_model_to_values(self, eod_freq, baseline_freq, sc, vs, fi_contrasts, fi_inf_values, fi_inf_slope, a_delta, a_tau, x0=None, init_simplex=None):
|
|
self.eod_freq = eod_freq
|
|
self.baseline_freq = baseline_freq
|
|
self.serial_correlation = sc
|
|
self.vector_strength = vs
|
|
self.fi_contrasts = fi_contrasts
|
|
self.f_infinities = fi_inf_values
|
|
self.f_infinities_slope = fi_inf_slope
|
|
self.a_delta = a_delta
|
|
self.a_tau = a_tau
|
|
|
|
return self.fit_model(x0, init_simplex)
|
|
|
|
def fit_model(self, x0=None, initial_simplex=None):
|
|
self.counter = 0
|
|
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, x0=x0, args=(self.a_tau, self.a_delta), method="Nelder-Mead")
|
|
# else:
|
|
fmin = minimize(fun=self.cost_function, x0=x0, args=(self.a_tau, self.a_delta), method="Nelder-Mead", options={"initial_simplex": initial_simplex})
|
|
|
|
return fmin, self.model.get_parameters()
|
|
|
|
|
|
if __name__ == '__main__':
|
|
main()
|