diff --git a/models/LIFACnoise.py b/models/LIFACnoise.py index 906f1c6..b738e5f 100644 --- a/models/LIFACnoise.py +++ b/models/LIFACnoise.py @@ -4,7 +4,9 @@ from models.AbstractModel import AbstractModel import numpy as np import functions as fu from numba import jit -import time +import helperFunctions as hF +from stimuli.SinusAmplitudeModulation import SinusAmplitudeModulationStimulus +from scipy.optimize import curve_fit class LifacNoiseModel(AbstractModel): @@ -64,7 +66,7 @@ class LifacNoiseModel(AbstractModel): return output_voltage, spiketimes - def simulate_fast(self, stimulus: AbstractStimulus, total_time_s): + def simulate_fast(self, stimulus: AbstractStimulus, total_time_s, time_start=0): v_zero = self.parameters["v_zero"] a_zero = self.parameters["a_zero"] @@ -77,10 +79,12 @@ class LifacNoiseModel(AbstractModel): mem_tau = self.parameters["mem_tau"] noise_strength = self.parameters["noise_strength"] - stimulus_array = stimulus.as_array(total_time_s, step_size) + stimulus_array = stimulus.as_array(time_start, total_time_s, step_size) + rectified_stimulus = rectify_stimulus_array(stimulus_array) parameters = np.array([v_zero, a_zero, step_size, threshold, v_base, delta_a, tau_a, v_offset, mem_tau, noise_strength]) - voltage_trace, adaption, spiketimes = simulate_fast(stimulus_array, total_time_s, parameters) + + voltage_trace, adaption, spiketimes = simulate_fast(rectified_stimulus, total_time_s, parameters) self.stimulus = stimulus self.voltage_trace = voltage_trace @@ -146,35 +150,53 @@ class LifacNoiseModel(AbstractModel): def get_model_copy(self): return LifacNoiseModel(self.parameters) - def calculate_baseline_markers(self, stimulus_freq=750): + def calculate_baseline_markers(self, base_stimulus_freq, max_lag=1): """ 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 = SinusAmplitudeModulationStimulus(base_stimulus_freq, 0, 0) + _, spiketimes = self.simulate_fast(base_stimulus, 30) + baseline_freq = hF.mean_freq_of_spiketimes_after_time_x(spiketimes, 5) + relative_spiketimes = np.array([s % (1 / base_stimulus_freq) for s in spiketimes]) + eod_durations = np.full((len(spiketimes)), 1 / base_stimulus_freq) + vector_strength = hF.__vector_strength__(relative_spiketimes, eod_durations) + serial_correlation = hF.calculate_serial_correlation(np.array(spiketimes), max_lag) + return baseline_freq, vector_strength, serial_correlation - pass - - def calculate_fi_markers(self, contrasts, ): + def calculate_fi_markers(self, contrasts, base_freq, modulation_frequency): """ calculates the fi markers f_infinity, f_infinity_slope for given contrasts based on simulated 2 seconds for each contrast - :return: + :return: f_inf_values_list, f_inf_slope """ + f_infinities = [] + for contrast in contrasts: + stimulus = SinusAmplitudeModulationStimulus(base_freq, contrast, modulation_frequency) + _, spiketimes = self.simulate_fast(stimulus, 0.5) + + if len(spiketimes) < 2: + f_infinities.append(0) + else: + f_infinity = hF.mean_freq_of_spiketimes_after_time_x(spiketimes, 0.4) + f_infinities.append(f_infinity) + + popt, pcov = curve_fit(fu.line, contrasts, f_infinities, maxfev=10000) -def stimulus_to_numpy_array(stimulus: AbstractStimulus, total_time_s, step_size): - total_time_points = int(total_time_s * 1000 / step_size) - stimulus_values = np.zeros(total_time_points) - for idx in range(len(stimulus_values)): - # rectified input: - stimulus_values[idx] = fu.rectify(stimulus.value_at_time_in_ms(step_size*idx)) + f_infinities_slope = popt[0] - return stimulus_values + return f_infinities, f_infinities_slope + + +@jit(nopython=True) +def rectify_stimulus_array(stimulus_array:np.ndarray): + return np.array([x if x > 0 else 0 for x in stimulus_array]) @jit(nopython=True)