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

separate FIcurve class, change fitter results structure

Browse Source
This commit is contained in:
a.ott 2020-05-12 14:14:55 +02:00
parent c1db288f97
commit 1f82f06209
7 changed files with 429 additions and 303 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
AdaptionCurrent.py
View File

@ -1,5 +1,5 @@
from FiCurve import FICurve from FiCurve import FICurve, get_fi_curve_class
from CellData import CellData from CellData import CellData
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
@ -13,7 +13,7 @@ class Adaption:
def __init__(self, cell_data: CellData, fi_curve: FICurve = None): def __init__(self, cell_data: CellData, fi_curve: FICurve = None):
self.cell_data = cell_data self.cell_data = cell_data
if fi_curve is None: if fi_curve is None:
self.fi_curve = FICurve(cell_data) self.fi_curve = get_fi_curve_class(cell_data, cell_data.get_fi_contrasts())
else: else:
self.fi_curve = fi_curve self.fi_curve = fi_curve
@ -51,7 +51,7 @@ class Adaption:
# start the actual fit: # start the actual fit:
try: try:
p0 = (self.fi_curve.f_zeros[i], tau, self.fi_curve.f_infinities[i]) p0 = (self.fi_curve.f_zero_frequencies[i], tau, self.fi_curve.f_inf_frequencies[i])
popt, pcov = curve_fit(fu.exponential_function, x_values, y_values, popt, pcov = curve_fit(fu.exponential_function, x_values, y_values,
p0=p0, maxfev=10000, bounds=([-np.inf, 0, -np.inf], [np.inf, np.inf, np.inf])) p0=p0, maxfev=10000, bounds=([-np.inf, 0, -np.inf], [np.inf, np.inf, np.inf]))
except RuntimeError: except RuntimeError:
@ -67,10 +67,10 @@ class Adaption:
self.exponential_fit_vars.append(popt) self.exponential_fit_vars.append(popt)
def __approximate_tau_for_exponential_fit(self, x_values, y_values, mean_freq_idx): def __approximate_tau_for_exponential_fit(self, x_values, y_values, mean_freq_idx):
if self.fi_curve.f_infinities[mean_freq_idx] < self.fi_curve.f_baselines[mean_freq_idx] * 0.95: if self.fi_curve.f_inf_frequencies[mean_freq_idx] < self.fi_curve.f_baseline_frequencies[mean_freq_idx] * 0.95:
test_val = [y > 0.65 * self.fi_curve.f_infinities[mean_freq_idx] for y in y_values] test_val = [y > 0.65 * self.fi_curve.f_inf_frequencies[mean_freq_idx] for y in y_values]
else: else:
test_val = [y < 0.65 * self.fi_curve.f_zeros[mean_freq_idx] for y in y_values] test_val = [y < 0.65 * self.fi_curve.f_zero_frequencies[mean_freq_idx] for y in y_values]
try: try:
idx = test_val.index(True) idx = test_val.index(True)
@ -85,13 +85,13 @@ class Adaption:
def __find_start_idx_for_exponential_fit(self, mean_freq_idx): def __find_start_idx_for_exponential_fit(self, mean_freq_idx):
time_axes = self.cell_data.get_time_axes_mean_frequencies()[mean_freq_idx] time_axes = self.cell_data.get_time_axes_mean_frequencies()[mean_freq_idx]
stimulus_start_idx = int((self.cell_data.get_stimulus_start() + time_axes[0]) / self.cell_data.get_sampling_interval()) stimulus_start_idx = int((self.cell_data.get_stimulus_start() + time_axes[0]) / self.cell_data.get_sampling_interval())
if self.fi_curve.f_infinities[mean_freq_idx] > self.fi_curve.f_baselines[mean_freq_idx] * 1.1: if self.fi_curve.f_inf_frequencies[mean_freq_idx] > self.fi_curve.f_baseline_frequencies[mean_freq_idx] * 1.1:
# start setting starting variables for the fit # start setting starting variables for the fit
# search for the start_index by searching for the max # search for the start_index by searching for the max
j = 0 j = 0
while True: while True:
try: try:
if self.cell_data.get_mean_isi_frequencies()[mean_freq_idx][stimulus_start_idx + j] == self.fi_curve.f_zeros[mean_freq_idx]: if self.cell_data.get_mean_isi_frequencies()[mean_freq_idx][stimulus_start_idx + j] == self.fi_curve.f_zero_frequencies[mean_freq_idx]:
start_idx = stimulus_start_idx + j start_idx = stimulus_start_idx + j
break break
except IndexError as e: except IndexError as e:
@ -99,7 +99,7 @@ class Adaption:
j += 1 j += 1
elif self.fi_curve.f_infinities[mean_freq_idx] < self.fi_curve.f_baselines[mean_freq_idx] * 0.9: elif self.fi_curve.f_inf_frequencies[mean_freq_idx] < self.fi_curve.f_baseline_frequencies[mean_freq_idx] * 0.9:
# start setting starting variables for the fit # start setting starting variables for the fit
# search for start by finding the end of the minimum # search for start by finding the end of the minimum
found_min = False found_min = False
@ -107,10 +107,10 @@ class Adaption:
nothing_to_fit = False nothing_to_fit = False
while True: while True:
if not found_min: if not found_min:
if self.cell_data.get_mean_isi_frequencies()[mean_freq_idx][stimulus_start_idx + j] == self.fi_curve.f_zeros[mean_freq_idx]: if self.cell_data.get_mean_isi_frequencies()[mean_freq_idx][stimulus_start_idx + j] == self.fi_curve.f_zero_frequencies[mean_freq_idx]:
found_min = True found_min = True
else: else:
if self.cell_data.get_mean_isi_frequencies()[mean_freq_idx][stimulus_start_idx + j + 1] > self.fi_curve.f_zeros[mean_freq_idx]: if self.cell_data.get_mean_isi_frequencies()[mean_freq_idx][stimulus_start_idx + j + 1] > self.fi_curve.f_zero_frequencies[mean_freq_idx]:
start_idx = stimulus_start_idx + j start_idx = stimulus_start_idx + j
break break
if j > 0.1 / self.cell_data.get_sampling_interval(): if j > 0.1 / self.cell_data.get_sampling_interval():
@ -133,7 +133,7 @@ class Adaption:
continue continue
tau_effs.append(self.exponential_fit_vars[i][1]) tau_effs.append(self.exponential_fit_vars[i][1])
f_infinity_slope = self.fi_curve.get_f_infinity_slope() f_infinity_slope = self.fi_curve.get_f_inf_slope()
# --- old way to calculate with the fi slope at middle of the fi curve # --- old way to calculate with the fi slope at middle of the fi curve
# fi_curve_slope = self.fi_curve.get_fi_curve_slope_of_straight() # fi_curve_slope = self.fi_curve.get_fi_curve_slope_of_straight()
# self.tau_real = np.median(tau_effs) * (fi_curve_slope / f_infinity_slope) # self.tau_real = np.median(tau_effs) * (fi_curve_slope / f_infinity_slope)
@ -174,7 +174,7 @@ class Adaption:
vars = self.exponential_fit_vars[i] vars = self.exponential_fit_vars[i]
fit_x = np.arange(0, 0.4, self.cell_data.get_sampling_interval()) fit_x = np.arange(0, 0.4, self.cell_data.get_sampling_interval())
plt.plot(fit_x, [fu.exponential_function(x, vars[0], vars[1], vars[2]) for x in fit_x]) plt.plot(fit_x, [fu.exponential_function(x, vars[0], vars[1], vars[2]) for x in fit_x])
plt.ylim([0, max(self.fi_curve.f_zeros[i], self.fi_curve.f_baselines[i])*1.1]) plt.ylim([0, max(self.fi_curve.f_zero_frequencies[i], self.fi_curve.f_baseline_frequencies[i])*1.1])
plt.xlabel("Time [s]") plt.xlabel("Time [s]")
plt.ylabel("Frequency [Hz]") plt.ylabel("Frequency [Hz]")

100
Baseline.py
View File

@ -31,15 +31,50 @@ class Baseline:
def get_interspike_intervals(self): def get_interspike_intervals(self):
raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS") raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS")
def plot_baseline(self, save_path=None): def plot_baseline(self, save_path=None, time_length=0.2):
raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS")
@staticmethod
def plot_baseline_given_data(time, eod, v1, spiketimes, sampling_interval, save_path=None, time_length=0.2):
""" """
plots the stimulus / eod, together with the v1, spiketimes and frequency plots the stimulus / eod, together with the v1, spiketimes and frequency
:return: :return:
""" """
raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS") length_data_points = int(time_length / sampling_interval)
start_idx = int(len(time) * 0.5 - length_data_points * 0.5)
start_idx = start_idx if start_idx >= 0 else 0
end_idx = int(len(time) * 0.5 + length_data_points * 0.5) + 1
end_idx = end_idx if end_idx <= len(time) else len(time)
spiketimes = np.array(spiketimes)
spiketimes_part = spiketimes[(spiketimes >= time[start_idx]) & (spiketimes < time[end_idx])]
def plot_inter_spike_interval_histogram(self, save_path=None): fig, axes = plt.subplots(3, 1, sharex="col", figsize=(12, 8))
isi = self.get_interspike_intervals() * 1000 # change unit to milliseconds fig.suptitle("Baseline middle part ({:.2f} seconds)".format(time_length))
axes[0].plot(time[start_idx:end_idx], eod[start_idx:end_idx])
axes[0].set_ylabel("Stimulus [mV]")
max_v1 = max(v1[start_idx:end_idx])
axes[1].plot(time[start_idx:end_idx], v1[start_idx:end_idx])
axes[1].plot(spiketimes_part, [max_v1 for _ in range(len(spiketimes_part))],
'o', color='orange')
axes[1].set_ylabel("V1-Trace [mV]")
t, f = hF.calculate_time_and_frequency_trace(spiketimes_part, sampling_interval)
axes[2].plot(t, f)
axes[2].set_ylabel("ISI-Frequency [Hz]")
axes[2].set_xlabel("Time [s]")
if save_path is not None:
plt.savefig(save_path + "baseline.png")
else:
plt.show()
plt.close()
def plot_interspike_interval_histogram(self, save_path=None):
isi = np.array(self.get_interspike_intervals()) * 1000 # change unit to milliseconds
maximum = max(isi) maximum = max(isi)
bins = np.arange(0, maximum * 1.01, 0.1) bins = np.arange(0, maximum * 1.01, 0.1)
@ -49,7 +84,7 @@ class Baseline:
plt.hist(isi, bins=bins) plt.hist(isi, bins=bins)
if save_path is not None: if save_path is not None:
plt.savefig(save_path) plt.savefig(save_path + "isi-histogram.png")
else: else:
plt.show() plt.show()
@ -60,10 +95,10 @@ class Baseline:
plt.xlabel("Lag") plt.xlabel("Lag")
plt.ylabel("Correlation") plt.ylabel("Correlation")
plt.ylim((-1, 1)) plt.ylim((-1, 1))
plt.plot(np.arange(1,max_lag+1, 1), self.get_serial_correlation(max_lag)) plt.plot(np.arange(1, max_lag+1, 1), self.get_serial_correlation(max_lag))
if save_path is not None: if save_path is not None:
plt.savefig(save_path) plt.savefig(save_path + "serial_correlation.png")
else: else:
plt.show() plt.show()
@ -83,7 +118,7 @@ class BaselineCellData(Baseline):
delay = self.data.get_delay() delay = self.data.get_delay()
sampling_interval = self.data.get_sampling_interval() sampling_interval = self.data.get_sampling_interval()
if delay < 0.1: if delay < 0.1:
warn("FICurve:__calculate_f_baseline__(): Quite short delay at the start.") warn("BaselineCellData:get_baseline_Frequency(): Quite short delay at the start.")
idx_start = int(0.025 / sampling_interval) idx_start = int(0.025 / sampling_interval)
idx_end = int((delay - 0.025) / sampling_interval) idx_end = int((delay - 0.025) / sampling_interval)
@ -136,30 +171,16 @@ class BaselineCellData(Baseline):
return isis return isis
def plot_baseline(self, save_path=None): def plot_baseline(self, save_path=None, time_length=0.2):
# eod, v1, spiketimes, frequency # eod, v1, spiketimes, frequency
times = self.data.get_base_traces(self.data.TIME) time = self.data.get_base_traces(self.data.TIME)[0]
eods = self.data.get_base_traces(self.data.EOD) eod = self.data.get_base_traces(self.data.EOD)[0]
v1_traces = self.data.get_base_traces(self.data.V1) v1_trace = self.data.get_base_traces(self.data.V1)[0]
spiketimes = self.data.get_base_spikes() spiketimes = self.data.get_base_spikes()[0]
fig, axes = plt.subplots(4, 1, sharex='col')
for i in range(len(times)): self.plot_baseline_given_data(time, eod, v1_trace, spiketimes,
axes[0].plot(times[i], eods[i]) self.data.get_sampling_interval(), save_path, time_length)
axes[1].plot(times[i], v1_traces[i])
axes[2].plot(spiketimes[i], [1 for i in range(len(spiketimes[i]))], 'o')
t, f = hF.calculate_time_and_frequency_trace(spiketimes[i], self.data.get_sampling_interval())
axes[3].plot(t, f)
if save_path is not None:
plt.savefig(save_path)
else:
plt.show()
plt.close()
class BaselineModel(Baseline): class BaselineModel(Baseline):
@ -200,24 +221,9 @@ class BaselineModel(Baseline):
def get_interspike_intervals(self): def get_interspike_intervals(self):
return np.diff(self.spiketimes) return np.diff(self.spiketimes)
def plot_baseline(self, save_path=None): def plot_baseline(self, save_path=None, time_length=0.2):
# eod, v1, spiketimes, frequency self.plot_baseline_given_data(self.time, self.eod, self.v1, self.spiketimes,
self.model.get_sampling_interval(), save_path, time_length)
fig, axes = plt.subplots(4, 1, sharex="col")
axes[0].plot(self.time, self.eod)
axes[1].plot(self.time, self.v1)
axes[2].plot(self.spiketimes, [1 for i in range(len(self.spiketimes))], 'o')
t, f = hF.calculate_time_and_frequency_trace(self.spiketimes, self.model.get_sampling_interval())
axes[3].plot(t, f)
if save_path is not None:
plt.savefig(save_path)
else:
plt.show()
plt.close()
def get_baseline_class(data, eod_freq=None) -> Baseline: def get_baseline_class(data, eod_freq=None) -> Baseline:

4
DataParserFactory.py
View File

@ -3,6 +3,7 @@ from os.path import isdir, exists
from warnings import warn from warnings import warn
import pyrelacs.DataLoader as Dl import pyrelacs.DataLoader as Dl
from models.AbstractModel import AbstractModel from models.AbstractModel import AbstractModel
import numpy as np
UNKNOWN = -1 UNKNOWN = -1
DAT_FORMAT = 0 DAT_FORMAT = 0
@ -82,11 +83,12 @@ class DatParser(AbstractParser):
return self.__get_traces__("BaselineActivity") return self.__get_traces__("BaselineActivity")
def get_baseline_spiketimes(self): def get_baseline_spiketimes(self):
# TODO change: reading from file -> detect from v1 trace
spiketimes = [] spiketimes = []
warn("Spiketimes don't fit time-wise to the baseline traces. Causes different vector strength angle per recording.") warn("Spiketimes don't fit time-wise to the baseline traces. Causes different vector strength angle per recording.")
for metadata, key, data in Dl.iload(self.baseline_file): for metadata, key, data in Dl.iload(self.baseline_file):
spikes = data[:, 0] spikes = np.array(data[:, 0]) / 1000 # timestamps are saved in ms -> conversion to seconds
spiketimes.append(spikes) spiketimes.append(spikes)
return spiketimes return spiketimes

387
FiCurve.py
View File

@ -1,5 +1,7 @@
from CellData import CellData from CellData import CellData
from models.LIFACnoise import LifacNoiseModel
from stimuli.SinusoidalStepStimulus import SinusoidalStepStimulus
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from warnings import warn from warnings import warn
@ -9,91 +11,66 @@ import helperFunctions as hF
class FICurve: class FICurve:
def __init__(self, cell_data: CellData, contrast: bool = True): def __init__(self, stimulus_values):
self.cell_data = cell_data self.stimulus_values = stimulus_values
self.using_contrast = contrast
if contrast:
self.stimulus_value = cell_data.get_fi_contrasts()
else:
self.stimulus_value = cell_data.get_fi_intensities()
self.f_zeros = [] self.f_baseline_frequencies = []
self.f_infinities = [] self.f_inf_frequencies = []
self.f_baselines = [] self.f_zero_frequencies = []
# f_max, f_min, k, x_zero
self.boltzmann_fit_vars = []
# increase, offset # increase, offset
self.f_infinity_fit = [] self.f_inf_fit = []
# f_max, f_min, k, x_zero
self.f_zero_fit = []
self.all_calculate_frequency_points() self.initialize()
self.f_infinity_fit = hF.fit_clipped_line(self.stimulus_value, self.f_infinities)
self.boltzmann_fit_vars = hF.fit_boltzmann(self.stimulus_value, self.f_zeros)
def all_calculate_frequency_points(self): def initialize(self):
mean_frequencies = self.cell_data.get_mean_isi_frequencies() self.calculate_all_frequency_points()
time_axes = self.cell_data.get_time_axes_mean_frequencies() self.f_inf_fit = hF.fit_clipped_line(self.stimulus_values, self.f_inf_frequencies)
stimulus_start = self.cell_data.get_stimulus_start() self.f_zero_fit = hF.fit_boltzmann(self.stimulus_values, self.f_zero_frequencies)
stimulus_duration = self.cell_data.get_stimulus_duration()
sampling_interval = self.cell_data.get_sampling_interval()
if len(mean_frequencies) == 0: def calculate_all_frequency_points(self):
warn("FICurve:all_calculate_frequency_points(): mean_frequencies is empty.\n" raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS")
"Was all_calculate_mean_isi_frequencies already called?")
for i in range(len(mean_frequencies)): def get_f_baseline_frequencies(self):
return self.f_baseline_frequencies
if time_axes[i][0] > self.cell_data.get_stimulus_start(): def get_f_inf_frequencies(self):
warn("TODO: Deal with to strongly cut frequency traces in cell data! ") return self.f_inf_frequencies
self.f_zeros.append(-1)
self.f_baselines.append(-1)
self.f_infinities.append(-1)
continue
f_zero = hF.detect_f_zero_in_frequency_trace(time_axes[i], mean_frequencies[i], def get_f_zero_frequencies(self):
stimulus_start, sampling_interval) return self.f_zero_frequencies
self.f_zeros.append(f_zero)
f_baseline = hF.detect_f_baseline_in_freq_trace(time_axes[i], mean_frequencies[i],
stimulus_start, sampling_interval)
self.f_baselines.append(f_baseline)
f_infinity = hF.detect_f_infinity_in_freq_trace(time_axes[i], mean_frequencies[i],
stimulus_start, stimulus_duration, sampling_interval)
self.f_infinities.append(f_infinity)
def get_f_zero_inverse_at_frequency(self, frequency): def get_f_inf_slope(self):
b_vars = self.boltzmann_fit_vars if len(self.f_inf_fit) > 0:
return fu.inverse_full_boltzmann(frequency, b_vars[0], b_vars[1], b_vars[2], b_vars[3]) return self.f_inf_fit[0]
def get_f_infinity_frequency_at_stimulus_value(self, stimulus_value):
infty_vars = self.f_infinity_fit
return fu.clipped_line(stimulus_value, infty_vars[0], infty_vars[1])
def get_f_infinity_slope(self): def get_f_zero_fit_slope_at_straight(self):
return self.f_infinity_fit[0] fit_vars = self.f_zero_fit
return fu.full_boltzmann_straight_slope(fit_vars[0], fit_vars[1], fit_vars[2], fit_vars[3])
def get_fi_curve_slope_at(self, stimulus_value): def get_f_zero_fit_slope_at_stimulus_value(self, stimulus_value):
fit_vars = self.boltzmann_fit_vars fit_vars = self.f_zero_fit
return fu.derivative_full_boltzmann(stimulus_value, fit_vars[0], fit_vars[1], fit_vars[2], fit_vars[3]) return fu.derivative_full_boltzmann(stimulus_value, fit_vars[0], fit_vars[1], fit_vars[2], fit_vars[3])
def get_fi_curve_slope_of_straight(self): def get_f_inf_frequency_at_stimulus_value(self, stimulus_value):
fit_vars = self.boltzmann_fit_vars return fu.clipped_line(stimulus_value, self.f_inf_fit[0], self.f_inf_fit[1])
return fu.full_boltzmann_straight_slope(fit_vars[0], fit_vars[1], fit_vars[2], fit_vars[3])
def get_f_zero_and_f_inf_intersection(self): def get_f_zero_and_f_inf_intersection(self):
x_values = np.arange(min(self.stimulus_value), max(self.stimulus_value), 0.0001) x_values = np.arange(min(self.stimulus_values), max(self.stimulus_values), 0.0001)
fit_vars = self.boltzmann_fit_vars fit_vars = self.f_zero_fit
f_zero = fu.full_boltzmann(x_values, fit_vars[0], fit_vars[1], fit_vars[2], fit_vars[3]) f_zero = fu.full_boltzmann(x_values, fit_vars[0], fit_vars[1], fit_vars[2], fit_vars[3])
f_inf = fu.clipped_line(x_values, self.f_infinity_fit[0], self.f_infinity_fit[1]) f_inf = fu.clipped_line(x_values, self.f_inf_fit[0], self.f_inf_fit[1])
intersection_indicies = np.argwhere(np.diff(np.sign(f_zero - f_inf))).flatten() intersection_indicies = np.argwhere(np.diff(np.sign(f_zero - f_inf))).flatten()
# print("fi-curve calc intersection:", intersection_indicies, x_values[intersection_indicies]) # print("fi-curve calc intersection:", intersection_indicies, x_values[intersection_indicies])
if len(intersection_indicies) > 1: if len(intersection_indicies) > 1:
f_baseline = np.median(self.f_baselines) f_baseline = np.median(self.f_baseline_frequencies)
best_dist = np.inf best_dist = np.inf
best_idx = -1 best_idx = -1
for idx in intersection_indicies: for idx in intersection_indicies:
dist = abs(fu.clipped_line(x_values[idx], self.f_infinity_fit[0], self.f_infinity_fit[1]) - f_baseline) dist = abs(fu.clipped_line(x_values[idx], self.f_inf_fit[0], self.f_inf_fit[1]) - f_baseline)
if dist < best_dist: if dist < best_dist:
best_dist = dist best_dist = dist
best_idx = idx best_idx = idx
@ -107,57 +84,291 @@ class FICurve:
def get_fi_curve_slope_at_f_zero_intersection(self): def get_fi_curve_slope_at_f_zero_intersection(self):
x = self.get_f_zero_and_f_inf_intersection() x = self.get_f_zero_and_f_inf_intersection()
fit_vars = self.boltzmann_fit_vars fit_vars = self.f_zero_fit
return fu.derivative_full_boltzmann(x, fit_vars[0], fit_vars[1], fit_vars[2], fit_vars[3]) return fu.derivative_full_boltzmann(x, fit_vars[0], fit_vars[1], fit_vars[2], fit_vars[3])
def plot_fi_curve(self, savepath: str = None, comp_f_baselines=None, comp_f_zeros=None, comp_f_infs=None): def plot_fi_curve(self, save_path=None):
min_x = min(self.stimulus_value) min_x = min(self.stimulus_values)
max_x = max(self.stimulus_value) max_x = max(self.stimulus_values)
step = (max_x - min_x) / 5000 step = (max_x - min_x) / 5000
x_values = np.arange(min_x, max_x, step) x_values = np.arange(min_x, max_x, step)
plt.plot(self.stimulus_value, self.f_baselines, color='blue', label='f_base') plt.plot(self.stimulus_values, self.f_baseline_frequencies, color='blue', label='f_base')
if comp_f_baselines is not None:
plt.plot(self.stimulus_value, comp_f_baselines, 'o', color='skyblue', label='comp_values base')
plt.plot(self.stimulus_value, self.f_infinities, 'o', color='green', label='f_inf') plt.plot(self.stimulus_values, self.f_inf_frequencies, 'o', color='green', label='f_inf')
plt.plot(x_values, [fu.clipped_line(x, self.f_infinity_fit[0], self.f_infinity_fit[1]) for x in x_values], plt.plot(x_values, [fu.clipped_line(x, self.f_inf_fit[0], self.f_inf_fit[1]) for x in x_values],
color='darkgreen', label='f_inf_fit') color='darkgreen', label='f_inf_fit')
if comp_f_infs is not None:
plt.plot(self.stimulus_value, comp_f_infs, 'o', color='lime', label='comp values f_inf')
plt.plot(self.stimulus_value, self.f_zeros, 'o', color='orange', label='f_zero') plt.plot(self.stimulus_values, self.f_zero_frequencies, 'o', color='orange', label='f_zero')
popt = self.boltzmann_fit_vars popt = self.f_zero_fit
plt.plot(x_values, [fu.full_boltzmann(x, popt[0], popt[1], popt[2], popt[3]) for x in x_values], plt.plot(x_values, [fu.full_boltzmann(x, popt[0], popt[1], popt[2], popt[3]) for x in x_values],
color='red', label='f_0_fit') color='red', label='f_0_fit')
if comp_f_zeros is not None:
plt.plot(self.stimulus_value, comp_f_zeros, 'o', color='wheat', label='comp_values f_zero')
plt.legend() plt.legend()
plt.ylabel("Frequency [Hz]") plt.ylabel("Frequency [Hz]")
if self.using_contrast: plt.xlabel("Stimulus value")
plt.xlabel("Stimulus contrast")
if save_path is None:
plt.show()
else: else:
plt.xlabel("Stimulus intensity [mv]") print("save")
if savepath is None: plt.savefig(save_path + "fi_curve.png")
plt.close()
@staticmethod
def plot_fi_curve_comparision(data_fi_curve, model_fi_curve, save_path=None):
min_x = min(min(data_fi_curve.stimulus_values), min(model_fi_curve.stimulus_values))
max_x = max(max(data_fi_curve.stimulus_values), max(model_fi_curve.stimulus_values))
step = (max_x - min_x) / 5000
x_values = np.arange(min_x, max_x+step, step)
fig, axes = plt.subplots(1, 3, sharex="all", sharey='all', figsize=(15, 6))
# plot baseline
data_origin = (data_fi_curve, model_fi_curve)
f_base_color = ("blue", "deepskyblue")
f_inf_color = ("green", "limegreen")
f_zero_color = ("red", "orange")
for i in range(2):
axes[i].plot(data_origin[i].stimulus_values, data_origin[i].get_f_baseline_frequencies(),
color=f_base_color[i], label='f_base')
axes[i].plot(data_origin[i].stimulus_values, data_origin[i].get_f_inf_frequencies(),
'o', color=f_inf_color[i], label='f_inf')
y_values = [fu.clipped_line(x, data_origin[i].f_inf_fit[0], data_origin[i].f_inf_fit[1]) for x in x_values]
axes[i].plot(x_values, y_values, color=f_inf_color[i], label='f_inf_fit')
axes[i].plot(data_origin[i].stimulus_values, data_origin[i].get_f_zero_frequencies(),
'o', color=f_zero_color[i], label='f_zero')
popt = data_origin[i].f_zero_fit
axes[i].plot(x_values, [fu.full_boltzmann(x, popt[0], popt[1], popt[2], popt[3]) for x in x_values],
color=f_zero_color[i], label='f_0_fit')
axes[i].set_xlabel("Stimulus value - contrast")
axes[i].legend()
axes[0].set_title("cell")
axes[0].set_ylabel("Frequency [Hz]")
axes[1].set_title("model")
median_baseline = np.median(data_fi_curve.get_f_baseline_frequencies())
axes[2].plot((min_x, max_x), (median_baseline, median_baseline), color=f_base_color[0], label="cell med base")
axes[2].plot(model_fi_curve.stimulus_values, model_fi_curve.get_f_baseline_frequencies(),
'o', color=f_base_color[1], label='model base')
y_values = [fu.clipped_line(x, data_fi_curve.f_inf_fit[0], data_fi_curve.f_inf_fit[1]) for x in x_values]
axes[2].plot(x_values, y_values, color=f_inf_color[0], label='f_inf_fit cell')
axes[2].plot(model_fi_curve.stimulus_values, model_fi_curve.get_f_inf_frequencies(),
'o', color=f_inf_color[1], label='f_inf model')
popt = data_fi_curve.f_zero_fit
axes[2].plot(x_values, [fu.full_boltzmann(x, popt[0], popt[1], popt[2], popt[3]) for x in x_values],
color=f_zero_color[0], label='f_0_fit cell')
axes[2].plot(model_fi_curve.stimulus_values, model_fi_curve.get_f_zero_frequencies(),
'o', color=f_zero_color[1], label='f_zero model')
axes[2].set_title("cell model comparision")
axes[2].set_xlabel("Stimulus value - contrast")
axes[2].legend()
if save_path is None:
plt.show() plt.show()
else: else:
print("save") print("save")
plt.savefig(savepath + "fi_curve.png") plt.savefig(save_path + "fi_curve_comparision.png")
plt.close() plt.close()
def plot_f_point_detections(self): def plot_f_point_detections(self, save_path=None):
raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS")
class FICurveCellData(FICurve):
def __init__(self, cell_data: CellData, stimulus_values):
self.cell_data = cell_data
super().__init__(stimulus_values)
def calculate_all_frequency_points(self):
mean_frequencies = self.cell_data.get_mean_isi_frequencies()
time_axes = self.cell_data.get_time_axes_mean_frequencies()
stimulus_start = self.cell_data.get_stimulus_start()
stimulus_duration = self.cell_data.get_stimulus_duration()
sampling_interval = self.cell_data.get_sampling_interval()
if len(mean_frequencies) == 0:
warn("FICurve:all_calculate_frequency_points(): mean_frequencies is empty.\n"
"Was all_calculate_mean_isi_frequencies already called?")
for i in range(len(mean_frequencies)):
if time_axes[i][0] > self.cell_data.get_stimulus_start():
raise ValueError("TODO: Deal with to strongly cut frequency traces in cell data! ")
# self.f_zero_frequencies.append(-1)
# self.f_baseline_frequencies.append(-1)
# self.f_inf_frequencies.append(-1)
# continue
f_zero = hF.detect_f_zero_in_frequency_trace(time_axes[i], mean_frequencies[i],
stimulus_start, sampling_interval)
self.f_zero_frequencies.append(f_zero)
f_baseline = hF.detect_f_baseline_in_freq_trace(time_axes[i], mean_frequencies[i],
stimulus_start, sampling_interval)
self.f_baseline_frequencies.append(f_baseline)
f_infinity = hF.detect_f_infinity_in_freq_trace(time_axes[i], mean_frequencies[i],
stimulus_start, stimulus_duration, sampling_interval)
self.f_inf_frequencies.append(f_infinity)
def get_f_zero_inverse_at_frequency(self, frequency):
# UNUSED
b_vars = self.f_zero_fit
return fu.inverse_full_boltzmann(frequency, b_vars[0], b_vars[1], b_vars[2], b_vars[3])
def get_f_infinity_frequency_at_stimulus_value(self, stimulus_value):
# UNUSED
infty_vars = self.f_inf_fit
return fu.clipped_line(stimulus_value, infty_vars[0], infty_vars[1])
# def get_fi_curve_slope_at(self, stimulus_value):
# fit_vars = self.f_zero_fit
# return fu.derivative_full_boltzmann(stimulus_value, fit_vars[0], fit_vars[1], fit_vars[2], fit_vars[3])
#
# def get_fi_curve_slope_of_straight(self):
# fit_vars = self.f_zero_fit
# return fu.full_boltzmann_straight_slope(fit_vars[0], fit_vars[1], fit_vars[2], fit_vars[3])
# def get_f_zero_and_f_inf_intersection(self):
# x_values = np.arange(min(self.stimulus_values), max(self.stimulus_values), 0.0001)
# fit_vars = self.f_zero_fit
# f_zero = fu.full_boltzmann(x_values, fit_vars[0], fit_vars[1], fit_vars[2], fit_vars[3])
# f_inf = fu.clipped_line(x_values, self.f_inf_fit[0], self.f_inf_fit[1])
#
# intersection_indicies = np.argwhere(np.diff(np.sign(f_zero - f_inf))).flatten()
# # print("fi-curve calc intersection:", intersection_indicies, x_values[intersection_indicies])
# if len(intersection_indicies) > 1:
# f_baseline = np.median(self.f_baseline_frequencies)
# best_dist = np.inf
# best_idx = -1
# for idx in intersection_indicies:
# dist = abs(fu.clipped_line(x_values[idx], self.f_inf_fit[0], self.f_inf_fit[1]) - f_baseline)
# if dist < best_dist:
# best_dist = dist
# best_idx = idx
#
# return x_values[best_idx]
#
# elif len(intersection_indicies) == 0:
# raise ValueError("No intersection found!")
# else:
# return x_values[intersection_indicies[0]]
# def get_fi_curve_slope_at_f_zero_intersection(self):
# x = self.get_f_zero_and_f_inf_intersection()
# fit_vars = self.f_zero_fit
# return fu.derivative_full_boltzmann(x, fit_vars[0], fit_vars[1], fit_vars[2], fit_vars[3])
# def plot_fi_curve(self, savepath: str = None, comp_f_baselines=None, comp_f_zeros=None, comp_f_infs=None):
# min_x = min(self.stimulus_values)
# max_x = max(self.stimulus_values)
# step = (max_x - min_x) / 5000
# x_values = np.arange(min_x, max_x, step)
#
# plt.plot(self.stimulus_values, self.f_baseline_frequencies, color='blue', label='f_base')
# if comp_f_baselines is not None:
# plt.plot(self.stimulus_values, comp_f_baselines, 'o', color='skyblue', label='comp_values base')
#
# plt.plot(self.stimulus_values, self.f_inf_frequencies, 'o', color='green', label='f_inf')
# plt.plot(x_values, [fu.clipped_line(x, self.f_inf_fit[0], self.f_inf_fit[1]) for x in x_values],
# color='darkgreen', label='f_inf_fit')
# if comp_f_infs is not None:
# plt.plot(self.stimulus_values, comp_f_infs, 'o', color='lime', label='comp values f_inf')
#
# plt.plot(self.stimulus_values, self.f_zero_frequencies, 'o', color='orange', label='f_zero')
# popt = self.f_zero_fit
# plt.plot(x_values, [fu.full_boltzmann(x, popt[0], popt[1], popt[2], popt[3]) for x in x_values],
# color='red', label='f_0_fit')
# if comp_f_zeros is not None:
# plt.plot(self.stimulus_values, comp_f_zeros, 'o', color='wheat', label='comp_values f_zero')
#
# plt.legend()
# plt.ylabel("Frequency [Hz]")
# plt.xlabel("Stimulus value")
#
# if savepath is None:
# plt.show()
# else:
# print("save")
# plt.savefig(savepath + "fi_curve.png")
# plt.close()
def plot_f_point_detections(self, save_path=None):
mean_frequencies = np.array(self.cell_data.get_mean_isi_frequencies()) mean_frequencies = np.array(self.cell_data.get_mean_isi_frequencies())
time_axes = self.cell_data.get_time_axes_mean_frequencies() time_axes = self.cell_data.get_time_axes_mean_frequencies()
for i in range(len(mean_frequencies)): for i in range(len(mean_frequencies)):
fig, axes = plt.subplots(1, 1, sharex="all") fig, axes = plt.subplots(1, 1, sharex="all")
axes.plot(time_axes[i], mean_frequencies[i], label="voltage") axes.plot(time_axes[i], mean_frequencies[i], label="voltage")
axes.plot((time_axes[i][0], time_axes[i][-1]), (self.f_zeros[i], self.f_zeros[i]), label="f_zero") axes.plot((time_axes[i][0], time_axes[i][-1]), (self.f_zero_frequencies[i], self.f_zero_frequencies[i]), label="f_zero")
axes.plot((time_axes[i][0], time_axes[i][-1]), (self.f_infinities[i], self.f_infinities[i]), '--', label="f_inf") axes.plot((time_axes[i][0], time_axes[i][-1]), (self.f_inf_frequencies[i], self.f_inf_frequencies[i]), '--', label="f_inf")
axes.plot((time_axes[i][0], time_axes[i][-1]), (self.f_baselines[i], self.f_baselines[i]), label="f_base") axes.plot((time_axes[i][0], time_axes[i][-1]), (self.f_baseline_frequencies[i], self.f_baseline_frequencies[i]), label="f_base")
axes.set_title(str(self.stimulus_value[i])) axes.set_title(str(self.stimulus_values[i]))
plt.legend() plt.legend()
plt.show() if save_path is None:
plt.show()
else:
plt.savefig(save_path + "GENERATE_NAMES.png")
plt.close()
class FICurveModel(FICurve):
def __init__(self, model, stimulus_values, eod_frequency):
self.eod_frequency = eod_frequency
self.model = model
super().__init__(stimulus_values)
def calculate_all_frequency_points(self):
stim_duration = 0.5
stim_start = 0.5
total_simulation_time = stim_duration + 2 * stim_start
sampling_interval = self.model.get_sampling_interval()
self.f_inf_frequencies = []
self.f_zero_frequencies = []
self.f_baseline_frequencies = []
for c in self.stimulus_values:
stimulus = SinusoidalStepStimulus(self.eod_frequency, c, stim_start, stim_duration)
_, spiketimes = self.model.simulate_fast(stimulus, total_simulation_time)
time, frequency = hF.calculate_time_and_frequency_trace(spiketimes, sampling_interval)
if len(spiketimes) < 10 or len(time) == 0 or min(time) > stim_start \
or max(time) < stim_start + stim_duration:
print("Too few spikes to calculate f_inf, f_0 and f_base")
self.f_inf_frequencies.append(0)
self.f_zero_frequencies.append(0)
self.f_baseline_frequencies.append(0)
continue
f_inf = hF.detect_f_infinity_in_freq_trace(time, frequency, stim_start, stim_duration, sampling_interval)
self.f_inf_frequencies.append(f_inf)
f_zero = hF.detect_f_zero_in_frequency_trace(time, frequency, stim_start, sampling_interval)
self.f_zero_frequencies.append(f_zero)
f_baseline = hF.detect_f_baseline_in_freq_trace(time, frequency, stim_start, sampling_interval)
self.f_baseline_frequencies.append(f_baseline)
def plot_f_point_detections(self, save_path=None):
raise NotImplementedError("TODO sorry... "
"The model version of the FiCurve class is still missing this implementation")
def get_fi_curve_class(data, stimulus_values, eod_freq=None) -> FICurve:
if isinstance(data, CellData):
return FICurveCellData(data, stimulus_values)
if isinstance(data, LifacNoiseModel):
if eod_freq is None:
raise ValueError("The FiCurveModel needs the eod variable to work")
return FICurveModel(data, stimulus_values, eod_freq)
raise ValueError("Unknown type: Cannot find corresponding Baseline class. Data was type:" + str(type(data)))

121
Fitter.py
View File

@ -3,7 +3,7 @@ from models.LIFACnoise import LifacNoiseModel
from stimuli.SinusoidalStepStimulus import SinusoidalStepStimulus from stimuli.SinusoidalStepStimulus import SinusoidalStepStimulus
from CellData import CellData, icelldata_of_dir from CellData import CellData, icelldata_of_dir
from Baseline import get_baseline_class from Baseline import get_baseline_class
from FiCurve import FICurve from FiCurve import get_fi_curve_class
from AdaptionCurrent import Adaption from AdaptionCurrent import Adaption
import helperFunctions as hF import helperFunctions as hF
import functions as fu import functions as fu
@ -40,19 +40,37 @@ def iget_start_parameters():
def run_with_real_data(): def run_with_real_data():
for cell_data in icelldata_of_dir("./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)
if not os.path.exists(results_path):
os.makedirs(results_path)
# plot cell images:
cell_save_path = results_path + "cell/"
if not os.path.exists(cell_save_path):
os.makedirs(cell_save_path)
data_baseline = get_baseline_class(cell_data)
data_baseline.plot_baseline(cell_save_path)
data_baseline.plot_interspike_interval_histogram(cell_save_path)
data_baseline.plot_serial_correlation(6, cell_save_path)
data_fi_curve = get_fi_curve_class(cell_data, cell_data.get_fi_contrasts())
data_fi_curve.plot_fi_curve(cell_save_path)
start_par_count = 0 start_par_count = 0
for start_parameters in iget_start_parameters(): for start_parameters in iget_start_parameters():
start_par_count += 1 start_par_count += 1
print("START PARAMETERS:", start_par_count) print("START PARAMETERS:", start_par_count)
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] + "/" parameter_set_path = results_path + "start_parameter_set_{}".format(start_par_count) + "/"
print("results at:", results_path)
results_path += "parameter_set_{}".format(start_par_count) + "/"
start_time = time.time() start_time = time.time()
fitter = Fitter() fitter = Fitter()
@ -61,54 +79,49 @@ def run_with_real_data():
print(fmin) print(fmin)
print(parameters) print(parameters)
end_time = time.time() end_time = time.time()
if not os.path.exists(parameter_set_path):
if not os.path.exists(results_path): os.makedirs(parameter_set_path)
os.makedirs(results_path) with open(parameter_set_path + "parameters_info.txt".format(start_par_count), "w") as file:
with open(results_path + "parameters_info.txt".format(start_par_count), "w") as file:
file.writelines(["start_parameters:\t" + str(start_parameters), file.writelines(["start_parameters:\t" + str(start_parameters),
"\nfinal_parameters:\t" + str(parameters), "\nfinal_parameters:\t" + str(parameters),
"\nfinal_fmin:\t" + str(fmin)]) "\nfinal_fmin:\t" + str(fmin)])
print('Fitting of cell took function took {:.3f} s'.format((end_time - start_time))) print('Fitting of cell took function took {:.3f} s'.format((end_time - start_time)))
# print(results_path) # print(results_path)
print_comparision_cell_model(cell_data, parameters, plot=True, savepath=results_path) print_comparision_cell_model(cell_data, data_baseline, data_fi_curve, parameters,
break plot=True, save_path=parameter_set_path)
# from Sounds import play_finished_sound # from Sounds import play_finished_sound
# play_finished_sound() # play_finished_sound()
pass pass
def print_comparision_cell_model(cell_data: CellData, parameters, plot=False, savepath=None): def print_comparision_cell_model(cell_data, data_baseline, data_fi_curve, parameters, plot=False, save_path=None):
res_model = LifacNoiseModel(parameters) model = LifacNoiseModel(parameters)
fi_curve = FICurve(cell_data) eod_frequency = cell_data.get_eod_frequency()
model_baseline = get_baseline_class(res_model, cell_data.get_eod_frequency()) model_baseline = get_baseline_class(model, eod_frequency)
m_bf = model_baseline.get_baseline_frequency() m_bf = model_baseline.get_baseline_frequency()
m_vs = model_baseline.get_vector_strength() m_vs = model_baseline.get_vector_strength()
m_sc = model_baseline.get_serial_correlation(1) m_sc = model_baseline.get_serial_correlation(1)
m_cv = model_baseline.get_coefficient_of_variation() m_cv = model_baseline.get_coefficient_of_variation()
_, m_f_zeros, m_f_infinities = res_model.calculate_fi_curve(fi_curve.stimulus_value, model_ficurve = get_fi_curve_class(model, cell_data.get_fi_contrasts(), eod_frequency)
cell_data.get_eod_frequency()) m_f_infinities = model_ficurve.get_f_inf_frequencies()
f_infinities_fit = hF.fit_clipped_line(fi_curve.stimulus_value, m_f_infinities) m_f_zeros = model_ficurve.get_f_zero_frequencies()
m_f_infinities_slope = f_infinities_fit[0] m_f_infinities_slope = model_ficurve.get_f_inf_slope()
m_f_zero_slope = model_ficurve.get_f_zero_fit_slope_at_straight()
f_zeros_fit = hF.fit_boltzmann(fi_curve.stimulus_value, m_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])
data_baseline = get_baseline_class(cell_data)
c_bf = data_baseline.get_baseline_frequency() c_bf = data_baseline.get_baseline_frequency()
c_vs = data_baseline.get_vector_strength() c_vs = data_baseline.get_vector_strength()
c_sc = data_baseline.get_serial_correlation(1) c_sc = data_baseline.get_serial_correlation(1)
c_cv = data_baseline.get_coefficient_of_variation() c_cv = data_baseline.get_coefficient_of_variation()
c_f_inf_slope = fi_curve.get_f_infinity_slope() c_f_inf_slope = data_fi_curve.get_f_inf_slope()
c_f_inf_values = fi_curve.f_infinities c_f_inf_values = data_fi_curve.f_inf_frequencies
c_f_zero_slope = fi_curve.get_fi_curve_slope_of_straight() c_f_zero_slope = data_fi_curve.get_f_zero_fit_slope_at_straight()
c_f_zero_values = fi_curve.f_zeros c_f_zero_values = data_fi_curve.f_zero_frequencies
print("EOD-frequency: {:.2f}".format(cell_data.get_eod_frequency())) print("EOD-frequency: {:.2f}".format(cell_data.get_eod_frequency()))
print("bf: cell - {:.2f} vs model {:.2f}".format(c_bf, m_bf)) 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("vs: cell - {:.2f} vs model {:.2f}".format(c_vs, m_vs))
@ -119,8 +132,8 @@ def print_comparision_cell_model(cell_data: CellData, parameters, plot=False, sa
print("f_zero_slope: cell - {:.2f} vs model {:.2f}".format(c_f_zero_slope, m_f_zero_slope)) 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) print("f zero values:\n cell -", c_f_zero_values, "\n model -", m_f_zeros)
if savepath is not None: if save_path is not None:
with open(savepath + "value_comparision.tsv", 'w') as value_file: with open(save_path + "value_comparision.tsv", 'w') as value_file:
value_file.write("Variable\tCell\tModel\n") 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("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("vector_strength\t{:.2f}\t{:.2f}\n".format(c_vs, m_vs))
@ -130,26 +143,13 @@ def print_comparision_cell_model(cell_data: CellData, parameters, plot=False, sa
value_file.write("f_zero_slope\t{:.2f}\t{:.2f}\n".format(c_f_zero_slope, m_f_zero_slope)) value_file.write("f_zero_slope\t{:.2f}\t{:.2f}\n".format(c_f_zero_slope, m_f_zero_slope))
if plot: if plot:
# plot cell images:
cell_save_path = savepath + "/cell/"
if not os.path.exists(cell_save_path):
os.makedirs(cell_save_path)
# data_baseline.plot_baseline(cell_save_path + "baseline.png")
data_baseline.plot_inter_spike_interval_histogram(cell_save_path + "isi-histogram.png")
data_baseline.plot_serial_correlation(6, cell_save_path + "serial_correlation.png")
# plot model images # plot model images
model_save_path = savepath + "/model/" model_baseline.plot_baseline(save_path)
if not os.path.exists(model_save_path): model_baseline.plot_interspike_interval_histogram(save_path)
os.makedirs(model_save_path) model_baseline.plot_serial_correlation(6, save_path)
# model_baseline.plot_baseline(model_save_path + "baseline.png")
model_baseline.plot_inter_spike_interval_histogram(model_save_path + "isi-histogram.png")
model_baseline.plot_serial_correlation(6, model_save_path + "serial_correlation.png")
if plot:
f_b, f_zero, f_inf = res_model.calculate_fi_curve(cell_data.get_fi_contrasts(), cell_data.get_eod_frequency())
fi_curve.plot_fi_curve(savepath=savepath, comp_f_baselines=f_b, comp_f_zeros=f_zero, comp_f_infs=f_inf) model_ficurve.plot_fi_curve(save_path)
model_ficurve.plot_fi_curve_comparision(data_fi_curve, model_ficurve, save_path)
class Fitter: class Fitter:
@ -196,14 +196,14 @@ class Fitter:
self.serial_correlation = data_baseline.get_serial_correlation(self.sc_max_lag) self.serial_correlation = data_baseline.get_serial_correlation(self.sc_max_lag)
self.coefficient_of_variation = data_baseline.get_coefficient_of_variation() self.coefficient_of_variation = data_baseline.get_coefficient_of_variation()
fi_curve = FICurve(data, contrast=True) fi_curve = get_fi_curve_class(data, data.get_fi_contrasts())
self.fi_contrasts = fi_curve.stimulus_value self.fi_contrasts = fi_curve.stimulus_values
self.f_inf_values = fi_curve.f_infinities self.f_inf_values = fi_curve.f_inf_frequencies
self.f_inf_slope = fi_curve.get_f_infinity_slope() self.f_inf_slope = fi_curve.get_f_inf_slope()
self.f_zero_values = fi_curve.f_zeros self.f_zero_values = fi_curve.f_zero_frequencies
self.f_zero_fit = fi_curve.boltzmann_fit_vars self.f_zero_fit = fi_curve.f_zero_fit
self.f_zero_slope = fi_curve.get_fi_curve_slope_of_straight() self.f_zero_slope = fi_curve.get_f_zero_fit_slope_at_straight()
# around 1/3 of the value at straight # 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()) # self.f_zero_slope = fi_curve.get_fi_curve_slope_at(fi_curve.get_f_zero_and_f_inf_intersection())
@ -216,7 +216,6 @@ class Fitter:
# 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_5(data, start_parameters) return self.fit_routine_5(data, start_parameters)
# return self.fit_model(fit_adaption=False)
def fit_routine_5(self, cell_data=None, start_parameters=None): def fit_routine_5(self, cell_data=None, start_parameters=None):
# [error_bf, error_vs, error_sc, error_cv, error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope] # [error_bf, error_vs, error_sc, error_cv, error_f_inf, error_f_inf_slope, error_f_zero, error_f_zero_slope]
@ -228,7 +227,7 @@ class Fitter:
x0 = np.array([start_parameters["mem_tau"], start_parameters["noise_strength"], x0 = np.array([start_parameters["mem_tau"], start_parameters["noise_strength"],
start_parameters["input_scaling"], self.tau_a, self.delta_a, start_parameters["dend_tau"]]) start_parameters["input_scaling"], self.tau_a, self.delta_a, start_parameters["dend_tau"]])
initial_simplex = create_init_simples(x0, search_scale=2) initial_simplex = create_init_simples(x0, search_scale=2)
error_weights = (0, 2, 3, 1, 1, 1, 0.5, 1) error_weights = (0, 2, 5, 1, 1, 1, 0.5, 1)
fmin = minimize(fun=self.cost_function_all, fmin = minimize(fun=self.cost_function_all,
args=(error_weights,), x0=x0, method="Nelder-Mead", args=(error_weights,), x0=x0, method="Nelder-Mead",
options={"initial_simplex": initial_simplex, "xatol": 0.001, "maxfev": 400, "maxiter": 400}) options={"initial_simplex": initial_simplex, "xatol": 0.001, "maxfev": 400, "maxiter": 400})

92
compareModelVsCell.py
View File

@ -1,92 +0,0 @@
from FiCurve import FICurve
from CellData import CellData
from models.LIFACnoise import LifacNoiseModel
import helperFunctions as hF
import functions as fu
import numpy as np
CELL_PATH = "data/2012-06-27-ah-invivo-1/"
# current parameters from fit with folus on improving model f0 response
# MODEL_PARAMETERS = {'step_size': 5e-05, 'mem_tau': 0.042542178602690675, 'v_base': 0, 'v_zero': 0, 'threshold': 1, 'v_offset': -111.328125, 'input_scaling': 388.9439738592248, 'delta_a': 0.05513136301255167, 'tau_a': 0.1017720885626184, 'a_zero': 2, 'noise_strength': 0.01740931732483443}
# parameters fit with fixed adaption focus on serial cor and f_inf_slope
MODEL_PARAMETERS = {'step_size': 5e-05, 'mem_tau': 0.03547683648372142, 'v_base': 0, 'v_zero': 0, 'threshold': 1, 'v_offset': -43.75, 'input_scaling': 162.97353975832954, 'delta_a': 0.024625305808413798, 'tau_a': 0.07632538029074364, 'a_zero': 2, 'noise_strength': 0.00408169739163286}
SAVE_PATH = "results/test/"
def main():
cell_data = CellData(CELL_PATH)
fi_curve = FICurve(cell_data)
model = LifacNoiseModel(MODEL_PARAMETERS)
print_characteristics(cell_data, fi_curve, model)
plot_fi_curve(cell_data, fi_curve, model)
def print_characteristics(cell_data: CellData, fi_curve: FICurve, model):
sc_max_lag = 1
eod_freq = cell_data.get_eod_frequency()
fi_contrasts = fi_curve.stimulus_value
cell_bf = cell_data.get_base_frequency()
cell_sc = cell_data.get_serial_correlation(sc_max_lag)
cell_vs = cell_data.get_vector_strength()
cell_f_inf_slope = fi_curve.get_f_infinity_slope()
cell_f_inf_values = fi_curve.f_infinities
cell_f_zero_slope = fi_curve.get_fi_curve_slope_of_straight()
cell_f_zero_values = fi_curve.f_zeros
# calculate Model characteristics:
baseline_freq, vector_strength, serial_correlation = model.calculate_baseline_markers(eod_freq, sc_max_lag)
f_baselines, f_zeros, f_infinities = model.calculate_fi_curve(fi_contrasts, eod_freq)
f_infinities_fit = hF.fit_clipped_line(fi_contrasts, f_infinities)
f_infinities_slope = f_infinities_fit[0]
f_zeros_fit = hF.fit_boltzmann(fi_contrasts, f_zeros)
f_zero_slope = fu.full_boltzmann_straight_slope(f_zeros_fit[0], f_zeros_fit[1], f_zeros_fit[2], f_zeros_fit[3])
print_value("Base frequency", cell_bf, baseline_freq)
for i in range(sc_max_lag):
print_value("Serial correlation lag {}".format(i+1), cell_sc[i], serial_correlation[i])
print_value("Vector strength", cell_vs, vector_strength)
print_value("f_inf slope", cell_f_inf_slope, f_infinities_slope)
print_f_values("f_inf values", cell_f_inf_values, f_infinities)
print_value("f_zero slope", cell_f_zero_slope, f_zero_slope)
print_f_values("f_zero values", cell_f_zero_values, f_zeros)
def print_value(name, cell_value, model_value):
print("{} - expected: {:.2f}, current: {:.2f}, %-error: {:.3f}".format(name,
cell_value, model_value, (model_value-cell_value)/cell_value))
def print_f_values(name, cell_values, model_values):
cell_values = np.array(cell_values)
model_values = np.array(model_values)
mean_perc_error = 0
for c,m in zip(cell_values, model_values):
mean_perc_error += (m-c)/ c
mean_perc_error = mean_perc_error/ len(cell_values)
print("{}:\nexpected:".format(name), np.around(cell_values), "\ncurrent: ", np.around(model_values),
"\nmean %-error: {:.3f}".format(mean_perc_error))
def plot_fi_curve(cell_data, fi_curve, model, save=False):
f_b, f_zero, f_inf = model.calculate_fi_curve(cell_data.get_fi_contrasts(), cell_data.get_eod_frequency())
if save:
fi_curve.plot_fi_curve(savepath=SAVE_PATH, comp_f_baselines=f_b, comp_f_zeros=f_zero, comp_f_infs=f_inf)
else:
fi_curve.plot_fi_curve(comp_f_baselines=f_b, comp_f_zeros=f_zero, comp_f_infs=f_inf)
if __name__ == '__main__':
main()

2
fit_lifacnoise.py
View File

@ -208,7 +208,7 @@ class Fitter:
self.coefficient_of_variation = data.get_coefficient_of_variation() self.coefficient_of_variation = data.get_coefficient_of_variation()
fi_curve = FICurve(data, contrast=True) fi_curve = FICurve(data, contrast=True)
self.fi_contrasts = fi_curve.stimulus_value self.fi_contrasts = fi_curve.stimulus_values
print("Fitter: fi-contrasts", self.fi_contrasts) print("Fitter: fi-contrasts", self.fi_contrasts)
self.f_infinities = fi_curve.f_infinities self.f_infinities = fi_curve.f_infinities
self.f_infinities_slope = fi_curve.get_f_infinity_slope() self.f_infinities_slope = fi_curve.get_f_infinity_slope()