diff --git a/doit.sh b/doit.sh new file mode 100644 index 0000000..999e2bc --- /dev/null +++ b/doit.sh @@ -0,0 +1,34 @@ +config="config_M_Sert_Cre_49.json" + +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter.py -p -e 1 -c $config +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter_svd.py -p -e 1 -c $config + +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter.py -p -e 2 -c $config +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter_svd.py -p -e 2 -c $config + +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter.py -p -e 3 -c $config +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter_svd.py -p -e 3 -c $config + +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter.py -p -e 4 -c $config +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter_svd.py -p -e 4 -c $config + +config="config_M_Sert_Cre_41.json" + +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter.py -p -e 1 -c $config +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter_svd.py -p -e 1 -c $config + +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter.py -p -e 2 -c $config +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter_svd.py -p -e 2 -c $config + +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter.py -p -e 3 -c $config +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter_svd.py -p -e 3 -c $config + +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter.py -p -e 4 -c $config +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter_svd.py -p -e 4 -c $config + +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter.py -p -e 5 -c $config +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter_svd.py -p -e 5 -c $config + +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter.py -p -e 6 -c $config +/data_1/davrot/P3.11/bin/python3 olivia_data_plotter_svd.py -p -e 6 -c $config + diff --git a/olivia_data_plotter.py b/olivia_data_plotter.py index b00dd06..0acf2eb 100644 --- a/olivia_data_plotter.py +++ b/olivia_data_plotter.py @@ -6,16 +6,27 @@ from functions.load_config import load_config from functions.get_trials import get_trials import h5py # type: ignore - +import torch +import scipy # type: ignore import argh +from functions.data_raw_loader import data_raw_loader -def main(*, experiment_id: int = 1, config_filename: str = "config.json") -> None: +def main( + *, + experiment_id: int = 4, + config_filename: str = "config.json", + highpass_freqency: float = 0.5, + lowpass_freqency: float = 10.0, + butter_worth_order: int = 4, + log_stage_name: str = "olivia", + plot_show: bool = True, +) -> None: mylogger = create_logger( save_logging_messages=True, display_logging_messages=True, - log_stage_name="test_xxx", + log_stage_name=log_stage_name, ) config = load_config(mylogger=mylogger, filename=config_filename) @@ -61,13 +72,142 @@ def main(*, experiment_id: int = 1, config_filename: str = "config.json") -> Non control = ratio_sequence[control_roi, :].mean(axis=0) s_darken = ratio_sequence[s_darken_roi, :].mean(axis=0) + max_value = max( + [ + control[config["skip_frames_in_the_beginning"] :].max(), + s_darken[config["skip_frames_in_the_beginning"] :].max(), + ] + ) + min_value = min( + [ + control[config["skip_frames_in_the_beginning"] :].min(), + s_darken[config["skip_frames_in_the_beginning"] :].min(), + ] + ) + + first_trial_id: int = int(get_trials(raw_data_path, experiment_id).min()) + ( + meta_channels, + meta_mouse_markings, + meta_recording_date, + meta_stimulation_times, + meta_experiment_names, + meta_trial_recording_duration, + meta_frame_time, + meta_mouse, + data, + ) = data_raw_loader( + raw_data_path=raw_data_path, + mylogger=mylogger, + experiment_id=experiment_id, + trial_id=first_trial_id, + device=torch.device("cpu"), + force_to_cpu_memory=True, + config=config, + ) + + idx = config["required_order"].index("acceptor") + acceptor = data[..., idx].mean(axis=0).mean(axis=0) + acceptor -= acceptor[config["skip_frames_in_the_beginning"] :].min() + acceptor /= acceptor[config["skip_frames_in_the_beginning"] :].max() + + acceptor_f0 = acceptor.clone() + acceptor_f0 *= max_value - min_value + acceptor_f0 += min_value + + b, a = scipy.signal.butter( + butter_worth_order, + lowpass_freqency, + btype="low", + output="ba", + fs=1.0 / meta_frame_time, + ) + control_f1 = scipy.signal.filtfilt(b, a, control) + s_darken_f1 = scipy.signal.filtfilt(b, a, s_darken) + + b, a = scipy.signal.butter( + butter_worth_order, + highpass_freqency, + btype="high", + output="ba", + fs=1.0 / meta_frame_time, + ) + control_f1 = scipy.signal.filtfilt(b, a, control_f1) + s_darken_f1 = scipy.signal.filtfilt(b, a, s_darken_f1) + + max_value = max( + [ + control_f1[config["skip_frames_in_the_beginning"] :].max(), + s_darken_f1[config["skip_frames_in_the_beginning"] :].max(), + ] + ) + min_value = min( + [ + control_f1[config["skip_frames_in_the_beginning"] :].min(), + s_darken_f1[config["skip_frames_in_the_beginning"] :].min(), + ] + ) + + acceptor_f1 = acceptor.clone() + acceptor_f1 *= max_value - min_value + acceptor_f1 += min_value + t = np.arange(0, control.shape[0]) / 100.0 - plt.plot(t, control, label="control") - plt.plot(t, s_darken, label="sDarken") + plt.figure(figsize=(10, 10)) + plt.subplot(2, 1, 1) + plt.plot( + t[config["skip_frames_in_the_beginning"] :], + acceptor_f0[config["skip_frames_in_the_beginning"] :], + color=(0.5, 0.5, 0.5), + label="light (acceptor)", + ) + + plt.plot( + t[config["skip_frames_in_the_beginning"] :], + control[config["skip_frames_in_the_beginning"] :], + label="control", + ) + plt.plot( + t[config["skip_frames_in_the_beginning"] :], + s_darken[config["skip_frames_in_the_beginning"] :], + label="sDarken", + ) + plt.title( + f"Experiment {experiment_id} {config['recoding_data']} {config['mouse_identifier']}" + ) + plt.legend() plt.xlabel("Time [sec]") - plt.show() + + plt.subplot(2, 1, 2) + + plt.plot( + t[config["skip_frames_in_the_beginning"] :], + acceptor_f1[config["skip_frames_in_the_beginning"] :], + color=(0.5, 0.5, 0.5), + label="light (acceptor)", + ) + + plt.plot( + t[config["skip_frames_in_the_beginning"] :], + control_f1[config["skip_frames_in_the_beginning"] :], + label=f"control ({highpass_freqency}Hz - {lowpass_freqency}Hz)", + ) + plt.plot( + t[config["skip_frames_in_the_beginning"] :], + s_darken_f1[config["skip_frames_in_the_beginning"] :], + label=f"sDarken ({highpass_freqency}Hz - {lowpass_freqency}Hz)", + ) + + plt.legend() + plt.xlabel("Time [sec]") + plt.savefig( + f"olivia_both_Exp{experiment_id}_{config['recoding_data']}_{config['mouse_identifier']}.png", + dpi=300, + ) + if plot_show: + plt.show() if __name__ == "__main__": diff --git a/olivia_data_plotter_svd.py b/olivia_data_plotter_svd.py index 207193a..956aaa9 100644 --- a/olivia_data_plotter_svd.py +++ b/olivia_data_plotter_svd.py @@ -7,16 +7,73 @@ from functions.load_config import load_config from functions.get_trials import get_trials import h5py # type: ignore import torch +import scipy # type: ignore import argh +from functions.data_raw_loader import data_raw_loader + +# def func(x, a, b, c, dt): +# return a * (x - dt) ** b + c -def main(*, experiment_id: int = 1, config_filename: str = "config.json") -> None: +# def fit( +# ratio_sequence: np.ndarray, +# t: np.ndarray, +# config: dict, +# ) -> tuple[np.ndarray, np.ndarray]: + +# data_min = ratio_sequence[config["skip_frames_in_the_beginning"] :].min() +# data_max = ratio_sequence[config["skip_frames_in_the_beginning"] :].max() +# b_min = 1.0 +# b_max = 3.0 + +# temp_1 = max([abs(data_min), abs(data_max)]) +# a_min = -temp_1 - 2 * abs(data_max - data_min) +# a_max = +temp_1 + 2 * abs(data_max - data_min) + +# try: +# popt, _ = scipy.optimize.curve_fit( +# f=func, +# xdata=t[config["skip_frames_in_the_beginning"] :], +# ydata=np.nan_to_num( +# ratio_sequence[config["skip_frames_in_the_beginning"] :] +# ), +# bounds=([a_min, b_min, a_min, -t[-1]], [a_max, b_max, a_max, t[-1]]), +# ) +# a: float | None = float(popt[0]) +# b: float | None = float(popt[1]) +# c: float | None = float(popt[2]) +# dt: float | None = float(popt[3]) + +# except ValueError: +# a = None +# b = None +# c = None +# dt = None + +# print(a, b, c, dt) + +# f1 = func(t, a, b, c, dt) +# ratio_sequence_f1 = ratio_sequence - f1 +# return ratio_sequence_f1, f1 + + +def main( + *, + experiment_id: int = 4, + config_filename: str = "config.json", + highpass_freqency: float = 0.5, + lowpass_freqency: float = 10.0, + butter_worth_order: int = 4, + log_stage_name: str = "olivia_svd", + scale_before_substraction: bool = True, + plot_show: bool = True, +) -> None: mylogger = create_logger( save_logging_messages=True, display_logging_messages=True, - log_stage_name="test_xxx", + log_stage_name=log_stage_name, ) config = load_config(mylogger=mylogger, filename=config_filename) @@ -58,10 +115,18 @@ def main(*, experiment_id: int = 1, config_filename: str = "config.json") -> Non rs_s_core, _, _ = torch.linalg.svd(torch.tensor(rs_s.T), full_matrices=False) rs_s_core = rs_s_core[:, 0].numpy() - rs_s_core -= rs_s_core.mean(keepdims=True) - rs_c_core -= rs_c_core.mean(keepdims=True) + rs_s_core -= rs_s_core[config["skip_frames_in_the_beginning"] :].mean( + keepdims=True + ) + rs_c_core -= rs_c_core[config["skip_frames_in_the_beginning"] :].mean( + keepdims=True + ) - rs_c_core *= (rs_s_core * rs_c_core).sum() / (rs_c_core**2).sum() + if scale_before_substraction: + rs_c_core *= ( + rs_s_core[config["skip_frames_in_the_beginning"] :] + * rs_c_core[config["skip_frames_in_the_beginning"] :] + ).sum() / (rs_c_core[config["skip_frames_in_the_beginning"] :] ** 2).sum() if i == 0: ratio_sequence = rs_s_core - rs_c_core @@ -72,10 +137,105 @@ def main(*, experiment_id: int = 1, config_filename: str = "config.json") -> Non t = np.arange(0, ratio_sequence.shape[0]) / 100.0 - plt.plot(t, ratio_sequence, label="sDarken - control") + # ratio_sequence_f1, f1 = fit( + # ratio_sequence=ratio_sequence, + # t=t, + # config=config, + # ) + + first_trial_id: int = int(get_trials(raw_data_path, experiment_id).min()) + ( + meta_channels, + meta_mouse_markings, + meta_recording_date, + meta_stimulation_times, + meta_experiment_names, + meta_trial_recording_duration, + meta_frame_time, + meta_mouse, + data, + ) = data_raw_loader( + raw_data_path=raw_data_path, + mylogger=mylogger, + experiment_id=experiment_id, + trial_id=first_trial_id, + device=torch.device("cpu"), + force_to_cpu_memory=True, + config=config, + ) + + b, a = scipy.signal.butter( + butter_worth_order, + lowpass_freqency, + btype="low", + output="ba", + fs=1.0 / meta_frame_time, + ) + ratio_sequence_f1 = scipy.signal.filtfilt(b, a, ratio_sequence) + + b, a = scipy.signal.butter( + butter_worth_order, + highpass_freqency, + btype="high", + output="ba", + fs=1.0 / meta_frame_time, + ) + ratio_sequence_f2 = scipy.signal.filtfilt(b, a, ratio_sequence_f1) + + idx = config["required_order"].index("acceptor") + acceptor = data[..., idx].mean(axis=0).mean(axis=0) + acceptor -= acceptor[config["skip_frames_in_the_beginning"] :].min() + acceptor /= acceptor[config["skip_frames_in_the_beginning"] :].max() + + acceptor *= ( + ratio_sequence_f2[config["skip_frames_in_the_beginning"] :].max() + - ratio_sequence_f2[config["skip_frames_in_the_beginning"] :].min() + ) + acceptor += ratio_sequence_f2[config["skip_frames_in_the_beginning"] :].min() + + plt.figure(figsize=(10, 10)) + plt.subplot(2, 1, 1) + plt.plot( + t[config["skip_frames_in_the_beginning"] :], + ratio_sequence[config["skip_frames_in_the_beginning"] :], + label=f"sDarken - control (scaled:{scale_before_substraction})", + ) + plt.plot( + t[config["skip_frames_in_the_beginning"] :], + ratio_sequence_f1[config["skip_frames_in_the_beginning"] :], + label=f"low pass {lowpass_freqency} Hz", + ) + plt.plot( + t[config["skip_frames_in_the_beginning"] :], + ratio_sequence_f2[config["skip_frames_in_the_beginning"] :], + label=f"high pass {highpass_freqency} Hz", + ) + plt.xlabel("Time [sec]") + plt.title( + f"Experiment {experiment_id} {config['recoding_data']} {config['mouse_identifier']}" + ) + plt.legend() + plt.subplot(2, 1, 2) + plt.plot( + t[config["skip_frames_in_the_beginning"] :], + acceptor[config["skip_frames_in_the_beginning"] :], + color=(0.5, 0.5, 0.5), + label="light (acceptor)", + ) + plt.plot( + t[config["skip_frames_in_the_beginning"] :], + ratio_sequence_f2[config["skip_frames_in_the_beginning"] :], + label=f"high pass {highpass_freqency} Hz", + ) plt.legend() plt.xlabel("Time [sec]") - plt.show() + + plt.savefig( + f"olivia_Exp{experiment_id}_{config['recoding_data']}_{config['mouse_identifier']}.png", + dpi=300, + ) + if plot_show: + plt.show() if __name__ == "__main__":