add function to test the effect of eod frequency and refractory period

variableEffect.py

@@ -8,12 +8,14 @@ import copy
 import os
 
 
-SEARCH_WIDTH = 1.1
+SEARCH_WIDTH = 3
-SEARCH_PRECISION = 1
+SEARCH_PRECISION = 40
 CONTRASTS = np.arange(-0.4, 0.45, 0.1)
 
 
 def main():
+    test_effect_of_two_variables()
+    quit()
     model_parameters1 = {'threshold': 1,
                          'step_size': 5e-05,
                          'a_zero': 2,
@@ -25,23 +27,110 @@ def main():
                          'tau_a': 0.05958195972016753,
                          'input_scaling': 119.81500448274554,
                          'dend_tau': 0.0027746086464721723,
-                         'v_offset': -24.21875}
+                         'v_offset': -24.21875,
+                         'refractory_period': 0.0006}
 
     model_parameters2 = {'v_offset': -15.234375, 'input_scaling': 64.94152780134829, 'step_size': 5e-05, 'a_zero': 2,
                          'threshold': 1, 'v_base': 0, 'delta_a': 0.04763179657857666, 'tau_a': 0.07891848949732623,
                          'mem_tau': 0.004828473985707999, 'noise_strength': 0.017132801387559883,
-                         'v_zero': 0, 'dend_tau': 0.0015230454266819539}
+                         'v_zero': 0, 'dend_tau': 0.0015230454266819539, 'refractory_period': 0.0006}
 
-    parameters_to_test = ["input_scaling", "dend_tau", "mem_tau", "noise_strength", "v_offset", "delta_a", "tau_a"]
+    parameters_to_test = ["input_scaling", "refractory_period", "dend_tau", "mem_tau", "noise_strength", "v_offset", "delta_a", "tau_a"]
+    # parameters_to_test = ["refractory_period", "input_scaling"]
     effect_data = []
     for p in parameters_to_test:
         print("Working on parameter " + p)
-        effect_data.append(test_parameter_effect(model_parameters2, p))
+        effect_data.append(test_parameter_effect(model_parameters2, p, 600))
 
     plot_effects(effect_data, "./figures/variable_effect/")
 
 
-def test_parameter_effect(model_parameters, test_parameter):
+def test_effect_of_two_variables():
+    eod_freqs = np.arange(100, 1001, 20)
+    ref_periods = np.arange(0, 0.00201, 0.0001)
+    variables = ("bf", "vs", "sc", "cv", "burst", "f_inf_s", "f_zero_s")
+    colorbar_labels = ("Frequency in Hz", "Vector strength", "serial correlation lag=1", "Coefficient of Variation",
+                       "Burstiness", "f_inf slope", "f_zero slope")
+
+    # base eod frequency would be 771!
+    base_parameters = {'step_size': 5e-05, 'mem_tau': 0.0076958612706114595, 'v_base': 0, 'v_zero': 0, 'threshold': 1,
+                       'v_offset': -37.5, 'input_scaling': 181.40702315746051, 'delta_a': 0.333391796423963,
+                       'tau_a': 0.17301586167067445, 'a_zero': 2, 'noise_strength': 0.017424670423939775,
+                       'dend_tau': 0.0037179224836952356, 'refractory_period': 0.0010602702699897444}
+
+    # base eod frequency would be 657!
+    base_parameters = {'refractory_period': 0.0008347981797599925, 'v_base': 0, 'v_zero': 0, 'a_zero': 2,
+                       'step_size': 5e-05, 'delta_a': 0.10570085698152036, 'threshold': 1,
+                       'input_scaling': 85.7818875779873, 'mem_tau': 0.01094261953657057, 'tau_a': 0.07741757133763925,
+                       'v_offset': -10.15625, 'noise_strength': 0.03080655781041302, 'dend_tau': 0.0013624430015225777}
+
+
+    effects = []
+    for eod_freq in eod_freqs:
+        effects_with_const_eod_freq = []
+        for ref_period in ref_periods:
+            model_parameters = copy.deepcopy(base_parameters)
+            model_parameters["refractory_period"] = ref_period
+            effects_with_const_eod_freq.append(test_model(model_parameters, eod_freq))
+
+        effects.append(effects_with_const_eod_freq)
+
+    if not os.path.exists("./figures/eod_and_ref_period_effect/"):
+        os.makedirs("./figures/eod_and_ref_period_effect/")
+
+    for x, variable in enumerate(variables):
+        matrix = np.zeros((len(eod_freqs), len(ref_periods)))
+        for i in range(len(eod_freqs)):
+            for j in range(len(ref_periods)):
+                matrix[i, j] = effects[i][j][variable]
+
+        fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(6, 6))
+
+        im = axes.imshow(matrix)
+        cbar = axes.figure.colorbar(im, ax=axes)
+        cbar.ax.set_ylabel(colorbar_labels[x], rotation=-90, va="bottom")
+
+        axes.set_title(variable)
+        axes.set_xlabel("Refractory periods in ms")
+        axes.set_ylabel("EOD frequency in Hz")
+
+        axes.set_xticks(np.arange(len(ref_periods)))
+        axes.set_yticks(np.arange(len(eod_freqs)))
+        # ... and label them with the respective list entries
+        axes.set_xticklabels(["{:.2f}".format(r*1000) for r in ref_periods])
+        axes.set_yticklabels(eod_freqs)
+        plt.setp(axes.get_xticklabels(), rotation=45, ha="right",
+                 rotation_mode="anchor")
+
+        axes.set_ylabel("Eod frequencies")
+        plt.tight_layout()
+
+        plt.savefig("./figures/eod_and_ref_period_effect/" + variable + ".png")
+        plt.close()
+
+
+def test_model(model_parameters, eod_freq):
+    model = LifacNoiseModel(model_parameters)
+    print(model.get_parameters())
+
+    fi_curve = FICurveModel(model, CONTRASTS, eod_freq, trials=10)
+    f_inf_s = fi_curve.get_f_inf_slope()
+    f_inf_v = fi_curve.get_f_inf_frequencies()
+    f_zero_s = fi_curve.get_f_zero_fit_slope_at_stimulus_value(0.1)
+    f_zero_v = fi_curve.get_f_zero_frequencies()
+
+    baseline = BaselineModel(model, eod_freq, trials=3)
+    bf = baseline.get_baseline_frequency()
+    vs = baseline.get_vector_strength()
+    sc = baseline.get_serial_correlation(1)[0]
+    cv = baseline.get_coefficient_of_variation()
+    burst = baseline.get_burstiness()
+
+    return {"f_inf_s": f_inf_s, "f_inf_v": f_inf_v, "f_zero_s": f_zero_s, "f_zero_v": f_zero_v,
+            "bf": bf, "vs": vs, "sc": sc, "cv": cv, "burst": burst}
+
+
+def test_parameter_effect(model_parameters, test_parameter, eod_freq):
     model_parameters = copy.deepcopy(model_parameters)
     start_value = model_parameters[test_parameter]
 
@@ -54,20 +143,19 @@ def test_parameter_effect(model_parameters, test_parameter):
     vs = []
     sc = []
     cv = []
+    burst = []
 
     f_inf_s = []
     f_inf_v = []
     f_zero_s = []
     f_zero_v = []
-    fi_curves = []
     broken_i = []
 
     for i in range(len(values)):
         model_parameters[test_parameter] = values[i]
         model = LifacNoiseModel(model_parameters)
 
-        fi_curve = FICurveModel(model, CONTRASTS, 600, trials=1)
+        fi_curve = FICurveModel(model, CONTRASTS, eod_freq, trials=10)
-        fi_curves.append(fi_curve)
         f_inf_s.append(fi_curve.get_f_inf_slope())
         f_inf_v.append(fi_curve.get_f_inf_frequencies())
         f_zero_s.append(fi_curve.get_f_zero_fit_slope_at_stimulus_value(0.1))
@@ -82,11 +170,12 @@ def test_parameter_effect(model_parameters, test_parameter):
 
         fi_curve.plot_f_point_detections(detection_save_path)
 
-        baseline = BaselineModel(model, 600, trials=1)
+        baseline = BaselineModel(model, eod_freq, trials=3)
         bf.append(baseline.get_baseline_frequency())
         vs.append(baseline.get_vector_strength())
         sc.append(baseline.get_serial_correlation(2))
         cv.append(baseline.get_coefficient_of_variation())
+        burst.append(baseline.get_burstiness())
 
     values = list(values)
     if len(broken_i) > 0:
@@ -94,22 +183,22 @@ def test_parameter_effect(model_parameters, test_parameter):
         for i in broken_i:
             del values[i]
 
-    return ParameterEffectData(fi_curves, values, test_parameter, bf, vs, sc, cv, f_inf_s, f_inf_v, f_zero_s, f_zero_v)
+    return ParameterEffectData(values, test_parameter, bf, vs, sc, cv, burst, f_inf_s, f_inf_v, f_zero_s, f_zero_v)
     # plot_effects(values, test_parameter, bf, vs, sc, cv, f_inf_s, f_inf_v, f_zero_s, f_zero_v)
 
 
 def plot_effects(par_effect_data_list, save_path=None):
 
-    names = ("bf", "vs", "sc", "cv", "f_inf_s", "f_inf_v", "f_zero_s", "f_zero_v", "f_zero_fit_x_0")
+    names = ("bf", "vs", "sc", "cv", "burstiness", "f_inf_s", "f_inf_v", "f_zero_s", "f_zero_v")
 
-    fig, axes = plt.subplots(len(names), len(par_effect_data_list), figsize=(32, 4*len(par_effect_data_list)), sharex="col")
+    fig, axes = plt.subplots(len(names), len(par_effect_data_list), figsize=(4*len(par_effect_data_list), 4*len(names)), sharex="col")
 
     for j in range(len(par_effect_data_list)):
         ped = par_effect_data_list[j]
 
-        ranges = ((0, max(ped.get_data("bf")) * 1.1), (0, 1), (-1, 1), (0, 1),
+        ranges = ((0, max(ped.get_data("bf")) * 1.1), (0, 1), (-1, 1), (0, 1), (0, 1),
                   (0, max(ped.get_data("f_inf_s")) * 1.1), (0, 800),
-                  (0, max(ped.get_data("f_zero_s")) * 1.1), (0, 10000), (-0.5, max(ped.get_data("f_zero_fit_x_0"))))
+                  (0, max(ped.get_data("f_zero_s")) * 1.1), (0, 10000))
         values = ped.values
 
         for i in range(len(names)):
@@ -137,10 +226,9 @@ def plot_effects(par_effect_data_list, save_path=None):
 
 
 class ParameterEffectData:
-    data_names = ("bf", "vs", "sc", "cv", "f_inf_s", "f_inf_v" "f_zero_s", "f_zero_v", "f_zero_fit_x_0")
+    data_names = ("bf", "vs", "sc", "cv", "burstiness","f_inf_s", "f_inf_v" "f_zero_s", "f_zero_v")
 
-    def __init__(self, fi_curves, values, test_parameter, bf, vs, sc, cv, f_inf_s, f_inf_v, f_zero_s, f_zero_v):
+    def __init__(self, values, test_parameter, bf, vs, sc, cv, burstiness, f_inf_s, f_inf_v, f_zero_s, f_zero_v):
-        self.fi_curves = fi_curves
         self.values = values
         self.test_parameter = test_parameter
         self.bf = bf
@@ -151,6 +239,7 @@ class ParameterEffectData:
         self.f_inf_v = f_inf_v
         self.f_zero_s = f_zero_s
         self.f_zero_v = f_zero_v
+        self.burstiness = burstiness
 
     def get_data(self, name):
         if name == "bf":
@@ -169,10 +258,8 @@ class ParameterEffectData:
             return self.f_zero_s
         elif name == "f_zero_v":
             return self.f_zero_v
-        elif name == "f_zero_fit_x_0":
+        elif name == "burstiness":
-            fits = [fi.f_zero_fit for fi in self.fi_curves]
+            return self.burstiness
-            x_zeros = [fit[3] for fit in fits]
-            return x_zeros
         else:
             raise ValueError("Unknown attribute name!")