diff --git a/analysis.py b/analysis.py index eb50555..bbf695e 100644 --- a/analysis.py +++ b/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) - # create_boxplots(errors) - # labels, corr_values, corrected_p_values = behaviour_correlations(fits_info, model_values=False) - # create_correlation_plot(labels, corr_values, corrected_p_values) + errors = calculate_percent_errors(fits_info) + create_boxplots(errors) + labels, corr_values, corrected_p_values = behaviour_correlations(fits_info, model_values=False) + create_correlation_plot(labels, corr_values, corrected_p_values) - # labels, corr_values, corrected_p_values = parameter_correlations(fits_info) - # create_correlation_plot(labels, corr_values, corrected_p_values) + labels, corr_values, corrected_p_values = parameter_correlations(fits_info) + 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