217 lines
6.5 KiB
Python
217 lines
6.5 KiB
Python
|
|
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()
|