API reference

Explainer and results

Counterfactual explainer for a tree-ensemble model.

model may be a native model object, a dump file path/dict, or an EnsembleIR. background fits the distance normalizers; alternatively pass normalizers explicitly (array or name->sigma dict).

Source code in src/treecf/api.py

class Explainer:
    """Counterfactual explainer for a tree-ensemble model.

    ``model`` may be a native model object, a dump file path/dict, or an
    ``EnsembleIR``. ``background`` fits the distance normalizers;
    alternatively pass ``normalizers`` explicitly (array or name->sigma dict).
    """

    _rust_cache: dict[str, object]  # marshaled Rust objects, created on first solve
    _prepared_trees: tuple[TreeArrays, ...]  # vectorized-verify arrays, created on first batch

    def __init__(
        self,
        model: object,
        background: FloatArray | None = None,
        constraints: Sequence[Constraint] = (),
        weights: dict[str, float] | None = None,
        normalizers: FloatArray | dict[str, float] | None = None,
        value_policy: dict[str, ValuePolicy] | None = None,
        plausibility: Plausibility | None = None,
    ) -> None:
        self.ir = model if isinstance(model, EnsembleIR) else parse_model(model)
        names = self.ir.feature_names
        self.compiled = compile_constraints(constraints, names)
        self.plausibility = plausibility
        if plausibility is not None:
            if plausibility.if_ir.n_features != self.ir.n_features:
                raise TreecfError("plausibility forest must share the model's feature space")
            if self.compiled.allow_missing:
                raise TreecfError(
                    "plausibility with AllowMissing is not supported "
                    "(isolation forests define no NaN routing)"
                )
        self.background = (
            None if background is None else np.asarray(background, dtype=np.float64)
        )
        self.sigma = _resolve_sigma(names, background, normalizers)
        self.weights = np.array([(weights or {}).get(name, 1.0) for name in names])
        self.value_policy = value_policy or {}
        for name, policy in self.value_policy.items():
            if name not in names:
                raise TreecfError(f"value_policy references unknown feature {name!r}")
            if isinstance(policy, str) and policy not in ("raw", "integer"):
                raise TreecfError(f"unknown value policy {policy!r} for {name!r}")

    def explain(
        self,
        x: FloatArray,
        target: Target,
        backend: str = "genetic",
        time_budget_s: float = 10.0,
        sparsity_weight: float = 0.0,
        seed: int | None = None,
    ) -> Counterfactual | Infeasible | dict[str, object]:
        """Search for a counterfactual (or one per band for ``Target.bands``).

        ``backend="genetic"`` runs the bundled Rust engine (default);
        ``backend="python"`` runs the reference numpy implementation of the
        same algorithm. Every result is float-verified before being returned.
        """
        x = np.asarray(x, dtype=np.float64)
        if self.plausibility is not None and np.isnan(x).any():
            raise TreecfError("plausibility with missing factual values is not supported")
        if backend in ("genetic", "genetic-rust"):
            rust = True
        elif backend == "python":
            rust = False
        else:
            raise TreecfError(f"unknown backend {backend!r}; use 'genetic' or 'python'")

        if target.bands_spec is not None:
            results: dict[str, object] = {}
            for name, interval in target.band_intervals(self.ir.link).items():
                results[name] = self._explain_genetic(
                    x, interval, time_budget_s, sparsity_weight, seed, rust=rust
                )
            return results
        interval = target.raw_interval(self.ir.link)
        return self._explain_genetic(
            x, interval, time_budget_s, sparsity_weight, seed, rust=rust
        )

    def explain_batch(
        self,
        X: FloatArray,
        target: Target,
        n_per_example: int = 1,
        diversity: str = "seeds",
        ids: Sequence[object] | None = None,
        backend: str = "genetic",
        time_budget_s: float = 10.0,
        sparsity_weight: float = 0.0,
        seed: int = 0,
        coalitions: Mapping[str, Sequence[str]] | None = None,
        include_full: bool = False,
    ) -> Any:
        """Mass-produce counterfactuals for a dataset; see ``treecf.batch``.

        ``n_per_example`` alternatives per row via ``diversity="seeds"`` (distinct
        change-sets from different seeds, best-effort) or ``"lever-blocking"``
        (freeze each plan's biggest lever; also records essential levers).
        ``diversity="coalitions"`` instead produces one plan per named feature
        group in ``coalitions`` per row (``n_per_example`` unused; see
        ``explain_coalitions``). The returned ``BatchResult`` supports
        save/load/for_id/to_frame.

        Solves run in parallel inside the Rust engine; ``time_budget_s`` is
        per solve, so a solve that hits its wall-clock budget while sharing
        cores may stop earlier than it would sequentially (results are
        otherwise identical to solving row by row).
        """
        from treecf.batch import explain_batch

        return explain_batch(
            self, X, target, n_per_example=n_per_example, diversity=diversity,
            ids=ids, backend=backend, time_budget_s=time_budget_s,
            sparsity_weight=sparsity_weight, seed=seed,
            coalitions=coalitions, include_full=include_full,
        )

    def explain_coalitions(
        self,
        x: FloatArray,
        target: Target,
        coalitions: Mapping[str, Sequence[str]],
        include_full: bool = False,
        backend: str = "genetic",
        time_budget_s: float = 10.0,
        sparsity_weight: float = 0.0,
        seed: int | None = None,
    ) -> dict[str, Counterfactual | Infeasible]:
        """One counterfactual per named feature coalition (opt-in mode).

        Each coalition is solved with every feature *outside* it frozen, so a
        plan only ever asks for changes within one group — grouped recourse
        instead of one plan that mixes unrelated levers. Coalitions may
        overlap; features in no coalition are never modified; an
        ``Infeasible`` for a coalition means that group alone cannot reach
        the target. ``include_full=True`` prepends an unrestricted baseline
        under the reserved key ``"(all levers)"``. One solve per coalition
        (milliseconds each); this mode is optional and never the default.
        """
        if target.bands_spec is not None:
            raise TreecfError(
                "Target.bands is not supported in explain_coalitions; loop bands explicitly"
            )
        normalized = _validate_coalitions(coalitions, self.ir.feature_names, include_full)
        results: dict[str, Counterfactual | Infeasible] = {}
        if include_full:
            results[_ALL_LEVERS] = self._explain_one(
                x, target, backend, time_budget_s, sparsity_weight, seed
            )
        for name, clone in self._coalition_explainers(normalized).items():
            results[name] = clone._explain_one(
                x, target, backend, time_budget_s, sparsity_weight, seed
            )
        return results

    def _explain_one(
        self,
        x: FloatArray,
        target: Target,
        backend: str,
        time_budget_s: float,
        sparsity_weight: float,
        seed: int | None,
    ) -> Counterfactual | Infeasible:
        """`explain` for a single-interval target, with the bands arm ruled out."""
        result = self.explain(
            x, target, backend=backend, time_budget_s=time_budget_s,
            sparsity_weight=sparsity_weight, seed=seed,
        )
        assert not isinstance(result, dict)  # bands are rejected by the callers
        return result

    def _coalition_explainers(
        self, coalitions: dict[str, tuple[str, ...]]
    ) -> dict[str, Explainer]:
        """One freeze-complement clone per coalition (Rust ensemble shared)."""
        names = self.ir.feature_names
        return {
            name: self._with_extra_freezes([f for f in names if f not in set(members)])
            for name, members in coalitions.items()
        }

    def _with_extra_freezes(self, features: Sequence[str]) -> Explainer:
        """Clone with additional Freeze constraints (lever-blocking, coalitions).

        ``AllowMissing`` on a newly frozen feature is dropped: a frozen value
        cannot transition to NaN, and keeping both would (correctly) fail
        constraint validation.
        """
        from treecf.constraints.objects import AllowMissing, Freeze

        frozen = set(features)
        kept = [
            c
            for c in self.compiled.constraints
            if not (isinstance(c, AllowMissing) and c.feature in frozen)
        ]
        clone = Explainer(
            self.ir,
            background=self.background,
            constraints=kept + [Freeze(f) for f in features],
            weights=dict(zip(self.ir.feature_names, self.weights.tolist(), strict=True)),
            normalizers=self.sigma,
            value_policy=self.value_policy,
            plausibility=self.plausibility,
        )
        # Same frozen IR -> the marshaled Rust ensembles are reusable; only the
        # constraints differ, so that cache entry is deliberately left out.
        parent_cache: dict[str, object] = getattr(self, "_rust_cache", {})
        clone._rust_cache = {
            key: parent_cache[key]
            for key in ("ensemble", "if_ensemble")
            if key in parent_cache
        }
        return clone

    def _explain_genetic(
        self,
        x: FloatArray,
        interval: tuple[float, float],
        time_budget_s: float,
        sparsity_weight: float,
        seed: int | None,
        rust: bool = True,
    ) -> Counterfactual | Infeasible:
        if rust:
            from treecf.backends.genetic_rust import solve_genetic_rust

            if not hasattr(self, "_rust_cache"):
                self._rust_cache = {}
            result = solve_genetic_rust(
                self.ir,
                x,
                interval,
                self.compiled,
                self.sigma,
                self.weights,
                lam=sparsity_weight,
                background=self.background,
                plausibility=self._plausibility_bound(),
                seed=seed,
                time_budget_s=time_budget_s,
                cache=self._rust_cache,
            )
        else:
            from treecf.backends.genetic import solve_genetic

            result = solve_genetic(
                self.ir,
                x,
                interval,
                self.compiled,
                self.sigma,
                self.weights,
                lam=sparsity_weight,
                background=self.background,
                plausibility=self._plausibility_bound(),
                seed=seed,
                time_budget_s=time_budget_s,
            )
        if result.x_cf is None:
            return Infeasible(reason="heuristic search exhausted (genetic backend)")
        return self._finalize_candidate(x, result.x_cf, interval, result.stats)

    def _finalize_candidate(
        self,
        x: FloatArray,
        x_cf: FloatArray,
        interval: tuple[float, float],
        stats: dict[str, object],
        score: float | None = None,
    ) -> Counterfactual | Infeasible:
        """Verify, snap, and package one solver candidate.

        ``score`` is an optional precomputed ``raw_score(self.ir, x_cf)``
        (e.g. from a vectorized batch evaluation); it is recomputed whenever
        value policies modify the candidate.
        """
        if score is None:
            score = raw_score(self.ir, x_cf)
        verification = self._verify(x, x_cf, interval, score=score)
        if verification is not None:  # defensive: the GA only returns checked individuals
            return Infeasible(reason=f"heuristic solution failed verification: {verification}")
        x_cf = self._prune_changes(x, x_cf, interval)
        final_cf, snapped = self._apply_value_policies(x, x_cf, interval)
        score = raw_score(self.ir, final_cf)
        return self._result(x, final_cf, "heuristic", stats, snapped, score=score)

    def _prune_changes(
        self, x: FloatArray, x_cf: FloatArray, interval: tuple[float, float]
    ) -> FloatArray:
        """Greedily revert changes that verification proves unnecessary.

        The search's revert-to-factual mutation is stochastic, so a stalled
        run can leave residual micro-changes that cross no decision threshold
        — they cost distance without moving the score. Reverting candidates
        one at a time (cheapest change first) can only lower the objective,
        and every kept revert is re-verified in float space, so the returned
        plan keeps all its guarantees.
        """
        allow = self.compiled.allow_missing

        def effort(j: int) -> float:
            source, dest = x[j], x_cf[j]
            if math.isnan(dest):
                delta = allow[j][0]
            elif math.isnan(source):
                delta = allow[j][1]
            else:
                delta = abs(dest - source)
            return float(self.weights[j] * delta / self.sigma[j])

        changed = [
            j
            for j in range(len(x))
            if (x[j] != x_cf[j]) and not (math.isnan(x[j]) and math.isnan(x_cf[j]))
        ]
        if len(changed) < 2:  # a single change is necessary by feasibility
            return x_cf
        candidate = x_cf.copy()
        for j in sorted(changed, key=effort):
            trial = candidate.copy()
            trial[j] = x[j]
            if self._verify(x, trial, interval) is None:
                candidate = trial
        return candidate

    def _prepared_tree_arrays(self) -> tuple[TreeArrays, ...]:
        if not hasattr(self, "_prepared_trees"):
            self._prepared_trees = prepare_tree_arrays(self.ir)
        return self._prepared_trees

    def _solve_batch(
        self,
        X: FloatArray,
        tasks: Sequence[tuple[int, int]],
        interval: tuple[float, float],
        time_budget_s: float,
        sparsity_weight: float,
    ) -> list[GeneticResult]:
        """Run independent seeded searches in one parallel Rust call."""
        from treecf.backends.genetic_rust import solve_genetic_batch_rust

        if not hasattr(self, "_rust_cache"):
            self._rust_cache = {}
        return solve_genetic_batch_rust(
            self.ir,
            X,
            tasks,
            interval,
            self.compiled,
            self.sigma,
            self.weights,
            lam=sparsity_weight,
            background=self.background,
            plausibility=self._plausibility_bound(),
            time_budget_s=time_budget_s,
            cache=self._rust_cache,
        )

    def _verify(
        self,
        x: FloatArray,
        x_cf: FloatArray,
        interval: tuple[float, float],
        score: float | None = None,
    ) -> str | None:
        """Float-space re-check of target and constraints. None = OK."""
        if score is None:
            score = raw_score(self.ir, x_cf)
        if not (interval[0] <= score <= interval[1]):
            return f"score {score} outside target {interval}"
        lo, hi, _frozen = self.compiled.instance_bounds(x)  # bounds anchor at the factual x
        lo = np.where(np.isnan(lo), -math.inf, lo)
        hi = np.where(np.isnan(hi), math.inf, hi)
        for j, value in enumerate(x_cf):
            if math.isnan(value):
                if not math.isnan(x[j]) and j not in self.compiled.allow_missing:
                    return f"feature {self.ir.feature_names[j]!r} became NaN without AllowMissing"
                continue
            if not (lo[j] <= value <= hi[j]):
                return f"feature {self.ir.feature_names[j]!r} violates its bounds"

        slack = 1e-9
        for lin in self.compiled.linears:
            values = [x_cf[j] for j in lin.indices]
            if any(math.isnan(v) for v in values):
                if lin.missing_policy == "satisfied":
                    continue
                return "Linear constraint references a missing value"
            total = sum(c * v for c, v in zip(lin.coefs, values, strict=True))
            ok = (
                total <= lin.rhs + slack
                if lin.op == "<="
                else total >= lin.rhs - slack
                if lin.op == ">="
                else abs(total - lin.rhs) <= slack
            )
            if not ok:
                return f"Linear constraint violated: {lin.coefficients} {lin.op} {lin.rhs}"
        for imp in self.compiled.implications:
            if x_cf[imp.cond_index] == imp.cond_value and x_cf[imp.cons_index] != imp.cons_value:
                return "Implies constraint violated"
        for group in self.compiled.onehot_groups:
            if sum(x_cf[j] for j in group) != 1.0:
                return "OneHot constraint violated"
        if self.plausibility is not None:
            score_anomaly = self.plausibility.anomaly_score(x_cf)
            if score_anomaly > self.plausibility.max_anomaly_score + 1e-12:
                return f"anomaly score {score_anomaly:.4f} exceeds plausibility bound"
        return None

    def _plausibility_bound(self) -> tuple[EnsembleIR, float] | None:
        if self.plausibility is None:
            return None
        return self.plausibility.if_ir, self.plausibility.min_total_path

    def _apply_value_policies(
        self, x: FloatArray, x_cf: FloatArray, interval: tuple[float, float]
    ) -> tuple[FloatArray, dict[str, bool]]:
        """Snap changed values per policy inside their cells; never break validity.

        The unsnapped ``x_cf`` is already verified, so reverting offending features
        one by one is guaranteed to terminate in a valid state.
        """
        applicable = [
            (j, name, self.value_policy[name])
            for j, name in enumerate(self.ir.feature_names)
            if name in self.value_policy
            and self.value_policy[name] != "raw"
            and not math.isnan(x_cf[j])
            and x_cf[j] != x[j]
        ]
        if not applicable:
            return x_cf, {}

        cells = feature_cells(self.ir)
        lo_b, hi_b, _ = self.compiled.instance_bounds(x)
        snapped: dict[str, bool] = {}
        candidate = x_cf.copy()
        for j, name, policy in applicable:
            cell = cells[j][cell_index(cells[j], x_cf[j])]
            value = _snap(x_cf[j], policy, cell.contains, float(lo_b[j]), float(hi_b[j]))
            if value is None:
                snapped[name] = False
            else:
                candidate[j] = value
                snapped[name] = True

        # Revert snapped features one at a time until the candidate verifies.
        order = [name for name in snapped if snapped[name]]
        while self._verify(x, candidate, interval) is not None and order:
            name = order.pop()
            j = self.ir.feature_names.index(name)
            candidate[j] = x_cf[j]
            snapped[name] = False
        if self._verify(x, candidate, interval) is not None:
            return x_cf, dict.fromkeys(snapped, False)
        return candidate, snapped

    def _result(
        self,
        x: FloatArray,
        x_cf: FloatArray,
        status: str,
        stats: dict[str, object],
        snapped: dict[str, bool] | None = None,
        score: float | None = None,
    ) -> Counterfactual:
        changes: dict[str, tuple[float, float]] = {}
        distance = 0.0
        for j, name in enumerate(self.ir.feature_names):
            x_nan, cf_nan = math.isnan(x[j]), math.isnan(x_cf[j])
            if (x[j] == x_cf[j]) or (x_nan and cf_nan):
                continue
            changes[name] = (float(x[j]), float(x_cf[j]))
            if cf_nan:  # value -> NaN priced by delta_miss
                delta = self.compiled.allow_missing[j][0]
            elif x_nan:  # NaN -> value priced by delta_from_miss
                delta = self.compiled.allow_missing[j][1]
            else:
                delta = abs(x_cf[j] - x[j])
            distance += self.weights[j] * delta / self.sigma[j]
        if score is None:
            score = raw_score(self.ir, x_cf)
        return Counterfactual(
            x_cf=x_cf,
            changes=changes,
            distance=float(distance),
            n_changed=len(changes),
            score_raw=score,
            score_prob=apply_link(Link.SIGMOID, score) if self.ir.link is Link.SIGMOID else None,
            proof=status,
            solver_stats=stats,
            snapped=snapped or {},
        )

explain(x, target, backend='genetic', time_budget_s=10.0, sparsity_weight=0.0, seed=None)

Search for a counterfactual (or one per band for Target.bands).

backend="genetic" runs the bundled Rust engine (default); backend="python" runs the reference numpy implementation of the same algorithm. Every result is float-verified before being returned.

Source code in src/treecf/api.py

def explain(
    self,
    x: FloatArray,
    target: Target,
    backend: str = "genetic",
    time_budget_s: float = 10.0,
    sparsity_weight: float = 0.0,
    seed: int | None = None,
) -> Counterfactual | Infeasible | dict[str, object]:
    """Search for a counterfactual (or one per band for ``Target.bands``).

    ``backend="genetic"`` runs the bundled Rust engine (default);
    ``backend="python"`` runs the reference numpy implementation of the
    same algorithm. Every result is float-verified before being returned.
    """
    x = np.asarray(x, dtype=np.float64)
    if self.plausibility is not None and np.isnan(x).any():
        raise TreecfError("plausibility with missing factual values is not supported")
    if backend in ("genetic", "genetic-rust"):
        rust = True
    elif backend == "python":
        rust = False
    else:
        raise TreecfError(f"unknown backend {backend!r}; use 'genetic' or 'python'")

    if target.bands_spec is not None:
        results: dict[str, object] = {}
        for name, interval in target.band_intervals(self.ir.link).items():
            results[name] = self._explain_genetic(
                x, interval, time_budget_s, sparsity_weight, seed, rust=rust
            )
        return results
    interval = target.raw_interval(self.ir.link)
    return self._explain_genetic(
        x, interval, time_budget_s, sparsity_weight, seed, rust=rust
    )

explain_batch(X, target, n_per_example=1, diversity='seeds', ids=None, backend='genetic', time_budget_s=10.0, sparsity_weight=0.0, seed=0, coalitions=None, include_full=False)

Mass-produce counterfactuals for a dataset; see treecf.batch.

n_per_example alternatives per row via diversity="seeds" (distinct change-sets from different seeds, best-effort) or "lever-blocking" (freeze each plan's biggest lever; also records essential levers). diversity="coalitions" instead produces one plan per named feature group in coalitions per row (n_per_example unused; see explain_coalitions). The returned BatchResult supports save/load/for_id/to_frame.

Solves run in parallel inside the Rust engine; time_budget_s is per solve, so a solve that hits its wall-clock budget while sharing cores may stop earlier than it would sequentially (results are otherwise identical to solving row by row).

Source code in src/treecf/api.py

def explain_batch(
    self,
    X: FloatArray,
    target: Target,
    n_per_example: int = 1,
    diversity: str = "seeds",
    ids: Sequence[object] | None = None,
    backend: str = "genetic",
    time_budget_s: float = 10.0,
    sparsity_weight: float = 0.0,
    seed: int = 0,
    coalitions: Mapping[str, Sequence[str]] | None = None,
    include_full: bool = False,
) -> Any:
    """Mass-produce counterfactuals for a dataset; see ``treecf.batch``.

    ``n_per_example`` alternatives per row via ``diversity="seeds"`` (distinct
    change-sets from different seeds, best-effort) or ``"lever-blocking"``
    (freeze each plan's biggest lever; also records essential levers).
    ``diversity="coalitions"`` instead produces one plan per named feature
    group in ``coalitions`` per row (``n_per_example`` unused; see
    ``explain_coalitions``). The returned ``BatchResult`` supports
    save/load/for_id/to_frame.

    Solves run in parallel inside the Rust engine; ``time_budget_s`` is
    per solve, so a solve that hits its wall-clock budget while sharing
    cores may stop earlier than it would sequentially (results are
    otherwise identical to solving row by row).
    """
    from treecf.batch import explain_batch

    return explain_batch(
        self, X, target, n_per_example=n_per_example, diversity=diversity,
        ids=ids, backend=backend, time_budget_s=time_budget_s,
        sparsity_weight=sparsity_weight, seed=seed,
        coalitions=coalitions, include_full=include_full,
    )

explain_coalitions(x, target, coalitions, include_full=False, backend='genetic', time_budget_s=10.0, sparsity_weight=0.0, seed=None)

One counterfactual per named feature coalition (opt-in mode).

Each coalition is solved with every feature outside it frozen, so a plan only ever asks for changes within one group — grouped recourse instead of one plan that mixes unrelated levers. Coalitions may overlap; features in no coalition are never modified; an Infeasible for a coalition means that group alone cannot reach the target. include_full=True prepends an unrestricted baseline under the reserved key "(all levers)". One solve per coalition (milliseconds each); this mode is optional and never the default.

Source code in src/treecf/api.py

def explain_coalitions(
    self,
    x: FloatArray,
    target: Target,
    coalitions: Mapping[str, Sequence[str]],
    include_full: bool = False,
    backend: str = "genetic",
    time_budget_s: float = 10.0,
    sparsity_weight: float = 0.0,
    seed: int | None = None,
) -> dict[str, Counterfactual | Infeasible]:
    """One counterfactual per named feature coalition (opt-in mode).

    Each coalition is solved with every feature *outside* it frozen, so a
    plan only ever asks for changes within one group — grouped recourse
    instead of one plan that mixes unrelated levers. Coalitions may
    overlap; features in no coalition are never modified; an
    ``Infeasible`` for a coalition means that group alone cannot reach
    the target. ``include_full=True`` prepends an unrestricted baseline
    under the reserved key ``"(all levers)"``. One solve per coalition
    (milliseconds each); this mode is optional and never the default.
    """
    if target.bands_spec is not None:
        raise TreecfError(
            "Target.bands is not supported in explain_coalitions; loop bands explicitly"
        )
    normalized = _validate_coalitions(coalitions, self.ir.feature_names, include_full)
    results: dict[str, Counterfactual | Infeasible] = {}
    if include_full:
        results[_ALL_LEVERS] = self._explain_one(
            x, target, backend, time_budget_s, sparsity_weight, seed
        )
    for name, clone in self._coalition_explainers(normalized).items():
        results[name] = clone._explain_one(
            x, target, backend, time_budget_s, sparsity_weight, seed
        )
    return results

Source code in src/treecf/api.py

@dataclass(frozen=True)
class Counterfactual:
    x_cf: FloatArray
    changes: dict[str, tuple[float, float]]
    distance: float
    n_changed: int
    score_raw: float
    score_prob: float | None
    proof: str  # "heuristic" — the genetic engine never claims optimality
    solver_stats: dict[str, object] = field(default_factory=dict)
    snapped: dict[str, bool] = field(default_factory=dict)  # value_policy outcome

Source code in src/treecf/api.py

@dataclass(frozen=True)
class Infeasible:
    reason: str

Batch production

Counterfactuals for a whole dataset, addressable by row id.

Source code in src/treecf/batch.py

@dataclass(frozen=True)
class BatchResult:
    """Counterfactuals for a whole dataset, addressable by row id."""

    feature_names: tuple[str, ...]
    diversity: str
    records: tuple[BatchRecord, ...]
    essential_levers: dict[object, list[str]] = field(default_factory=dict)

    def __len__(self) -> int:
        return len(self.records)

    def __iter__(self) -> Iterator[BatchRecord]:
        return iter(self.records)

    def for_id(self, row_id: object) -> list[BatchRecord]:
        return [r for r in self.records if r.id == row_id]

    def save(self, path: str | os.PathLike[str]) -> None:
        data = {
            "feature_names": list(self.feature_names),
            "diversity": self.diversity,
            "essential_levers": {str(k): v for k, v in self.essential_levers.items()},
            "essential_lever_ids": [encode_floats(k) for k in self.essential_levers],
            "records": [
                {
                    "id": record.id,
                    "k": record.k,
                    "feasible": record.feasible,
                    "x_cf": None if record.x_cf is None else encode_floats(record.x_cf),
                    "changes": {
                        name: encode_floats(list(pair))
                        for name, pair in record.changes.items()
                    },
                    "distance": record.distance,
                    "n_changed": record.n_changed,
                    "score_raw": record.score_raw,
                    "score_prob": record.score_prob,
                    "seed": record.seed,
                    "blocked_lever": record.blocked_lever,
                    "coalition": record.coalition,
                }
                for record in self.records
            ],
        }
        with open(path, "w", encoding="utf-8") as fh:
            json.dump(data, fh)

    @classmethod
    def load(cls, path: str | os.PathLike[str]) -> BatchResult:
        with open(path, encoding="utf-8") as fh:
            data = json.load(fh)
        records = []
        for raw in data["records"]:
            records.append(
                BatchRecord(
                    id=raw["id"],
                    k=int(raw["k"]),
                    feasible=bool(raw["feasible"]),
                    x_cf=(
                        None
                        if raw["x_cf"] is None
                        else np.asarray(decode_floats(raw["x_cf"]), dtype=np.float64)
                    ),
                    changes={
                        name: tuple(decode_floats(pair))
                        for name, pair in raw["changes"].items()
                    },
                    distance=raw["distance"],
                    n_changed=raw["n_changed"],
                    score_raw=raw["score_raw"],
                    score_prob=raw["score_prob"],
                    seed=raw["seed"],
                    blocked_lever=raw["blocked_lever"],
                    coalition=raw.get("coalition"),  # absent in pre-coalition files
                )
            )
        essential_ids = [decode_floats(k) for k in data.get("essential_lever_ids", [])]
        essential_values = list(data.get("essential_levers", {}).values())
        return cls(
            feature_names=tuple(data["feature_names"]),
            diversity=data["diversity"],
            records=tuple(records),
            essential_levers=dict(zip(essential_ids, essential_values, strict=True)),
        )

    def to_frame(self) -> Any:
        """One row per (id, k), wide ``cf_<feature>`` columns (pandas, lazy import)."""
        try:
            import pandas as pd
        except ImportError as exc:  # pragma: no cover - exercised without pandas
            raise TreecfError("to_frame() requires pandas: pip install pandas") from exc
        rows = []
        for record in self.records:
            row: dict[str, object] = {
                "id": record.id,
                "k": record.k,
                "feasible": record.feasible,
                "distance": record.distance,
                "n_changed": record.n_changed,
                "score_raw": record.score_raw,
                "score_prob": record.score_prob,
                "seed": record.seed,
                "blocked_lever": record.blocked_lever,
                "coalition": record.coalition,
                "changed_features": sorted(record.changes),
            }
            for j, name in enumerate(self.feature_names):
                row[f"cf_{name}"] = (
                    float(record.x_cf[j]) if record.x_cf is not None else math.nan
                )
            rows.append(row)
        return pd.DataFrame(rows)

to_frame()

One row per (id, k), wide cf_<feature> columns (pandas, lazy import).

Source code in src/treecf/batch.py

def to_frame(self) -> Any:
    """One row per (id, k), wide ``cf_<feature>`` columns (pandas, lazy import)."""
    try:
        import pandas as pd
    except ImportError as exc:  # pragma: no cover - exercised without pandas
        raise TreecfError("to_frame() requires pandas: pip install pandas") from exc
    rows = []
    for record in self.records:
        row: dict[str, object] = {
            "id": record.id,
            "k": record.k,
            "feasible": record.feasible,
            "distance": record.distance,
            "n_changed": record.n_changed,
            "score_raw": record.score_raw,
            "score_prob": record.score_prob,
            "seed": record.seed,
            "blocked_lever": record.blocked_lever,
            "coalition": record.coalition,
            "changed_features": sorted(record.changes),
        }
        for j, name in enumerate(self.feature_names):
            row[f"cf_{name}"] = (
                float(record.x_cf[j]) if record.x_cf is not None else math.nan
            )
        rows.append(row)
    return pd.DataFrame(rows)

One counterfactual (or the infeasibility marker) for one dataset row.

Source code in src/treecf/batch.py

@dataclass(frozen=True)
class BatchRecord:
    """One counterfactual (or the infeasibility marker) for one dataset row."""

    id: object
    k: int
    feasible: bool
    x_cf: FloatArray | None
    changes: dict[str, tuple[float, float]]
    distance: float | None
    n_changed: int | None
    score_raw: float | None
    score_prob: float | None
    seed: int | None = None  # diversity="seeds": the seed that produced this plan
    blocked_lever: str | None = None  # diversity="lever-blocking": the frozen lever
    coalition: str | None = None  # diversity="coalitions": the group this plan may touch

Targets

Closed interval target, expressed in raw-score or probability space.

Target.bands builds a named ladder of intervals (rating grades); Explainer.explain then returns one result per band.

Source code in src/treecf/targets.py

@dataclass(frozen=True)
class Target:
    """Closed interval target, expressed in raw-score or probability space.

    ``Target.bands`` builds a named ladder of intervals (rating grades);
    ``Explainer.explain`` then returns one result per band.
    """

    space: str  # "raw" | "probability"
    lo: float
    hi: float
    bands_spec: tuple[tuple[str, float, float], ...] | None = None

    @classmethod
    def raw(
        cls,
        range: tuple[float, float] | None = None,
        op: str | None = None,
        value: float | None = None,
    ) -> Target:
        lo, hi = _interval_from(range, op, value, lo_limit=-math.inf, hi_limit=math.inf)
        return cls(space="raw", lo=lo, hi=hi)

    @classmethod
    def probability(
        cls,
        range: tuple[float, float] | None = None,
        op: str | None = None,
        value: float | None = None,
    ) -> Target:
        lo, hi = _interval_from(range, op, value, lo_limit=0.0, hi_limit=1.0)
        if not (0.0 <= lo < hi <= 1.0):
            raise TargetError(f"probability interval [{lo}, {hi}] must lie within [0, 1]")
        return cls(space="probability", lo=lo, hi=hi)

    @classmethod
    def bands(
        cls, bands: dict[str, tuple[float, float]], space: str = "probability"
    ) -> Target:
        if space not in ("raw", "probability"):
            raise TargetError("bands space must be 'raw' or 'probability'")
        if not bands:
            raise TargetError("bands must contain at least one named interval")
        spec = []
        for name, (lo, hi) in bands.items():
            if not lo < hi:
                raise TargetError(f"band {name!r}: empty interval [{lo}, {hi}]")
            if space == "probability" and not (0.0 <= lo < hi <= 1.0):
                raise TargetError(f"band {name!r} must lie within [0, 1]")
            spec.append((name, float(lo), float(hi)))
        first = spec[0]
        return cls(space=space, lo=first[1], hi=first[2], bands_spec=tuple(spec))

    def raw_interval(self, link: Link) -> tuple[float, float]:
        """Interval [L, U] on the raw score; probability targets require the SIGMOID link."""
        if self.space == "raw":
            return self.lo, self.hi
        if link is not Link.SIGMOID:
            raise TargetError(
                "probability target requires a SIGMOID-link model; "
                "use Target.raw for identity-link outputs"
            )
        return _logit(self.lo), _logit(self.hi)

    def band_intervals(self, link: Link) -> dict[str, tuple[float, float]]:
        assert self.bands_spec is not None
        out: dict[str, tuple[float, float]] = {}
        for name, lo, hi in self.bands_spec:
            single = Target(space=self.space, lo=lo, hi=hi)
            out[name] = single.raw_interval(link)
        return out

raw_interval(link)

Interval [L, U] on the raw score; probability targets require the SIGMOID link.

Source code in src/treecf/targets.py

def raw_interval(self, link: Link) -> tuple[float, float]:
    """Interval [L, U] on the raw score; probability targets require the SIGMOID link."""
    if self.space == "raw":
        return self.lo, self.hi
    if link is not Link.SIGMOID:
        raise TargetError(
            "probability target requires a SIGMOID-link model; "
            "use Target.raw for identity-link outputs"
        )
    return _logit(self.lo), _logit(self.hi)

Constraints

Parse "2*a - b <= 3"-style sugar into a canonical Linear object.

When feature_names is given, identifiers are validated immediately; otherwise validation happens later in compile_constraints.

Source code in src/treecf/constraints/parser.py

def constraint(text: str, feature_names: Sequence[str] | None = None) -> Linear:
    """Parse ``"2*a - b <= 3"``-style sugar into a canonical Linear object.

    When ``feature_names`` is given, identifiers are validated immediately;
    otherwise validation happens later in ``compile_constraints``.
    """
    tokens = _tokenize(text)
    parser = _Parser(text, tokens, feature_names)
    return parser.parse()

Canonical constraint objects. Frozen dataclasses; validation at compile time.

M1 subset: Freeze, Monotone, Range. Linear, Implies, OneHot, AllowMissing arrive in M2.

AllowMissing dataclass

NaN is a feasible counterfactual value for this feature.

delta_miss prices the value<->NaN transition; pass delta_from_miss for an asymmetric NaN->value cost (defaults to delta_miss).

Source code in src/treecf/constraints/objects.py

@dataclass(frozen=True)
class AllowMissing:
    """NaN is a feasible counterfactual value for this feature.

    ``delta_miss`` prices the value<->NaN transition; pass ``delta_from_miss``
    for an asymmetric NaN->value cost (defaults to ``delta_miss``).
    """

    feature: str
    delta_miss: float
    delta_from_miss: float | None = None

Equals dataclass

Binary-feature equality (used standalone or inside Implies).

Source code in src/treecf/constraints/objects.py

@dataclass(frozen=True)
class Equals:
    """Binary-feature equality (used standalone or inside Implies)."""

    feature: str
    value: float

Freeze dataclass

The feature is immutable: the counterfactual keeps the factual value.

Source code in src/treecf/constraints/objects.py

@dataclass(frozen=True)
class Freeze:
    """The feature is immutable: the counterfactual keeps the factual value."""

    feature: str

Implies dataclass

If condition holds then consequence must hold; binary features only (v0.1).

Source code in src/treecf/constraints/objects.py

@dataclass(frozen=True)
class Implies:
    """If `condition` holds then `consequence` must hold; binary features only (v0.1)."""

    condition: Equals
    consequence: Equals

Linear dataclass

Linear inter-feature constraint: sum(coef * feature) op rhs.

missing_policy resolves the constraint when a referenced feature is NaN in the counterfactual: "satisfied" (vacuously true, the default), "violated"/"forbid_missing" (the counterfactual may not use NaN there).

Source code in src/treecf/constraints/objects.py

@dataclass(frozen=True)
class Linear:
    """Linear inter-feature constraint: sum(coef * feature) op rhs.

    ``missing_policy`` resolves the constraint when a referenced feature is NaN
    in the counterfactual: "satisfied" (vacuously true, the default),
    "violated"/"forbid_missing" (the counterfactual may not use NaN there).
    """

    coefficients: dict[str, float]
    op: str  # "<=" | ">=" | "=="
    rhs: float
    missing_policy: str = "satisfied"

Monotone dataclass

The feature may only move in one direction from the factual value.

Source code in src/treecf/constraints/objects.py

@dataclass(frozen=True)
class Monotone:
    """The feature may only move in one direction from the factual value."""

    feature: str
    direction: str  # "increase" | "decrease"

OneHot dataclass

The listed binary columns sum to exactly one.

Source code in src/treecf/constraints/objects.py

@dataclass(frozen=True)
class OneHot:
    """The listed binary columns sum to exactly one."""

    features: tuple[str, ...]

Range dataclass

Hard domain bounds for the counterfactual value (inclusive).

Source code in src/treecf/constraints/objects.py

@dataclass(frozen=True)
class Range:
    """Hard domain bounds for the counterfactual value (inclusive)."""

    feature: str
    lo: float
    hi: float

Mining

Source code in src/treecf/mining.py

def suggest_constraints(
    X: FloatArray,
    feature_names: Sequence[str] | None = None,
    min_support: float = 1.0,
    top_k: int = 50,
    report_threshold: float = 0.999,
    include_ranges: bool = False,
) -> SuggestionSet:
    X = np.asarray(X, dtype=np.float64)
    n, p = X.shape
    names = list(feature_names) if feature_names is not None else [f"f{i}" for i in range(p)]
    present = ~np.isnan(X)

    binary = [
        j
        for j in range(p)
        if present[:, j].any() and np.isin(X[present[:, j], j], (0.0, 1.0)).all()
    ]
    binary_set = set(binary)

    suggestions: list[SuggestedConstraint] = []
    findings: list[DataQualityFinding] = []

    # --- pairwise order / equality (O(p^2 n), vectorized per anchor column) ---
    order_edges: dict[tuple[int, int], SuggestedConstraint] = {}
    equal_pairs: list[tuple[int, int]] = []
    for a in range(p):
        for b in range(a + 1, p):
            both = present[:, a] & present[:, b]
            n_both = int(both.sum())
            if n_both == 0:
                continue
            va, vb = X[both, a], X[both, b]
            viol_ab = int((va > vb).sum())  # violations of a <= b
            viol_ba = int((vb > va).sum())
            if viol_ab == 0 and viol_ba == 0:
                equal_pairs.append((a, b))
                suggestions.append(
                    SuggestedConstraint(
                        constraint=Linear({names[a]: 1.0, names[b]: -1.0}, op="==", rhs=0.0),
                        kind="equality",
                        support=1.0,
                        n_rows_checked=n_both,
                        n_violations=0,
                        rationale=f"{names[a]} == {names[b]} on every co-present row; "
                        "usually a redundant feature, not a constraint to impose",
                    )
                )
                continue
            for lo_idx, hi_idx, viol in ((a, b, viol_ab), (b, a, viol_ba)):
                support = 1.0 - viol / n_both
                if support >= min_support:
                    order_edges[(lo_idx, hi_idx)] = SuggestedConstraint(
                        constraint=Linear(
                            {names[lo_idx]: 1.0, names[hi_idx]: -1.0}, op="<=", rhs=0.0
                        ),
                        kind="order",
                        support=support,
                        n_rows_checked=n_both,
                        n_violations=viol,
                        evidence=_violation_evidence(X, present, lo_idx, hi_idx),
                        rationale=_order_rationale(names[lo_idx], names[hi_idx]),
                    )
                elif support >= report_threshold:
                    findings.append(
                        DataQualityFinding(
                            kind="near_invariant",
                            description=f"{names[lo_idx]} <= {names[hi_idx]} holds on "
                            f"{support:.4%} of rows — likely an ETL defect",
                            support=support,
                            n_rows_checked=n_both,
                            n_violations=viol,
                            evidence=_violation_evidence(X, present, lo_idx, hi_idx),
                        )
                    )

    # equality-class collapse, then transitive reduction of the <= graph
    representative = _union_find(p, equal_pairs)
    rep_edges = {
        (representative[a], representative[b])
        for (a, b) in order_edges
        if representative[a] != representative[b]
    }
    reduced = transitive_reduction(rep_edges)
    for (a, b), suggestion in order_edges.items():
        edge = (representative[a], representative[b])
        if edge in reduced and edge[0] != edge[1]:
            suggestions.append(suggestion)
            reduced.discard(edge)  # one edge per class pair

    # --- binary implications A=1 => B=1 ---
    for a in binary:
        a_is_one = present[:, a] & (X[:, a] == 1.0)
        if not a_is_one.any():
            continue
        for b in binary:
            if a == b:
                continue
            checked = a_is_one & present[:, b]
            if not checked.any():
                continue
            if (X[checked, b] == 1.0).all():
                suggestions.append(
                    SuggestedConstraint(
                        constraint=Implies(Equals(names[a], 1.0), Equals(names[b], 1.0)),
                        kind="implication",
                        support=1.0,
                        n_rows_checked=int(checked.sum()),
                        n_violations=0,
                        rationale=f"{names[a]}=1 always co-occurs with {names[b]}=1",
                    )
                )

    # --- one-hot groups: exclusivity components with row sum == 1 ---
    complete_binary = [j for j in binary if present[:, j].all()]
    for component in _exclusivity_components(X, complete_binary):
        if len(component) < 2:
            continue
        if np.all(X[:, component].sum(axis=1) == 1.0):
            suggestions.append(
                SuggestedConstraint(
                    constraint=OneHot(tuple(names[j] for j in component)),
                    kind="onehot",
                    support=1.0,
                    n_rows_checked=n,
                    n_violations=0,
                    rationale="binary columns with row sum identically 1",
                )
            )

    # --- missingness links miss(A) => miss(B) ---
    for a in range(p):
        miss_a = ~present[:, a]
        if not miss_a.any():
            continue
        for b in range(p):
            if a == b or present[:, b].all():
                continue
            if (~present[miss_a, b]).all():
                both_ways = bool((~present[~present[:, b], a]).all())
                suggestions.append(
                    SuggestedConstraint(
                        constraint=None,
                        kind="missing_link",
                        support=1.0,
                        n_rows_checked=int(miss_a.sum()),
                        n_violations=0,
                        rationale=(
                            f"miss({names[a]}) {'<=>' if both_ways else '=>'} miss({names[b]}); "
                            "consider joint AllowMissing / missing_policy"
                        ),
                    )
                )
            if bool((~present[~present[:, b], a]).all()):
                break  # symmetric link already reported from this anchor

    # --- integer-valuedness -> value_policy suggestion ---
    for j in range(p):
        col = X[present[:, j], j]
        if len(col) and j not in binary_set and np.all(col == np.round(col)):
            suggestions.append(
                SuggestedConstraint(
                    constraint=None,
                    kind="integer",
                    support=1.0,
                    n_rows_checked=len(col),
                    n_violations=0,
                    rationale=(
                        f"{names[j]} is integer-valued; "
                        f'value_policy={{"{names[j]}": "integer"}}'
                    ),
                )
            )

    if include_ranges:
        for j in range(p):
            col = X[present[:, j], j]
            if len(col):
                lo, hi = np.percentile(col, [1, 99])
                pad = 0.1 * (hi - lo)
                suggestions.append(
                    SuggestedConstraint(
                        constraint=None,
                        kind="range",
                        support=1.0,
                        n_rows_checked=len(col),
                        n_violations=0,
                        rationale=f"observed 1-99% range [{lo:.4g}, {hi:.4g}] padded by {pad:.4g}",
                    )
                )

    suggestions.sort(key=_rank_key, reverse=True)
    return SuggestionSet(suggestions=tuple(suggestions[:top_k]), findings=tuple(findings))

Source code in src/treecf/mining.py

@dataclass(frozen=True)
class SuggestedConstraint:
    constraint: Constraint | None  # None for advisory kinds (missing_link, integer)
    kind: str  # "order" | "equality" | "implication" | "onehot" | "missing_link" | "integer"
    support: float
    n_rows_checked: int
    n_violations: int
    evidence: list[dict[str, object]] = field(default_factory=list)
    rationale: str = ""

    def as_code(self) -> str:
        tail = f"  # support={self.support:.4f}, n={self.n_rows_checked}"
        if self.kind == "order" and isinstance(self.constraint, Linear):
            coeffs = self.constraint.coefficients
            smaller = max(coeffs, key=lambda k: coeffs[k])
            larger = min(coeffs, key=lambda k: coeffs[k])
            return f'constraint("{smaller} <= {larger}")' + tail
        if self.kind == "equality" and isinstance(self.constraint, Linear):
            a, b = list(self.constraint.coefficients)
            return f'# equality: {a} == {b} — likely a redundant feature' + tail
        if self.kind == "implication" and isinstance(self.constraint, Implies):
            c = self.constraint
            return (
                f'Implies(Equals("{c.condition.feature}", {c.condition.value:g}), '
                f'Equals("{c.consequence.feature}", {c.consequence.value:g}))' + tail
            )
        if self.kind == "onehot" and isinstance(self.constraint, OneHot):
            inner = ", ".join(f'"{f}"' for f in self.constraint.features)
            return f"OneHot(({inner}))" + tail
        return f"# {self.kind}: {self.rationale}" + tail

Source code in src/treecf/mining.py

@dataclass(frozen=True)
class DataQualityFinding:
    kind: str  # "near_invariant"
    description: str
    support: float
    n_rows_checked: int
    n_violations: int
    evidence: list[dict[str, object]] = field(default_factory=list)

Plausibility

Source code in src/treecf/plausibility.py

@dataclass(frozen=True)
class Plausibility:
    if_ir: EnsembleIR
    max_anomaly_score: float

    @classmethod
    def isolation_forest(
        cls, model_or_ir: object, max_anomaly_score: float = 0.55
    ) -> Plausibility:
        if not 0.0 < max_anomaly_score < 1.0:
            raise TreecfError("max_anomaly_score must lie in (0, 1)")
        if isinstance(model_or_ir, EnsembleIR):
            if_ir = model_or_ir
        else:
            from treecf.ir.parsers.sklearn import parse_isolation_forest

            if_ir = parse_isolation_forest(model_or_ir)
        return cls(if_ir=if_ir, max_anomaly_score=max_anomaly_score)

    @property
    def normalizer(self) -> float:
        """c(n) for the forest's subsample size."""
        from treecf.ir.parsers.sklearn import _avg_path

        return _avg_path(float(self.if_ir.meta["max_samples"]))  # type: ignore[arg-type]

    @property
    def min_total_path(self) -> float:
        """Feasibility bound: sum_t h_t(x') >= -T * c(n) * log2(theta)."""
        n_trees = len(self.if_ir.trees)
        return -n_trees * self.normalizer * math.log2(self.max_anomaly_score)

    def anomaly_score(self, x: FloatArray) -> float:
        total = raw_score(self.if_ir, np.asarray(x, dtype=np.float64))
        mean_path = total / len(self.if_ir.trees)
        return float(2.0 ** (-mean_path / self.normalizer))

min_total_path property

Feasibility bound: sum_t h_t(x') >= -T * c(n) * log2(theta).

normalizer property

c(n) for the forest's subsample size.

Visualization

Counterfactual visualizations. matplotlib lives behind the [viz] extra.

plot_alternatives(results, explainer=None, ax=None)

Overlaid dumbbells: every alternative plan's changes for one instance.

Accepts a sequence of Counterfactual objects or feasible BatchRecord entries, or a mapping of outcomes as returned by explain_coalitions (keys become legend labels; Infeasible values are skipped). Each plan keeps one color across all its changes — meant for a handful of alternatives for the same row (at most 10). With explainer, changes are plotted as standardized deltas from the factual (Δ/σ), so features of different scales share one axis; without, raw values are shown with gray factual dots.

Source code in src/treecf/viz.py

def plot_alternatives(results: Any, explainer: Any = None, ax: Any = None) -> Any:
    """Overlaid dumbbells: every alternative plan's changes for one instance.

    Accepts a sequence of ``Counterfactual`` objects or feasible
    ``BatchRecord`` entries, or a mapping of outcomes as returned by
    ``explain_coalitions`` (keys become legend labels; ``Infeasible`` values
    are skipped). Each plan keeps one color across all its changes — meant
    for a handful of alternatives for the same row (at most 10). With
    ``explainer``, changes are plotted as standardized deltas from the
    factual (Δ/σ), so features of different scales share one axis; without,
    raw values are shown with gray factual dots.
    """
    plt = _import_pyplot()
    plans = _plans_with_labels(results)
    if not plans:
        raise TreecfError("no feasible plans to plot")
    if len(plans) > 10:
        raise TreecfError("plot_alternatives compares at most 10 plans")
    sigma: dict[str, float] = {}
    if explainer is not None:
        sigma = {
            name: float(s)
            for name, s in zip(explainer.ir.feature_names, explainer.sigma, strict=True)
        }
    frequency: dict[str, int] = {}
    for _, plan in plans:
        for name in plan.changes:
            frequency[name] = frequency.get(name, 0) + 1
    features = sorted(frequency, key=lambda name: (-frequency[name], name))
    slots = {name: i for i, name in enumerate(features)}

    if ax is None:
        _, ax = plt.subplots(figsize=(7, 0.8 * max(2, len(features))))
    step = min(0.18, 0.7 / len(plans))
    for p, (plan_name, plan) in enumerate(plans):
        color = f"C{p}"
        offset = (p - (len(plans) - 1) / 2) * step
        base = plan_name if plan_name is not None else f"plan {p + 1}"
        label: str | None = f"{base} (J={plan.distance:.3g})"
        for name, (source, dest) in plan.changes.items():
            y = slots[name] + offset
            if math.isnan(dest) or math.isnan(source):
                anchor = source if math.isnan(dest) else dest
                if explainer is not None:
                    anchor = 0.0
                ax.plot([anchor], [y], "o", color=color, markersize=5, label=label)
                ax.annotate(
                    "-> NaN" if math.isnan(dest) else "NaN ->",
                    xy=(anchor, y), xytext=(6, 0), textcoords="offset points",
                    va="center", color="tab:red", fontsize=9,
                )
            else:
                if explainer is not None:
                    start, end = 0.0, (dest - source) / sigma[name]
                else:
                    start, end = source, dest
                ax.plot([start, end], [y, y], "-", color=color, alpha=0.5, zorder=1)
                ax.plot([start], [y], "o", color="tab:gray", markersize=4)
                ax.plot([end], [y], "o", color=color, markersize=5, label=label)
            label = None  # one legend entry per plan
    if explainer is not None:
        ax.axvline(0.0, color="0.6", linestyle="--", linewidth=1)
        ax.set_xlabel("standardized change from factual (Δ/σ)")
    else:
        ax.set_xlabel("feature value (gray = factual)")
    ax.set_yticks(range(len(features)), features)
    ax.invert_yaxis()
    ax.set_title(f"{len(plans)} alternative plan(s) for one instance")
    ax.legend(loc="best")
    return ax

plot_changes(cf, ax=None)

Dumbbell chart of per-feature changes (from -> to); NaN transitions annotated.

Source code in src/treecf/viz.py

def plot_changes(cf: Counterfactual, ax: Any = None) -> Any:
    """Dumbbell chart of per-feature changes (from -> to); NaN transitions annotated."""
    plt = _import_pyplot()
    if ax is None:
        _, ax = plt.subplots(figsize=(7, 0.6 * max(2, len(cf.changes))))
    names = list(cf.changes)
    labeled = False
    for i, name in enumerate(names):
        source, target = cf.changes[name]
        if math.isnan(target) or math.isnan(source):
            anchor = source if math.isnan(target) else target
            ax.plot([anchor], [i], "o", color="tab:gray")
            ax.annotate(
                "-> NaN" if math.isnan(target) else "NaN ->",
                xy=(anchor, i),
                xytext=(6, 0),
                textcoords="offset points",
                va="center",
                color="tab:red",
            )
            continue
        ax.plot([source, target], [i, i], "-", color="tab:gray", zorder=1)
        factual_label = None if labeled else "factual"
        cf_label = None if labeled else "counterfactual"
        ax.plot([source], [i], "o", color="tab:gray", label=factual_label)
        ax.plot([target], [i], "o", color="tab:blue", label=cf_label)
        labeled = True
    ax.set_yticks(range(len(names)), names)
    ax.set_xlabel("feature value")
    ax.set_title(f"{cf.n_changed} change(s), distance {cf.distance:.3g} ({cf.proof})")
    if labeled:
        ax.legend(loc="best")
    return ax

plot_counterfactuals(results, ax=None)

Changed-feature matrix comparing diverse counterfactuals.

Source code in src/treecf/viz.py

def plot_counterfactuals(results: Sequence[Counterfactual], ax: Any = None) -> Any:
    """Changed-feature matrix comparing diverse counterfactuals."""
    plt = _import_pyplot()
    features = sorted({name for cf in results for name in cf.changes})
    if ax is None:
        _, ax = plt.subplots(figsize=(1.0 + 0.8 * len(features), 0.8 + 0.5 * len(results)))
    matrix = [[1.0 if f in cf.changes else 0.0 for f in features] for cf in results]
    ax.imshow(matrix, cmap="Blues", aspect="auto", vmin=0.0, vmax=1.0)
    ax.set_xticks(range(len(features)), features, rotation=45, ha="right")
    ax.set_yticks(
        range(len(results)),
        [f"#{i + 1} (J={cf.distance:.3g})" for i, cf in enumerate(results)],
    )
    ax.set_title("changed features per counterfactual")
    return ax

plot_effort(explainer, cf, ax=None)

Cost-space companion: how the distance J splits across the changes.

Source code in src/treecf/viz.py

def plot_effort(explainer: Any, cf: Counterfactual, ax: Any = None) -> Any:
    """Cost-space companion: how the distance J splits across the changes."""
    plt = _import_pyplot()
    contributions = sorted(
        _change_effort(explainer, cf.changes).items(), key=lambda pair: pair[1], reverse=True
    )

    if ax is None:
        _, ax = plt.subplots(figsize=(7, 0.6 * max(2, len(contributions)) + 0.8))
    labels = [name for name, _ in contributions]
    efforts = [effort for _, effort in contributions]
    ax.barh(range(len(labels)), efforts, color="tab:blue", height=0.6)
    for i, effort in enumerate(efforts):
        ax.annotate(
            f"{effort:.3g}", xy=(effort, i), xytext=(4, 0),
            textcoords="offset points", va="center", fontsize=9,
        )
    ax.set_yticks(range(len(labels)), labels)
    ax.invert_yaxis()
    ax.set_xlabel("effort contribution (w·|Δ|/σ)")
    ax.set_title(f"where the effort goes — total J = {cf.distance:.3g}")
    return ax

plot_ladder(bands_result, ax=None)

Cost of reaching each rating band (Target.bands): the price of every grade.

Source code in src/treecf/viz.py

def plot_ladder(bands_result: Mapping[str, object], ax: Any = None) -> Any:
    """Cost of reaching each rating band (Target.bands): the price of every grade."""
    plt = _import_pyplot()
    if ax is None:
        _, ax = plt.subplots(figsize=(1.5 + 0.9 * len(bands_result), 4))
    names = list(bands_result)
    heights = []
    for name in names:
        outcome = bands_result[name]
        heights.append(outcome.distance if isinstance(outcome, Counterfactual) else 0.0)
    bars = ax.bar(names, heights, color="tab:blue")
    for bar, name in zip(bars, names, strict=True):
        outcome = bands_result[name]
        if isinstance(outcome, Infeasible):
            ax.text(
                bar.get_x() + bar.get_width() / 2,
                0.02,
                "infeasible",
                ha="center",
                va="bottom",
                rotation=90,
                color="tab:red",
            )
    ax.set_xticks(range(len(names)), names)
    ax.set_ylabel("distance J")
    ax.set_title("cost of reaching each band")
    return ax

plot_tradeoff(results, target=None, ax=None)

Cost vs achieved score for alternative plans of one instance.

One dot per plan: x = distance J, y = the achieved probability (sigmoid models) or raw score. target draws the interval bounds the plans had to reach. Accepts a sequence of Counterfactual objects or feasible BatchRecord entries, or a mapping as returned by explain_coalitions (keys label the dots; Infeasible skipped).

Source code in src/treecf/viz.py

def plot_tradeoff(results: Any, target: Any = None, ax: Any = None) -> Any:
    """Cost vs achieved score for alternative plans of one instance.

    One dot per plan: x = distance J, y = the achieved probability (sigmoid
    models) or raw score. ``target`` draws the interval bounds the plans had
    to reach. Accepts a sequence of ``Counterfactual`` objects or feasible
    ``BatchRecord`` entries, or a mapping as returned by
    ``explain_coalitions`` (keys label the dots; ``Infeasible`` skipped).
    """
    plt = _import_pyplot()
    plans = _plans_with_labels(results)
    if not plans:
        raise TreecfError("no feasible plans to plot")
    prob_space = all(plan.score_prob is not None for _, plan in plans)

    if ax is None:
        _, ax = plt.subplots(figsize=(6, 4))
    for p, (plan_name, plan) in enumerate(plans):
        score = plan.score_prob if prob_space else plan.score_raw
        ax.plot([plan.distance], [score], "o", color=f"C{p}", markersize=8)
        ax.annotate(
            plan_name if plan_name is not None else f"{p + 1}",
            xy=(plan.distance, score), xytext=(6, 4),
            textcoords="offset points", fontsize=9,
        )
    if target is not None:
        for bound in _target_bounds(target, prob_space):
            ax.axhline(bound, color="tab:red", linewidth=1)
    ax.set_xlabel("distance J (effort)")
    ax.set_ylabel("model probability" if prob_space else "raw score")
    ax.set_title("what each plan costs, and what it buys")
    return ax

plot_waterfall(explainer, cf, target=None, ax=None)

SHAP-style waterfall: exact score deltas of the counterfactual's changes.

Starts at the factual score, applies the changes one at a time (largest single effect first), each bar being the EXACT score delta from that change (recomputed through the IR — endpoints are exact; per-bar attribution is sequential and therefore order-dependent, like any sequential decomposition). Sigmoid-link models are plotted in probability space.

Source code in src/treecf/viz.py

def plot_waterfall(explainer: Any, cf: Counterfactual, target: Any = None, ax: Any = None) -> Any:
    """SHAP-style waterfall: exact score deltas of the counterfactual's changes.

    Starts at the factual score, applies the changes one at a time (largest
    single effect first), each bar being the EXACT score delta from that change
    (recomputed through the IR — endpoints are exact; per-bar attribution is
    sequential and therefore order-dependent, like any sequential decomposition).
    Sigmoid-link models are plotted in probability space.
    """
    import numpy as np

    from treecf.ir.evaluate import apply_link, raw_score
    from treecf.ir.model import Link

    plt = _import_pyplot()
    ir = explainer.ir
    index = {name: j for j, name in enumerate(ir.feature_names)}

    x = cf.x_cf.copy()
    for name, (source, _) in cf.changes.items():
        x[index[name]] = source

    def single_delta(name: str) -> float:
        probe = x.copy()
        probe[index[name]] = cf.changes[name][1]
        return raw_score(ir, probe) - raw_score(ir, x)

    order = sorted(cf.changes, key=lambda f: abs(single_delta(f)), reverse=True)

    sigmoid = ir.link is Link.SIGMOID
    to_display = (lambda s: apply_link(Link.SIGMOID, s)) if sigmoid else (lambda s: s)

    current = x.copy()
    scores = [to_display(raw_score(ir, current))]
    for name in order:
        current[index[name]] = cf.changes[name][1]
        scores.append(to_display(raw_score(ir, current)))

    if ax is None:
        _, ax = plt.subplots(figsize=(7, 0.7 * max(2, len(order)) + 1))
    for i, _name in enumerate(order):
        before, after = scores[i], scores[i + 1]
        delta = after - before
        color = "tab:blue" if delta < 0 else "tab:orange"
        ax.barh(i, delta, left=before, color=color, height=0.6)
        ax.plot([after, after], [i, i + 1], color="0.6", linestyle=":", linewidth=1)
        ax.annotate(
            f"{delta:+.4g}",
            xy=(max(before, after), i),
            xytext=(4, 0),
            textcoords="offset points",
            va="center",
            fontsize=9,
        )
    ax.axvline(scores[0], color="0.4", linestyle="--", linewidth=1)
    ax.text(scores[0], -0.55, f"f(x) = {scores[0]:.4g}", ha="center", va="top", fontsize=9)
    ax.axvline(scores[-1], color="tab:green", linestyle="--", linewidth=1)
    ax.text(
        scores[-1], len(order) - 0.3, f"f(x') = {scores[-1]:.4g}",
        ha="center", va="bottom", fontsize=9, color="tab:green",
    )
    if target is not None:
        for bound in target.raw_interval(ir.link):
            if np.isfinite(bound):
                ax.axvline(to_display(bound), color="tab:red", linewidth=1)
    ax.set_yticks(range(len(order)), order)
    ax.invert_yaxis()  # largest effect on top, like SHAP
    ax.set_xlabel("model probability" if sigmoid else "raw score")
    ax.set_title("what moves the score (sequential, exact)")
    if sigmoid:
        low = min(0.0, min(scores))
        high = max(1.0, max(scores))
        ax.set_xlim(low - 0.02, min(high + 0.05, 1.05))
    return ax

Batch-level counterfactual visualizations. matplotlib lives behind the [viz] extra.

Every function consumes a BatchResult. k=0 (the default) keeps each row's best plan; k=None keeps every feasible plan, so shares are per plan, not per row.

plot_batch_deltas(batch, explainer=None, k=0, top_n=10, ax=None)

Strip plot of actual deltas (to − from) per feature, top-N most-changed.

One jittered dot per plan, a median tick per feature; NaN transitions are counted in a per-feature annotation instead of plotted. With explainer, deltas are divided by the per-feature normalizer sigma so features of different scales share one axis.

Source code in src/treecf/viz_batch.py

def plot_batch_deltas(
    batch: BatchResult,
    explainer: Any = None,
    k: int | None = 0,
    top_n: int = 10,
    ax: Any = None,
) -> Any:
    """Strip plot of actual deltas (to − from) per feature, top-N most-changed.

    One jittered dot per plan, a median tick per feature; NaN transitions are
    counted in a per-feature annotation instead of plotted. With ``explainer``,
    deltas are divided by the per-feature normalizer sigma so features of
    different scales share one axis.
    """
    plt = _import_pyplot()
    import numpy as np

    if explainer is not None and tuple(explainer.ir.feature_names) != batch.feature_names:
        raise TreecfError("explainer and batch describe different feature spaces")
    selected = _select_records(batch, k)
    sigma = {name: 1.0 for name in batch.feature_names}
    if explainer is not None:
        sigma = dict(zip(batch.feature_names, (float(s) for s in explainer.sigma), strict=True))
    deltas: dict[str, list[float]] = {}
    nan_counts: Counter[str] = Counter()
    totals: Counter[str] = Counter()
    for record in selected:
        for name, (source, dest) in record.changes.items():
            totals[name] += 1
            if math.isnan(source) or math.isnan(dest):
                nan_counts[name] += 1
            else:
                deltas.setdefault(name, []).append((dest - source) / sigma[name])
    order = sorted(totals, key=lambda name: (-totals[name], name))[:top_n]

    if ax is None:
        _, ax = plt.subplots(figsize=(7, 0.6 * max(2, len(order))))
    rng = np.random.default_rng(0)  # fixed jitter: figures stay deterministic
    for i, name in enumerate(order):
        values = deltas.get(name, [])
        if values:
            jitter = rng.uniform(-0.15, 0.15, len(values))
            ax.plot(values, i + jitter, "o", color="tab:blue", alpha=0.6, markersize=4)
            ax.plot([float(np.median(values))], [i], "|", color="tab:orange", markersize=14)
        if nan_counts[name]:
            ax.annotate(
                f"→NaN ×{nan_counts[name]}", xy=(1.0, i),
                xycoords=("axes fraction", "data"), xytext=(-4, 0),
                textcoords="offset points", ha="right", va="center",
                color="tab:red", fontsize=9,
            )
    ax.axvline(0.0, color="0.6", linestyle="--", linewidth=1)
    ax.set_yticks(range(len(order)), order)
    ax.invert_yaxis()
    ax.set_xlabel("delta (to − from)" if explainer is None else "standardized delta (Δ/σ)")
    ax.set_title(f"how far the levers move ({len(selected)} plan(s))")
    return ax

plot_batch_levers(batch, k=0, normalize=True, top_n=20, show_essential=True, ax=None)

Horizontal stacked bars: share of plans changing each feature, by direction.

Increases, decreases, and NaN transitions stack per feature, ordered by how often the feature is used. For diversity="lever-blocking" results, features recorded as essential levers are annotated with their count.

Source code in src/treecf/viz_batch.py

def plot_batch_levers(
    batch: BatchResult,
    k: int | None = 0,
    normalize: bool = True,
    top_n: int = 20,
    show_essential: bool = True,
    ax: Any = None,
) -> Any:
    """Horizontal stacked bars: share of plans changing each feature, by direction.

    Increases, decreases, and NaN transitions stack per feature, ordered by how
    often the feature is used. For ``diversity="lever-blocking"`` results,
    features recorded as essential levers are annotated with their count.
    """
    plt = _import_pyplot()
    selected = _select_records(batch, k)
    increase: Counter[str] = Counter()
    decrease: Counter[str] = Counter()
    to_nan: Counter[str] = Counter()
    for record in selected:
        for name, (source, dest) in record.changes.items():
            if math.isnan(source) or math.isnan(dest):
                to_nan[name] += 1
            elif dest > source:
                increase[name] += 1
            else:
                decrease[name] += 1
    total = increase + decrease + to_nan
    order = sorted(total, key=lambda name: (-total[name], name))[:top_n]
    scale = 1.0 / len(selected) if normalize else 1.0

    if ax is None:
        _, ax = plt.subplots(figsize=(7, 0.6 * max(2, len(order))))
    positions = range(len(order))
    left = [0.0] * len(order)
    parts = [("increase", increase, "tab:orange"), ("decrease", decrease, "tab:blue"),
             ("NaN", to_nan, "tab:gray")]
    for label, counter, color in parts:
        widths = [counter[name] * scale for name in order]
        if not any(widths):
            continue
        ax.barh(positions, widths, left=left, height=0.6, color=color, label=label)
        left = [acc + w for acc, w in zip(left, widths, strict=True)]

    essential: Counter[str] = Counter()
    if show_essential and batch.diversity == "lever-blocking":
        essential = Counter(
            lever for levers in batch.essential_levers.values() for lever in levers
        )
    for i, name in enumerate(order):
        if essential[name]:
            ax.annotate(
                f"essential ×{essential[name]}", xy=(left[i], i), xytext=(4, 0),
                textcoords="offset points", va="center", color="tab:red", fontsize=9,
            )
    ax.set_yticks(positions, order)
    ax.invert_yaxis()
    ax.set_xlabel("fraction of plans" if normalize else "plans")
    ax.set_title(f"levers used across {len(selected)} plan(s)")
    ax.legend(loc="best")
    return ax

plot_batch_matrix(batch, explainer=None, k=0, sort_rows=True, max_row_labels=30, ax=None)

Plans × features heatmap: binary changes, or effort-shaded with an explainer.

With explainer, each cell shows the change's effort w·|Δ|/σ (NaN legs priced via AllowMissing); without, cells mark changed features like plot_counterfactuals. Rows sort by distance; columns by how often the feature is changed.

Source code in src/treecf/viz_batch.py

def plot_batch_matrix(
    batch: BatchResult,
    explainer: Any = None,
    k: int | None = 0,
    sort_rows: bool = True,
    max_row_labels: int = 30,
    ax: Any = None,
) -> Any:
    """Plans × features heatmap: binary changes, or effort-shaded with an explainer.

    With ``explainer``, each cell shows the change's effort ``w·|Δ|/σ`` (NaN
    legs priced via ``AllowMissing``); without, cells mark changed features
    like ``plot_counterfactuals``. Rows sort by distance; columns by how often
    the feature is changed.
    """
    plt = _import_pyplot()
    import numpy as np

    if explainer is not None and tuple(explainer.ir.feature_names) != batch.feature_names:
        raise TreecfError("explainer and batch describe different feature spaces")
    selected = _select_records(batch, k)
    if sort_rows:
        selected.sort(key=lambda record: record.distance or 0.0)
    frequency = Counter(name for record in selected for name in record.changes)
    features = sorted(frequency, key=lambda name: (-frequency[name], name))

    matrix = np.zeros((len(selected), len(features)))
    for i, record in enumerate(selected):
        row_values = (
            {name: 1.0 for name in record.changes}
            if explainer is None
            else _change_effort(explainer, record.changes)
        )
        for jf, name in enumerate(features):
            matrix[i, jf] = row_values.get(name, 0.0)

    if ax is None:
        height = 0.8 + min(0.3 * max(2, len(selected)), 6.0)
        _, ax = plt.subplots(figsize=(1.0 + 0.8 * len(features), height))
    if explainer is None:
        vmax = 1.0
    else:
        # robust ceiling: one extreme change must not wash out the rest
        positive = matrix[matrix > 0]
        vmax = max(float(np.percentile(positive, 95)) if positive.size else 0.0, 1e-12)
    ax.imshow(matrix, cmap="Blues", aspect="auto", vmin=0.0, vmax=vmax)
    ax.set_xticks(range(len(features)), features, rotation=45, ha="right")
    if len(selected) <= max_row_labels:
        labels = [
            f"{r.id} (J={r.distance:.3g})" + (f" k={r.k}" if k is None else "")
            for r in selected
        ]
        ax.set_yticks(range(len(selected)), labels)
    else:
        ax.set_yticks([])
        ax.set_ylabel(f"{len(selected)} plans")
    ax.set_title(
        "effort per change (w·|Δ|/σ)" if explainer is not None else "changed features per plan"
    )
    return ax

plot_batch_summary(batch, k=0, axs=None)

Three-panel batch overview: plan cost, sparsity, and feasibility.

Creates its own figure when axs is None and returns the array of three axes (unlike the single-axes functions, which return one ax).

Source code in src/treecf/viz_batch.py

def plot_batch_summary(batch: BatchResult, k: int | None = 0, axs: Any = None) -> Any:
    """Three-panel batch overview: plan cost, sparsity, and feasibility.

    Creates its own figure when ``axs`` is None and returns the array of three
    axes (unlike the single-axes functions, which return one ``ax``).
    """
    plt = _import_pyplot()
    ids_all = {record.id for record in batch.records}
    if not ids_all:
        raise TreecfError("empty batch")
    ids_ok = {record.id for record in batch.records if record.feasible}
    selected = [r for r in batch.records if r.feasible and (k is None or r.k == k)]

    own_figure = axs is None
    if own_figure:
        _, axs = plt.subplots(1, 3, figsize=(11, 3.2), constrained_layout=True)
    distances = [record.distance for record in selected if record.distance is not None]
    if distances:
        axs[0].hist(distances, bins="auto", color="tab:blue")
    else:
        axs[0].text(
            0.5, 0.5, "no feasible plans", ha="center", va="center",
            transform=axs[0].transAxes, color="tab:red",
        )
    axs[0].set_xlabel("distance J")
    axs[0].set_title("plan cost")

    sparsity = Counter(record.n_changed for record in selected)
    counts = sorted((n, c) for n, c in sparsity.items() if n is not None)
    if counts:
        axs[1].bar([n for n, _ in counts], [c for _, c in counts], color="tab:blue")
        axs[1].set_xticks([n for n, _ in counts])
    axs[1].set_xlabel("features changed")
    axs[1].set_title("sparsity")

    axs[2].bar(
        ["feasible", "infeasible"],
        [len(ids_ok), len(ids_all) - len(ids_ok)],
        color=["tab:blue", "tab:red"],
    )
    axs[2].set_title(f"{len(ids_ok) / len(ids_all):.0%} of rows solvable")

    if own_figure:
        axs[0].figure.suptitle(
            f"batch summary — {len(ids_all)} rows, diversity={batch.diversity!r}"
        )
    return axs