ReferenceDiscovery

Discovery

Causal discovery: learn causal structure from data. Import from ergodic.discovery.

Two ingredients drive the searches. A conditional-independence test (CITest) is a class constructed with data, which makes it the stateful oracle the constraint algorithms query: Linear (partial correlation), Discrete (G squared, natively categorical), Mixed (degenerate Gaussian for mixed columns), GCM (any regression Learner becomes a test), and DSeparationOracle as ground truth. A decomposable Score answers local(node, parents) for the score-based searches: BIC (Gaussian), DiscreteBIC, BDeu, and MixedBIC (degenerate Gaussian). Algorithms are thin consumers: PC wraps a test and returns a CPDAG, FCI wraps a test and returns a PAG (latent confounders allowed), MMPC wraps a test and returns an undirected skeleton screen, GES wraps a score and returns a CPDAG, AStar wraps a score and returns the exact optimum, and HillClimb is the fast greedy baseline. discover is the one-line facade, and bootstrap() reruns any method on resamples and returns edge frequencies (BootstrapResult), the first soft output. EdgeProbabilities holds a model's directed-edge probabilities (or bootstrap frequencies) over named variables, and decode() turns one into a DAG with an explicit, versioned decoder (DecodeResult). Background knowledge (including a skeleton restriction through restrict_to) enters through DomainKnowledge.

Time series and panels get their own methods: PCMCI reads a TimeSeriesData or PanelData through a lagged design matrix (lagged_tabular()) and returns a window graph over nodes like "X[t-1]" (TimeSeriesDiscoveryResult), and granger() is the predictive baseline, a pair-level table on purpose rather than a graph.

Also exported from ergodic: DecodeResult, EdgeProbabilities, decode, discover, granger.

BIC#

class BIC(data: TabularData | None = None, *, penalty: float = 1.0) -> None

The Gaussian BIC, decomposable and score-equivalent.

Nodes must be single numeric columns (continuous, datetime, or binary); scoring a wider categorical block raises TypeError lazily, at the first query against that node. Categorical parents are fine: their one-hot columns enter the regression and each column counts toward k.

BIC.penalty#

property BIC.penalty: float

The multiplier on the complexity term.

BIC.local_from_arrays#

def BIC.local_from_arrays(
    y: npt.NDArray[np.float64],
    X: npt.NDArray[np.float64],
) -> float

The same local by direct least squares on (n, 1) y and (n, k) X.

Agrees with the covariance path to numerical precision; the tests hold the two to 1e-8.

BIC.is_bound#

property BIC.is_bound: bool

Whether the score is bound to data and can be evaluated.

BIC.data#

property BIC.data: TabularData | None

The dataset this score is bound to, after dropping incomplete rows.

BIC.n#

property BIC.n: int

Rows the score is computed from (after dropping incomplete rows); 0 before binding.

BIC.variables#

property BIC.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

BIC.bind#

def BIC.bind(data: TabularData) -> Self

A new instance of this score bound to data; never mutates self.

BIC.local#

def BIC.local(node: str, parents: Iterable[str] = ()) -> float

The cached local score of node given parents.

BIC.local_diff#

def BIC.local_diff(node: str, parents: Iterable[str], add: str) -> float

The gain from adding add to the parents of node.

The default is two locals; an implementation may override it with a cheaper incremental form.

BIC.enumerate_locals#

def BIC.enumerate_locals(
    node: str,
    candidates: Iterable[str],
    max_parents: int | None = None,
) -> dict[frozenset[str], float]

Every candidate parent set up to max_parents, scored.

The access pattern of the exact solvers (A*, dynamic programming, ILP): all the local scores up front under an in-degree cap.

BIC.graph_score#

def BIC.graph_score(parents_of: dict[str, frozenset[str]]) -> float

The decomposable total: the sum of locals over a parent map.

BIC.diagnostics#

def BIC.diagnostics() -> dict[str, float]

Run counters: evaluations and cache hits.

BIC.local_from_codes#

def BIC.local_from_codes(
    y: npt.NDArray[np.int64],
    y_card: int,
    zs: tuple[npt.NDArray[np.int64], ...],
    z_cards: tuple[int, ...],
) -> float

The local score from codes: y is (n,) integer codes, zs the parents'.

BOSS#

class BOSS(
    score: Score | str = 'auto',
    *,
    use_bes: bool = True,
    restarts: int = 0,
    seed: int | None = None,
) -> None

Best order score search over a decomposable score.

fit(data, knowledge) returns a DiscoveryResult whose graph carries the CPDAG face when the final PDAG is a valid CPDAG and the mixed face otherwise (forced knowledge orientations can leave the class; diagnostics["cpdag_valid"] says which happened).

BOSS.fit#

def BOSS.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Sweep permutations to convergence and return the class found.

DECODERS#

DECODERS: mappingproxy = mappingproxy({'greedy_v1': 'greedy pair-MAP: best-scoring edges first, skip any that closes a cycle', 'exact_map_v1': 'exact MAP of the factorized pair score, A* over topological orders', 'expected_count_v1': "the MAP order's best K edges, K the rounded expected edge count", 'soft_f1_v1': "the MAP order's best K edges, K from a label-free soft-F1 proxy", 'pair_sampler_v1': "sample each pair's state, add edges in random order, drop cycle-closing ones", 'perturb_and_map_v1': 'Gumbel-perturb the pair states, then solve exact MAP per sample'})

The registered decoders: versioned id to a one-line summary (read-only).

FCI#

class FCI(
    test: CITest | str = 'auto',
    *,
    alpha: float = 0.05,
    max_cond_size: int | None = None,
) -> None

FCI over a conditional-independence test, returning a PAG.

fit(data, knowledge) returns a DiscoveryResult whose graph carries the PAG face when the result is a valid complete PAG and the mixed face otherwise (noisy tests or knowledge-seeded marks). With a perfect oracle and no knowledge the output is exactly the canonical PAG: pag_of_mag(mag_of_dag(truth)).

FCI.fit#

def FCI.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Run FCI and return the discovered graph with its provenance.

GCM#

class GCM(data: TabularData | None = None, *, learner: Any = None) -> None

The generalized covariance measure over a regression learner.

Regresses x and y on z with a fresh clone of the learner per endpoint, forms the elementwise residual product, and reports sqrt(n) * mean(R) / sd(R) against the standard normal, two-sided. Residuals are in-sample, as in the paper's basic form. With an empty conditioning set the endpoints are centered instead, so the statistic is a plain covariance test.

GCM.statistic#

def GCM.statistic(
    x: npt.NDArray[np.float64],
    y: npt.NDArray[np.float64],
    z: npt.NDArray[np.float64],
) -> tuple[float, float, float | None]

Return (statistic, p_value, dof); signature depends on needs.

numeric: statistic(x, y, z) with x (n, dx), y (n, dy), z (n, dz) float blocks. codes: statistic(x, y, z, cards) with x, y (n,) integer codes, z a tuple of (n,) code columns, and cards the matching (card_x, card_y, cards_z).

GCM.is_bound#

property GCM.is_bound: bool

Whether the test is bound to data and can answer queries.

GCM.data#

property GCM.data: TabularData | None

The dataset this test is bound to, after dropping incomplete rows.

GCM.n#

property GCM.n: int

Rows the oracle answers from (after dropping incomplete rows).

GCM.variables#

property GCM.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

GCM.bind#

def GCM.bind(data: TabularData) -> Self

A new instance of this test bound to data; never mutates self.

GCM.test#

def GCM.test(x: str, y: str, given: Iterable[str] = ()) -> CIResult

The cached conditional-independence query, symmetric in x and y.

GCM.independent#

def GCM.independent(
    x: str,
    y: str,
    given: Iterable[str] = (),
    *,
    alpha: float,
) -> bool

The decision the algorithms read: fail to reject at alpha.

GCM.diagnostics#

def GCM.diagnostics() -> dict[str, float]

Run counters: total queries, largest conditioning set, cache hits.

GES#

class GES(
    score: Score | str = 'auto',
    *,
    max_degree: int | None = None,
) -> None

Greedy equivalence search over a decomposable score.

fit(data, knowledge) returns a DiscoveryResult whose graph carries the CPDAG face when the final PDAG is a valid CPDAG and the mixed face otherwise (forced knowledge orientations can leave the class; diagnostics["cpdag_valid"] says which happened).

GES.fit#

def GES.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Run the forward and backward phases and return the class found.

MMPC#

class MMPC(
    test: CITest | str = 'auto',
    *,
    alpha: float = 0.05,
    max_cond_size: int | None = 3,
) -> None

The max-min parents-and-children screen over a conditional-independence test.

MMPC is a cheap skeleton screen to run before an expensive search: it returns the undirected adjacencies the data supports, and DomainKnowledge().restrict_to(result.graph) turns them into a support that PC, GES, or AStar then searches inside, with far fewer tests or score calls. Adjacencies are kept under the AND rule (each endpoint must sit in the other's parents-children set), the stricter of the two symmetry choices, so the screen errs toward sparsity.

Knowledge enters in three modes. Prune: a pair forbidden in both directions is never a candidate. Protect: required and confounded pairs enter every CPC untested, survive the backward sweep, and pass the AND rule unconditionally. Validate: the finished skeleton is checked at the adjacency level only, because a skeleton has no directions, so is_consistent (which demands required edges be directed) cannot apply; orientation constraints (tiers, single-direction forbids, roots, sinks) do not constrain an undirected adjacency and are left to the search that runs inside the support.

fit(data, knowledge) returns a DiscoveryResult whose graph holds one tail-tail edge per surviving pair, with the mixed face.

MMPC.fit#

def MMPC.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Run MMPC and return the undirected skeleton with its provenance.

PC#

class PC(
    test: CITest | str = 'auto',
    *,
    alpha: float = 0.05,
    max_cond_size: int | None = None,
    colliders: str = 'standard',
) -> None

PC-stable over a conditional-independence test.

fit(data, knowledge) returns a DiscoveryResult whose graph carries the CPDAG face when the oriented PDAG is a valid CPDAG and the mixed face otherwise.

PC.fit#

def PC.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Run PC-stable and return the discovered graph with its provenance.

PCMCI#

class PCMCI(
    test: CITest | str = 'auto',
    *,
    alpha: float = 0.05,
    max_lag: int = 1,
    max_cond_size: int | None = None,
) -> None

PCMCI+ over a conditional-independence test, for time-series and panel data.

fit(data) takes TimeSeriesData or PanelData (panel windows stack per entity) and returns a TimeSeriesDiscoveryResult whose graph lives on the window nodes v[t-k].

fit(data, knowledge) reads a DomainKnowledge over the process names, not the window nodes, so a constraint on (x, y) applies at every lag, the contemporaneous slice included:

• A directional forbid x -> y means the link (x, y, lag) is never screened or MCI-tested at any lag and the lag-zero orientation x -> y is banned; the contemporaneous edge itself may survive when y -> x stays permitted (only the direction dies). A pair forbidden in both directions is an adjacency ban: the contemporaneous pair is removed without testing (and records no sepset) and both lagged directions are pruned. An explicit forbid(x, x) removes the self-memory link x[t-1] -> x[t].
• Tiers and order: a later variable never causes an earlier one at any lag, pruned exactly as a forbid; additionally a surviving contemporaneous edge between ordered processes is seeded earlier -> later before the collider phase, and a data-driven arrowhead that disagrees follows the orientation-conflict convention already in this file (the edge stays undirected and is counted in diagnostics["orientation_conflicts"]).
• Roots take no incoming cross-variable link at any lag and sinks no outgoing one; self-memory x[t-1] -> x[t] survives both, since a process can carry its own past without being caused by anything else. A skeleton restriction (restrict_to) bans pairs outside its support the same way, self-memory exempt.
• required (x, y) names no lag, so it cannot be forced; it is validated after the run instead, and fit raises ValueError when no compatible link (a lagged x -> y, an oriented lag-zero x -> y, or a still-undirected lag-zero edge) survived.
• confound has no analogue in a window graph (it carries no bidirected edges) and raises NotImplementedError.

diagnostics["n_knowledge_pruned"] counts what the knowledge took off the table: lagged links never screened, contemporaneous pairs removed outright, and lag-zero direction bans, one count each.

PCMCI.fit#

def PCMCI.fit(
    data: TimeSeriesData | PanelData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> TimeSeriesDiscoveryResult

Run PCMCI+ and return the window graph with its provenance.

knowledge constrains processes, so every constraint applies at every lag; the class docstring lists the mappings. confound raises here (a window graph has no bidirected edges), and a required edge that survived at no lag raises after the run.

RFCI#

class RFCI(
    test: CITest | str = 'auto',
    *,
    alpha: float = 0.05,
    max_cond_size: int | None = None,
) -> None

RFCI over a conditional-independence test: FCI's marks at PC's price.

fit(data, knowledge) returns a DiscoveryResult whose graph carries the PAG face when the output happens to be a valid complete PAG and the mixed face with diagnostics["pag_valid"] = 0.0 otherwise. Unlike FCI, the mixed face is an expected honest outcome here even under a perfect oracle: the RFCI skeleton can keep edges only the possible-d-sep phase would remove, and the result then stands for more than one Markov equivalence class. Domain knowledge enters exactly as in FCI: prune (a pair forbidden in both directions is never tested), protect (required and confounded pairs are never tested or removed), orient through ancestry (the same seeding helpers), validate at the end.

RFCI.fit#

def RFCI.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Run RFCI and return the discovered graph with its provenance.

AStar#

class AStar(
    score: Score | str = 'auto',
    *,
    max_parents: int | None = None,
    node_limit: int = 20,
) -> None

Exact structure search: the highest-BIC DAG, found with A*.

AStar.fit#

def AStar.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Search for a highest-scoring DAG and return it with its provenance.

BDeu#

class BDeu(data: TabularData | None = None, *, ess: float = 10.0) -> None

The Bayesian Dirichlet equivalent uniform score, decomposable.

Every column must be discrete (categorical, ordinal, or binary): the score runs on integer codes, so binding to data with a continuous column raises TypeError. Score-equivalent for any fixed ess, but the value of ess shapes the result; see the module docstring.

BDeu.ess#

property BDeu.ess: float

The equivalent sample size of the Dirichlet prior.

BDeu.local_from_codes#

def BDeu.local_from_codes(
    y: npt.NDArray[np.int64],
    y_card: int,
    zs: tuple[npt.NDArray[np.int64], ...],
    z_cards: tuple[int, ...],
) -> float

The Dirichlet log marginal likelihood from integer codes.

BDeu.is_bound#

property BDeu.is_bound: bool

Whether the score is bound to data and can be evaluated.

BDeu.data#

property BDeu.data: TabularData | None

The dataset this score is bound to, after dropping incomplete rows.

BDeu.n#

property BDeu.n: int

Rows the score is computed from (after dropping incomplete rows); 0 before binding.

BDeu.variables#

property BDeu.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

BDeu.bind#

def BDeu.bind(data: TabularData) -> Self

A new instance of this score bound to data; never mutates self.

BDeu.local#

def BDeu.local(node: str, parents: Iterable[str] = ()) -> float

The cached local score of node given parents.

BDeu.local_diff#

def BDeu.local_diff(node: str, parents: Iterable[str], add: str) -> float

The gain from adding add to the parents of node.

The default is two locals; an implementation may override it with a cheaper incremental form.

BDeu.enumerate_locals#

def BDeu.enumerate_locals(
    node: str,
    candidates: Iterable[str],
    max_parents: int | None = None,
) -> dict[frozenset[str], float]

Every candidate parent set up to max_parents, scored.

The access pattern of the exact solvers (A*, dynamic programming, ILP): all the local scores up front under an in-degree cap.

BDeu.graph_score#

def BDeu.graph_score(parents_of: dict[str, frozenset[str]]) -> float

The decomposable total: the sum of locals over a parent map.

BDeu.diagnostics#

def BDeu.diagnostics() -> dict[str, float]

Run counters: evaluations and cache hits.

BDeu.local_from_arrays#

def BDeu.local_from_arrays(
    y: npt.NDArray[np.float64],
    X: npt.NDArray[np.float64],
) -> float

The local score from arrays: y is (n, dy), X is (n, k).

BootstrapResult#

class BootstrapResult(
    method: str,
    n_runs: int,
    n_failed: int,
    edges: Mapping[frozenset[str], float],
    marks: Mapping[str, float],
    variables: tuple[str, ...],
) -> None

Edge frequencies from rerunning one discovery method on bootstrap resamples.

edges maps an unordered node pair to the fraction of successful runs in which the pair was adjacent, whatever the marks. marks maps a canonical edge glyph (glyph(), for example "X --> Y" or "X <-> Y") to its frequency, so one pair can spread its mass over several glyphs; the glyph frequencies for a pair sum to its adjacency frequency. Both mappings hold only pairs and glyphs that appeared at least once, with keys in sorted order. Frequencies are over the n_runs successful replicates, not over what was attempted. A temporal method keys them on the window nodes, "X[t-1] --> Y[t]".

BootstrapResult.frequency#

def BootstrapResult.frequency(a: str, b: str) -> float

The adjacency frequency of the pair, 0.0 when it never appeared.

BootstrapResult.summary#

def BootstrapResult.summary(threshold: float = 0.5) -> MixedGraph

The stability graph: pairs at or above threshold, one edge each.

A qualifying pair gets its modal glyph, the orientation it wore most often (ties broken by the sorted glyph string, so the result is deterministic). The graph keeps the mixed kind: a vote over graphs carries no equivalence-class guarantee, even when every replicate was a CPDAG. Nodes are all the data's variables, so a variable that never earned a stable edge shows up isolated. threshold must be in (0, 1].

CIResult#

class CIResult(
    statistic: float,
    p_value: float,
    dof: float | None = None,
) -> None

One conditional-independence query: the statistic, its p-value, and dof.

CITest#

class CITest(data: TabularData | None = None) -> None

Base for conditional-independence tests; see the module docstring.

Subclasses set name and needs, may take options in __init__ (set them before calling super().__init__(data), since binding may read them), and implement statistic. The optional _prepare hook runs once after binding, for sufficient statistics.

CITest.is_bound#

property CITest.is_bound: bool

Whether the test is bound to data and can answer queries.

CITest.data#

property CITest.data: TabularData | None

The dataset this test is bound to, after dropping incomplete rows.

CITest.n#

property CITest.n: int

Rows the oracle answers from (after dropping incomplete rows).

CITest.variables#

property CITest.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

CITest.bind#

def CITest.bind(data: TabularData) -> Self

A new instance of this test bound to data; never mutates self.

CITest.test#

def CITest.test(x: str, y: str, given: Iterable[str] = ()) -> CIResult

The cached conditional-independence query, symmetric in x and y.

CITest.independent#

def CITest.independent(
    x: str,
    y: str,
    given: Iterable[str] = (),
    *,
    alpha: float,
) -> bool

The decision the algorithms read: fail to reject at alpha.

CITest.diagnostics#

def CITest.diagnostics() -> dict[str, float]

Run counters: total queries, largest conditioning set, cache hits.

CITest.statistic#

def CITest.statistic(*args: Any) -> tuple[float, float, float | None]

Return (statistic, p_value, dof); signature depends on needs.

numeric: statistic(x, y, z) with x (n, dx), y (n, dy), z (n, dz) float blocks. codes: statistic(x, y, z, cards) with x, y (n,) integer codes, z a tuple of (n,) code columns, and cards the matching (card_x, card_y, cards_z).

DSeparationOracle#

class DSeparationOracle(graph: MixedGraph) -> None

The ground-truth oracle: answers queries by m-separation in a known graph.

Quacks like any bound test (p-value 1.0 when separated, 0.0 when connected), so every constraint algorithm can be checked for exactness and run in teaching mode. Needs no data; bind returns self.

Latent and selection nodes are not query-able variables: the oracle exposes the observed margin, which is what a constraint algorithm searching over data would see. Separation is still computed on the full graph, latents included, and every query conditions on the selection nodes, because data sampled under selection only ever shows the distribution given them.

DSeparationOracle.is_bound#

property DSeparationOracle.is_bound: bool

Always True: the oracle answers from its graph and needs no data.

DSeparationOracle.graph#

property DSeparationOracle.graph: MixedGraph

The graph whose d-separations (m-separations) the oracle reports.

DSeparationOracle.bind#

def DSeparationOracle.bind(data: TabularData) -> Self

A new instance of this test bound to data; never mutates self.

DSeparationOracle.data#

property DSeparationOracle.data: TabularData | None

The dataset this test is bound to, after dropping incomplete rows.

DSeparationOracle.n#

property DSeparationOracle.n: int

Rows the oracle answers from (after dropping incomplete rows).

DSeparationOracle.variables#

property DSeparationOracle.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

DSeparationOracle.test#

def DSeparationOracle.test(
    x: str,
    y: str,
    given: Iterable[str] = (),
) -> CIResult

The cached conditional-independence query, symmetric in x and y.

DSeparationOracle.independent#

def DSeparationOracle.independent(
    x: str,
    y: str,
    given: Iterable[str] = (),
    *,
    alpha: float,
) -> bool

The decision the algorithms read: fail to reject at alpha.

DSeparationOracle.diagnostics#

def DSeparationOracle.diagnostics() -> dict[str, float]

Run counters: total queries, largest conditioning set, cache hits.

DSeparationOracle.statistic#

def DSeparationOracle.statistic(
    *args: Any,
) -> tuple[float, float, float | None]

Return (statistic, p_value, dof); signature depends on needs.

numeric: statistic(x, y, z) with x (n, dx), y (n, dy), z (n, dz) float blocks. codes: statistic(x, y, z, cards) with x, y (n,) integer codes, z a tuple of (n,) code columns, and cards the matching (card_x, card_y, cards_z).

DecodeResult#

class DecodeResult(
    graph: MixedGraph,
    method: str,
    test: str | None = None,
    score: str | None = None,
    alpha: float | None = None,
    n: int = 0,
    diagnostics: Mapping[str, Any] = <factory>,
    sepsets: Mapping[frozenset[str], frozenset[str]] | None = None,
    decoder: str = '',
    candidates: tuple[MixedGraph, ...] = (),
    probabilities: EdgeProbabilities | None = None,
) -> None

A DAG decoded from edge probabilities, with how it was decoded.

A DiscoveryResult whose graph is the decoded DAG (kind dag): the single decoded graph for a deterministic decoder, the consensus (medoid) of the candidates for a sampler. method is "decode:<decoder id>" (a Welt run decoded locally keeps "welt:<model>").

diagnostics always labels the graph: decoded is True, native is False, and it records decoder, parameters, seed, samples, rejected_cycles (nominated edges skipped because they would close a directed cycle), order_skipped (skipped only because they would break the order knowledge) and search_expansions (each one per candidate), threshold_source, edge_budget, objective (the factorized pair log score of graph), n_edges, source_kind, source_provenance, knowledge_applied and note ("not a graph posterior; orientation not certified").

Decoder#

class Decoder(
    method: str = 'greedy',
    threshold: float | None = None,
    samples: int = 1,
    seed: int | None = None,
    max_nodes: int = 20,
) -> None

A decoder and its settings, checked up front and reusable.

Decoder("greedy", threshold=0.3).decode(probabilities) is decode(probabilities, "greedy", threshold=0.3). Pass one to discover() as decoder= to fix every setting before a paid run.

Decoder.stochastic#

property Decoder.stochastic: bool

Whether this decoder samples candidates (and so needs a seed).

Decoder.decode#

def Decoder.decode(
    probabilities: EdgeProbabilities,
    *,
    knowledge: DomainKnowledge | None = None,
) -> DecodeResult

Decode probabilities with these settings; see decode().

DirectLiNGAM#

class DirectLiNGAM(prune: str = 'bic') -> None

DirectLiNGAM over continuous tabular data.

fit(data, knowledge) returns a DiscoveryResult whose graph is a DAG: with independent non-Gaussian noises the whole orientation is identified, so there is no equivalence class to report. The class consumes neither a CI test nor a decomposable score (ingredient is None); the data is its only ingredient.

DirectLiNGAM.fit#

def DirectLiNGAM.fit(
    data: TabularData,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Find the causal order, prune the parents, and return the DAG.

DiscoveryResult#

class DiscoveryResult(
    graph: MixedGraph,
    method: str,
    test: str | None = None,
    score: str | None = None,
    alpha: float | None = None,
    n: int = 0,
    diagnostics: Mapping[str, Any] = <factory>,
    sepsets: Mapping[frozenset[str], frozenset[str]] | None = None,
) -> None

What a discovery run found, and how.

graph is the headline. test or score names the ingredient the algorithm consumed (one of them is None). diagnostics carries the run counters (tests by conditioning size, cache hits, score evaluations); a welt:<model> run carries its Welt provenance there instead (ids, versions and the native edge-score dataframe, see discover()). sepsets maps an unordered removed pair to the set that separated it, when a constraint method produced one.

Discrete#

class Discrete(data: TabularData | None = None) -> None

The G squared conditional-independence test on contingency tables.

Degrees of freedom are adjusted for sparsity, bnlearn-style: each nonempty stratum of the conditioning configuration contributes (card_x - 1) * (card_y - 1), and configurations with no observed rows contribute nothing. The unadjusted product over the full configuration space would overstate the dof badly on sparse tables and make the test blind. When the dof comes out zero (a single-level endpoint, or no rows), the query is vacuous and the result is (0.0, 1.0, 0.0).

Discrete.statistic#

def Discrete.statistic(
    x: npt.NDArray[np.int64],
    y: npt.NDArray[np.int64],
    z: tuple[npt.NDArray[np.int64], ...],
    cards: tuple[int, int, tuple[int, ...]],
) -> tuple[float, float, float | None]

Return (statistic, p_value, dof); signature depends on needs.

numeric: statistic(x, y, z) with x (n, dx), y (n, dy), z (n, dz) float blocks. codes: statistic(x, y, z, cards) with x, y (n,) integer codes, z a tuple of (n,) code columns, and cards the matching (card_x, card_y, cards_z).

Discrete.is_bound#

property Discrete.is_bound: bool

Whether the test is bound to data and can answer queries.

Discrete.data#

property Discrete.data: TabularData | None

The dataset this test is bound to, after dropping incomplete rows.

Discrete.n#

property Discrete.n: int

Rows the oracle answers from (after dropping incomplete rows).

Discrete.variables#

property Discrete.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

Discrete.bind#

def Discrete.bind(data: TabularData) -> Self

A new instance of this test bound to data; never mutates self.

Discrete.test#

def Discrete.test(x: str, y: str, given: Iterable[str] = ()) -> CIResult

The cached conditional-independence query, symmetric in x and y.

Discrete.independent#

def Discrete.independent(
    x: str,
    y: str,
    given: Iterable[str] = (),
    *,
    alpha: float,
) -> bool

The decision the algorithms read: fail to reject at alpha.

Discrete.diagnostics#

def Discrete.diagnostics() -> dict[str, float]

Run counters: total queries, largest conditioning set, cache hits.

DiscreteBIC#

class DiscreteBIC(
    data: TabularData | None = None,
    *,
    penalty: float = 1.0,
) -> None

The multinomial BIC, decomposable and score-equivalent.

Every column must be discrete (categorical, ordinal, or binary): the score runs on integer codes, so binding to data with a continuous column raises TypeError. Parameter counts agree across a Markov equivalence class, so Markov-equivalent DAGs get equal totals, exactly as with the Gaussian BIC.

DiscreteBIC.penalty#

property DiscreteBIC.penalty: float

The multiplier on the complexity term.

DiscreteBIC.local_from_codes#

def DiscreteBIC.local_from_codes(
    y: npt.NDArray[np.int64],
    y_card: int,
    zs: tuple[npt.NDArray[np.int64], ...],
    z_cards: tuple[int, ...],
) -> float

The multinomial local from integer codes; see the module docstring.

DiscreteBIC.is_bound#

property DiscreteBIC.is_bound: bool

Whether the score is bound to data and can be evaluated.

DiscreteBIC.data#

property DiscreteBIC.data: TabularData | None

The dataset this score is bound to, after dropping incomplete rows.

DiscreteBIC.n#

property DiscreteBIC.n: int

Rows the score is computed from (after dropping incomplete rows); 0 before binding.

DiscreteBIC.variables#

property DiscreteBIC.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

DiscreteBIC.bind#

def DiscreteBIC.bind(data: TabularData) -> Self

A new instance of this score bound to data; never mutates self.

DiscreteBIC.local#

def DiscreteBIC.local(node: str, parents: Iterable[str] = ()) -> float

The cached local score of node given parents.

DiscreteBIC.local_diff#

def DiscreteBIC.local_diff(
    node: str,
    parents: Iterable[str],
    add: str,
) -> float

The gain from adding add to the parents of node.

The default is two locals; an implementation may override it with a cheaper incremental form.

DiscreteBIC.enumerate_locals#

def DiscreteBIC.enumerate_locals(
    node: str,
    candidates: Iterable[str],
    max_parents: int | None = None,
) -> dict[frozenset[str], float]

Every candidate parent set up to max_parents, scored.

The access pattern of the exact solvers (A*, dynamic programming, ILP): all the local scores up front under an in-degree cap.

DiscreteBIC.graph_score#

def DiscreteBIC.graph_score(parents_of: dict[str, frozenset[str]]) -> float

The decomposable total: the sum of locals over a parent map.

DiscreteBIC.diagnostics#

def DiscreteBIC.diagnostics() -> dict[str, float]

Run counters: evaluations and cache hits.

DiscreteBIC.local_from_arrays#

def DiscreteBIC.local_from_arrays(
    y: npt.NDArray[np.float64],
    X: npt.NDArray[np.float64],
) -> float

The local score from arrays: y is (n, dy), X is (n, k).

EdgeProbabilities#

class EdgeProbabilities(
    variables: tuple[str, ...],
    values: npt.NDArray[np.float64],
    kind: ProbabilityKind = 'independent',
    calibrated: bool = False,
    provenance: Mapping[str, Any] = <factory>,
) -> None

A matrix of directed-edge probabilities over named variables.

values[i, j] is the number for an edge from variables[i] (the source row) to variables[j] (the target column). The diagonal is zero, every entry is finite and in [0, 1], and the stored matrix is a read-only copy, so the caller's array is never shared or changed.

kind says how to read the numbers:

• "independent": one marginal per directed edge, each on its own. Both directions of a pair may be high at once (AVICI and CDFM return these).
• "pair_state": each unordered pair has three exclusive states (no edge, i -> j, j -> i) and the two directed cells are two of them, so P[i, j] + P[j, i] <= 1; the rest is the no-edge mass.
• "frequency": how often each directed edge appeared over resampled runs (see from_bootstrap()). One edge per pair per run, so the pair simplex holds here too.

calibrated is False unless the source states the numbers match observed edge frequencies; ergodic never sets it on its own. provenance is free-form, read-only metadata (model, version, job id, threshold, ...).

None of this is a posterior over graphs. The marginals do not account for acyclicity or for dependence between edges; a decoded graph built from them is a modelling choice, not the model's claim.

EdgeProbabilities.to_pair_states#

def EdgeProbabilities.to_pair_states() -> EdgeProbabilities

The pair-state reading: three exclusive states per unordered pair.

For "independent" marginals p = P[i, j] and q = P[j, i], the two edges are treated as independent coin flips and conditioned on "not both" (a DAG cannot hold a two-cycle): no edge (1 - p)(1 - q), i -> j p(1 - q), j -> i q(1 - p), renormalized. When both are exactly 1 the limit is undecided, so the pair becomes (0, 0.5, 0.5). This keeps each model's direction ranking and never looks at a true graph. It is the conversion Ergodic's research decoders were evaluated with.

"pair_state" values are returned unchanged (the same object). "frequency" values already satisfy the pair simplex (one edge per pair per run), so they are read as pair states as they are: a new object of kind "pair_state" with the same numbers.

Returns: Edge probabilities of kind "pair_state", with the conversion recorded in provenance["pair_state_conversion"].

EdgeProbabilities.existence#

def EdgeProbabilities.existence() -> npt.NDArray[np.float64]

Skeleton probabilities: the chance each pair is adjacent, either way.

Read from the pair-state view (to_pair_states()), so it is P[i, j] + P[j, i] there: for "independent" marginals that is (p(1 - q) + q(1 - p)) / (1 - pq). The matrix is symmetric with a zero diagonal and read-only.

EdgeProbabilities.threshold#

def EdgeProbabilities.threshold(t: float = 0.5) -> MixedGraph

The directed edges whose value is above t, as a diagnostic graph.

Keeps every edge with values[i, j] > t (strictly, as AVICI and CDFM threshold natively). The result is a mixed-kind graph over every variable and may hold directed cycles. A pair whose two directions both pass cannot be drawn as two edges (one edge per pair), so it shows as an undirected X --- Y edge: adjacent, direction unresolved. This is not a DAG and not a decoded graph; use decode() for one.

Returns: A mixed MixedGraph.

EdgeProbabilities.to_frame#

def EdgeProbabilities.to_frame() -> pd.DataFrame

The matrix as a dataframe: index source, columns target.

EdgeProbabilities.from_frame#

classmethod EdgeProbabilities.from_frame(
    frame: pd.DataFrame,
    *,
    kind: ProbabilityKind = 'independent',
    calibrated: bool = False,
    provenance: Mapping[str, Any] | None = None,
) -> EdgeProbabilities

Build from a square dataframe: rows the source, columns the target.

The index gives the variable order; the columns must hold the same names (in any order) and are reordered to match. Labels are read as strings.

EdgeProbabilities.from_bootstrap#

classmethod EdgeProbabilities.from_bootstrap(
    result: BootstrapResult,
    *,
    undirected: Literal['split', 'drop'] = 'split',
) -> EdgeProbabilities

Directed-edge frequencies from a bootstrap() run.

Each canonical glyph frequency lands on its directed cell: "X --> Y" on [X, Y] and "X <-- Y" on [Y, X]. An undirected X --- Y (a CPDAG edge, "one of the two directions") is split half and half under undirected="split" and dropped under "drop". Bidirected and circle-marked glyphs say something other than "this direction" and are always dropped; the dropped mass is recorded in provenance["dropped_frequency"]. Frequencies measure resampling stability, not probability.

Returns: Edge probabilities of kind "frequency".

GRaSP#

class GRaSP(score: Score | str = 'auto', *, depth: int = 3) -> None

Greedy relaxations of the sparsest permutation over a decomposable score.

fit(data, knowledge) returns a DiscoveryResult whose graph carries the CPDAG face when the final PDAG is a valid CPDAG and the mixed face otherwise (forced knowledge orientations can leave the class; diagnostics["cpdag_valid"] says which happened).

GRaSP.fit#

def GRaSP.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Walk tucks to convergence and return the class found.

Glasso#

class Glasso(
    lam: float | None = None,
    *,
    gamma: float = 0.5,
    n_lambdas: int = 100,
    tol: float = 0.001,
    max_iter: int = 100,
) -> None

The graphical-lasso screen: an undirected skeleton from a sparse precision.

Like MMPC, this is a cheap screen to run before an expensive search: DomainKnowledge().restrict_to(result.graph) turns the recovered support into a constraint that PC, GES, or AStar then searches inside, with far fewer tests or score calls. The honesty note: a precision zero is a conditional independence under Gaussianity, so the graph is a skeleton source, not causal structure; it says which adjacencies the data supports, never which way they point.

The solve is FHT 2008 block coordinate descent on the correlation matrix (the graph, not the scale, is the deliverable), with the diagonal unpenalized (the scikit-learn convention; the paper's differs, see the module docstring) and the shrunk initialization that keeps every column subproblem positive definite.

Knowledge enters through the penalty, and only its adjacency content binds, because the screen is undirected. A pair forbidden in both directions (an adjacency ban, including pairs outside a restrict_to support) is a hard zero, excluded from the updates entirely. A required pair (either direction) is unpenalized, so it survives every lambda the data gives it any conditional association at. Directional knowledge (tiers, one-way forbids, roots, sinks) does not constrain an undirected adjacency and is left to the search that runs inside the support.

fit(data, knowledge) needs a TabularData whose variables are all continuous (a Gaussian model) and returns a DiscoveryResult whose graph holds one tail-tail edge per nonzero precision pair, on the mixed face. The chosen penalty lives in the diagnostics, with the EBIC value when the path ran.

Glasso.fit#

def Glasso.fit(
    data: TabularData,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Run the screen and return the undirected skeleton with its provenance.

GrangerResult#

class GrangerResult(
    pvalues: Mapping[tuple[str, str], float],
    significant: frozenset[tuple[str, str]],
    max_lag: int,
    alpha: float,
    n: int,
) -> None

Granger tests for every ordered pair: a predictive table, not a graph.

pvalues maps (cause, effect) to the likelihood-ratio p-value; significant holds the pairs with p <= alpha. n counts the stacked window rows the regressions used.

HillClimb#

class HillClimb(
    score: Score | str = 'auto',
    *,
    tabu: int = 10,
    max_iter: int = 10000,
    max_parents: int | None = None,
) -> None

Greedy DAG search: add, remove, and reverse moves with a tabu ring.

Steepest ascent from the empty (or knowledge-seeded) graph, returning the best state seen. n_iterations in the diagnostics counts applied moves, so it always equals n_adds + n_removes + n_reverses.

HillClimb.fit#

def HillClimb.fit(
    data: TabularData | None = None,
    knowledge: DomainKnowledge | None = None,
) -> DiscoveryResult

Climb to a high-scoring DAG and return it with its provenance.

Linear#

class Linear(data: TabularData | None = None) -> None

Partial correlation through the Fisher z transform.

Residualizes x and y on [1, z], computes the correlation r of the residuals, and reports sqrt(n - dz - 3) * atanh(r) against the standard normal, two-sided. The sign of r survives in the statistic; the p-value does not depend on it.

Linear.statistic#

def Linear.statistic(
    x: npt.NDArray[np.float64],
    y: npt.NDArray[np.float64],
    z: npt.NDArray[np.float64],
) -> tuple[float, float, float | None]

Return (statistic, p_value, dof); signature depends on needs.

numeric: statistic(x, y, z) with x (n, dx), y (n, dy), z (n, dz) float blocks. codes: statistic(x, y, z, cards) with x, y (n,) integer codes, z a tuple of (n,) code columns, and cards the matching (card_x, card_y, cards_z).

Linear.is_bound#

property Linear.is_bound: bool

Whether the test is bound to data and can answer queries.

Linear.data#

property Linear.data: TabularData | None

The dataset this test is bound to, after dropping incomplete rows.

Linear.n#

property Linear.n: int

Rows the oracle answers from (after dropping incomplete rows).

Linear.variables#

property Linear.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

Linear.bind#

def Linear.bind(data: TabularData) -> Self

A new instance of this test bound to data; never mutates self.

Linear.test#

def Linear.test(x: str, y: str, given: Iterable[str] = ()) -> CIResult

The cached conditional-independence query, symmetric in x and y.

Linear.independent#

def Linear.independent(
    x: str,
    y: str,
    given: Iterable[str] = (),
    *,
    alpha: float,
) -> bool

The decision the algorithms read: fail to reject at alpha.

Linear.diagnostics#

def Linear.diagnostics() -> dict[str, float]

Run counters: total queries, largest conditioning set, cache hits.

Mixed#

class Mixed(data: TabularData | None = None) -> None

Partial canonical correlations with Bartlett's chi-square.

Residualizes the x and y blocks on [1, z], takes the reduced QR of each residual block, and reads the canonical correlations rho_i off the singular values of Qx.T @ Qy (clipped below 1). Bartlett's statistic -(n - dz - (dx + dy + 3) / 2) * sum(log(1 - rho_i^2)) is referred to a chi-square with dx * dy degrees of freedom.

With scalar endpoints (dx == dy == 1) the single canonical correlation is the absolute partial correlation, so Mixed agrees with Linear asymptotically. A zero-width block (a single-level categorical) carries no information; that query is vacuous and returns (0.0, 1.0, 0.0).

Mixed.statistic#

def Mixed.statistic(
    x: npt.NDArray[np.float64],
    y: npt.NDArray[np.float64],
    z: npt.NDArray[np.float64],
) -> tuple[float, float, float | None]

Return (statistic, p_value, dof); signature depends on needs.

numeric: statistic(x, y, z) with x (n, dx), y (n, dy), z (n, dz) float blocks. codes: statistic(x, y, z, cards) with x, y (n,) integer codes, z a tuple of (n,) code columns, and cards the matching (card_x, card_y, cards_z).

Mixed.is_bound#

property Mixed.is_bound: bool

Whether the test is bound to data and can answer queries.

Mixed.data#

property Mixed.data: TabularData | None

The dataset this test is bound to, after dropping incomplete rows.

Mixed.n#

property Mixed.n: int

Rows the oracle answers from (after dropping incomplete rows).

Mixed.variables#

property Mixed.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

Mixed.bind#

def Mixed.bind(data: TabularData) -> Self

A new instance of this test bound to data; never mutates self.

Mixed.test#

def Mixed.test(x: str, y: str, given: Iterable[str] = ()) -> CIResult

The cached conditional-independence query, symmetric in x and y.

Mixed.independent#

def Mixed.independent(
    x: str,
    y: str,
    given: Iterable[str] = (),
    *,
    alpha: float,
) -> bool

The decision the algorithms read: fail to reject at alpha.

Mixed.diagnostics#

def Mixed.diagnostics() -> dict[str, float]

Run counters: total queries, largest conditioning set, cache hits.

MixedBIC#

class MixedBIC(
    data: TabularData | None = None,
    *,
    penalty: float = 1.0,
) -> None

The degenerate-Gaussian BIC, decomposable and score-equivalent.

Any column kind is a valid node: a continuous node is its single column and a categorical node is its one-hot block (first level dropped), scored jointly. local = -(n/2) log det Sigma_{Y|X} - penalty * p/2 * log n with p = dy dx + dy (dy + 1) / 2 for a dy-column node and dx total parent columns.

MixedBIC.local_from_arrays#

def MixedBIC.local_from_arrays(
    y: npt.NDArray[np.float64],
    X: npt.NDArray[np.float64],
) -> float

The same local by direct least squares on (n, dy) y and (n, k) X.

Agrees with the covariance path to numerical precision; the tests hold the two to 1e-8.

MixedBIC.penalty#

property MixedBIC.penalty: float

The multiplier on the complexity term.

MixedBIC.is_bound#

property MixedBIC.is_bound: bool

Whether the score is bound to data and can be evaluated.

MixedBIC.data#

property MixedBIC.data: TabularData | None

The dataset this score is bound to, after dropping incomplete rows.

MixedBIC.n#

property MixedBIC.n: int

Rows the score is computed from (after dropping incomplete rows); 0 before binding.

MixedBIC.variables#

property MixedBIC.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

MixedBIC.bind#

def MixedBIC.bind(data: TabularData) -> Self

A new instance of this score bound to data; never mutates self.

MixedBIC.local#

def MixedBIC.local(node: str, parents: Iterable[str] = ()) -> float

The cached local score of node given parents.

MixedBIC.local_diff#

def MixedBIC.local_diff(
    node: str,
    parents: Iterable[str],
    add: str,
) -> float

The gain from adding add to the parents of node.

The default is two locals; an implementation may override it with a cheaper incremental form.

MixedBIC.enumerate_locals#

def MixedBIC.enumerate_locals(
    node: str,
    candidates: Iterable[str],
    max_parents: int | None = None,
) -> dict[frozenset[str], float]

Every candidate parent set up to max_parents, scored.

The access pattern of the exact solvers (A*, dynamic programming, ILP): all the local scores up front under an in-degree cap.

MixedBIC.graph_score#

def MixedBIC.graph_score(parents_of: dict[str, frozenset[str]]) -> float

The decomposable total: the sum of locals over a parent map.

MixedBIC.diagnostics#

def MixedBIC.diagnostics() -> dict[str, float]

Run counters: evaluations and cache hits.

MixedBIC.local_from_codes#

def MixedBIC.local_from_codes(
    y: npt.NDArray[np.int64],
    y_card: int,
    zs: tuple[npt.NDArray[np.int64], ...],
    z_cards: tuple[int, ...],
) -> float

The local score from codes: y is (n,) integer codes, zs the parents'.

Score#

class Score(data: TabularData | None = None) -> None

Base for decomposable scores; see the module docstring.

Subclasses set name, needs, and the property flags, may take options in __init__ (set them before calling super().__init__(data)), and implement the method matching needs: local_from_arrays for "numeric", local_from_codes for "codes". _prepare runs once after binding for sufficient statistics; a native score may also override _local_value to skip array assembly entirely (the covariance path in BIC).

Score.is_bound#

property Score.is_bound: bool

Whether the score is bound to data and can be evaluated.

Score.data#

property Score.data: TabularData | None

The dataset this score is bound to, after dropping incomplete rows.

Score.n#

property Score.n: int

Rows the score is computed from (after dropping incomplete rows); 0 before binding.

Score.variables#

property Score.variables: tuple[str, ...]

The variables this oracle answers about, sorted.

Score.bind#

def Score.bind(data: TabularData) -> Self

A new instance of this score bound to data; never mutates self.

Score.local#

def Score.local(node: str, parents: Iterable[str] = ()) -> float

The cached local score of node given parents.

Score.local_diff#

def Score.local_diff(node: str, parents: Iterable[str], add: str) -> float

The gain from adding add to the parents of node.

The default is two locals; an implementation may override it with a cheaper incremental form.

Score.enumerate_locals#

def Score.enumerate_locals(
    node: str,
    candidates: Iterable[str],
    max_parents: int | None = None,
) -> dict[frozenset[str], float]

Every candidate parent set up to max_parents, scored.

The access pattern of the exact solvers (A*, dynamic programming, ILP): all the local scores up front under an in-degree cap.

Score.graph_score#

def Score.graph_score(parents_of: dict[str, frozenset[str]]) -> float

The decomposable total: the sum of locals over a parent map.

Score.diagnostics#

def Score.diagnostics() -> dict[str, float]

Run counters: evaluations and cache hits.

Score.local_from_arrays#

def Score.local_from_arrays(
    y: npt.NDArray[np.float64],
    X: npt.NDArray[np.float64],
) -> float

The local score from arrays: y is (n, dy), X is (n, k).

Score.local_from_codes#

def Score.local_from_codes(
    y: npt.NDArray[np.int64],
    y_card: int,
    zs: tuple[npt.NDArray[np.int64], ...],
    z_cards: tuple[int, ...],
) -> float

The local score from codes: y is (n,) integer codes, zs the parents'.

TimeSeriesDiscoveryResult#

class TimeSeriesDiscoveryResult(
    graph: MixedGraph,
    method: str,
    test: str | None = None,
    score: str | None = None,
    alpha: float | None = None,
    n: int = 0,
    diagnostics: Mapping[str, Any] = <factory>,
    sepsets: Mapping[frozenset[str], frozenset[str]] | None = None,
    max_lag: int = 0,
    links: Mapping[tuple[str, str, int], float] = <factory>,
) -> None

A DiscoveryResult over window nodes, with the per-link p-values.

links maps (source, target, lag) to the p-value of the test that kept the link. A lag of one or more is the lagged edge source[t-lag] -> target[t]. A lag of zero is a contemporaneous edge: an oriented edge is keyed (source, target), an undirected edge is keyed once with the pair in sorted order.

TimeSeriesDiscoveryResult.summary_edges#

def TimeSeriesDiscoveryResult.summary_edges() -> frozenset[tuple[str, str]]

Variable-level pairs implied by the links.

(x, y) is present when any lagged link x -> y exists or a contemporaneous edge x -> y is oriented in the window graph. An undirected contemporaneous pair contributes both orders, and a lagged self-link (autocorrelation) contributes (x, x).

bootstrap#

def bootstrap(
    data: TabularData | TimeSeriesData | PanelData,
    *,
    method: str = 'pc',
    n_bootstrap: int = 100,
    seed: int = 0,
    knowledge: DomainKnowledge | None = None,
    block_length: int | None = None,
    **kwargs: Any,
) -> BootstrapResult

Rerun a discovery method on bootstrap resamples and count what survives.

Answers the question a single run cannot: how stable is this edge under resampling? Each replicate draws a resample, runs discover(sample, method=method, knowledge=knowledge, **kwargs), and the result aggregates adjacency and glyph frequencies over the runs that finished. Extra keyword arguments flow to the algorithm unchanged (alpha, max_cond_size, score=..., max_lag=...), so any registered method, iid or temporal, sits behind the same call.

An iid method on a TabularData draws data.n_rows rows with replacement, the classic bootstrap; block_length does not apply there and raises. A temporal method ("pcmci") never resamples rows, which would shred the serial dependence it reads. On a TimeSeriesData it runs a moving-block bootstrap: ceil(T / L) blocks of L consecutive rows (block_length, default max(ceil(sqrt(T)), 2 * (max_lag + 1))) are drawn with replacement and enter the replicate as separate entities of a panel, "block_0", "block_1", ... with fresh time stamps. Blocks-as-entities is the point: gluing blocks back into one series would build lagged windows that straddle two unrelated fragments and let fake adjacencies in, while per-entity pooling never reads across a block seam. On a PanelData it runs the cluster bootstrap instead, drawing n_entities entities whole with replacement (relabeled in draw order so repeated draws stay distinct); block_length does not apply to a panel and raises. One honesty note: block resampling treats the series as stationary at the block scale, so a trend or a regime change inside the sampling window leaks into the frequencies. An iid method on temporal data still refuses, inside the replicate, exactly as discover() does.

The frequencies measure resampling variability, not Bayesian uncertainty; they are not posterior probabilities over graphs (ergodic has no graph posteriors). The run is deterministic: the same data, seed, and settings give the same result. A temporal method keys the frequencies on window-graph glyphs like "X[t-1] --> Y[t]".

A replicate that raises ValueError or RuntimeError (degenerate resamples happen: a constant column, a knowledge conflict on noisy data) counts in n_failed and is skipped; frequencies are over the successful runs only. When every replicate fails, the wrapper raises RuntimeError chaining the last error.

decode#

def decode(
    probabilities: EdgeProbabilities,
    method: str = 'greedy',
    *,
    threshold: float | None = None,
    knowledge: DomainKnowledge | None = None,
    samples: int = 1,
    seed: int | None = None,
    max_nodes: int = 20,
) -> DecodeResult

Decode edge probabilities into a DAG with a named, versioned decoder.

The input is read as pair states (to_pair_states()) and every decoder works on the edge weight w(i -> j) = log P(i -> j) - log P(no edge). The decoders:

• "greedy" (greedy_v1, the default): for each pair nominate its better direction when w > 0, or, with threshold=t, when the input value of that direction is above t; add nominations by decreasing weight and skip any that would close a cycle. On independent marginals, w > 0 is p > 0.5, AVICI's native rule; threshold lets a model's own cut (CDFM's adaptive one) choose the edges instead.
• "exact_map" (exact_map_v1): the DAG maximizing the factorized pair score, solved exactly by A* over topological orders.
• "expected_count" / "soft_f1" (expected_count_v1, soft_f1_v1): the best K edges compatible with the exact MAP order, K the rounded expected edge count or the size that maximizes a label-free soft-F1 proxy. They may keep edges a MAP decoder drops.
• "pair_sampler" / "perturb_and_map" (pair_sampler_v1, perturb_and_map_v1): samples candidate DAGs, from independent pair-state draws with random-order cycle rejection, or from exact MAP under Gumbel-perturbed pair states. graph is their medoid (the candidate with the smallest total pair-state distance to the others, earliest on an exact tie; the research posterior_medoid can pick another tied candidate). Candidates are not posterior samples.

Ties are broken by variable order, so a decoder is deterministic given its inputs (and seed). The input is never changed and the global numpy random state is never touched.

Knowledge constrains the decode. Forbidden edges, tiers, roots, sinks and a skeleton restriction remove directed edges (their weight becomes minus infinity). Order (order, tiers) is enforced as one time order: the decoded graph has a topological order consistent with every order statement and every edge at once. That is exactly the tier reading, and stricter than violations() for separate order pairs (it only asks that no later variable be an ancestor of an earlier one), so a DAG that satisfies each pair alone but no single order is ruled out. Required edges go in first, whatever their weight, and must not form a cycle. confound is refused: a DAG cannot hold a latent confounder. The decoded graph is checked against the knowledge afterwards.

Returns: A DecodeResult: graph a DAG over every variable, with the decoder, its parameters and the diagnostics recorded.

discover#

def discover(
    data: TabularData | TimeSeriesData | PanelData,
    method: str = 'pc',
    *,
    test: Any = None,
    score: Any = None,
    knowledge: DomainKnowledge | None = None,
    alpha: float = 0.05,
    **kwargs: Any,
) -> DiscoveryResult

Learn causal structure from data with the named method.

test feeds a constraint method (a CITest instance built without data, since discover binds it to data, or an alias like "linear"; a bound test runs through the algorithm class, as in PC(test=GCM(data)).fit()); score feeds a score-based method (a Score instance or an alias like "bic", "bdeu"); both default to "auto", picked from the column schema. Each algorithm declares which ingredient it consumes. Extra keyword arguments pass to the algorithm (max_cond_size, max_parents, node_limit, max_degree, tabu, max_lag).

The iid methods read a TabularData and reject temporal shapes (rows are not exchangeable); method="pcmci" reads a TimeSeriesData or PanelData and returns a window graph over lagged nodes like "X[t-1]".

A few methods consume neither ingredient because the machinery is the algorithm: "lingam" (regression and residual independence) and "glasso" (a sparse precision estimate). They reject test= and score=, and alpha does not apply (their own knobs, like lam or prune, pass through the keyword arguments).

method="welt:<model>" (for example "welt:arrow") runs a hosted Welt model through discover(), one paid Welt job. It needs pip install "ergodic[welt]", takes no test, score or alpha, and passes only client, seed, timeout and decoder through the keyword arguments; anything else raises before any request. knowledge is refused unless decoder= asks for local decoding of the model's edge scores (decode()), which is also what lets a non-DAG model such as "welt:avici" return a (decoded) DAG.

granger#

def granger(
    data: TimeSeriesData,
    *,
    max_lag: int = 1,
    alpha: float = 0.05,
) -> GrangerResult

Granger-test every ordered pair of variables in data.

For each effect Y, the full regression puts Y[t] on an intercept plus every variable's lags 1..max_lag; the restricted one drops the candidate cause's lags. All variables must be continuous. See the module docstring for what this does and does not mean.

lagged_node#

def lagged_node(name: str, lag: int) -> str

The window-column name for name at lag steps back.

Lag zero reads "X[t]", lag k reads "X[t-k]". The name must carry no whitespace, since the graph layer parses edges on whitespace.

lagged_tabular#

def lagged_tabular(
    data: TimeSeriesData | PanelData,
    max_lag: int,
) -> TabularData

Stack the windows of data into a tabular frame of lagged columns.

Row t holds every variable at lags 0..max_lag; the column for variable v at lag k is named lagged_node(v, k) and holds v at time t - k. Column types carry over: continuous stays continuous, categorical keeps its categories, ordinal keeps its level order.

A time series needs at least max_lag + 8 rows. Panel windows are built within each entity (sorted entity order) and stacked; entities with fewer than max_lag + 2 rows are skipped with a warning, and it is an error when no rows survive. Unequal time steps warn too: lags count rows, not calendar time.

parse_lagged#

def parse_lagged(node: str) -> tuple[str, int]

Invert lagged_node(): "X[t-2]" becomes ("X", 2).