Adaptive sampling algorithms from structural optimization can improve the efficiency of machine learning surrogate models for simulation problems.

Adaptive sampling algorithms originally developed for structural/mechanical optimization (e.g., space-filling designs, expected improvement, variance-based sequential sampling) can reduce the number of expensive simulation evaluations required to train a machine learning surrogate model to a target accuracy threshold by ≥30% compared to passive (random or grid) sampling strategies, across at least two distinct simulation problem classes (e.g., computational fluid dynamics, finite element analysis) and at least two surrogate model architectures (e.g., Gaussian Process, neural network). Formally: given a fixed accuracy budget ε (e.g., RMSE ≤ 5% of output range), adaptive sampling reaches ε with N_adaptive evaluations such that N_adaptive ≤ 0.70 × N_passive, with statistical significance p < 0.05 over ≥10 independent trials per configuration.

MINIMUM VIABLE TEST (MVT) — 4-benchmark, 3-algorithm, 2-surrogate factorial design:

BENCHMARK PROBLEMS (4 total): B1: Branin-Hoo function (d=2, analytical, cheap proxy for smooth structural response) B2: Hartmann-6 function (d=6, multimodal, proxy for nonlinear FEA response) B3: OpenFOAM lid-driven cavity CFD (d=4: Re, mesh density, boundary conditions × 2) — real simulation, ~15 min/run B4: FEniCS linear elasticity beam deflection (d=5: geometry, material, load parameters) — real simulation, ~8 min/run

ADAPTIVE ALGORITHMS (3, all from structural optimization literature): A1: Bayesian Optimization with Expected Improvement (EI) + GP surrogate (Mockus 1978, widely used in structural opt.) A2: LOLA-Voronoi sequential space-filling (Crombecq et al. 2011, developed for structural simulation surrogates) A3: FLOLA-Voronoi with gradient-enhanced exploration (van der Herten et al. 2015)

BASELINE: A0: Latin Hypercube Sampling (LHS) — standard passive baseline

SURROGATE MODELS (2): S1: Gaussian Process (Matérn 5/2 kernel, GPy library) S2: Neural Network (3-layer MLP, 64-128-64 units, ReLU, Adam optimizer)

FACTORIAL DESIGN: 4 benchmarks × 4 algorithms (3 adaptive + 1 baseline) × 2 surrogates × 10 independent trials = 320 experimental runs

EVALUATION METRIC: Number of simulation evaluations N* to reach RMSE ≤ 5% of output range on held-out uniform test grid (N_test = 500 for analytical, 200 for real simulations). Secondary: RMSE at fixed N ∈ {50, 100, 200, 500}.

# ============================================================
# ADAPTIVE SAMPLING EVP — IMPLEMENTATION SKETCH
# ============================================================
# Dependencies: botorch, gpytorch, smt, fenics, PyFOAM, scipy, numpy, h5py

import numpy as np
from scipy.stats.qmc import LatinHypercube
from botorch.models import SingleTaskGP
from botorch.acquisition import ExpectedImprovement
from botorch.optim import optimize_acqf
from smt.sampling_methods import LHS as SMT_LHS
from smt.surrogate_models import KRG  # for LOLA-Voronoi integration
import h5py, time

# ---- CONFIGURATION ----
BENCHMARKS = {
    'branin': {'d': 2, 'bounds': [(-5,10),(0,15)], 'func': branin_fn, 'sim': False},
    'hartmann6': {'d': 6, 'bounds': [(0,1)]*6, 'func': hartmann6_fn, 'sim': False},
    'openfoam_cavity': {'d': 4, 'bounds': [(100,5000),(20,80),(0.1,1.0),(0.1,1.0)],
                        'func': openfoam_wrapper, 'sim': True, 'cost_min': 900},
    'fenics_beam': {'d': 5, 'bounds': [(0.5,2.0),(0.05,0.2),(70e9,210e9),(1e3,50e3),(0.25,0.45)],
                    'func': fenics_wrapper, 'sim': True, 'cost_min': 480},
}
ALGORITHMS = ['lhs_baseline', 'ei_bo', 'lola_voronoi', 'flola_voronoi']
SURROGATES = ['gp', 'mlp']
N_TRIALS = 10
N_INIT = 10
N_MAX = 500
EPSILON_FRAC = 0.05  # 5% of output range

# ---- SIMULATION WRAPPERS ----
def openfoam_wrapper(x: np.ndarray) -> float:
    """
    x: [Re, mesh_size, inlet_u, inlet_v]
    Returns: L2 norm of velocity field (scalar QoI)
    Calls PyFOAM to modify blockMeshDict, run icoFoam, extract postProcessing data
    """
    Re, mesh, u_in, v_in = x
    config = write_openfoam_config(Re=Re, mesh=int(mesh), u=(u_in, v_in))
    result = run_openfoam(config, timeout=1200)  # 20 min timeout
    if result.converged:
        return extract_velocity_norm(result.path)
    else:
        raise SimulationError(f"OpenFOAM non-convergence at x={x}")

def fenics_wrapper(x: np.ndarray) -> float:
    """
    x: [length, width, E, load, nu]
    Returns: max deflection (m)
    """
    L, W, E, F, nu = x
    mesh = create_beam_mesh(L, W)
    u = solve_elasticity(mesh, E=E, nu=nu, load=F)
    return float(u.vector().max())

# ---- SURROGATE TRAINING ----
def train_surrogate(X: np.ndarray, y: np.ndarray, stype: str):
    if stype == 'gp':
        from gpytorch.kernels import MaternKernel, ScaleKernel
        model = SingleTaskGP(torch.tensor(X), torch.tensor(y).unsqueeze(-1))
        mll = ExactMarginalLogLikelihood(model.likelihood, model)
        fit_gpytorch_mll(mll)
        return model
    elif stype == 'mlp':
        from sklearn.neural_network import MLPRegressor
        model = MLPRegressor(hidden_layer_sizes=(64,128,64), activation='relu',
                             max_iter=2000, random_state=42)
        model.fit(X, y)
        return model

def evaluate_surrogate(model, X_test: np.ndarray, y_test: np.ndarray,
                       stype: str, y_range: float) -> float:
    if stype == 'gp':
        with torch.no_grad():
            pred = model(torch.tensor(X_test)).mean.numpy()
    else:
        pred = model.predict(X_test)
    rmse = np.sqrt(np.mean((pred - y_test)**2))
    return rmse / y_range  # normalized RMSE

# ---- ADAPTIVE SAMPLING LOOPS ----
def run_ei_bo(benchmark: dict, surrogate_type: str, trial: int) -> dict:
    """Bayesian Optimization with Expected Improvement"""
    bounds = torch.tensor(benchmark['bounds']).T.float()
    X = lhs_init(benchmark['d'], N_INIT, benchmark['bounds'])
    y = np.array([benchmark['func'](x) for x in X])
    results = {'N': [], 'rmse': [], 'overhead_frac': []}

    for n in range(N_INIT, N_MAX):
        t0 = time.time()
        model = train_surrogate(X, y, 'gp')  # GP always for acquisition
        best_f = torch.tensor(y.min())
        EI = ExpectedImprovement(model, best_f=best_f)
        x_new, _ = optimize_acqf(EI, bounds=bounds, q=1, num_restarts=10,
                                  raw_samples=512)
        t_overhead = time.time() - t0

        t_sim = time.time()
        y_new = benchmark['func'](x_new.numpy().flatten())
        t_sim_elapsed = time.time() - t_sim

        X = np.vstack([X, x_new.numpy()])
        y = np.append(y, y_new)

        surrogate = train_surrogate(X, y, surrogate_type)
        rmse_norm = evaluate_surrogate(surrogate, X_test, y_test,
                                       surrogate_type, y_range)
        results['N'].append(n+1)
        results['rmse'].append(rmse_norm)
        results['overhead_frac'].append(t_overhead / (t_overhead + t_sim_elapsed))

        if rmse_norm <= EPSILON_FRAC:
            results['N_star'] = n + 1
            break
    return results

def run_lola_voronoi(benchmark: dict, surrogate_type: str, trial: int) -> dict:
    """LOLA-Voronoi sequential space-filling (via SMT)"""
    from smt.sampling_methods import LolaVoronoi
    sampling = LolaVoronoi(xlimits=np.array(benchmark['bounds']))
    X = sampling(N_INIT)
    y = np.array([benchmark['func'](x) for x in X])
    results = {'N': [], 'rmse': [], 'overhead_frac': []}

    for n in range(N_INIT, N_MAX):
        t0 = time.time()
        sampling.update(X, y)
        x_new = sampling(1)
        t_overhead = time.time() - t0

        t_sim = time.time()
        y_new = benchmark['func'](x_new.flatten())
        t_sim_elapsed = time.time() - t_sim

        X = np.vstack([X, x_new])
        y = np.append(y, y_new)

        surrogate = train_surrogate(X, y, surrogate_type)
        rmse_norm = evaluate_surrogate(surrogate, X_test, y_test,
                                       surrogate_type, y_range)
        results['N'].append(n+1)
        results['rmse'].append(rmse_norm)
        results['overhead_frac'].append(t_overhead / (t_overhead + t_sim_elapsed))

        if rmse_norm <= EPSILON_FRAC:
            results['N_star'] = n + 1
            break
    return results

# ---- MAIN EXPERIMENT LOOP ----
def run_full_experiment():
    all_results = {}
    for bname, benchmark in BENCHMARKS.items():
        # Generate held-out test grid
        X_test, y_test = generate_test_grid(benchmark, N=500 if not benchmark['sim'] else 200)
        y_range = y_test.max() - y_test.min()

        for algo in ALGORITHMS:
            for stype in SURROGATES:
                trial_results = []
                for trial in range(N_TRIALS):
                    np.random.seed(trial * 42)
                    if algo == 'lhs_baseline':
                        res = run_lhs_baseline(benchmark, stype, trial)
                    elif algo == 'ei_bo':
                        res = run_ei_bo(benchmark, stype, trial)
                    elif algo == 'lola_voronoi':
                        res = run_lola_voronoi(benchmark, stype, trial)
                    elif algo == 'flola_voronoi':
                        res = run_flola_voronoi(benchmark, stype, trial)
                    trial_results.append(res)

                key = f"{bname}_{algo}_{stype}"
                all_results[key] = aggregate_trial_results(trial_results)

    # Statistical analysis
    efficiency_ratios = compute_efficiency_ratios(all_results)
    p_values = run_wilcoxon_tests(all_results, correction='bonferroni')
    overhead_stats = compute_overhead_stats(all_results)

    # Save results
    with h5py.File('evp_results.h5', 'w') as f:
        save_results_hdf5(f, all_results, efficiency_ratios, p_values)

    generate_figures(efficiency_ratios, p_values, overhead_stats)
    return all_results, efficiency_ratios, p_values

# ---- ABORT CHECKPOINT LOGIC ----
def check_abort_conditions(results_so_far: dict, checkpoint_n: int) -> bool:
    """Called after each benchmark completes"""
    if checkpoint_n >= 2:  # After first 2 benchmarks
        mean_R = np.mean([r.get('efficiency_ratio', 1.0)
                          for r in results_so_far.values()])
        if mean_R < 0.95:
            print(f"ABORT: Mean efficiency ratio {mean_R:.3f} < 0.95 after {checkpoint_n} benchmarks")
            return True
        convergence_failures = sum(1 for r in results_so_far.values()
                                   if r.get('sim_failure_rate', 0) > 0.20)
        if convergence_failures > 0:
            print(f"ABORT: Simulation failure rate > 20% in {convergence_failures} conditions")
            return True
    return False

if __name__ == '__main__':
    results, ratios, pvals = run_full_experiment()
    print_summary_table(ratios, pvals)