From 3cdce5775f8a2dc0e74280bb38513e1146e34387 Mon Sep 17 00:00:00 2001 From: Javier Garcia Ordonez Date: Tue, 29 Sep 2026 16:15:38 +0200 Subject: [PATCH 1/4] feat(data+evaluate): port incubator CV diagnostics, coverage, convergence fits, scale selection (T48np) Clean re-port of the still-unpromoted NIH-incubator tooling (CONFIDENTIAL/nih-merck-validation) onto the log-scale/api era: - data: select_variable_scale, auto_select_distributions - evaluate: compute_coverage, compute_cv_diagnostics, fit_convergence_exponential(_asymptotic), bootstrap+draws+diagnostics redesign of compute_cv_convergence, mean_signed_error in CV metrics - api CVAccuracyMetrics absorbs mean_signed_error - propagate_uq_with_uncertainty deliberately NOT ported: the api-era propagate_manual_uq_with_uncertainty already covers it Tests: NIH versions of test_funs_evaluate_cv_stats.py (superset) and test_dakota_funs_data_processing.py additions, ruff-normalized. (cherry picked from commit 35d2b9933aa26bffdbcc65ed1878eeac9a5313a3) --- src/itis_sumo/api/types.py | 1 + src/itis_sumo/data/__init__.py | 4 + src/itis_sumo/data/funs_data_processing.py | 166 ++++++++- src/itis_sumo/evaluate/__init__.py | 8 + src/itis_sumo/evaluate/funs_evaluate.py | 410 +++++++++++++++++++-- tests/test_dakota_funs_data_processing.py | 91 +++++ tests/test_funs_evaluate_cv_stats.py | 370 +++++++++++++++++++ 7 files changed, 1018 insertions(+), 32 deletions(-) diff --git a/src/itis_sumo/api/types.py b/src/itis_sumo/api/types.py index 6eeec11..2d38b8c 100644 --- a/src/itis_sumo/api/types.py +++ b/src/itis_sumo/api/types.py @@ -228,6 +228,7 @@ class CVAccuracyMetrics: sum_abs: float mean_abs: float max_abs: float + mean_signed_error: float seed: int diff --git a/src/itis_sumo/data/__init__.py b/src/itis_sumo/data/__init__.py index 1c8482c..275140f 100644 --- a/src/itis_sumo/data/__init__.py +++ b/src/itis_sumo/data/__init__.py @@ -1,6 +1,7 @@ """Training-data parsing, filtering, sampling grids, and statistics.""" from itis_sumo.data.funs_data_processing import ( + auto_select_distributions, compute_correlation_indices, create_grid_samples, create_manual_uq_samples, @@ -16,6 +17,7 @@ load_data, process_input_file, sanitize_varnames, + select_variable_scale, ) from itis_sumo.data.funs_dataset_diagnostics import ( DatasetDiagnostics, @@ -29,6 +31,7 @@ "OutlierSummary", "VariableDiagnostics", "analyze_dataset", + "auto_select_distributions", "compute_correlation_indices", "create_grid_samples", "create_manual_uq_samples", @@ -44,4 +47,5 @@ "load_data", "process_input_file", "sanitize_varnames", + "select_variable_scale", ] diff --git a/src/itis_sumo/data/funs_data_processing.py b/src/itis_sumo/data/funs_data_processing.py index 891c58c..d45a991 100644 --- a/src/itis_sumo/data/funs_data_processing.py +++ b/src/itis_sumo/data/funs_data_processing.py @@ -5,7 +5,7 @@ import re from collections.abc import Callable, Mapping, Sequence from pathlib import Path -from typing import Literal, TypeVar, overload +from typing import Any, Literal, TypeVar, overload import numpy as np import pandas as pd @@ -582,11 +582,175 @@ def create_manual_uq_samples( elif dist_type == "constant": value = dist_info["value"] samples[var] = [float(value)] * num_samples + # elif dist_type == "lognormal": + # mean = dist_info["mean"] + # std = dist_info["std"] + # samples[var] = float(rng.lognormal(mean, std, (num_samples,))) + # DEV NOTE (internal only, not user-facing): the supported taxonomy is + # deliberately just constant/uniform/normal (each possibly fit in + # log-space via a "log_"-prefixed column, see process_input_file's + # make_log and select_variable_scale/auto_select_distributions below) + # -- no variable we've seen in practice needed anything else. + # Exponential is a plausible future addition (e.g. `{"distribution": + # "exponential", "rate": ...}` -> rng.exponential) but don't surface + # it in docs or auto-selection until there's a concrete case for it. else: raise ValueError(f"Unsupported distribution type: {dist_type}") return samples +def select_variable_scale( + values: np.ndarray | list[float], alpha: float = 0.05 +) -> dict[str, Any]: + """Decide whether a variable is better modeled as normal or uniform, in + linear or log space, by testing all scale/distribution combinations that + apply and picking the one least contradicted by the data. + + Every variable this module samples is assumed to be one of exactly three + shapes -- constant, uniform, or normal (see `create_manual_uq_samples`) + -- each possibly fit in linear or log space (log space achieved by + pre-transforming the data, e.g. `process_input_file(..., make_log=True)`; + a normal fit in log space is the same thing as a lognormal fit in linear + space). This picks among those four non-trivial combinations + automatically instead of assuming one. + + Decision rule: run `scipy.stats.shapiro` (normal fit) and + `scipy.stats.kstest` against a fitted uniform (uniform fit), on both the + raw values and their log (log skipped if any value is <= 0), then pick + whichever of the (up to) four candidates has the highest p-value -- the + hypothesis least rejected by the data. `p_value` is the resulting + quantification of fit quality; `confident` is `p_value > alpha`, i.e. + whether the *winning* candidate is itself a statistically defensible fit + (not rejected at the `alpha` significance level) rather than merely the + least-bad of four options that are all a poor match for the data -- + signaling a call worth a human look rather than a clean choice. + + A variable with (near-)zero variance short-circuits to + `{"distribution": "constant", ...}` without running either test -- + Shapiro-Wilk and the KS test are both undefined/degenerate on constant + data, and constant is the correct answer regardless. + + Args: + values: Raw (linear-space) sample values for one variable. + alpha: Significance level the winning candidate's `p_value` must + clear for `confident` to be `True`. + + Returns: + `{"scale": "linear"|"log", "distribution": "constant"|"uniform"|"normal", + "p_value": float, "confident": bool, "candidates": {...}}` -- + `candidates` keyed `"{scale}_{distribution}"` (log entries omitted if + log isn't defined), each `{"test", "statistic", "p_value"}`. Constant + variables return `"value"` instead of `"p_value"`/`"candidates"`. + """ + from scipy.stats import kstest, shapiro + + arr = np.asarray(values, dtype=float) + if np.ptp(arr) == 0: + return { + "scale": "linear", + "distribution": "constant", + "value": float(arr[0]), + "confident": True, + "candidates": {}, + } + + candidates: dict[str, dict[str, Any]] = {} + scales = {"linear": arr} + if np.all(arr > 0): + scales["log"] = np.log(arr) + + for scale_name, x in scales.items(): + w_stat, w_p = shapiro(x) + candidates[f"{scale_name}_normal"] = { + "test": "shapiro", + "statistic": float(w_stat), + "p_value": float(w_p), + } + lo, span = float(x.min()), float(x.max() - x.min()) + ks_stat, ks_p = kstest(x, "uniform", args=(lo, span)) + candidates[f"{scale_name}_uniform"] = { + "test": "kstest", + "statistic": float(ks_stat), + "p_value": float(ks_p), + } + + best_key = max(candidates, key=lambda k: candidates[k]["p_value"]) + best_p_value = candidates[best_key]["p_value"] + scale, distribution = best_key.split("_") + + return { + "scale": scale, + "distribution": distribution, + "p_value": best_p_value, + "confident": best_p_value > alpha, + "candidates": candidates, + } + + +def auto_select_distributions( + df: pd.DataFrame, columns: list[str], alpha: float = 0.05 +) -> tuple[dict[str, dict[str, float | str]], dict[str, Any]]: + """Run `select_variable_scale` over each of `columns` and assemble the + result directly into a `distributions` dict usable by + `create_manual_uq_samples`/`evaluate_sobol_indices`, alongside the raw + per-variable diagnostics. + + A column selected for log space is renamed with a `"log_"` prefix in the + returned `distributions` dict (matching `process_input_file`'s + `make_log` convention) -- callers must fit/evaluate the corresponding + surrogate against log-transformed data for that variable to match. + + This is an auto-selected *default*: the returned `distributions` dict is + a plain dict like any hand-written one, so a caller can inspect + `diagnostics[var]` (e.g. print/log it, or surface it however their + application does) and override any entry before passing `distributions` + on to downstream sampling/evaluation. + + Args: + df: DataFrame with one column per variable, raw (linear-space) values. + columns: Column names to select a distribution for. + alpha: Forwarded to `select_variable_scale`. + + Returns: + `(distributions, diagnostics)` -- `distributions` keyed by the + (possibly `log_`-prefixed) variable name, `diagnostics` keyed by the + original column name and holding each variable's full + `select_variable_scale` result (including `"confident"`). + """ + distributions: dict[str, dict[str, float | str]] = {} + diagnostics: dict[str, Any] = {} + + for col in columns: + values = df[col].to_numpy(dtype=float) + decision = select_variable_scale(values, alpha=alpha) + diagnostics[col] = decision + + if decision["distribution"] == "constant": + distributions[col] = { + "distribution": "constant", + "value": decision["value"], + } + continue + + scale = decision["scale"] + x = np.log(values) if scale == "log" else values + key = f"log_{col}" if scale == "log" else col + if decision["distribution"] == "normal": + distributions[key] = { + "distribution": "normal", + "mean": float(np.mean(x)), + "std": float(np.std(x, ddof=1)), + } + else: + distributions[key] = { + "distribution": "uniform", + "min": float(np.min(x)), + "max": float(np.max(x)), + } + + return distributions, diagnostics + + T = TypeVar("T") diff --git a/src/itis_sumo/evaluate/__init__.py b/src/itis_sumo/evaluate/__init__.py index 856d59a..8088679 100644 --- a/src/itis_sumo/evaluate/__init__.py +++ b/src/itis_sumo/evaluate/__init__.py @@ -1,8 +1,10 @@ """End-to-end evaluations on top of Dakota runs.""" from itis_sumo.evaluate.funs_evaluate import ( + compute_coverage, compute_cv_accuracy_metrics, compute_cv_convergence, + compute_cv_diagnostics, compute_paired_ttest, evaluate_sobol_indices, evaluate_sumo, @@ -11,6 +13,8 @@ evaluate_sumo_manual_crossvalidation, evaluate_sumo_on_grid, export_sumo_model, + fit_convergence_exponential, + fit_convergence_exponential_asymptotic, import_sumo_model, perform_moga_optimization, propagate_uq, @@ -18,8 +22,10 @@ ) __all__ = [ + "compute_coverage", "compute_cv_accuracy_metrics", "compute_cv_convergence", + "compute_cv_diagnostics", "compute_paired_ttest", "evaluate_sobol_indices", "evaluate_sumo", @@ -28,6 +34,8 @@ "evaluate_sumo_manual_crossvalidation", "evaluate_sumo_on_grid", "export_sumo_model", + "fit_convergence_exponential", + "fit_convergence_exponential_asymptotic", "import_sumo_model", "perform_moga_optimization", "propagate_uq", diff --git a/src/itis_sumo/evaluate/funs_evaluate.py b/src/itis_sumo/evaluate/funs_evaluate.py index 86a3842..b1dfbea 100644 --- a/src/itis_sumo/evaluate/funs_evaluate.py +++ b/src/itis_sumo/evaluate/funs_evaluate.py @@ -8,6 +8,8 @@ import numpy as np import pandas as pd +from scipy.optimize import curve_fit +from scipy.special import erfinv from scipy.stats import ttest_rel from sklearn.model_selection import KFold @@ -36,6 +38,7 @@ sanitize_varnames, scale_distribution, ) +from itis_sumo.preprocess.data_preprocessor import DataPreprocessor _logger = logging.getLogger(__name__) @@ -175,7 +178,7 @@ def propagate_manual_uq_with_uncertainty( input_vars: list[str], output_response: str, distributions: dict[str, dict[str, float | str]], - preprocessor, + preprocessor: DataPreprocessor | None, num_samples: int, n_histograms: int = 100, seed: int = 42, @@ -193,6 +196,13 @@ def propagate_manual_uq_with_uncertainty( Flask route rather than ``propagate_uq`` (Dakota-native, normal-only, no predictive-uncertainty injection) -- see test_metamodeling_analytical.py. + If ``preprocessor`` is given, samples are drawn/transformed/evaluated in + the preprocessor's mapped space and results are mapped back via + ``inverse_transform``; if omitted, ``input_vars``/``output_response`` are + used directly (sanitized for Dakota) and the returned samples live in the + file's own space -- the pathway the NIH worked example uses to propagate + in a hand-built log-transformed space without a ``DataPreprocessor``. + Returns: A ``(n_histograms, num_samples)`` array of propagated output samples, already inverse-transformed to the caller's original units. @@ -212,14 +222,21 @@ def propagate_manual_uq_with_uncertainty( SAMPLES_FILE = run_dir / "manual_uq_samples.csv" df_samples.to_csv(SAMPLES_FILE, index=False) - df_samples_transformed = preprocessor.transform(df_samples) + if preprocessor is not None: + df_samples_transformed = preprocessor.transform(df_samples) + mapped_input_vars = [ + preprocessor.input_variables[var].mapped_name for var in input_vars + ] + mapped_response = preprocessor.output_variables[output_response].mapped_name + else: + # No preprocessor: the caller's units ARE the training file's space, + # so the samples pass through untransformed (NIH-worked-example path). + df_samples_transformed = df_samples + mapped_input_vars = sanitize_varnames(input_vars) + mapped_response = sanitize_varnames(output_response) PROCESSED_SAMPLES_FILE = run_dir / "manual_uq_samples_processed.csv" df_samples_transformed.to_csv(PROCESSED_SAMPLES_FILE, sep=" ", index=False) - mapped_input_vars = [ - preprocessor.input_variables[var].mapped_name for var in input_vars - ] - mapped_response = preprocessor.output_variables[output_response].mapped_name results = evaluate_sumo( run_dir, PROCESSED_TRAINING_FILE, @@ -245,6 +262,8 @@ def propagate_manual_uq_with_uncertainty( r = np.sqrt(2) * erfinv(rng.uniform(-1 + 1e-10, 1 - 1e-10, size=num_samples)) all_results_transformed[i, :] = prediction + r * uncertainty + if preprocessor is None: + return all_results_transformed all_samples_dict = {mapped_response: all_results_transformed.flatten().tolist()} all_samples_original = preprocessor.inverse_transform(all_samples_dict) return np.asarray(all_samples_original[output_response]).reshape( @@ -592,12 +611,16 @@ def evaluate_sumo_manual_crossvalidation( def compute_cv_accuracy_metrics( actual: list[float] | np.ndarray, predicted: list[float] | np.ndarray ) -> dict[str, float]: - """Compute RMSE/MAE/sum-abs/max-abs directly from paired CV actual/predicted values. + """Compute RMSE/MAE/sum-abs/max-abs/bias directly from paired CV actual/predicted values. Unlike `_parse_crossvalidation_outputlogs`, this does not depend on parsing Dakota's stdout (which `evaluate_sumo_crossvalidation` no longer captures) - it derives the same metrics straight from the actual/predicted arrays already produced by `evaluate_sumo_manual_crossvalidation`. + + `mean_signed_error` (bias) is `mean(actual - predicted)`, signed rather than + absolute - a systematic over/under-prediction can net out near zero in `mean_abs` + if it isn't consistent in sign, so this is a distinct diagnostic, not a duplicate. """ actual_arr = np.asarray(actual, dtype=float) predicted_arr = np.asarray(predicted, dtype=float) @@ -622,6 +645,7 @@ def compute_cv_accuracy_metrics( "sum_abs": float("nan"), "mean_abs": float("nan"), "max_abs": float("nan"), + "mean_signed_error": float("nan"), } residuals = actual_arr - predicted_arr abs_residuals = np.abs(residuals) @@ -630,6 +654,7 @@ def compute_cv_accuracy_metrics( "sum_abs": float(np.sum(abs_residuals)), "mean_abs": float(np.mean(abs_residuals)), "max_abs": float(np.max(abs_residuals)), + "mean_signed_error": float(np.mean(residuals)), } @@ -676,6 +701,201 @@ def _convergence_subset_sizes( return unique_sizes +def _tukey_outlier_mask(residuals: np.ndarray) -> np.ndarray: + """Boolean mask, True where `residuals` fall outside the Tukey IQR fence + (Q1 - 1.5*IQR, Q3 + 1.5*IQR). Flags, does not drop -- see `compute_cv_diagnostics` + and V17kb: unfiltered MAE stays the primary metric, this only feeds + `mean_abs_filtered`/`n_outliers` (a secondary diagnostic). + """ + q1, q3 = np.percentile(residuals, [25, 75]) + iqr = q3 - q1 + lower, upper = q1 - 1.5 * iqr, q3 + 1.5 * iqr + return (residuals < lower) | (residuals > upper) + + +def _paired_cohens_d(actual: np.ndarray, predicted: np.ndarray) -> float: + """Standardized paired-difference effect size (Cohen's dz) for actual vs. predicted. + + dz = mean(actual - predicted) / std(actual - predicted, ddof=1). Unlike the + paired t-test's p-value, dz is not confounded by sample size -- it is the + convergence-stabilization signal V17kb designates as primary, with p-value + evolution kept as a secondary diagnostic (power grows with N regardless of + whether the underlying bias is shrinking). + """ + diff = actual - predicted + diff_std = np.std(diff, ddof=1) + return float(np.mean(diff) / diff_std) if diff_std > 0 else float("nan") + + +def compute_coverage( + actual: list[float] | np.ndarray, + predicted: list[float] | np.ndarray, + predicted_std: list[float] | np.ndarray, + levels: tuple[float, ...] = (0.6827, 0.95, 0.9973), +) -> dict[str, Any]: + """Empirical vs. nominal prediction-interval coverage (V19cz): for each + `levels` entry (default 1sigma/2sigma/3sigma two-sided Gaussian mass), + what fraction of `actual` fall within `predicted +- z*predicted_std`, + z = sqrt(2)*erfinv(level) (same two-sided-z convention as `create_manual_uq_samples`). + Tests whether the surrogate's OWN reported uncertainty + (`{output}_std_hat`, V8df) is trustworthy, not just whether point predictions + are close -- a persistent (non-shrinking-with-N) gap between empirical and + nominal coverage flags a miscalibrated variance model, not a data-shortage. + + Points with NaN/non-positive `predicted_std` are dropped (Dakota can omit + variances.dat for a fold, B23) -- coverage is reported over however many + points have a usable std, `n_points` says how many that was. + """ + actual_arr = np.asarray(actual, dtype=float) + predicted_arr = np.asarray(predicted, dtype=float) + std_arr = np.asarray(predicted_std, dtype=float) + if not (actual_arr.shape == predicted_arr.shape == std_arr.shape): + raise ValueError( + f"actual (shape {actual_arr.shape}), predicted (shape {predicted_arr.shape}), " + f"predicted_std (shape {std_arr.shape}) must have the same shape" + ) + valid = ( + ~np.isnan(actual_arr) + & ~np.isnan(predicted_arr) + & ~np.isnan(std_arr) + & (std_arr > 0) + ) + actual_arr = actual_arr[valid] + predicted_arr = predicted_arr[valid] + std_arr = std_arr[valid] + + n_points = int(actual_arr.size) + if n_points == 0: + return { + "levels": list(levels), + "empirical": [float("nan")] * len(levels), + "n_points": 0, + } + + z_scores = (actual_arr - predicted_arr) / std_arr + empirical = [] + for level in levels: + z_crit = np.sqrt(2) * erfinv(level) + empirical.append(float(np.mean(np.abs(z_scores) <= z_crit))) + return {"levels": list(levels), "empirical": empirical, "n_points": n_points} + + +def compute_cv_diagnostics( + actual: list[float] | np.ndarray, + predicted: list[float] | np.ndarray, + predicted_std: list[float] | np.ndarray | None = None, +) -> dict[str, Any]: + """One-pass statistical bundle for a CV subset (V16wq): everything + `_cv_subset_diagnostics` needs from a single Dakota CV rerun -- accuracy metrics, + Tukey-flagged outlier count + filtered MAE, paired t-test, and Cohen's dz effect + size -- so a caller never has to rerun CV just to derive a different metric. + Also returns the NaN-filtered `actual`/`predicted` pairs themselves so downstream + N=50 diagnostic plots (QQ, diagonal, error-vs-mean) can reuse them without a + second CV rerun. + + `predicted_std` (V19cz), if given, is the per-point GP predictive std + (`{output}_std_hat`) aligned with `actual`/`predicted` before NaN-filtering; + it is filtered by the same actual/predicted validity mask and returned + verbatim (it may still carry its own NaNs from folds that never wrote + variances.dat -- `compute_coverage` drops those independently) so a caller + can pool it across bootstrap draws for coverage without a second CV rerun. + """ + actual_arr = np.asarray(actual, dtype=float) + predicted_arr = np.asarray(predicted, dtype=float) + if actual_arr.shape != predicted_arr.shape: + raise ValueError( + f"actual (shape {actual_arr.shape}) and predicted (shape {predicted_arr.shape}) " + "must have the same shape" + ) + std_arr = ( + np.asarray(predicted_std, dtype=float) + if predicted_std is not None + else np.full(actual_arr.shape, np.nan) + ) + if std_arr.shape != actual_arr.shape: + raise ValueError( + f"predicted_std (shape {std_arr.shape}) must match actual (shape {actual_arr.shape})" + ) + # See compute_cv_accuracy_metrics: NaN entries are dropped CV rows (B22), not real data. + valid = ~np.isnan(actual_arr) & ~np.isnan(predicted_arr) + actual_arr = actual_arr[valid] + predicted_arr = predicted_arr[valid] + std_arr = std_arr[valid] + + metrics = compute_cv_accuracy_metrics(actual_arr, predicted_arr) + base = { + **metrics, + "mean_abs_filtered": float("nan"), + "n_outliers": 0, + "ttest_statistic": float("nan"), + "ttest_p_value": float("nan"), + "cohens_d": float("nan"), + "actual": actual_arr.tolist(), + "predicted": predicted_arr.tolist(), + "predicted_std": std_arr.tolist(), + } + # paired t-test/effect-size/outlier-fence all need >=2 points; below that, + # leave them NaN rather than raise (compute_paired_ttest requires >=2) -- + # matches compute_cv_accuracy_metrics's own NaN-over-crash convention. + if actual_arr.size < 2: + return base + + residuals = actual_arr - predicted_arr + outlier_mask = _tukey_outlier_mask(residuals) + kept = ~outlier_mask + ttest = compute_paired_ttest(actual_arr, predicted_arr) + return { + **base, + "mean_abs_filtered": float(np.mean(np.abs(residuals[kept]))) + if kept.any() + else float("nan"), + "n_outliers": int(np.sum(outlier_mask)), + "ttest_statistic": ttest["statistic"], + "ttest_p_value": ttest["p_value"], + "cohens_d": _paired_cohens_d(actual_arr, predicted_arr), + } + + +def _cv_subset_diagnostics( + run_dir: Path, + training_file: Path, + input_vars: list[str], + output_response: str, + N_CROSS_VALIDATION: int, + n: int, + keep_idxs: list[int] | None, + tag: str, +) -> dict[str, Any]: + if keep_idxs is not None: + subset_file = process_input_file( + training_file, + columns_to_keep=input_vars + [output_response], + suffix=f"convergence_{tag}", + keep_idxs=keep_idxs, + ) + else: + subset_file = process_input_file( + training_file, + columns_to_keep=input_vars + [output_response], + suffix=f"convergence_{tag}", + filter_N_samples=n, + ) + subset_run_dir = run_dir / f"convergence_{tag}" + os.makedirs(subset_run_dir, exist_ok=True) + result = evaluate_sumo_manual_crossvalidation( + subset_run_dir, + subset_file, + input_vars, + output_response, + N_CROSS_VALIDATION=min(N_CROSS_VALIDATION, n), + ) + return compute_cv_diagnostics( + result[output_response], + result[output_response + "_hat"], + result[output_response + "_std_hat"], + ) + + def compute_cv_convergence( run_dir: Path, training_file: Path, @@ -684,43 +904,171 @@ def compute_cv_convergence( N_CROSS_VALIDATION: int = 5, min_samples: int = 5, max_points: int = 5, -) -> list[dict[str, float]]: + n_bootstrap: int = 1, + seed: int = 42, +) -> list[dict[str, Any]]: """Rerun manual K-fold CV at increasing training-sample-count subsets. Reuses `evaluate_sumo_manual_crossvalidation` (the same compute path - `/sumo_cross_validation` already runs) on the first `n` rows of `training_file` for - each subset size, deriving RMSE via `compute_cv_accuracy_metrics` at each step. - Subset sizes are evenly spaced between `min_samples` and the full sample count, - capped at `max_points` to bound the number of extra Dakota reruns (⊥ single-N - snapshot only). Returns a `{n_samples, metric}` series for accuracy-vs-N plotting. + `/sumo_cross_validation` already runs), deriving RMSE via + `compute_cv_accuracy_metrics` at each step. Subset sizes are evenly spaced + between `min_samples` and the full sample count, capped at `max_points` to + bound the number of extra Dakota reruns. Returns a + `{n_samples, metric, metric_std, n_bootstrap, draws, diagnostics}` series for + accuracy-vs-N plotting -- `draws` is the raw (pre-aggregation) RMSE per + bootstrap draw at that size (length 1 when `n_bootstrap<=1` or at the + full-sample point), for callers that want the individual points rather + than just the mean/std (e.g. fitting a trend across the pooled draws). + `diagnostics` is the full `compute_cv_diagnostics` bundle per draw (V16wq) -- + MAE/Tukey-outlier-count/paired-ttest/Cohen's-d/raw actual+predicted -- so a + caller wanting a different metric never needs a second CV rerun. + + `n_bootstrap` (default 1, i.e. a single deterministic first-`n`-rows + subset per size) draws that many *distinct* random subsets (no + replacement) of each subset size below the full sample count and reports + the mean/std RMSE across them via a seeded `np.random.Generator` -- one + arbitrary subset is a noisy point estimate of "does more data help"; + bootstrapping estimates the trend with its own spread. The full-sample + point always runs once (`metric_std=0.0`): there is only one subset of + that size, so bootstrapping it is meaningless. """ n_total = len(load_data(training_file)) subset_sizes = _convergence_subset_sizes(n_total, min_samples, max_points) + rng = np.random.default_rng(seed) series = [] for n in subset_sizes: - subset_file = process_input_file( - training_file, - columns_to_keep=input_vars + [output_response], - filter_N_samples=n, - suffix=f"convergence_{n}", + if n >= n_total or n_bootstrap <= 1: + diag = _cv_subset_diagnostics( + run_dir, + training_file, + input_vars, + output_response, + N_CROSS_VALIDATION, + n, + None, + str(n), + ) + rmse = diag["root_mean_squared"] + series.append( + { + "n_samples": n, + "metric": rmse, + "metric_std": 0.0, + "n_bootstrap": 1, + "draws": [rmse], + "diagnostics": [diag], + } + ) + continue + + draw_diagnostics = [ + _cv_subset_diagnostics( + run_dir, + training_file, + input_vars, + output_response, + N_CROSS_VALIDATION, + n, + rng.choice(n_total, size=n, replace=False).tolist(), + f"{n}_draw{draw}", + ) + for draw in range(n_bootstrap) + ] + draw_rmses = [diag["root_mean_squared"] for diag in draw_diagnostics] + # every draw can come back NaN (e.g. too few points/fold for the surrogate to + # train at all) -- nanmean/nanstd on an all-NaN slice warn but still return NaN, + # so short-circuit instead of letting that warning through on an expected case. + valid_rmses = [r for r in draw_rmses if not np.isnan(r)] + series.append( + { + "n_samples": n, + "metric": float(np.mean(valid_rmses)) if valid_rmses else float("nan"), + "metric_std": float(np.std(valid_rmses)) + if valid_rmses + else float("nan"), + "n_bootstrap": n_bootstrap, + "draws": draw_rmses, + "diagnostics": draw_diagnostics, + } ) - subset_run_dir = run_dir / f"convergence_{n}" - os.makedirs(subset_run_dir, exist_ok=True) - n_folds = min(N_CROSS_VALIDATION, n) - result = evaluate_sumo_manual_crossvalidation( - subset_run_dir, - subset_file, - input_vars, - output_response, - N_CROSS_VALIDATION=n_folds, + + return series + + +def fit_convergence_exponential( + n_samples: list[float], values: list[float] +) -> dict[str, float]: + """Fit `y = a * exp(-b * n)` to pooled convergence data via nonlinear least + squares, reporting the fit's `r_squared`. + + Meant to run on every individual bootstrap draw as its own `(n, value)` + pair (e.g. `compute_cv_convergence`'s `draws` per size, flattened; see + NIH worked example § convergence), not just the per-size mean -- pooling + the raw draws means `r_squared` reflects fit quality against the actual + spread, not an already-smoothed curve. `nan` entries in `values` (e.g. an + all-NaN bootstrap draw) are dropped before fitting. + """ + n_arr = np.asarray(n_samples, dtype=float) + y_arr = np.asarray(values, dtype=float) + mask = ~np.isnan(y_arr) + n_arr, y_arr = n_arr[mask], y_arr[mask] + if n_arr.size < 3: + raise ValueError( + "Need at least 3 valid (n_samples, value) pairs to fit a 2-parameter exponential" ) - metrics = compute_cv_accuracy_metrics( - result[output_response], result[output_response + "_hat"] + + def model(n: np.ndarray, a: float, b: float) -> np.ndarray: + return a * np.exp(-b * n) + + a0 = y_arr[np.argmin(n_arr)] + p0 = (a0 if a0 != 0 else 1.0, 1.0 / max(float(n_arr.max()), 1.0)) + (a, b), _ = curve_fit(model, n_arr, y_arr, p0=p0, maxfev=10000) + + y_pred = model(n_arr, a, b) + ss_res = float(np.sum((y_arr - y_pred) ** 2)) + ss_tot = float(np.sum((y_arr - np.mean(y_arr)) ** 2)) + r_squared = 1.0 - ss_res / ss_tot if ss_tot > 0 else float("nan") + return {"a": float(a), "b": float(b), "r_squared": r_squared} + + +def fit_convergence_exponential_asymptotic( + n_samples: list[float], values: list[float] +) -> dict[str, float]: + """Fit `y = a * exp(-b * n) + c` to pooled convergence data (V18wp). + + `fit_convergence_exponential`'s 2-parameter `a*exp(-b*n)` model decays to + a ZERO asymptote, which is the right shape for a quantity defined to be 0 + at the reference point (e.g. % deviation from a final value) but the + WRONG shape for a raw error metric (MAE/RMSE in physical units): even at + N->infinity a surrogate keeps some irreducible error, so the curve should + flatten at a positive floor `c`, not at 0. Otherwise identical to + `fit_convergence_exponential` (nonlinear least squares, `nan` values + dropped, needs >=3 points -- one more than the 2-param fit since there + are 3 parameters to identify). + """ + n_arr = np.asarray(n_samples, dtype=float) + y_arr = np.asarray(values, dtype=float) + mask = ~np.isnan(y_arr) + n_arr, y_arr = n_arr[mask], y_arr[mask] + if n_arr.size < 4: + raise ValueError( + "Need at least 4 valid (n_samples, value) pairs to fit a 3-parameter exponential" ) - series.append({"n_samples": n, "metric": metrics["root_mean_squared"]}) - return series + def model(n: np.ndarray, a: float, b: float, c: float) -> np.ndarray: + return a * np.exp(-b * n) + c + + c0 = y_arr[np.argmax(n_arr)] + a0 = y_arr[np.argmin(n_arr)] - c0 + p0 = (a0 if a0 != 0 else 1.0, 1.0 / max(float(n_arr.max()), 1.0), c0) + (a, b, c), _ = curve_fit(model, n_arr, y_arr, p0=p0, maxfev=10000) + + y_pred = model(n_arr, a, b, c) + ss_res = float(np.sum((y_arr - y_pred) ** 2)) + ss_tot = float(np.sum((y_arr - np.mean(y_arr)) ** 2)) + r_squared = 1.0 - ss_res / ss_tot if ss_tot > 0 else float("nan") + return {"a": float(a), "b": float(b), "c": float(c), "r_squared": r_squared} def evaluate_sumo( diff --git a/tests/test_dakota_funs_data_processing.py b/tests/test_dakota_funs_data_processing.py index a311da2..e3371cc 100644 --- a/tests/test_dakota_funs_data_processing.py +++ b/tests/test_dakota_funs_data_processing.py @@ -9,6 +9,7 @@ from itis_sumo.data.funs_data_processing import ( _filter_data, _parse_data, + auto_select_distributions, create_manual_uq_samples, get_bounds_uniform_distribution, get_bounds_uniform_distributions, @@ -18,6 +19,7 @@ is_dominated, load_data, sanitize_varnames, + select_variable_scale, ) # --- sanitize_varnames ------------------------------------------------------- @@ -505,3 +507,92 @@ def test_filter_data_both_filters_raises(): df = pd.DataFrame({"a": [1, 2, 3]}) with pytest.raises(AssertionError, match="only one of"): _filter_data(df, filter_highest_N=1, filter_N_samples=1) + + +# --- select_variable_scale / auto_select_distributions ----------------------- + + +class TestSelectVariableScale: + def test_constant_short_circuits_without_running_tests(self): + result = select_variable_scale([5.0] * 20) + assert result == { + "scale": "linear", + "distribution": "constant", + "value": 5.0, + "confident": True, + "candidates": {}, + } + + def test_recovers_normal_linear_on_normal_data(self): + rng = np.random.default_rng(2) + values = rng.normal(loc=10.0, scale=2.0, size=100) + result = select_variable_scale(values) + assert result["scale"] == "linear" + assert result["distribution"] == "normal" + assert result["confident"] + assert set(result["candidates"]) == { + "linear_normal", + "linear_uniform", + "log_normal", + "log_uniform", + } + + def test_recovers_uniform_linear_on_uniform_data(self): + rng = np.random.default_rng(1) + values = rng.uniform(low=1.0, high=2.0, size=500) + result = select_variable_scale(values) + assert result["scale"] == "linear" + assert result["distribution"] == "uniform" + + def test_recovers_normal_log_on_lognormal_data(self): + rng = np.random.default_rng(2) + values = rng.lognormal(mean=0.0, sigma=0.5, size=500) + result = select_variable_scale(values) + assert result["scale"] == "log" + assert result["distribution"] == "normal" + + def test_log_candidates_omitted_when_values_nonpositive(self): + values = np.array([-1.0, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0]) + result = select_variable_scale(values) + assert "log_normal" not in result["candidates"] + assert "log_uniform" not in result["candidates"] + assert result["scale"] == "linear" + + def test_p_value_matches_best_candidate(self): + rng = np.random.default_rng(3) + values = rng.normal(loc=0.0, scale=1.0, size=200) + result = select_variable_scale(values) + best_key = f"{result['scale']}_{result['distribution']}" + assert result["p_value"] == result["candidates"][best_key]["p_value"] + + +class TestAutoSelectDistributions: + def test_builds_usable_distributions_dict_with_log_prefix(self): + rng = np.random.default_rng(4) + df = pd.DataFrame( + { + "x_normal": rng.normal(10.0, 2.0, 300), + "x_lognormal": rng.lognormal(0.0, 0.5, 300), + "x_constant": [3.0] * 300, + } + ) + distributions, diagnostics = auto_select_distributions(df, list(df.columns)) + + assert distributions["x_normal"]["distribution"] == "normal" + assert "log_x_lognormal" in distributions + assert distributions["log_x_lognormal"]["distribution"] == "normal" + assert distributions["x_constant"] == {"distribution": "constant", "value": 3.0} + assert set(diagnostics) == {"x_normal", "x_lognormal", "x_constant"} + + # the resulting dict must be directly usable by create_manual_uq_samples + samples = create_manual_uq_samples( + list(distributions.keys()), distributions, num_samples=5, seed=1 + ) + assert set(samples) == set(distributions) + + def test_diagnostics_carry_confidence_and_p_value_per_variable(self): + rng = np.random.default_rng(5) + df = pd.DataFrame({"x": rng.normal(0.0, 1.0, 200)}) + _, diagnostics = auto_select_distributions(df, ["x"]) + assert "confident" in diagnostics["x"] + assert "p_value" in diagnostics["x"] diff --git a/tests/test_funs_evaluate_cv_stats.py b/tests/test_funs_evaluate_cv_stats.py index 177d780..0463252 100644 --- a/tests/test_funs_evaluate_cv_stats.py +++ b/tests/test_funs_evaluate_cv_stats.py @@ -10,9 +10,13 @@ from itis_sumo.evaluate.funs_evaluate import ( _convergence_subset_sizes, + compute_coverage, compute_cv_accuracy_metrics, compute_cv_convergence, + compute_cv_diagnostics, compute_paired_ttest, + fit_convergence_exponential, + fit_convergence_exponential_asymptotic, ) @@ -104,6 +108,160 @@ def test_nan_pair_from_dropped_cv_row_is_excluded(self): assert result["p_value"] == pytest.approx(1.0) +class TestComputeCvDiagnostics: + def test_returns_rmse_mae_ttest_cohens_d_in_one_call(self): + """§T19df/V16wq: a single call surfaces every metric the convergence + study needs -- no second CV rerun required for a different metric.""" + actual = [1.0, 2.0, 3.0, 4.0] + predicted = [2.0, 1.0, 4.0, 3.0] # residuals -1,+1,-1,+1 -> mean 0 + diag = compute_cv_diagnostics(actual, predicted) + assert diag["root_mean_squared"] == pytest.approx(1.0) + assert diag["mean_abs"] == pytest.approx(1.0) + assert diag["ttest_statistic"] == pytest.approx(0.0) + assert diag["ttest_p_value"] == pytest.approx(1.0) + assert diag["cohens_d"] == pytest.approx(0.0) + assert diag["actual"] == pytest.approx(actual) + assert diag["predicted"] == pytest.approx(predicted) + + def test_tukey_flags_extreme_residual(self): + """One wildly-off pair among otherwise-tight residuals gets flagged and + excluded from the filtered MAE, without touching the unfiltered `mean_abs` + (V17kb: unfiltered MAE stays primary, filtering is a secondary view).""" + actual = [1.0, 2.0, 3.0, 4.0, 5.0, 100.0] + predicted = [1.1, 1.9, 3.1, 3.9, 5.1, 0.0] # last pair: residual +100 + diag = compute_cv_diagnostics(actual, predicted) + assert diag["n_outliers"] == 1 + assert diag["mean_abs_filtered"] < diag["mean_abs"] + + def test_no_systematic_bias_gives_small_effect_size(self): + rng = np.random.default_rng(42) + actual = rng.normal(loc=0.0, scale=1.0, size=200) + predicted = actual + rng.normal(loc=0.0, scale=0.01, size=200) + diag = compute_cv_diagnostics(actual, predicted) + assert abs(diag["cohens_d"]) < 0.2 + + def test_constant_bias_gives_large_effect_size_regardless_of_n(self): + """Cohen's d (unlike the t-test p-value) isn't confounded by sample size -- + the same constant bias yields ~the same effect size at N=20 and N=200.""" + rng = np.random.default_rng(7) + for n in (20, 200): + actual = rng.normal(loc=0.0, scale=1.0, size=n) + predicted = actual + 2.0 + diag = compute_cv_diagnostics(actual, predicted) + assert abs(diag["cohens_d"]) > 0.8 + + def test_fewer_than_two_points_returns_nan_diagnostics_not_raise(self): + diag = compute_cv_diagnostics([1.0], [2.0]) + assert np.isnan(diag["ttest_p_value"]) + assert np.isnan(diag["cohens_d"]) + assert diag["n_outliers"] == 0 + + def test_shape_mismatch_raises(self): + with pytest.raises(ValueError, match="same shape"): + compute_cv_diagnostics([1.0, 2.0], [1.0, 2.0, 3.0]) + + def test_predicted_std_passed_through_when_provided(self): + """T24bn/V19cz: predicted_std rides along in the bundle so downstream + coverage/calibration checks don't need a second CV rerun.""" + actual = [1.0, 2.0, 3.0, 4.0] + predicted = [2.0, 1.0, 4.0, 3.0] + predicted_std = [0.5, 0.6, 0.7, 0.8] + diag = compute_cv_diagnostics(actual, predicted, predicted_std) + assert diag["predicted_std"] == pytest.approx(predicted_std) + + def test_predicted_std_defaults_to_nan_when_omitted(self): + """Back-compat: existing 2-arg callers keep working, predicted_std + just comes back as NaN instead of a real value.""" + actual = [1.0, 2.0, 3.0, 4.0] + predicted = [2.0, 1.0, 4.0, 3.0] + diag = compute_cv_diagnostics(actual, predicted) + assert len(diag["predicted_std"]) == len(actual) + assert all(np.isnan(v) for v in diag["predicted_std"]) + + def test_predicted_std_filtered_by_same_nan_mask_as_actual_predicted(self): + """B22-style: a dropped CV row's NaN in actual/predicted must also drop + the corresponding predicted_std entry, keeping all three arrays aligned.""" + actual = [1.0, 2.0, 3.0, 4.0] + predicted = [2.0, 1.0, 4.0, float("nan")] + predicted_std = [0.5, 0.6, 0.7, 0.8] + diag = compute_cv_diagnostics(actual, predicted, predicted_std) + assert diag["predicted_std"] == pytest.approx([0.5, 0.6, 0.7]) + + +class TestComputeCoverage: + def test_well_calibrated_gaussian_matches_nominal_levels(self): + rng = np.random.default_rng(0) + n = 5000 + actual = rng.normal(loc=0.0, scale=1.0, size=n) + predicted = np.zeros(n) + predicted_std = np.ones(n) + result = compute_coverage(actual, predicted, predicted_std) + assert result["n_points"] == n + for empirical, nominal in zip(result["empirical"], result["levels"]): + assert empirical == pytest.approx(nominal, abs=0.02) + + def test_overconfident_std_undercounts_coverage(self): + """A predicted_std much smaller than the actual residual spread should + show up as empirical coverage well below the nominal level -- the + exact overconfidence pattern V19cz exists to catch.""" + rng = np.random.default_rng(1) + n = 2000 + actual = rng.normal(loc=0.0, scale=1.0, size=n) + predicted = np.zeros(n) + predicted_std = np.full(n, 0.2) # much narrower than the true scale=1.0 + result = compute_coverage(actual, predicted, predicted_std, levels=(0.6827,)) + assert result["empirical"][0] < 0.6827 - 0.05 + + def test_shape_mismatch_raises(self): + with pytest.raises(ValueError, match="same shape"): + compute_coverage([1.0, 2.0], [1.0, 2.0, 3.0], [1.0, 1.0, 1.0]) + + def test_nan_and_nonpositive_std_excluded(self): + actual = [1.0, 2.0, 3.0, 4.0] + predicted = [1.0, 2.0, 3.0, 4.0] + predicted_std = [1.0, float("nan"), 0.0, -1.0] + result = compute_coverage(actual, predicted, predicted_std, levels=(0.6827,)) + assert result["n_points"] == 1 + + def test_empty_after_filtering_returns_nan_not_raise(self): + result = compute_coverage([1.0], [1.0], [float("nan")], levels=(0.6827,)) + assert result["n_points"] == 0 + assert np.isnan(result["empirical"][0]) + + +class TestFitConvergenceExponentialAsymptotic: + def test_recovers_known_params_with_nonzero_floor(self): + """V18wp: raw error metrics have an irreducible-error floor as N->infinity, + unlike the old zero-asymptote fit -- this is why the 3-parameter model exists.""" + a_true, b_true, c_true = 2.0, 0.1, 0.5 + n = np.linspace(10, 100, 40) + y = a_true * np.exp(-b_true * n) + c_true + + fit = fit_convergence_exponential_asymptotic(n.tolist(), y.tolist()) + + assert fit["a"] == pytest.approx(a_true, rel=1e-2) + assert fit["b"] == pytest.approx(b_true, rel=1e-2) + assert fit["c"] == pytest.approx(c_true, rel=1e-2) + assert fit["r_squared"] == pytest.approx(1.0, abs=1e-6) + + def test_nan_values_are_dropped_before_fitting(self): + a_true, b_true, c_true = 3.0, 0.1, 0.2 + n = np.linspace(5, 50, 10) + y = a_true * np.exp(-b_true * n) + c_true + y_with_nans = y.tolist() + y_with_nans[2] = float("nan") + y_with_nans[7] = float("nan") + + fit = fit_convergence_exponential_asymptotic(n.tolist(), y_with_nans) + + assert fit["a"] == pytest.approx(a_true, rel=1e-2) + assert fit["c"] == pytest.approx(c_true, rel=1e-2) + + def test_raises_below_four_valid_points(self): + with pytest.raises(ValueError): + fit_convergence_exponential_asymptotic([10.0, 20.0, 30.0], [1.0, 0.5, 0.3]) + + class TestConvergenceSubsetSizes: def test_below_minimum_returns_empty(self): assert _convergence_subset_sizes(n_total=3, min_samples=5, max_points=5) == [] @@ -169,6 +327,14 @@ def fake_manual_cv( assert call_sizes[-1] == n_total for point in series: assert point["metric"] == pytest.approx(0.1) + # V16wq: diagnostics bundle rides along, one entry per draw, no rerun + assert len(point["diagnostics"]) == point["n_bootstrap"] + assert point["diagnostics"][0]["root_mean_squared"] == pytest.approx(0.1) + # T24bn/V19cz: std_hat rides along through the convergence series too + assert all( + v == pytest.approx(0.0) + for v in point["diagnostics"][0]["predicted_std"] + ) def test_empty_series_when_below_minimum(self, tmp_path): training_file = tmp_path / "df_processed_jobs.dat" @@ -181,3 +347,207 @@ def test_empty_series_when_below_minimum(self, tmp_path): tmp_path, training_file, ["x1"], "y", min_samples=5 ) assert series == [] + + +class TestComputeCvConvergenceBootstrap: + def test_bootstrap_draws_multiple_distinct_subsets_below_full_size( + self, tmp_path, monkeypatch + ): + training_file = tmp_path / "df_processed_jobs.dat" + n_total = 10 + rng = np.random.default_rng(1) + x1 = rng.uniform(-1, 1, n_total) + y = x1 * 2.0 + with open(training_file, "w") as f: + f.write("x1 y\n") + f.writelines(f"{xi} {yi}\n" for xi, yi in zip(x1, y)) + + seen_subsets = [] + + def fake_manual_cv( + run_dir, subset_file, input_vars, output_response, N_CROSS_VALIDATION=5 + ): + import pandas as pd + + df = pd.read_csv(subset_file, sep=" ") + seen_subsets.append(tuple(sorted(df["x1"].round(6).tolist()))) + actual = df[output_response].astype(float).tolist() + predicted = [v + 0.1 for v in actual] + return { + output_response: actual, + output_response + "_hat": predicted, + output_response + "_std_hat": [0.0] * len(actual), + } + + monkeypatch.setattr( + "itis_sumo.evaluate.funs_evaluate.evaluate_sumo_manual_crossvalidation", + fake_manual_cv, + ) + + series = compute_cv_convergence( + tmp_path, + training_file, + ["x1"], + "y", + min_samples=5, + max_points=2, + n_bootstrap=4, + ) + + # sizes: [5, 10] -- 4 bootstrap draws at n=5, exactly 1 run at n=10 (full set) + assert [point["n_samples"] for point in series] == [5, 10] + assert len(seen_subsets) == 4 + 1 + + small_point, full_point = series + assert small_point["n_bootstrap"] == 4 + assert full_point["n_bootstrap"] == 1 + assert full_point["metric_std"] == pytest.approx(0.0) + + # the 4 draws at n=5 aren't all the same subset of rows + small_subsets = set(seen_subsets[:4]) + assert len(small_subsets) > 1 + + # raw per-draw values are exposed for pooled downstream analysis + assert len(small_point["draws"]) == 4 + assert small_point["draws"] == pytest.approx([0.1] * 4) + assert full_point["draws"] == pytest.approx([0.1]) + + # V16wq: diagnostics bundle length matches draws length at each size + assert len(small_point["diagnostics"]) == 4 + assert len(full_point["diagnostics"]) == 1 + + def test_bootstrap_is_deterministic_given_seed(self, tmp_path, monkeypatch): + training_file = tmp_path / "df_processed_jobs.dat" + n_total = 10 + rng = np.random.default_rng(2) + x1 = rng.uniform(-1, 1, n_total) + y = x1 * 2.0 + with open(training_file, "w") as f: + f.write("x1 y\n") + f.writelines(f"{xi} {yi}\n" for xi, yi in zip(x1, y)) + + def fake_manual_cv( + run_dir, subset_file, input_vars, output_response, N_CROSS_VALIDATION=5 + ): + import pandas as pd + + df = pd.read_csv(subset_file, sep=" ") + actual = df[output_response].astype(float).tolist() + predicted = [v + 0.1 for v in actual] + return { + output_response: actual, + output_response + "_hat": predicted, + output_response + "_std_hat": [0.0] * len(actual), + } + + monkeypatch.setattr( + "itis_sumo.evaluate.funs_evaluate.evaluate_sumo_manual_crossvalidation", + fake_manual_cv, + ) + + series_a = compute_cv_convergence( + tmp_path / "a", + training_file, + ["x1"], + "y", + min_samples=5, + max_points=2, + n_bootstrap=4, + seed=7, + ) + series_b = compute_cv_convergence( + tmp_path / "b", + training_file, + ["x1"], + "y", + min_samples=5, + max_points=2, + n_bootstrap=4, + seed=7, + ) + assert series_a == series_b + + def test_default_n_bootstrap_matches_pre_bootstrap_behavior( + self, tmp_path, monkeypatch + ): + """n_bootstrap defaults to 1: single deterministic first-n-rows subset, + same as the pre-bootstrap implementation -- back-compat for existing callers.""" + training_file = tmp_path / "df_processed_jobs.dat" + n_total = 10 + rng = np.random.default_rng(3) + x1 = rng.uniform(-1, 1, n_total) + y = x1 * 2.0 + with open(training_file, "w") as f: + f.write("x1 y\n") + f.writelines(f"{xi} {yi}\n" for xi, yi in zip(x1, y)) + + seen_subsets = [] + + def fake_manual_cv( + run_dir, subset_file, input_vars, output_response, N_CROSS_VALIDATION=5 + ): + import pandas as pd + + df = pd.read_csv(subset_file, sep=" ") + seen_subsets.append(df["x1"].round(6).tolist()) + actual = df[output_response].astype(float).tolist() + predicted = [v + 0.1 for v in actual] + return { + output_response: actual, + output_response + "_hat": predicted, + output_response + "_std_hat": [0.0] * len(actual), + } + + monkeypatch.setattr( + "itis_sumo.evaluate.funs_evaluate.evaluate_sumo_manual_crossvalidation", + fake_manual_cv, + ) + + series = compute_cv_convergence( + tmp_path, training_file, ["x1"], "y", min_samples=5, max_points=2 + ) + + assert len(seen_subsets) == 2 # one run per size, no bootstrap fan-out + assert seen_subsets[0] == pytest.approx(x1[:5].round(6).tolist()) + assert all( + point["n_bootstrap"] == 1 and point["metric_std"] == 0.0 for point in series + ) + + +class TestFitConvergenceExponential: + def test_recovers_known_params_on_noiseless_data(self): + a_true, b_true = 5.0, 0.05 + n = np.linspace(10, 100, 30) + y = a_true * np.exp(-b_true * n) + + fit = fit_convergence_exponential(n.tolist(), y.tolist()) + + assert fit["a"] == pytest.approx(a_true, rel=1e-3) + assert fit["b"] == pytest.approx(b_true, rel=1e-3) + assert fit["r_squared"] == pytest.approx(1.0, abs=1e-6) + + def test_nan_values_are_dropped_before_fitting(self): + a_true, b_true = 3.0, 0.1 + n = np.linspace(5, 50, 10) + y = a_true * np.exp(-b_true * n) + y_with_nans = y.tolist() + y_with_nans[2] = float("nan") + y_with_nans[7] = float("nan") + + fit = fit_convergence_exponential(n.tolist(), y_with_nans) + + assert fit["a"] == pytest.approx(a_true, rel=1e-3) + assert fit["b"] == pytest.approx(b_true, rel=1e-3) + + def test_noisy_data_gives_high_but_imperfect_r_squared(self): + rng = np.random.default_rng(0) + n = np.linspace(10, 100, 40) + y = 5.0 * np.exp(-0.05 * n) + rng.normal(0, 0.05, size=n.size) + + fit = fit_convergence_exponential(n.tolist(), y.tolist()) + + assert 0.5 < fit["r_squared"] < 1.0 + + def test_raises_below_three_valid_points(self): + with pytest.raises(ValueError): + fit_convergence_exponential([10.0, 20.0], [1.0, 0.5]) From 362a426b234972a4578c6ad270a35c4d8b494ec6 Mon Sep 17 00:00:00 2001 From: Javier Garcia Ordonez Date: Tue, 29 Sep 2026 16:38:31 +0200 Subject: [PATCH 2/4] feat(data): raw-value Tukey outlier detector + real analyze_dataset (T49np) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The 'does not exist yet, anywhere' gap from BRANCH_CONSOLIDATION.md §3: - detect_raw_outliers: Tukey fence computed in the variable's SELECTED scale (geometric fence for log-selected variables -- arithmetic fences over-flag the benign right tail of multiplicative data), flag-only, NaN-excluded, fences reported in the values' own units - analyze_dataset wired for real (was T18ry scaffold): auto scale/ distribution + per-column outliers + optional candidate detail, dataclasses.asdict() JSON round-trip for the flaskapi boundary Flip tests: geometric-vs-arithmetic over-flagging guard; ln<->log10 affine log-basis moves zero flags (V46np). (cherry picked from commit 21a5a233e5e90f30f5a4da56dd6ab2adbd35d027) --- src/itis_sumo/data/__init__.py | 2 + .../data/funs_dataset_diagnostics.py | 111 ++++++++++-- tests/test_funs_dataset_diagnostics.py | 166 +++++++++++++++++- 3 files changed, 265 insertions(+), 14 deletions(-) diff --git a/src/itis_sumo/data/__init__.py b/src/itis_sumo/data/__init__.py index 275140f..247f3c8 100644 --- a/src/itis_sumo/data/__init__.py +++ b/src/itis_sumo/data/__init__.py @@ -24,6 +24,7 @@ OutlierSummary, VariableDiagnostics, analyze_dataset, + detect_raw_outliers, ) __all__ = [ @@ -36,6 +37,7 @@ "create_grid_samples", "create_manual_uq_samples", "create_samples_along_axes", + "detect_raw_outliers", "extract_predictions_along_axes", "extract_predictions_gridpoints", "get_bounds_uniform_distribution", diff --git a/src/itis_sumo/data/funs_dataset_diagnostics.py b/src/itis_sumo/data/funs_dataset_diagnostics.py index e6bc9f8..e168093 100644 --- a/src/itis_sumo/data/funs_dataset_diagnostics.py +++ b/src/itis_sumo/data/funs_dataset_diagnostics.py @@ -1,20 +1,25 @@ -"""Dataset diagnostics schema: the narrow, versioned entrypoint flaskapi/mmux_vite +"""Dataset diagnostics: the narrow, versioned entrypoint flaskapi/mmux_vite consumes for scale/distribution/outlier information (SPEC.md V16qf, T18ry). -Schema only for now. The underlying detection logic (`select_variable_scale`, -`auto_select_distributions`, a raw-value outlier detector) is not yet promoted -from the confidential incubator branch — see BRANCH_CONSOLIDATION.md §3 for -per-piece readiness. `analyze_dataset` raises `NotImplementedError` until that -logic is wired in. +Detection is stitched from pieces promoted off the confidential incubator +branch: `select_variable_scale`/`auto_select_distributions` (T48np) for the +scale/distribution half, and `detect_raw_outliers` (this module, V46np) for +the outlier-surfacing half — the raw-value Tukey detector BRANCH_CONSOLIDATION.md +§3 recorded as "new work, not a port". Flag-only, like `_tukey_outlier_mask` +does for CV residuals (V17kb): nothing here drops data, callers decide. """ from __future__ import annotations +from collections.abc import Sequence from dataclasses import dataclass from typing import Any, Literal +import numpy as np import pandas as pd +from itis_sumo.data.funs_data_processing import auto_select_distributions + Scale = Literal["linear", "log"] Distribution = Literal["constant", "uniform", "normal"] @@ -43,6 +48,66 @@ class DatasetDiagnostics: detail: dict[str, Any] | None = None +def detect_raw_outliers( + values: Sequence[float] | np.ndarray, + *, + scale: Scale = "linear", + k: float = 1.5, +) -> OutlierSummary: + """Flag raw column values outside the Tukey IQR fence (Q1 - k*IQR, + Q3 + k*IQR) — the raw-column counterpart of `evaluate`'s private + `_tukey_outlier_mask` (which operates on CV residuals only, V17kb). + Flags, never drops: the caller decides what a flag means. + + `scale` is the space the fence is *computed* in (V46np): `"log"` requires + strictly-positive values, computes the fence in ln-space and returns + `fence_low`/`fence_high` mapped back (`exp`) into the values' own units. + That makes the log fence the GEOMETRIC analogue `(q1·(q1/q3)^k, + q3·(q3/q1)^k)` — deliberately NOT the arithmetic fence re-expressed, + which is the point: linear-space fences over-flag the benign right tail + of multiplicative/lognormal quantities (the flip test flags ~4x the + injected count), so imposing them on a log-selected variable manufactures + phantom outliers. Quantiles DO commute with the affine log-basis change + (ln = ln10·log10), so the log base can never change which rows are + flagged — asserted, not assumed. + + NaN entries are excluded from the quartile estimate (nanpercentile) and + are never flagged (comparison against a finite fence is False for NaN). + A zero-IQR column yields coincident fences; only values strictly outside + are flagged (so a constant column flags nothing). + """ + arr = np.asarray(values, dtype=float) + if scale not in ("linear", "log"): + raise ValueError(f"Unknown scale: {scale!r}") + finite = np.isfinite(arr) + if scale == "log": + if not np.all(arr[finite] > 0): + raise ValueError( + "detect_raw_outliers(scale='log') requires strictly positive " + "values; non-positive rows are a scale-selection error" + ) + with np.errstate(divide="ignore"): + x = np.where(finite, np.log(np.where(finite, arr, 1.0)), np.nan) + else: + x = arr + + q1, q3 = np.nanpercentile(x, [25, 75]) + iqr = q3 - q1 + lower, upper = q1 - k * iqr, q3 + k * iqr + flagged = finite & ((x < lower) | (x > upper)) + + if scale == "log": + fence_low, fence_high = float(np.exp(lower)), float(np.exp(upper)) + else: + fence_low, fence_high = float(lower), float(upper) + return OutlierSummary( + count=int(np.sum(flagged)), + indices=[int(i) for i in np.nonzero(flagged)[0]], + fence_low=fence_low, + fence_high=fence_high, + ) + + def analyze_dataset( df: pd.DataFrame, input_cols: list[str], @@ -54,6 +119,13 @@ def analyze_dataset( on a training dataset — the only itis-sumo call flaskapi should need for this feature (V16qf). + Scale/distribution come from `auto_select_distributions` (the same + decision the incubator's NIH/Merck studies used); outliers from + `detect_raw_outliers` computed in each variable's own selected scale + (V46np). Detection is per column of `df` as given — if the caller keeps + a hand-built `log_`-prefixed column, that column's raw (log-space) + values are what get inspected. + Args: df: training data, one column per variable. input_cols: names of `df` columns to treat as inputs. @@ -66,8 +138,27 @@ def analyze_dataset( Returns: `DatasetDiagnostics`, JSON-serializable via `dataclasses.asdict()`. """ - raise NotImplementedError( - "analyze_dataset is scaffolded (schema only, T18ry) — scale/distribution/" - "outlier detection logic is not yet promoted from feat/nih-in-silico-example; " - "see BRANCH_CONSOLIDATION.md §3" + all_cols = list(input_cols) + list(output_cols) + _, decisions = auto_select_distributions(df, all_cols, alpha=alpha) + + def per_column(col: str) -> VariableDiagnostics: + decision = decisions[col] + scale: Scale = ( + decision["scale"] if decision["distribution"] != "constant" else "linear" + ) + values = df[col].to_numpy(dtype=float) + return VariableDiagnostics( + name=col, + scale=scale, + distribution=decision["distribution"], + confident=bool(decision["confident"]), + outliers=detect_raw_outliers(values, scale=scale), + ) + + return DatasetDiagnostics( + inputs={col: per_column(col) for col in input_cols}, + outputs={col: per_column(col) for col in output_cols}, + detail={col: dict(decisions[col]) for col in all_cols} + if include_detail + else None, ) diff --git a/tests/test_funs_dataset_diagnostics.py b/tests/test_funs_dataset_diagnostics.py index a80e660..ec2a449 100644 --- a/tests/test_funs_dataset_diagnostics.py +++ b/tests/test_funs_dataset_diagnostics.py @@ -1,10 +1,168 @@ +import dataclasses +import json + +import numpy as np import pandas as pd import pytest from itis_sumo.data import analyze_dataset +from itis_sumo.data.funs_dataset_diagnostics import detect_raw_outliers + +# --- detect_raw_outliers ------------------------------------------------------ + + +def test_flags_single_extreme_value(): + values = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 100.0] + summary = detect_raw_outliers(values) + assert summary.count == 1 + assert summary.indices == [9] + assert summary.fence_high < 100.0 + + +def test_constant_column_flags_nothing_with_coincident_fences(): + summary = detect_raw_outliers([2.5] * 20) + assert summary.count == 0 + assert summary.indices == [] + assert summary.fence_low == summary.fence_high == 2.5 + + +def test_nan_rows_never_flagged_and_excluded_from_quartiles(): + clean = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0] + with_nan = [1.0, 2.0, np.nan, 3.0, 4.0, 5.0, np.nan, 6.0] + assert detect_raw_outliers(with_nan).indices == detect_raw_outliers(clean).indices + # a NaN is never itself flagged + spiky = [1.0, 2.0, 3.0, np.nan, 50.0] + assert 3 not in detect_raw_outliers(spiky).indices + + +def test_log_scale_requires_strictly_positive_values(): + with pytest.raises(ValueError, match="strictly positive"): + detect_raw_outliers([1.0, 0.0, 2.0], scale="log") + with pytest.raises(ValueError, match="strictly positive"): + detect_raw_outliers([-1.0, 2.0], scale="log") + + +def test_unknown_scale_raises(): + with pytest.raises(ValueError, match="Unknown scale"): + detect_raw_outliers([1.0, 2.0], scale="sqrt") # ty: ignore[invalid-argument-type] + + +class TestScaleAwareFences: + """V46np: the fence is computed in the variable's SELECTED scale. For a + log-scale (multiplicative) variable that means the geometric fence — + `exp(k-th Tukey fence of ln values)` — not the arithmetic fence in + disguise; and the log-basis choice (ln vs log10, an affine change on the + log variable) must never move a single flag.""" + + def test_log_fence_is_the_geometric_analogue(self): + rng = np.random.default_rng(0) + values = rng.lognormal(mean=1.0, sigma=1.2, size=400) + q1_ln, q3_ln = np.percentile(np.log(values), [25, 75]) + iqr_ln = q3_ln - q1_ln + summary = detect_raw_outliers(values, scale="log") + assert summary.fence_low == pytest.approx( + np.exp(q1_ln - 1.5 * iqr_ln), rel=1e-12 + ) + assert summary.fence_high == pytest.approx( + np.exp(q3_ln + 1.5 * iqr_ln), rel=1e-12 + ) + # geometric (multiplicative) form, to within quantile-interpolation noise: + q1, q3 = np.percentile(values, [25, 75]) + assert summary.fence_low == pytest.approx(q1 * (q1 / q3) ** 1.5, rel=1e-4) + assert summary.fence_high == pytest.approx(q3 * (q3 / q1) ** 1.5, rel=1e-4) + + def test_linear_fence_on_lognormal_data_overflags_vs_geometric_fence(self): + # the flip test: 3 injected outliers among benign multiplicative spread. + rng = np.random.default_rng(1) + values = rng.lognormal(mean=1.0, sigma=1.2, size=200) + values[rng.choice(200, size=3, replace=False)] *= 50.0 + linear = detect_raw_outliers(values, scale="linear") + log = detect_raw_outliers(values, scale="log") + assert log.count >= 3 # the injected outliers are found + assert linear.count > log.count + 3 # linear manufactures extras + + @pytest.mark.parametrize("seed", [0, 1, 2]) + def test_log_base_choice_moves_no_flags(self, seed): + rng = np.random.default_rng(seed) + values = rng.lognormal(mean=1.0, sigma=1.2, size=200) + values[rng.choice(200, size=3, replace=False)] *= 50.0 + ln_space = detect_raw_outliers(values, scale="log") + log10_space = detect_raw_outliers(np.log10(values), scale="linear") + assert ln_space.count == log10_space.count + assert ln_space.indices == log10_space.indices + assert ln_space.fence_low == pytest.approx( + 10.0**log10_space.fence_low, rel=1e-9 + ) + assert ln_space.fence_high == pytest.approx( + 10.0**log10_space.fence_high, rel=1e-9 + ) + + def test_log_fences_reported_in_original_units(self): + rng = np.random.default_rng(3) + values = rng.lognormal(mean=0.0, sigma=0.5, size=300) + summary = detect_raw_outliers(values, scale="log") + assert summary.fence_low > 0.0 + assert summary.fence_low < np.median(values) < summary.fence_high + assert summary.count <= max( + 1, len(values) // 100 + ) # benign spread, tail-sized flags only + + +# --- analyze_dataset ---------------------------------------------------------- + + +def _synthetic_df(n=200, seed=7): + rng = np.random.default_rng(seed) + return pd.DataFrame( + { + "x_normal": rng.normal(10.0, 2.0, n), + "x_lognormal": rng.lognormal(0.0, 0.5, n), + "x_constant": [3.0] * n, + "y_out": rng.normal(0.0, 1.0, n), + } + ) + + +def test_analyze_dataset_returns_diagnostics_for_every_column(): + df = _synthetic_df() + diag = analyze_dataset( + df, input_cols=["x_normal", "x_lognormal", "x_constant"], output_cols=["y_out"] + ) + assert set(diag.inputs) == {"x_normal", "x_lognormal", "x_constant"} + assert set(diag.outputs) == {"y_out"} + assert diag.outputs["y_out"].scale in ("linear", "log") + assert diag.outputs["y_out"].outliers is not None + assert diag.inputs["x_constant"].distribution == "constant" + assert diag.inputs["x_constant"].confident is True + + +def test_analyze_dataset_selects_log_scale_for_lognormal_column(): + df = _synthetic_df() + diag = analyze_dataset(df, input_cols=["x_lognormal"], output_cols=["y_out"]) + assert diag.inputs["x_lognormal"].scale == "log" + assert diag.inputs["x_lognormal"].distribution == "normal" + + +def test_analyze_dataset_surfaces_injected_raw_outliers(): + df = _synthetic_df() + df.loc[10, "x_lognormal"] *= 100.0 + diag = analyze_dataset(df, input_cols=["x_lognormal"], output_cols=["y_out"]) + summary = diag.inputs["x_lognormal"].outliers + assert summary is not None + assert summary.count >= 1 + assert 10 in summary.indices + +def test_analyze_dataset_detail_only_with_flag_and_json_serializable(): + df = _synthetic_df(n=100, seed=8) + lean = analyze_dataset(df, input_cols=["x_normal"], output_cols=["y_out"]) + assert lean.detail is None -def test_analyze_dataset_raises_not_implemented_until_detection_logic_lands(): - df = pd.DataFrame({"x": [1.0, 2.0, 3.0], "y": [4.0, 5.0, 6.0]}) - with pytest.raises(NotImplementedError, match="T18ry"): - analyze_dataset(df, input_cols=["x"], output_cols=["y"]) + full = analyze_dataset( + df, input_cols=["x_normal"], output_cols=["y_out"], include_detail=True + ) + assert full.detail is not None + assert "candidates" in full.detail["x_normal"] + # the whole point of the plain-dataclass schema: crosses the flaskapi boundary + payload = json.dumps(dataclasses.asdict(full)) + assert "x_normal" in payload From 47fdcd79d983fbfe927253e56f97d69c750921af Mon Sep 17 00:00:00 2001 From: Javier Garcia Ordonez Date: Thu, 1 Oct 2026 16:13:44 +0200 Subject: [PATCH 3/4] feat(report): itis_sumo.report facade for automated report notebooks (T54rp) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Report code (marimo report notebooks, tissue-conductivity LHS instance) reached into itis_sumo.{api,data,evaluate.funs_evaluate} internals in parallel; every engine refactor then forced coordinated report edits -- the report-authoring class of V16qf violations. - new itis_sumo.report: flat curated surface (diagnostics, CV validation/calibration incl. T48np coverage/bootstrap/convergence fits, sensitivity/correlation, manual-UQ propagation); re-exports only, no logic - contract test pins __all__ completeness + object identity to the source modules (tests/test_report_surface.py) - SPEC: V49rp facade rule; T48np/T49np/T54rp land rows; T18ry unblocked (analyze_dataset real w/ this promotion); §I report row Co-authored-by: Qwen3.8 Flash Next NVFP4 (AWS vLLM) --- SPEC.md | 14 +++++++-- src/itis_sumo/report/__init__.py | 51 ++++++++++++++++++++++++++++++++ tests/test_report_surface.py | 37 +++++++++++++++++++++++ 3 files changed, 100 insertions(+), 2 deletions(-) create mode 100644 src/itis_sumo/report/__init__.py create mode 100644 tests/test_report_surface.py diff --git a/SPEC.md b/SPEC.md index d8e2d61..19ac2e8 100644 --- a/SPEC.md +++ b/SPEC.md @@ -60,7 +60,8 @@ Consumer contract (grill 2026-08-18): itis-sumo owns ALL surrogate machinery end - CLI: `itis-sumo validate` - python (consumer-facing) → `itis_sumo.api` = THE module flaskapi/mmux_vite imports; nothing else. Tabular in, typed results out, taxonomy errors (T22ax) - python (in-package/headless) → `itis_sumo.{core,config,data,sampling,evaluate,preprocess,utils}` public funcs; `core` now re-exports the E1 model store, `evaluate` now re-exports `export_sumo_model`/`import_sumo_model` (were reachable only via `funs_evaluate`, breaking V16qf) -- python (forthcoming) → `analyze_dataset` narrow diagnostics entrypoint (T18ry) — its output IS the override payload for the preprocessing defaults AND the auto-`domain` source, ⊥ a side feature +- python (diagnostics) → `analyze_dataset` narrow entrypoint (T18ry — unblocked by the T48np/T49np promotion, ships next tag) — its output IS the override payload for the preprocessing defaults AND the auto-`domain` source, ⊥ a side feature +- python (report authoring) → `itis_sumo.report` = THE flat surface automated report notebooks import (V49rp, T54rp): diagnostics ∧ CV validation/calibration ∧ sensitivity/correlation ∧ manual-UQ propagation; re-exports only, ⊥ logic lives there - artifacts: run_dir `{dakota_stdout.txt,dakota_stderr.txt,*.dat,*.sps/.alg}`; models dir `{id}.metadata.json` (E1) ## §V @@ -114,6 +115,12 @@ V45ls: ∀ value-producing `itis_sumo.api` entry point → `scale` ! be threaded V46rn: log scale ! compose w/ ∀ uncertainty distribution it is defined on: `uniform` → log-uniform draw ∧ `normal` → μ/σ parameterize the ln-space distribution (raw draw = `exp(N(μ,σ))`, lognormal in original units, positive by construction; the surrogate sees exactly `N(μ,σ)` in its training space — same "distribution describes what the model sees" contract as log-uniform); ⊥ a sampler path rejecting log⊗normal or drawing it linear; `constant` ⊥ composes w/ log (rejection stays) V47st: scale/spec surface input guards ! fail loud at the boundary, never mid-engine: `preprocessing` override maps ! be exact-cover at ∀ entry points (⊥ unknown-column override silently ignored — session path's rule extends to `compute_correlations`/`generate_lhs_samples`/`generate_grid_samples`); required scale maps ! be indexed by required key (⊥ per-variable `.get(..., "linear")` silent default); log-uniform bounds ! validated present ∧ ordered at the api boundary (→ `SumoInputError`, not `_run_engine`'s `SumoEngineError`); ∀ bounded spec (`DomainSpec`/`DistributionSpec`) ! reject `minimum ≥ maximum` at construction V48tr: Sobol' 2nd-order pairs ! come from the exact joint-pair estimator `S_ij = Var(E[Y|X_i,X_j])/V − S_i − S_j` (independent C stream ∧ U^ij/V^ij mixed designs, cost `n·(2+d+d(d−1))`, EXACT ∀ d); ⊥ the retired B26nc identity that inferred `O(d²)` pairs from `d` first/total gaps (underdetermined d≥4: pair sums collapse→0/negative); ∀ estimator quantity ! use pooled-A/B mean removal ∧ pooled variance ∧ centered products → translation-invariant under `f→f+c` ∧ identical to displayed scipy `saltelli_2010` point estimator; order masses `M1=ΣᵢSᵢ`/`M2=Σ_{i DatasetDiagnostics` narrow entrypoint (scale/distribution auto-selection + outlier surfacing, plain dataclasses, JSON-serializable via `dataclasses.asdict()`) as the sole flaskapi-facing dataset-diagnostics API, replacing any per-function flaskapi orchestration; BLOCKED — building blocks `select_variable_scale`/`auto_select_distributions` (+ a new raw-value outlier detector, reusing `_tukey_outlier_mask`'s IQR technique) currently exist only on confidential incubator branch `feat/nih-in-silico-example`, not `develop` — needs individual promotion first, same promotion rule as T15mn/E1|§C,V16qf,I +T18ry|x|design+implement `analyze_dataset(df, input_cols, output_cols, alpha=0.05, include_detail=False) -> DatasetDiagnostics` narrow entrypoint (scale/distribution auto-selection + outlier surfacing, plain dataclasses, JSON-serializable via `dataclasses.asdict()`) as the sole flaskapi-facing dataset-diagnostics API, replacing any per-function flaskapi orchestration; unblocked by T48np's promotion of `select_variable_scale`/`auto_select_distributions` + T49np's new raw-value outlier detector|§C,V16qf,I T19kp|.|headless notebook, post-alpha — add runnable notebook counterpart to docs getting-started flow once alpha docs publish is stable|T10le T20hm|.|CI workflow changes ! activate shared validation jobs and regression-check detector classification|V18rs T21vk|✓|VOCAB section (above) + `docs/reference/glossary.md` (nav-registered, strict build green); vocabulary enforced on the api layer by `tests/test_api_contract.py::TestPublicSurface`. Sweep of the OLDER modules' docstrings still outstanding → folded into T24cm|V19cn @@ -181,6 +188,9 @@ T50vb|x|`DomainSpec`/`DistributionSpec` ! reject `minimum ≥ maximum` at constr T51pq|x|port mmux_vite T31rb (PR#649): exact arbitrary-d 2nd-order Sobol' pair estimator + order masses — `_saltelli_abc` 3-stream design builder ∧ `_saltelli_pair_designs` (U/V mixed designs), `_sobol_algebra` (scipy `saltelli_2010` parity pinned by test, translation-invariant), `_sobol_joint_bootstrap` (shared-row CIs over ∀ index+mass), zero-variance → null masses; replaces the V9gh closed-form identity + the `d==1` special case + the runtime `scipy.stats.sobol_indices` call; public surface gains `OrderMasses` dataclass ∧ `SobolResult.order_contributions`; analytic benchmarks additive d=8 / pair-interaction d=5 / Ishigami / pair-quadratic d=10 drive the PRODUCTION helpers with 3·CI-scaled tolerances; B26nc degeneracy regression kept|V48tr,V9gh,V22rs,B20qt T52xx|x|land V26dd domain⊥distribution split (ahead of the handle transformation, per owner call 2026-09-30 — no defer): `evaluate_sobol`/`SumoSession.sobol` take optional `domains: Mapping[str, DomainSpec]` (unknown names → `SumoInputError`; boxes auto-inferred from observed bounds + detected scale; columns constant in the samples pinned ∧ echoed in `fixed`); Saltelli draws become uniform/log-uniform over the box (V44ls) — the `distributions` parameter, log⊗normal Sobol sampling ∧ the `mean±3σ` conflation die; engine `evaluate_sobol_indices` consumes a domain-`sampling` map (`{"minimum","maximum","log_scale"?}` | `{"value"}`); `SobolResult`: `distributions` → `domains` + `fixed` + `effective_config` (V21pf); NORMAL flip matrix ⊥ Sobol row (its flip rides DOMAIN boxes, uniform matrix); docs/examples/V&V aligned|V26dd,V44ls,V45ls,V21pf,V47st T53xx|x|explicit factor pinning for Sobol: `evaluate_sobol`/`SumoSession.sobol` accept `fixed: Mapping[str, float]` (unknown names rejected ∧ boxed∧fixed overlap rejected ∧ finite ∧ log-positive); pins echo in `SobolResult.fixed` beside auto-inferred constants — the domain-vocabulary freeze a consumer needs to hold a factor without distribution parameters (mmux FE Sobol-panel migration, prelude); test: pinned factor zero indices ∧ surrogate keeps its real spread|V26dd,V47st +T48np|x|port the NIH incubator's still-unpromoted diagnostics tooling onto this line (`CONFIDENTIAL/nih-merck-validation` → here, clean re-port not rebase): data `select_variable_scale`/`auto_select_distributions`; evaluate `compute_coverage`/`compute_cv_diagnostics`/`fit_convergence_exponential{,_asymptotic}` + bootstrap/`draws`/`diagnostics` redesign of `compute_cv_convergence` + `mean_signed_error`; NIH test files ride along; NIH's duplicate `propagate_uq_with_uncertainty` ⊥ ported — `propagate_manual_uq_with_uncertainty` made its superset instead (optional-`preprocessor` None-branch), callers adapt; api `CVAccuracyMetrics` absorbs the new `mean_signed_error` field; `ty`-clean (heterogeneous result dicts `dict[str, Any]` per house precedent)|V16wq,V17kb,V18wp,V19cz +T49np|x|new raw-value outlier detector `detect_raw_outliers(values, scale=, k=)` in `data/funs_dataset_diagnostics.py` (Tukey fence computed in the variable's selected scale, flag-only, NaN-safe) + wire `analyze_dataset` for real (scale/distribution via `auto_select_distributions`, outliers in each variable's own scale, `include_detail` candidates, JSON round-trip) — closes T18ry and the "raw-value outlier detection does not exist yet, anywhere" gap in BRANCH_CONSOLIDATION §3; flip-tests: geometric-vs-arithmetic over-flagging guard + ln↔log10 zero-flag-move|V46np,V16qf,T18ry +T54rp|x|new `itis_sumo.report` façade (V49rp): flat curated surface for automated report notebooks — api workflow verbs (cross_validate/evaluate_sobol/evaluate_correlations/…) ∧ `analyze_dataset` ∧ the T48np diagnostics (coverage, CV bundle, convergence bootstrap/fits) ∧ `propagate_manual_uq_with_uncertainty`; re-exports only, ⊥ logic; contract test pins completeness + object identity so report code never reaches into engine internals|V49rp,V16qf ## §B id|date|cause|fix diff --git a/src/itis_sumo/report/__init__.py b/src/itis_sumo/report/__init__.py new file mode 100644 index 0000000..84e695d --- /dev/null +++ b/src/itis_sumo/report/__init__.py @@ -0,0 +1,51 @@ +"""Report-authoring surface (V49rp): the curated entrypoints automated report +notebooks (marimo, T51np pattern; tissue-conductivity LHS instance, T52rv) +build on — the report-facing counterpart of `itis_sumo.api` (which serves the +flaskapi/mmux_vite request-response consumer, V16qf). + +One flat import (`from itis_sumo.report import ...`) for the whole authoring +arc: dataset diagnostics -> CV validation & calibration (coverage, Cohen's dz, +bootstrap convergence sweep, asymptotic error fits) -> sensitivity & correlation +-> manual UQ propagation. Re-exports only: no logic lives here, so engine +refactors stay invisible to report code as long as this surface holds. +""" + +from itis_sumo.api import ( + DEFAULT_SEED, + DistributionSpec, + compute_correlations, + cross_validate, + evaluate_along_axes, + evaluate_correlations, + evaluate_sobol, +) +from itis_sumo.data import analyze_dataset +from itis_sumo.evaluate import ( + compute_coverage, + compute_cv_accuracy_metrics, + compute_cv_convergence, + compute_cv_diagnostics, + compute_paired_ttest, + fit_convergence_exponential, + fit_convergence_exponential_asymptotic, +) +from itis_sumo.evaluate.funs_evaluate import propagate_manual_uq_with_uncertainty + +__all__ = [ + "DEFAULT_SEED", + "DistributionSpec", + "analyze_dataset", + "compute_correlations", + "compute_coverage", + "compute_cv_accuracy_metrics", + "compute_cv_convergence", + "compute_cv_diagnostics", + "compute_paired_ttest", + "cross_validate", + "evaluate_along_axes", + "evaluate_correlations", + "evaluate_sobol", + "fit_convergence_exponential", + "fit_convergence_exponential_asymptotic", + "propagate_manual_uq_with_uncertainty", +] diff --git a/tests/test_report_surface.py b/tests/test_report_surface.py new file mode 100644 index 0000000..3b77745 --- /dev/null +++ b/tests/test_report_surface.py @@ -0,0 +1,37 @@ +"""`itis_sumo.report` surface contract (V49rp, T54rp): the flat re-export set +automated report notebooks import against stays complete and identical to the +underlying api/data/evaluate objects (the facade must not fork behavior).""" + +import importlib + +import pytest + +from itis_sumo import report + +#: modules the facade is allowed to source from, in lookup order +_SOURCE_MODULES = ( + "itis_sumo.api", + "itis_sumo.data", + "itis_sumo.evaluate", + "itis_sumo.evaluate.funs_evaluate", +) + + +def test_report_all_importable_and_sorted() -> None: + names = report.__all__ + assert names == sorted(names), "__all__ must stay sorted (ruff RUF022 precedent)" + for name in names: + assert hasattr(report, name), f"itis_sumo.report missing {name}" + + +@pytest.mark.parametrize("name", report.__all__) +def test_reexports_are_identical_objects(name: str) -> None: + """Facade members must BE the underlying functions, not copies/wrappers — + so a fix in the engine lands in reports without a second promotion step.""" + facade_obj = getattr(report, name) + origins = [ + module + for module in _SOURCE_MODULES + if getattr(importlib.import_module(module), name, None) is facade_obj + ] + assert origins, f"{name} is not the same object as in any source module" From ca867ee81c0e3cd998e7c646fa1bcbff3d0d8e46 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Thu, 1 Oct 2026 14:14:25 +0000 Subject: [PATCH 4/4] chore: bump alpha version to 0.1.0a11 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 445b296..e5244ea 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "uv_build" [project] name = "itis-sumo" -version = "0.1.0a10" +version = "0.1.0a11" description = "Surrogate Modeling functionality for IT'IS Foundation / ZMT Modeling Intelligence suite: build, evaluate, cross-validate surrogates + UQ + sampling, headless or embedded" readme = "README.md" # T16mo rung 1 (1.5.9->1.5.11): 1.5.11 ships cp313 wheels, so the <3.13