diff --git a/geci_loader.py b/geci_loader.py index 5ee5884..8d530f5 100644 --- a/geci_loader.py +++ b/geci_loader.py @@ -5,10 +5,26 @@ from jsmin import jsmin # type:ignore import argh from functions.get_trials import get_trials from functions.get_experiments import get_experiments +import scipy # type: ignore + + +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 def loader( - filename: str = "config_M_Sert_Cre_49.json", fpath: str = "/data_1/hendrik/gevi" + filename: str = "config_M_Sert_Cre_49.json", + fpath: str = "/data_1/hendrik/gevi", + 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 = True, + fit_power: bool = False, # True => -ax^b ; False => exp(-b) ) -> None: if os.path.isfile(filename) is False: @@ -25,9 +41,26 @@ def loader( config["raw_path"], ) + 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" + + 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 + + roi_control: np.ndarray = np.load(roi_control_path) + roi_darken: np.ndarray = np.load(roi_sdarken_path) + experiments = get_experiments(raw_data_path).numpy() n_exp = experiments.shape[0] + first_run: bool = True + for i_exp in range(0, n_exp): trials = get_trials(raw_data_path, experiments[i_exp]).numpy() n_tri = trials.shape[0] @@ -46,18 +79,83 @@ def loader( tmp = np.load(tmp_fname) tmp_data_sequence = tmp["data_donor"] + tmp_data_sequence = tmp_data_sequence[:, :, skip_timesteps:] tmp_light_signal = tmp["data_acceptor"] + tmp_light_signal = tmp_light_signal[:, :, skip_timesteps:] - if (i_exp == 0) and (i_tri == 0): + if first_run: mask = tmp["mask"] new_shape = [n_exp, *tmp_data_sequence.shape] data_sequence = np.zeros(new_shape) light_signal = np.zeros(new_shape) + first_run = False - # Here you might want to use the exp fit and removal... - data_sequence[i_exp] += tmp_data_sequence / n_tri - light_signal[i_exp] += tmp_light_signal / n_tri + if remove_fit: + roi_control *= mask + assert roi_control.sum() > 0, "ROI control empty" + roi_darken *= mask + assert roi_darken.sum() > 0, "ROI sDarken empty" + if remove_fit: + combined_matrix = (roi_darken + roi_control) > 0 + idx = np.where(combined_matrix) + for idx_pos in range(0, idx[0].shape[0]): + + temp = tmp_data_sequence[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 + tmp_data_sequence[idx[0][idx_pos], idx[1][idx_pos], :] = temp + + data_sequence[i_exp] += tmp_data_sequence + 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) diff --git a/geci_plot.py b/geci_plot.py index 8fd9ab4..e12ac17 100644 --- a/geci_plot.py +++ b/geci_plot.py @@ -19,11 +19,11 @@ def func_exp(x, a, b, c): def plot( filename: str = "config_M_Sert_Cre_49.json", experiment: int = 4, - skip_timesteps: int = 100, + skip_timesteps: int = 0, # 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 = True, + remove_fit: bool = False, fit_power: bool = False, # True => -ax^b ; False => exp(-b) ) -> None: @@ -154,6 +154,7 @@ def plot( if pattern is not None: temp -= pattern + a_exp[idx[0][idx_pos], idx[1][idx_pos], :] = temp darken = a_exp[roi_darken > 0.5, :].sum(axis=0) / (roi_darken > 0.5).sum() lighten = a_exp[roi_control > 0.5, :].sum(axis=0) / (roi_control > 0.5).sum()