feat: SOTA AutoML engine, serve-safe GenAI, and publishable-grade benchmarks

.github/workflows/ci.yml

@@ -35,6 +35,9 @@ jobs:
 
   test:
     runs-on: ubuntu-latest
+    env:
+      OMP_NUM_THREADS: "1"
+      OPENBLAS_NUM_THREADS: "1"
     strategy:
       matrix:
         python-version: ["3.13"]

CONTRIBUTING.md

@@ -32,6 +32,12 @@ uv run pytest                          # tests (integration/nightly excluded by
 - **`-m integration`** runs network/heavy tests (OpenML, HuggingFace downloads).
 - **`-m nightly`** runs long-running suites (full AMLB, GPU).
 
+> **OpenMP note (macOS):** the suite must run with the OpenMP/BLAS pools pinned to one thread.
+> `tests/conftest.py` sets `OMP_NUM_THREADS=1` / `OPENBLAS_NUM_THREADS=1` before any ML import
+> (CI sets the same), because xgboost / lightgbm / scikit-learn each bundle a libomp on macOS and
+> the duplicate runtimes deadlock in an OpenMP join barrier (`__kmp_join_barrier`) at teardown.
+> If you run scripts/benchmarks that import the ML stack outside pytest, export the same variables.
+
 ## Conventions
 
 - **CalVer** `YY.MM.PATCH`; the version lives in `src/fireflyframework_datascience/_version.py`.

README.md

@@ -68,12 +68,17 @@ science actually delivers business value:
 - **Production-ready** — serving, data validation, lineage and real benchmarks are built in, not bolted on.
 
 **Proven, not promised — unbiased and significance-tested.** Under **nested cross-validation** (no
-selection bias), Firefly's AutoML significantly beats a single LogisticRegression (Δ +0.029, *p* = 0.046)
-and a single XGBoost (Δ +0.030, ***p* = 7.5e-6**), and is statistically on par with RandomForest —
-adapting per dataset, up to **+0.15** on non-linear `phoneme`. With a real LLM (`claude-haiku-4-5`),
-governed GenAI feature engineering adds a **significant +0.021** lift on a linear model (*p* = 0.0039) by
-rediscovering a withheld driver (`revenue = price × units`) from the schema — and the cost/benefit gate
-guarantees it never regresses, at **< $0.01**. Every number is reproducible — see the
+selection bias) against five systems, Firefly's AutoML earns the **best mean rank (1.60 of 6, Friedman
+*p* = 0.034)** — it beats a **tuned** XGBoost at a matched trial budget on **5/5 datasets** (one-sided
+Wilcoxon *p* = 0.031), an untuned XGBoost 5/5 (*p* = 0.031), and FLAML on 4/5 at matched wall-clock,
+adapting per dataset up to **+0.15** on non-linear `phoneme`. At *n* = 5 datasets the power is limited and
+the Holm-adjusted pairwise tests are not significant — we report the per-dataset statistics in full rather
+than lean on one headline. With a real LLM (`claude-haiku-4-5`), governed GenAI feature engineering adds a
+**significant +0.021** lift on a linear model (*p* = 0.0039) by rediscovering a withheld driver
+(`revenue = price × units`) from the schema — and the cost/benefit gate
+keeps only measured train-CV wins (its no-regression guarantee is scoped to the train-CV metric it
+optimizes; holdout results are reported as measured), at a token cost taken from metering. Every number is
+reproducible — see the
 [benchmark results](benchmarks/RESULTS.md).
 
 > 📄 **The whole story in one document:** **[The Complete Guide (PDF)](docs/brief/firefly-datascience-complete-guide.pdf)**

benchmarks/_stats.py

@@ -0,0 +1,104 @@
+# Copyright 2026 Firefly Software Foundation.
+"""Publishable-grade comparison statistics for the benchmark harnesses.
+
+The unit of analysis is the DATASET (Demsar, JMLR 2006): callers aggregate fold scores to one value
+per dataset per system BEFORE anything here runs. ``paired_wilcoxon`` compares two systems across
+datasets (one-sided); ``friedman_holm`` compares k systems (Friedman test, mean ranks, Holm-corrected
+two-sided pairwise Wilcoxon post-hocs). Pooling correlated fold deltas is exactly the error this
+module exists to prevent.
+"""
+
+from __future__ import annotations
+
+from collections.abc import Mapping, Sequence
+from dataclasses import dataclass
+from itertools import combinations
+from typing import Any
+
+
+@dataclass(frozen=True)
+class PairedResult:
+    """Paired comparison of two systems across datasets (one aggregated value per dataset)."""
+
+    mean_delta: float
+    wins: int
+    ties: int
+    losses: int
+    p_value: float
+    n: int
+
+
+def paired_wilcoxon(a: Sequence[float], b: Sequence[float], *, band: float = 1e-3) -> PairedResult:
+    """One-sided Wilcoxon signed-rank test that system ``a`` beats system ``b`` across datasets.
+
+    ``a`` and ``b`` hold one aggregated score per dataset, in the same dataset order. Win/tie/loss
+    uses the practical-equivalence ``band``: ``|delta| <= band`` counts as a tie.
+    """
+    deltas = [x - y for x, y in zip(a, b, strict=True)]
+    if not deltas:
+        raise ValueError("paired_wilcoxon needs at least one dataset")
+    wins = sum(d > band for d in deltas)
+    losses = sum(d < -band for d in deltas)
+    ties = len(deltas) - wins - losses
+    if all(abs(d) < 1e-12 for d in deltas):
+        p_value = 1.0  # identical systems: no evidence of a difference (scipy would raise)
+    else:
+        from scipy.stats import wilcoxon
+
+        p_value = float(wilcoxon(deltas, alternative="greater").pvalue)
+    return PairedResult(
+        mean_delta=sum(deltas) / len(deltas), wins=wins, ties=ties, losses=losses, p_value=p_value, n=len(deltas)
+    )
+
+
+def holm(p_values: Sequence[float]) -> list[float]:
+    """Holm step-down correction; returns adjusted p-values in the input order."""
+    order = sorted(range(len(p_values)), key=lambda i: p_values[i])
+    m = len(p_values)
+    adjusted = [0.0] * m
+    running = 0.0
+    for rank, i in enumerate(order):
+        running = max(running, (m - rank) * p_values[i])
+        adjusted[i] = min(1.0, running)
+    return adjusted
+
+
+def friedman_holm(scores: Mapping[str, Sequence[float]]) -> dict[str, Any]:
+    """Compare k systems on the same datasets: Friedman + mean ranks + Holm-corrected pairwise Wilcoxons.
+
+    ``scores`` maps system name -> one aggregated score per dataset (same dataset order everywhere).
+    Returns ``{"friedman_p": float, "mean_ranks": {system: rank}, "pairwise": {(a, b): holm_p}}``
+    with rank 1 = best; ``friedman_p`` is NaN when fewer than three systems are compared.
+    """
+    systems = list(scores)
+    if len(systems) < 2:
+        raise ValueError("friedman_holm needs at least two systems")
+    n_datasets = len(scores[systems[0]])
+    if any(len(scores[s]) != n_datasets for s in systems):
+        raise ValueError("every system needs one score per dataset (equal-length sequences)")
+
+    import numpy as np
+    from scipy.stats import friedmanchisquare, rankdata, wilcoxon
+
+    matrix = np.asarray([list(scores[s]) for s in systems], dtype=float)  # (k systems, n datasets)
+    ranks = np.vstack([rankdata(-matrix[:, j]) for j in range(n_datasets)])  # (n, k); rank 1 = best
+    mean_ranks = {s: float(ranks[:, i].mean()) for i, s in enumerate(systems)}
+
+    friedman_p = float("nan")
+    if len(systems) >= 3:
+        friedman_p = float(friedmanchisquare(*[matrix[i] for i in range(len(systems))]).pvalue)
+
+    pairs = list(combinations(systems, 2))
+    raw: list[float] = []
+    for sys_a, sys_b in pairs:
+        deltas = matrix[systems.index(sys_a)] - matrix[systems.index(sys_b)]
+        raw.append(1.0 if bool(np.all(np.abs(deltas) < 1e-12)) else float(wilcoxon(deltas).pvalue))
+    adjusted = holm(raw)
+    return {
+        "friedman_p": friedman_p,
+        "mean_ranks": mean_ranks,
+        "pairwise": dict(zip(pairs, adjusted, strict=True)),
+    }
+
+
+__all__ = ["PairedResult", "friedman_holm", "holm", "paired_wilcoxon"]

benchmarks/ablation_eval.py

@@ -0,0 +1,152 @@
+# Copyright 2026 Firefly Software Foundation.
+"""Component ablation on real OpenML-CC18 tasks: does each new engine piece pay its way?
+
+Arms (identical data, splits, seeds and trainer portfolio — only the engine knobs differ):
+
+    legacy           untuned defaults, no portfolio, no ensemble, no class-weighting (the old engine)
+    tuned            + Optuna search at the shared trial budget
+    tuned+portfolio  + zero-shot config portfolio seeding
+    tuned+ensemble   + OOF-measured greedy weighted ensemble (no portfolio)
+    full             everything on (the shipped defaults)
+
+Each arm is fit on the same train split and scored on the same held-out split (ROC-AUC). Statistics
+follow Demsar (2006): one score per dataset per arm, Wilcoxon across datasets for full-vs-each-arm,
+and Friedman + Holm-corrected pairwise tests across all arms. Requires the ``data`` extra + network.
+
+Run:  python benchmarks/ablation_eval.py
+"""
+
+from __future__ import annotations
+
+import sys
+import time
+from collections.abc import Mapping, Sequence
+from pathlib import Path
+from typing import Any
+
+from fireflyframework_datascience.automl import AutoML
+from fireflyframework_datascience.datasets import Dataset
+
+sys.path.insert(0, str(Path(__file__).resolve().parent))
+from _stats import friedman_holm, paired_wilcoxon  # noqa: E402  (sibling module on sys.path)
+
+# CPU-feasible OpenML-CC18 binary tasks (id -> short label); mixed types, all small enough for laptops.
+CC18_ABLATION: dict[int, str] = {
+    31: "credit-g",
+    37: "diabetes",
+    1049: "pc4",
+    1464: "blood-transfusion",
+    1480: "ilpd",
+    1494: "qsar-biodeg",
+    1510: "wdbc",
+    40994: "climate-model",
+}
+ARMS: tuple[str, ...] = ("legacy", "tuned", "tuned+portfolio", "tuned+ensemble", "full")
+SEED = 0
+
+
+class _NoPortfolio:
+    """Portfolio stub returning no configs — disables zero-shot seeding for an ablation arm."""
+
+    def configs_for(self, trainer: str, *, n_rows: int, n_features: int) -> list[dict[str, Any]]:
+        return []
+
+
+def build_arm(arm: str, *, cv: int, n_trials: int, trainers: Sequence[Any] | None = None) -> AutoML:
+    """One AutoML instance per arm; every arm shares cv/seed/trainers so only the knobs differ."""
+    from fireflyframework_datascience.search.adapters import DefaultSearchPolicy
+
+    common: dict[str, Any] = {"cv": cv, "random_state": SEED}
+    if trainers is not None:
+        common["trainers"] = list(trainers)
+    if arm == "legacy":
+        return AutoML(
+            search_policy=DefaultSearchPolicy(),
+            n_trials=1,
+            ensemble=False,
+            portfolio=_NoPortfolio(),
+            auto_class_weight=False,
+            **common,
+        )
+    if arm == "tuned":
+        return AutoML(n_trials=n_trials, ensemble=False, portfolio=_NoPortfolio(), auto_class_weight=False, **common)
+    if arm == "tuned+portfolio":
+        return AutoML(n_trials=n_trials, ensemble=False, auto_class_weight=False, **common)
+    if arm == "tuned+ensemble":
+        return AutoML(n_trials=n_trials, ensemble="auto", portfolio=_NoPortfolio(), auto_class_weight=False, **common)
+    if arm == "full":
+        return AutoML(n_trials=n_trials, ensemble="auto", auto_class_weight=True, **common)
+    raise ValueError(f"Unknown ablation arm {arm!r}")
+
+
+def run_ablation(
+    datasets: Sequence[Dataset],
+    *,
+    cv: int = 5,
+    n_trials: int = 20,
+    test_size: float = 0.25,
+    trainers: Sequence[Any] | None = None,
+    arms: Sequence[str] = ARMS,
+) -> dict[str, Any]:
+    """Fit every arm on identical train/test splits; return holdout ROC-AUC + wall-clock per arm/dataset."""
+    scores: dict[str, dict[str, float]] = {arm: {} for arm in arms}
+    seconds: dict[str, dict[str, float]] = {arm: {} for arm in arms}
+    for dataset in datasets:
+        train, test = dataset.train_test_split(test_size=test_size, random_state=SEED)
+        for arm in arms:
+            automl = build_arm(arm, cv=cv, n_trials=n_trials, trainers=trainers)
+            started = time.perf_counter()
+            outcome = automl.fit(train, metric="roc_auc")
+            elapsed = time.perf_counter() - started
+            scores[arm][dataset.name] = float(outcome.evaluate(test).primary_value)
+            seconds[arm][dataset.name] = round(elapsed, 2)
+    return {"scores": scores, "seconds": seconds}
+
+
+def summarize(
+    scores: Mapping[str, Mapping[str, float]], seconds: Mapping[str, Mapping[str, float]] | None = None
+) -> str:
+    """Markdown report: per-dataset table, mean wall-clock, and per-dataset (Demsar-correct) statistics."""
+    arms = list(scores)
+    dataset_names = list(next(iter(scores.values())))
+    lines = ["# AutoML component ablation (holdout ROC-AUC · identical splits, folds and seeds)", ""]
+    lines.append("| dataset | " + " | ".join(arms) + " |")
+    lines.append("|---|" + "---:|" * len(arms))
+    for name in dataset_names:
+        lines.append("| " + " | ".join([name, *[f"{scores[arm][name]:.4f}" for arm in arms]]) + " |")
+    if seconds is not None:
+        mean_secs = [f"{sum(seconds[arm].values()) / len(seconds[arm]):.1f}s" for arm in arms]
+        lines.append("| mean wall-clock | " + " | ".join(mean_secs) + " |")
+    if len(dataset_names) >= 2 and len(arms) >= 2:
+        lines += ["", "## Statistics — unit of analysis is the dataset (one score each)", ""]
+        full = arms[-1]
+        for arm in arms[:-1]:
+            res = paired_wilcoxon([scores[full][d] for d in dataset_names], [scores[arm][d] for d in dataset_names])
+            lines.append(
+                f"- **{full} vs {arm}**: mean Δ={res.mean_delta:+.4f}, "
+                f"W/T/L={res.wins}/{res.ties}/{res.losses}, one-sided Wilcoxon p={res.p_value:.4g} (n={res.n})"
+            )
+        multi = friedman_holm({arm: [scores[arm][d] for d in dataset_names] for arm in arms})
+        ranks = ", ".join(f"{arm}={multi['mean_ranks'][arm]:.2f}" for arm in arms)
+        lines += ["", f"- Friedman p={multi['friedman_p']:.4g} · mean ranks (1=best): {ranks}"]
+        for (first, second), p in multi["pairwise"].items():
+            lines.append(f"  - Holm-adjusted Wilcoxon {first} vs {second}: p={p:.4g}")
+    return "\n".join(lines)
+
+
+def load_suite(dataset_ids: Sequence[int] | None = None) -> list[Dataset]:
+    """Load the CC18 ablation datasets (the facade owns string-target encoding — no hand-encoding here)."""
+    from fireflyframework_datascience.datasets.adapters import OpenMLDatasetLoader
+
+    loader = OpenMLDatasetLoader()
+    return [loader.load(f"openml:{data_id}") for data_id in (dataset_ids or list(CC18_ABLATION))]
+
+
+def main() -> None:
+    print("Firefly DataScience — component ablation on OpenML-CC18 (Tier-1)\n")
+    report = run_ablation(load_suite())
+    print(summarize(report["scores"], report["seconds"]))
+
+
+if __name__ == "__main__":
+    main()

benchmarks/amlb_benchmark.py

@@ -54,7 +54,7 @@ def run_amlb(dataset_ids: list[int] | None = None, *, cv: int = 5) -> list[Bench
         dataset = _encode_classification_target(loader.load(f"openml:{data_id}"))
         train, test = dataset.train_test_split(test_size=0.25, random_state=0)
         started = time.perf_counter()
-        outcome = AutoML(cv=cv).fit(train)
+        outcome = AutoML(cv=cv, n_trials=1).fit(train)  # untuned default arm (tuned arms live in the ablation)
         elapsed = time.perf_counter() - started
         holdout = outcome.evaluate(test)
         results.append(

benchmarks/automl_benchmark.py

@@ -44,7 +44,7 @@ def run_suite(datasets: list[str] | None = None, *, cv: int = 3, test_size: floa
         dataset = loader.load(name)
         train, test = dataset.train_test_split(test_size=test_size, random_state=0)
         started = time.perf_counter()
-        outcome = AutoML(cv=cv).fit(train)
+        outcome = AutoML(cv=cv, n_trials=1).fit(train)  # untuned default arm (tuned arms live in the ablation)
         elapsed = time.perf_counter() - started
         holdout = outcome.evaluate(test)
         results.append(