Add LUCID: regime-adaptive deconfounding for latent confounders

.github/linters/.codespellrc

@@ -6,10 +6,12 @@
 # ignore-words-list: legitimate domain acronyms codespell mistakes for typos
 # (GES = Greedy Equivalence Search; FPR = false-positive rate) and identifiers
 # that intentionally misspell a Python reserved word to use it as a name
-# (lamda = lambda, a regularization-coefficient parameter name in sigkci_gpu.py).
+# (lamda = lambda, a regularization-coefficient parameter name in sigkci_gpu.py),
+# plus "retuned" -- a real word ("re-tuned") that codespell mistakes for
+# "returned" (routed_deconf.py's research-history docstring).
 [codespell]
 # docs/_build is generated Sphinx output (HTML + a minified search index) and
 # is gitignored, so CI never scans it -- but a local run over the repo root
 # does, and reports spurious hits from the minified JS.
 skip = *.ipynb,./docs/_build
-ignore-words-list = ges,fpr,lamda
+ignore-words-list = ges,fpr,lamda,retuned

CHANGELOG.md

@@ -11,6 +11,15 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 
 ### Added
 
+* **LUCID** (`causalts.confounders`) — regime-adaptive deconfounding for causal discovery
+  under latent confounders. `run_lucid(df, max_lag)` diagnoses whether latent confounding
+  is sparse or pervasive from the residual spectrum (against a Marchenko–Pastur no-factor
+  null) and applies the matching correction, returning a `LucidResult`. Every
+  `CausalResult` now also exposes `.deconfound()`, `.tetrad_filter()`, and `.pds_filter()`
+  so LUCID (or a fixed-strategy comparator) can be applied to a graph already discovered
+  with any algorithm, without re-running the skeleton search. See the new
+  [Unobserved Confounders (LUCID)](examples/latent_confounder_detection) tutorial and the
+  `causalts.confounders` API page.
 * `corrplot(..., diag="glyph")` — renders the diagonal as an ordinary cell,
   using the same `method` and colormap as the rest of the matrix. Intended for
   *directed* matrices (a cause→effect adjacency or an edge-stability matrix),

causalts/__init__.py

@@ -37,17 +37,28 @@
 
 if TYPE_CHECKING:  # pragma: no cover - import for type checkers only
     from .cedar.legacy import SYPI  # noqa: F401
+    from .confounders import (  # noqa: F401
+        LucidResult,
+        deconfound,
+        routed_deconfound,
+        run_lucid,
+    )
     from .grace.gated_discovery import (  # noqa: F401
         run_cdnots_gated,
         run_stability_selection,
     )
     from .grace.result import GraceResult  # noqa: F401
 
 # Served on demand so that ``import causalts`` stays cheap: cedar.legacy pulls
-# in dcor + statsmodels, grace pulls in pytorch-lightning, and neither is on the
-# common discovery path.  Maps attribute name -> module that defines it.
+# in dcor + statsmodels, grace pulls in pytorch-lightning, confounders pulls in
+# statsmodels + scikit-learn's cluster/linear_model, and none is on the common
+# discovery path.  Maps attribute name -> module that defines it.
 _LAZY_ATTRS = {
     "SYPI": "causalts.cedar.legacy",
+    "LucidResult": "causalts.confounders",
+    "deconfound": "causalts.confounders",
+    "routed_deconfound": "causalts.confounders",
+    "run_lucid": "causalts.confounders",
     "run_cdnots_gated": "causalts.grace.gated_discovery",
     "run_stability_selection": "causalts.grace.gated_discovery",
     "GraceResult": "causalts.grace.result",

causalts/confounders/__init__.py

@@ -0,0 +1,39 @@
+# Copyright 2025 Bloomberg Finance L.P.
+# SPDX-License-Identifier: GPL-3.0-or-later
+"""LUCID: regime-adaptive deconfounding for time-series causal discovery.
+
+Latent confounders leave different statistical fingerprints depending on their
+structure, and no single correction handles all of them. LUCID infers the confounding
+regime from the data with a Marchenko-Pastur spectral router, then applies the strategy
+matched to that regime.
+
+The main entry point is :func:`run_lucid`, which returns a
+:class:`~causalts.confounders.result.LucidResult` carrying the graph alongside the
+router's diagnostics. :func:`routed_deconfound` is the lower-level array-returning form.
+Every result object also exposes ``.deconfound()``, which applies LUCID to an
+already-discovered graph without re-running the skeleton search.
+
+See :mod:`causalts.confounders.routed_deconf` for details.
+"""
+
+from .result import LucidResult
+from .routed_deconf import (
+    deconfound,
+    mp_factor_count,
+    pds_filter,
+    routed_deconfound,
+    run_lucid,
+    spectral_gap,
+    tetrad_filter,
+)
+
+__all__ = [
+    "run_lucid",
+    "LucidResult",
+    "routed_deconfound",
+    "deconfound",
+    "spectral_gap",
+    "mp_factor_count",
+    "pds_filter",
+    "tetrad_filter",
+]

causalts/confounders/result.py

@@ -0,0 +1,109 @@
+# Copyright 2025 Bloomberg Finance L.P.
+# SPDX-License-Identifier: GPL-3.0-or-later
+"""Result object for LUCID."""
+
+from __future__ import annotations
+
+import numpy as np
+
+from ..result import CausalResult
+
+
+class LucidResult(CausalResult):
+    """Result from :func:`~causalts.confounders.run_lucid`.
+
+    Inherits plotting and DoWhy bridge methods from :class:`CausalResult`.
+
+    Attributes
+    ----------
+    cg_tig : np.ndarray, shape (d, d, max_lag+1)
+        Binary adjacency. ``cg_tig[cause, effect, lag] == 1`` means
+        cause(t-lag) -> effect(t).
+    var_names : list[str]
+        Variable names.
+    regime : str
+        Confounding regime chosen by the router: ``"sparse"``, ``"pervasive"``, or
+        ``"sf"`` (an intermediate hub/scale-free regime).
+    spectral_ratio : float or None
+        The router statistic ``R`` -- the share of residual correlation mass carried by
+        the top ``router_k`` eigenvalues. ``None`` for routers that do not compute it.
+    tau : float or None
+        Routing threshold ``R`` was compared against. ``R > tau`` selects the pervasive
+        branch.
+    router : str
+        Which router ran (``"auto"``, ``"spectral"`` or ``"mp"``).
+    n_factors : int or None
+        Number of pervasive latent factors implied by the Marchenko-Pastur edge
+        (:func:`~causalts.confounders.mp_factor_count`). ``0`` means no pervasive
+        factor was detected.
+    factor_loadings : np.ndarray or None, shape (n_factors, d)
+        Leading right singular vectors of the VAR residuals, one row per detected
+        factor, columns aligned with ``var_names``.
+
+        .. warning::
+           These are **descriptive, not identified**. Factor directions are recovered
+           only up to rotation, so a large entry means "this variable carries dominant
+           shared variation", *not* "this variable has an identified latent parent".
+           LUCID does not identify latent variables; it corrects for their footprint.
+    info : dict
+        Full router/pipeline diagnostics as returned by
+        :func:`~causalts.confounders.routed_deconfound` with ``return_info=True``.
+    runtime : float or None
+        Wall-clock seconds for the LUCID call.
+    """
+
+    def __init__(
+        self,
+        graph: np.ndarray,
+        df,
+        var_names: list[str],
+        *,
+        info: dict | None = None,
+        n_factors: int | None = None,
+        factor_loadings: np.ndarray | None = None,
+        runtime: float | None = None,
+    ):
+        self.cg_tig = graph
+        self.var_names = list(var_names)
+        self._df = df
+        self._scm_cache = {}
+        self.info = dict(info or {})
+        self.regime = self.info.get("regime")
+        self.spectral_ratio = self.info.get("R")
+        self.tau = self.info.get("tau")
+        self.router = self.info.get("router")
+        self.n_factors = n_factors
+        self.factor_loadings = factor_loadings
+        self.runtime = runtime
+
+    def plot(self, **kwargs):
+        """Plot the deconfounded causal graph."""
+        kwargs.setdefault("show_colorbar", False)
+        return super().plot(**kwargs)
+
+    def top_factor_variables(self, factor: int = 0, n: int = 5):
+        """Variables loading most strongly on one factor direction, by ``|loading|``.
+
+        Descriptive only -- see the ``factor_loadings`` warning above.
+
+        Returns a list of ``(var_name, loading)`` pairs, largest ``|loading|`` first.
+        """
+        if self.factor_loadings is None:
+            return []
+        if not 0 <= factor < self.factor_loadings.shape[0]:
+            raise IndexError(
+                f"factor {factor} out of range (n_factors={self.factor_loadings.shape[0]})"
+            )
+        row = self.factor_loadings[factor]
+        order = np.argsort(np.abs(row))[::-1][:n]
+        return [(self.var_names[i], float(row[i])) for i in order]
+
+    def __repr__(self):
+        d = len(self.var_names)
+        edges = int((self.cg_tig == 1).sum())
+        bits = [f"regime={self.regime!r}"]
+        if self.spectral_ratio is not None and self.tau is not None:
+            bits.append(f"R={self.spectral_ratio:.3f} vs tau={self.tau:.3f}")
+        if self.n_factors is not None:
+            bits.append(f"n_factors={self.n_factors}")
+        return f"LucidResult(d={d}, edges={edges}, " + ", ".join(bits) + ")"