diff --git a/pyspi/config.yaml b/pyspi/config.yaml index 5bb5855..fd215af 100644 --- a/pyspi/config.yaml +++ b/pyspi/config.yaml @@ -1391,3 +1391,15 @@ - orth: True log: True absolute: True + + # Interdependence score + InterDependenceScore: + labels: + - unsigned + - undirected + - nonlinear + dependencies: + configs: # default params + - terms: 6 + pnorm: 'max' + bandwidth: 0.5 diff --git a/pyspi/fast_config.yaml b/pyspi/fast_config.yaml index e05f1cf..9fbe31e 100644 --- a/pyspi/fast_config.yaml +++ b/pyspi/fast_config.yaml @@ -1291,3 +1291,16 @@ - orth: True log: True absolute: True + + # Interdependence score + InterDependenceScore: + labels: + - unsigned + - undirected + - nonlinear + dependencies: + configs: # default params + - terms: 6 + pnorm: 'max' + bandwidth: 0.5 + \ No newline at end of file diff --git a/pyspi/lib/ids/LICENSE.txt b/pyspi/lib/ids/LICENSE.txt new file mode 100644 index 0000000..7f1e5c2 --- /dev/null +++ b/pyspi/lib/ids/LICENSE.txt @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2025 Adityanarayanan Radhakrishnan + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/pyspi/lib/ids/__init__.py b/pyspi/lib/ids/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/pyspi/lib/ids/dependence.py b/pyspi/lib/ids/dependence.py new file mode 100644 index 0000000..7402132 --- /dev/null +++ b/pyspi/lib/ids/dependence.py @@ -0,0 +1,31 @@ +""" +Interdependence Score (IDS) computation. +Based on the work by Adityanarayanan Radhakrishnan (MIT License) +Original: https://github.com/aradha/interdependence_scores +Modified for use in pyspi package. +""" +from .numpy_dependence import compute_IDS_numpy + +def compute_IDS(X, Y=None, num_terms=6, p_norm='max', + p_val=False, num_tests=100, bandwidth_term=1/2): + """Compute IDS between all pairs of variables in X (or between X and Y). + + This is a modified version of the implementation from: + https://github.com/aradha/interdependence_scores + + Original author: Adityanarayanan Radhakrishnan + License: MIT (see LICENSE.txt) + + Parameters: + X: np.ndarray or torch.Tensor + Y: np.ndarray or torch.Tensor (optional) + num_terms: Number of terms for Taylor series approximation (optional) + p_norm: String 'max' if using IDS-max. 1 or 2 for IDS-1, IDS-2, respectively. (optional) + p_val: Boolean. Indicates whether to compute p-values using permutation tests + num_tests: Number of permutation tests if p_val=True + bandwidth_term: Constant term in Gaussian kernel + Returns: + IDS matrix, p-value matrix (if p_val=True) + """ + return compute_IDS_numpy(X, Y=Y, num_terms=num_terms, p_norm=p_norm, + p_val=p_val, num_tests=num_tests, bandwidth_term=bandwidth_term) diff --git a/pyspi/lib/ids/numpy_dependence.py b/pyspi/lib/ids/numpy_dependence.py new file mode 100644 index 0000000..09f3e5d --- /dev/null +++ b/pyspi/lib/ids/numpy_dependence.py @@ -0,0 +1,99 @@ +""" +Interdependence Score (IDS) computation. +Based on the work by Adityanarayanan Radhakrishnan (MIT License) +Original: https://github.com/aradha/interdependence_scores +""" +import numpy as np +import math +import sys +from tqdm import tqdm + +SEED = 1717 +np.random.seed(SEED) + +EPSILON = sys.float_info.epsilon + +def transform(y, num_terms=6, bandwidth_term=1/2): + B = bandwidth_term + exp = np.exp(-B * y**2) + terms = [] + for i in range(num_terms): + terms.append(exp * (y)**i / math.sqrt(math.factorial(i) *1.)) + y_ = np.concatenate(terms, axis=-1) + return y_ + +def center(X): + return X - np.mean(X, axis=0, keepdims=True) + + +def compute_p_val(C, X, Y=None, num_terms=6, p_norm='max', n_tests=100, bandwidth_term=1/2): + + gt = C + count = 0 + + n, dx = X.shape + for i in tqdm(range(n_tests)): + + # Used to shuffle data + random_noise = np.random.normal(size=(n, dx)) + permutations = np.argsort(random_noise, axis=0) + X_permuted = X[permutations, np.arange(dx)[None, :]] + + if Y is not None: + n, dy = Y.shape + random_noise = np.random.normal(size=(n, dy)) + permutations = np.argsort(random_noise, axis=0) + Y_permuted = Y[permutations, np.arange(dy)[None, :]] + null = compute_IDS_numpy(X_permuted, Y=Y_permuted, num_terms=num_terms, + p_norm=p_norm, bandwidth_term=bandwidth_term) + else: + null = compute_IDS_numpy(X_permuted, Y=Y, num_terms=num_terms, + p_norm=p_norm, bandwidth_term=bandwidth_term) + + + count += np.where(null > gt, 1, 0) + + p_vals = count / n_tests + return p_vals + + +def compute_IDS_numpy(X, Y=None, num_terms=6, p_norm='max', + p_val=False, num_tests=100, bandwidth_term=1/2): + n, dx = X.shape + X_t = transform(X, num_terms=num_terms, bandwidth_term=bandwidth_term) + X_t = center(X_t) + + if Y is not None: + _, dy = Y.shape + Y_t = transform(Y, num_terms=num_terms, bandwidth_term=bandwidth_term) + Y_t = center(Y_t) + cov = X_t.T @ Y_t + X_std = np.sqrt(np.sum(X_t**2, axis=0)) + Y_std = np.sqrt(np.sum(Y_t**2, axis=0)) + correlations = cov / (X_std.reshape(-1, 1) + EPSILON) + C = correlations / (Y_std.reshape(1, -1) + EPSILON) + C = C.reshape(num_terms, dx, num_terms, dy) + else: + C = np.corrcoef(X_t.T) + C = C.reshape(num_terms, dx, num_terms, dx) + + C = np.nan_to_num(C, nan=0, posinf=0, neginf=0) + C = np.abs(C) + + if p_norm == 'max': + C = np.amax(C, axis=(0, 2)) + elif p_norm == 2: + C = C**2 + C = np.mean(C, axis=0) + C = np.mean(C, axis=1) + C = np.sqrt(C) + elif p_norm == 1: + C = np.mean(C, axis=0) + C = np.mean(C, axis=1) + + if p_val: + p_vals = compute_p_val(C, X, Y=Y, num_terms=num_terms, p_norm=p_norm, + n_tests=num_tests, bandwidth_term=bandwidth_term) + return C, p_vals + else: + return C diff --git a/pyspi/statistics/misc.py b/pyspi/statistics/misc.py index a99a48d..b885b49 100644 --- a/pyspi/statistics/misc.py +++ b/pyspi/statistics/misc.py @@ -8,6 +8,7 @@ from sklearn.metrics import mean_squared_error from sklearn import linear_model import mne.connectivity as mnec +from pyspi.lib.ids.dependence import compute_IDS from pyspi.base import ( Directed, @@ -147,7 +148,7 @@ def bivariate(self, data, i=None, j=None): class PowerEnvelopeCorrelation(Undirected, Unsigned): - humanname = "Power envelope correlation" + name = "Power envelope correlation" identifier = "pec" labels = ["unsigned", "misc", "undirected"] @@ -173,3 +174,28 @@ def multivariate(self, data): ) np.fill_diagonal(adj, np.nan) return adj + +class InterDependenceScore(Undirected, Unsigned): + name = "Interdependence score" + identifier = "ids" + labels = ["unsigned", "misc", "undirected", "nonlinear"] + + def __init__( + self, + terms=6, + pnorm='max', + bandwidth=0.5 + ): + self._num_terms = terms + self._p_norm = pnorm + self._bandwidth_term = bandwidth + + + @parse_multivariate + def multivariate(self, data): + # reshape for the compute_IDS function which expects shape (obs, proc) + z = np.squeeze(data.to_numpy(), axis=2).T + ids = compute_IDS(z, num_terms=self._num_terms, p_norm=self._p_norm, + bandwidth_term=self._bandwidth_term) + return ids + \ No newline at end of file diff --git a/tests/CML7_benchmark_tables.pkl b/tests/CML7_benchmark_tables.pkl index f23b1df..23d5a47 100644 Binary files a/tests/CML7_benchmark_tables.pkl and b/tests/CML7_benchmark_tables.pkl differ