Adding combo dFC methods

dFC_method_pairs.txt

@@ -0,0 +1,9 @@
+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
+

pydfc/dfc_methods/__init__.py

@@ -1,6 +1,11 @@
 """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 .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
@@ -12,6 +17,9 @@
 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_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
 from .discrete_hmm import HMM_DISC
@@ -31,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
@@ -44,10 +54,15 @@
 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_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
@@ -104,6 +119,11 @@
     "HMM_DISC",
     "EXPONENTIAL_WINDOW",
     "ADAPTIVE_EXPONENTIAL_WINDOW",
+    "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",
@@ -120,12 +140,17 @@
     "PRECISION_SHRINKAGE_WINDOW",
     "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",
     "COPULA_TAIL_DEPENDENCE",
     "OJA_SUBSPACE_CONNECTIVITY",
     "GRAPH_DIFFUSION_COACTIVATION",
     "SLIDING_WINDOW",
+    "SLIDING_TIME_FREQ_COMBO",
     "TEMPORAL_DERIVATIVE_MULTIPLICATION",
     "TIME_FREQ",
     "WINDOWLESS",
@@ -150,6 +175,11 @@
     "POSITIVE_NEGATIVE_ASYMMETRY",
     "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",

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

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