From dacde1f1b62874f224d2a9f29bdff6dd79cddf65 Mon Sep 17 00:00:00 2001 From: "a.ott" Date: Mon, 16 Mar 2020 17:17:05 +0100 Subject: [PATCH] add tests for sinusoidalStep, add FiCurve plot possibilities (temp), add initial simplex generation --- FiCurve.py | 12 ++- fit_lifacnoise.py | 126 ++++++++++++++++-------- functions.py | 22 ++--- generalTests.py | 13 ++- helperFunctions.py | 54 ++++++++++ models/LIFACnoise.py | 60 ++++++++--- stimuli/SinusoidalStepStimulus.py | 2 +- unittests/testSinusoidalStepStimulus.py | 117 ++++++++++++++++++++++ 8 files changed, 339 insertions(+), 67 deletions(-) create mode 100644 unittests/testSinusoidalStepStimulus.py diff --git a/FiCurve.py b/FiCurve.py index e81cc76..0376809 100644 --- a/FiCurve.py +++ b/FiCurve.py @@ -202,22 +202,30 @@ class FICurve: fit_vars = self.boltzmann_fit_vars 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): + 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_value) max_x = max(self.stimulus_value) step = (max_x - min_x) / 5000 x_values = np.arange(min_x, max_x, step) plt.plot(self.stimulus_value, self.f_baselines, 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='lime', label='f_inf') + plt.plot(self.stimulus_value, self.f_infinities, '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], 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') popt = self.boltzmann_fit_vars 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_value, comp_f_zeros, 'o', color='wheat', label='comp_values f_zero') + plt.legend() plt.ylabel("Frequency [Hz]") diff --git a/fit_lifacnoise.py b/fit_lifacnoise.py index 7f080ed..9dc30aa 100644 --- a/fit_lifacnoise.py +++ b/fit_lifacnoise.py @@ -12,8 +12,8 @@ import matplotlib.pyplot as plt def main(): - # run_test_with_fixed_model() - # quit() + #run_test_with_fixed_model() + #quit() # fitter = Fitter() # fmin, params = fitter.fit_model_to_values(700, 1400, [-0.3], 1, [0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3], [1370, 1380, 1390, 1400, 1410, 1420, 1430], 100, 0.02, 0.01) @@ -34,8 +34,8 @@ def run_with_real_data(): print('Fitting of cell took function took {:.3f} s'.format((end_time - start_time))) res_model = LifacNoiseModel(parameters) - m_bf, m_vs, m_sc = res_model.calculate_baseline_markers(celldata.get_base_frequency()) - m_f_values, m_f_slope = res_model.calculate_fi_markers(celldata.get_fi_contrasts(), celldata.get_base_frequency(), 10) + m_bf, m_vs, m_sc = res_model.calculate_baseline_markers(celldata.get_eod_frequency()) + m_f_values, m_f_slope = res_model.calculate_fi_markers(celldata.get_fi_contrasts(), celldata.get_eod_frequency()) c_bf = celldata.get_base_frequency() c_vs = celldata.get_vector_strength() @@ -49,31 +49,64 @@ def run_with_real_data(): 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) + f_b, f_zero, f_inf = res_model.calculate_fi_curve(celldata.get_fi_contrasts(), celldata.get_eod_frequency()) + + fi_curve.plot_fi_curve(comp_f_baselines=f_b, comp_f_zeros=f_zero, comp_f_infs=f_inf) + break + pass def run_test_with_fixed_model(): - a_tau = 10 + a_tau = 0.1 a_delta = 0.08 - parameters = {'mem_tau': 5, 'delta_a': a_delta, 'input_scaling': 100, - 'v_offset': 80, 'threshold': 1, 'v_base': 0, 'step_size': 0.00005, 'tau_a': a_tau, - 'a_zero': 0, 'v_zero': 0, 'noise_strength': 0.5} + parameters = {"mem_tau": 0.015, + "v_base": 0, + "v_zero": 0, + "threshold": 1, + "v_offset": -10, + "input_scaling": 60, + "delta_a": a_delta, + "tau_a": a_tau, + "a_zero": 10, + "noise_strength": 0.05, + "step_size": 0.00005} model = LifacNoiseModel(parameters) eod_freq = 750 contrasts = np.arange(0.5, 1.51, 0.1) - modulation_freq = 10 baseline_freq, vector_strength, serial_correlation = model.calculate_baseline_markers(eod_freq) - f_infinities, f_infinities_slope = model.calculate_fi_markers(contrasts, eod_freq, modulation_freq) - + f_infinities, f_infinities_slope = model.calculate_fi_markers(contrasts, eod_freq) + print("Baseline freq:{:.2f}\nVector strength: {:.3f}\nSerial cor:".format(baseline_freq, vector_strength), serial_correlation) + print("FI-Curve\nSlope: {:.2f}\nValues:".format(f_infinities_slope), f_infinities) fitter = Fitter() - fmin, fit_parameters = fitter.fit_model_to_values(eod_freq, baseline_freq, serial_correlation, vector_strength, contrasts, f_infinities, f_infinities_slope, a_delta, a_tau) - print("calculated parameters:") - print(fit_parameters) + starting_conditions = [[0.05, 0.02, 50], [0.01, 0.08, 75], [0.02, 0.03, 70], + [0.02, 0.03, 70]] + for x0 in starting_conditions: + init_simplex = create_init_simples(x0) + + fmin, fit_parameters = fitter.fit_model_to_values(eod_freq, baseline_freq, serial_correlation, vector_strength, contrasts, f_infinities, f_infinities_slope, a_delta, a_tau, x0, init_simplex) + print(x0) + print(fmin) + print("calculated parameters:") + print(fit_parameters) + + print("ref parameters:") + print(parameters) - print("ref parameters:") - print(parameters) + +def create_init_simples(x0, search_scale=2): + dim = len(x0) + simplex = [[x0[0]/search_scale], [x0[0]*search_scale]] + for i in range(1, dim, 1): + for vertex in simplex: + vertex.append(x0[i]*search_scale) + new_vertex = list(x0[:i]) + new_vertex.append(x0[i]/search_scale) + simplex.append(new_vertex) + + return simplex class Fitter: @@ -118,30 +151,32 @@ class Fitter: self.f_infinities_slope = fi_curve.get_f_infinity_slope() f_zero_slope = fi_curve.get_fi_curve_slope_of_straight() - self.a_delta = (f_zero_slope / self.f_infinities_slope) + self.a_delta = (f_zero_slope / self.f_infinities_slope) / 1000 adaption = Adaption(data, fi_curve) - self.a_tau = adaption.get_tau_real() + self.a_tau = adaption.get_tau_real() *2 + print() # mem_tau, (threshold?), (v_offset), noise_strength, input_scaling - def cost_function(self, X, tau_a=0.001, delta_a=3, error_scaling=()): + def cost_function(self, X, tau_a=0.1, delta_a=0.08, error_scaling=()): # set model parameters to the given ones: self.model.set_variable("mem_tau", X[0]) self.model.set_variable("noise_strength", X[1]) self.model.set_variable("input_scaling", X[2]) - self.model.set_variable("tau_a", X[3]) - self.model.set_variable("delta_a", X[4]) - # self.model.set_variable("tau_a", tau_a) - # self.model.set_variable("delta_a", delta_a) + #self.model.set_variable("tau_a", X[3]) + #self.model.set_variable("delta_a", X[4]) + self.model.set_variable("tau_a", tau_a) + self.model.set_variable("delta_a", delta_a) # minimize the difference in baseline_freq first by fitting v_offset # v_offset = self.__fit_v_offset_to_baseline_frequency__() base_stimulus = SinusoidalStepStimulus(self.eod_freq, 0) - - v_offset = self.model.find_v_offset(self.baseline_freq, base_stimulus) + test_model = self.model.get_model_copy() + test_model.set_variable("noise_strength", 0) + v_offset = test_model.find_v_offset(self.baseline_freq, base_stimulus) self.model.set_variable("v_offset", v_offset) - baseline_freq, vector_strength, serial_correlation = self.model.calculate_baseline_markers(self.baseline_freq, self.sc_max_lag) + baseline_freq, vector_strength, serial_correlation = self.model.calculate_baseline_markers(self.eod_freq, self.sc_max_lag) # only eod with amplitude 1 and no modulation # _, spiketimes = self.model.simulate_fast(base_stimulus, 30) @@ -168,21 +203,22 @@ class Fitter: f_infinity = hF.mean_freq_of_spiketimes_after_time_x(spiketimes, 0.5) f_infinities.append(f_infinity) - popt, pcov = curve_fit(fu.line, self.fi_contrasts, f_infinities, maxfev=10000) + popt, pcov = curve_fit(fu.clipped_line, self.fi_contrasts, f_infinities, maxfev=10000) f_infinities_slope = popt[0] - if self.counter == 600: - bf, vs, sc = self.model.calculate_baseline_markers(self.eod_freq) error_bf = abs((baseline_freq - self.baseline_freq) / self.baseline_freq) error_vs = abs((vector_strength - self.vector_strength) / self.vector_strength) error_sc = abs((serial_correlation[0] - self.serial_correlation[0]) / self.serial_correlation[0]) - error_f_inf_slope = abs((f_infinities_slope - self.f_infinities_slope) / self.f_infinities_slope) + error_f_inf_slope = abs((f_infinities_slope - self.f_infinities_slope) / self.f_infinities_slope) * 4 #print("vs:", vector_strength, self.vector_strength) #print("sc", serial_correlation[0], self.serial_correlation[0]) #print("f slope:", f_infinities_slope, self.f_infinities_slope) error_f_inf = 0 for i in range(len(f_infinities)): - error_f_inf += abs((f_infinities[i] - self.f_infinities[i]) / f_infinities[i]) + f_inf = f_infinities[i] + if f_inf <= 0: + f_inf = 1 + error_f_inf += abs((f_infinities[i] - self.f_infinities[i]) / f_inf) error_f_inf = error_f_inf / len(f_infinities) @@ -190,16 +226,21 @@ class Fitter: errors = [error_bf, error_vs, error_sc, error_f_inf_slope, error_f_inf] self.counter += 1 - if self.counter % 10 == 0: - print("Cost function run times:", self.counter, "error sum:", sum(errors), errors) + if self.counter % 100 == 0: + print("\nCost function run times: {:}\n".format(self.counter), + "Total error: {:.4f}\n".format(sum(errors)), + "Baseline frequency - expected: {:.0f}, current: {:.0f}, error: {:.3f}\n".format(self.baseline_freq, baseline_freq, error_bf), + "Vector strength - expected: {:.2f}, current: {:.2f}, error: {:.3f}\n".format(self.vector_strength, vector_strength, error_vs), + "Serial correlation - expected: {:.2f}, current: {:.2f}, error: {:.3f}\n".format(self.serial_correlation[0], serial_correlation[0], error_sc), + "f-infinity slope - expected: {:.0f}, current: {:.0f}, error: {:.3f}\n".format(self.f_infinities_slope, f_infinities_slope, error_f_inf_slope), + "f-infinity values:\nexpected:", np.around(self.f_infinities), "\ncurrent: ", np.around(f_infinities), "\nerror: {:.3f}".format(error_f_inf)) return error_bf + error_vs + error_sc + error_f_inf_slope + error_f_inf def fit_model_to_data(self, data: CellData): self.calculate_needed_values_from_data(data) - self.counter = 0 return self.fit_model() - def fit_model_to_values(self, eod_freq, baseline_freq, sc, vs, fi_contrasts, fi_inf_values, fi_inf_slope, a_delta, a_tau): + def fit_model_to_values(self, eod_freq, baseline_freq, sc, vs, fi_contrasts, fi_inf_values, fi_inf_slope, a_delta, a_tau, x0=None, init_simplex=None): self.eod_freq = eod_freq self.baseline_freq = baseline_freq self.serial_correlation = sc @@ -210,12 +251,17 @@ class Fitter: self.a_delta = a_delta self.a_tau = a_tau - return self.fit_model() + return self.fit_model(x0, init_simplex) - def fit_model(self): - x0 = np.array([0.02, 3, 75, 0.05, 20]) - init_simplex = np.array([np.array([2, 1, 10]), np.array([40, 10, 140]), np.array([20, 20, 70]), np.array([150, 1, 200])]) - fmin = minimize(fun=self.cost_function, x0=x0, args=(self.a_tau, self.a_delta), method="Nelder-Mead") #, options={"initial_simplex": init_simplex}) + def fit_model(self, x0=None, initial_simplex=None): + self.counter = 0 + if x0 is None: + x0 = np.array([0.02, 0.03, 70]) + if initial_simplex is None: + initial_simplex = create_init_simples(x0) + # fmin = minimize(fun=self.cost_function, x0=x0, args=(self.a_tau, self.a_delta), method="Nelder-Mead") + # else: + fmin = minimize(fun=self.cost_function, x0=x0, args=(self.a_tau, self.a_delta), method="Nelder-Mead", options={"initial_simplex": initial_simplex}) return fmin, self.model.get_parameters() diff --git a/functions.py b/functions.py index 3459db0..d094c20 100644 --- a/functions.py +++ b/functions.py @@ -11,6 +11,17 @@ def inverse_line(x, m, c): return (x-c)/m +def clipped_line(x, m, c): + return np.clip(m*x + c, 0, None) + + +def inverse_clipped_line(x, a, b): + if clipped_line(x, a, b) == 0: + raise ValueError("Value undefined in inverse_clipped_line.") + + return (x-a)/b + + def exponential_function(x, a, b, c): return (a-c)*np.exp(-x/b)+c @@ -38,17 +49,6 @@ def inverse_full_boltzmann(x, f_max, f_min, k, x_zero): return -(np.log((f_max-f_min) / (x - f_min) - 1) / k) + x_zero -def clipped_line(x, m, c): - return np.clip(m*x + c, 0, None) - - -def inverse_clipped_line(x, a, b): - if clipped_line(x, a, b) == 0: - raise ValueError("Value undefined in inverse_clipped_line.") - - return (x-a)/b - - @jit(nopython=True) # useful in less that 1000x10000 calls (1000 tests with 10k data points) def rectify(x): if x < 0: diff --git a/generalTests.py b/generalTests.py index 6139142..0a58126 100644 --- a/generalTests.py +++ b/generalTests.py @@ -195,11 +195,20 @@ def rectify_stimulus_array(stimulus_array: np.ndarray): if __name__ == '__main__': - + # X = [0.05, 0.02, 50, 0.1, 0.03] model = LifacNoiseModel() + # model.set_variable("mem_tau", X[0]) + # model.set_variable("noise_strength", X[1]) + # model.set_variable("input_scaling", X[2]) + # model.set_variable("tau_a", X[3]) + # model.set_variable("delta_a", X[4]) stim = SinusoidalStepStimulus(700, 0.2, start_time=1, duration=1) bf, vs, sc = model.calculate_baseline_markers(700) - print(bf, vs, "\n", sc) + print("baseline freq:{:.2f}\nVector strength: {:.3f}\nSerial cor:".format(bf, vs), sc) + contrasts = np.arange(-0.3, 0.31, 0.05) + model.calculate_fi_curve(contrasts, 700) + f_infinities, slope = model.calculate_fi_markers(contrasts, 700) + print("FI-Curve\nSlope: {:.2f}\nValues:".format(slope), f_infinities) plot_model_during_stimulus(model, stim, 3) quit() diff --git a/helperFunctions.py b/helperFunctions.py index 68f1dcd..d926007 100644 --- a/helperFunctions.py +++ b/helperFunctions.py @@ -398,3 +398,57 @@ def __vector_strength__(relative_spike_times: np.ndarray, eod_durations: np.ndar return vs + +def detect_f_zero_in_frequency_trace(time, frequency, stimulus_start, sampling_interval, peak_buffer_percent=0.05, buffer=0.025): + 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)] + min_before = min(freq_before) + max_before = max(freq_before) + mean_before = np.mean(freq_before) + + # time where the f-zero is searched in + start_idx = int((stimulus_start-0.1*buffer) / sampling_interval) + end_idx = int((stimulus_start + buffer) / sampling_interval) + + min_during_start_of_stim = min(frequency[start_idx:end_idx]) + max_during_start_of_stim = max(frequency[start_idx:end_idx]) + + if abs(mean_before-min_during_start_of_stim) > abs(max_during_start_of_stim-mean_before): + f_zero = min_during_start_of_stim + else: + f_zero = max_during_start_of_stim + + peak_buffer = (max_before - min_before) * peak_buffer_percent + if min_before - peak_buffer <= f_zero <= max_before + peak_buffer: + end_idx = start_idx + int((end_idx-start_idx)/2) + f_zero = np.mean(frequency[start_idx:end_idx]) + + # import matplotlib.pyplot as plt + # plt.plot(time, frequency) + # plt.plot(time[start_idx:end_idx], [f_zero for i in range(end_idx-start_idx)]) + # plt.show() + + return f_zero + + +def detect_f_infinity_in_freq_trace(time, frequency, stimulus_start, stimulus_duration, sampling_interval, length=0.1, buffer=0.025): + stimulus_end_time = stimulus_start + stimulus_duration - time[0] + + start_idx = int((stimulus_end_time - length - buffer) / sampling_interval) + end_idx = int((stimulus_end_time - buffer) / sampling_interval) + + return np.mean(frequency[start_idx:end_idx]) + + +def detect_f_baseline_in_freq_trace(time, frequency, stimulus_start, sampling_interval, buffer=0.025): + stim_start = stimulus_start - time[0] + + if stim_start < 0.1: + warn("FICurve:__calculate_f_baseline__(): Quite short delay at the start.") + + start_idx = int(buffer/sampling_interval) + end_idx = int((stim_start-buffer)/sampling_interval) + f_baseline = np.mean(frequency[start_idx:end_idx]) + + return f_baseline \ No newline at end of file diff --git a/models/LIFACnoise.py b/models/LIFACnoise.py index cedefb0..e0b0c6c 100644 --- a/models/LIFACnoise.py +++ b/models/LIFACnoise.py @@ -21,7 +21,7 @@ class LifacNoiseModel(AbstractModel): "input_scaling": 60, "delta_a": 0.08, "tau_a": 0.1, - "a_zero": 10, + "a_zero": 2, "noise_strength": 0.05, "step_size": 0.00005} @@ -157,14 +157,14 @@ class LifacNoiseModel(AbstractModel): def get_model_copy(self): return LifacNoiseModel(self.parameters) - def calculate_baseline_markers(self, base_stimulus_freq, max_lag=1): + def calculate_baseline_markers(self, 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 = SinusoidalStepStimulus(base_stimulus_freq, 0) + base_stimulus = SinusoidalStepStimulus(stimulus_freq, 0) _, spiketimes = self.simulate_fast(base_stimulus, 30) time_x = 5 baseline_freq = hF.mean_freq_of_spiketimes_after_time_x(spiketimes, time_x) @@ -181,7 +181,7 @@ class LifacNoiseModel(AbstractModel): return baseline_freq, vector_strength, serial_correlation - def calculate_fi_markers(self, contrasts, base_freq, modulation_frequency): + def calculate_fi_markers(self, contrasts, stimulus_freq): """ calculates the fi markers f_infinity, f_infinity_slope for given contrasts based on simulated 2 seconds for each contrast @@ -190,7 +190,7 @@ class LifacNoiseModel(AbstractModel): f_infinities = [] for contrast in contrasts: - stimulus = SinusoidalStepStimulus(base_freq, contrast) + stimulus = SinusoidalStepStimulus(stimulus_freq, contrast) _, spiketimes = self.simulate_fast(stimulus, 1) f_infinity = hF.mean_freq_of_spiketimes_after_time_x(spiketimes, 0.3) @@ -202,7 +202,45 @@ class LifacNoiseModel(AbstractModel): return f_infinities, f_infinities_slope - def find_v_offset(self, goal_baseline_frequency, base_stimulus, threshold=10, border=50000): + def calculate_fi_curve(self, contrasts, stimulus_freq): + stim_duration = 1 + stim_start = 1 + sampling_interval = self.get_sampling_interval() + f_infinities = [] + f_zeros = [] + f_baselines = [] + import matplotlib.pyplot as plt + for c in contrasts: + stimulus = SinusoidalStepStimulus(stimulus_freq, c, stim_start, stim_duration) + _, spiketimes = self.simulate_fast(stimulus, 3.5) + + time, frequency = hF.calculate_time_and_frequency_trace(spiketimes, sampling_interval) + + f_inf = hF.detect_f_infinity_in_freq_trace(time, frequency, stim_start, stim_duration, sampling_interval) + f_infinities.append(f_inf) + + f_zero = hF.detect_f_zero_in_frequency_trace(time, frequency, stim_start, sampling_interval) + f_zeros.append(f_zero) + + f_baseline = hF.detect_f_baseline_in_freq_trace(time, frequency, stim_start, sampling_interval) + f_baselines.append(f_baseline) + + + # fig, axes = plt.subplots(2, 1, sharex="all") + # stim_time = np.arange(0,3.5, sampling_interval) + # axes[0].set_title("Contrast: " + str(c)) + # axes[0].plot(stim_time, [stimulus.value_at_time_in_s(t) for t in stim_time]) # stimulus.as_array(0, 3.5, sampling_interval)) + # + # axes[1].plot(time, frequency) + # axes[1].plot((time[0], time[-1]), (f_inf, f_inf), label="inf") + # axes[1].plot((time[0], time[-1]), (f_zero, f_zero), label="zero") + # axes[1].plot((time[0], time[-1]), (f_baseline, f_baseline), label="base") + # plt.legend() + # plt.show() + + return f_baselines, f_zeros, f_infinities + + def find_v_offset(self, goal_baseline_frequency, base_stimulus, threshold=2, border=50000): test_model = self.get_model_copy() simulation_length = 5 @@ -218,8 +256,8 @@ class LifacNoiseModel(AbstractModel): current_v_offset += v_search_step_size current_freq = test_v_offset(test_model, current_v_offset, base_stimulus, simulation_length) - if current_v_offset == 0: - return -1000 + # if current_v_offset == 0: + # return -1000 lower_bound = current_v_offset - v_search_step_size upper_bound = current_v_offset @@ -235,8 +273,7 @@ def binary_search_base_freq(model: LifacNoiseModel, base_stimulus, goal_frequenc counter += 1 middle = upper_bound - (upper_bound - lower_bound)/2 frequency = test_v_offset(model, middle, base_stimulus, simulation_length) - if counter > 100: - print("meep") + # print('{:.1f}, {:.1f}, {:.1f}, {:.1f} vs {:.1f} '.format(lower_bound, middle, upper_bound, frequency, goal_frequency)) if abs(frequency - goal_frequency) < threshold: return middle @@ -244,7 +281,8 @@ def binary_search_base_freq(model: LifacNoiseModel, base_stimulus, goal_frequenc lower_bound = middle elif frequency > goal_frequency: upper_bound = middle - elif abs(upper_bound-lower_bound) < 0.01 * threshold: + elif abs(upper_bound-lower_bound) < 0.001: + warn("Search was stopped no value was found!") return 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)) diff --git a/stimuli/SinusoidalStepStimulus.py b/stimuli/SinusoidalStepStimulus.py index bd383b9..6271cc7 100644 --- a/stimuli/SinusoidalStepStimulus.py +++ b/stimuli/SinusoidalStepStimulus.py @@ -15,7 +15,7 @@ class SinusoidalStepStimulus(AbstractStimulus): def value_at_time_in_s(self, time_point): carrier = np.sin(2 * np.pi * self.frequency * time_point) - if time_point < self.start_time or time_point > self.start_time + self.duration: + if time_point > self.start_time and time_point < self.start_time + self.duration: return self.amplitude * carrier * self.contrast return self.amplitude * carrier diff --git a/unittests/testSinusoidalStepStimulus.py b/unittests/testSinusoidalStepStimulus.py new file mode 100644 index 0000000..09d819e --- /dev/null +++ b/unittests/testSinusoidalStepStimulus.py @@ -0,0 +1,117 @@ +from stimuli.SinusoidalStepStimulus import SinusoidalStepStimulus +import unittest +import numpy as np +import helperFunctions as hF +import matplotlib.pyplot as plt +from warnings import warn + + +class SinusoidalStepStimulusTester(unittest.TestCase): + + base_frequencies = [0, 10, 100, 1000] + contrasts = [0, 0.5, 1, 1.5] + modulation_frequencies = [0, 5, 10, 100] + step_sizes = [1, 0.5, 0.00005] + time_starts = [0, 2, -2] + durations = [2] + + def setUp(self): + pass + + def tearDown(self): + pass + + def test_consistency_base_frequency(self): + contrast = 0.1 + mod_freq = 5 + time_start = -1 + duration = 10 + step_size = 0.00005 + for base_freq in self.base_frequencies: + stimulus = SinusoidalStepStimulus(base_freq, contrast, 0, 8) + self.assertTrue(array_and_time_points_equal(stimulus, time_start, duration, step_size), + msg="Stimulus values inconsistent with base freq: {:.2f}".format(base_freq)) + + def test_consistency_contrast(self): + base_freq = 700 + mod_freq = 5 + time_start = -1 + duration = 10 + step_size = 0.00005 + for contrast in self.contrasts: + stimulus = SinusoidalStepStimulus(base_freq, contrast, 0, 8) + self.assertTrue(array_and_time_points_equal(stimulus, time_start, duration, step_size), + msg="Stimulus values inconsistent with contrast: {:.2f}".format(contrast)) + + def test_consistency_modulation_frequency(self): + contrast = 0.1 + base_freq = 700 + time_start = -1 + duration = 10 + step_size = 0.00005 + for mod_freq in self.modulation_frequencies: + stimulus = SinusoidalStepStimulus(base_freq, contrast, 0, 1) + self.assertTrue(array_and_time_points_equal(stimulus, time_start, duration, step_size), + msg="Stimulus values inconsistent with mod freq: {:.2f}".format(mod_freq)) + + def test_consistency_step_size(self): + contrast = 0.1 + base_freq = 700 + time_start = -1 + duration = 10 + mod_freq = 10 + for step_size in self.step_sizes: + stimulus = SinusoidalStepStimulus(base_freq, contrast, 0, 8) + self.assertTrue(array_and_time_points_equal(stimulus, time_start, duration, step_size), + msg="Stimulus values inconsistent with step_size: {:.3f}ms".format(step_size)*1000) + + def test_consistency_time_start(self): + contrast = 0.1 + base_freq = 700 + mod_freq = 10 + duration = 10 + step_size = 0.00005 + for time_start in self.time_starts: + stimulus = SinusoidalStepStimulus(base_freq, contrast, 0, 8) + self.assertTrue(array_and_time_points_equal(stimulus, time_start, duration, step_size), + msg="Stimulus values inconsistent when the time starts at: {:.2f}s".format(time_start)) + + +def array_and_time_points_equal(stimulus, start, duration, step_size): + precision = 5 + array = np.around(stimulus.as_array(start, duration, step_size), precision) + time = np.arange(start, start+duration, step_size) + for i, time_point in enumerate(time): + value = stimulus.value_at_time_in_s(time_point) + if array[i] != np.round(value, precision): + stim_per_point = [] + for t in time: + stim_per_point.append(stimulus.value_at_time_in_s(t)) + + stim_per_point = np.around(np.array(stim_per_point), precision) + fig, axes = plt.subplots(2, 1, sharex="all") + axes[0].plot(time, array, label="array") + axes[0].plot(time, stim_per_point, label="individual") + axes[0].set_title("stimulus values") + axes[0].legend() + + axes[1].plot(time, np.array(stim_per_point)-array) + axes[1].set_title("difference") + plt.show() + + return False + + # stim_per_point = [] + # for t in time: + # stim_per_point.append(stimulus.value_at_time_in_s(t)) + # + # stim_per_point = np.around(np.array(stim_per_point), precision) + # fig, axes = plt.subplots(1, 1, sharex="all") + # axes.plot(time, array, label="array") + # axes.plot(time, stim_per_point, label="individual") + # axes.set_title("stimulus values") + # axes.legend() + # + # plt.show() + + return True