diff --git a/analysis.py b/analysis.py new file mode 100644 index 0000000..eb50555 --- /dev/null +++ b/analysis.py @@ -0,0 +1,216 @@ + +import argparse +import os +import numpy as np +import matplotlib.pyplot as plt +from scipy.stats import pearsonr + +from ModelFit import get_best_fit + + +def main(): + # parser = argparse.ArgumentParser() + # parser.add_argument("dir", help="folder containing the cell folders with the fit results") + # args = parser.parse_args() + + dir_path = "results/invivo_results/" # args.dir + + # if not os.path.isdir(dir_path): + # print("Argument dir is not a directory.") + # parser.print_usage() + # exit(0) + + 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) + + # 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)) + + pass + + +def calculate_percent_errors(fits_info): + errors = {} + + for cell in sorted(fits_info.keys()): + for behaviour in fits_info[cell][1].keys(): + if behaviour not in errors.keys(): + errors[behaviour] = [] + + if fits_info[cell][2][behaviour] == 0: + if fits_info[cell][1][behaviour] == 0: + errors[behaviour].append(0) + else: + print("Cannot calc % error if reference is 0") + continue + errors[behaviour].append((fits_info[cell][1][behaviour] - fits_info[cell][2][behaviour]) / fits_info[cell][2][behaviour]) + return errors + + +def get_parameter_values(fits_info): + par_keys = sorted(["input_scaling", "delta_a", "mem_tau", "noise_strength", + "refractory_period", "tau_a", "v_offset", "dend_tau"]) + parameter_values = {} + for cell in sorted(fits_info.keys()): + for par in par_keys: + if par not in parameter_values.keys(): + parameter_values[par] = [] + + parameter_values[par].append(fits_info[cell][0][par]) + + return parameter_values + + +def get_behaviour_values(fits_info): + behaviour_values_cell = {} + behaviour_values_model = {} + for cell in sorted(fits_info.keys()): + for behaviour in fits_info[cell][1].keys(): + if behaviour not in behaviour_values_cell.keys(): + behaviour_values_cell[behaviour] = [] + behaviour_values_model[behaviour] = [] + + behaviour_values_model[behaviour].append(fits_info[cell][1][behaviour]) + behaviour_values_cell[behaviour].append(fits_info[cell][2][behaviour]) + + return behaviour_values_cell, behaviour_values_model + + +def behaviour_correlations(fits_info, model_values=True): + bv_cell, bv_model = get_behaviour_values(fits_info) + if model_values: + behaviour_values = bv_model + else: + behaviour_values = bv_cell + + labels = sorted(behaviour_values.keys()) + corr_values = np.zeros((len(labels), len(labels))) + p_values = np.ones((len(labels), len(labels))) + + for i in range(len(labels)): + for j in range(len(labels)): + c, p = pearsonr(behaviour_values[labels[i]], behaviour_values[labels[j]]) + corr_values[i, j] = c + p_values[i, j] = p + + corrected_p_values = p_values * sum(range(len(labels))) + + return labels, corr_values, corrected_p_values + + +def parameter_correlations(fits_info): + parameter_values = get_parameter_values(fits_info) + + labels = sorted(parameter_values.keys()) + corr_values = np.zeros((len(labels), len(labels))) + p_values = np.ones((len(labels), len(labels))) + + for i in range(len(labels)): + for j in range(len(labels)): + c, p = pearsonr(parameter_values[labels[i]], parameter_values[labels[j]]) + corr_values[i, j] = c + p_values[i, j] = p + + corrected_p_values = p_values * sum(range(len(labels))) + + return labels, corr_values, corrected_p_values + + +def get_fit_info(folder): + fits_info = {} + + for item in os.listdir(folder): + cell_folder = os.path.join(folder, item) + + results = get_best_fit(cell_folder) + cell_behaviour, model_behaviour = results.get_behaviour_values() + fits_info[item] = [results.get_final_parameters(), model_behaviour, cell_behaviour] + + return fits_info + + +def create_correlation_plot(labels, correlations, p_values): + + cleaned_cors = np.zeros(correlations.shape) + + for i in range(correlations.shape[0]): + for j in range(correlations.shape[1]): + if abs(p_values[i, j]) < 0.05: + cleaned_cors[i, j] = correlations[i, j] + + fig, ax = plt.subplots() + im = ax.imshow(cleaned_cors, vmin=-1, vmax=1) + + cbar = ax.figure.colorbar(im, ax=ax) + cbar.ax.set_ylabel("Correlation coefficient", rotation=-90, va="bottom") + + # We want to show all ticks... + ax.set_xticks(np.arange(len(labels))) + ax.set_yticks(np.arange(len(labels))) + # ... and label them with the respective list entries + ax.set_xticklabels(labels) + ax.set_yticklabels(labels) + + # Rotate the tick labels and set their alignment. + plt.setp(ax.get_xticklabels(), rotation=45, ha="right", + rotation_mode="anchor") + + # Loop over data dimensions and create text annotations. + for i in range(len(labels)): + for j in range(len(labels)): + text = ax.text(j, i, "{:.2f}".format(correlations[i, j]), + ha="center", va="center", color="w") + + fig.tight_layout() + plt.show() + + +def create_boxplots(errors): + + labels = ["{}_n:{}".format(k, len(errors[k])) for k in sorted(errors.keys())] + + y_values = [errors[k] for k in sorted(errors.keys())] + + plt.boxplot(y_values) + plt.xticks(np.arange(1, len(y_values)+1, 1), labels, rotation=45) + plt.tight_layout() + plt.show() + plt.close() + + pass + + +def create_parameter_distributions(par_values): + + fig, axes = plt.subplots(4, 2) + + if len(par_values.keys()) != 8: + print("not eight parameters") + + labels = sorted(par_values.keys()) + axes_flat = axes.flatten() + for i, l in enumerate(labels): + min_v = min(par_values[l]) * 0.95 + max_v = max(par_values[l]) * 1.05 + step = (max_v - min_v) / 15 + bins = np.arange(min_v, max_v+step, step) + axes_flat[i].hist(par_values[l], bins=bins) + axes_flat[i].set_title(l) + + plt.tight_layout() + plt.show() + plt.close() + + +def create_behaviour_distributions(cell_b_values, model_b_values): + pass + + +if __name__ == '__main__': + main() diff --git a/fit_statistics.py b/fit_statistics.py deleted file mode 100644 index ea825ba..0000000 --- a/fit_statistics.py +++ /dev/null @@ -1,64 +0,0 @@ - -from ModelFit import get_best_fit -import matplotlib.pyplot as plt - -import os -import numpy as np - - -def get_all_cell_folders(dir): - paths = [] - for item in os.listdir(dir): - cell_path = os.path.join(dir, item) - if os.path.isdir(cell_path): - paths.append(cell_path) - - return paths - - -def main(): - - values_models = {} - values_cells = {} - - for cell in get_all_cell_folders("results/invivo_results/"): - best_fit = get_best_fit(cell) - c_behaviour, m_behaviour = best_fit.get_behaviour_values() - - for key in c_behaviour.keys(): - if key not in values_cells.keys(): - values_cells[key] = [] - if key not in values_models.keys(): - values_models[key] = [] - - values_cells[key].append(c_behaviour[key]) - values_models[key].append(m_behaviour[key]) - - percentage_error_data = [] - labels = [] - for key in values_models.keys(): - errors = [] - for i in range(len(values_models[key])): - if values_cells[key][i] == 0: - if values_models[key][i] == 0: - errors.append(0) - else: - print("Cannot calc % error if reference is 0") - continue - errors.append((values_models[key][i] / values_cells[key][i]) -1) - - percentage_error_data.append(errors) - labels.append(key) - - plt.boxplot(percentage_error_data) - plt.xticks(range(1, len(percentage_error_data)+1, 1), labels=labels, rotation=45) - plt.show() - plt.close() - - # plt.plot([np.median(x) for x in percentage_error_data], 'o') - # plt.xticks(range(len(percentage_error_data)), labels=labels, rotation=45) - # plt.show() - - -if __name__ == '__main__': - main()