diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index fc4a7c0..3fbaf4c 100644 --- a/CONTRIBUTING.md +++ b/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"]: diff --git a/README.md b/README.md index 779f8d8..3c98337 100644 --- a/README.md +++ b/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 ) ``` diff --git a/causationentropy/core/discovery.py b/causationentropy/core/discovery.py index ae5d442..bfe842a 100644 --- a/causationentropy/core/discovery.py +++ b/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.") diff --git a/causationentropy/tests/test_discovery.py b/causationentropy/tests/test_discovery.py index d681a99..c2426d5 100644 --- a/causationentropy/tests/test_discovery.py +++ b/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."""