add IDS SPI + config YAML

Author: joshuabmoore

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

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
+        

pyspi/lib/ids/__init__.py

@@ -0,0 +1,22 @@
+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).
+    
+    Taken from the implementation in: https://github.com/aradha/interdependence_scores
+    
+
+    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)

pyspi/lib/ids/dependence.py

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

pyspi/lib/ids/numpy_dependence.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
+