diff --git a/CHANGELOG.md b/CHANGELOG.md index 376275c..7fd7abb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,6 +10,7 @@ All notable changes to this project will be documented in this file. - Surrogate accuracy reporting: `fit_surrogate()`/`SurrogateModel` (and `fit_regime_surrogate()`/`RegimeSurrogate` via passthrough) now compute a uniform k-fold cross-validated `cv_r2`/`cv_rmse` per observable for both the `gp` and `rf` backends. Warns (`warn_below_r2`, default threshold `0.0`) at fit time and in `RegimeSurrogate.recommend()` when an observable's accuracy is too low to trust (#114). - `sensitivity_from_table()`: post-hoc Sobol/Morris sensitivity from an already-collected `ResultsTable`, by fitting a cheap surrogate over it (`fit_surrogate`) and running `screen()`'s existing machinery against the surrogate instead of a fresh simulator. Unlike a marginal Spearman correlation, this correctly detects non-monotonic (e.g. U-shaped) factor effects. Returns a `TableSensitivity` with `importance` indices and the surrogate's `surrogate_cv_r2` (#114) so callers can judge whether to trust the result (#113). - `sobol_indices()`: like `screen(method="sobol")`, but returns both S1 (first-order) and ST (total-order) per observable instead of discarding ST. `screen()` itself is unchanged for backward compatibility; `ST - S1` is the standard way to detect interaction effects that a first-order-only view misses entirely (#120). +- Replicate averaging in `run_adaptive()` and `screen()`/`sobol_indices()` (`n_reps`, #122): the same `run_grid(..., n_reps=N)` convention (#112), now applied consistently across every entry point that repeatedly evaluates a simulator/`run_fn`. `run_adaptive` detects an opt-in `rep` keyword on `Simulator.generate` and averages each trial's objective(s) over `n_reps` draws before Optuna sees them; `screen()`/`sobol_indices()` detect the same convention on the bare `run_fn` callable they take instead. Without this, adaptive optimization or sensitivity screening against a stochastic simulator could select a "best" config or report an "important" factor that's actually just a single lucky/unlucky data draw. Default `n_reps=1` preserves prior behavior. ## [0.2.0] — 2026-08-17 diff --git a/src/trade_study/design.py b/src/trade_study/design.py index 5db69e4..cfc1111 100644 --- a/src/trade_study/design.py +++ b/src/trade_study/design.py @@ -5,6 +5,7 @@ from __future__ import annotations +import inspect from dataclasses import dataclass, field from enum import Enum from itertools import product @@ -417,12 +418,57 @@ def _rejection_sample( return accepted +def _run_fn_accepts_rep(run_fn: Callable[..., dict[str, float]]) -> bool: + """Whether ``run_fn`` opts into the ``rep`` convention. + + Mirrors :func:`~trade_study.runner._generate_accepts_rep`'s convention + for :meth:`Simulator.generate`, applied to the bare callable + ``screen()``/``sobol_indices()`` take instead: a ``run_fn`` that wants + per-replicate stochasticity under ``n_reps>1`` may accept an optional + keyword-only ``rep`` parameter. Detected via introspection so + ``run_fn``s written before replication support existed keep working + unmodified. + + Returns: + True if ``run_fn``'s signature includes a ``rep`` parameter. + """ + try: + sig = inspect.signature(run_fn) + except (TypeError, ValueError): + return False + return "rep" in sig.parameters + + +def _evaluate_averaged( + run_fn: Callable[..., dict[str, float]], + cfg: dict[str, Any], + n_reps: int, + *, + supports_rep: bool, +) -> dict[str, float]: + """Evaluate ``run_fn`` at ``cfg``, averaging over ``n_reps`` replicates. + + Returns: + Mapping from observable name to its mean over ``n_reps`` calls (a + single call's result, unchanged, when ``n_reps == 1``). + """ + if n_reps == 1: + return run_fn(cfg, rep=0) if supports_rep else run_fn(cfg) + rep_scores = [ + run_fn(cfg, rep=rep) if supports_rep else run_fn(cfg) for rep in range(n_reps) + ] + return { + name: float(np.mean([s[name] for s in rep_scores])) for name in rep_scores[0] + } + + def screen( run_fn: Callable[[dict[str, Any]], dict[str, float]], factors: list[Factor], *, method: str = "morris", n_trajectories: int = 100, + n_reps: int = 1, seed: int = 42, ) -> dict[str, NDArray[np.floating[Any]]]: """Screen factors for influence on observables via SALib. @@ -435,6 +481,14 @@ def screen( n_trajectories: Number of Morris trajectories. For Sobol, this controls the base sample size *N*; the total number of model evaluations is *N* x (num_vars + 2). + n_reps: Replicate draws averaged into each sampled point's scores + before SALib sees them (#122), the same convention + ``run_grid(..., n_reps>1)`` uses (#112) -- ``run_fn`` opts in + via an optional keyword-only ``rep`` parameter, detected via + introspection. Without it, a single stochastic draw per point + can produce importance estimates that are artifacts of that + one draw rather than the factor's real effect. Default 1 (a + single draw per point, today's behavior). seed: Random seed. Returns: @@ -442,13 +496,16 @@ def screen( (mu_star for Morris, S1 for Sobol), one value per factor. Raises: - ValueError: If *method* is unknown or no continuous factors are - provided. + ValueError: If *method* is unknown, no continuous factors are + provided, or ``n_reps`` is less than 1. """ continuous = [f for f in factors if f.factor_type == FactorType.CONTINUOUS] if not continuous: msg = "Screening requires at least one continuous factor" raise ValueError(msg) + if n_reps < 1: + msg = f"n_reps must be >= 1; got {n_reps}" + raise ValueError(msg) problem: dict[str, Any] = { "num_vars": len(continuous), @@ -457,9 +514,9 @@ def screen( } if method == "morris": - return _screen_morris(run_fn, problem, n_trajectories, seed) + return _screen_morris(run_fn, problem, n_trajectories, n_reps, seed) if method == "sobol": - return _screen_sobol(run_fn, problem, n_trajectories, seed) + return _screen_sobol(run_fn, problem, n_trajectories, n_reps, seed) msg = f"Unknown screening method: {method!r}" raise ValueError(msg) @@ -469,6 +526,7 @@ def _screen_morris( run_fn: Callable[[dict[str, Any]], dict[str, float]], problem: dict[str, Any], n_trajectories: int, + n_reps: int, seed: int, ) -> dict[str, NDArray[np.floating[Any]]]: """Morris elementary-effects screening. @@ -480,11 +538,12 @@ def _screen_morris( from SALib.sample import morris as morris_sample # type: ignore[import-untyped] param_values = morris_sample.sample(problem, n_trajectories, seed=seed) + supports_rep = _run_fn_accepts_rep(run_fn) results_by_obs: dict[str, list[float]] = {} for row in param_values: cfg = dict(zip(problem["names"], row, strict=True)) - scores = run_fn(cfg) + scores = _evaluate_averaged(run_fn, cfg, n_reps, supports_rep=supports_rep) for obs_name, val in scores.items(): results_by_obs.setdefault(obs_name, []).append(val) @@ -505,6 +564,7 @@ def _sobol_sample_and_evaluate( run_fn: Callable[[dict[str, Any]], dict[str, float]], problem: dict[str, Any], n_samples: int, + n_reps: int, seed: int, ) -> dict[str, list[float]]: """Draw a Saltelli design and evaluate ``run_fn`` at every point. @@ -519,11 +579,12 @@ def _sobol_sample_and_evaluate( from SALib.sample import sobol as sobol_sample param_values = sobol_sample.sample(problem, n_samples, seed=seed) + supports_rep = _run_fn_accepts_rep(run_fn) results_by_obs: dict[str, list[float]] = {} for row in param_values: cfg = dict(zip(problem["names"], row, strict=True)) - scores = run_fn(cfg) + scores = _evaluate_averaged(run_fn, cfg, n_reps, supports_rep=supports_rep) for obs_name, val in scores.items(): results_by_obs.setdefault(obs_name, []).append(val) return results_by_obs @@ -533,6 +594,7 @@ def _screen_sobol( run_fn: Callable[[dict[str, Any]], dict[str, float]], problem: dict[str, Any], n_samples: int, + n_reps: int, seed: int, ) -> dict[str, NDArray[np.floating[Any]]]: """Sobol variance-based sensitivity analysis. @@ -542,7 +604,9 @@ def _screen_sobol( """ from SALib.analyze import sobol as sobol_analyze - results_by_obs = _sobol_sample_and_evaluate(run_fn, problem, n_samples, seed) + results_by_obs = _sobol_sample_and_evaluate( + run_fn, problem, n_samples, n_reps, seed + ) importance: dict[str, NDArray[np.floating[Any]]] = {} for obs_name, vals in results_by_obs.items(): @@ -561,6 +625,7 @@ def sobol_indices( factors: list[Factor], *, n_samples: int = 100, + n_reps: int = 1, seed: int = 42, ) -> dict[str, tuple[NDArray[np.floating[Any]], NDArray[np.floating[Any]]]]: """Sobol first- and total-order sensitivity indices (#120). @@ -579,6 +644,10 @@ def sobol_indices( varied (as in ``screen()``). n_samples: Base sample size *N* for the Saltelli design; total evaluations are *N* x (2 x num_vars + 2). + n_reps: Replicate draws averaged into each sampled point's scores + before SALib sees them (#122) -- see ``screen()``'s ``n_reps`` + for the full rationale and the ``run_fn`` opt-in convention. + Default 1 (a single draw per point, today's behavior). seed: Random seed. Returns: @@ -587,7 +656,8 @@ def sobol_indices( appear in ``factors``). Raises: - ValueError: If no continuous factors are provided. + ValueError: If no continuous factors are provided, or ``n_reps`` + is less than 1. """ from SALib.analyze import sobol as sobol_analyze @@ -595,13 +665,18 @@ def sobol_indices( if not continuous: msg = "Screening requires at least one continuous factor" raise ValueError(msg) + if n_reps < 1: + msg = f"n_reps must be >= 1; got {n_reps}" + raise ValueError(msg) problem: dict[str, Any] = { "num_vars": len(continuous), "names": [f.name for f in continuous], "bounds": [list(f.bounds) for f in continuous if f.bounds is not None], } - results_by_obs = _sobol_sample_and_evaluate(run_fn, problem, n_samples, seed) + results_by_obs = _sobol_sample_and_evaluate( + run_fn, problem, n_samples, n_reps, seed + ) indices: dict[str, tuple[NDArray[np.floating[Any]], NDArray[np.floating[Any]]]] = {} for obs_name, vals in results_by_obs.items(): diff --git a/src/trade_study/runner.py b/src/trade_study/runner.py index e1148df..5503030 100644 --- a/src/trade_study/runner.py +++ b/src/trade_study/runner.py @@ -183,6 +183,7 @@ def run_adaptive( observables: list[Observable], *, n_trials: int = 100, + n_reps: int = 1, seed: int = 42, ) -> ResultsTable: """Run adaptive multi-objective optimization via optuna. @@ -193,15 +194,35 @@ def run_adaptive( factors: Factor definitions (from design module). observables: Observable definitions. n_trials: Number of optuna trials. + n_reps: Replicate draws averaged into each trial's objective + values before Optuna sees them (#122), the same convention + ``run_grid(..., n_reps>1)`` uses (#112) -- a simulator opts in + via an optional keyword-only ``rep`` parameter on + ``generate``, detected via introspection. Without it, NSGA-II + can select a "best" config that's just a lucky draw for one + data realization rather than one that's robust across draws. + Simulators that don't accept ``rep`` see identical draws + repeated ``n_reps`` times, which only adds redundant compute. + Default 1 (a single draw per trial, today's behavior). seed: Random seed. Returns: - ResultsTable with scored results. + ResultsTable with scored results, one row per optuna trial (each + row the mean of ``n_reps`` replicate draws). + + Raises: + ValueError: If ``n_reps`` is less than 1. """ import optuna as _optuna from .design import FactorType + if n_reps < 1: + msg = f"n_reps must be >= 1; got {n_reps}" + raise ValueError(msg) + + supports_rep = _generate_accepts_rep(world) + directions_str = [ "minimize" if o.direction == Direction.MINIMIZE else "maximize" for o in observables @@ -229,10 +250,15 @@ def objective(trial: optuna.trial.Trial) -> tuple[float, ...]: FactorType.DISCRETE, }: config[f.name] = trial.suggest_categorical(f.name, f.levels) - truth, observations = world.generate(config) - scores = scorer.score(truth, observations, config) + rep_scores: list[dict[str, float]] = [] + for rep in range(n_reps): + if supports_rep: + truth, observations = world.generate(config, rep=rep) # type: ignore[call-arg] + else: + truth, observations = world.generate(config) + rep_scores.append(scorer.score(truth, observations, config)) return tuple( - scores.get(name, float("nan")) * w + float(np.mean([s.get(name, float("nan")) for s in rep_scores])) * w for name, w in zip(obs_names, obs_weights, strict=True) ) diff --git a/tests/test_design.py b/tests/test_design.py index 58d7c17..4201958 100644 --- a/tests/test_design.py +++ b/tests/test_design.py @@ -2,6 +2,7 @@ from __future__ import annotations +import importlib from typing import Any import numpy as np @@ -17,6 +18,8 @@ sobol_indices, ) +_design = importlib.import_module("trade_study.design") + # --------------------------------------------------------------------------- # Fixtures # --------------------------------------------------------------------------- @@ -334,6 +337,91 @@ def test_screen_rejects_no_continuous() -> None: screen(lambda _c: {"y": 0.0}, factors) +# --------------------------------------------------------------------------- +# screen/sobol_indices — replicate averaging (#122) +# --------------------------------------------------------------------------- + + +def test_run_fn_accepts_rep_detects_rep_param() -> None: + def with_rep(_c: dict[str, Any], *, rep: int = 0) -> dict[str, float]: + return {"y": float(rep)} + + def without_rep(_c: dict[str, Any]) -> dict[str, float]: + return {"y": 0.0} + + assert _design._run_fn_accepts_rep(with_rep) is True # ruff: ignore[private-member-access] + assert _design._run_fn_accepts_rep(without_rep) is False # ruff: ignore[private-member-access] + + +def test_evaluate_averaged_single_call_when_n_reps_one() -> None: + calls: list[int] = [] + + def run_fn(_c: dict[str, Any]) -> dict[str, float]: + calls.append(1) + return {"y": 5.0} + + result = _design._evaluate_averaged(run_fn, {}, 1, supports_rep=False) # ruff: ignore[private-member-access] + assert result == {"y": 5.0} + assert len(calls) == 1 + + +def test_evaluate_averaged_averages_across_reps() -> None: + def run_fn(_c: dict[str, Any], *, rep: int = 0) -> dict[str, float]: + return {"y": float(rep)} + + # average of reps 0..3 is 1.5 + result = _design._evaluate_averaged(run_fn, {}, 4, supports_rep=True) # ruff: ignore[private-member-access] + assert result == pytest.approx({"y": 1.5}) + + +def test_evaluate_averaged_repeats_identical_draw_when_run_fn_ignores_rep() -> None: + calls: list[int] = [] + + def run_fn(_c: dict[str, Any]) -> dict[str, float]: + calls.append(1) + return {"y": 7.0} + + result = _design._evaluate_averaged(run_fn, {}, 3, supports_rep=False) # ruff: ignore[private-member-access] + assert result == pytest.approx({"y": 7.0}) + assert len(calls) == 3 + + +def test_screen_rejects_invalid_n_reps(continuous_factors: list[Factor]) -> None: + with pytest.raises(ValueError, match="n_reps must be >= 1"): + screen(_linear_model, continuous_factors, n_trajectories=5, n_reps=0) + + +def test_screen_n_reps_passes_incrementing_rep_to_run_fn( + continuous_factors: list[Factor], +) -> None: + seen_reps: set[int] = set() + + def run_fn(cfg: dict[str, Any], *, rep: int = 0) -> dict[str, float]: + seen_reps.add(rep) + return {"y": cfg["alpha"]} + + screen(run_fn, continuous_factors, n_trajectories=5, n_reps=3, seed=0) + assert seen_reps == {0, 1, 2} + + +def test_sobol_indices_rejects_invalid_n_reps(continuous_factors: list[Factor]) -> None: + with pytest.raises(ValueError, match="n_reps must be >= 1"): + sobol_indices(_linear_model, continuous_factors, n_samples=8, n_reps=0) + + +def test_sobol_indices_n_reps_passes_incrementing_rep_to_run_fn( + continuous_factors: list[Factor], +) -> None: + seen_reps: set[int] = set() + + def run_fn(cfg: dict[str, Any], *, rep: int = 0) -> dict[str, float]: + seen_reps.add(rep) + return {"y": cfg["alpha"]} + + sobol_indices(run_fn, continuous_factors, n_samples=8, n_reps=3, seed=0) + assert seen_reps == {0, 1, 2} + + # --------------------------------------------------------------------------- # screen — Sobol (#76) # --------------------------------------------------------------------------- diff --git a/tests/test_runner.py b/tests/test_runner.py index 7dd2a80..5297b8e 100644 --- a/tests/test_runner.py +++ b/tests/test_runner.py @@ -333,6 +333,71 @@ def test_run_adaptive_deterministic_seed( np.testing.assert_allclose(r1.scores, r2.scores) +# --------------------------------------------------------------------------- +# run_adaptive replicate averaging (#122) +# --------------------------------------------------------------------------- + + +def test_run_adaptive_rejects_non_positive_n_reps( + world: _ToySimulator, + scorer: _ToyScorer, + observables: list[Observable], +) -> None: + factors = [Factor("alpha", FactorType.CONTINUOUS, bounds=(0.0, 1.0))] + with pytest.raises(ValueError, match="n_reps must be >= 1"): + run_adaptive(world, scorer, factors, observables, n_trials=5, n_reps=0) + + +def test_run_adaptive_n_reps_default_matches_n_reps_one( + world: _ToySimulator, + scorer: _ToyScorer, + observables: list[Observable], +) -> None: + """n_reps defaults to 1: scores unchanged from before #122.""" + factors = [Factor("alpha", FactorType.CONTINUOUS, bounds=(0.0, 1.0))] + default = run_adaptive(world, scorer, factors, observables, n_trials=10, seed=3) + explicit = run_adaptive( + world, scorer, factors, observables, n_trials=10, n_reps=1, seed=3 + ) + np.testing.assert_allclose(default.scores, explicit.scores) + + +def test_run_adaptive_non_rep_aware_simulator_repeats_identically( + world: _ToySimulator, + scorer: _ToyScorer, + observables: list[Observable], +) -> None: + """A simulator without a rep parameter sees the same draw n_reps times.""" + factors = [Factor("alpha", FactorType.CONTINUOUS, bounds=(0.0, 1.0))] + single = run_adaptive(world, scorer, factors, observables, n_trials=10, seed=5) + replicated = run_adaptive( + world, scorer, factors, observables, n_trials=10, n_reps=4, seed=5 + ) + np.testing.assert_allclose(single.scores, replicated.scores) + + +def test_run_adaptive_rep_aware_simulator_averages_across_reps( + rep_world: _RepAwareSimulator, + rep_scorer: _RepSensitiveScorer, + observables: list[Observable], +) -> None: + """A rep-aware simulator's objective is the mean over n_reps draws. + + _RepAwareSimulator offsets alpha by ``0.01 * rep``; at n_reps=3 the + averaged cost should reflect the mean offset (0.01) rather than a + single draw's (0 for rep=0 alone). + """ + factors = [Factor("alpha", FactorType.CONTINUOUS, bounds=(0.0, 1.0))] + result = run_adaptive( + rep_world, rep_scorer, factors, observables, n_trials=1, n_reps=3, seed=0 + ) + cost_idx = result.observable_names.index("cost") + suggested_alpha = result.configs[0]["alpha"] + # mean offset over rep=0,1,2 is 0.01*(0+1+2)/3. + expected_alpha = suggested_alpha + 0.01 * (0 + 1 + 2) / 3 + assert result.scores[0, cost_idx] == pytest.approx(expected_alpha * 10.0) + + # --------------------------------------------------------------------------- # Progress callback (#77) # ---------------------------------------------------------------------------