Add configurable random state

CONTRIBUTING.md

@@ -315,10 +315,11 @@ def discover_network(
     k_means: int = 5,
     n_shuffles: int = 200,
     n_jobs: int = -1,
+    random_state=42,
 ) -> nx.MultiDiGraph:
     """Main discovery interface."""
 
-    rng = np.random.default_rng(42)
+    rng = np.random.default_rng(random_state)
 
     # Validate method
     if method not in ["standard", "alternative", "information_lasso", "lasso", "your_method"]:

README.md

@@ -99,7 +99,8 @@ network = discover_network(
     max_lag=5,                  # Maximum time lag to consider
     alpha_forward=0.05,         # Forward selection significance
     alpha_backward=0.05,        # Backward elimination significance
-    n_shuffles=200              # Permutation test iterations
+    n_shuffles=200,             # Permutation test iterations
+    random_state=42,            # Permutation-test seed; pass None or another int for independent replicates
 )
 ```
 

causationentropy/core/discovery.py

@@ -3,6 +3,7 @@
 Email: kslote@clarkson.edu
 version = 1.1.0
 """
+
 import copy
 from typing import Dict, Tuple, Union
 
@@ -28,6 +29,7 @@ def discover_network(
     k_means: int = 5,
     n_shuffles: int = 200,
     n_jobs=-1,
+    random_state: Union[int, np.random.Generator, None] = 42,
 ) -> nx.MultiDiGraph:
     r"""
     Infer a causal graph via Optimal Causation Entropy (oCSE).
@@ -98,6 +100,13 @@ def discover_network(
         provide more accurate p-value estimates but increase computational cost.
     n_jobs : int, default=-1
         Number of parallel jobs for computation. -1 uses all available processors.
+    random_state : int, numpy.random.Generator, or None, default=42
+        Controls the permutation-test random stream. An integer seed is
+        converted with ``numpy.random.default_rng`` and is reproducible across
+        calls. The default of 42 preserves historical bit-reproducible
+        behavior. Pass a different integer or ``None`` for independent
+        replicates (``None`` draws entropy from the OS). A
+        ``numpy.random.Generator`` is used as-is and advanced in place.
 
     Returns
     -------
@@ -140,6 +149,9 @@ def discover_network(
     >>>
     >>> # Discover causal network
     >>> G = discover_network(data, max_lag=3, alpha_forward=0.01)
+    >>>
+    >>> # Independent shuffle-test replicate
+    >>> G_rep = discover_network(data, max_lag=3, random_state=0)
 
     References
     ----------
@@ -148,7 +160,7 @@ def discover_network(
 
     .. [2] Schreiber, T. Measuring information transfer. Physical Review Letters 85, 461 (2000).
     """
-    rng = np.random.default_rng(42)
+    rng = np.random.default_rng(random_state)
 
     if method not in ["standard", "alternative", "information_lasso", "lasso"]:
         raise NotImplementedError(f"discover_network: method={method} not supported.")

causationentropy/tests/test_discovery.py

@@ -8,6 +8,7 @@
 from causationentropy.core.discovery import (
     discover_network,
     lasso_optimal_causation_entropy,
+    shuffle_test,
 )
 
 
@@ -496,6 +497,108 @@ def test_pvalue_calculation_correctness(self):
             pytest.skip("No X0->X1 edges found - test setup may need adjustment")
 
 
+def _edge_signature(G):
+    """Comparable edge payload for reproducibility checks."""
+    return sorted(
+        (u, v, k, d.get("lag"), d.get("cmi"), d.get("p_value"))
+        for u, v, k, d in G.edges(keys=True, data=True)
+    )
+
+
+def _permutation_rows(X, seed, n_shuffles):
+    """Replay the row shuffles shuffle_test draws from an integer seed."""
+    rng = np.random.default_rng(seed)
+    return [X[rng.permutation(len(X)), :].copy() for _ in range(n_shuffles)]
+
+
+class TestDiscoverNetworkRandomState:
+    """Tests for user-controllable permutation-test randomness."""
+
+    def test_same_random_state_identical_results(self):
+        """Two calls with the same integer seed produce identical graphs."""
+        rng = np.random.default_rng(0)
+        data = rng.normal(size=(40, 2))
+
+        G1 = discover_network(data, max_lag=1, n_shuffles=15, random_state=7)
+        G2 = discover_network(data, max_lag=1, n_shuffles=15, random_state=7)
+
+        assert set(G1.nodes()) == set(G2.nodes())
+        assert _edge_signature(G1) == _edge_signature(G2)
+
+    def test_different_seeds_change_shuffle_sequence(self):
+        """Different seeds change the permutation sequence, not necessarily the graph."""
+        X = np.arange(16, dtype=float).reshape(16, 1)
+        Y = np.zeros((16, 1))
+        n_shuffles = 4
+        captured = []
+
+        def fake_cmi(X_perm, Y_arg, Z, **kwargs):
+            captured.append(np.asarray(X_perm).copy())
+            return 0.0
+
+        with patch(
+            "causationentropy.core.discovery.conditional_mutual_information",
+            side_effect=fake_cmi,
+        ):
+            shuffle_test(X, Y, None, 0.1, alpha=0.05, n_shuffles=n_shuffles, rng=0)
+            first = captured[:]
+            captured.clear()
+            shuffle_test(X, Y, None, 0.1, alpha=0.05, n_shuffles=n_shuffles, rng=1)
+            second = captured[:]
+
+        expected_0 = _permutation_rows(X, 0, n_shuffles)
+        expected_1 = _permutation_rows(X, 1, n_shuffles)
+        assert len(first) == n_shuffles
+        assert len(second) == n_shuffles
+        for got, expected in zip(first, expected_0):
+            np.testing.assert_array_equal(got, expected)
+        for got, expected in zip(second, expected_1):
+            np.testing.assert_array_equal(got, expected)
+        assert any(not np.array_equal(a, b) for a, b in zip(expected_0, expected_1))
+
+    def test_default_random_state_matches_seed_42(self):
+        """Omitting random_state is the same as random_state=42."""
+        rng = np.random.default_rng(1)
+        data = rng.normal(size=(40, 2))
+
+        G_default = discover_network(data, max_lag=1, n_shuffles=15)
+        G_42 = discover_network(data, max_lag=1, n_shuffles=15, random_state=42)
+        G_default_again = discover_network(data, max_lag=1, n_shuffles=15)
+
+        assert _edge_signature(G_default) == _edge_signature(G_42)
+        assert _edge_signature(G_default) == _edge_signature(G_default_again)
+
+    def test_numpy_generator_is_accepted(self):
+        """Equivalent Generators produce identical results; a Generator is consumed."""
+        rng = np.random.default_rng(2)
+        data = rng.normal(size=(40, 2))
+
+        G1 = discover_network(
+            data, max_lag=1, n_shuffles=15, random_state=np.random.default_rng(99)
+        )
+        G2 = discover_network(
+            data, max_lag=1, n_shuffles=15, random_state=np.random.default_rng(99)
+        )
+        assert _edge_signature(G1) == _edge_signature(G2)
+
+        shared = np.random.default_rng(99)
+        discover_network(data, max_lag=1, n_shuffles=10, random_state=shared)
+        next_after_one_call = int(shared.integers(2**63))
+        shared = np.random.default_rng(99)
+        discover_network(data, max_lag=1, n_shuffles=10, random_state=shared)
+        discover_network(data, max_lag=1, n_shuffles=10, random_state=shared)
+        next_after_two_calls = int(shared.integers(2**63))
+        assert next_after_one_call != next_after_two_calls
+
+    def test_random_state_none_runs(self):
+        """random_state=None is a valid independent-entropy setting."""
+        rng = np.random.default_rng(3)
+        data = rng.normal(size=(30, 2))
+        G = discover_network(data, max_lag=1, n_shuffles=10, random_state=None)
+        assert isinstance(G, nx.MultiDiGraph)
+        assert set(G.nodes()) == {"X0", "X1"}
+
+
 class TestLassoOptimalCausationEntropy:
     """Test LASSO-based variable selection for causal discovery."""