From db2723f4084c7c2ae23507f530722843912f7e4a Mon Sep 17 00:00:00 2001 From: Achille-V <72287888+Achille-V@users.noreply.github.com> Date: Sat, 18 Jul 2026 21:06:34 -0400 Subject: [PATCH 1/2] Adding pari-combo methods --- dFC_method_pairs.txt | 6 + pydfc/dfc_methods/__init__.py | 14 ++ pydfc/dfc_methods/adaptive_dcc_combo.py | 145 ++++++++++++++++ .../dfc_methods/adaptive_multiscale_combo.py | 141 ++++++++++++++++ .../adaptive_random_fourier_combo.py | 140 ++++++++++++++++ pydfc/dfc_methods/dcc_multiscale_combo.py | 144 ++++++++++++++++ pydfc/dfc_methods/dcc_time_freq_combo.py | 157 ++++++++++++++++++ pydfc/dfc_methods/random_fourier_dcc_combo.py | 145 ++++++++++++++++ .../random_fourier_time_freq_combo.py | 124 ++++++++++++++ pydfc/multi_analysis.py | 28 ++++ tests/test_validation/dfc_method_wrappers.py | 7 + 11 files changed, 1051 insertions(+) create mode 100644 dFC_method_pairs.txt create mode 100644 pydfc/dfc_methods/adaptive_dcc_combo.py create mode 100644 pydfc/dfc_methods/adaptive_multiscale_combo.py create mode 100644 pydfc/dfc_methods/adaptive_random_fourier_combo.py create mode 100644 pydfc/dfc_methods/dcc_multiscale_combo.py create mode 100644 pydfc/dfc_methods/dcc_time_freq_combo.py create mode 100644 pydfc/dfc_methods/random_fourier_dcc_combo.py create mode 100644 pydfc/dfc_methods/random_fourier_time_freq_combo.py diff --git a/dFC_method_pairs.txt b/dFC_method_pairs.txt new file mode 100644 index 0000000..7ff4f6c --- /dev/null +++ b/dFC_method_pairs.txt @@ -0,0 +1,6 @@ +AdaptiveExponentialWindow and DCCConnectivity +AdaptiveExponentialWindow and MultiscaleWindow +RandomFourierDependence and DCCConnectivity +RandomFourierDependence and Time-Freq +DCCConnectivity and MultiscaleWindow +DCCConnectivity and Time-Freq diff --git a/pydfc/dfc_methods/__init__.py b/pydfc/dfc_methods/__init__.py index 30b88ba..3145f6e 100644 --- a/pydfc/dfc_methods/__init__.py +++ b/pydfc/dfc_methods/__init__.py @@ -1,6 +1,9 @@ """The :mod:`pydfc.dfc_methods` contains dFC methods objects.""" from .adaptive_exponential_window import ADAPTIVE_EXPONENTIAL_WINDOW +from .adaptive_dcc_combo import ADAPTIVE_DCC_COMBO +from .adaptive_multiscale_combo import ADAPTIVE_MULTISCALE_COMBO +from .adaptive_random_fourier_combo import ADAPTIVE_RANDOM_FOURIER_COMBO from .agglomerative_states import AGGLOMERATIVE_STATES from .amplitude_envelope_correlation import AMPLITUDE_ENVELOPE_CORRELATION from .base_dfc_method import BaseDFCMethod @@ -12,6 +15,8 @@ from .copula_tail_dependence import COPULA_TAIL_DEPENDENCE from .curvature_correlation import CURVATURE_CORRELATION from .dcc_connectivity import DCC_CONNECTIVITY +from .dcc_multiscale_combo import DCC_MULTISCALE_COMBO +from .dcc_time_freq_combo import DCC_TIME_FREQ_COMBO from .derivative_weighted_window import DERIVATIVE_WEIGHTED_WINDOW from .differential_coactivation import DIFFERENTIAL_COACTIVATION from .discrete_hmm import HMM_DISC @@ -44,6 +49,8 @@ from .precision_shrinkage_window import PRECISION_SHRINKAGE_WINDOW from .quantum_mutual_information import QUANTUM_MUTUAL_INFORMATION from .random_fourier_dependence import RANDOM_FOURIER_DEPENDENCE +from .random_fourier_dcc_combo import RANDOM_FOURIER_DCC_COMBO +from .random_fourier_time_freq_combo import RANDOM_FOURIER_TIME_FREQ_COMBO from .recurrence_kernel_dependence import RECURRENCE_KERNEL_DEPENDENCE from .reservoir_echo_state import RESERVOIR_ECHO_STATE from .robust_sliding_window import ROBUST_SLIDING_WINDOW @@ -104,6 +111,9 @@ "HMM_DISC", "EXPONENTIAL_WINDOW", "ADAPTIVE_EXPONENTIAL_WINDOW", + "ADAPTIVE_DCC_COMBO", + "ADAPTIVE_MULTISCALE_COMBO", + "ADAPTIVE_RANDOM_FOURIER_COMBO", "MULTISCALE_WINDOW", "EDGE_COACTIVATION", "PHASE_LOCKING_WINDOW", @@ -120,6 +130,8 @@ "PRECISION_SHRINKAGE_WINDOW", "RECURRENCE_KERNEL_DEPENDENCE", "RANDOM_FOURIER_DEPENDENCE", + "RANDOM_FOURIER_DCC_COMBO", + "RANDOM_FOURIER_TIME_FREQ_COMBO", "SPECTRAL_STATES", "EVENT_SYNCHRONIZATION", "COPULA_TAIL_DEPENDENCE", @@ -150,6 +162,8 @@ "POSITIVE_NEGATIVE_ASYMMETRY", "STATE_SPACE_NEIGHBORHOOD", "DCC_CONNECTIVITY", + "DCC_MULTISCALE_COMBO", + "DCC_TIME_FREQ_COMBO", "PERSISTENT_HOMOLOGY", "TIME_REVERSAL_ASYMMETRY", "QUANTUM_MUTUAL_INFORMATION", diff --git a/pydfc/dfc_methods/adaptive_dcc_combo.py b/pydfc/dfc_methods/adaptive_dcc_combo.py new file mode 100644 index 0000000..4bf78db --- /dev/null +++ b/pydfc/dfc_methods/adaptive_dcc_combo.py @@ -0,0 +1,145 @@ +""" +Adaptive DCC dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class ADAPTIVE_DCC_COMBO(BaseDFCMethod): + """DCC conditional correlation with signal-change-adaptive innovation weight.""" + + MEASURE_NAME = "AdaptiveDcccombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "garch_lambda", + "dcc_b", + "min_periods", + "alpha_min", + "alpha_max", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["garch_lambda"] is None: + self.params["garch_lambda"] = 0.94 + if self.params["dcc_b"] is None: + self.params["dcc_b"] = 0.85 + if self.params["min_periods"] is None: + self.params["min_periods"] = 10 + if self.params["alpha_min"] is None: + self.params["alpha_min"] = 0.02 + if self.params["alpha_max"] is None: + self.params["alpha_max"] = 0.20 + + @property + def measure_name(self): + return self.params["measure_name"] + + @staticmethod + def _corr_from_cov(covariance): + variance = np.diag(covariance) + scale = np.sqrt(np.outer(variance, variance)) + corr = np.divide( + covariance, + scale, + out=np.zeros_like(covariance, dtype=float), + where=scale > 0, + ) + corr = 0.5 * (corr + corr.T) + corr[np.diag_indices_from(corr)] = 1 + corr[np.isnan(corr)] = 0 + return np.clip(corr, -1.0, 1.0) + + @staticmethod + def _nearest_spd(matrix): + matrix = 0.5 * (matrix + matrix.T) + vals, vecs = np.linalg.eigh(matrix) + return vecs @ np.diag(np.maximum(vals, 1e-6)) @ vecs.T + + def _standardized_residuals(self, time_series): + n_regions, n_time = time_series.shape + lam = float(self.params["garch_lambda"]) + eps = time_series - time_series.mean(axis=1, keepdims=True) + sigma2 = np.zeros((n_regions, n_time)) + sigma2[:, 0] = np.maximum(np.var(eps, axis=1), 1e-8) + for tr in range(1, n_time): + sigma2[:, tr] = lam * sigma2[:, tr - 1] + (1.0 - lam) * eps[:, tr - 1] ** 2 + return eps / np.sqrt(np.maximum(sigma2, 1e-12)) + + def _adaptive_alpha(self, diff_energy, baseline, tr): + alpha_min = float(self.params["alpha_min"]) + alpha_max = float(self.params["alpha_max"]) + if not 0 <= alpha_min <= alpha_max <= 1: + raise ValueError( + "alpha_min and alpha_max must satisfy " + "0 <= alpha_min <= alpha_max <= 1." + ) + local_energy = diff_energy[tr] / baseline + adapt = local_energy / (1.0 + local_energy) + return alpha_min + (alpha_max - alpha_min) * adapt + + def dFC(self, time_series, Fs): + min_periods = int(self.params["min_periods"]) + beta = float(self.params["dcc_b"]) + if not 0 <= beta < 1: + raise ValueError("dcc_b must satisfy 0 <= dcc_b < 1.") + + Z = self._standardized_residuals(time_series) + Q_bar = self._nearest_spd(np.corrcoef(Z)) + Q_bar[np.isnan(Q_bar)] = 0 + Q_bar[np.diag_indices_from(Q_bar)] = 1 + Q = Q_bar.copy() + + diff_energy = np.zeros(time_series.shape[1]) + diff_energy[1:] = np.mean(np.abs(np.diff(time_series, axis=1)), axis=0) + baseline = np.median(diff_energy[1 : max(min_periods, 2)]) + 1e-8 + + FCSs = [] + TR_array = [] + for tr in range(1, time_series.shape[1]): + alpha = self._adaptive_alpha(diff_energy, baseline, tr) + z = Z[:, tr - 1] + Q = (1.0 - alpha) * ((1.0 - beta) * Q_bar + beta * Q) + alpha * np.outer(z, z) + R = self._corr_from_cov(Q) + if tr >= min_periods - 1: + FCSs.append(R) + TR_array.append(tr) + baseline = 0.98 * baseline + 0.02 * max(diff_energy[tr], 1e-8) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/dfc_methods/adaptive_multiscale_combo.py b/pydfc/dfc_methods/adaptive_multiscale_combo.py new file mode 100644 index 0000000..c9aca10 --- /dev/null +++ b/pydfc/dfc_methods/adaptive_multiscale_combo.py @@ -0,0 +1,141 @@ +""" +Adaptive multiscale window dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class ADAPTIVE_MULTISCALE_COMBO(BaseDFCMethod): + """Multiscale recent-window correlations with adaptive exponential weights.""" + + MEASURE_NAME = "AdaptiveMultiscalecombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "windows", + "min_periods", + "alpha_min", + "alpha_max", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["windows"] is None: + self.params["windows"] = [15, 30, 60] + if self.params["min_periods"] is None: + self.params["min_periods"] = 10 + if self.params["alpha_min"] is None: + self.params["alpha_min"] = 0.02 + if self.params["alpha_max"] is None: + self.params["alpha_max"] = 0.35 + + @property + def measure_name(self): + return self.params["measure_name"] + + @staticmethod + def _samples(value, Fs, minimum=1): + return max(int(round(value * Fs)), minimum) + + @staticmethod + def _corr_from_cov(covariance): + variance = np.diag(covariance) + scale = np.sqrt(np.outer(variance, variance)) + corr = np.divide( + covariance, + scale, + out=np.zeros_like(covariance, dtype=float), + where=scale > 0, + ) + corr = 0.5 * (corr + corr.T) + corr[np.diag_indices_from(corr)] = 1 + corr[np.isnan(corr)] = 0 + return np.clip(corr, -1.0, 1.0) + + def _weighted_corr(self, samples, weights): + weights = np.maximum(np.asarray(weights, dtype=float), 0) + weights_sum = np.sum(weights) + if weights_sum <= 0: + weights = np.ones(samples.shape[1], dtype=float) + weights_sum = np.sum(weights) + weights = weights / weights_sum + mean = np.sum(samples * weights[None, :], axis=1, keepdims=True) + centered = samples - mean + covariance = (centered * weights[None, :]) @ centered.T + return self._corr_from_cov(covariance) + + def _adaptive_alpha(self, diff_energy, baseline, tr): + alpha_min = float(self.params["alpha_min"]) + alpha_max = float(self.params["alpha_max"]) + if not 0 <= alpha_min <= alpha_max <= 1: + raise ValueError( + "alpha_min and alpha_max must satisfy " + "0 <= alpha_min <= alpha_max <= 1." + ) + local_energy = diff_energy[tr] / baseline + adapt = local_energy / (1.0 + local_energy) + return alpha_min + (alpha_max - alpha_min) * adapt + + def dFC(self, time_series, Fs): + windows = [self._samples(window, Fs, minimum=2) for window in self.params["windows"]] + min_periods = int(self.params["min_periods"]) + min_start = max(min(min(windows), min_periods), 2) + diff_energy = np.zeros(time_series.shape[1]) + diff_energy[1:] = np.mean(np.abs(np.diff(time_series, axis=1)), axis=0) + baseline = np.median(diff_energy[1 : max(min_start, 2)]) + 1e-8 + FCSs = [] + TR_array = [] + + for tr in range(min_start - 1, time_series.shape[1]): + alpha = self._adaptive_alpha(diff_energy, baseline, tr) + matrices = [] + scale_weights = [] + for window in windows: + start = max(0, tr - window + 1) + segment = time_series[:, start : tr + 1] + if segment.shape[1] < 2: + continue + age = segment.shape[1] - 1 - np.arange(segment.shape[1]) + weights = alpha * np.power(1.0 - alpha, age) + matrices.append(self._weighted_corr(segment, weights)) + scale_weights.append(np.sqrt(segment.shape[1])) + FCSs.append(np.average(np.array(matrices), axis=0, weights=scale_weights)) + TR_array.append(tr) + baseline = 0.98 * baseline + 0.02 * max(diff_energy[tr], 1e-8) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/dfc_methods/adaptive_random_fourier_combo.py b/pydfc/dfc_methods/adaptive_random_fourier_combo.py new file mode 100644 index 0000000..6173d4a --- /dev/null +++ b/pydfc/dfc_methods/adaptive_random_fourier_combo.py @@ -0,0 +1,140 @@ +""" +Adaptive random Fourier feature dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class ADAPTIVE_RANDOM_FOURIER_COMBO(BaseDFCMethod): + """Nonlinear random-feature dependence with volatility-adaptive forgetting.""" + + MEASURE_NAME = "AdaptiveRandomFouriercombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "min_periods", + "n_random_features", + "alpha_min", + "alpha_max", + "random_state", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + + if self.params["min_periods"] is None: + self.params["min_periods"] = 10 + if self.params["n_random_features"] is None: + self.params["n_random_features"] = 32 + if self.params["alpha_min"] is None: + self.params["alpha_min"] = 0.02 + if self.params["alpha_max"] is None: + self.params["alpha_max"] = 0.35 + + @property + def measure_name(self): + return self.params["measure_name"] + + def _feature_map(self, samples): + rng = np.random.default_rng(self.params["random_state"]) + n_features = max(int(self.params["n_random_features"]), 2) + omega = rng.normal(size=n_features) + phase = rng.uniform(0, 2 * np.pi, size=n_features) + return np.sqrt(2.0 / n_features) * np.cos( + samples[:, :, None] * omega + phase + ) + + def _instant_dependence(self, phi): + phi = phi - np.mean(phi, axis=1, keepdims=True) + norm = np.linalg.norm(phi, axis=1, keepdims=True) + phi = np.divide(phi, norm, out=np.zeros_like(phi), where=norm > 0) + instant = phi @ phi.T + instant = 0.5 * (instant + instant.T) + instant[np.diag_indices_from(instant)] = 1 + instant[np.isnan(instant)] = 0 + return instant + + def _adaptive_alpha(self, diff_energy, baseline, alpha_min, alpha_max, tr): + local_energy = diff_energy[tr] / baseline + adapt = local_energy / (1.0 + local_energy) + return alpha_min + (alpha_max - alpha_min) * adapt + + def dFC(self, time_series, Fs): + min_periods = int(self.params["min_periods"]) + alpha_min = float(self.params["alpha_min"]) + alpha_max = float(self.params["alpha_max"]) + if not 0 <= alpha_min <= alpha_max <= 1: + raise ValueError( + "alpha_min and alpha_max must satisfy " + "0 <= alpha_min <= alpha_max <= 1." + ) + + diff_energy = np.zeros(time_series.shape[1]) + diff_energy[1:] = np.mean(np.abs(np.diff(time_series, axis=1)), axis=0) + baseline = np.median(diff_energy[1 : max(min_periods, 2)]) + 1e-8 + features = self._feature_map(time_series) + n_regions = time_series.shape[0] + dependence = np.eye(n_regions) + FCSs = [] + TR_array = [] + + for tr in range(time_series.shape[1]): + alpha = self._adaptive_alpha( + diff_energy=diff_energy, + baseline=baseline, + alpha_min=alpha_min, + alpha_max=alpha_max, + tr=tr, + ) + instant = self._instant_dependence(features[:, tr, :]) + dependence = (1.0 - alpha) * dependence + alpha * instant + dependence = 0.5 * (dependence + dependence.T) + dependence[np.diag_indices_from(dependence)] = 1 + dependence[np.isnan(dependence)] = 0 + if tr >= min_periods - 1: + FCSs.append(dependence.copy()) + TR_array.append(tr) + baseline = 0.98 * baseline + 0.02 * max(diff_energy[tr], 1e-8) + + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert ( + len(time_series.subj_id_lst) == 1 + ), "this function takes only one subject as input." + assert ( + type(time_series) is TIME_SERIES + ), "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/dfc_methods/dcc_multiscale_combo.py b/pydfc/dfc_methods/dcc_multiscale_combo.py new file mode 100644 index 0000000..2a19838 --- /dev/null +++ b/pydfc/dfc_methods/dcc_multiscale_combo.py @@ -0,0 +1,144 @@ +""" +Multiscale DCC dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class DCC_MULTISCALE_COMBO(BaseDFCMethod): + """Average DCC conditional correlations over multiple recent scales.""" + + MEASURE_NAME = "DccMultiscalecombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "windows", + "garch_lambda", + "dcc_a", + "dcc_b", + "min_periods", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["windows"] is None: + self.params["windows"] = [15, 30, 60] + if self.params["garch_lambda"] is None: + self.params["garch_lambda"] = 0.94 + if self.params["dcc_a"] is None: + self.params["dcc_a"] = 0.05 + if self.params["dcc_b"] is None: + self.params["dcc_b"] = 0.85 + if self.params["min_periods"] is None: + self.params["min_periods"] = 10 + + @property + def measure_name(self): + return self.params["measure_name"] + + @staticmethod + def _samples(value, Fs, minimum=1): + return max(int(round(value * Fs)), minimum) + + @staticmethod + def _corr_from_cov(covariance): + variance = np.diag(covariance) + scale = np.sqrt(np.outer(variance, variance)) + corr = np.divide( + covariance, + scale, + out=np.zeros_like(covariance, dtype=float), + where=scale > 0, + ) + corr = 0.5 * (corr + corr.T) + corr[np.diag_indices_from(corr)] = 1 + corr[np.isnan(corr)] = 0 + return np.clip(corr, -1.0, 1.0) + + @staticmethod + def _nearest_spd(matrix): + matrix = 0.5 * (matrix + matrix.T) + vals, vecs = np.linalg.eigh(matrix) + return vecs @ np.diag(np.maximum(vals, 1e-6)) @ vecs.T + + def _standardized_residuals(self, time_series): + n_regions, n_time = time_series.shape + lam = float(self.params["garch_lambda"]) + eps = time_series - time_series.mean(axis=1, keepdims=True) + sigma2 = np.zeros((n_regions, n_time)) + sigma2[:, 0] = np.maximum(np.var(eps, axis=1), 1e-8) + for tr in range(1, n_time): + sigma2[:, tr] = lam * sigma2[:, tr - 1] + (1.0 - lam) * eps[:, tr - 1] ** 2 + return eps / np.sqrt(np.maximum(sigma2, 1e-12)) + + def _window_qbar(self, Z, tr, window): + start = max(0, tr - window + 1) + segment = Z[:, start : tr + 1] + if segment.shape[1] < 2: + return np.eye(Z.shape[0]) + qbar = np.corrcoef(segment) + qbar[np.isnan(qbar)] = 0 + qbar[np.diag_indices_from(qbar)] = 1 + return self._nearest_spd(qbar) + + def dFC(self, time_series, Fs): + windows = [self._samples(window, Fs, minimum=2) for window in self.params["windows"]] + min_periods = int(self.params["min_periods"]) + a = float(self.params["dcc_a"]) + b = float(self.params["dcc_b"]) + if not 0 <= a < 1 or not 0 <= b < 1 or a + b >= 1: + raise ValueError("dcc_a and dcc_b must be non-negative and sum to less than 1.") + + Z = self._standardized_residuals(time_series) + Qs = [self._window_qbar(Z, 0, window) for window in windows] + FCSs = [] + TR_array = [] + for tr in range(1, time_series.shape[1]): + matrices = [] + scale_weights = [] + z = Z[:, tr - 1] + for idx, window in enumerate(windows): + qbar = self._window_qbar(Z, tr, window) + Qs[idx] = (1.0 - a - b) * qbar + a * np.outer(z, z) + b * Qs[idx] + matrices.append(self._corr_from_cov(Qs[idx])) + scale_weights.append(np.sqrt(min(window, tr + 1))) + if tr >= min_periods - 1: + FCSs.append(np.average(np.array(matrices), axis=0, weights=scale_weights)) + TR_array.append(tr) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/dfc_methods/dcc_time_freq_combo.py b/pydfc/dfc_methods/dcc_time_freq_combo.py new file mode 100644 index 0000000..9610975 --- /dev/null +++ b/pydfc/dfc_methods/dcc_time_freq_combo.py @@ -0,0 +1,157 @@ +""" +Time-frequency DCC dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class DCC_TIME_FREQ_COMBO(BaseDFCMethod): + """DCC conditional correlations of local time-frequency power innovations.""" + + MEASURE_NAME = "DccTimeFreqcombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "tf_window", + "freq_bins", + "garch_lambda", + "dcc_a", + "dcc_b", + "min_periods", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["tf_window"] is None: + self.params["tf_window"] = 30 + if self.params["freq_bins"] is None: + self.params["freq_bins"] = 12 + if self.params["garch_lambda"] is None: + self.params["garch_lambda"] = 0.94 + if self.params["dcc_a"] is None: + self.params["dcc_a"] = 0.05 + if self.params["dcc_b"] is None: + self.params["dcc_b"] = 0.85 + if self.params["min_periods"] is None: + self.params["min_periods"] = 10 + + @property + def measure_name(self): + return self.params["measure_name"] + + @staticmethod + def _samples(value, Fs, minimum=1): + return max(int(round(value * Fs)), minimum) + + @staticmethod + def _corr_from_cov(covariance): + variance = np.diag(covariance) + scale = np.sqrt(np.outer(variance, variance)) + corr = np.divide( + covariance, + scale, + out=np.zeros_like(covariance, dtype=float), + where=scale > 0, + ) + corr = 0.5 * (corr + corr.T) + corr[np.diag_indices_from(corr)] = 1 + corr[np.isnan(corr)] = 0 + return np.clip(corr, -1.0, 1.0) + + @staticmethod + def _nearest_spd(matrix): + matrix = 0.5 * (matrix + matrix.T) + vals, vecs = np.linalg.eigh(matrix) + return vecs @ np.diag(np.maximum(vals, 1e-6)) @ vecs.T + + def _spectral_power(self, time_series, tr, window_samples): + start = max(0, tr - window_samples + 1) + segment = time_series[:, start : tr + 1] + segment = segment - segment.mean(axis=1, keepdims=True) + if segment.shape[1] < window_samples: + pad = np.zeros((segment.shape[0], window_samples - segment.shape[1])) + segment = np.hstack([pad, segment]) + taper = np.hanning(window_samples) + if not np.any(taper): + taper = np.ones(window_samples) + spectra = np.abs(np.fft.rfft(segment * taper[None, :], axis=1))[:, 1:] + n_bins = min(int(self.params["freq_bins"]), spectra.shape[1]) + power = np.sqrt(np.mean(spectra[:, :n_bins] ** 2, axis=1)) + return np.log1p(power) + + def _tf_innovations(self, time_series, Fs): + min_periods = int(self.params["min_periods"]) + window_samples = self._samples(self.params["tf_window"], Fs, minimum=max(min_periods, 4)) + features = np.array([ + self._spectral_power(time_series, tr, window_samples) + for tr in range(time_series.shape[1]) + ]).T + return features - features.mean(axis=1, keepdims=True) + + def _standardize(self, features): + n_regions, n_time = features.shape + lam = float(self.params["garch_lambda"]) + sigma2 = np.zeros((n_regions, n_time)) + sigma2[:, 0] = np.maximum(np.var(features, axis=1), 1e-8) + for tr in range(1, n_time): + sigma2[:, tr] = lam * sigma2[:, tr - 1] + (1.0 - lam) * features[:, tr - 1] ** 2 + return features / np.sqrt(np.maximum(sigma2, 1e-12)) + + def dFC(self, time_series, Fs): + min_periods = int(self.params["min_periods"]) + a = float(self.params["dcc_a"]) + b = float(self.params["dcc_b"]) + if not 0 <= a < 1 or not 0 <= b < 1 or a + b >= 1: + raise ValueError("dcc_a and dcc_b must be non-negative and sum to less than 1.") + + Z = self._standardize(self._tf_innovations(time_series, Fs)) + Q_bar = self._nearest_spd(np.corrcoef(Z)) + Q_bar[np.isnan(Q_bar)] = 0 + Q_bar[np.diag_indices_from(Q_bar)] = 1 + Q = Q_bar.copy() + FCSs = [] + TR_array = [] + for tr in range(1, time_series.shape[1]): + z = Z[:, tr - 1] + Q = (1.0 - a - b) * Q_bar + a * np.outer(z, z) + b * Q + R = self._corr_from_cov(Q) + if tr >= min_periods - 1: + FCSs.append(R) + TR_array.append(tr) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/dfc_methods/random_fourier_dcc_combo.py b/pydfc/dfc_methods/random_fourier_dcc_combo.py new file mode 100644 index 0000000..1cd51e5 --- /dev/null +++ b/pydfc/dfc_methods/random_fourier_dcc_combo.py @@ -0,0 +1,145 @@ +""" +Random Fourier DCC dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class RANDOM_FOURIER_DCC_COMBO(BaseDFCMethod): + """DCC smoothing of nonlinear random Fourier feature dependence.""" + + MEASURE_NAME = "RandomFourierDcccombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "garch_lambda", + "dcc_a", + "dcc_b", + "min_periods", + "n_random_features", + "random_state", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["garch_lambda"] is None: + self.params["garch_lambda"] = 0.94 + if self.params["dcc_a"] is None: + self.params["dcc_a"] = 0.05 + if self.params["dcc_b"] is None: + self.params["dcc_b"] = 0.85 + if self.params["min_periods"] is None: + self.params["min_periods"] = 10 + if self.params["n_random_features"] is None: + self.params["n_random_features"] = 32 + + @property + def measure_name(self): + return self.params["measure_name"] + + @staticmethod + def _corr_from_cov(covariance): + variance = np.diag(covariance) + scale = np.sqrt(np.outer(variance, variance)) + corr = np.divide( + covariance, + scale, + out=np.zeros_like(covariance, dtype=float), + where=scale > 0, + ) + corr = 0.5 * (corr + corr.T) + corr[np.diag_indices_from(corr)] = 1 + corr[np.isnan(corr)] = 0 + return np.clip(corr, -1.0, 1.0) + + @staticmethod + def _nearest_spd(matrix): + matrix = 0.5 * (matrix + matrix.T) + vals, vecs = np.linalg.eigh(matrix) + return vecs @ np.diag(np.maximum(vals, 1e-6)) @ vecs.T + + def _standardized_residuals(self, time_series): + n_regions, n_time = time_series.shape + lam = float(self.params["garch_lambda"]) + eps = time_series - time_series.mean(axis=1, keepdims=True) + sigma2 = np.zeros((n_regions, n_time)) + sigma2[:, 0] = np.maximum(np.var(eps, axis=1), 1e-8) + for tr in range(1, n_time): + sigma2[:, tr] = lam * sigma2[:, tr - 1] + (1.0 - lam) * eps[:, tr - 1] ** 2 + return eps / np.sqrt(np.maximum(sigma2, 1e-12)) + + def _feature_map(self, samples): + rng = np.random.default_rng(self.params["random_state"]) + n_features = max(int(self.params["n_random_features"]), 2) + omega = rng.normal(size=n_features) + phase = rng.uniform(0, 2 * np.pi, size=n_features) + return np.sqrt(2.0 / n_features) * np.cos(samples[:, :, None] * omega + phase) + + @staticmethod + def _instant_dependence(phi): + phi = phi - np.mean(phi, axis=1, keepdims=True) + norm = np.linalg.norm(phi, axis=1, keepdims=True) + phi = np.divide(phi, norm, out=np.zeros_like(phi), where=norm > 0) + instant = phi @ phi.T + instant = 0.5 * (instant + instant.T) + instant[np.diag_indices_from(instant)] = 1 + instant[np.isnan(instant)] = 0 + return instant + + def dFC(self, time_series, Fs): + min_periods = int(self.params["min_periods"]) + a = float(self.params["dcc_a"]) + b = float(self.params["dcc_b"]) + if not 0 <= a < 1 or not 0 <= b < 1 or a + b >= 1: + raise ValueError("dcc_a and dcc_b must be non-negative and sum to less than 1.") + + Z = self._standardized_residuals(time_series) + features = self._feature_map(Z) + instants = np.array([self._instant_dependence(features[:, tr, :]) for tr in range(Z.shape[1])]) + Q_bar = self._nearest_spd(np.mean(instants, axis=0)) + Q = Q_bar.copy() + FCSs = [] + TR_array = [] + for tr in range(1, Z.shape[1]): + Q = (1.0 - a - b) * Q_bar + a * instants[tr - 1] + b * Q + R = self._corr_from_cov(Q) + if tr >= min_periods - 1: + FCSs.append(R) + TR_array.append(tr) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/dfc_methods/random_fourier_time_freq_combo.py b/pydfc/dfc_methods/random_fourier_time_freq_combo.py new file mode 100644 index 0000000..c47569d --- /dev/null +++ b/pydfc/dfc_methods/random_fourier_time_freq_combo.py @@ -0,0 +1,124 @@ +""" +Random Fourier time-frequency dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class RANDOM_FOURIER_TIME_FREQ_COMBO(BaseDFCMethod): + """Random Fourier dependence on local time-frequency spectral profiles.""" + + MEASURE_NAME = "RandomFourierTimeFreqcombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "tf_window", + "min_periods", + "freq_bins", + "n_random_features", + "random_state", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["tf_window"] is None: + self.params["tf_window"] = 30 + if self.params["min_periods"] is None: + self.params["min_periods"] = 10 + if self.params["freq_bins"] is None: + self.params["freq_bins"] = 12 + if self.params["n_random_features"] is None: + self.params["n_random_features"] = 32 + + @property + def measure_name(self): + return self.params["measure_name"] + + @staticmethod + def _samples(value, Fs, minimum=1): + return max(int(round(value * Fs)), minimum) + + def _spectral_profiles(self, time_series, tr, window_samples): + start = max(0, tr - window_samples + 1) + segment = time_series[:, start : tr + 1] + segment = segment - segment.mean(axis=1, keepdims=True) + if segment.shape[1] < window_samples: + pad = np.zeros((segment.shape[0], window_samples - segment.shape[1])) + segment = np.hstack([pad, segment]) + taper = np.hanning(window_samples) + if not np.any(taper): + taper = np.ones(window_samples) + spectra = np.abs(np.fft.rfft(segment * taper[None, :], axis=1))[:, 1:] + n_bins = min(int(self.params["freq_bins"]), spectra.shape[1]) + profiles = spectra[:, :n_bins] + scale = np.linalg.norm(profiles, axis=1, keepdims=True) + return np.divide(profiles, scale, out=np.zeros_like(profiles), where=scale > 0) + + def _feature_projector(self, n_input): + rng = np.random.default_rng(self.params["random_state"]) + n_features = max(int(self.params["n_random_features"]), 2) + omega = rng.normal(size=(n_input, n_features)) + phase = rng.uniform(0, 2 * np.pi, size=n_features) + return omega, phase + + @staticmethod + def _dependence(phi): + phi = phi - np.mean(phi, axis=1, keepdims=True) + norm = np.linalg.norm(phi, axis=1, keepdims=True) + phi = np.divide(phi, norm, out=np.zeros_like(phi), where=norm > 0) + matrix = phi @ phi.T + matrix = 0.5 * (matrix + matrix.T) + matrix[np.diag_indices_from(matrix)] = 1 + matrix[np.isnan(matrix)] = 0 + return np.clip(matrix, -1.0, 1.0) + + def dFC(self, time_series, Fs): + min_periods = int(self.params["min_periods"]) + window_samples = self._samples(self.params["tf_window"], Fs, minimum=max(min_periods, 4)) + first_profiles = self._spectral_profiles(time_series, min_periods - 1, window_samples) + omega, phase = self._feature_projector(first_profiles.shape[1]) + FCSs = [] + TR_array = [] + + for tr in range(min_periods - 1, time_series.shape[1]): + profiles = self._spectral_profiles(time_series, tr, window_samples) + phi = np.sqrt(2.0 / omega.shape[1]) * np.cos(profiles @ omega + phase) + FCSs.append(self._dependence(phi)) + TR_array.append(tr) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/multi_analysis.py b/pydfc/multi_analysis.py index fa1a269..252a01c 100644 --- a/pydfc/multi_analysis.py +++ b/pydfc/multi_analysis.py @@ -162,6 +162,14 @@ def create_measure_obj(self, MEASURES_name_lst, **params): if MEASURES_name == "AdaptiveExponentialWindow": measure = ADAPTIVE_EXPONENTIAL_WINDOW(**params) + ###### ADAPTIVE DCC COMBO ###### + if MEASURES_name == "AdaptiveDcccombo": + measure = ADAPTIVE_DCC_COMBO(**params) + + ###### ADAPTIVE MULTISCALE COMBO ###### + if MEASURES_name == "AdaptiveMultiscalecombo": + measure = ADAPTIVE_MULTISCALE_COMBO(**params) + ###### MULTISCALE WINDOW ###### if MEASURES_name == "MultiscaleWindow": measure = MULTISCALE_WINDOW(**params) @@ -202,6 +210,26 @@ def create_measure_obj(self, MEASURES_name_lst, **params): if MEASURES_name == "RandomFourierDependence": measure = RANDOM_FOURIER_DEPENDENCE(**params) + ###### RANDOM FOURIER DCC COMBO ###### + if MEASURES_name == "RandomFourierDcccombo": + measure = RANDOM_FOURIER_DCC_COMBO(**params) + + ###### RANDOM FOURIER TIME FREQ COMBO ###### + if MEASURES_name == "RandomFourierTimeFreqcombo": + measure = RANDOM_FOURIER_TIME_FREQ_COMBO(**params) + + ###### ADAPTIVE RANDOM FOURIER COMBO ###### + if MEASURES_name == "AdaptiveRandomFouriercombo": + measure = ADAPTIVE_RANDOM_FOURIER_COMBO(**params) + + ###### DCC MULTISCALE COMBO ###### + if MEASURES_name == "DccMultiscalecombo": + measure = DCC_MULTISCALE_COMBO(**params) + + ###### DCC TIME FREQ COMBO ###### + if MEASURES_name == "DccTimeFreqcombo": + measure = DCC_TIME_FREQ_COMBO(**params) + ###### EVENT SYNCHRONIZATION ###### if MEASURES_name == "EventSynchronization": measure = EVENT_SYNCHRONIZATION(**params) diff --git a/tests/test_validation/dfc_method_wrappers.py b/tests/test_validation/dfc_method_wrappers.py index f575a23..4a837c8 100644 --- a/tests/test_validation/dfc_method_wrappers.py +++ b/tests/test_validation/dfc_method_wrappers.py @@ -10,6 +10,8 @@ "Time-Freq": ["tf", "wtc"], "ExponentialWindow": ["ew"], "AdaptiveExponentialWindow": ["aew"], + "AdaptiveDcccombo": ["adcccombo"], + "AdaptiveMultiscalecombo": ["amscombo"], "MultiscaleWindow": ["msw"], "EdgeCoactivation": ["eca"], "PhaseLockingWindow": ["plv"], @@ -21,6 +23,9 @@ "PrecisionShrinkageWindow": ["psw"], "RecurrenceKernelDependence": ["rkd"], "RandomFourierDependence": ["rfd"], + "AdaptiveRandomFouriercombo": ["arfc"], + "RandomFourierDcccombo": ["rfdcccombo"], + "RandomFourierTimeFreqcombo": ["rftfcombo"], "EventSynchronization": ["event"], "CopulaTailDependence": ["ctd"], "OjaSubspaceConnectivity": ["oja"], @@ -61,6 +66,8 @@ "PositiveNegativeAsymmetryFC": ["pnafc"], "StateSpaceNeighborhoodFC": ["ssnfc", "knnfc"], "DCCConnectivity": ["dcc"], + "DccMultiscalecombo": ["dccmscombo"], + "DccTimeFreqcombo": ["dcctfcombo"], "PersistentHomologyFC": ["phfc", "tda"], "TimeReversalAsymmetryFC": ["trac"], "QuantumMutualInformationFC": ["qmic"], From 4e7db7a843d47b04751dc434d82877a79659b33e Mon Sep 17 00:00:00 2001 From: Achille-V <72287888+Achille-V@users.noreply.github.com> Date: Mon, 20 Jul 2026 14:35:11 -0400 Subject: [PATCH 2/2] Adding 9 new combo methods --- dFC_method_pairs.txt | 15 +- pydfc/dfc_methods/__init__.py | 16 +++ .../adaptive_sliding_window_combo.py | 126 +++++++++++++++++ pydfc/dfc_methods/adaptive_time_freq_combo.py | 130 ++++++++++++++++++ pydfc/dfc_methods/dcc_sliding_window_combo.py | 124 +++++++++++++++++ .../multiscale_sliding_window_combo.py | 103 ++++++++++++++ .../dfc_methods/multiscale_time_freq_combo.py | 110 +++++++++++++++ .../random_fourier_multiscale_combo.py | 108 +++++++++++++++ .../random_fourier_sliding_window_combo.py | 107 ++++++++++++++ pydfc/dfc_methods/sliding_time_freq_combo.py | 101 ++++++++++++++ 10 files changed, 934 insertions(+), 6 deletions(-) create mode 100644 pydfc/dfc_methods/adaptive_sliding_window_combo.py create mode 100644 pydfc/dfc_methods/adaptive_time_freq_combo.py create mode 100644 pydfc/dfc_methods/dcc_sliding_window_combo.py create mode 100644 pydfc/dfc_methods/multiscale_sliding_window_combo.py create mode 100644 pydfc/dfc_methods/multiscale_time_freq_combo.py create mode 100644 pydfc/dfc_methods/random_fourier_multiscale_combo.py create mode 100644 pydfc/dfc_methods/random_fourier_sliding_window_combo.py create mode 100644 pydfc/dfc_methods/sliding_time_freq_combo.py diff --git a/dFC_method_pairs.txt b/dFC_method_pairs.txt index 7ff4f6c..ae89fda 100644 --- a/dFC_method_pairs.txt +++ b/dFC_method_pairs.txt @@ -1,6 +1,9 @@ -AdaptiveExponentialWindow and DCCConnectivity -AdaptiveExponentialWindow and MultiscaleWindow -RandomFourierDependence and DCCConnectivity -RandomFourierDependence and Time-Freq -DCCConnectivity and MultiscaleWindow -DCCConnectivity and Time-Freq +AdaptiveExponentialWindow and SlidingWindow +AdaptiveExponentialWindow and Time-Freq +RandomFourierDependence and MultiscaleWindow +RandomFourierDependence and Sliding Window +DCCConnectivity and SlidingWindow +MultiscaleWindow and SlidingWindow +MultiscaleWindow and Time-Freq +SlidingWindow and Time-Freq + diff --git a/pydfc/dfc_methods/__init__.py b/pydfc/dfc_methods/__init__.py index 3145f6e..b1cc78c 100644 --- a/pydfc/dfc_methods/__init__.py +++ b/pydfc/dfc_methods/__init__.py @@ -4,6 +4,8 @@ from .adaptive_dcc_combo import ADAPTIVE_DCC_COMBO from .adaptive_multiscale_combo import ADAPTIVE_MULTISCALE_COMBO from .adaptive_random_fourier_combo import ADAPTIVE_RANDOM_FOURIER_COMBO +from .adaptive_sliding_window_combo import ADAPTIVE_SLIDING_WINDOW_COMBO +from .adaptive_time_freq_combo import ADAPTIVE_TIME_FREQ_COMBO from .agglomerative_states import AGGLOMERATIVE_STATES from .amplitude_envelope_correlation import AMPLITUDE_ENVELOPE_CORRELATION from .base_dfc_method import BaseDFCMethod @@ -16,6 +18,7 @@ from .curvature_correlation import CURVATURE_CORRELATION from .dcc_connectivity import DCC_CONNECTIVITY from .dcc_multiscale_combo import DCC_MULTISCALE_COMBO +from .dcc_sliding_window_combo import DCC_SLIDING_WINDOW_COMBO from .dcc_time_freq_combo import DCC_TIME_FREQ_COMBO from .derivative_weighted_window import DERIVATIVE_WEIGHTED_WINDOW from .differential_coactivation import DIFFERENTIAL_COACTIVATION @@ -36,6 +39,8 @@ from .markov_smoothed_kmeans_states import MARKOV_SMOOTHED_KMEANS_STATES from .minibatch_kmeans_states import MINIBATCH_KMEANS_STATES from .multiscale_window import MULTISCALE_WINDOW +from .multiscale_sliding_window_combo import MULTISCALE_SLIDING_WINDOW_COMBO +from .multiscale_time_freq_combo import MULTISCALE_TIME_FREQ_COMBO from .mutual_compression import MUTUAL_COMPRESSION from .nmf_states import NMF_STATES from .oja_subspace_connectivity import OJA_SUBSPACE_CONNECTIVITY @@ -50,11 +55,14 @@ from .quantum_mutual_information import QUANTUM_MUTUAL_INFORMATION from .random_fourier_dependence import RANDOM_FOURIER_DEPENDENCE from .random_fourier_dcc_combo import RANDOM_FOURIER_DCC_COMBO +from .random_fourier_multiscale_combo import RANDOM_FOURIER_MULTISCALE_COMBO +from .random_fourier_sliding_window_combo import RANDOM_FOURIER_SLIDING_WINDOW_COMBO from .random_fourier_time_freq_combo import RANDOM_FOURIER_TIME_FREQ_COMBO from .recurrence_kernel_dependence import RECURRENCE_KERNEL_DEPENDENCE from .reservoir_echo_state import RESERVOIR_ECHO_STATE from .robust_sliding_window import ROBUST_SLIDING_WINDOW from .sliding_window import SLIDING_WINDOW +from .sliding_time_freq_combo import SLIDING_TIME_FREQ_COMBO from .sliding_window_clustr import SLIDING_WINDOW_CLUSTR from .sparse_coactivation_code import SPARSE_COACTIVATION_CODE from .spectral_similarity import SPECTRAL_SIMILARITY @@ -114,6 +122,8 @@ "ADAPTIVE_DCC_COMBO", "ADAPTIVE_MULTISCALE_COMBO", "ADAPTIVE_RANDOM_FOURIER_COMBO", + "ADAPTIVE_SLIDING_WINDOW_COMBO", + "ADAPTIVE_TIME_FREQ_COMBO", "MULTISCALE_WINDOW", "EDGE_COACTIVATION", "PHASE_LOCKING_WINDOW", @@ -131,6 +141,8 @@ "RECURRENCE_KERNEL_DEPENDENCE", "RANDOM_FOURIER_DEPENDENCE", "RANDOM_FOURIER_DCC_COMBO", + "RANDOM_FOURIER_MULTISCALE_COMBO", + "RANDOM_FOURIER_SLIDING_WINDOW_COMBO", "RANDOM_FOURIER_TIME_FREQ_COMBO", "SPECTRAL_STATES", "EVENT_SYNCHRONIZATION", @@ -138,6 +150,7 @@ "OJA_SUBSPACE_CONNECTIVITY", "GRAPH_DIFFUSION_COACTIVATION", "SLIDING_WINDOW", + "SLIDING_TIME_FREQ_COMBO", "TEMPORAL_DERIVATIVE_MULTIPLICATION", "TIME_FREQ", "WINDOWLESS", @@ -163,7 +176,10 @@ "STATE_SPACE_NEIGHBORHOOD", "DCC_CONNECTIVITY", "DCC_MULTISCALE_COMBO", + "DCC_SLIDING_WINDOW_COMBO", "DCC_TIME_FREQ_COMBO", + "MULTISCALE_SLIDING_WINDOW_COMBO", + "MULTISCALE_TIME_FREQ_COMBO", "PERSISTENT_HOMOLOGY", "TIME_REVERSAL_ASYMMETRY", "QUANTUM_MUTUAL_INFORMATION", diff --git a/pydfc/dfc_methods/adaptive_sliding_window_combo.py b/pydfc/dfc_methods/adaptive_sliding_window_combo.py new file mode 100644 index 0000000..5a8bb19 --- /dev/null +++ b/pydfc/dfc_methods/adaptive_sliding_window_combo.py @@ -0,0 +1,126 @@ +""" +Adaptive sliding-window dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class ADAPTIVE_SLIDING_WINDOW_COMBO(BaseDFCMethod): + """Sliding-window correlations with adaptive within-window recency weights.""" + + MEASURE_NAME = "AdaptiveSlidingWindowcombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "W", + "n_overlap", + "tapered_window", + "window_std", + "alpha_min", + "alpha_max", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["W"] is None: + self.params["W"] = 44 + if self.params["n_overlap"] is None: + self.params["n_overlap"] = 0.5 + if self.params["tapered_window"] is None: + self.params["tapered_window"] = True + if self.params["alpha_min"] is None: + self.params["alpha_min"] = 0.02 + if self.params["alpha_max"] is None: + self.params["alpha_max"] = 0.35 + + @property + def measure_name(self): + return self.params["measure_name"] + + @staticmethod + def _corr_from_cov(covariance): + variance = np.diag(covariance) + scale = np.sqrt(np.outer(variance, variance)) + corr = np.divide(covariance, scale, out=np.zeros_like(covariance), where=scale > 0) + corr = 0.5 * (corr + corr.T) + corr[np.diag_indices_from(corr)] = 1 + corr[np.isnan(corr)] = 0 + return np.clip(corr, -1.0, 1.0) + + def _weighted_corr(self, samples, weights): + weights = np.maximum(np.asarray(weights, dtype=float), 0) + if np.sum(weights) <= 0: + weights = np.ones(samples.shape[1], dtype=float) + weights = weights / np.sum(weights) + centered = samples - np.sum(samples * weights[None, :], axis=1, keepdims=True) + covariance = (centered * weights[None, :]) @ centered.T + return self._corr_from_cov(covariance) + + def _alpha(self, diff_energy, baseline, tr): + alpha_min = float(self.params["alpha_min"]) + alpha_max = float(self.params["alpha_max"]) + if not 0 <= alpha_min <= alpha_max <= 1: + raise ValueError("alpha_min and alpha_max must satisfy 0 <= alpha_min <= alpha_max <= 1.") + adapt = (diff_energy[tr] / baseline) / (1.0 + diff_energy[tr] / baseline) + return alpha_min + (alpha_max - alpha_min) * adapt + + def dFC(self, time_series, Fs): + n_time = time_series.shape[1] + W = max(int(round(float(self.params["W"]) * Fs)), 2) + W = min(W, n_time) + step = max(int(round((1.0 - float(self.params["n_overlap"])) * W)), 1) + diff_energy = np.zeros(n_time) + diff_energy[1:] = np.mean(np.abs(np.diff(time_series, axis=1)), axis=0) + baseline = np.median(diff_energy[1 : max(W, 2)]) + 1e-8 + FCSs = [] + TR_array = [] + for start in range(0, n_time - W + 1, step): + end = start + W + tr = int((start + end) / 2) + alpha = self._alpha(diff_energy, baseline, min(tr, n_time - 1)) + age = W - 1 - np.arange(W) + weights = alpha * np.power(1.0 - alpha, age) + if self.params["tapered_window"]: + taper = np.hanning(W) + if np.any(taper): + weights = weights * taper + FCSs.append(self._weighted_corr(time_series[:, start:end], weights)) + TR_array.append(tr) + baseline = 0.98 * baseline + 0.02 * max(diff_energy[min(tr, n_time - 1)], 1e-8) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/dfc_methods/adaptive_time_freq_combo.py b/pydfc/dfc_methods/adaptive_time_freq_combo.py new file mode 100644 index 0000000..dd985b6 --- /dev/null +++ b/pydfc/dfc_methods/adaptive_time_freq_combo.py @@ -0,0 +1,130 @@ +""" +Adaptive time-frequency dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class ADAPTIVE_TIME_FREQ_COMBO(BaseDFCMethod): + """Local spectral-profile dependence smoothed with adaptive forgetting.""" + + MEASURE_NAME = "AdaptiveTimeFreqcombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "tf_window", + "freq_bins", + "min_periods", + "alpha_min", + "alpha_max", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["tf_window"] is None: + self.params["tf_window"] = 30 + if self.params["freq_bins"] is None: + self.params["freq_bins"] = 12 + if self.params["min_periods"] is None: + self.params["min_periods"] = 10 + if self.params["alpha_min"] is None: + self.params["alpha_min"] = 0.02 + if self.params["alpha_max"] is None: + self.params["alpha_max"] = 0.35 + + @property + def measure_name(self): + return self.params["measure_name"] + + @staticmethod + def _samples(value, Fs, minimum=1): + return max(int(round(float(value) * Fs)), minimum) + + @staticmethod + def _dependence(features): + features = features - np.mean(features, axis=1, keepdims=True) + scale = np.linalg.norm(features, axis=1, keepdims=True) + features = np.divide(features, scale, out=np.zeros_like(features), where=scale > 0) + matrix = features @ features.T + matrix = 0.5 * (matrix + matrix.T) + matrix[np.diag_indices_from(matrix)] = 1 + matrix[np.isnan(matrix)] = 0 + return np.clip(matrix, -1.0, 1.0) + + def _spectral_profiles(self, time_series, tr, window_samples): + start = max(0, tr - window_samples + 1) + segment = time_series[:, start : tr + 1] + segment = segment - segment.mean(axis=1, keepdims=True) + if segment.shape[1] < window_samples: + segment = np.hstack([np.zeros((segment.shape[0], window_samples - segment.shape[1])), segment]) + taper = np.hanning(window_samples) + if not np.any(taper): + taper = np.ones(window_samples) + spectra = np.abs(np.fft.rfft(segment * taper[None, :], axis=1))[:, 1:] + n_bins = min(int(self.params["freq_bins"]), spectra.shape[1]) + return np.log1p(spectra[:, :n_bins]) + + def _alpha(self, diff_energy, baseline, tr): + alpha_min = float(self.params["alpha_min"]) + alpha_max = float(self.params["alpha_max"]) + if not 0 <= alpha_min <= alpha_max <= 1: + raise ValueError("alpha_min and alpha_max must satisfy 0 <= alpha_min <= alpha_max <= 1.") + ratio = diff_energy[tr] / baseline + return alpha_min + (alpha_max - alpha_min) * (ratio / (1.0 + ratio)) + + def dFC(self, time_series, Fs): + min_periods = int(self.params["min_periods"]) + window_samples = self._samples(self.params["tf_window"], Fs, minimum=max(min_periods, 4)) + diff_energy = np.zeros(time_series.shape[1]) + diff_energy[1:] = np.mean(np.abs(np.diff(time_series, axis=1)), axis=0) + baseline = np.median(diff_energy[1 : max(min_periods, 2)]) + 1e-8 + smoothed = np.eye(time_series.shape[0]) + FCSs = [] + TR_array = [] + for tr in range(time_series.shape[1]): + profiles = self._spectral_profiles(time_series, tr, window_samples) + alpha = self._alpha(diff_energy, baseline, tr) + smoothed = (1.0 - alpha) * smoothed + alpha * self._dependence(profiles) + smoothed = 0.5 * (smoothed + smoothed.T) + smoothed[np.diag_indices_from(smoothed)] = 1 + if tr >= min_periods - 1: + FCSs.append(smoothed.copy()) + TR_array.append(tr) + baseline = 0.98 * baseline + 0.02 * max(diff_energy[tr], 1e-8) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/dfc_methods/dcc_sliding_window_combo.py b/pydfc/dfc_methods/dcc_sliding_window_combo.py new file mode 100644 index 0000000..46e87d2 --- /dev/null +++ b/pydfc/dfc_methods/dcc_sliding_window_combo.py @@ -0,0 +1,124 @@ +""" +DCC sliding-window dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class DCC_SLIDING_WINDOW_COMBO(BaseDFCMethod): + """Sliding-window summaries of DCC conditional correlation trajectories.""" + + MEASURE_NAME = "DccSlidingWindowcombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "W", + "n_overlap", + "garch_lambda", + "dcc_a", + "dcc_b", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["W"] is None: + self.params["W"] = 44 + if self.params["n_overlap"] is None: + self.params["n_overlap"] = 0.5 + if self.params["garch_lambda"] is None: + self.params["garch_lambda"] = 0.94 + if self.params["dcc_a"] is None: + self.params["dcc_a"] = 0.05 + if self.params["dcc_b"] is None: + self.params["dcc_b"] = 0.85 + + @property + def measure_name(self): + return self.params["measure_name"] + + @staticmethod + def _corr_from_cov(covariance): + variance = np.diag(covariance) + scale = np.sqrt(np.outer(variance, variance)) + corr = np.divide(covariance, scale, out=np.zeros_like(covariance), where=scale > 0) + corr = 0.5 * (corr + corr.T) + corr[np.diag_indices_from(corr)] = 1 + corr[np.isnan(corr)] = 0 + return np.clip(corr, -1.0, 1.0) + + @staticmethod + def _nearest_spd(matrix): + matrix = 0.5 * (matrix + matrix.T) + vals, vecs = np.linalg.eigh(matrix) + return vecs @ np.diag(np.maximum(vals, 1e-6)) @ vecs.T + + def _dcc_path(self, time_series): + a = float(self.params["dcc_a"]) + b = float(self.params["dcc_b"]) + if not 0 <= a < 1 or not 0 <= b < 1 or a + b >= 1: + raise ValueError("dcc_a and dcc_b must be non-negative and sum to less than 1.") + lam = float(self.params["garch_lambda"]) + eps = time_series - time_series.mean(axis=1, keepdims=True) + sigma2 = np.zeros_like(eps) + sigma2[:, 0] = np.maximum(np.var(eps, axis=1), 1e-8) + for tr in range(1, eps.shape[1]): + sigma2[:, tr] = lam * sigma2[:, tr - 1] + (1.0 - lam) * eps[:, tr - 1] ** 2 + Z = eps / np.sqrt(np.maximum(sigma2, 1e-12)) + Q_bar = self._nearest_spd(np.corrcoef(Z)) + Q_bar[np.isnan(Q_bar)] = 0 + Q_bar[np.diag_indices_from(Q_bar)] = 1 + Q = Q_bar.copy() + matrices = [self._corr_from_cov(Q)] + for tr in range(1, Z.shape[1]): + Q = (1.0 - a - b) * Q_bar + a * np.outer(Z[:, tr - 1], Z[:, tr - 1]) + b * Q + matrices.append(self._corr_from_cov(Q)) + return np.array(matrices) + + def dFC(self, time_series, Fs): + path = self._dcc_path(time_series) + n_time = time_series.shape[1] + W = min(max(int(round(float(self.params["W"]) * Fs)), 2), n_time) + step = max(int(round((1.0 - float(self.params["n_overlap"])) * W)), 1) + FCSs = [] + TR_array = [] + for start in range(0, n_time - W + 1, step): + end = start + W + FCSs.append(np.mean(path[start:end], axis=0)) + TR_array.append(int((start + end) / 2)) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/dfc_methods/multiscale_sliding_window_combo.py b/pydfc/dfc_methods/multiscale_sliding_window_combo.py new file mode 100644 index 0000000..09dff66 --- /dev/null +++ b/pydfc/dfc_methods/multiscale_sliding_window_combo.py @@ -0,0 +1,103 @@ +""" +Multiscale sliding-window dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class MULTISCALE_SLIDING_WINDOW_COMBO(BaseDFCMethod): + """Sliding-window correlations averaged over multiple centered window sizes.""" + + MEASURE_NAME = "MultiscaleSlidingWindowcombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "windows", + "n_overlap", + "tapered_window", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["windows"] is None: + self.params["windows"] = [15, 30, 60] + if self.params["n_overlap"] is None: + self.params["n_overlap"] = 0.5 + if self.params["tapered_window"] is None: + self.params["tapered_window"] = True + + @property + def measure_name(self): + return self.params["measure_name"] + + @staticmethod + def _samples(value, Fs, minimum=1): + return max(int(round(float(value) * Fs)), minimum) + + @staticmethod + def _corr(samples): + corr = np.corrcoef(samples) + corr[np.isnan(corr)] = 0 + corr[np.diag_indices_from(corr)] = 1 + return np.clip(0.5 * (corr + corr.T), -1.0, 1.0) + + def dFC(self, time_series, Fs): + n_time = time_series.shape[1] + windows = [min(self._samples(window, Fs, minimum=2), n_time) for window in self.params["windows"]] + base_window = max(windows) + step = max(int(round((1.0 - float(self.params["n_overlap"])) * base_window)), 1) + FCSs = [] + TR_array = [] + for center in range(base_window // 2, n_time - (base_window - base_window // 2) + 1, step): + matrices = [] + weights = [] + for window in windows: + start = max(0, center - window // 2) + end = min(n_time, start + window) + start = max(0, end - window) + segment = time_series[:, start:end] + if self.params["tapered_window"]: + taper = np.hanning(segment.shape[1]) + if np.any(taper): + segment = segment * taper[None, :] + matrices.append(self._corr(segment)) + weights.append(np.sqrt(segment.shape[1])) + FCSs.append(np.average(np.array(matrices), axis=0, weights=weights)) + TR_array.append(center) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/dfc_methods/multiscale_time_freq_combo.py b/pydfc/dfc_methods/multiscale_time_freq_combo.py new file mode 100644 index 0000000..fe1737d --- /dev/null +++ b/pydfc/dfc_methods/multiscale_time_freq_combo.py @@ -0,0 +1,110 @@ +""" +Multiscale time-frequency dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class MULTISCALE_TIME_FREQ_COMBO(BaseDFCMethod): + """Time-frequency spectral-profile dependence averaged across window scales.""" + + MEASURE_NAME = "MultiscaleTimeFreqcombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "windows", + "freq_bins", + "min_periods", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["windows"] is None: + self.params["windows"] = [15, 30, 60] + if self.params["freq_bins"] is None: + self.params["freq_bins"] = 12 + if self.params["min_periods"] is None: + self.params["min_periods"] = 10 + + @property + def measure_name(self): + return self.params["measure_name"] + + @staticmethod + def _samples(value, Fs, minimum=1): + return max(int(round(float(value) * Fs)), minimum) + + @staticmethod + def _dependence(features): + features = features - np.mean(features, axis=1, keepdims=True) + scale = np.linalg.norm(features, axis=1, keepdims=True) + features = np.divide(features, scale, out=np.zeros_like(features), where=scale > 0) + matrix = features @ features.T + matrix = 0.5 * (matrix + matrix.T) + matrix[np.diag_indices_from(matrix)] = 1 + matrix[np.isnan(matrix)] = 0 + return np.clip(matrix, -1.0, 1.0) + + def _profiles(self, time_series, tr, window): + start = max(0, tr - window + 1) + segment = time_series[:, start : tr + 1] + segment = segment - segment.mean(axis=1, keepdims=True) + if segment.shape[1] < window: + segment = np.hstack([np.zeros((segment.shape[0], window - segment.shape[1])), segment]) + taper = np.hanning(window) + if not np.any(taper): + taper = np.ones(window) + spectra = np.abs(np.fft.rfft(segment * taper[None, :], axis=1))[:, 1:] + n_bins = min(int(self.params["freq_bins"]), spectra.shape[1]) + return np.log1p(spectra[:, :n_bins]) + + def dFC(self, time_series, Fs): + windows = [self._samples(window, Fs, minimum=4) for window in self.params["windows"]] + min_periods = max(int(self.params["min_periods"]), min(windows)) + FCSs = [] + TR_array = [] + for tr in range(min_periods - 1, time_series.shape[1]): + matrices = [] + weights = [] + for window in windows: + matrices.append(self._dependence(self._profiles(time_series, tr, window))) + weights.append(np.sqrt(min(window, tr + 1))) + FCSs.append(np.average(np.array(matrices), axis=0, weights=weights)) + TR_array.append(tr) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/dfc_methods/random_fourier_multiscale_combo.py b/pydfc/dfc_methods/random_fourier_multiscale_combo.py new file mode 100644 index 0000000..d27b32c --- /dev/null +++ b/pydfc/dfc_methods/random_fourier_multiscale_combo.py @@ -0,0 +1,108 @@ +""" +Random Fourier multiscale dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class RANDOM_FOURIER_MULTISCALE_COMBO(BaseDFCMethod): + """Random Fourier feature dependence averaged across recent scales.""" + + MEASURE_NAME = "RandomFourierMultiscalecombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "windows", + "min_periods", + "n_random_features", + "random_state", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["windows"] is None: + self.params["windows"] = [15, 30, 60] + if self.params["min_periods"] is None: + self.params["min_periods"] = 10 + if self.params["n_random_features"] is None: + self.params["n_random_features"] = 32 + + @property + def measure_name(self): + return self.params["measure_name"] + + @staticmethod + def _samples(value, Fs, minimum=1): + return max(int(round(float(value) * Fs)), minimum) + + def _projector(self): + rng = np.random.default_rng(self.params["random_state"]) + n_features = max(int(self.params["n_random_features"]), 2) + return rng.normal(size=n_features), rng.uniform(0, 2 * np.pi, size=n_features) + + @staticmethod + def _dependence(phi): + phi = phi.reshape(phi.shape[0], -1) + phi = phi - np.mean(phi, axis=1, keepdims=True) + norm = np.linalg.norm(phi, axis=1, keepdims=True) + phi = np.divide(phi, norm, out=np.zeros_like(phi), where=norm > 0) + matrix = phi @ phi.T + matrix = 0.5 * (matrix + matrix.T) + matrix[np.diag_indices_from(matrix)] = 1 + matrix[np.isnan(matrix)] = 0 + return np.clip(matrix, -1.0, 1.0) + + def dFC(self, time_series, Fs): + windows = [self._samples(window, Fs, minimum=2) for window in self.params["windows"]] + min_periods = max(int(self.params["min_periods"]), min(windows)) + omega, phase = self._projector() + FCSs = [] + TR_array = [] + for tr in range(min_periods - 1, time_series.shape[1]): + matrices = [] + weights = [] + for window in windows: + start = max(0, tr - window + 1) + segment = time_series[:, start : tr + 1] + phi = np.sqrt(2.0 / len(omega)) * np.cos(segment[:, :, None] * omega + phase) + matrices.append(self._dependence(phi)) + weights.append(np.sqrt(segment.shape[1])) + FCSs.append(np.average(np.array(matrices), axis=0, weights=weights)) + TR_array.append(tr) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/dfc_methods/random_fourier_sliding_window_combo.py b/pydfc/dfc_methods/random_fourier_sliding_window_combo.py new file mode 100644 index 0000000..5e5c707 --- /dev/null +++ b/pydfc/dfc_methods/random_fourier_sliding_window_combo.py @@ -0,0 +1,107 @@ +""" +Random Fourier sliding-window dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class RANDOM_FOURIER_SLIDING_WINDOW_COMBO(BaseDFCMethod): + """Sliding-window dependence after nonlinear random Fourier projection.""" + + MEASURE_NAME = "RandomFourierSlidingWindowcombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "W", + "n_overlap", + "tapered_window", + "n_random_features", + "random_state", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["W"] is None: + self.params["W"] = 44 + if self.params["n_overlap"] is None: + self.params["n_overlap"] = 0.5 + if self.params["tapered_window"] is None: + self.params["tapered_window"] = True + if self.params["n_random_features"] is None: + self.params["n_random_features"] = 32 + + @property + def measure_name(self): + return self.params["measure_name"] + + def _projector(self): + rng = np.random.default_rng(self.params["random_state"]) + n_features = max(int(self.params["n_random_features"]), 2) + return rng.normal(size=n_features), rng.uniform(0, 2 * np.pi, size=n_features) + + @staticmethod + def _dependence(phi): + phi = phi.reshape(phi.shape[0], -1) + phi = phi - np.mean(phi, axis=1, keepdims=True) + norm = np.linalg.norm(phi, axis=1, keepdims=True) + phi = np.divide(phi, norm, out=np.zeros_like(phi), where=norm > 0) + matrix = phi @ phi.T + matrix = 0.5 * (matrix + matrix.T) + matrix[np.diag_indices_from(matrix)] = 1 + matrix[np.isnan(matrix)] = 0 + return np.clip(matrix, -1.0, 1.0) + + def dFC(self, time_series, Fs): + n_time = time_series.shape[1] + W = min(max(int(round(float(self.params["W"]) * Fs)), 2), n_time) + step = max(int(round((1.0 - float(self.params["n_overlap"])) * W)), 1) + omega, phase = self._projector() + FCSs = [] + TR_array = [] + for start in range(0, n_time - W + 1, step): + end = start + W + segment = time_series[:, start:end] + if self.params["tapered_window"]: + taper = np.hanning(W) + if np.any(taper): + segment = segment * taper[None, :] + phi = np.sqrt(2.0 / len(omega)) * np.cos(segment[:, :, None] * omega + phase) + FCSs.append(self._dependence(phi)) + TR_array.append(int((start + end) / 2)) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC diff --git a/pydfc/dfc_methods/sliding_time_freq_combo.py b/pydfc/dfc_methods/sliding_time_freq_combo.py new file mode 100644 index 0000000..7938a3f --- /dev/null +++ b/pydfc/dfc_methods/sliding_time_freq_combo.py @@ -0,0 +1,101 @@ +""" +Sliding time-frequency dFC method. +""" + +import time + +import numpy as np + +from ..dfc import DFC +from ..time_series import TIME_SERIES +from .base_dfc_method import BaseDFCMethod + + +class SLIDING_TIME_FREQ_COMBO(BaseDFCMethod): + """Sliding-window dependence of local spectral profiles.""" + + MEASURE_NAME = "SlidingTimeFreqcombo" + + def __init__(self, **params): + self.logs_ = "" + self.TPM = [] + self.FCS_ = [] + self.FCS_fit_time_ = None + self.dFC_assess_time_ = None + self.params_name_lst = [ + "measure_name", + "is_state_based", + "W", + "n_overlap", + "freq_bins", + "tapered_window", + "normalization", + "num_select_nodes", + "num_time_point", + "Fs_ratio", + "noise_ratio", + "num_realization", + "session", + ] + self.params = {} + for params_name in self.params_name_lst: + self.params[params_name] = params.get(params_name, None) + self.params["measure_name"] = self.MEASURE_NAME + self.params["is_state_based"] = False + if self.params["W"] is None: + self.params["W"] = 44 + if self.params["n_overlap"] is None: + self.params["n_overlap"] = 0.5 + if self.params["freq_bins"] is None: + self.params["freq_bins"] = 12 + if self.params["tapered_window"] is None: + self.params["tapered_window"] = True + + @property + def measure_name(self): + return self.params["measure_name"] + + @staticmethod + def _dependence(features): + features = features - np.mean(features, axis=1, keepdims=True) + scale = np.linalg.norm(features, axis=1, keepdims=True) + features = np.divide(features, scale, out=np.zeros_like(features), where=scale > 0) + matrix = features @ features.T + matrix = 0.5 * (matrix + matrix.T) + matrix[np.diag_indices_from(matrix)] = 1 + matrix[np.isnan(matrix)] = 0 + return np.clip(matrix, -1.0, 1.0) + + def dFC(self, time_series, Fs): + n_time = time_series.shape[1] + W = min(max(int(round(float(self.params["W"]) * Fs)), 4), n_time) + step = max(int(round((1.0 - float(self.params["n_overlap"])) * W)), 1) + FCSs = [] + TR_array = [] + for start in range(0, n_time - W + 1, step): + end = start + W + segment = time_series[:, start:end] + segment = segment - segment.mean(axis=1, keepdims=True) + if self.params["tapered_window"]: + taper = np.hanning(W) + if np.any(taper): + segment = segment * taper[None, :] + spectra = np.abs(np.fft.rfft(segment, axis=1))[:, 1:] + n_bins = min(int(self.params["freq_bins"]), spectra.shape[1]) + FCSs.append(self._dependence(np.log1p(spectra[:, :n_bins]))) + TR_array.append(int((start + end) / 2)) + return np.array(FCSs), np.array(TR_array) + + def estimate_FCS(self, time_series): + return self + + def estimate_dFC(self, time_series): + assert len(time_series.subj_id_lst) == 1, "this function takes only one subject as input." + assert type(time_series) is TIME_SERIES, "time_series must be of TIME_SERIES class." + time_series = self.manipulate_time_series4dFC(time_series) + tic = time.time() + FCSs, TR_array = self.dFC(time_series=time_series.data, Fs=time_series.Fs) + self.set_dFC_assess_time(time.time() - tic) + dFC = DFC(measure=self) + dFC.set_dFC(FCSs=FCSs, TR_array=TR_array, TS_info=time_series.info_dict) + return dFC