updated clingo function with new parameters

gunfolds/estimation/linear_model.py

@@ -5,6 +5,8 @@
 import numpy as np
 from progressbar import ProgressBar, Percentage
 from scipy import linalg, optimize
+import scipy.sparse as sp
+from scipy.sparse.linalg import eigs
 from statsmodels.tsa.api import VAR
 from sympy.matrices import SparseMatrix
 
@@ -418,11 +420,11 @@ def init():
         if distribution == 'flat':
             x = np.ones(len(edges[0]))
         elif distribution == 'flatsigned':
-            x = np.sign(np.randn(len(edges[0])))*np.ones(len(edges[0]))
+            x = np.sign(np.random.randn(len(edges[0]))) * np.ones(len(edges[0]))
         elif distribution == 'beta':
             x = np.random.beta(0.5, 0.5, len(edges[0]))*3-1.5
         elif distribution == 'normal':
-            x = np.randn(len(edges[0]))
+            x = np.random.randn(len(edges[0]))
         elif distribution == 'uniform':
             x = np.sign(np.randn(len(edges[0])))*np.rand(len(edges[0]))
         else:
@@ -1083,4 +1085,92 @@ def randomSVARs(n, repeats=100, rate=2, density=0.1, th=0.09,
             'bidirected': B
             }
 
+def check_matrix_powers(W, powers, threshold):
+    """
+     Check if the powers of a matrix W preserve a threshold for non-zero entries.
+
+     :param W: The input matrix.
+     :type W: array_like
+
+     :param powers: List of powers to check.
+     :type powers: iterable
+
+     :param threshold: Threshold for non-zero entries.
+     :type threshold: float
+
+     :return: True if all powers of the matrix preserve the threshold for non-zero entries, False otherwise.
+     :rtype: bool
+     """
+    for n in powers:
+        W_n = np.linalg.matrix_power(W, n)
+        non_zero_indices = np.nonzero(W_n)
+        if (np.abs(W_n[non_zero_indices]) < threshold).any():
+            return False
+    return True
+
+
+def create_stable_weighted_matrix(
+    A,
+    threshold=0.1,
+    powers=[1, 2, 3, 4],
+    max_attempts=1000,
+    damping_factor=0.99,
+    random_state=None,
+):
+    """
+    Create a stable weighted matrix with a specified spectral radius.
+
+    :param A: The input matrix.
+    :type A: array_like
+
+    :param threshold: Threshold for non-zero entries preservation. Default is 0.1.
+    :type threshold: float, optional
+
+    :param powers: List of powers to check in the stability condition. Default is [1, 2, 3, 4].
+    :type powers: iterable, optional
+
+    :param max_attempts: Maximum attempts to create a stable matrix. Default is 1000.
+    :type max_attempts: int, optional
+
+    :param damping_factor: Damping factor for scaling the matrix. Default is 0.99.
+    :type damping_factor: float, optional
+
+    :param random_state: Random seed for reproducibility. Default is None.
+    :type random_state: int or None, optional
+
+    :return: A stable weighted matrix.
+    :rtype: array_like
+
+    :raises ValueError: If unable to create a stable matrix after the maximum attempts.
+    """
+    np.random.seed(
+        random_state
+    )  # Set random seed for reproducibility if provided
+    attempts = 0
+
+    while attempts < max_attempts:
+        # Generate a random matrix with the same sparsity pattern as A
+        random_weights = np.random.randn(*A.shape)
+        weighted_matrix = A * random_weights
+
+        # Convert to sparse format for efficient eigenvalue computation
+        weighted_sparse = sp.csr_matrix(weighted_matrix)
+
+        # Compute the largest eigenvalue in magnitude
+        eigenvalues, _ = eigs(weighted_sparse, k=1, which="LM")
+        max_eigenvalue = np.abs(eigenvalues[0])
+
+        # Scale the matrix so that the spectral radius is slightly less than 1
+        if max_eigenvalue > 0:
+            weighted_matrix *= damping_factor / max_eigenvalue
+            # Check if the powers of the matrix preserve the threshold for non-zero entries of A
+            if check_matrix_powers(weighted_matrix, powers, threshold):
+                return weighted_matrix
+
+        attempts += 1
+
+    raise ValueError(
+        f"Unable to create a matrix satisfying the condition after {max_attempts} attempts."
+    )
+
 # option #3

gunfolds/estimation/modifiedgc.py

@@ -0,0 +1,92 @@
+from gunfolds.utils import graphkit as gk
+from gunfolds import conversions as conv
+from itertools import combinations
+import numpy as np
+import statsmodels.api as sm
+import copy
+
+# Written by John Cook
+# https://github.com/neuroneural/Undersampled_Graph_Estimation/blob/master/dbnestimation/Tools/grangercausality.py
+# Modified to return edge weights by Sergey Plis
+
+
+def gc(data, pval=0.05, bidirected=True):
+    """
+    :param data: time series data
+    :type data: numpy matrix
+
+    :param pval: p-value for statistical significance
+    :type pval: float
+
+    :param bidirected: flag to indicate bidirectionality
+    :type bidirected: boolean
+
+    :returns: modified graph and matrices D and B
+    :rtype: tuple
+    """
+    n = data.shape[0]
+    # Initialize D and B matrices with zeros
+    D = np.zeros((n, n))
+    B = np.zeros((n, n))
+
+    data = np.asarray(np.r_[data[:, :-1], data[:, 1:]])
+
+    def FastPCtest(y, x):
+        y = y - 1
+        x = x - 1
+        nodelist = list(range(n))
+        nodelist.remove(x)
+        nodelist.remove(y)
+        yd = data[n + int(y), :].T
+        Xalt = data[:n]
+        Xnull = np.vstack([Xalt, data[n + x, :]]).T
+        Xalt = Xalt.T
+        Xnull = sm.add_constant(Xnull)
+        Xalt = sm.add_constant(Xalt)
+        estnull = sm.OLS(yd, Xnull).fit()
+        estalt = sm.OLS(yd, Xalt).fit()
+        diff = sm.stats.anova_lm(estalt, estnull)
+        # Return F-statistic and condition with the correct indices
+        Fval = diff.iloc[1, 4]
+        return Fval, diff.iloc[1, 5] > pval
+
+    def grangertest(y, x):
+        y = data[n + int(y) - 1, :].T
+        Xnull = data[:n, :].T
+        Xalt = np.vstack([data[: (x - 1), :], data[x:n, :]]).T
+        Xnull = sm.add_constant(Xnull)
+        Xalt = sm.add_constant(Xalt)
+        estnull = sm.OLS(y, Xnull).fit()
+        estalt = sm.OLS(y, Xalt).fit()
+        diff = sm.stats.anova_lm(estalt, estnull)
+        # Return F-statistic and condition with the correct indices
+        Fval = diff.iloc[1, 4]  # F statistic comparing to previous model
+        return Fval, diff.iloc[1, 5] > pval  # PR(>F) is the p-value
+
+    # Assuming gk.superclique and conv.ian2g are defined elsewhere
+    new_g = gk.superclique(n)
+    g = conv.ian2g(new_g)
+
+    for i in range(1, n + 1):
+        for j in range(1, n + 1):
+            Fval, condition = grangertest(j, i)
+            if condition:
+                sett = copy.deepcopy(g[str(i)][str(j)])
+                sett.remove((0, 1))
+                g[str(i)][str(j)] = sett
+            else:
+                D[i - 1, j - 1] = Fval
+
+    biedges = combinations(range(1, n + 1), 2)
+    for i, j in biedges:
+        Fval, condition = FastPCtest(j, i)
+        if bidirected and condition:
+            sett = copy.deepcopy(g[str(i)][str(j)])
+            sett.remove((2, 0))
+            g[str(i)][str(j)] = sett
+            g[str(j)][str(i)] = sett
+        else:
+            B[i - 1, j - 1] = Fval
+            B[j - 1, i - 1] = Fval
+
+    return conv.dict_format_converter(g), D, B