From 814e223134c9a5ef652b1286d4930eff9bb6a645 Mon Sep 17 00:00:00 2001 From: Jan Grewe Date: Tue, 29 Sep 2020 17:03:31 +0200 Subject: [PATCH] add chirpsize to the loops --- chirps_as_probing_signals.md | 2 + nix_util.py | 9 ++- plots.py | 2 +- punit_responses.py | 31 +++++----- response_discriminability.py | 113 ++++++++++++++++++++++++----------- 5 files changed, 104 insertions(+), 53 deletions(-) diff --git a/chirps_as_probing_signals.md b/chirps_as_probing_signals.md index 4943425..c36f5f4 100644 --- a/chirps_as_probing_signals.md +++ b/chirps_as_probing_signals.md @@ -57,3 +57,5 @@ Won't do, this is trivial?! * Raab et al show this is also the case with rises. * Check role of AFRs and rises in Tallarovic et al, Hupe et al. * we actually do not observe chirps without stimulation + +### correlated discrimunation results with different cell characteristics. diff --git a/nix_util.py b/nix_util.py index e212413..1c82414 100644 --- a/nix_util.py +++ b/nix_util.py @@ -14,6 +14,7 @@ def sort_blocks(nix_file): block_map = {} contrasts = [] deltafs = [] + chirp_sizes = [] conditions = [] for b in nix_file.blocks: if "baseline" not in b.name.lower(): @@ -27,10 +28,14 @@ def sort_blocks(nix_file): dltf = float(name_parts[5]) if dltf not in deltafs: deltafs.append(dltf) - block_map[(cntrst, dltf, cndtn)] = b + chirpsize = int(name_parts[7]) + if chirpsize not in chirp_sizes: + chirp_sizes.append(chirpsize) + + block_map[(cntrst, dltf, chirpsize, cndtn)] = b else: block_map["baseline"] = b - return block_map, contrasts, deltafs, conditions + return block_map, contrasts, deltafs, chirp_sizes, conditions def get_spikes(block): diff --git a/plots.py b/plots.py index fa2ef23..f8621b6 100644 --- a/plots.py +++ b/plots.py @@ -308,7 +308,7 @@ def performance_plot(args): errors[i] = np.std(df.auc[(df.kernel_width == kernel_width) & (df.contrast == c) & (df.df == d) & (df.detection_task == t)]) df_ax.errorbar(dfs, performances, yerr=errors, fmt=".-", label="%.2f" % c) df_ax.set_ylim([0.25, 1.0]) - df_ax.set_ylabel("performance", fontsize=8) + df_ax.set_ylabel("performance", fontsize=8, rotation=180) df_ax.set_xlabel(r"$\Delta_f$ [Hz]", fontsize=8) df_ax.hlines(0.5, dfs[0], dfs[-1], color="k", ls="--", lw=0.2) df_ax.legend(fontsize=7, ncol=4, frameon=False, loc="lower center", mode="expand", handlelength=1.0, handletextpad=0.25) diff --git a/punit_responses.py b/punit_responses.py index 698ae13..68ec34b 100644 --- a/punit_responses.py +++ b/punit_responses.py @@ -149,7 +149,7 @@ def simulate_responses(stimulus_params, model_params, repeats=10, deltaf=20): del cell_params["cell"] del cell_params["EODf"] conditions = ["other", "self"] - + chirp_size = stimulus_params["chirp_size"] pre_time, pre_stim = get_pre_stimulus(stimulus_params["eodfs"]["self"], dt=model_params["deltat"]) for contrast in stimulus_params["contrasts"]: params = stimulus_params.copy() @@ -159,7 +159,7 @@ def simulate_responses(stimulus_params, model_params, repeats=10, deltaf=20): for condition in conditions: print("\tcontrast: %s, condition: %s" %(contrast, condition), " "*10, end="\r") - block_name = "contrast_%.3f_condition_%s_deltaf_%i" %(contrast, condition, deltaf) + block_name = "contrast_%.3f_condition_%s_deltaf_%i_chirpsize_%i" %(contrast, condition, deltaf, chirp_size) params["condition"] = condition time, self_signal, self_freq, other_signal, other_freq = get_signals(**params) full_signal = (self_signal + other_signal) @@ -177,7 +177,7 @@ def simulate_responses(stimulus_params, model_params, repeats=10, deltaf=20): sp = simulate(np.hstack((pre_stim, self_signal)), **cell_params) no_other_spikes.append(sp[sp > pre_time[-1]] - pre_time[-1]) if condition == "self": - name = "contrast_%.3f_condition_no-other_deltaf_%i" %(contrast, deltaf) + name = "contrast_%.3f_condition_no-other_deltaf_%i_chirpsize_%i" %(contrast, deltaf, chirp_size) save(filename, name, params, cell_params, self_signal, None, self_freq, None, self_signal, no_other_spikes) save(filename, block_name, params, cell_params, self_signal, other_signal, self_freq, other_freq, full_signal, spikes) print("\n") @@ -185,6 +185,7 @@ def simulate_responses(stimulus_params, model_params, repeats=10, deltaf=20): def simulate_cell(cell_id, models): deltafs = [-200, -100, -50, -20, -10, -5, 5, 10, 20, 50, 100, 200] # Hz, difference frequency between self and other + chirp_sizes = [40, 60, 100] stimulus_params = { "eodfs": {"self": 0.0, "other": 0.0}, # eod frequency in Hz, to be overwritten "contrasts": [20, 10, 5, 2.5, 1.25, 0.625, 0.3125], "chirp_size": 100, # Hz, frequency excursion @@ -200,23 +201,25 @@ def simulate_cell(cell_id, models): save_baseline_response(filename, "baseline response", baseline_spikes, model_params) print("Cell: %s" % model_params["cell"]) - for deltaf in deltafs: - stimulus_params["eodfs"] = {"self": model_params["EODf"], "other": model_params["EODf"] + deltaf} - stimulus_params["dt"] = model_params["deltat"] - - print("\t Deltaf: %i" % deltaf) - chirp_times = np.arange(stimulus_params["chirp_duration"], - stimulus_params["duration"] - stimulus_params["chirp_duration"], - 1./stimulus_params["chirp_frequency"]) - stimulus_params["chirp_times"] = chirp_times - simulate_responses(stimulus_params, model_params, repeats=25, deltaf=deltaf) + for cs in chirp_sizes: + for deltaf in deltafs: + stimulus_params["eodfs"] = {"self": model_params["EODf"], "other": model_params["EODf"] + deltaf} + stimulus_params["dt"] = model_params["deltat"] + stimulus_params["chirp_size"] = cs + + print("\t Deltaf: %i" % deltaf) + chirp_times = np.arange(stimulus_params["chirp_duration"], + stimulus_params["duration"] - stimulus_params["chirp_duration"], + 1./stimulus_params["chirp_frequency"]) + stimulus_params["chirp_times"] = chirp_times + simulate_responses(stimulus_params, model_params, repeats=25, deltaf=deltaf) def main(): models = load_models("models.csv") num_cores = multiprocessing.cpu_count() - 6 - Parallel(n_jobs=num_cores)(delayed(simulate_cell)(cell_id, models) for cell_id in range(len(models[:10]))) + Parallel(n_jobs=num_cores)(delayed(simulate_cell)(cell_id, models) for cell_id in range(len(models[:20]))) if __name__ == "__main__": diff --git a/response_discriminability.py b/response_discriminability.py index 9d9173c..a2e3c08 100644 --- a/response_discriminability.py +++ b/response_discriminability.py @@ -60,17 +60,32 @@ def get_firing_rate(block_map, df, contrast, condition, kernel_width=0.0005): return time, rates, spikes +def foreign_fish_detection_beat(block_map, df, cs, all_contrasts, all_conditions, kernel_width=0.0005, cell_name="", store_roc=False): + """Tries to detect the presence of a foreign fish by estimating the discriminability of the responses during the beat + versus the responses without another fish beeing there, i.e. the baseline activity. + Applies a ROC analysis to the response segments between chirps. Calculates a) the distances between the baseline responses and + b) distances between the baseline and beat responses. Tests whether distances in b) are larger than a) + Args: + block_map ([type]): maps nix blocks to combination of stimulus parameters + df ([type]): the difference frequency that should be used + cs ([type]): ths chirpsize that should be used + all_contrasts ([type]): list of all used contrasts + all_conditions ([type]): list of all chirp conditions, i.e. self, other, or no-other + kernel_width (float, optional): std of Gaussian kernel. Defaults to 0.0005. + cell_name (str, optional): name of the cell. Defaults to "". + store_roc (bool, optional): if true the full false positives and true positives will be returned leads to huge file sizes!. Defaults to False. - -def foreign_fish_detection_beat(block_map, df, all_contrasts, all_conditions, kernel_width=0.0005, cell_name="", store_roc=False): + Returns: + list of dictionaries: the results, auc is the area under the curve, i.e. the discrimination performance in the range [0, 1]. The 'detection_task' is 'beat' + """ detection_performances = [] for contrast in all_contrasts: print(" " * 50, end="\r") print("Contrast: %.3f" % contrast, end="\r") - no_other_block = block_map[(contrast, df, "no-other")] - self_block = block_map[(contrast, df, "self")] + no_other_block = block_map[(contrast, df, cs, "no-other")] + self_block = block_map[(contrast, df, cs, "self")] # get some metadata assuming they are all the same for each condition, which they should duration, dt, _, chirp_duration, chirp_times = get_chirp_metadata(self_block) @@ -111,26 +126,53 @@ def foreign_fish_detection_beat(block_map, df, all_contrasts, all_conditions, ke temp2 = np.ones_like(valid_distances_comparison) group = np.hstack((temp1, temp2)) - score = np.hstack((valid_distances_baseline, valid_distances_comparison)) - fpr, tpr, _ = roc_curve(group, score, pos_label=1) + score = np.hstack((valid_distances_baseline, valid_distances_comparison)) auc = roc_auc_score(group, score) if store_roc: - detection_performances.append({"cell": cell_name, "detection_task": "beat", "contrast": contrast, "df": df, "kernel_width": kernel_width, "auc": auc, "true_positives": tpr, "false_positives": fpr}) + fpr, tpr, _ = roc_curve(group, score, pos_label=1) + detection_performances.append({"cell": cell_name, "detection_task": "beat", "contrast": contrast, "df": df, "kernel_width": kernel_width, "chirpsize": cs, "auc": auc, "true_positives": tpr, "false_positives": fpr}) else: - detection_performances.append({"cell": cell_name, "detection_task": "beat", "contrast": contrast, "df": df, "kernel_width": kernel_width, "auc": auc}) + detection_performances.append({"cell": cell_name, "detection_task": "beat", "contrast": contrast, "df": df, "kernel_width": kernel_width, "chirpsize": cs, "auc": auc}) print("\n") return detection_performances -def foreign_fish_detection_chirp(block_map, df, all_contrasts, all_conditions, kernel_width=0.0005, cell_name="", store_roc=False): +def foreign_fish_detection_chirp(block_map, df, cs, all_contrasts, all_conditions, kernel_width=0.0005, cell_name="", store_roc=False): + """Tries to detect the presence of a foreign fish by estimating the discriminability of the chirp + responses in the presence of another fish versus the responses without another fish beeing around. + + Applies a ROC analysis to the response segments containing the chirp. Does two discrimination tests: + 1) compares the responses to self-chirping alone to the responses to self-chriping in company. + 2) compares the responess to other-chirping to the response during the beat. + + Tests the assumptions that the distances a) between the self-chriping alone and self-chriping in company + are larger than the distances within the the self-chirping alone condition and b) the distances between + other-chirping in company and no one is chirping in company (i.e. beat) are larger than the distances + within the beat responses. + + Args: + block_map ([type]): maps nix blocks to combination of stimulus parameters + df ([type]): the difference frequency that should be used + cs ([type]): ths chirpsize that should be used + all_contrasts ([type]): list of all used contrasts + all_conditions ([type]): list of all chirp conditions, i.e. self, other, or no-other + kernel_width (float, optional): std of Gaussian kernel. Defaults to 0.0005. + cell_name (str, optional): name of the cell. Defaults to "". + store_roc (bool, optional): if true the full false positives and true positives will be returned leads to huge file sizes!. Defaults to False. + + Returns: + list of dictionaries: the results, auc is the area under the curve, i.e. the discrimination performance in the range [0, 1]. + The 'detection_task' is either "self vs soliloquy" for 1) or "other vs quietness" for 2) + + """ detection_performances = [] for contrast in all_contrasts: print(" " * 50, end="\r") print("Contrast: %.3f" % contrast, end="\r") - no_other_block = block_map[(contrast, df, "no-other")] - self_block = block_map[(contrast, df, "self")] - other_block = block_map[(contrast, df, "self")] + no_other_block = block_map[(contrast, df, cs, "no-other")] + self_block = block_map[(contrast, df, cs, "self")] + other_block = block_map[(contrast, df, cs, "self")] # get some metadata assuming they are all the same for each condition, which they should duration, dt, _, chirp_duration, chirp_times = get_chirp_metadata(self_block) @@ -191,54 +233,52 @@ def foreign_fish_detection_chirp(block_map, df, all_contrasts, all_conditions, k group = np.hstack((no_other_temp, self_vs_alone_temp)) score = np.hstack((valid_no_other_distances, valid_self_vs_alone_distances)) - fpr, tpr, _ = roc_curve(group, score, pos_label=1) auc = roc_auc_score(group, score) if store_roc: - detection_performances.append({"cell": cell_name, "detection_task": "self vs soliloquy", "contrast": contrast, "df": df, "kernel_width": kernel_width, "auc": auc, "true_positives": tpr, "false_positives": fpr}) + fpr, tpr, _ = roc_curve(group, score, pos_label=1) + detection_performances.append({"cell": cell_name, "detection_task": "self vs soliloquy", "contrast": contrast, "df": df, "kernel_width": kernel_width, "chirpsize": cs, "auc": auc, "true_positives": tpr, "false_positives": fpr}) else: - detection_performances.append({"cell": cell_name, "detection_task": "self vs soliloquy", "contrast": contrast, "df": df, "kernel_width": kernel_width, "auc": auc}) + detection_performances.append({"cell": cell_name, "detection_task": "self vs soliloquy", "contrast": contrast, "df": df, "kernel_width": kernel_width, "chirpsize": cs, "auc": auc}) group = np.hstack((silence_temp, other_vs_silence_temp)) score = np.hstack((valid_silence_distances, valid_other_vs_silence_distances)) - fpr, tpr, _ = roc_curve(group, score, pos_label=1) auc = roc_auc_score(group, score) if store_roc: - detection_performances.append({"cell": cell_name, "detection_task": "other vs quietness", "contrast": contrast, "df": df, "kernel_width": kernel_width, "auc": auc, "true_positives": tpr, "false_positives": fpr}) + fpr, tpr, _ = roc_curve(group, score, pos_label=1) + detection_performances.append({"cell": cell_name, "detection_task": "other vs quietness", "contrast": contrast, "df": df, "kernel_width": kernel_width, "chirpsize": cs, "auc": auc, "true_positives": tpr, "false_positives": fpr}) else: - detection_performances.append({"cell": cell_name, "detection_task": "other vs quietness", "contrast": contrast, "df": df, "kernel_width": kernel_width, "auc": auc}) + detection_performances.append({"cell": cell_name, "detection_task": "other vs quietness", "contrast": contrast, "df": df, "kernel_width": kernel_width, "chirpsize": cs, "auc": auc}) print("\n") return detection_performances -def foreign_fish_detection(block_map, all_dfs, all_contrasts, all_conditions, current_df=None, cell_name="", store_roc=False): - dfs = [current_df] if current_df is not None else all_dfs +def foreign_fish_detection(block_map, all_dfs, all_contrasts, all_conditions, all_chirpsizes, current_df=None, current_chirpsize=None, cell_name="", store_roc=False): + dfs = [current_df] if current_df is not None else all_dfs + chirp_sizes = [current_chirpsize] if current_chirpsize is not None else all_chirpsizes kernels = [0.00025, 0.0005, 0.001, 0.0025] result_dicts = [] - for df in dfs: - for kw in kernels: - print("df: %i, kernel: %.4f" % (df, kw)) - print("Foreign fish detection during beat:") - result_dicts.extend(foreign_fish_detection_beat(block_map, df, all_contrasts, all_conditions, kw, cell_name, store_roc)) - print("Foreign fish detection during chirp:") - result_dicts.extend(foreign_fish_detection_chirp(block_map, df, all_contrasts, all_conditions, kw, cell_name, store_roc)) + for cs in chirp_sizes: + for df in dfs: + for kw in kernels: + print("cs: %i Hz, df: %i Hz, kernel: %.4fs" % (cs, df, kw)) + print("Foreign fish detection during beat:") + result_dicts.extend(foreign_fish_detection_beat(block_map, df, cs, all_contrasts, all_conditions, kw, cell_name, store_roc)) + print("Foreign fish detection during chirp:") + result_dicts.extend(foreign_fish_detection_chirp(block_map, df, cs, all_contrasts, all_conditions, kw, cell_name, store_roc)) return result_dicts -def estimate_chirp_phase(am, chirp_times): - - pass - - -def process_cell(filename, dfs=[], contrasts=[], conditions=[]): +def process_cell(filename): print(filename) nf = nix.File.open(filename, nix.FileMode.ReadOnly) - block_map, all_contrasts, all_dfs, all_conditions = sort_blocks(nf) + block_map, all_contrasts, all_dfs, all_chirpsizes, all_conditions = sort_blocks(nf) if "baseline" in block_map.keys(): baseline_spikes = read_baseline(block_map["baseline"]) else: print("ERROR: no baseline data for file %s!" % filename) - results = foreign_fish_detection(block_map, all_dfs, all_contrasts, all_conditions, current_df=None, + results = foreign_fish_detection(block_map, all_dfs, all_contrasts, all_conditions, all_chirpsizes, + current_df=None, current_chirpsize=None, cell_name=filename.split(os.path.sep)[-1].split(".nix")[0], store_roc=False) nf.close() return results @@ -253,7 +293,8 @@ def main(): for pr in processed_list: results.extend(pr) df = pd.DataFrame(results) - df.to_csv(os.path.join(data_folder, "discimination_results.csv"), sep=";") + df.to_csv(os.path.join(data_folder, "discimination_results2.csv"), sep=";") + if __name__ == "__main__": main() \ No newline at end of file