diff --git a/causationentropy/core/discovery.py b/causationentropy/core/discovery.py index b658d7f..1a6f2ae 100644 --- a/causationentropy/core/discovery.py +++ b/causationentropy/core/discovery.py @@ -5,12 +5,13 @@ """ import copy +import math from typing import Dict, Tuple, Union import networkx as nx import numpy as np import pandas as pd -from sklearn.linear_model import Lasso, LassoLarsIC +from sklearn.linear_model import Lasso, LassoCV, LassoLarsIC from causationentropy.core.information.conditional_mutual_information import ( conditional_mutual_information, @@ -31,6 +32,7 @@ def discover_network( n_jobs=-1, random_state: Union[int, np.random.Generator, None] = 42, only_return_significant: bool = True, + screen_retention: float = 0.40, ) -> nx.MultiDiGraph: r""" Infer a causal graph via Optimal Causation Entropy (oCSE). @@ -75,6 +77,8 @@ def discover_network( - 'standard': Uses initial conditioning set of lagged target variables - 'alternative': No initial conditioning set - 'information_lasso': Information-theoretic variant with LASSO regularization + - 'information_screened': Information-LASSO proposal plus conditional rescue + followed by restricted standard oCSE refinement - 'lasso': Pure LASSO-based selection information : str, default='gaussian' Information measure estimator type. Options: @@ -115,6 +119,10 @@ def discover_network( ``significant=False`` edge attribute. This is intended for small networks only (e.g. delay-analysis plots), since it runs a CMI + shuffle test for all ``n * max_lag`` candidates per target. + screen_retention : float, default=0.40 + Soft target fraction of external lagged candidates retained by the + information_screened method before restricted oCSE refinement. + Ignored by other methods. Appended to preserve existing positional calls. Returns ------- @@ -174,7 +182,13 @@ def discover_network( """ rng = np.random.default_rng(random_state) - if method not in ["standard", "alternative", "information_lasso", "lasso"]: + if method not in [ + "standard", + "alternative", + "information_lasso", + "information_screened", + "lasso", + ]: raise NotImplementedError(f"discover_network: method={method} not supported.") supported_information_types = ["gaussian", "knn", "kde", "geometric_knn", "poisson"] if information not in supported_information_types: @@ -218,13 +232,21 @@ def discover_network( print(f"Estimating edges for node {i} ({var_names[i]})") Y = Y_all[:, [i]] # shape: (T - max_lag, 1) + candidate_ids = list(range(len(feature_names))) + output_base = None + report_candidate_ids = None if method == "standard": Z_init = [] for tau in range(1, max_lag + 1): Z_init.append(series[max_lag - tau : T - tau, i]) # lagged Y_i Z_init = np.column_stack(Z_init) # shape: (T - max_lag, max_lag) - S = standard_optimal_causation_entropy( - X_lagged, + # Target-history lags are already in Z_init; do not retest them. + candidate_ids = [ + idx for idx, (source, _lag) in enumerate(feature_names) if source != i + ] + X_candidates = X_lagged[:, candidate_ids] + S_local = standard_optimal_causation_entropy( + X_candidates, Y, Z_init, rng, @@ -236,6 +258,9 @@ def discover_network( k_means, bandwidth, ) + S = [candidate_ids[idx] for idx in S_local] + output_base = Z_init + report_candidate_ids = candidate_ids if method == "alternative": S = alternative_optimal_causation_entropy( X_lagged, @@ -250,7 +275,39 @@ def discover_network( bandwidth, ) if method == "information_lasso": - S = information_lasso_optimal_causation_entropy(X_lagged, Y, rng) + S = information_lasso_optimal_causation_entropy( + X_lagged, + Y, + rng, + information=information, + metric=metric, + k_means=k_means, + bandwidth=bandwidth, + ) + if method == "information_screened": + Z_init = np.column_stack( + [series[max_lag - tau : T - tau, i] for tau in range(1, max_lag + 1)] + ) + candidate_ids = [ + idx for idx, (source, _lag) in enumerate(feature_names) if source != i + ] + screened_local = information_screened_optimal_causation_entropy( + X_lagged[:, candidate_ids], + Y, + Z_init, + rng, + retention=screen_retention, + alpha_forward=alpha_forward, + alpha_backward=alpha_backward, + n_shuffles=n_shuffles, + information=information, + metric=metric, + k_means=k_means, + bandwidth=bandwidth, + ) + S = [candidate_ids[int(local_idx)] for local_idx in screened_local] + output_base = Z_init + report_candidate_ids = candidate_ids if method == "lasso": S = lasso_optimal_causation_entropy(X_lagged, Y, rng) for s in S: @@ -260,9 +317,17 @@ def discover_network( X_predictor = X_lagged[:, [s]] # predictor at this lag Y_target = Y # target variable - # Conditioning set: all other selected predictors for this target + # Conditioning set: all other selected predictors for this target. + # Standard and screened oCSE keep the target-history baseline used + # throughout forward/backward refinement. other_selected = [idx for idx in S if idx != s] - Z_cond = X_lagged[:, other_selected] if other_selected else None + Z_selected = X_lagged[:, other_selected] if other_selected else None + if output_base is None: + Z_cond = Z_selected + elif Z_selected is None: + Z_cond = output_base + else: + Z_cond = np.hstack((output_base, Z_selected)) # Compute conditional mutual information cmi = conditional_mutual_information( @@ -308,25 +373,36 @@ def discover_network( significant=True, ) - selected_sets.append((i, Y, list(S))) + selected_sets.append((i, Y, list(S), output_base, report_candidate_ids)) # Report the remaining candidates only after every target is done, so the # extra shuffle tests below don't change the random stream used by the # forward/backward selection of later targets. if not only_return_significant: - for i, Y, S in selected_sets: + for i, Y, S, output_base, report_candidate_ids in selected_sets: # Report every tested candidate, including insignificant links, # so delay-analysis plots can show CMI versus lag. Conditioning # matches the significant edges above (selected set, minus the # candidate itself when it is selected). selected = set(S) n_features = len(feature_names) - for cand in range(n_features): + candidates_to_report = ( + report_candidate_ids + if report_candidate_ids is not None + else range(n_features) + ) + for cand in candidates_to_report: if cand in selected: continue src_var, src_lag = feature_names[cand] X_predictor = X_lagged[:, [cand]] - Z_cond = X_lagged[:, S] if S else None + Z_selected = X_lagged[:, S] if S else None + if output_base is None: + Z_cond = Z_selected + elif Z_selected is None: + Z_cond = output_base + else: + Z_cond = np.hstack((output_base, Z_selected)) cmi = conditional_mutual_information( X_predictor, @@ -431,6 +507,7 @@ def standard_optimal_causation_entropy( metric, k_means, bandwidth, + Z_init=Z_init, ) return S @@ -499,14 +576,40 @@ def alternative_optimal_causation_entropy( def information_lasso_optimal_causation_entropy( - X, Y, rng, criterion="bic", max_lambda=100, cross_val=10, information="gaussian" + X, + Y, + rng, + criterion="bic", + max_lambda=100, + cross_val=10, + information="gaussian", + metric="euclidean", + k_means=5, + bandwidth="silverman", ): - """ - Execute information-theoretic variant of oCSE with LASSO regularization. + r""" + Select predictors with information-weighted LASSO. + + Each candidate first receives a non-negative information weight from its + conditional mutual information with the target (with no additional + conditioning set in this standalone pathway): + + .. math:: + + w_j = \frac{I(X_j; Y)}{\sum_k I(X_k; Y)}. - This method combines information-theoretic causal discovery with LASSO regularization - to handle high-dimensional predictor spaces. The approach balances causal relationship - strength with model complexity. + The weighted-LASSO objective is + + .. math:: + + \min_{\boldsymbol{\beta}} + \frac{1}{2n}\|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\|_2^2 + + \lambda \sum_j \frac{1}{w_j}|\beta_j|. + + It is solved with the standard scikit-learn LASSO solvers by scaling + feature j by w_j. This is the usual reparameterization of a weighted + l1 penalty. Candidates with zero information weight receive an infinite + effective penalty and cannot enter the model. Parameters ---------- @@ -515,29 +618,166 @@ def information_lasso_optimal_causation_entropy( Y : array-like of shape (T, 1) Target variable column. rng : numpy.random.Generator - Random number generator. + Random number generator retained for API consistency. The current + weighted-LASSO solvers are deterministic. criterion : str, default='bic' - Information criterion for model selection ('bic' or 'aic'). + Information criterion used by LassoLarsIC when there are enough + samples relative to predictors. max_lambda : int, default=100 - Maximum number of LASSO iterations. + Maximum solver iterations for LassoLarsIC. LassoCV uses at least 1000 + iterations to avoid premature convergence in high dimensions. cross_val : int, default=10 - Cross-validation folds (currently unused). + Number of cross-validation folds used by LassoCV in the + high-dimensional regime. information : str, default='gaussian' - Information measure estimator type. + Information estimator passed to conditional_mutual_information. + metric : str, default='euclidean' + Distance metric used by k-NN information estimators. + k_means : int, default=5 + Number of neighbors used by k-NN information estimators. + bandwidth : str or float, default='silverman' + Bandwidth used by KDE information estimators. Returns ------- S : list of int - Indices of selected predictor variables. + Indices of predictors with non-zero weighted-LASSO coefficients. Notes ----- - This is a simplified implementation that delegates to LASSO. Future versions - will incorporate information-theoretic weighting into the regularization. + LassoLarsIC is used when the information criterion is identifiable from + the sample size. When predictors are at least as numerous as the available + samples, LassoCV uses the existing cross_val parameter to choose the + regularization strength instead of falling back to a fixed default penalty. + """ + del rng # Retained in the public signature for consistency with other methods. + + n_features = X.shape[1] + information_values = np.empty(n_features, dtype=float) + + for j in range(n_features): + value = conditional_mutual_information( + X[:, [j]], + Y, + None, + method=information, + metric=metric, + k=k_means, + bandwidth=bandwidth, + ) + if not np.isfinite(value): + raise ValueError( + "Information-LASSO received a non-finite information weight " + f"for predictor {j}." + ) + information_values[j] = max(0.0, value) + + total_information = information_values.sum() + if total_information <= 0: + return [] + + weights = information_values / total_information + X_weighted = X * weights.reshape(1, -1) + + if X.shape[0] > n_features + 1: + model = LassoLarsIC(criterion=criterion, max_iter=max_lambda).fit( + X_weighted, Y.flatten() + ) + else: + cv_folds = min(cross_val, X.shape[0]) + if cv_folds < 2: + raise ValueError( + "Information-LASSO requires at least two samples for " + "high-dimensional cross-validation." + ) + model = LassoCV(cv=cv_folds, max_iter=max(1000, max_lambda)).fit( + X_weighted, Y.flatten() + ) + + return np.flatnonzero(model.coef_ != 0).tolist() + + +def information_screened_optimal_causation_entropy( + X, + Y, + Z_init, + rng, + retention=0.40, + alpha_forward=0.05, + alpha_backward=0.05, + n_shuffles=200, + information="gaussian", + metric="euclidean", + k_means=5, + bandwidth="silverman", +): + r"""Screen with Information-LASSO plus conditional rescue, then refine. + + X contains only external candidate predictors for one target. Target + history is supplied separately through Z_init and is preserved in both the + conditional-rescue and restricted-oCSE conditioning sets. """ + if not 0 < retention <= 1: + raise ValueError("retention must be in (0, 1].") - # This is a simplified implementation - needs proper information-theoretic weighting - return lasso_optimal_causation_entropy(X, Y, rng, criterion, max_lambda, cross_val) + n_features = X.shape[1] + if n_features == 0: + return [] + + endpoint = information_lasso_optimal_causation_entropy( + X, + Y, + rng, + information=information, + metric=metric, + k_means=k_means, + bandwidth=bandwidth, + ) + endpoint = [int(j) for j in endpoint] + selected = set(endpoint) + + target_size = int(math.ceil(retention * n_features)) + if len(selected) < n_features: + target_size = max(target_size, len(selected) + 1) + target_size = min(n_features, target_size) + + if len(selected) < target_size: + Z_rescue = np.hstack((Z_init, X[:, endpoint])) if endpoint else Z_init + scored = [] + for j in range(n_features): + if j in selected: + continue + score = conditional_mutual_information( + X[:, [j]], + Y, + Z_rescue, + method=information, + metric=metric, + k=k_means, + bandwidth=bandwidth, + ) + if not np.isfinite(score): + score = -np.inf + scored.append((float(score), int(j))) + + scored.sort(key=lambda item: (-item[0], item[1])) + selected.update(j for _score, j in scored[: target_size - len(selected)]) + + screened = sorted(selected) + refined_local = standard_optimal_causation_entropy( + X[:, screened], + Y, + Z_init, + rng, + alpha1=alpha_forward, + alpha2=alpha_backward, + n_shuffles=n_shuffles, + information=information, + metric=metric, + k_means=k_means, + bandwidth=bandwidth, + ) + return [screened[int(local_idx)] for local_idx in refined_local] def lasso_optimal_causation_entropy( @@ -819,6 +1059,7 @@ def backward( metric="euclidean", k_means=5, bandwidth="silverman", + Z_init=None, ): r""" Backward elimination phase of optimal Causation Entropy. @@ -826,13 +1067,13 @@ def backward( This function performs backward elimination to remove spurious causal relationships identified during forward selection. For each predictor selected in the forward phase, it tests whether the predictor maintains statistical significance when conditioned on - all other selected predictors. + all other selected predictors and any fixed initial conditioning set. For each predictor :math:`X_j` in the selected set, the test evaluates: .. math:: - I(X_j^{(t)}; Y^{(t+\tau)} | \mathbf{S}_{-j}^{(t)}) > \text{threshold} + I(X_j^{(t)}; Y^{(t+\tau)} | \mathbf{Z}_{init}, \mathbf{S}_{-j}^{(t)}) > \text{threshold} where :math:`\mathbf{S}_{-j}^{(t)}` represents all selected predictors except :math:`X_j`. @@ -853,6 +1094,10 @@ def backward( information : str, default='gaussian' Information measure estimator type. + Z_init : array-like of shape (T, p), optional + Fixed conditioning baseline retained during every backward test. + The default None preserves unconditioned alternative-oCSE behavior. + Returns ------- S_final : list of int @@ -872,8 +1117,14 @@ def backward( S = copy.deepcopy(S_init) # working copy for j in rng.permutation(S_init): - # conditioning set Z = S \ {j} - Z = X_full[:, [k for k in S if k != j]] if len(S) > 1 else None + other_ids = [k for k in S if k != j] + Z_selected = X_full[:, other_ids] if other_ids else None + if Z_init is None: + Z = Z_selected + elif Z_selected is None: + Z = Z_init + else: + Z = np.hstack((Z_init, Z_selected)) Xj = X_full[:, [j]] cmij = conditional_mutual_information( diff --git a/causationentropy/tests/test_discovery.py b/causationentropy/tests/test_discovery.py index 936ae7b..50044fa 100644 --- a/causationentropy/tests/test_discovery.py +++ b/causationentropy/tests/test_discovery.py @@ -88,7 +88,13 @@ def test_discover_network_valid_methods(self): """Test that all valid methods are accepted.""" data = np.random.normal(0, 1, (30, 2)) - valid_methods = ["standard", "alternative", "information_lasso", "lasso"] + valid_methods = [ + "standard", + "alternative", + "information_lasso", + "information_screened", + "lasso", + ] for method in valid_methods: G = discover_network(data, method=method, max_lag=1, n_shuffles=10) assert isinstance(G, nx.MultiDiGraph) diff --git a/causationentropy/tests/test_information_lasso.py b/causationentropy/tests/test_information_lasso.py new file mode 100644 index 0000000..fbd4dbd --- /dev/null +++ b/causationentropy/tests/test_information_lasso.py @@ -0,0 +1,149 @@ +from unittest.mock import MagicMock, patch + +import numpy as np +import pytest + +from causationentropy.core.discovery import ( + discover_network, + information_lasso_optimal_causation_entropy, +) + + +def test_information_lasso_scales_features_by_information_weights(): + """CMI weights must alter the sparse-selection design matrix.""" + X = np.arange(10, dtype=float).reshape(5, 2) + 1.0 + Y = np.arange(5, dtype=float).reshape(-1, 1) + rng = np.random.default_rng(0) + + fitted = MagicMock() + fitted.fit.return_value = fitted + fitted.coef_ = np.array([1.0, 0.0]) + + with patch( + "causationentropy.core.discovery.conditional_mutual_information", + side_effect=[3.0, 1.0], + ) as cmi, patch( + "causationentropy.core.discovery.LassoLarsIC", + return_value=fitted, + ): + selected = information_lasso_optimal_causation_entropy( + X, + Y, + rng, + information="kde", + metric="cityblock", + k_means=7, + bandwidth=0.5, + ) + + assert selected == [0] + np.testing.assert_allclose( + fitted.fit.call_args.args[0], + X * np.array([0.75, 0.25]), + ) + for call in cmi.call_args_list: + assert call.args[2] is None + assert call.kwargs == { + "method": "kde", + "metric": "cityblock", + "k": 7, + "bandwidth": 0.5, + } + + +def test_information_lasso_uses_cross_validation_when_high_dimensional(): + """p >= n uses the existing cross_val parameter instead of fixed-alpha LASSO.""" + X = np.arange(20, dtype=float).reshape(4, 5) + 1.0 + Y = np.arange(4, dtype=float).reshape(-1, 1) + rng = np.random.default_rng(0) + + fitted = MagicMock() + fitted.fit.return_value = fitted + fitted.coef_ = np.array([0.0, 1.0, 0.0, 0.0, -1.0]) + + with patch( + "causationentropy.core.discovery.conditional_mutual_information", + side_effect=[5.0, 4.0, 3.0, 2.0, 1.0], + ), patch( + "causationentropy.core.discovery.LassoCV", + return_value=fitted, + ) as lasso_cv: + selected = information_lasso_optimal_causation_entropy( + X, + Y, + rng, + max_lambda=77, + cross_val=10, + ) + + assert selected == [1, 4] + lasso_cv.assert_called_once_with(cv=4, max_iter=1000) + + +def test_information_lasso_returns_empty_when_all_information_is_zero(): + """No information signal should not fall back to ordinary LASSO.""" + X = np.ones((6, 3)) + Y = np.arange(6, dtype=float).reshape(-1, 1) + + with patch( + "causationentropy.core.discovery.conditional_mutual_information", + return_value=0.0, + ), patch("causationentropy.core.discovery.LassoLarsIC") as lars, patch( + "causationentropy.core.discovery.LassoCV" + ) as lasso_cv: + selected = information_lasso_optimal_causation_entropy( + X, + Y, + np.random.default_rng(0), + ) + + assert selected == [] + lars.assert_not_called() + lasso_cv.assert_not_called() + + +def test_information_lasso_rejects_nonfinite_information_weight(): + """Estimator failures should surface rather than silently changing selection.""" + X = np.ones((6, 2)) + Y = np.arange(6, dtype=float).reshape(-1, 1) + + with patch( + "causationentropy.core.discovery.conditional_mutual_information", + side_effect=[1.0, np.nan], + ): + with pytest.raises(ValueError, match="non-finite information weight"): + information_lasso_optimal_causation_entropy( + X, + Y, + np.random.default_rng(0), + ) + + +def test_discover_network_forwards_information_lasso_estimator_options(): + """The public API estimator settings must reach Information-LASSO.""" + data = np.arange(36, dtype=float).reshape(12, 3) + + with patch( + "causationentropy.core.discovery.information_lasso_optimal_causation_entropy", + return_value=[], + ) as information_lasso: + graph = discover_network( + data, + method="information_lasso", + information="kde", + max_lag=1, + metric="cityblock", + bandwidth=0.4, + k_means=9, + n_shuffles=2, + ) + + assert graph.number_of_nodes() == 3 + assert information_lasso.call_count == 3 + for call in information_lasso.call_args_list: + assert call.kwargs == { + "information": "kde", + "metric": "cityblock", + "k_means": 9, + "bandwidth": 0.4, + } diff --git a/causationentropy/tests/test_information_screened_ocse.py b/causationentropy/tests/test_information_screened_ocse.py new file mode 100644 index 0000000..ff878ef --- /dev/null +++ b/causationentropy/tests/test_information_screened_ocse.py @@ -0,0 +1,421 @@ +import inspect +from unittest.mock import Mock + +import networkx as nx +import numpy as np +import pytest + +import causationentropy.core.discovery as discovery + + +def test_information_screened_rescue_then_maps_refined_support(monkeypatch): + X = np.column_stack( + [ + np.linspace(0.0, 1.0, 20), + np.linspace(1.0, 2.0, 20), + np.linspace(2.0, 3.0, 20), + np.linspace(3.0, 4.0, 20), + ] + ) + Y = np.linspace(-1.0, 1.0, 20).reshape(-1, 1) + Z_init = np.ones((20, 1)) + + monkeypatch.setattr( + discovery, + "information_lasso_optimal_causation_entropy", + lambda *args, **kwargs: [0], + ) + + rescue_conditioning = [] + + def fake_cmi(x, _y, z, **kwargs): + rescue_conditioning.append(z.copy()) + # Candidate 2 is the strongest excluded rescue candidate. + first = float(x[0, 0]) + return {1.0: 1.0, 2.0: 5.0, 3.0: 2.0}.get(first, 0.0) + + monkeypatch.setattr(discovery, "conditional_mutual_information", fake_cmi) + + captured = {} + + def fake_standard(x, _y, z_init, _rng, **kwargs): + captured["screened"] = x.copy() + captured["z_init"] = z_init.copy() + # Return the second column of the screened matrix. + return [1] + + monkeypatch.setattr( + discovery, + "standard_optimal_causation_entropy", + fake_standard, + ) + + selected = discovery.information_screened_optimal_causation_entropy( + X, + Y, + Z_init, + np.random.default_rng(0), + retention=0.5, + ) + + assert selected == [2] + np.testing.assert_allclose(captured["screened"], X[:, [0, 2]]) + np.testing.assert_allclose(captured["z_init"], Z_init) + + expected_rescue_conditioning = np.hstack((Z_init, X[:, [0]])) + assert len(rescue_conditioning) == 3 + for z in rescue_conditioning: + np.testing.assert_allclose(z, expected_rescue_conditioning) + + +def test_information_screened_rescue_uses_history_when_proposal_is_empty( + monkeypatch, +): + X = np.column_stack( + [ + np.linspace(0.0, 1.0, 20), + np.linspace(1.0, 2.0, 20), + ] + ) + Y = np.linspace(-1.0, 1.0, 20).reshape(-1, 1) + Z_init = np.column_stack( + [ + np.linspace(4.0, 5.0, 20), + np.linspace(5.0, 6.0, 20), + ] + ) + rescue_conditioning = [] + + monkeypatch.setattr( + discovery, + "information_lasso_optimal_causation_entropy", + lambda *args, **kwargs: [], + ) + + def fake_cmi(_x, _y, z, **kwargs): + rescue_conditioning.append(z.copy()) + return 1.0 + + monkeypatch.setattr(discovery, "conditional_mutual_information", fake_cmi) + monkeypatch.setattr( + discovery, + "standard_optimal_causation_entropy", + lambda *args, **kwargs: [], + ) + + discovery.information_screened_optimal_causation_entropy( + X, + Y, + Z_init, + np.random.default_rng(0), + retention=0.5, + ) + + assert len(rescue_conditioning) == 2 + for z in rescue_conditioning: + np.testing.assert_allclose(z, Z_init) + + +def test_information_screened_retention_is_soft_and_preserves_proposal(monkeypatch): + X = np.column_stack( + [ + np.linspace(0.0, 1.0, 24), + np.linspace(1.0, 2.0, 24), + np.linspace(2.0, 3.0, 24), + np.linspace(3.0, 4.0, 24), + np.linspace(4.0, 5.0, 24), + ] + ) + Y = np.linspace(-1.0, 1.0, 24).reshape(-1, 1) + Z_init = np.linspace(10.0, 11.0, 24).reshape(-1, 1) + + # The proposal already exceeds the nominal 40% target (3 > ceil(0.4 * 5)=2). + # Preserve the proposal and leave one rescue slot rather than truncating it. + monkeypatch.setattr( + discovery, + "information_lasso_optimal_causation_entropy", + lambda *args, **kwargs: [0, 1, 2], + ) + + rescue_conditioning = [] + + def fake_cmi(x, _y, z, **kwargs): + rescue_conditioning.append(z.copy()) + # Candidate 4 wins the single rescue slot over candidate 3. + return 2.0 if np.isclose(x[0, 0], 4.0) else 1.0 + + monkeypatch.setattr(discovery, "conditional_mutual_information", fake_cmi) + + captured = {} + + def fake_standard(x, _y, z_init, _rng, **kwargs): + captured["screened"] = x.copy() + captured["z_init"] = z_init.copy() + return [] + + monkeypatch.setattr( + discovery, + "standard_optimal_causation_entropy", + fake_standard, + ) + + selected = discovery.information_screened_optimal_causation_entropy( + X, + Y, + Z_init, + np.random.default_rng(0), + retention=0.40, + ) + + assert selected == [] + np.testing.assert_allclose(captured["screened"], X[:, [0, 1, 2, 4]]) + np.testing.assert_allclose(captured["z_init"], Z_init) + + expected_conditioning = np.hstack((Z_init, X[:, [0, 1, 2]])) + assert len(rescue_conditioning) == 2 + for z in rescue_conditioning: + np.testing.assert_allclose(z, expected_conditioning) + + +def test_information_screened_rejects_invalid_retention(): + X = np.ones((20, 2)) + Y = np.ones((20, 1)) + Z = np.ones((20, 1)) + + for value in (0.0, -0.1, 1.1): + with pytest.raises(ValueError, match="retention"): + discovery.information_screened_optimal_causation_entropy( + X, + Y, + Z, + np.random.default_rng(0), + retention=value, + ) + + +def test_discover_information_screened_uses_external_candidates(monkeypatch): + captured = [] + + def fake_screened(X, Y, Z_init, rng, **kwargs): + captured.append((X.shape, Y.shape, Z_init.shape, kwargs["retention"])) + return [] + + monkeypatch.setattr( + discovery, + "information_screened_optimal_causation_entropy", + fake_screened, + ) + + data = np.random.default_rng(0).normal(size=(40, 3)) + graph = discovery.discover_network( + data, + method="information_screened", + max_lag=2, + screen_retention=0.35, + n_shuffles=3, + random_state=0, + ) + + assert len(graph) == 3 + assert len(captured) == 3 + for x_shape, y_shape, z_shape, retention in captured: + assert x_shape == (38, 4) # (n-1) * max_lag external candidates + assert y_shape == (38, 1) + assert z_shape == (38, 2) + assert retention == 0.35 + + +def test_backward_preserves_initial_conditioning_for_screened_refinement(monkeypatch): + rng = np.random.default_rng(7) + X = rng.normal(size=(50, 3)) + Y = rng.normal(size=(50, 1)) + Z_init = rng.normal(size=(50, 2)) + seen = [] + + def fake_cmi(_x, _y, z, **kwargs): + seen.append(z.copy()) + return 1.0 + + monkeypatch.setattr(discovery, "conditional_mutual_information", fake_cmi) + monkeypatch.setattr( + discovery, + "shuffle_test", + lambda *args, **kwargs: { + "Threshold": 0.0, + "Value": 1.0, + "Pass": True, + "P_value": 0.0, + }, + ) + + selected = discovery.backward( + X, + Y, + [0, 1], + np.random.default_rng(11), + n_shuffles=3, + Z_init=Z_init, + ) + + assert selected == [0, 1] + assert len(seen) == 2 + for z in seen: + np.testing.assert_allclose(z[:, :2], Z_init) + + +def test_standard_discovery_does_not_retest_target_history(monkeypatch): + calls = [] + + def fake_standard(X, Y, Z_init, rng, *args, **kwargs): + calls.append((X.copy(), Z_init.copy())) + return [] + + monkeypatch.setattr( + discovery, + "standard_optimal_causation_entropy", + fake_standard, + ) + + data = np.arange(120, dtype=float).reshape(40, 3) + discovery.discover_network( + data, + method="standard", + max_lag=2, + n_shuffles=3, + ) + + assert len(calls) == 3 + for X_candidates, Z_init in calls: + assert X_candidates.shape == (38, 4) + assert Z_init.shape == (38, 2) + for candidate in X_candidates.T: + assert not any( + np.array_equal(candidate, conditioned) for conditioned in Z_init.T + ) + + +def test_discover_network_preserves_legacy_positional_arguments(): + """Old calls must still bind n_jobs, random_state and reporting correctly.""" + legacy = [ + "data", + "method", + "information", + "max_lag", + "alpha_forward", + "alpha_backward", + "metric", + "bandwidth", + "k_means", + "n_shuffles", + "n_jobs", + "random_state", + "only_return_significant", + ] + assert list(inspect.signature(discovery.discover_network).parameters) == ( + legacy + ["screen_retention"] + ) + data = np.random.default_rng(17).normal(size=(40, 3)) + positional = discovery.discover_network( + data, + "standard", + "gaussian", + 1, + 0.05, + 0.05, + None, + "silverman", + 5, + 5, + 1, + 17, + False, + ) + named = discovery.discover_network( + data, + method="standard", + information="gaussian", + max_lag=1, + n_shuffles=5, + n_jobs=1, + random_state=17, + only_return_significant=False, + ) + assert nx.utils.graphs_equal(positional, named) + assert positional.number_of_edges() == 6 + + +@pytest.mark.parametrize("method", ["standard", "information_screened"]) +@pytest.mark.parametrize("only_significant", [True, False]) +def test_reporting_preserves_each_targets_history( + monkeypatch, method, only_significant +): + """Deferred reporting must not reuse the final target's baseline or mask.""" + data = np.random.default_rng(33).normal(size=(40, 3)) + calls = [] + + def check_conditioning(x, y, z, method, metric, k, bandwidth): + target = next(i for i in range(3) if np.array_equal(y[:, 0], data[2:, i])) + expected = np.column_stack((data[1:-1, target], data[:-2, target])) + np.testing.assert_array_equal(z[:, :2], expected) + assert not any(np.array_equal(x[:, 0], col) for col in expected.T) + calls.append(target) + return 0.5 + + monkeypatch.setattr( + discovery, "standard_optimal_causation_entropy", Mock(return_value=[0]) + ) + monkeypatch.setattr( + discovery, + "information_screened_optimal_causation_entropy", + Mock(return_value=[0]), + ) + monkeypatch.setattr(discovery, "conditional_mutual_information", check_conditioning) + monkeypatch.setattr(discovery, "shuffle_test", Mock(return_value={"P_value": 0.01})) + graph = discovery.discover_network( + data, + method=method, + max_lag=2, + only_return_significant=only_significant, + ) + per_target = 1 if only_significant else 4 + assert graph.number_of_edges() == 3 * per_target + assert all(calls.count(target) == per_target for target in range(3)) + assert all(source != target for source, target in graph.edges()) + + +def test_standard_passes_history_to_backward(monkeypatch): + """The standard method, including restricted refinement, forwards its baseline.""" + data = np.random.default_rng(0).normal(size=(40, 4)) + backward = Mock(return_value=[0]) + monkeypatch.setattr(discovery, "standard_forward", Mock(return_value=[0])) + monkeypatch.setattr(discovery, "backward", backward) + discovery.standard_optimal_causation_entropy( + data[:, :2], + data[:, [2]], + data[:, [3]], + np.random.default_rng(1), + ) + np.testing.assert_array_equal(backward.call_args.kwargs["Z_init"], data[:, [3]]) + + +@pytest.mark.parametrize("method", ["standard", "information_screened"]) +def test_reporting_mode_keeps_selected_edges_identical(method): + """Compose the history correction with deferred reporting's RNG guarantee.""" + rng = np.random.default_rng(12) + data = rng.normal(size=(120, 3)) + data[1:, 1] = 0.9 * data[:-1, 0] + 0.1 * rng.normal(size=119) + selected = discovery.discover_network(data, method=method, max_lag=1, n_shuffles=19) + all_edges = discovery.discover_network( + data, + method=method, + max_lag=1, + n_shuffles=19, + only_return_significant=False, + ) + retained = nx.MultiDiGraph() + retained.add_nodes_from(all_edges.nodes(data=True)) + for source, target, attrs in all_edges.edges(data=True): + if attrs["significant"]: + attrs = {key: value for key, value in attrs.items() if key != "significant"} + retained.add_edge(source, target, **attrs) + assert nx.utils.graphs_equal(selected, retained) diff --git a/docs/source/theory/index.rst b/docs/source/theory/index.rst index 65dad02..195a60c 100644 --- a/docs/source/theory/index.rst +++ b/docs/source/theory/index.rst @@ -12,6 +12,7 @@ including mathematical foundations, algorithmic details, and methodological comp methods/standard_oce methods/alternative_oce methods/information_lasso + methods/information_screened_ocse methods/lasso_methods information_theory statistical_foundations @@ -65,4 +66,5 @@ Next Steps - :doc:`methods/standard_oce`: Detailed explanation of the standard oCSE algorithm - :doc:`methods/alternative_oce`: Alternative oCSE formulation without initial conditioning - :doc:`methods/information_lasso`: Information-theoretic variant with LASSO regularization -- :doc:`methods/lasso_methods`: Pure LASSO-based approaches for comparison \ No newline at end of file +- :doc:`methods/information_screened_ocse`: Conservative information screening followed by restricted oCSE +- :doc:`methods/lasso_methods`: Pure LASSO-based approaches for comparison diff --git a/docs/source/theory/methods/information_lasso.rst b/docs/source/theory/methods/information_lasso.rst index d335930..5a47851 100644 --- a/docs/source/theory/methods/information_lasso.rst +++ b/docs/source/theory/methods/information_lasso.rst @@ -44,18 +44,30 @@ where: Information-Theoretic Weights ============================= -The weights :math:`w_j` are derived from conditional mutual information measures: +The broader info-LASSO family can derive :math:`w_j` from conditional +mutual information. The standalone Path A implementation in this package is +the empty-conditioning specialization of that family. -Base Weights ------------ +General Conditional Form +------------------------ -For each potential predictor :math:`X_j^{(t-\tau)}`: +For each potential predictor :math:`X_j^{(t-\tau)}`, a general conditioned +weight can be written as: .. math:: w_{j,\tau} = \frac{I(X_j^{(t-\tau)}; X_i^{(t)} | \mathbf{Z}_i)}{\sum_{k,\tau'} I(X_k^{(t-\tau')}; X_i^{(t)} | \mathbf{Z}_i)} -This normalizes the conditional mutual information values to create relative importance weights. +For the current standalone ``method="information_lasso"`` pathway, +:math:`\mathbf{Z}_i = \varnothing`. Therefore + +.. math:: + + I(X_j; Y | \varnothing) = I(X_j; Y), + +and the implemented base weights are marginal-information weights. Conditioned +screening with a non-empty target-history set is handled separately by the +screened-oCSE pathway rather than by Path A. Adaptive Weighting ------------------ @@ -175,6 +187,40 @@ Two-Stage Implementation selected = np.where(np.abs(beta_original) > 1e-6)[0] return selected, beta_original +Current Standalone Implementation +--------------------------------- + +The public method="information_lasso" pathway implements the two-stage +weighted-LASSO approach above. + +For each lagged candidate, the standalone pathway computes + +.. math:: + + w_j = \frac{I(X_j;Y | \varnothing)}{\sum_k I(X_k;Y | \varnothing)} + = \frac{I(X_j;Y)}{\sum_k I(X_k;Y)}. + +In code, the empty conditioning set is represented by ``Z=None``. This is an +intentional scope choice for standalone Path A, not an implicit use of the +target-history conditioning set. The normalized information weights are then +used as predictor-specific penalty weights. Equivalently, scaling column +j of the design matrix by w_j lets the existing LASSO solvers optimize + +.. math:: + + \frac{1}{2n}\|\mathbf{y}-\mathbf{X}\boldsymbol{\beta}\|_2^2 + + \lambda \sum_j \frac{1}{w_j}|\beta_j|. + +Candidates with zero information receive an infinite effective penalty and are +not selected. When there are enough samples relative to predictors, +LassoLarsIC selects the regularization strength using AIC or BIC. In the +high-dimensional regime, LassoCV uses the existing cross_val parameter. + +The information estimator selected through discover_network (for example +Gaussian, KDE, or k-NN) is also used to construct these weights. Iterative +reweighting and the separate Information-LASSO screening plus exact oCSE +refinement pathway are not part of this standalone implementation. + Adaptive Implementation ---------------------- diff --git a/docs/source/theory/methods/information_screened_ocse.rst b/docs/source/theory/methods/information_screened_ocse.rst new file mode 100644 index 0000000..46f026f --- /dev/null +++ b/docs/source/theory/methods/information_screened_ocse.rst @@ -0,0 +1,146 @@ +========================== +Information-Screened oCSE +========================== + +``method="information_screened"`` uses Information-LASSO as a proposal stage, +adds a conditional-information rescue pass, and then runs standard oCSE only +inside the retained candidate set. + +The screen is deliberately conservative. Its purpose is not to estimate the +final parent set directly; it is to reduce the candidate universe while +preserving variables needed by the downstream oCSE refinement. + +Method +====== + +For a target variable :math:`Y`, let: + +- :math:`C` be the external lagged-candidate set; +- :math:`H_Y` be the lagged target-history conditioning set; +- :math:`S_0` be the Information-LASSO proposal; +- :math:`\widehat{W}` be the retained working set. + +The proposal stage first computes + +.. math:: + + S_0 = \operatorname{InfoLasso}(C \rightarrow Y). + +Candidates outside :math:`S_0` are then ranked by conditional mutual +information + +.. math:: + + r_j = \widehat{I}(X_j; Y \mid H_Y \cup S_0). + +The target of this rescue pass is a parent-superset working set rather than an +exact graph estimate. In particular, the screening objective is + +.. math:: + + \Pr\!\left(\operatorname{Pa}(Y) \subseteq \widehat{W}\right) \rightarrow 1, + +while keeping :math:`|\widehat{W}|` substantially smaller than :math:`|C|`. +The final graph decision remains the responsibility of standard oCSE: + +.. math:: + + \widehat{\operatorname{Pa}}(Y) + = \operatorname{oCSE}(Y; \widehat{W}, H_Y). + +Retention Is a Soft Target +========================== + +``screen_retention`` is a target fraction, not a hard cap. Let +:math:`p = |C|` and + +.. math:: + + B_0 = \left\lceil \text{screen\_retention} \cdot p \right\rceil. + +The implementation never truncates the Information-LASSO proposal. When there +are candidates outside :math:`S_0`, it also keeps at least one rescue slot so +that the screened pathway cannot silently collapse back to the proposal-only +path. The realized target size is therefore + +.. math:: + + B = \min\!\left(p,\; \max(B_0, |S_0| + 1)\right) + +when :math:`|S_0| < p`, and :math:`B=p` otherwise. + +This asymmetry is intentional. A false positive at the screening stage only +increases downstream computation, whereas a false negative removes a candidate +that restricted oCSE can no longer recover. Consequently, preserving the +proposal and leaving a rescue channel open takes priority over enforcing an +exact candidate count. + +The resulting working set is + +.. math:: + + \widehat{W} + = S_0 \cup \operatorname{Top}_{B-|S_0|}\{r_j : j \notin S_0\}. + +Interpretation +============== + +A useful way to view the method is as a two-level search: + +.. code-block:: none + + all external lagged candidates + | + v + Information-LASSO proposal + | + v + conditional rescue given target history + proposal + | + v + conservative working set + | + v + standard oCSE refinement + | + v + causal graph + +The screening stage therefore trades a small amount of extra candidate +retention for protection against proposal-stage misses. The downstream oCSE +stage still performs the significance-based causal selection. + +Current Scope +============= + +Standard-path conditioning is corrected independently of screening; see +:doc:`standard_oce`. Both full and restricted searches must use that same +baseline for a meaningful comparison. Existing positional arguments retain +their order; pass the new retention control by name. + +Evaluation distinguishes three quantities: + +* True-parent recall: the fraction of generating external parent-lag pairs in + the screen. Parent-complete target rate separately measures whether every + generating parent of a target survives, excluding targets with no parents. +* Full-oCSE-support coverage: the fraction of the finite-sample full-oCSE + selection in the screen. This is an agreement diagnostic, not ground truth. +* Final graph precision/recall/F1 after restricted refinement, plus total CMI + evaluations and runtime including screening. + +Measure these over sample sizes and realized retention fractions. A +proposal-only ablation removes rescue but changes candidate count; a +budget-matched marginal-fill ablation retains the proposal and fills to the +same size without conditioning on the proposal. Report both to distinguish +the value of conditional rescue from simply keeping more candidates. + +The current implementation does not claim a finite-sample sure-screening +theorem and does not guarantee that every causal parent is retained for every +finite dataset. The relevant empirical quantity is therefore the +retention/coverage frontier: how much the candidate universe can be reduced +before parent or full-oCSE-support coverage begins to degrade. + +Likewise, the rescue score is univariate conditional CMI for each excluded +candidate given :math:`H_Y \cup S_0`. Higher-order synergy searches are outside +the present scope because they would increase screening cost and change the +computational objective of this pathway. diff --git a/docs/source/theory/methods/standard_oce.rst b/docs/source/theory/methods/standard_oce.rst index d516f49..dd0b0a0 100644 --- a/docs/source/theory/methods/standard_oce.rst +++ b/docs/source/theory/methods/standard_oce.rst @@ -7,6 +7,32 @@ of the causation entropy framework. This method begins with an initial condition typically consisting of lagged values of the target variable, and systematically builds the causal predictor set through forward selection and backward elimination phases. +Conditioning Consistency Correction +================================== + +The standard-path changes accompanying information-screened oCSE are a separate +correctness correction, not a screening operation. Both the full and restricted +standard-oCSE searches use the following contract: + +* Target-history lags are fixed background variables, not selectable candidates. +* Forward selection, backward elimination, and reported edge CMI/p-values all + preserve that same target-history baseline. +* With ``only_return_significant=False``, all external candidates are reported; + target self-lags remain excluded. Background arrays and candidate indices are + saved separately for each target before deferred edge reporting. + +Previously, standard discovery could retest a lag already present in its +conditioning set, then drop the initial history during backward elimination +and edge reporting. Correcting this can change selected edges, CMI values, +p-values, self-lag output, and seeded results. It is an intentional observable +correction to the default method, not a claim of numerical backward compatibility. +The alternative method still starts without a history baseline. + +``screen_retention`` is appended after the existing ``discover_network`` +parameters, preserving their positional binding. Comparisons of full versus +screened oCSE must use the same corrected standard-path implementation; otherwise +conditioning corrections could be mistaken for screening improvements. + Mathematical Foundation ======================= @@ -333,4 +359,4 @@ rigorous statistical testing. However, users should be aware of its computationa requirements and the assumptions inherent in the initial conditioning approach. The method is particularly well-suited for time series with strong temporal dependencies -where controlling for autoregressive effects is crucial for accurate causal inference. \ No newline at end of file +where controlling for autoregressive effects is crucial for accurate causal inference.