diff --git a/config_M3879M_2021-10-05.json b/config_M3879M_2021-10-05.json index 87ae628..da3980a 100644 --- a/config_M3879M_2021-10-05.json +++ b/config_M3879M_2021-10-05.json @@ -5,6 +5,8 @@ "raw_path": "raw", "export_path": "output_M3879M_2021-10-05", "ref_image_path": "ref_images_M3879M_2021-10-05", + "heartbeat_remove": true, // if gevi must be true; geci: who knows... + "gevi": true, // true => gevi, false => geci // Ratio Sequence "classical_ratio_mode": true, // true: a/d false: 1+a-d // Regression diff --git a/config_M_Sert_Cre_41.json b/config_M_Sert_Cre_41.json index 668bd6c..3675a3d 100644 --- a/config_M_Sert_Cre_41.json +++ b/config_M_Sert_Cre_41.json @@ -5,6 +5,8 @@ "raw_path": "raw", "export_path": "output_M_Sert_Cre_41", "ref_image_path": "ref_images_M_Sert_Cre_41", + "heartbeat_remove": false, + "gevi": false, // true => gevi, false => geci // Ratio Sequence "classical_ratio_mode": true, // true: a/d false: 1+a-d // Regression diff --git a/config_M_Sert_Cre_42.json b/config_M_Sert_Cre_42.json new file mode 100644 index 0000000..47d0ab6 --- /dev/null +++ b/config_M_Sert_Cre_42.json @@ -0,0 +1,62 @@ +{ + "basic_path": "/data_1/hendrik", + "recoding_data": "2023-07-18", + "mouse_identifier": "M_Sert_Cre_42", + "raw_path": "raw", + "export_path": "output_M_Sert_Cre_42", + "ref_image_path": "ref_images_M_Sert_Cre_42", + "heartbeat_remove": false, + "gevi": false, // true => gevi, false => geci + // Ratio Sequence + "classical_ratio_mode": true, // true: a/d false: 1+a-d + // Regression + //"target_camera_acceptor": "acceptor", + "target_camera_acceptor": "", + "regressor_cameras_acceptor": [ + "oxygenation", + "volume" + ], + "target_camera_donor": "donor", + "regressor_cameras_donor": [ + // "oxygenation", + "volume" + ], + // binning + "binning_enable": true, + "binning_at_the_end": false, + "binning_kernel_size": 4, + "binning_stride": 4, + "binning_divisor_override": 1, + // alignment + "alignment_batch_size": 200, + "rotation_stabilization_threshold_factor": 3.0, // >= 1.0 + "rotation_stabilization_threshold_border": 0.9, // <= 1.0 + // Heart beat detection + "lower_freqency_bandpass": 5.0, // Hz + "upper_freqency_bandpass": 14.0, // Hz + "heartbeat_filtfilt_chuck_size": 10, + // Gauss smear + "gauss_smear_spatial_width": 8, + "gauss_smear_temporal_width": 0.1, + "gauss_smear_use_matlab_mask": false, + // LED Ramp on + "skip_frames_in_the_beginning": 100, // Frames + // PyTorch + "dtype": "float32", + "force_to_cpu": false, + // Save + "save_as_python": true, // produces .npz files (compressed) + "save_as_matlab": false, // produces .hd5 file (compressed) + // Save extra information + "save_alignment": false, + "save_heartbeat": false, + "save_factors": false, + "save_regression_coefficients": false, + // Not important parameter + "required_order": [ + "acceptor", + "donor", + "oxygenation", + "volume" + ] +} diff --git a/config_M_Sert_Cre_45.json b/config_M_Sert_Cre_45.json new file mode 100644 index 0000000..e28c337 --- /dev/null +++ b/config_M_Sert_Cre_45.json @@ -0,0 +1,62 @@ +{ + "basic_path": "/data_1/hendrik", + "recoding_data": "2023-07-18", + "mouse_identifier": "M_Sert_Cre_45", + "raw_path": "raw", + "export_path": "output_M_Sert_Cre_45", + "ref_image_path": "ref_images_M_Sert_Cre_45", + "heartbeat_remove": false, + "gevi": false, // true => gevi, false => geci + // Ratio Sequence + "classical_ratio_mode": true, // true: a/d false: 1+a-d + // Regression + //"target_camera_acceptor": "acceptor", + "target_camera_acceptor": "", + "regressor_cameras_acceptor": [ + "oxygenation", + "volume" + ], + "target_camera_donor": "donor", + "regressor_cameras_donor": [ + // "oxygenation", + "volume" + ], + // binning + "binning_enable": true, + "binning_at_the_end": false, + "binning_kernel_size": 4, + "binning_stride": 4, + "binning_divisor_override": 1, + // alignment + "alignment_batch_size": 200, + "rotation_stabilization_threshold_factor": 3.0, // >= 1.0 + "rotation_stabilization_threshold_border": 0.9, // <= 1.0 + // Heart beat detection + "lower_freqency_bandpass": 5.0, // Hz + "upper_freqency_bandpass": 14.0, // Hz + "heartbeat_filtfilt_chuck_size": 10, + // Gauss smear + "gauss_smear_spatial_width": 8, + "gauss_smear_temporal_width": 0.1, + "gauss_smear_use_matlab_mask": false, + // LED Ramp on + "skip_frames_in_the_beginning": 100, // Frames + // PyTorch + "dtype": "float32", + "force_to_cpu": false, + // Save + "save_as_python": true, // produces .npz files (compressed) + "save_as_matlab": false, // produces .hd5 file (compressed) + // Save extra information + "save_alignment": false, + "save_heartbeat": false, + "save_factors": false, + "save_regression_coefficients": false, + // Not important parameter + "required_order": [ + "acceptor", + "donor", + "oxygenation", + "volume" + ] +} diff --git a/config_M_Sert_Cre_46.json b/config_M_Sert_Cre_46.json new file mode 100644 index 0000000..21db1d5 --- /dev/null +++ b/config_M_Sert_Cre_46.json @@ -0,0 +1,62 @@ +{ + "basic_path": "/data_1/hendrik", + "recoding_data": "2023-03-16", + "mouse_identifier": "M_Sert_Cre_46", + "raw_path": "raw", + "export_path": "output_M_Sert_Cre_46", + "ref_image_path": "ref_images_M_Sert_Cre_46", + "heartbeat_remove": false, + "gevi": false, // true => gevi, false => geci + // Ratio Sequence + "classical_ratio_mode": true, // true: a/d false: 1+a-d + // Regression + //"target_camera_acceptor": "acceptor", + "target_camera_acceptor": "", + "regressor_cameras_acceptor": [ + "oxygenation", + "volume" + ], + "target_camera_donor": "donor", + "regressor_cameras_donor": [ + // "oxygenation", + "volume" + ], + // binning + "binning_enable": true, + "binning_at_the_end": false, + "binning_kernel_size": 4, + "binning_stride": 4, + "binning_divisor_override": 1, + // alignment + "alignment_batch_size": 200, + "rotation_stabilization_threshold_factor": 3.0, // >= 1.0 + "rotation_stabilization_threshold_border": 0.9, // <= 1.0 + // Heart beat detection + "lower_freqency_bandpass": 5.0, // Hz + "upper_freqency_bandpass": 14.0, // Hz + "heartbeat_filtfilt_chuck_size": 10, + // Gauss smear + "gauss_smear_spatial_width": 8, + "gauss_smear_temporal_width": 0.1, + "gauss_smear_use_matlab_mask": false, + // LED Ramp on + "skip_frames_in_the_beginning": 100, // Frames + // PyTorch + "dtype": "float32", + "force_to_cpu": false, + // Save + "save_as_python": true, // produces .npz files (compressed) + "save_as_matlab": false, // produces .hd5 file (compressed) + // Save extra information + "save_alignment": false, + "save_heartbeat": false, + "save_factors": false, + "save_regression_coefficients": false, + // Not important parameter + "required_order": [ + "acceptor", + "donor", + "oxygenation", + "volume" + ] +} diff --git a/config_M_Sert_Cre_49.json b/config_M_Sert_Cre_49.json index 1b0e58e..2621f5e 100644 --- a/config_M_Sert_Cre_49.json +++ b/config_M_Sert_Cre_49.json @@ -5,6 +5,8 @@ "raw_path": "raw", "export_path": "output_M_Sert_Cre_49", "ref_image_path": "ref_images_M_Sert_Cre_49", + "heartbeat_remove": false, + "gevi": false, // true => gevi, false => geci // Ratio Sequence "classical_ratio_mode": true, // true: a/d false: 1+a-d // Regression diff --git a/geci_loader.py b/geci_loader.py new file mode 100644 index 0000000..95826ae --- /dev/null +++ b/geci_loader.py @@ -0,0 +1,71 @@ +import numpy as np +import os +import argh + + +# mouse:int = 0, 1, 2, 3, 4 +def loader(mouse:int = 0, fpath:str = "/data_1/hendrik/gevi") -> None: + + mouse_name = [ + 'M_Sert_Cre_41', + 'M_Sert_Cre_42', + 'M_Sert_Cre_45', + 'M_Sert_Cre_46', + 'M_Sert_Cre_49' + ] + + n_tris = [ + [15, 15, 30, 30, 30, 30,], # 0 in cond 7 + [15, 15, 30, 30, 30, 30,], # 0 in cond 7 + [15, 15, 30, 30, 30, 30,], # 0 in cond 7 + [20, 40, 20, 20,], # 0, 0, 0 in cond 5-7 + [20, 40, 20, 20,], # 0, 0, 0 in cond 5-7 + ] + + # 41, 42, 45, 46, 49: + # "1": "control", + # "2": "visual control grating 100 1s", + # "3": "optical Stimulation 20Hz 50% 5 Intervals", + # "4": "optical Stimulation 20Hz 100% 5 Intervals", + # "5": "optical Stimulation 20Hz 50% and grating 100", + # "6": "optical Stimulation 20Hz 100% and grating 100", + # "7": "grating 3s" + + lbs = [ + ['control', 'visual control', 'op20 50 5', 'op20 100 5', 'op20 50 grat', 'op20 100 grat'], + ['control', 'visual control', 'op20 50 5', 'op20 100 5', 'op20 50 grat', 'op20 100 grat'], + ['control', 'visual control', 'op20 50 5', 'op20 100 5', 'op20 50 grat', 'op20 100 grat'], + ['control', 'visual control', 'op20 50 5', 'op20 100 5'], + ['control', 'visual control', 'op20 50 5', 'op20 100 5'] + ] + + n_exp = len(n_tris[mouse]) + + for i_exp in range(n_exp): + n_tri = n_tris[mouse][i_exp] + for i_tri in range(n_tri): + + experiment_name: str = f"Exp{i_exp+1:03d}_Trial{i_tri+1:03d}" + tmp_fname = os.path.join( + fpath, "output_" + mouse_name[mouse], experiment_name + "_acceptor_donor.npz" + ) + print(f'Processing file "{tmp_fname}"...') + tmp = np.load(tmp_fname) + + tmp_data_sequence = tmp["data_donor"] + tmp_light_signal = tmp["data_acceptor"] + + if (i_exp == 0) and (i_tri == 0): + mask = tmp["mask"] + new_shape = [n_exp, *tmp_data_sequence.shape] + data_sequence = np.zeros(new_shape) + light_signal = np.zeros(new_shape) + data_sequence[i_exp] += tmp_data_sequence / n_tri + light_signal[i_exp] += tmp_light_signal / n_tri + + np.save("dsq_" + mouse_name[mouse], data_sequence) + np.save("lsq_" + mouse_name[mouse], light_signal) + np.save("msq_" + mouse_name[mouse], mask) + +if __name__ == "__main__": + argh.dispatch_command(loader) diff --git a/geci_plot.py b/geci_plot.py new file mode 100644 index 0000000..f7b45d6 --- /dev/null +++ b/geci_plot.py @@ -0,0 +1,209 @@ +import numpy as np +import matplotlib.pyplot as plt +import h5py # type: ignore +import argh +import scipy # type: ignore +import json +import os +from jsmin import jsmin # type:ignore +from functions.get_trials import get_trials + + +def func_pow(x, a, b, c): + return -a * x**b + c + + +def func_exp(x, a, b, c): + return a * np.exp(-x / b) + c + + +# mouse: int = 0, 1, 2, 3, 4 +def plot( + filename: str = "config_M_Sert_Cre_49.json", + experiment: int = 4, + skip_timesteps: int = 100, + # If there is no special ROI... Get one! This is just a backup + roi_control_path: str = "/data_1/hendrik/2023-03-15/ROI_control.mat", + roi_sdraken_path: str = "/data_1/hendrik/2023-03-15/ROI_sDarken.mat", + remove_fit: bool = True, + fit_power: bool = False, # True => -ax^b ; False => exp(-b) +) -> None: + + # lbs = [ + # [ + # "control", + # "visual control", + # "op20 50 5", + # "op20 100 5", + # "op20 50 grat", + # "op20 100 grat", + # ], + # [ + # "control", + # "visual control", + # "op20 50 5", + # "op20 100 5", + # "op20 50 grat", + # "op20 100 grat", + # ], + # [ + # "control", + # "visual control", + # "op20 50 5", + # "op20 100 5", + # "op20 50 grat", + # "op20 100 grat", + # ], + # ["control", "visual control", "op20 50 5", "op20 100 5"], + # ["control", "visual control", "op20 50 5", "op20 100 5"], + # ] + + if os.path.isfile(filename) is False: + print(f"{filename} is missing") + exit() + + with open(filename, "r") as file: + config = json.loads(jsmin(file.read())) + + raw_data_path: str = os.path.join( + config["basic_path"], + config["recoding_data"], + config["mouse_identifier"], + config["raw_path"], + ) + + if os.path.isdir(raw_data_path) is False: + print(f"ERROR: could not find raw directory {raw_data_path}!!!!") + exit() + + trials = get_trials(raw_data_path, experiment).numpy() + assert trials.shape[0] > 0 + + with open(os.path.join(raw_data_path, f"Exp{experiment:03d}_Trial{trials[0]:03d}_Part001_meta.txt"), "r") as file: + metadata = json.loads(jsmin(file.read())) + + experiment_names = metadata['sessionMetaData']['experimentNames'][str(experiment)] + + roi_path: str = os.path.join(config["basic_path"], config["recoding_data"]) + roi_control_mat: str = os.path.join(roi_path, "ROI_control.mat") + roi_sdarken_mat: str = os.path.join(roi_path, "ROI_sDarken.mat") + + if os.path.isdir(roi_control_mat): + roi_control_path = roi_control_mat + + if os.path.isdir(roi_sdarken_mat): + roi_sdraken_path = roi_sdarken_mat + + print("Load data...") + data = np.load("dsq_" + config["mouse_identifier"] + ".npy", mmap_mode="r") + + print("Load light signal...") + light = np.load("lsq_" + config["mouse_identifier"] + ".npy", mmap_mode="r") + + print("Load mask...") + mask = np.load("msq_" + config["mouse_identifier"] + ".npy") + + hf = h5py.File(roi_control_path, "r") + roi_lighten = np.array(hf["roi"]).T + roi_lighten *= mask + + hf = h5py.File(roi_sdraken_path, "r") + roi_darken = np.array(hf["roi"]).T + roi_darken *= mask + + plt.figure(1) + a_show = data[experiment - 1, :, :, 1000].copy() + a_show[(roi_darken + roi_lighten) < 0.5] = np.nan + plt.imshow(a_show) + 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_lighten) > 0.5] = np.nan + plt.imshow(a_dontshow) + plt.title(f"{config["mouse_identifier"]} -- Experiment: {experiment}") + plt.show(block=False) + + plt.figure(3) + light_exp = light[experiment - 1, :, :, skip_timesteps:].copy() + light_exp[(roi_darken + roi_lighten) < 0.5, :] = 0.0 + light_signal = light_exp.mean(axis=(0, 1)) + light_signal -= light_signal.min() + light_signal /= light_signal.max() + + a_exp = data[experiment - 1, :, :, skip_timesteps:].copy() + + if remove_fit: + combined_matrix = (roi_darken + roi_lighten) > 0 + idx = np.where(combined_matrix) + for idx_pos in range(0, idx[0].shape[0]): + temp = a_exp[idx[0][idx_pos], idx[1][idx_pos], :] + temp -= temp.mean() + + data_time = np.arange(0, temp.shape[0], dtype=np.float32) + skip_timesteps + data_time /= 100.0 + + data_min = temp.min() + data_max = temp.max() + data_delta = data_max - data_min + a_min = data_min - data_delta + b_min = 0.01 + a_max = data_max + data_delta + if fit_power: + b_max = 10.0 + else: + b_max = 100.0 + c_min = data_min - data_delta + c_max = data_max + data_delta + + try: + if fit_power: + popt, _ = scipy.optimize.curve_fit( + f=func_pow, + xdata=data_time, + ydata=np.nan_to_num(temp), + bounds=([a_min, b_min, c_min], [a_max, b_max, c_max]), + ) + pattern: np.ndarray | None = func_pow(data_time, *popt) + else: + popt, _ = scipy.optimize.curve_fit( + f=func_exp, + xdata=data_time, + ydata=np.nan_to_num(temp), + bounds=([a_min, b_min, c_min], [a_max, b_max, c_max]), + ) + pattern = func_exp(data_time, *popt) + + assert pattern is not None + pattern -= pattern.mean() + + scale = (temp * pattern).sum() / (pattern**2).sum() + pattern *= scale + + except ValueError: + print(f"Fit failed: Position ({idx[0][idx_pos]}, {idx[1][idx_pos]}") + pattern = None + + if pattern is not None: + temp -= pattern + + darken = a_exp[roi_darken > 0.5, :].sum(axis=0) / (roi_darken > 0.5).sum() + lighten = a_exp[roi_lighten > 0.5, :].sum(axis=0) / (roi_lighten > 0.5).sum() + + light_signal *= darken.max() - darken.min() + light_signal += darken.min() + + time_axis = np.arange(0, lighten.shape[-1], dtype=np.float32) + skip_timesteps + time_axis /= 100.0 + + 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.legend() + plt.show() + + +if __name__ == "__main__": + argh.dispatch_command(plot) diff --git a/stage_4_process.py b/stage_4_process.py index 909731f..f046ea9 100644 --- a/stage_4_process.py +++ b/stage_4_process.py @@ -99,7 +99,7 @@ def process_trial( 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(f"CUDA memory: {free_mem // 1024} MByte") mylogger.info(f"Data shape: {data.shape}") mylogger.info("-==- Done -==-") @@ -266,9 +266,9 @@ def process_trial( batch_size=config["alignment_batch_size"], fill_value=-100.0, ) - mylogger.info(f"Rotation: {round(float(angle_refref[0]),2)} degree") + mylogger.info(f"Rotation: {round(float(angle_refref[0]), 2)} degree") mylogger.info( - f"Translation: {round(float(tvec_refref[0]),1)} x {round(float(tvec_refref[1]),1)} pixel" + f"Translation: {round(float(tvec_refref[0]), 1)} x {round(float(tvec_refref[1]), 1)} pixel" ) if config["save_alignment"]: @@ -285,7 +285,7 @@ def process_trial( np.save(temp_path, tvec_refref.cpu()) mylogger.info("Moving & rotating the oxygenation ref image") - ref_image_oxygenation = tv.transforms.functional.affine( + ref_image_oxygenation = tv.transforms.functional.affine( # type: ignore img=ref_image_oxygenation.unsqueeze(0), angle=-float(angle_refref), translate=[0, 0], @@ -295,7 +295,7 @@ def process_trial( fill=-100.0, ) - ref_image_oxygenation = tv.transforms.functional.affine( + ref_image_oxygenation = tv.transforms.functional.affine( # type: ignore img=ref_image_oxygenation, angle=0, translate=[tvec_refref[1], tvec_refref[0]], @@ -313,8 +313,8 @@ def process_trial( volume_index: int = config["required_order"].index("volume") mylogger.info("Rotate acceptor") - data[acceptor_index, ...] = tv.transforms.functional.affine( - img=data[acceptor_index, ...], + data[acceptor_index, ...] = tv.transforms.functional.affine( # type: ignore + img=data[acceptor_index, ...], # type: ignore angle=-float(angle_refref), translate=[0, 0], scale=1.0, @@ -324,7 +324,7 @@ def process_trial( ) mylogger.info("Translate acceptor") - data[acceptor_index, ...] = tv.transforms.functional.affine( + data[acceptor_index, ...] = tv.transforms.functional.affine( # type: ignore img=data[acceptor_index, ...], angle=0, translate=[tvec_refref[1], tvec_refref[0]], @@ -335,7 +335,7 @@ def process_trial( ) mylogger.info("Rotate oxygenation") - data[oxygenation_index, ...] = tv.transforms.functional.affine( + data[oxygenation_index, ...] = tv.transforms.functional.affine( # type: ignore img=data[oxygenation_index, ...], angle=-float(angle_refref), translate=[0, 0], @@ -346,7 +346,7 @@ def process_trial( ) mylogger.info("Translate oxygenation") - data[oxygenation_index, ...] = tv.transforms.functional.affine( + data[oxygenation_index, ...] = tv.transforms.functional.affine( # type: ignore img=data[oxygenation_index, ...], angle=0, translate=[tvec_refref[1], tvec_refref[0]], @@ -359,7 +359,7 @@ def process_trial( mylogger.info("Perform rotation between donor and volume and its ref images") mylogger.info("for all frames and then rotate all the data accordingly") - perform_donor_volume_rotation + ( data[acceptor_index, ...], data[donor_index, ...], @@ -381,9 +381,9 @@ def process_trial( mylogger.info( f"angles: " - f"min {round(float(angle_donor_volume.min()),2)} " - f"max {round(float(angle_donor_volume.max()),2)} " - f"mean {round(float(angle_donor_volume.mean()),2)} " + f"min {round(float(angle_donor_volume.min()), 2)} " + f"max {round(float(angle_donor_volume.max()), 2)} " + f"mean {round(float(angle_donor_volume.mean()), 2)} " ) if config["save_alignment"]: @@ -417,15 +417,15 @@ def process_trial( mylogger.info( f"translation dim 0: " - f"min {round(float(tvec_donor_volume[:,0].min()),1)} " - f"max {round(float(tvec_donor_volume[:,0].max()),1)} " - f"mean {round(float(tvec_donor_volume[:,0].mean()),1)} " + f"min {round(float(tvec_donor_volume[:, 0].min()), 1)} " + f"max {round(float(tvec_donor_volume[:, 0].max()), 1)} " + f"mean {round(float(tvec_donor_volume[:, 0].mean()), 1)} " ) mylogger.info( f"translation dim 1: " - f"min {round(float(tvec_donor_volume[:,1].min()),1)} " - f"max {round(float(tvec_donor_volume[:,1].max()),1)} " - f"mean {round(float(tvec_donor_volume[:,1].mean()),1)} " + f"min {round(float(tvec_donor_volume[:, 1].min()), 1)} " + f"max {round(float(tvec_donor_volume[:, 1].max()), 1)} " + f"mean {round(float(tvec_donor_volume[:, 1].mean()), 1)} " ) if config["save_alignment"]: @@ -471,172 +471,183 @@ def process_trial( sample_frequency: float = 1.0 / meta_frame_time - mylogger.info("Extract heartbeat from volume signal") - heartbeat_ts: torch.Tensor = bandpass( - data=data[volume_index, ...].movedim(0, -1).clone(), - low_frequency=config["lower_freqency_bandpass"], - high_frequency=config["upper_freqency_bandpass"], - fs=sample_frequency, - filtfilt_chuck_size=config["heartbeat_filtfilt_chuck_size"], - ) - heartbeat_ts = heartbeat_ts.flatten(start_dim=0, end_dim=-2) - mask_flatten: torch.Tensor = mask_positve.flatten(start_dim=0, end_dim=-1) + if config["gevi"]: + assert config["heartbeat_remove"] - heartbeat_ts = heartbeat_ts[mask_flatten, :] - - heartbeat_ts = heartbeat_ts.movedim(0, -1) - heartbeat_ts -= heartbeat_ts.mean(dim=0, keepdim=True) - - try: - volume_heartbeat, _, _ = torch.linalg.svd(heartbeat_ts, full_matrices=False) - except torch.cuda.OutOfMemoryError: - mylogger.info("torch.cuda.OutOfMemoryError: Fallback to cpu") - volume_heartbeat_cpu, _, _ = torch.linalg.svd( - heartbeat_ts.cpu(), full_matrices=False - ) - volume_heartbeat = volume_heartbeat_cpu.to(heartbeat_ts.data, copy=True) - del volume_heartbeat_cpu - - volume_heartbeat = volume_heartbeat[:, 0] - volume_heartbeat -= volume_heartbeat[ - config["skip_frames_in_the_beginning"] : - ].mean() - - del heartbeat_ts - - 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 config["save_heartbeat"]: - temp_path = os.path.join( - config["export_path"], experiment_name + "_volume_heartbeat.npy" - ) - mylogger.info(f"Save volume heartbeat to {temp_path}") - np.save(temp_path, volume_heartbeat.cpu()) - mylogger.info("-==- Done -==-") - - volume_heartbeat = volume_heartbeat.unsqueeze(0).unsqueeze(0) - norm_volume_heartbeat = ( - volume_heartbeat[..., config["skip_frames_in_the_beginning"] :] ** 2 - ).sum(dim=-1) - - heartbeat_coefficients: torch.Tensor = torch.zeros( - (data.shape[0], data.shape[-2], data.shape[-1]), - dtype=data.dtype, - device=data.device, - ) - for i in range(0, data.shape[0]): - y = bandpass( - data=data[i, ...].movedim(0, -1).clone(), + if config["heartbeat_remove"]: + mylogger.info("Extract heartbeat from volume signal") + heartbeat_ts: torch.Tensor = bandpass( + data=data[volume_index, ...].movedim(0, -1).clone(), low_frequency=config["lower_freqency_bandpass"], high_frequency=config["upper_freqency_bandpass"], fs=sample_frequency, filtfilt_chuck_size=config["heartbeat_filtfilt_chuck_size"], - )[..., config["skip_frames_in_the_beginning"] :] - y -= y.mean(dim=-1, keepdim=True) - - heartbeat_coefficients[i, ...] = ( - volume_heartbeat[..., config["skip_frames_in_the_beginning"] :] * y - ).sum(dim=-1) / norm_volume_heartbeat - - heartbeat_coefficients[i, ...] *= mask_positve.type( - dtype=heartbeat_coefficients.dtype ) - del y + heartbeat_ts = heartbeat_ts.flatten(start_dim=0, end_dim=-2) + mask_flatten: torch.Tensor = mask_positve.flatten(start_dim=0, end_dim=-1) - if config["save_heartbeat"]: - temp_path = os.path.join( - config["export_path"], experiment_name + "_heartbeat_coefficients.npy" + heartbeat_ts = heartbeat_ts[mask_flatten, :] + + heartbeat_ts = heartbeat_ts.movedim(0, -1) + heartbeat_ts -= heartbeat_ts.mean(dim=0, keepdim=True) + + try: + volume_heartbeat, _, _ = torch.linalg.svd(heartbeat_ts, full_matrices=False) + except torch.cuda.OutOfMemoryError: + mylogger.info("torch.cuda.OutOfMemoryError: Fallback to cpu") + volume_heartbeat_cpu, _, _ = torch.linalg.svd( + heartbeat_ts.cpu(), full_matrices=False + ) + volume_heartbeat = volume_heartbeat_cpu.to(heartbeat_ts.data, copy=True) + del volume_heartbeat_cpu + + volume_heartbeat = volume_heartbeat[:, 0] + volume_heartbeat -= volume_heartbeat[ + config["skip_frames_in_the_beginning"] : + ].mean() + + del heartbeat_ts + + 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 config["save_heartbeat"]: + temp_path = os.path.join( + config["export_path"], experiment_name + "_volume_heartbeat.npy" + ) + mylogger.info(f"Save volume heartbeat to {temp_path}") + np.save(temp_path, volume_heartbeat.cpu()) + mylogger.info("-==- Done -==-") + + volume_heartbeat = volume_heartbeat.unsqueeze(0).unsqueeze(0) + norm_volume_heartbeat = ( + volume_heartbeat[..., config["skip_frames_in_the_beginning"] :] ** 2 + ).sum(dim=-1) + + heartbeat_coefficients: torch.Tensor = torch.zeros( + (data.shape[0], data.shape[-2], data.shape[-1]), + dtype=data.dtype, + device=data.device, ) - mylogger.info(f"Save heartbeat coefficients to {temp_path}") - np.save(temp_path, heartbeat_coefficients.cpu()) - mylogger.info("-==- Done -==-") + for i in range(0, data.shape[0]): + y = bandpass( + data=data[i, ...].movedim(0, -1).clone(), + low_frequency=config["lower_freqency_bandpass"], + high_frequency=config["upper_freqency_bandpass"], + fs=sample_frequency, + filtfilt_chuck_size=config["heartbeat_filtfilt_chuck_size"], + )[..., config["skip_frames_in_the_beginning"] :] + y -= y.mean(dim=-1, keepdim=True) - mylogger.info("Remove heart beat from data") - data -= heartbeat_coefficients.unsqueeze(1) * volume_heartbeat.unsqueeze(0).movedim( - -1, 1 - ) - mylogger.info("-==- Done -==-") + heartbeat_coefficients[i, ...] = ( + volume_heartbeat[..., config["skip_frames_in_the_beginning"] :] * y + ).sum(dim=-1) / norm_volume_heartbeat - donor_heartbeat_factor = heartbeat_coefficients[donor_index, ...].clone() - acceptor_heartbeat_factor = heartbeat_coefficients[acceptor_index, ...].clone() - del heartbeat_coefficients + heartbeat_coefficients[i, ...] *= mask_positve.type( + dtype=heartbeat_coefficients.dtype + ) + del y - 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)] + if config["save_heartbeat"]: + temp_path = os.path.join( + config["export_path"], experiment_name + "_heartbeat_coefficients.npy" + ) + mylogger.info(f"Save heartbeat coefficients to {temp_path}") + np.save(temp_path, heartbeat_coefficients.cpu()) + mylogger.info("-==- Done -==-") + + mylogger.info("Remove heart beat from data") + data -= heartbeat_coefficients.unsqueeze(1) * volume_heartbeat.unsqueeze( + 0 + ).movedim(-1, 1) + mylogger.info("-==- Done -==-") + + 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") + + 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 + + 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 -==-") + + 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("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("Divide by mean over time") + data /= data[:, config["skip_frames_in_the_beginning"] :, ...].nanmean( + dim=1, + keepdim=True, ) - 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("-==- Done -==-") - 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()) - - 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("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("Divide by mean over time") - data /= data[:, config["skip_frames_in_the_beginning"] :, ...].nanmean( - dim=1, - keepdim=True, - ) data = data.nan_to_num(nan=0.0) - mylogger.info("-==- Done -==-") - mylogger.info("Preparation for regression -- Gauss smear") spatial_width = float(config["gauss_smear_spatial_width"]) @@ -669,7 +680,7 @@ def process_trial( 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(f"CUDA memory: {free_mem // 1024} MByte") overwrite_fft_gauss: None | torch.Tensor = None for i in range(0, data_filtered.shape[0]): @@ -703,7 +714,7 @@ def process_trial( 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(f"CUDA memory: {free_mem // 1024} MByte") mylogger.info("-==- Done -==-") mylogger.info("Preperation for Regression") @@ -747,6 +758,9 @@ def process_trial( mylogger.info("-==- Done -==-") else: dual_signal_mode = False + target_id = config["required_order"].index("acceptor") + data_acceptor = data[target_id, ...].clone() + data_acceptor[mask_negative, :] = 0.0 if len(config["target_camera_donor"]) > 0: mylogger.info("Regression Donor") @@ -781,6 +795,9 @@ def process_trial( mylogger.info("-==- Done -==-") else: dual_signal_mode = False + target_id = config["required_order"].index("donor") + data_donor = data[target_id, ...].clone() + data_donor[mask_negative, :] = 0.0 del data del data_filtered @@ -791,24 +808,119 @@ def process_trial( 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(f"CUDA memory: {free_mem // 1024} MByte") + + # ##################### + + if config["gevi"]: + assert dual_signal_mode + else: + assert dual_signal_mode is False + + if dual_signal_mode is False: + + mylogger.info("mono signal model") + + mylogger.info("Remove nan") + data_acceptor = torch.nan_to_num(data_acceptor, nan=0.0) + data_donor = torch.nan_to_num(data_donor, nan=0.0) + mylogger.info("-==- Done -==-") + + if config["binning_enable"] and config["binning_at_the_end"]: + mylogger.info("Binning of data") + mylogger.info( + ( + f"kernel_size={int(config['binning_kernel_size'])}, " + f"stride={int(config['binning_stride'])}, " + "divisor_override=None" + ) + ) + + data_acceptor = binning( + data_acceptor.unsqueeze(-1), + kernel_size=int(config["binning_kernel_size"]), + stride=int(config["binning_stride"]), + divisor_override=None, + ).squeeze(-1) + + data_donor = binning( + data_donor.unsqueeze(-1), + kernel_size=int(config["binning_kernel_size"]), + stride=int(config["binning_stride"]), + divisor_override=None, + ).squeeze(-1) + + mask_positve = ( + binning( + mask_positve.unsqueeze(-1).unsqueeze(-1).type(dtype=dtype), + kernel_size=int(config["binning_kernel_size"]), + stride=int(config["binning_stride"]), + divisor_override=None, + ) + .squeeze(-1) + .squeeze(-1) + ) + mask_positve = (mask_positve > 0).type(torch.bool) + + if config["save_as_python"]: + + temp_path = os.path.join( + config["export_path"], experiment_name + "_acceptor_donor.npz" + ) + mylogger.info(f"Save data donor and acceptor and mask to {temp_path}") + np.savez_compressed( + temp_path, + data_acceptor=data_acceptor.cpu(), + data_donor=data_donor.cpu(), + mask=mask_positve.cpu(), + ) + + if config["save_as_matlab"]: + temp_path = os.path.join( + config["export_path"], experiment_name + "_acceptor_donor.hd5" + ) + mylogger.info(f"Save data donor and acceptor and mask to {temp_path}") + file_handle = h5py.File(temp_path, "w") + + mask_positve = mask_positve.movedim(0, -1) + data_acceptor = data_acceptor.movedim(1, -1).movedim(0, -1) + data_donor = data_donor.movedim(1, -1).movedim(0, -1) + _ = file_handle.create_dataset( + "mask", + data=mask_positve.type(torch.uint8).cpu(), + compression="gzip", + compression_opts=9, + ) + _ = file_handle.create_dataset( + "data_acceptor", + data=data_acceptor.cpu(), + compression="gzip", + compression_opts=9, + ) + _ = file_handle.create_dataset( + "data_donor", + data=data_donor.cpu(), + compression="gzip", + compression_opts=9, + ) + mylogger.info("Reminder: How to read with matlab:") + mylogger.info(f"mask = h5read('{temp_path}','/mask');") + mylogger.info(f"data_acceptor = h5read('{temp_path}','/data_acceptor');") + mylogger.info(f"data_donor = h5read('{temp_path}','/data_donor');") + file_handle.close() + return + # ##################### mylogger.info("Calculate ratio sequence") - if dual_signal_mode: - if config["classical_ratio_mode"]: - mylogger.info("via acceptor / donor") - ratio_sequence: torch.Tensor = data_acceptor / data_donor - mylogger.info("via / mean over time") - ratio_sequence /= ratio_sequence.mean(dim=-1, keepdim=True) - else: - mylogger.info("via 1.0 + acceptor - donor") - ratio_sequence = 1.0 + data_acceptor - data_donor + + if config["classical_ratio_mode"]: + mylogger.info("via acceptor / donor") + ratio_sequence: torch.Tensor = data_acceptor / data_donor + mylogger.info("via / mean over time") + ratio_sequence /= ratio_sequence.mean(dim=-1, keepdim=True) else: - mylogger.info("mono signal model") - if len(config["target_camera_donor"]) > 0: - ratio_sequence = data_donor.clone() - else: - ratio_sequence = data_acceptor.clone() + mylogger.info("via 1.0 + acceptor - donor") + ratio_sequence = 1.0 + data_acceptor - data_donor mylogger.info("Remove nan") ratio_sequence = torch.nan_to_num(ratio_sequence, nan=0.0)