From 11c37c9f2a41aaf0c7e7efc547e804d5e8201830 Mon Sep 17 00:00:00 2001 From: "a.ott" Date: Mon, 11 May 2020 13:56:56 +0200 Subject: [PATCH] seperate baseline attribute calculation --- Baseline.py | 231 +++++++++++++++++++++++++++++++++++++++++++ CellData.py | 26 ----- Fitter.py | 38 ++++--- models/LIFACnoise.py | 75 ++++++-------- 4 files changed, 283 insertions(+), 87 deletions(-) create mode 100644 Baseline.py diff --git a/Baseline.py b/Baseline.py new file mode 100644 index 0000000..fe17151 --- /dev/null +++ b/Baseline.py @@ -0,0 +1,231 @@ + +from CellData import CellData +from models.LIFACnoise import LifacNoiseModel +from stimuli.SinusoidalStepStimulus import SinusoidalStepStimulus +import helperFunctions as hF +import numpy as np +from warnings import warn +import matplotlib.pyplot as plt + + +class Baseline: + + def __init__(self): + self.baseline_frequency = -1 + self.serial_correlation = [] + self.vector_strength = -1 + self.coefficient_of_variation = -1 + + def get_baseline_frequency(self): + raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS") + + def get_serial_correlation(self, max_lag): + raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS") + + def get_vector_strength(self): + raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS") + + def get_coefficient_of_variation(self): + raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS") + + def get_interspike_intervals(self): + raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS") + + def plot_baseline(self, save_path=None): + """ + plots the stimulus / eod, together with the v1, spiketimes and frequency + :return: + """ + raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS") + + def plot_inter_spike_interval_histogram(self, save_path=None): + isi = self.get_interspike_intervals() * 1000 # change unit to milliseconds + maximum = max(isi) + bins = np.arange(0, maximum * 1.01, 0.1) + + plt.title('Baseline ISIs') + plt.xlabel('ISI in ms') + plt.ylabel('Count') + plt.hist(isi, bins=bins) + + if save_path is not None: + plt.savefig(save_path) + else: + plt.show() + + plt.close() + + def plot_serial_correlation(self, max_lag, save_path=None): + plt.title("Baseline Serial correlation") + plt.xlabel("Lag") + plt.ylabel("Correlation") + + plt.plot(np.arange(1,max_lag+1, 1), self.get_serial_correlation(max_lag)) + + if save_path is not None: + plt.savefig(save_path) + else: + plt.show() + + plt.close() + + +class BaselineCellData(Baseline): + + def __init__(self, cell_data: CellData): + super().__init__() + self.data = cell_data + + def get_baseline_frequency(self): + if self.baseline_frequency == -1: + base_freqs = [] + for freq in self.data.get_mean_isi_frequencies(): + delay = self.data.get_delay() + sampling_interval = self.data.get_sampling_interval() + if delay < 0.1: + warn("FICurve:__calculate_f_baseline__(): Quite short delay at the start.") + + idx_start = int(0.025 / sampling_interval) + idx_end = int((delay - 0.025) / sampling_interval) + base_freqs.append(np.mean(freq[idx_start:idx_end])) + + self.baseline_frequency = np.median(base_freqs) + + return self.baseline_frequency + + def get_vector_strength(self): + if self.vector_strength == -1: + times = self.data.get_base_traces(self.data.TIME) + eods = self.data.get_base_traces(self.data.EOD) + v1_traces = self.data.get_base_traces(self.data.V1) + self.vector_strength = hF.calculate_vector_strength_from_v1_trace(times, eods, v1_traces) + + return self.vector_strength + + def get_serial_correlation(self, max_lag): + if len(self.serial_correlation) != max_lag: + serial_cors = [] + for spiketimes in self.data.get_base_spikes(): + sc = hF.calculate_serial_correlation(spiketimes, max_lag) + serial_cors.append(sc) + serial_cors = np.array(serial_cors) + mean_sc = np.mean(serial_cors, axis=0) + + self.serial_correlation = mean_sc + return self.serial_correlation + + def get_coefficient_of_variation(self): + if self.coefficient_of_variation == -1: + spiketimes = self.data.get_base_spikes() + # CV (stddev of ISI divided by mean ISI (np.diff(spiketimes)) + cvs = [] + for st in spiketimes: + st = np.array(st) + cvs.append(hF.calculate_coefficient_of_variation(st)) + + self.coefficient_of_variation = np.mean(cvs) + return self.coefficient_of_variation + + def get_interspike_intervals(self): + spiketimes = self.data.get_base_spikes() + # CV (stddev of ISI divided by mean ISI (np.diff(spiketimes)) + isis = [] + for st in spiketimes: + st = np.array(st) + isis.extend(np.diff(st)) + + return isis + + def plot_baseline(self, save_path=None): + # eod, v1, spiketimes, frequency + + times = self.data.get_base_traces(self.data.TIME) + eods = self.data.get_base_traces(self.data.EOD) + v1_traces = self.data.get_base_traces(self.data.V1) + spiketimes = self.data.get_base_spikes() + + fig, axes = plt.subplots(4, 1, sharex="True") + + for i in range(len(times)): + axes[0].plot(times[i], eods[i]) + axes[1].plot(times[i], v1_traces[i]) + axes[2].plot(spiketimes, [1]*len(spiketimes), '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): + + simulation_time = 30 + + def __init__(self, model: LifacNoiseModel, eod_frequency): + super().__init__() + self.model = model + self.eod_frequency = eod_frequency + self.stimulus = SinusoidalStepStimulus(eod_frequency, 0) + self.v1, self.spiketimes = model.simulate_fast(self.stimulus, self.simulation_time) + self.eod = self.stimulus.as_array(0, self.simulation_time, model.get_sampling_interval()) + self.time = np.arange(0, self.simulation_time, model.get_sampling_interval()) + + def get_baseline_frequency(self): + if self.baseline_frequency == -1: + self.baseline_frequency = hF.calculate_mean_isi_freq(self.spiketimes) + + return self.baseline_frequency + + def get_vector_strength(self): + if self.vector_strength == -1: + self.vector_strength = hF.calculate_vector_strength_from_spiketimes(self.time, self.eod, self.spiketimes, + self.model.get_sampling_interval()) + return self.vector_strength + + def get_serial_correlation(self, max_lag): + if len(self.serial_correlation) != max_lag: + self.serial_correlation = hF.calculate_serial_correlation(self.spiketimes, max_lag) + return self.serial_correlation + + def get_coefficient_of_variation(self): + if self.coefficient_of_variation == -1: + self.coefficient_of_variation = hF.calculate_coefficient_of_variation(self.spiketimes) + return self.coefficient_of_variation + + def get_interspike_intervals(self): + return np.diff(self.spiketimes) + + def plot_baseline(self, save_path=None): + # eod, v1, spiketimes, frequency + + fig, axes = plt.subplots(4, 1, sharex="True") + + axes[0].plot(self.time, self.eod) + axes[1].plot(self.time, self.v1) + axes[2].plot(self.spiketimes, [1]*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: + if isinstance(data, CellData): + return BaselineCellData(data) + if isinstance(data, LifacNoiseModel): + if eod_freq is None: + raise ValueError("The EOD frequency is needed for the BaselineModel Class.") + return BaselineModel(data, eod_freq) + + raise ValueError("Unknown type: Cannot find corresponding Baseline class. data was type:" + str(type(data))) diff --git a/CellData.py b/CellData.py index 741885f..837ecd2 100644 --- a/CellData.py +++ b/CellData.py @@ -149,32 +149,6 @@ class CellData: def get_after_stimulus_duration(self) -> float: return self.recording_times[3] - def get_vector_strength(self): - times = self.get_base_traces(self.TIME) - eods = self.get_base_traces(self.EOD) - v1_traces = self.get_base_traces(self.V1) - return hf.calculate_vector_strength_from_v1_trace(times, eods, v1_traces) - - def get_serial_correlation(self, max_lag): - serial_cors = [] - for spiketimes in self.get_base_spikes(): - sc = hf.calculate_serial_correlation(spiketimes, max_lag) - serial_cors.append(sc) - serial_cors = np.array(serial_cors) - mean_sc = np.mean(serial_cors, axis=0) - - return mean_sc - - def get_coefficient_of_variation(self): - spiketimes = self.get_base_spikes() - # CV (stddev of ISI divided by mean ISI (np.diff(spiketimes)) - cvs = [] - for st in spiketimes: - st = np.array(st) - cvs.append(hf.calculate_coefficient_of_variation(st)) - - return np.mean(cvs) - def get_eod_frequency(self): eods = self.get_base_traces(self.EOD) sampling_interval = self.get_sampling_interval() diff --git a/Fitter.py b/Fitter.py index 7a28695..c0a358c 100644 --- a/Fitter.py +++ b/Fitter.py @@ -2,6 +2,7 @@ from models.LIFACnoise import LifacNoiseModel from stimuli.SinusoidalStepStimulus import SinusoidalStepStimulus from CellData import CellData, icelldata_of_dir +from Baseline import get_baseline_class from FiCurve import FICurve from AdaptionCurrent import Adaption import helperFunctions as hF @@ -52,7 +53,6 @@ def run_with_real_data(): 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) @@ -80,11 +80,16 @@ def run_with_real_data(): pass -def print_comparision_cell_model(cell_data, parameters, plot=False, savepath=None): +def print_comparision_cell_model(cell_data: CellData, parameters, plot=False, savepath=None): res_model = LifacNoiseModel(parameters) fi_curve = FICurve(cell_data) - m_bf, m_vs, m_sc, m_cv = res_model.calculate_baseline_markers(cell_data.get_eod_frequency()) + model_baseline = get_baseline_class(res_model, cell_data.get_eod_frequency()) + 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() + 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) @@ -93,10 +98,12 @@ def print_comparision_cell_model(cell_data, parameters, plot=False, savepath=Non 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) - c_cv = cell_data.get_coefficient_of_variation() + 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_f_slope = fi_curve.get_f_infinity_slope() c_f_values = fi_curve.f_infinities @@ -156,10 +163,11 @@ class Fitter: 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() - self.vector_strength = data.get_vector_strength() - self.serial_correlation = data.get_serial_correlation(self.sc_max_lag) - self.coefficient_of_variation = data.get_coefficient_of_variation() + data_baseline = get_baseline_class(data) + self.baseline_freq = data_baseline.get_baseline_frequency() + self.vector_strength = data_baseline.get_vector_strength() + self.serial_correlation = data_baseline.get_serial_correlation(self.sc_max_lag) + self.coefficient_of_variation = data_baseline.get_coefficient_of_variation() fi_curve = FICurve(data, contrast=True) self.fi_contrasts = fi_curve.stimulus_value @@ -347,9 +355,11 @@ class Fitter: return sum(error_list) def calculate_errors(self, error_weights=None): - baseline_freq, vector_strength, serial_correlation, coefficient_of_variation =\ - self.base_model.calculate_baseline_markers(self.eod_freq, self.sc_max_lag) - # print("baseline features calculated!") + model_baseline = get_baseline_class(self.base_model, self.eod_freq) + baseline_freq = model_baseline.get_baseline_frequency() + vector_strength = model_baseline.get_vector_strength() + serial_correlation = model_baseline.get_serial_correlation(self.sc_max_lag) + coefficient_of_variation = model_baseline.get_coefficient_of_variation() # f_infinities, f_infinities_slope = self.base_model.calculate_fi_markers(self.fi_contrasts, self.eod_freq) f_baselines, f_zeros, f_infinities = self.base_model.calculate_fi_curve(self.fi_contrasts, self.eod_freq) diff --git a/models/LIFACnoise.py b/models/LIFACnoise.py index c8ff10a..3010b29 100644 --- a/models/LIFACnoise.py +++ b/models/LIFACnoise.py @@ -1,4 +1,3 @@ - from stimuli.AbstractStimulus import AbstractStimulus from models.AbstractModel import AbstractModel import numpy as np @@ -10,6 +9,7 @@ from scipy.optimize import curve_fit from warnings import warn import matplotlib.pyplot as plt + class LifacNoiseModel(AbstractModel): # all times in milliseconds # possible mem_res: 100 * 1000000 exact value unknown in p-units @@ -56,7 +56,8 @@ class LifacNoiseModel(AbstractModel): for i in range(1, len(time), 1): time_point = time[i] # rectified input: - stimulus_strength = self._calculate_input_voltage_step(input_voltage[i-1], fu.rectify(stimulus.value_at_time_in_s(time_point))) + stimulus_strength = self._calculate_input_voltage_step(input_voltage[i - 1], + fu.rectify(stimulus.value_at_time_in_s(time_point))) v_next = self._calculate_voltage_step(current_v, stimulus_strength - current_a) a_next = self._calculate_adaption_step(current_a) @@ -97,7 +98,9 @@ class LifacNoiseModel(AbstractModel): def _calculate_input_voltage_step(self, current_i, rectified_input): # input_voltage[i] = input_voltage[i - 1] + (-input_voltage[i - 1] + rectified_stimulus_array[i] * input_scaling) / dend_tau - return current_i + ((-current_i + rectified_input * self.parameters["input_scaling"]) / self.parameters["dend_tau"]) * self.parameters["step_size"] + return current_i + ( + (-current_i + rectified_input * self.parameters["input_scaling"]) / self.parameters["dend_tau"]) * \ + self.parameters["step_size"] def simulate_fast(self, stimulus: AbstractStimulus, total_time_s, time_start=0): @@ -115,7 +118,9 @@ class LifacNoiseModel(AbstractModel): dend_tau = self.parameters["dend_tau"] rectified_stimulus = rectify_stimulus_array(stimulus.as_array(time_start, total_time_s, step_size)) - 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]) + 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]) voltage_trace, adaption, spiketimes, input_voltage = simulate_fast(rectified_stimulus, total_time_s, parameters) @@ -168,31 +173,6 @@ class LifacNoiseModel(AbstractModel): def get_model_copy(self): return LifacNoiseModel(self.parameters) - def calculate_baseline_markers(self, stimulus_freq, max_lag=1, simulation_time=30): - """ - calculates the baseline markers baseline frequency, vector strength and serial correlation - based on simulated 30 seconds with a standard Sinusoidal stimulus with the given frequency - - :return: baseline_freq, vs, sc - """ - base_stimulus = SinusoidalStepStimulus(stimulus_freq, 0) - _, spiketimes = self.simulate_fast(base_stimulus, simulation_time) - time_x = 5 - baseline_freq = hF.mean_freq_of_spiketimes_after_time_x(spiketimes, time_x) - - if baseline_freq < 1: - return baseline_freq, 0, [0]*max_lag - - else: - time_trace = np.arange(0, 30, self.get_sampling_interval()) - stimulus_array = base_stimulus.as_array(0, 30, self.get_sampling_interval()) - - vector_strength = hF.calculate_vector_strength_from_spiketimes(time_trace, stimulus_array, spiketimes, self.get_sampling_interval()) - serial_correlation = hF.calculate_serial_correlation(np.array(spiketimes), max_lag) - coeffient_of_variation = hF.calculate_coefficient_of_variation(np.array(spiketimes)) - - return baseline_freq, vector_strength, serial_correlation, coeffient_of_variation - def calculate_fi_markers(self, contrasts, stimulus_freq): """ calculates the fi markers f_infinity, f_infinity_slope for given contrasts @@ -204,9 +184,10 @@ class LifacNoiseModel(AbstractModel): f_infinities = [] for contrast in contrasts: stimulus = SinusoidalStepStimulus(stimulus_freq, contrast, stimulus_start, stimulus_duration) - _, spiketimes = self.simulate_fast(stimulus, stimulus_start*2+stimulus_duration) + _, spiketimes = self.simulate_fast(stimulus, stimulus_start * 2 + stimulus_duration) time, freq = hF.calculate_time_and_frequency_trace(spiketimes, self.get_sampling_interval()) - f_inf = hF.detect_f_infinity_in_freq_trace(time, freq, stimulus_start, stimulus_duration, self.get_sampling_interval()) + f_inf = hF.detect_f_infinity_in_freq_trace(time, freq, stimulus_start, stimulus_duration, + self.get_sampling_interval()) f_infinities.append(f_inf) popt = hF.fit_clipped_line(contrasts, f_infinities) @@ -217,7 +198,6 @@ class LifacNoiseModel(AbstractModel): def calculate_fi_curve(self, contrasts, stimulus_freq): - stim_duration = 0.5 stim_start = 0.5 total_simulation_time = stim_duration + 2 * stim_start @@ -240,7 +220,8 @@ class LifacNoiseModel(AbstractModel): # plt.plot(frequency) # plt.show() - if len(spiketimes) < 10 or len(time) == 0 or min(time) > stim_start or max(time) < stim_start+stim_duration: + 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") f_infinities.append(0) f_zeros.append(0) @@ -290,16 +271,18 @@ class LifacNoiseModel(AbstractModel): lower_bound = current_v_offset - v_search_step_size upper_bound = current_v_offset - return binary_search_base_freq(test_model, base_stimulus, goal_baseline_frequency, simulation_length, lower_bound, upper_bound, threshold) + return binary_search_base_freq(test_model, base_stimulus, goal_baseline_frequency, simulation_length, + lower_bound, upper_bound, threshold) -def binary_search_base_freq(model: LifacNoiseModel, base_stimulus, goal_frequency, simulation_length, lower_bound, upper_bound, threshold): +def binary_search_base_freq(model: LifacNoiseModel, base_stimulus, goal_frequency, simulation_length, lower_bound, + upper_bound, threshold): counter = 0 if threshold <= 0: raise ValueError("binary_search_base_freq() - LifacNoiseModel: threshold is not allowed to be negative!") while True: counter += 1 - middle = upper_bound - (upper_bound - lower_bound)/2 + middle = upper_bound - (upper_bound - lower_bound) / 2 frequency = test_v_offset(model, middle, base_stimulus, simulation_length) # print('{:.1f}, {:.1f}, {:.1f}, {:.1f} vs {:.1f} '.format(lower_bound, middle, upper_bound, frequency, goal_frequency)) @@ -311,10 +294,11 @@ def binary_search_base_freq(model: LifacNoiseModel, base_stimulus, goal_frequenc elif frequency > goal_frequency: upper_bound = middle else: - print('lower bound: {:.1f}, middle: {:.1f}, upper_bound: {:.1f}, frequency: {:.1f} vs goal: {:.1f} '.format(lower_bound, middle, upper_bound, frequency, goal_frequency)) + print('lower bound: {:.1f}, middle: {:.1f}, upper_bound: {:.1f}, frequency: {:.1f} vs goal: {:.1f} '.format( + lower_bound, middle, upper_bound, frequency, goal_frequency)) raise ValueError("binary_search_base_freq() - LifacNoiseModel: Goal frequency might be nan?") - if abs(upper_bound-lower_bound) < 0.0001: + if abs(upper_bound - lower_bound) < 0.0001: warn("Search was stopped no value was found!") return middle @@ -373,18 +357,15 @@ def simulate_fast(rectified_stimulus_array, total_time_s, parameters: np.ndarray noise_value = np.random.normal() noise = noise_strength * noise_value / np.sqrt(step_size) - input_voltage[i] = input_voltage[i - 1] + ((-input_voltage[i - 1] + rectified_stimulus_array[i]) / dend_tau) * 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 + input_voltage[i] = input_voltage[i - 1] + ( + (-input_voltage[i - 1] + rectified_stimulus_array[i]) / dend_tau) * 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 output_voltage[i] > threshold: output_voltage[i] = v_base - spiketimes.append(i*step_size) + spiketimes.append(i * step_size) adaption[i] += delta_a / tau_a return output_voltage, adaption, spiketimes, input_voltage - - - - -