From 804f94ad0b2ad8db551728dad859b715e62a60ec Mon Sep 17 00:00:00 2001 From: "a.ott" Date: Sun, 3 Jan 2021 13:11:05 +0100 Subject: [PATCH] stuff --- CellData.py | 3 + DataParserFactory.py | 6 ++ lines_of_code.py | 33 +++++++++ sam_experiments.py | 161 ++++++++++++++++++++++++++----------------- 4 files changed, 139 insertions(+), 64 deletions(-) create mode 100644 lines_of_code.py diff --git a/CellData.py b/CellData.py index 9bb6213..d2c3cb9 100644 --- a/CellData.py +++ b/CellData.py @@ -94,6 +94,9 @@ class CellData: def get_cell_name(self): return os.path.basename(self.data_path) + def has_sam_recordings(self): + return self.parser.has_sam_recordings() + def get_baseline_length(self): return self.parser.get_baseline_length() diff --git a/DataParserFactory.py b/DataParserFactory.py index 36a1f80..952540c 100644 --- a/DataParserFactory.py +++ b/DataParserFactory.py @@ -19,6 +19,9 @@ class AbstractParser: def get_baseline_length(self): raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS") + def has_sam_recordings(self): + raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS") + def get_fi_curve_contrasts(self): raise NotImplementedError("NOT YET OVERRIDDEN FROM ABSTRACT CLASS") @@ -70,6 +73,9 @@ class DatParser(AbstractParser): self.fi_recording_times = [] self.sampling_interval = -1 + def has_sam_recordings(self): + return exists(self.sam_file) + def get_baseline_length(self): lengths = [] for metadata, key, data in Dl.iload(self.baseline_file): diff --git a/lines_of_code.py b/lines_of_code.py new file mode 100644 index 0000000..d30a4f6 --- /dev/null +++ b/lines_of_code.py @@ -0,0 +1,33 @@ + +import os + + +def count_lines_folder(folder): + lines_of_code = 0 + files = 0 + for file in os.listdir(folder): + + if os.path.isdir(file): + continue + if not file.endswith(".py"): + continue + # print(file) + files += 1 + with open(os.path.join(folder, file)) as file: + lines_of_code += len(file.readlines()) + return lines_of_code, files + + +total_lines = 0 +total_files = 0 + +folders = [".", "tests/", "models/", "introduction/", "stimuli/"] + +for folder in folders: + lines, files = count_lines_folder(folder) + print(folder, files, lines) + total_lines += lines + total_files += files + +print("Total lines of code:", total_lines) +print("Total files with code:", total_files) \ No newline at end of file diff --git a/sam_experiments.py b/sam_experiments.py index e6a5e68..1e6561d 100644 --- a/sam_experiments.py +++ b/sam_experiments.py @@ -9,10 +9,13 @@ import helperFunctions as hF from CellData import CellData from ModelFit import ModelFit, get_best_fit import os +import shutil def main(): - sam_analysis("results/final_2/2011-10-25-ad-invivo-1/") + run_sam_analysis_for_all_cells("results/final_2") + + # sam_analysis("results/final_2/2011-10-25-ad-invivo-1/") # plot_traces_with_spiketimes() # plot_mean_of_cuts() @@ -27,6 +30,21 @@ def main(): test_model_response(model, eod_freq, 0.1, np.arange(5, 2500, 5)) +def run_sam_analysis_for_all_cells(folder): + count = 0 + for item in os.listdir(folder): + cell_folder = os.path.join(folder, item) + fit = get_best_fit(cell_folder, use_comparable_error=False) + cell_data = fit.get_cell_data() + + if cell_data.has_sam_recordings(): + count += 1 + # print("Fit quality:", fit.get_fit_routine_error()) + sam_analysis(cell_folder) + print(count) + + + def test_model_response(model: LifacNoiseModel, eod_freq, contrast, modulation_frequencies): stds = [] @@ -182,11 +200,11 @@ def sam_analysis(fit_path): # TODO problem of cutting the pdf as in some cases the pdf is shorter than 1 modulation frequency period! # length info wrong ? always at least one period? - if 1/mod_freq > durations[0] / 4: - print("skipped mod_freq: {}".format(mod_freq)) - print("Duration: {} while mod_freq period: {:.2f}".format(durations[0], 1/mod_freq)) - print("Maybe long enough duration? unique durations:", u_durations) - continue + # if 1/mod_freq > durations[0] / 4: + # print("skipped mod_freq: {}".format(mod_freq)) + # print("Duration: {} while mod_freq period: {:.2f}".format(durations[0], 1/mod_freq)) + # print("Maybe long enough duration? unique durations:", u_durations) + # continue mfreq_data = {} cell_means = [] model_means = [] @@ -196,24 +214,32 @@ def sam_analysis(fit_path): for i in range(len(delta_freqs)): if delta_freqs[i] != mod_freq: continue - + if len(spiketimes[i]) == 0: + print("No spiketimes found at index!") + continue if len(spiketimes[i]) > 1: print("There are more spiketimes in one 'point'! Only the first was used! ") + + spikes = spiketimes[i][0] cell_pdf = spiketimes_calculate_pdf(spikes, step_size) - cell_cuts = cut_pdf_into_periods(cell_pdf, 1/mod_freq, step_size, factor=1.0) + cell_cuts = cut_pdf_into_periods(cell_pdf, 1/mod_freq, step_size) cell_mean = np.mean(cell_cuts, axis=0) cell_means.append(cell_mean) stimulus = SAM(eod_freq, contrasts[i] / 100, mod_freq) - v1, spikes_model = model.simulate(stimulus, durations[i] * 4) + v1, spikes_model = model.simulate(stimulus, 10) model_pdf = spiketimes_calculate_pdf(spikes_model, step_size) - model_cuts = cut_pdf_into_periods(model_pdf, 1/mod_freq, step_size, factor=1.0) + model_cuts = cut_pdf_into_periods(model_pdf, 1/mod_freq, step_size) model_mean = np.mean(model_cuts, axis=0) model_means.append(model_mean) + min_length = min(min([len(cm) for cm in cell_means]), min([len(mm) for mm in model_means])) + for i in range(len(cell_means)): + cell_means[i] = cell_means[i][:min_length] + model_means[i] = model_means[i][:min_length] final_cell_mean = np.mean(cell_means, axis=0) final_model_mean = np.mean(model_means, axis=0) cell_stds.append(np.std(final_cell_mean)) @@ -225,53 +251,53 @@ def sam_analysis(fit_path): final_model_mean_phase_corrected = np.roll(final_model_mean, approx_offset) # PLOT EVERY MOD FREQ - fig, axes = plt.subplots(1, 5, figsize=(15, 5), sharex=True) - for c in cell_means: - axes[0].plot(c, color="grey", alpha=0.2) - axes[0].plot(np.mean(cell_means, axis=0), color="black") - axes[0].set_title("Cell response") - axis_cell = axes[0].axis() - - for m in model_means: - axes[1].plot(m, color="grey", alpha=0.2) - axes[1].plot(np.mean(model_means, axis=0), color="black") - axes[1].set_title("Model response") - axis_model = axes[1].axis() - ylim_top = max(axis_cell[3], axis_model[3]) - axes[1].set_ylim(0, ylim_top) - axes[0].set_ylim(0, ylim_top) - axes[2].set_ylim(0, ylim_top) - - - axes[2].plot(final_cell_mean, label="cell") - axes[2].plot(final_model_mean, label="model") - axes[2].plot(final_model_mean_phase_corrected, label="model p-cor") - axes[2].legend() - axes[2].set_title("cell-model overlapped") - axes[3].plot((final_model_mean - final_cell_mean) / final_cell_mean, label="normal") - axes[3].plot((final_model_mean_phase_corrected- final_cell_mean) / final_cell_mean, label="phase cor") - axes[3].set_title("rel. error") - axes[3].legend() - axes[4].plot(final_model_mean - final_cell_mean, label="normal") - axes[4].plot(final_model_mean_phase_corrected - final_cell_mean, label="phase cor") - axes[4].set_title("abs. error (Hz)") - axes[4].legend() - - fig.suptitle("modulation frequency: {}".format(mod_freq)) - - # plt.tight_layout() - plt.show() - plt.close() + # fig, axes = plt.subplots(1, 5, figsize=(15, 5), sharex=True) + # for c in cell_means: + # axes[0].plot(c, color="grey", alpha=0.2) + # axes[0].plot(np.mean(cell_means, axis=0), color="black") + # axes[0].set_title("Cell response") + # axis_cell = axes[0].axis() + # + # for m in model_means: + # axes[1].plot(m, color="grey", alpha=0.2) + # axes[1].plot(np.mean(model_means, axis=0), color="black") + # axes[1].set_title("Model response") + # axis_model = axes[1].axis() + # ylim_top = max(axis_cell[3], axis_model[3]) + # axes[1].set_ylim(0, ylim_top) + # axes[0].set_ylim(0, ylim_top) + # axes[2].set_ylim(0, ylim_top) + # + # axes[2].plot(final_cell_mean, label="cell") + # axes[2].plot(final_model_mean, label="model") + # axes[2].plot(final_model_mean_phase_corrected, label="model p-cor") + # axes[2].legend() + # axes[2].set_title("cell-model overlapped") + # axes[3].plot((final_model_mean - final_cell_mean) / final_cell_mean, label="normal") + # axes[3].plot((final_model_mean_phase_corrected- final_cell_mean) / final_cell_mean, label="phase cor") + # axes[3].set_title("rel. error") + # axes[3].legend() + # axes[4].plot(final_model_mean - final_cell_mean, label="normal") + # axes[4].plot(final_model_mean_phase_corrected - final_cell_mean, label="phase cor") + # axes[4].set_title("abs. error (Hz)") + # axes[4].legend() + # + # fig.suptitle("modulation frequency: {}".format(mod_freq)) + # + # # plt.tight_layout() + # # plt.show() + # plt.close() fig, ax = plt.subplots(1, 1) - ax.plot(u_delta_freqs, cell_stds, label="cell stds") - ax.plot(u_delta_freqs, model_stds, label="model stds") + ax.plot(u_delta_freqs[-len(cell_stds):], cell_stds, label="cell stds") + ax.plot(u_delta_freqs[-len(model_stds):], model_stds, label="model stds") ax.set_title("response modulation depth") ax.set_xlabel("Modulation frequency") ax.set_ylabel("STD") ax.legend() - plt.show() + plt.savefig("figures/sam/" + cell_data.get_cell_name() + ".png") + # plt.show() plt.close() @@ -335,14 +361,16 @@ def approximate_axon_delay_in_idx(cell_data, model): cell_pdf = spiketimes_calculate_pdf(spikes, step_size) - cell_cuts = cut_pdf_into_periods(cell_pdf, 1/mod_freq, step_size, factor=1.0) + cell_cuts = cut_pdf_into_periods(cell_pdf, 1/mod_freq, step_size) + if len(cell_cuts) == 0: + continue cell_mean = np.mean(cell_cuts, axis=0) cell_means.append(cell_mean) stimulus = SAM(eod_freq, contrasts[i] / 100, mod_freq) v1, spikes_model = model.simulate(stimulus, durations[i] * 4) model_pdf = spiketimes_calculate_pdf(spikes_model, step_size) - model_cuts = cut_pdf_into_periods(model_pdf, 1/mod_freq, step_size, factor=1.0) + model_cuts = cut_pdf_into_periods(model_pdf, 1/mod_freq, step_size) model_mean = np.mean(model_cuts, axis=0) model_means.append(model_mean) @@ -355,7 +383,10 @@ def approximate_axon_delay_in_idx(cell_data, model): axon_delays.append(offset) mean_delay = np.mean(axon_delays) - return int(round(mean_delay)) + if np.isnan(mean_delay): + return 0 + else: + return int(round(mean_delay)) def generate_pdf(model, stimulus, trials=4, sim_length=3, kernel_width=0.005): @@ -393,25 +424,27 @@ def spiketimes_calculate_pdf(spikes, step_size, kernel_width=0.001): return rate -def cut_pdf_into_periods(pdf, period, step_size, factor=1.5): +def cut_pdf_into_periods(pdf, period, step_size, factor=0.0): if period < 0: - print("cut_pdf_into_periods(): Period was negative! Absolute value taken to continue") + # print("cut_pdf_into_periods(): Period was negative! Absolute value taken to continue") period = abs(period) - if period / step_size > len(pdf): - return [pdf] - - idx_period_length = int(period/float(step_size)) - offset_per_step = period/float(step_size) - idx_period_length - cut_length = int(period / float(step_size) * factor) - num_of_cuts = int(len(pdf) / (idx_period_length+offset_per_step)) + idx_period_length = int(period / float(step_size)) + offset_per_step = period / float(step_size) - idx_period_length + cut_length = idx_period_length + int(factor * idx_period_length) + num_of_cuts = int(len(pdf) / (idx_period_length + offset_per_step)) if len(pdf) - (num_of_cuts * idx_period_length + (num_of_cuts * offset_per_step)) < cut_length - idx_period_length: num_of_cuts -= 1 - if num_of_cuts <= 1: - raise RuntimeError("Probability density function to short to cut.") + if idx_period_length * 0.9 > len(pdf): + return [] + # raise RuntimeError("SAM stimulus is too short for the given mod freq period.") + + if cut_length > len(pdf) or num_of_cuts < 1: + return [pdf] + cuts = np.zeros((num_of_cuts-1, cut_length)) for i in np.arange(1, num_of_cuts, 1): offset_correction = int(offset_per_step * i)