From 87044903f5541ba0c86ad9e3fccaf40979df3103 Mon Sep 17 00:00:00 2001 From: bloergl <56762962+bloergl@users.noreply.github.com> Date: Wed, 18 Dec 2024 18:43:12 +0100 Subject: [PATCH] Add files via upload --- geci_loader.py | 17 +++--- geci_plot.py | 36 ++++++------ stage_4_process.py | 140 +++++++++++++++++++++++---------------------- 3 files changed, 100 insertions(+), 93 deletions(-) diff --git a/geci_loader.py b/geci_loader.py index 8d530f5..a8f4da1 100644 --- a/geci_loader.py +++ b/geci_loader.py @@ -18,7 +18,7 @@ def func_exp(x, a, b, c): def loader( filename: str = "config_M_Sert_Cre_49.json", - fpath: str = "/data_1/hendrik/gevi", + fpath: str|None = None, skip_timesteps: int = 100, # If there is no special ROI... Get one! This is just a backup roi_control_path_default: str = "roi_controlM_Sert_Cre_49.npy", @@ -27,6 +27,9 @@ def loader( fit_power: bool = False, # True => -ax^b ; False => exp(-b) ) -> None: + if fpath is None: + fpath = os.getcwd() + if os.path.isfile(filename) is False: print(f"{filename} is missing") exit() @@ -42,8 +45,8 @@ def loader( ) if remove_fit: - roi_control_path: str = f"roi_control{config["mouse_identifier"]}.npy" - roi_sdarken_path: str = f"roi_sdarken{config["mouse_identifier"]}.npy" + roi_control_path: str = f"roi_control{config['mouse_identifier']}.npy" + roi_sdarken_path: str = f"roi_sdarken{config['mouse_identifier']}.npy" if os.path.isfile(roi_control_path) is False: print(f"Using replacement RIO: {roi_control_path_default}") @@ -72,7 +75,7 @@ def loader( ) tmp_fname = os.path.join( fpath, - "output_" + config["mouse_identifier"], + config["export_path"], experiment_name + "_acceptor_donor.npz", ) print(f'Processing file "{tmp_fname}"...') @@ -156,9 +159,9 @@ def loader( light_signal[i_exp] += tmp_light_signal data_sequence[i_exp] /= n_tri light_signal[i_exp] /= n_tri - np.save("dsq_" + config["mouse_identifier"], data_sequence) - np.save("lsq_" + config["mouse_identifier"], light_signal) - np.save("msq_" + config["mouse_identifier"], mask) + np.save(os.path.join(fpath, config["export_path"], "dsq_" + config["mouse_identifier"]), data_sequence) + np.save(os.path.join(fpath, config["export_path"], "lsq_" + config["mouse_identifier"]), light_signal) + np.save(os.path.join(fpath, config["export_path"], "msq_" + config["mouse_identifier"]), mask) if __name__ == "__main__": diff --git a/geci_plot.py b/geci_plot.py index a5e4717..8d32ab4 100644 --- a/geci_plot.py +++ b/geci_plot.py @@ -1,3 +1,5 @@ +# %% + import numpy as np import matplotlib.pyplot as plt import argh @@ -18,15 +20,16 @@ def func_exp(x, a, b, c): # mouse: int = 0, 1, 2, 3, 4 def plot( filename: str = "config_M_Sert_Cre_49.json", + fpath: str | None = None, experiment: int = 4, skip_timesteps: int = 100, - # If there is no special ROI... Get one! This is just a backup - roi_control_path_default: str = "roi_controlM_Sert_Cre_49.npy", - roi_sdarken_path_default: str = "roi_sdarkenM_Sert_Cre_49.npy", remove_fit: bool = False, fit_power: bool = False, # True => -ax^b ; False => exp(-b) ) -> None: + if fpath is None: + fpath = os.getcwd() + if os.path.isfile(filename) is False: print(f"{filename} is missing") exit() @@ -45,30 +48,25 @@ def plot( print(f"ERROR: could not find raw directory {raw_data_path}!!!!") exit() - with open(f"meta_{config["mouse_identifier"]}_exp{experiment:03d}.json", "r") as file: + with open(f"meta_{config['mouse_identifier']}_exp{experiment:03d}.json", "r") as file: metadata = json.loads(jsmin(file.read())) experiment_names = metadata['sessionMetaData']['experimentNames'][str(experiment)] - roi_control_path: str = f"roi_control{config["mouse_identifier"]}.npy" - roi_sdarken_path: str = f"roi_sdarken{config["mouse_identifier"]}.npy" + roi_control_path: str = f"roi_control{config['mouse_identifier']}.npy" + roi_sdarken_path: str = f"roi_sdarken{config['mouse_identifier']}.npy" - if os.path.isfile(roi_control_path) is False: - print(f"Using replacement RIO: {roi_control_path_default}") - roi_control_path = roi_control_path_default - - if os.path.isfile(roi_sdarken_path) is False: - print(f"Using replacement RIO: {roi_sdarken_path_default}") - roi_sdarken_path = roi_sdarken_path_default + assert os.path.isfile(roi_control_path) + assert os.path.isfile(roi_sdarken_path) print("Load data...") - data = np.load("dsq_" + config["mouse_identifier"] + ".npy", mmap_mode="r") + data = np.load(os.path.join(fpath, config["export_path"], "dsq_" + config["mouse_identifier"] + ".npy"), mmap_mode="r") print("Load light signal...") - light = np.load("lsq_" + config["mouse_identifier"] + ".npy", mmap_mode="r") + light = np.load(os.path.join(fpath, config["export_path"], "lsq_" + config["mouse_identifier"] + ".npy"), mmap_mode="r") print("Load mask...") - mask = np.load("msq_" + config["mouse_identifier"] + ".npy") + mask = np.load(os.path.join(fpath, config["export_path"], "msq_" + config["mouse_identifier"] + ".npy")) roi_control = np.load(roi_control_path) roi_control *= mask @@ -82,14 +80,14 @@ def plot( a_show = data[experiment - 1, :, :, 1000].copy() a_show[(roi_darken + roi_control) < 0.5] = np.nan plt.imshow(a_show) - plt.title(f"{config["mouse_identifier"]} -- Experiment: {experiment}") + plt.title(f"{config['mouse_identifier']} -- Experiment: {experiment}") plt.show(block=False) plt.figure(2) a_dontshow = data[experiment - 1, :, :, 1000].copy() a_dontshow[(roi_darken + roi_control) > 0.5] = np.nan plt.imshow(a_dontshow) - plt.title(f"{config["mouse_identifier"]} -- Experiment: {experiment}") + plt.title(f"{config['mouse_identifier']} -- Experiment: {experiment}") plt.show(block=False) plt.figure(3) @@ -174,7 +172,7 @@ def plot( plt.plot(time_axis, light_signal, c="k", label="light") plt.plot(time_axis, darken, label="sDarken") plt.plot(time_axis, lighten, label="control") - plt.title(f"{config["mouse_identifier"]} -- Experiment: {experiment} ({experiment_names})") + plt.title(f"{config['mouse_identifier']} -- Experiment: {experiment} ({experiment_names})") plt.legend() plt.show() diff --git a/stage_4_process.py b/stage_4_process.py index b8e4382..2bdf57d 100644 --- a/stage_4_process.py +++ b/stage_4_process.py @@ -1,3 +1,5 @@ +#%% + import numpy as np import torch import torchvision as tv # type: ignore @@ -568,84 +570,88 @@ def process_trial( data -= heartbeat_coefficients.unsqueeze(1) * volume_heartbeat.unsqueeze( 0 ).movedim(-1, 1) + # data_herzlos = data.clone() mylogger.info("-==- Done -==-") - donor_heartbeat_factor = heartbeat_coefficients[donor_index, ...].clone() - acceptor_heartbeat_factor = heartbeat_coefficients[acceptor_index, ...].clone() - del heartbeat_coefficients + if config["gevi"]: # UDO scaling performed! + donor_heartbeat_factor = heartbeat_coefficients[donor_index, ...].clone() + acceptor_heartbeat_factor = heartbeat_coefficients[acceptor_index, ...].clone() + del heartbeat_coefficients - if device != torch.device("cpu"): - torch.cuda.empty_cache() - mylogger.info("Empty CUDA cache") - free_mem = cuda_total_memory - max( - [ - torch.cuda.memory_reserved(device), - torch.cuda.memory_allocated(device), - ] - ) - mylogger.info(f"CUDA memory: {free_mem // 1024} MByte") + if device != torch.device("cpu"): + torch.cuda.empty_cache() + mylogger.info("Empty CUDA cache") + free_mem = cuda_total_memory - max( + [ + torch.cuda.memory_reserved(device), + torch.cuda.memory_allocated(device), + ] + ) + mylogger.info(f"CUDA memory: {free_mem // 1024} MByte") - mylogger.info("Calculate scaling factor for donor and acceptor") - donor_factor: torch.Tensor = ( - donor_heartbeat_factor + acceptor_heartbeat_factor - ) / (2 * donor_heartbeat_factor) - acceptor_factor: torch.Tensor = ( - donor_heartbeat_factor + acceptor_heartbeat_factor - ) / (2 * acceptor_heartbeat_factor) + mylogger.info("Calculate scaling factor for donor and acceptor") + donor_factor: torch.Tensor = ( + donor_heartbeat_factor + acceptor_heartbeat_factor + ) / (2 * donor_heartbeat_factor) + acceptor_factor: torch.Tensor = ( + donor_heartbeat_factor + acceptor_heartbeat_factor + ) / (2 * acceptor_heartbeat_factor) - del donor_heartbeat_factor - del acceptor_heartbeat_factor + del donor_heartbeat_factor + del acceptor_heartbeat_factor - if config["save_factors"]: - temp_path = os.path.join( - config["export_path"], experiment_name + "_donor_factor.npy" - ) - mylogger.info(f"Save donor factor to {temp_path}") - np.save(temp_path, donor_factor.cpu()) + if config["save_factors"]: + temp_path = os.path.join( + config["export_path"], experiment_name + "_donor_factor.npy" + ) + mylogger.info(f"Save donor factor to {temp_path}") + np.save(temp_path, donor_factor.cpu()) - temp_path = os.path.join( - config["export_path"], experiment_name + "_acceptor_factor.npy" - ) - mylogger.info(f"Save acceptor factor to {temp_path}") - np.save(temp_path, acceptor_factor.cpu()) - mylogger.info("-==- Done -==-") + temp_path = os.path.join( + config["export_path"], experiment_name + "_acceptor_factor.npy" + ) + mylogger.info(f"Save acceptor factor to {temp_path}") + np.save(temp_path, acceptor_factor.cpu()) + mylogger.info("-==- Done -==-") - mylogger.info("Scale acceptor to heart beat amplitude") - mylogger.info("Calculate mean") - mean_values_acceptor = data[ - acceptor_index, config["skip_frames_in_the_beginning"] :, ... - ].nanmean(dim=0, keepdim=True) + mylogger.info("Scale acceptor to heart beat amplitude") + mylogger.info("Calculate mean") + mean_values_acceptor = data[ + acceptor_index, config["skip_frames_in_the_beginning"] :, ... + ].nanmean(dim=0, keepdim=True) - mylogger.info("Remove mean") - data[acceptor_index, ...] -= mean_values_acceptor - mylogger.info("Apply acceptor_factor and mask") - data[acceptor_index, ...] *= acceptor_factor.unsqueeze( - 0 - ) * mask_positve.unsqueeze(0) - mylogger.info("Add mean") - data[acceptor_index, ...] += mean_values_acceptor - mylogger.info("-==- Done -==-") + mylogger.info("Remove mean") + data[acceptor_index, ...] -= mean_values_acceptor + mylogger.info("Apply acceptor_factor and mask") + data[acceptor_index, ...] *= acceptor_factor.unsqueeze( + 0 + ) * mask_positve.unsqueeze(0) + mylogger.info("Add mean") + data[acceptor_index, ...] += mean_values_acceptor + mylogger.info("-==- Done -==-") - mylogger.info("Scale donor to heart beat amplitude") - mylogger.info("Calculate mean") - mean_values_donor = data[ - donor_index, config["skip_frames_in_the_beginning"] :, ... - ].nanmean(dim=0, keepdim=True) - mylogger.info("Remove mean") - data[donor_index, ...] -= mean_values_donor - mylogger.info("Apply donor_factor and mask") - data[donor_index, ...] *= donor_factor.unsqueeze(0) * mask_positve.unsqueeze(0) - mylogger.info("Add mean") - data[donor_index, ...] += mean_values_donor - mylogger.info("-==- Done -==-") + mylogger.info("Scale donor to heart beat amplitude") + mylogger.info("Calculate mean") + mean_values_donor = data[ + donor_index, config["skip_frames_in_the_beginning"] :, ... + ].nanmean(dim=0, keepdim=True) + mylogger.info("Remove mean") + data[donor_index, ...] -= mean_values_donor + mylogger.info("Apply donor_factor and mask") + data[donor_index, ...] *= donor_factor.unsqueeze(0) * mask_positve.unsqueeze(0) + mylogger.info("Add mean") + data[donor_index, ...] += mean_values_donor + mylogger.info("-==- Done -==-") + else: + mylogger.info("GECI does not require acceptor/donor scaling, skipping!") + mylogger.info("-==- Done -==-") - mylogger.info("Divide by mean over time") - data /= data[:, config["skip_frames_in_the_beginning"] :, ...].nanmean( - dim=1, - keepdim=True, - ) - - mylogger.info("-==- Done -==-") + mylogger.info("Divide by mean over time") + data /= data[:, config["skip_frames_in_the_beginning"] :, ...].nanmean( + dim=1, + keepdim=True, + ) + mylogger.info("-==- Done -==-") data = data.nan_to_num(nan=0.0) mylogger.info("Preparation for regression -- Gauss smear")