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

to test remotely

Browse Source
This commit is contained in:
alexanderott 2020-04-16 15:48:42 +02:00
parent f9d4838b71
commit 0855596623
4 changed files with 237 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

225
Fitter.py
View File

@ -12,60 +12,104 @@ import time
import os
def main():
run_with_real_data()
SAVE_PATH_PREFIX = ""
def run_with_real_data():
for cell_data in icelldata_of_dir("./data/"):
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)
def main():
run_with_real_data()
start_time = time.time()
fitter = Fitter()
fmin, parameters = fitter.fit_model_to_data(cell_data)
print(fmin)
print(parameters)
end_time = time.time()
def iget_start_parameters(mem_tau_list=None, input_scaling_list=None, noise_strength_list=None, dend_tau_list=None):
# mem_tau, input_scaling, noise_strength, dend_tau,
# expand by tau_a, delta_a ?
if mem_tau_list is None:
mem_tau_list = [0.01]
if input_scaling_list is None:
input_scaling_list = [40, 60, 80]
if noise_strength_list is None:
noise_strength_list = [0.02, 0.06]
if dend_tau_list is None:
dend_tau_list = [0.001, 0.002]
if not os.path.exists(results_path):
os.makedirs(results_path)
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}
with open(results_path + "fit_parameters.txt", "w") as file:
file.writelines([str(parameters)])
results_path += "fit_routine_3_"
print('Fitting of cell took function took {:.3f} s'.format((end_time - start_time)))
print_comparision_cell_model(cell_data, parameters, plot=True, savepath=results_path)
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
if start_par_count <= 4:
continue
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),
"final_parameters:\t" + str(parameters),
"final_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 = res_model.calculate_baseline_markers(cell_data.get_eod_frequency())
m_f_values, m_f_slope = res_model.calculate_fi_markers(cell_data.get_fi_contrasts(), 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)
fi_curve = FICurve(cell_data)
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("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)
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())
@ -108,7 +152,7 @@ class Fitter:
# counts how often the cost_function was called
self.counter = 0
def fit_model_to_data(self, data: CellData):
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()
@ -123,18 +167,20 @@ class Fitter:
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()
self.f_zero_slope = fi_curve.get_fi_curve_slope_at(fi_curve.get_f_zero_and_f_inf_intersection()) # 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()) # around 1/3 of the value at straight
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))
# print("delta_a: {:.3f}".format(self.delta_a), "tau_a: {:.3f}".format(self.tau_a))
return self.fit_routine_3(data)
return self.fit_routine_4(data, start_parameters)
# return self.fit_model(fit_adaption=False)
def fit_routine_1(self, cell_data=None):
global SAVE_PATH_PREFIX
SAVE_PATH_PREFIX = "fit_routine_1_"
# 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, noise_strength, input_scaling
@ -149,7 +195,6 @@ class Fitter:
print("##### After step 1: (fixed adaption)")
print_comparision_cell_model(cell_data, res_parameters_step1)
self.counter = 0
x0 = np.array([res_parameters_step1["mem_tau"], res_parameters_step1["noise_strength"],
res_parameters_step1["input_scaling"], res_parameters_step1["tau_a"],
@ -167,6 +212,8 @@ class Fitter:
return fmin_step2, res_parameters_step2
def fit_routine_2(self, cell_data=None):
global SAVE_PATH_PREFIX
SAVE_PATH_PREFIX = "fit_routine_2_"
# 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, noise_strength, input_scaling
@ -181,6 +228,8 @@ class Fitter:
return fmin, res_parameters
def fit_routine_3(self, cell_data=None):
global SAVE_PATH_PREFIX
SAVE_PATH_PREFIX = "fit_routine_3_"
# 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, noise_strength, input_scaling, dend_tau
@ -194,6 +243,64 @@ class Fitter:
return fmin, res_parameters
def fit_routine_4(self, cell_data=None, start_parameters=None):
global SAVE_PATH_PREFIX
SAVE_PATH_PREFIX = "fit_routine_4_"
# 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["input_scaling"], start_parameters["dend_tau"]])
initial_simplex = create_init_simples(x0, search_scale=2)
error_weights = (0, 1, 5, 1, 2, 0, 0)
fmin = minimize(fun=self.cost_function_with_fixed_adaption_with_dend_tau_no_noise,
args=(self.tau_a, self.delta_a, error_weights), x0=x0, method="Nelder-Mead",
options={"initial_simplex": initial_simplex, "xatol": 0.001})
res_parameters = fmin.x
print_comparision_cell_model(cell_data, self.base_model.get_parameters())
self.counter = 0
x0 = np.array([res_parameters[0], res_parameters[1], self.tau_a,
self.delta_a, res_parameters[2]])
initial_simplex = create_init_simples(x0, search_scale=2)
error_weights = (0, 0, 0, 0, 0, 4, 2)
fmin = minimize(fun=self.cost_function_all_without_noise,
args=(error_weights,), x0=x0, method="Nelder-Mead",
options={"initial_simplex": initial_simplex, "xatol": 0.001})
res_parameters = fmin.x
print(fmin)
print_comparision_cell_model(cell_data, self.base_model.get_parameters())
self.counter = 0
x0 = np.array([res_parameters[0],
res_parameters[1], self.tau_a,
self.delta_a, res_parameters[2]])
initial_simplex = create_init_simples(x0, search_scale=2)
error_weights = (0, 0, 1, 0, 0, 5, 2)
fmin = minimize(fun=self.cost_function_all_without_noise,
args=(error_weights,), x0=x0, method="Nelder-Mead",
options={"initial_simplex": initial_simplex, "xatol": 0.001})
res_parameters = self.base_model.get_parameters()
print_comparision_cell_model(cell_data, self.base_model.get_parameters())
# noise_strength = 0.03
# self.counter = 0
# x0 = np.array([res_parameters["mem_tau"], noise_strength,
# res_parameters["input_scaling"], res_parameters["tau_a"],
# res_parameters["delta_a"], res_parameters["dend_tau"]])
# initial_simplex = create_init_simples(x0, search_scale=2)
# error_weights = (0, 2, 2, 1, 1, 3, 2)
# fmin = minimize(fun=self.cost_function_all,
# args=(error_weights,), x0=x0, method="Nelder-Mead",
# options={"initial_simplex": initial_simplex, "xatol": 0.001})
# res_parameters = self.base_model.get_parameters()
return fmin, self.base_model.get_parameters()
def fit_model(self, x0=None, initial_simplex=None, fit_adaption=False):
self.counter = 0
@ -220,6 +327,27 @@ class Fitter:
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
@ -276,6 +404,27 @@ class Fitter:
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
@ -329,7 +478,7 @@ class Fitter:
error = sum(error_list)
self.counter += 1
if self.counter % 200 == 0:
if self.counter % 200 == 0 and False: # TODO currently shut off!
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(
@ -354,8 +503,8 @@ def calculate_f_values_error(fit, reference):
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 = 0
error += abs((fit[i] - reference[i]) / (reference[i] + constant))
constant = 50
error += abs((fit[i] - reference[i])+constant) / (abs(reference[i]) + constant)
norm_error = error / len(reference)

40
Sounds.py Normal file
View File

@ -0,0 +1,40 @@
from pyaudio import PyAudio
import numpy as np
BITRATE = 20000
LENGTH = 0.5
def play_finished_sound():
global BITRATE
global LENGTH
frequency = 261
num_of_frames = int(BITRATE*LENGTH)
frames = np.arange(0, num_of_frames, 1)
wave_data_numeric = np.sin(frames / ((BITRATE / frequency) / np.pi)) * 127 + 128
wave_data_numeric = wave_data_numeric.astype(int)
wave_data_chr = "".join([chr(x) for x in wave_data_numeric])
rest_frames = num_of_frames % BITRATE
rest = [chr(128)]*rest_frames
wave_data_chr.join(rest)
p = PyAudio()
stream = p.open(
format=p.get_format_from_width(1),
channels=1,
rate=BITRATE,
output=True,
)
stream.write(wave_data_chr)
stream.stop_stream()
stream.close()
p.terminate()
if __name__ == '__main__':
play_finished_sound()

4
helperFunctions.py
View File

@ -427,6 +427,10 @@ def detect_f_zero_in_frequency_trace(time, frequency, stimulus_start, sampling_i
stimulus_start = stimulus_start - time[0] # time start is generally != 0 and != delay
freq_before = frequency[int(buffer/sampling_interval):int((stimulus_start - buffer) / sampling_interval)]
if len(freq_before) < 3:
print("mäh")
min_before = min(freq_before)
max_before = max(freq_before)
mean_before = np.mean(freq_before)

6
models/LIFACnoise.py
View File

@ -227,6 +227,12 @@ class LifacNoiseModel(AbstractModel):
# if c == contrasts[0] or c == contrasts[-1]:
# plt.plot(frequency)
# plt.show()
if len(time) == 0 or time[0] >= stim_start or len(spiketimes) < 5:
f_infinities.append(0)
f_zeros.append(0)
f_baselines.append(0)
continue
f_inf = hF.detect_f_infinity_in_freq_trace(time, frequency, stim_start, stim_duration, sampling_interval)
f_infinities.append(f_inf)