Alexander/P-unit_model
2
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

add option to fit without dendTau or ref_period

Browse Source
This commit is contained in:
a.ott 2020-09-11 17:21:47 +02:00
parent 4c0bde5833
commit 5af4117d4c
3 changed files with 256 additions and 28 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

199
Fitter.py
View File

@ -136,9 +136,16 @@ class Fitter:
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")
# tau mins:
tau_min = 0.001
for i in (0, 3, 5):
if X[i] < tau_min:
print("tried too small tau value")
return 1000 + abs(X[i] - tau_min) * 10000
for i in (1, 2, 4, 6):
if X[i] < 0:
print("tried negative parameter value")
return 1000 + abs(X[i]) * 10000
if X[6] > 1.05/self.eod_freq: # refractory period shouldn't be larger than one eod period
@ -170,6 +177,184 @@ class Fitter:
self.best_parameters_found = X
return sum(error_list)
def fit_routine_no_dend_tau(self, start_parameters, error_weights=None):
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["refractory_period"]])
initial_simplex = create_init_simples(x0, search_scale=3)
# 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]
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})
return fmin, self.base_model.get_parameters()
def cost_function_no_dend_tau(self, X, error_weights=None):
# tau mins:
tau_min = 0.001
for i in (0, 3):
if X[i] < tau_min:
print("tried too small tau value")
return 1000 + abs(X[i] - tau_min) * 10000
for i in (1, 2, 4, 5):
if X[i] < 0:
print("tried negative parameter value")
return 1000 + abs(X[i]) * 10000
if X[6] > 1.05/self.eod_freq: # refractory period shouldn't be larger than one eod period
print("tried too large ref period")
return 1000 + abs(X[6]) * 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("refractory_period", 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)
# 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_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
return sum(error_list)
def fit_routine_no_ref_period(self, start_parameters, error_weights=None):
self.counter = 0
# fit only v_offset, mem_tau, input_scaling, dend_tau
self.base_model.set_variable("refractory_period", 0)
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"]])
initial_simplex = create_init_simples(x0, search_scale=3)
# 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]
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})
return fmin, self.base_model.get_parameters()
def cost_function_no_ref_period(self, X, error_weights=None):
# tau mins:
tau_min = 0.001
for i in (0, 3, 5):
if X[i] < tau_min:
print("tried too small tau value")
return 1000 + abs(X[i] - tau_min) * 10000
for i in (1, 2, 4):
if X[i] < 0:
print("tried negative parameter value")
return 1000 + abs(X[i]) * 10000
if X[6] > 1.05/self.eod_freq: # refractory period shouldn't be larger than one eod period
print("tried too large ref period")
return 1000 + abs(X[6]) * 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])
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_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
return sum(error_list)
def fit_routine_no_dend_tau_and_no_ref_period(self, start_parameters, error_weights=None):
self.counter = 0
# fit only v_offset, mem_tau, input_scaling, dend_tau
self.base_model.set_variable("refractory_period", 0)
self.base_model.parameters["dend_tau"] = 0
x0 = np.array([start_parameters["mem_tau"], start_parameters["noise_strength"],
start_parameters["input_scaling"], start_parameters["tau_a"], start_parameters["delta_a"]])
initial_simplex = create_init_simples(x0, search_scale=3)
# 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]
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})
return fmin, self.base_model.get_parameters()
def cost_function_no_dend_tau_and_no_ref_period(self, X, error_weights=None):
# tau mins:
tau_min = 0.001
for i in (0, 3):
if X[i] < tau_min:
print("tried too small tau value")
return 1000 + abs(X[i] - tau_min) * 10000
for i in (1, 2, 4):
if X[i] < 0:
print("tried negative parameter value")
return 1000 + abs(X[i]) * 10000
if X[6] > 1.05/self.eod_freq: # refractory period shouldn't be larger than one eod period
print("tried too large ref period")
return 1000 + abs(X[6]) * 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])
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_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
return sum(error_list)
def calculate_errors(self, error_weights=None, model=None):
if model is None:
model = self.base_model
@ -202,7 +387,7 @@ class Fitter:
error_vs = abs((vector_strength - self.vector_strength) / 0.01)
error_cv = abs((coefficient_of_variation - self.coefficient_of_variation) / 0.05)
error_bursty = (abs(burstiness - self.burstiness) / 0.2)
error_hist = np.mean((isi_bins - self.isi_bins) ** 2) / 600
error_hist = np.sqrt(np.mean((isi_bins - self.isi_bins) ** 2)) / 10
# print("error hist: {:.2f}".format(error_hist))
# print("Burstiness: cell {:.2f}, model: {:.2f}, error: {:.2f}".format(self.burstiness, burstiness, error_bursty))
@ -211,7 +396,7 @@ class Fitter:
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_slope = abs((f_infinities_slope - self.f_inf_slope) / abs(self.f_inf_slope+1)) * 25
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)
@ -219,7 +404,7 @@ class Fitter:
/ abs(self.f_zero_slope_at_straight+1)
error_f_zero = calculate_list_error(f_zeros, self.f_zero_values) / 10
error_f0_curve = self.calculate_f0_curve_error(model, fi_curve_model) / 1000
error_f0_curve = self.calculate_f0_curve_error(model, fi_curve_model) / 20
error_list = [error_vs, error_sc, error_cv, error_hist, error_bursty,
error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope_at_straight, error_f0_curve]
@ -290,7 +475,7 @@ class Fitter:
elif len(model_curve) < len(cell_curve):
cell_curve = cell_curve[:len(model_curve)]
error_f0_curve = np.mean((model_curve - cell_curve) ** 2)
error_f0_curve = np.sqrt(np.mean((model_curve - cell_curve) ** 2))
return error_f0_curve

54
models/LIFACnoise.py
View File

@ -128,8 +128,10 @@ class LifacNoiseModel(AbstractModel):
parameters = np.array(
[v_zero, a_zero, step_size, threshold, v_base, delta_a, tau_a, v_offset, mem_tau, noise_strength,
time_start, input_scaling, dend_tau, ref_period])
voltage_trace, adaption, spiketimes, input_voltage = simulate_fast(rectified_stimulus, total_time_s, parameters)
if dend_tau >= step_size:
voltage_trace, adaption, spiketimes, input_voltage = simulate_fast(rectified_stimulus, total_time_s, parameters)
else:
voltage_trace, adaption, spiketimes, input_voltage = simulate_fast_no_dend_tau(rectified_stimulus, total_time_s, parameters)
self.stimulus = stimulus
self.input_voltage = input_voltage
@ -304,3 +306,51 @@ def simulate_fast(rectified_stimulus_array, total_time_s, parameters: np.ndarray
adaption[i] += delta_a / tau_a
return output_voltage, adaption, spiketimes, input_voltage
@jit(nopython=True)
def simulate_fast_no_dend_tau(rectified_stimulus_array, total_time_s, parameters: np.ndarray):
v_zero = parameters[0]
a_zero = parameters[1]
step_size = parameters[2]
threshold = parameters[3]
v_base = parameters[4]
delta_a = parameters[5]
tau_a = parameters[6]
v_offset = parameters[7]
mem_tau = parameters[8]
noise_strength = parameters[9]
time_start = parameters[10]
input_scaling = parameters[11]
dend_tau = parameters[12]
ref_period = parameters[13]
time = np.arange(time_start, total_time_s, step_size)
length = len(time)
output_voltage = np.zeros(length)
adaption = np.zeros(length)
input_voltage = rectified_stimulus_array
spiketimes = []
output_voltage[0] = v_zero
adaption[0] = a_zero
for i in range(1, len(time), 1):
noise_value = np.random.normal()
noise = noise_strength * noise_value / np.sqrt(step_size)
output_voltage[i] = output_voltage[i - 1] + ((v_base - output_voltage[i - 1] + v_offset + (
input_voltage[i] * input_scaling) - adaption[i - 1] + noise) / mem_tau) * step_size
adaption[i] = adaption[i - 1] + ((-adaption[i - 1]) / tau_a) * step_size
if len(spiketimes) > 0 and time[i] - spiketimes[-1] < ref_period + step_size/2:
output_voltage[i] = v_base
if output_voltage[i] > threshold:
output_voltage[i] = v_base
spiketimes.append((i * step_size) + time_start)
adaption[i] += delta_a / tau_a
return output_voltage, adaption, spiketimes, input_voltage

31
run_Fitter.py
View File

@ -4,7 +4,7 @@ from CellData import CellData
from Baseline import get_baseline_class
from FiCurve import get_fi_curve_class
from Fitter import Fitter
from ModelFit import ModelFit
from ModelFit import get_best_fit, ModelFit
import time
import os
@ -14,8 +14,8 @@ from helperFunctions import plot_errors
import multiprocessing as mp
SAVE_DIRECTORY = "./results/final_2/"
SAVE_DIRECTORY_BEST = "./results/final_2_best/"
SAVE_DIRECTORY = "./results/final_3/"
SAVE_DIRECTORY_BEST = "./results/final_3_best/"
# [vs, sc, cv, isi_hist, bursty, f_inf, f_inf_slope, f_zero, f_zero_slope, f0_curve]
ERROR_WEIGHTS = (1, 1, 1, 1, 1, 1, 1, 1, 0, 1)
@ -45,7 +45,7 @@ def test_single_cell(path):
for i, p in enumerate(start_parameters):
fitter = Fitter()
fitter.set_data_reference_values(cell_data)
fmin, res_par = fitter.fit_routine(p, ERROR_WEIGHTS)
fmin, res_par = fitter.fit_routine_no_ref_period(p, ERROR_WEIGHTS)
cell_path = os.path.split(cell_data.get_data_path())[-1]
@ -63,17 +63,19 @@ def fit_cell_base(parameters):
fmin, res_par = fitter.fit_routine(parameters[2], ERROR_WEIGHTS)
cell_data = parameters[0]
cell_path = os.path.split(cell_data.get_data_path())[-1]
cell_name = os.path.split(cell_data.get_data_path())[-1]
error = fitter.calculate_errors(model=LifacNoiseModel(res_par))
save_path = SAVE_DIRECTORY + "/" + cell_path + "/start_parameter_{:}_err_{:.2f}/".format(parameters[1], sum(error))
save_path = SAVE_DIRECTORY + "/" + cell_name + "/start_parameter_{:}_err_{:.2f}/".format(parameters[1], sum(error))
save_fitting_run_info(parameters[0], res_par, parameters[2], plot=True, save_path=save_path)
plot_errors(fitter.errors, save_path)
fit = ModelFit(save_path)
fit.generate_master_plot(save_path)
time2 = time.time()
del fitter
print("Time taken for " + cell_path +
print("Time taken for " + cell_name +
"\n and start parameters ({:}): {:.2f}s thread time".format(parameters[1]+1, time2 - time1) +
"\n error: {:.2f}".format(sum(error)))
@ -97,19 +99,10 @@ def fit_cell_parallel(cell_data, start_parameters):
del pool
del cell_data
save_master_plot(save_path_cell)
def save_master_plot(save_path_cell):
best_fit = None
min_err = np.inf
for fit in os.listdir(save_path_cell):
cur_fit = ModelFit(os.path.join(save_path_cell, fit))
if cur_fit.comparable_error() < min_err:
min_err = cur_fit.comparable_error()
best_fit = cur_fit
best_fit = get_best_fit(save_path_cell)
best_fit.generate_master_plot(SAVE_DIRECTORY_BEST)
best_fit.generate_master_plot(SAVE_DIRECTORY)
def iget_start_parameters():