P-unit_model/run_Fitter.py
2020-09-21 00:05:45 +02:00

232 lines
9.6 KiB
Python

from models.LIFACnoise import LifacNoiseModel
from CellData import CellData
from Baseline import get_baseline_class
from FiCurve import get_fi_curve_class
from Fitter import Fitter
from ModelFit import get_best_fit, ModelFit
import time
import os
import argparse
import numpy as np
from helperFunctions import plot_errors
import multiprocessing as mp
SAVE_DIRECTORY = "./results/ref_and_tau/no_ref_period/"
SAVE_DIRECTORY_BEST = "./results/ref_and_tau/nrp_best/"
# SAVE_DIRECTORY = "./results/ref_and_tau/no_dend_tau/"
# SAVE_DIRECTORY_BEST = "./results/ref_and_tau/ndt_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)
def main():
parser = argparse.ArgumentParser()
parser.add_argument("--cell", help="folder (with .dat files) containing the cell data")
args = parser.parse_args()
if args.cell is not None:
cell_data = CellData(args.cell)
cell_name = os.path.split(cell_data.get_data_path())[-1]
if os.path.exists(SAVE_DIRECTORY + "/" + cell_name + "/"):
print(cell_name, "already done")
return
start_parameters = [p for p in iget_start_parameters()]
fit_cell_parallel(cell_data, start_parameters)
quit()
# test_single_cell("data/invivo/2010-11-08-al-invivo-1/")
test_single_cell("data/final/2012-12-21-am-invivo-1/")
# start_parameters = [p for p in iget_start_parameters()]
# cell_data = CellData("data/invivo_bursty/2014-03-19-ae-invivo-1/")
# fit_cell_parallel(cell_data, start_parameters)
def test_single_cell(path):
cell_data = CellData(path)
start_parameters = [p for p in iget_start_parameters()]
for i, p in enumerate(start_parameters):
fitter = Fitter()
fitter.set_data_reference_values(cell_data)
fmin, res_par = fitter.fit_routine_no_dend_tau_and_no_ref_period(p, ERROR_WEIGHTS)
cell_path = 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(i, sum(error))
save_fitting_run_info(cell_data, res_par, p, plot=True, save_path=save_path)
print("Done with start parameters {}".format(str(i)))
def fit_cell_base(parameters):
# parameter = (cell_data, start_parameter_index, start_parameter)
time1 = time.time()
fitter = Fitter()
fitter.set_data_reference_values(parameters[0])
fmin, res_par = fitter.fit_routine(parameters[2], ERROR_WEIGHTS)
cell_data = parameters[0]
cell_name = os.path.split(cell_data.get_data_path())[-1]
error = fitter.calculate_errors(model=LifacNoiseModel(res_par))
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_name +
"\n and start parameters ({:}): {:.2f}s thread time".format(parameters[1]+1, time2 - time1) +
"\n error: {:.2f}".format(sum(error)))
def fit_cell_parallel(cell_data, start_parameters):
cell_path = os.path.basename(cell_data.get_data_path())
save_path_cell = os.path.join(SAVE_DIRECTORY, cell_data.get_cell_name())
print(cell_path)
core_count = mp.cpu_count()
pool = mp.Pool(core_count - 1)
parameters = []
for i, p in enumerate(start_parameters):
parameters.append((cell_data, i, p))
time1 = time.time()
pool.map(fit_cell_base, parameters)
time2 = time.time()
print("Time taken for all start parameters ({:}): {:.2f}s".format(len(start_parameters), time2-time1))
del pool
del cell_data
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():
# mem_tau, input_scaling, noise_strength, dend_tau,
# expand by tau_a, delta_a ?
mem_tau_list = [0.001]
input_scaling_list = [80]
noise_strength_list = [0.01]
dend_tau_list = [0.002]
delta_a_list = [0.01, 0.03, 0.065]
tau_a_list = [0.02, 0.04]
ref_time_list = [0.00065, 0.0012]
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:
for delta_a in delta_a_list:
for tau_a in tau_a_list:
for ref_time in ref_time_list:
yield {"mem_tau": mem_tau, "input_scaling": input_scaling,
"noise_strength": noise_strength, "dend_tau": dend_tau,
"delta_a": delta_a, "tau_a": tau_a, "refractory_period": ref_time}
def save_fitting_run_info(cell_data, parameters, start_parameters, plot=False, save_path=None):
if save_path is not None:
if not os.path.exists(save_path):
os.makedirs(save_path)
if save_path is None:
return
with open(save_path + "parameters_info.txt", "w") as file:
file.writelines(["start_parameters:\t" + str(start_parameters),
"\nfinal_parameters:\t" + str(parameters)])
model = LifacNoiseModel(parameters)
eod_frequency = cell_data.get_eod_frequency()
data_baseline = get_baseline_class(cell_data)
c_bf = data_baseline.get_baseline_frequency()
c_vs = data_baseline.get_vector_strength()
c_sc = data_baseline.get_serial_correlation(1)
c_cv = data_baseline.get_coefficient_of_variation()
c_burst = data_baseline.get_burstiness()
data_fi_curve = get_fi_curve_class(cell_data, cell_data.get_fi_contrasts(), save_dir=cell_data.get_data_path())
c_f_inf_slope = data_fi_curve.get_f_inf_slope()
c_f_inf_values = data_fi_curve.f_inf_frequencies
c_f_zero_slope = data_fi_curve.get_f_zero_fit_slope_at_straight()
c_f_zero_values = data_fi_curve.f_zero_frequencies
if c_f_inf_slope < 0:
contrasts = np.array(cell_data.get_fi_contrasts()) * -1
else:
contrasts = np.array(cell_data.get_fi_contrasts())
model_baseline = get_baseline_class(model, eod_frequency, trials=15)
m_bf = model_baseline.get_baseline_frequency()
m_vs = model_baseline.get_vector_strength()
m_sc = model_baseline.get_serial_correlation(1)
m_cv = model_baseline.get_coefficient_of_variation()
m_burst = model_baseline.get_burstiness()
model_ficurve = get_fi_curve_class(model, cell_data.get_fi_contrasts(), eod_frequency, trials=60)
m_f_infinities = model_ficurve.get_f_inf_frequencies()
m_f_zeros = model_ficurve.get_f_zero_frequencies()
m_f_infinities_slope = model_ficurve.get_f_inf_slope()
m_f_zero_slope = model_ficurve.get_f_zero_fit_slope_at_straight()
np.save(os.path.join(save_path, "model_fi_inf_values.npy"), np.array(m_f_infinities))
np.save(os.path.join(save_path, "cell_fi_inf_values.npy"), np.array(c_f_inf_values))
np.save(os.path.join(save_path, "model_fi_zero_values.npy"), np.array(m_f_zeros))
np.save(os.path.join(save_path, "cell_fi_zero_values.npy"), np.array(c_f_zero_values))
with open(os.path.join(save_path, "cell_data_path.txt"), "w") as f:
path = cell_data.get_data_path() + "\n"
f.write(path)
if c_f_inf_slope < 0:
model_ficurve.stimulus_values = contrasts * -1
# print("EOD-frequency: {:.2f}".format(cell_data.get_eod_frequency()))
# 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("cv: cell - {:.2f} vs model {:.2f}".format(c_cv, m_cv))
# print("f_inf_slope: cell - {:.2f} vs model {:.2f}".format(c_f_inf_slope, m_f_infinities_slope))
# print("f infinity values:\n cell -", c_f_inf_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 -", m_f_zeros)
if save_path is not None:
with open(save_path + "value_comparision.tsv", 'w') as value_file:
value_file.write("Variable\tCell\tModel\n")
value_file.write("baseline_frequency\t{:.2f}\t{:.2f}\n".format(c_bf, m_bf))
value_file.write("vector_strength\t{:.2f}\t{:.2f}\n".format(c_vs, m_vs))
value_file.write("serial_correlation\t{:.2f}\t{:.2f}\n".format(c_sc[0], m_sc[0]))
value_file.write("Burstiness\t{:.2f}\t{:.2f}\n".format(c_burst, m_burst))
value_file.write("coefficient_of_variation\t{:.2f}\t{:.2f}\n".format(c_cv, m_cv))
value_file.write("f_inf_slope\t{:.2f}\t{:.2f}\n".format(c_f_inf_slope, m_f_infinities_slope))
value_file.write("f_zero_slope\t{:.2f}\t{:.2f}\n".format(c_f_zero_slope, m_f_zero_slope))
if plot:
# plot model images
model_baseline.plot_baseline(save_path)
model_baseline.plot_interspike_interval_histogram(save_path)
model_baseline.plot_isi_histogram_comparision(data_baseline.get_interspike_intervals(),
model_baseline.get_interspike_intervals(), save_path)
model_baseline.plot_serial_correlation(6, save_path)
model_ficurve.plot_fi_curve(save_path)
model_ficurve.plot_fi_curve_comparision(data_fi_curve, model_ficurve, save_path)
if __name__ == '__main__':
main()