add tests for sinusoidalStep, add FiCurve plot possibilities (temp), add initial simplex generation

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]")

fit_lifacnoise.py

@@ -12,8 +12,8 @@ import matplotlib.pyplot as plt
 
 
 def main():
-    # run_test_with_fixed_model()
+    #run_test_with_fixed_model()
-    # quit()
+    #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_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_base_frequency(), 10)
+        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,
+    parameters = {"mem_tau": 0.015,
-                  'v_offset': 80, 'threshold': 1, 'v_base': 0, 'step_size': 0.00005, 'tau_a': a_tau,
+                      "v_base": 0,
-                  'a_zero': 0, 'v_zero': 0, 'noise_strength': 0.5}
+                      "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)
+    starting_conditions = [[0.05, 0.02, 50], [0.01, 0.08, 75], [0.02, 0.03, 70],
-    print("calculated parameters:")
+                           [0.02, 0.03, 70]]
-    print(fit_parameters)
+    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("tau_a", X[3])
-        self.model.set_variable("delta_a", X[4])
+        #self.model.set_variable("delta_a", X[4])
-        # self.model.set_variable("tau_a", tau_a)
+        self.model.set_variable("tau_a", tau_a)
-        # self.model.set_variable("delta_a", delta_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)
-
+        test_model = self.model.get_model_copy()
-        v_offset = self.model.find_v_offset(self.baseline_freq, base_stimulus)
+        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:
+        if self.counter % 100 == 0:
-            print("Cost function run times:", self.counter, "error sum:", sum(errors), errors)
+            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):
+    def fit_model(self, x0=None, initial_simplex=None):
-        x0 = np.array([0.02, 3, 75, 0.05, 20])
+        self.counter = 0
-        init_simplex = np.array([np.array([2, 1, 10]), np.array([40, 10, 140]), np.array([20, 20, 70]), np.array([150, 1, 200])])
+        if x0 is None:
-        fmin = minimize(fun=self.cost_function, x0=x0, args=(self.a_tau, self.a_delta), method="Nelder-Mead")  #, options={"initial_simplex": init_simplex})
+            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()
 

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:

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()

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

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:
+        # if current_v_offset == 0:
-            return -1000
+        #     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))

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

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