add behaviour distributions beginning of sensitivity analysis

analysis.py

@@ -6,6 +6,10 @@ import matplotlib.pyplot as plt
 from scipy.stats import pearsonr
 
 from ModelFit import get_best_fit
+from Baseline import BaselineModel
+from FiCurve import FICurveModel
+from CellData import CellData
+from stimuli.SinusoidalStepStimulus import SinusoidalStepStimulus
 
 
 def main():
@@ -19,19 +23,20 @@ def main():
     #     print("Argument dir is not a directory.")
     #     parser.print_usage()
     #     exit(0)
-
+    # sensitivity_analysis(dir_path, max_models=2)
     fits_info = get_fit_info(dir_path)
 
-    # errors = calculate_percent_errors(fits_info)
+    errors = calculate_percent_errors(fits_info)
-    # create_boxplots(errors)
+    create_boxplots(errors)
-    # labels, corr_values, corrected_p_values = behaviour_correlations(fits_info, model_values=False)
+    labels, corr_values, corrected_p_values = behaviour_correlations(fits_info, model_values=False)
-    # create_correlation_plot(labels, corr_values, corrected_p_values)
+    create_correlation_plot(labels, corr_values, corrected_p_values)
 
-    # labels, corr_values, corrected_p_values = parameter_correlations(fits_info)
+    labels, corr_values, corrected_p_values = parameter_correlations(fits_info)
-    # create_correlation_plot(labels, corr_values, corrected_p_values)
+    create_correlation_plot(labels, corr_values, corrected_p_values)
 
     create_parameter_distributions(get_parameter_values(fits_info))
-
+    cell_b, model_b = get_behaviour_values(fits_info)
+    create_behaviour_distributions(cell_b, model_b)
     pass
 
 
@@ -183,8 +188,6 @@ def create_boxplots(errors):
     plt.show()
     plt.close()
 
-    pass
-
 
 def create_parameter_distributions(par_values):
 
@@ -209,6 +212,98 @@ def create_parameter_distributions(par_values):
 
 
 def create_behaviour_distributions(cell_b_values, model_b_values):
+    fig, axes = plt.subplots(4, 2)
+
+    labels = sorted(cell_b_values.keys())
+    axes_flat = axes.flatten()
+    for i, l in enumerate(labels):
+        min_v = min(min(cell_b_values[l]), min(model_b_values[l])) * 0.95
+        max_v = max(max(cell_b_values[l]), max(model_b_values[l])) * 1.05
+        limit = 50000
+        if max_v > limit:
+            print("For {} the max value was limited to {},  {} values were excluded!".format(l, limit, np.sum(np.array(cell_b_values[l]) > limit)))
+            max_v = limit
+        step = (max_v - min_v) / 15
+        bins = np.arange(min_v, max_v + step, step)
+        axes_flat[i].hist(cell_b_values[l], bins=bins, alpha=0.5)
+        axes_flat[i].hist(model_b_values[l], bins=bins, alpha=0.5)
+        axes_flat[i].set_title(l)
+
+    plt.tight_layout()
+    plt.show()
+    plt.close()
+
+    pass
+
+
+def sensitivity_analysis(dir_path, par_range=(0.5, 1.6, 0.1), contrast_range=(-0.3, 0.4, 0.1), parameters=None, behaviours=None, max_models=None):
+    models = []
+    eods = []
+    base_freqs = []
+
+    count = 0
+    for item in sorted(os.listdir(dir_path)):
+        count += 1
+        if max_models is not None and count > max_models:
+            break
+        cell_folder = os.path.join(dir_path, item)
+
+        results = get_best_fit(cell_folder)
+
+        models.append(results.get_model())
+        eods.append(CellData(results.get_cell_path()).get_eod_frequency())
+        cell, model = results.get_behaviour_values()
+        base_freqs.append(cell["baseline_frequency"])
+
+    if parameters is None:
+        parameters = ["input_scaling", "delta_a", "mem_tau", "noise_strength",
+                      "refractory_period", "tau_a", "dend_tau"]
+
+    if behaviours is None:
+        behaviours = ["burstiness", "coefficient_of_variation", "serial_correlation",
+                      "vector_strength", "f_inf_slope", "f_zero_slope", "f_zero_middle"]
+
+    model_behaviour_responses = []
+
+    contrasts = np.arange(contrast_range[0], contrast_range[1], contrast_range[2])
+    factors = np.arange(par_range[0], par_range[1], par_range[2])
+    for model, eod, base_freq in zip(models, eods, base_freqs):
+        par_responses = {}
+        for par in parameters:
+            par_responses[par] = {}
+            for b in behaviours:
+                par_responses[par][b] = np.zeros(len(factors))
+            for i, factor in enumerate(factors):
+
+                model_copy = model.get_model_copy()
+                model_copy.set_variable(par, model.get_parameters()[par] * factor)
+                print("{} at {}, ({} of {})".format(par, model.get_parameters()[par] * factor, i+1, len(factors)))
+                base_stimulus = SinusoidalStepStimulus(eod, 0)
+                v_offset = model_copy.find_v_offset(base_freq, base_stimulus)
+                model_copy.set_variable("v_offset", v_offset)
+
+                baseline = BaselineModel(model_copy, eod, trials=3)
+                print(baseline.get_baseline_frequency())
+                if "burstiness" in behaviours:
+                    par_responses[par]["burstiness"][i] = baseline.get_burstiness()
+                if "coefficient_of_variation" in behaviours:
+                    par_responses[par]["coefficient_of_variation"][i] = baseline.get_coefficient_of_variation()
+                if "serial_correlation" in behaviours:
+                    par_responses[par]["serial_correlation"][i] = baseline.get_serial_correlation(1)[0]
+                if "vector_strength" in behaviours:
+                    par_responses[par]["vector_strength"][i] = baseline.get_vector_strength()
+
+                fi_curve = FICurveModel(model_copy, contrasts, eod, trials=10)
+
+                if "f_inf_slope" in behaviours:
+                    par_responses[par]["f_inf_slope"][i] = fi_curve.get_f_inf_slope()
+                if "f_zero_slope" in behaviours:
+                    par_responses[par]["f_zero_slope"][i] = fi_curve.get_f_zero_fit_slope_at_straight()
+                if "f_zero_middle" in behaviours:
+                    par_responses[par]["f_zero_middle"][i] = fi_curve.f_zero_fit[3]
+
+        model_behaviour_responses.append(par_responses)
+
     pass