PyApprox

https://sandialabs.github.io/pyapprox/
Admin
PyApprox is an open-source Python toolkit for surrogate modeling, sensitivity analysis, Bayesian...
Tokens:1,051,107
Snippets:6,675
Trust Score:9.7
Update:1 month ago
Show doc for...
Tokens:
Raw
### Shared setup: imports, benchmark, and optimizer

Source: https://sandialabs.github.io/pyapprox/notebooks/multifidelity_estimation_cookbook.ipynb

Sets up the common environment for all examples: imports, random seed, SLSQP optimizer, and the five-model polynomial benchmark. It prints the number of models and QoIs and the cost of each model.

```python
import numpy as np
import matplotlib.pyplot as plt
from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox_benchmarks.statest import (
    PolynomialEnsembleBenchmark,
)
from pyapprox.statest.statistics import MultiOutputMean
from pyapprox.statest.mc_estimator import MCEstimator
from pyapprox.statest import (
    MLMCEstimator, MFMCEstimator, GMFEstimator, GRDEstimator, GISEstimator,
)
from pyapprox.statest.acv import ACVAllocator, default_allocator_factory
from pyapprox.statest.acv.base import FittedACVEstimator
from pyapprox.statest.acv.search import ACVSearch
from pyapprox.statest.acv.strategies import (
    FixedRecursionStrategy, TreeDepthRecursionStrategy,
)
from pyapprox.statest.allocation import MCAllocator, CVAllocator
from pyapprox.statest.plotting import (
    plot_allocation, plot_estimator_variance_reductions,
)
from pyapprox.optimization.minimize.scipy.slsqp import ScipySLSQPOptimizer

bkd = NumpyBkd()
np.random.seed(42)

# SLSQP is more robust than the default trust-constr optimizer for ACV allocation
optimizer = ScipySLSQPOptimizer(maxiter=200)
allocator_factory = lambda est: default_allocator_factory(est, optimizer=optimizer)

benchmark = PolynomialEnsembleBenchmark(bkd, nmodels=5)
models    = benchmark.problem().models()        # f_0 (HF) ... f_4 (cheapest)
variable  = benchmark.problem().prior()
costs     = benchmark.problem().costs()
nqoi      = models[0].nqoi()
nmodels   = len(models)

costs_np = bkd.to_numpy(costs)
print(f"{nmodels} models, {nqoi} QoI(s)")
for a, c in enumerate(costs_np):
    print(f"  model {a}: cost = {c:.5f}")
```

--------------------------------

### Setup imports and backend

Source: https://sandialabs.github.io/pyapprox/boed_kl_nonlinear_usage.html

Imports required modules and sets up the backend and random seed for the Lotka-Volterra OED example.

```python
import numpy as np
import matplotlib.pyplot as plt
from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox_benchmarks.expdesign.lotka_volterra import LotkaVolterraPredictionOEDProblem
from pyapprox.expdesign.data import OEDDataGenerator
from pyapprox.expdesign import (
    GaussianOEDInnerLoopLikelihood,
    create_kl_oed_objective_from_data,
    RelaxedKLOEDSolver,
    RelaxedOEDConfig,
)

np.random.seed(1)
bkd = NumpyBkd()
```

--------------------------------

### Setup: Import PyApprox Hamiltonian surrogate and integrator classes

Source: https://sandialabs.github.io/pyapprox/learning_hamiltonians_usage.html

Imports the required PyApprox classes for Hamiltonian surrogates, derivative-matching fitters, implicit steppers, and the Numpy backend. This setup is needed before constructing any surrogate or running the tutorial's examples.

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox.probability import UniformMarginal
from pyapprox.surrogates.affine.basis import OrthonormalPolynomialBasis
from pyapprox.surrogates.affine.expansions import BasisExpansion
from pyapprox.surrogates.affine.indices import compute_hyperbolic_indices
from pyapprox.surrogates.affine.univariate import create_bases_1d
from pyapprox.surrogates.dynamical_systems import (
    BatchedBoundODEResidual,
    SnapshotDataset,
    FixedPoissonVariableHamiltonianSurrogate,
    FixedPoissonVariableHamiltonianDerivativeMatchingFitter,
    VariablePoissonFixedHamiltonianSurrogate,
    VariablePoissonFixedHamiltonianDerivativeMatchingFitter,
)
from pyapprox.ode.implicit_steppers.backward_euler import BackwardEulerAdjoint
from pyapprox.ode.implicit_steppers.crank_nicolson import CrankNicolsonAdjoint
from pyapprox.ode.implicit_steppers.implicit_midpoint import (
    ImplicitMidpointAdjoint,
)
from pyapprox.ode.implicit_steppers.integrator import TimeIntegrator
from pyapprox.util.rootfinding.newton import NewtonSolver

bkd = NumpyBkd()
```

--------------------------------

### Setup for MeshKLE and DataDrivenKLE

Source: https://sandialabs.github.io/pyapprox/kle_mesh_and_data_driven.html

Initializes shared setup: imports, random seed, backend, kernel, and mesh coordinates with trapezoid quadrature weights. Used as the foundation for both MeshKLE and DataDrivenKLE examples.

```python
import numpy as np
np.random.seed(42)
import matplotlib.pyplot as plt

from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox.surrogates.kernels import ExponentialKernel
from pyapprox.surrogates.kle import MeshKLE, DataDrivenKLE

bkd = NumpyBkd()
rng = np.random.default_rng(0)

# Shared setup
npts = 150
x = np.linspace(0, 2, npts)
mesh_coords = bkd.array(x[None, :])

# Trapezoid quadrature weights
dx = x[1] - x[0]
w = np.ones(npts) * dx; w[0] /= 2; w[-1] /= 2
quad_weights = bkd.array(w)
```

--------------------------------

### Shared Setup: Imports and PolynomialEnsembleBenchmark

Source: https://sandialabs.github.io/pyapprox/multifidelity_estimation_cookbook.html

Shared setup for all multi-fidelity examples. Imports PyApprox estimators, allocators, and plotting utilities; creates a five-model PolynomialEnsembleBenchmark; prints model costs. Requires the pyapprox_benchmarks package and a NumpyBkd backend.

```python
import numpy as np
import matplotlib.pyplot as plt
from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox_benchmarks.statest import (
    PolynomialEnsembleBenchmark,
)
from pyapprox.statest.statistics import MultiOutputMean
from pyapprox.statest.mc_estimator import MCEstimator
from pyapprox.statest import (
    MLMCEstimator, MFMCEstimator, GMFEstimator, GRDEstimator, GISEstimator,
)
from pyapprox.statest.acv import ACVAllocator, default_allocator_factory
from pyapprox.statest.acv.base import FittedACVEstimator
from pyapprox.statest.acv.search import ACVSearch
from pyapprox.statest.acv.strategies import (
    FixedRecursionStrategy, TreeDepthRecursionStrategy,
)
from pyapprox.statest.allocation import MCAllocator, CVAllocator
from pyapprox.statest.plotting import (
    plot_allocation, plot_estimator_variance_reductions,
)
from pyapprox.optim.minimize.scipy.slsqp import ScipySLSQPOptimizer

bkd = NumpyBkd()
np.random.seed(42)

# SLSQP is more robust than the default trust-constr optimizer for ACV allocation
optimizer = ScipySLSQPOptimizer(maxiter=200)
allocator_factory = lambda est: default_allocator_factory(est, optimizer=optimizer)

benchmark = PolynomialEnsembleBenchmark(bkd, nmodels=5)
models    = benchmark.problem().models()        # f_0 (HF) ... f_4 (cheapest)
variable  = benchmark.problem().prior()
costs     = benchmark.problem().costs()
nqoi      = models[0].nqoi()
nmodels   = len(models)

costs_np = bkd.to_numpy(costs)
print(f"{nmodels} models, {nqoi} QoI(s)")
for a, c in enumerate(costs_np):
    print(f"  model {a}: cost = {c:.5f}")
```

--------------------------------

### Initialize beam example environment

Source: https://sandialabs.github.io/pyapprox/notebooks/experimental_design_intro.ipynb

Set up the random seed, backend, and beam length for the cantilever beam example. This is a hidden setup cell that initializes the computational environment.

```python
#| echo: false
import numpy as np
np.random.seed(42)
import matplotlib.pyplot as plt
from pyapprox.util.backends.numpy import NumpyBkd

bkd = NumpyBkd()
rng = np.random.default_rng(42)
L = 100.0
```

--------------------------------

### Setup: imports, marginal, and test samples

Source: https://sandialabs.github.io/pyapprox/local_adaptive_vs_combination.html

Sets up the environment for the tutorial: imports PyApprox modules, defines a uniform marginal on [0,1], sets a maximum point budget of 1000, and generates 4000 random test samples in 2D. This setup is required before constructing the adaptive sparse grids and running the convergence experiments.

```python
import numpy as np
import matplotlib.pyplot as plt

from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox.surrogates.sparsegrids import (
    SingleFidelityAdaptiveSparseGridFitter,
    TensorProductSubspaceFactory,
)
from pyapprox.surrogates.sparsegrids.basis_factory import (
    ClenshawCurtisLagrangeFactory,
    PiecewiseFactory,
)
from pyapprox.surrogates.sparsegrids.basis.hierarchical_basis_1d import (
    HierarchicalBasis1D,
)
from pyapprox.surrogates.sparsegrids.hierarchical.hierarchical_fitter import (
    SingleFidelityHierarchicalFitter,
)
from pyapprox.surrogates.affine.indices import (
    ClenshawCurtisGrowthRule,
    MaxLevelCriteria,
)
from pyapprox.probability.univariate.uniform import UniformMarginal
from pyapprox.interface.functions.fromcallable.function import (
    FunctionFromCallable,
)

bkd = NumpyBkd()
np.random.seed(0)

marginal = UniformMarginal(0.0, 1.0, bkd)
max_pts = 1000

rng = np.random.default_rng(0)
test_samples = bkd.asarray(rng.uniform(0, 1, (2, 4000)))
```

--------------------------------

### Setup for Bin-Based Sensitivity Analysis with Ishigami Benchmark

Source: https://sandialabs.github.io/pyapprox/bin_based_sensitivity.html

Imports required libraries, initializes the random seed, loads the Ishigami benchmark, and prints the function and domain. This setup is required before running bin-based sensitivity analysis.

```python
import numpy as np
import matplotlib.pyplot as plt
from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox_benchmarks.sensitivity import IshigamiBenchmark
from pyapprox.sensitivity.variance_based import BinBasedSensitivityAnalysis

bkd = NumpyBkd()
np.random.seed(42)

# Load the Ishigami benchmark
benchmark = IshigamiBenchmark(bkd)
func = benchmark.problem().function()
prior = benchmark.problem().prior()

print(f"Ishigami function: f(x) = sin(x1) + 7*sin^2(x2) + 0.1*x3^4*sin(x1)")
print(f"Domain: x_i ~ U[-pi, pi], i = 1,2,3")
```

--------------------------------

### Setup: Import PyApprox modules and create NumpyBkd backend

Source: https://sandialabs.github.io/pyapprox/derivative_matching_usage.html

Imports the required PyApprox modules and creates a NumpyBkd backend instance for the tutorial. This setup is required before any of the subsequent steps.

```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox.probability import UniformMarginal
from pyapprox.surrogates.affine.basis import OrthonormalPolynomialBasis
from pyapprox.surrogates.affine.expansions import BasisExpansion
from pyapprox.surrogates.affine.expansions.fitters.least_squares import (
    LeastSquaresFitter,
)
from pyapprox.surrogates.affine.indices import compute_hyperbolic_indices
from pyapprox.surrogates.affine.univariate import create_bases_1d
from pyapprox.surrogates.dynamical_systems import (
    BatchedBoundODEResidual,
    SnapshotDataset,
)
from pyapprox.ode.implicit_steppers.backward_euler import BackwardEulerAdjoint
from pyapprox.ode.implicit_steppers.integrator import TimeIntegrator
from pyapprox.util.rootfinding.newton import NewtonSolver

bkd = NumpyBkd()
```

--------------------------------

### Set up beam model and posterior for DRAM tutorial

Source: https://sandialabs.github.io/pyapprox/notebooks/dram_mcmc.ipynb

Set up the cantilever beam model with 2 KLE terms, define the tip-deflection wrapper, generate synthetic observations, and construct the posterior. Also prepares a 1D beam with 1 KLE term for the Delayed Rejection illustration.

```python
#| echo: false
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox_benchmarks.pde.cantilever_beam import (
    build_cantilever_beam_1d,
)
from pyapprox.probability.likelihood import (
    DiagonalGaussianLogLikelihood,
    ModelBasedLogLikelihood,
)
from pyapprox.interface.functions.fromcallable.function import (
    FunctionFromCallable,
)
from pyapprox.inverse.posterior.log_unnormalized import (
    LogUnNormalizedPosterior,
)
from pyapprox.inverse.sampling.metropolis import (
    MetropolisHastingsSampler,
    AdaptiveMetropolisSampler,
)
from pyapprox.inverse.sampling.dram import (
    DelayedRejectionAdaptiveMetropolis,
)

bkd = NumpyBkd()

# 1D FEM beam with 2 KLE terms for spatially varying EI(x)
problem_2d = build_cantilever_beam_1d(bkd, num_kle_terms=2)
beam_2d = problem_2d.function()   # nvars=2, nqoi=3
prior_2d = problem_2d.prior()     # IndependentJoint of 2 x N(0,1)

# Tip-only wrapper (first QoI is tip deflection)
tip_model_2d = FunctionFromCallable(
    nqoi=1, nvars=2,
    fun=lambda xi: beam_2d(xi)[0:1, :],
    bkd=bkd,
)

# True KLE parameters and synthetic data
xi_true_2d = bkd.asarray([[0.5], [-0.3]])
xi1_true = float(xi_true_2d[0, 0])
xi2_true = float(xi_true_2d[1, 0])
delta_true_2d = float(tip_model_2d(xi_true_2d)[0, 0])
sigma_noise_2d = 0.1 * abs(delta_true_2d)
noise_lik_2d = DiagonalGaussianLogLikelihood(
    bkd.asarray([sigma_noise_2d**2]), bkd,
)
model_lik_2d = ModelBasedLogLikelihood(tip_model_2d, noise_lik_2d, bkd)

np.random.seed(7)
y_obs_2d = float(model_lik_2d.rvs(xi_true_2d))
noise_lik_2d.set_observations(bkd.asarray([[y_obs_2d]]))

# Posterior
posterior_2d = LogUnNormalizedPosterior(
    model_fn=tip_model_2d, likelihood=noise_lik_2d, prior=prior_2d, bkd=bkd,
)

# 1D FEM beam with 1 KLE term (for the DR illustration)
problem_1d = build_cantilever_beam_1d(bkd, num_kle_terms=1)
beam_1d = problem_1d.function()   # nvars=1, nqoi=3
prior_1d = problem_1d.prior()

tip_model_1d = FunctionFromCallable(
    nqoi=1, nvars=1,
    fun=lambda xi: beam_1d(xi)[0:1, :],
    bkd=bkd,
)

xi_true_1d = 0.5

```

--------------------------------

### Setup: import NumPy, Matplotlib, and NumpyBkd

Source: https://sandialabs.github.io/pyapprox/multifidelity_sparse_grids.html

Sets up the NumPy backend and random seed for the multi-fidelity sparse grid tutorial. Required before building the model hierarchy and discrepancy approximations.

```python
import numpy as np
import matplotlib.pyplot as plt
from pyapprox.util.backends.numpy import NumpyBkd

bkd = NumpyBkd()
np.random.seed(42)
```

--------------------------------

### Import pyapprox piecewise interpolation classes

Source: https://sandialabs.github.io/pyapprox/notebooks/piecewise_interpolation.ipynb

Imports the required pyapprox classes and sets up the NumPy backend. This setup is needed for the subsequent piecewise interpolation examples.

```python
import numpy as np
np.random.seed(42)
import matplotlib.pyplot as plt

from pyapprox.util.backends.numpy import NumpyBkd

from pyapprox.surrogates.affine.univariate import (
    LagrangeBasis1D,
    PiecewiseLinear,
    PiecewiseQuadratic,
    PiecewiseCubic,
)
from pyapprox.surrogates.affine.univariate.globalpoly import (
    ClenshawCurtisQuadratureRule,
)

bkd = NumpyBkd()
```

--------------------------------

### Setup: imports and shared test data

Source: https://sandialabs.github.io/pyapprox/local_basis_sparse_grids.html

Imports required PyApprox modules and sets up a shared backend, random seed, marginal distribution, and test samples for the sparse grid examples.

```python
import numpy as np
import matplotlib.pyplot as plt

from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox.surrogates.sparsegrids import (
    IsotropicSparseGridFitter,
    SingleFidelityAdaptiveSparseGridFitter,
    TensorProductSubspaceFactory,
)
from pyapprox.surrogates.sparsegrids.basis_factory import (
    ClenshawCurtisLagrangeFactory,
    PiecewiseFactory,
)
from pyapprox.surrogates.affine.indices import (
    ClenshawCurtisGrowthRule,
    CubicNestedGrowthRule,
    MaxLevelCriteria,
)
from pyapprox.probability.univariate.uniform import UniformMarginal

bkd = NumpyBkd()
np.random.seed(42)

marginal = UniformMarginal(-1.0, 1.0, bkd)

# Shared test sets
rng = np.random.default_rng(0)
test_samples = bkd.asarray(rng.uniform(-1, 1, (2, 4000)))
```

--------------------------------

### Setup imports and backend for Leja sequence tutorial

Source: https://sandialabs.github.io/pyapprox/leja_sequences.html

Imports the PyApprox modules needed for Leja sequence and quadrature examples, sets a random seed for reproducibility, and creates a NumPy backend instance.

```python
import numpy as np
np.random.seed(42)
import matplotlib.pyplot as plt

from pyapprox.util.backends.numpy import NumpyBkd

from pyapprox.surrogates.affine.univariate import (
    JacobiPolynomial1D,
    LegendrePolynomial1D,
    GaussQuadratureRule,
    LagrangeBasis1D,
)
from pyapprox.surrogates.affine.univariate.globalpoly import (
    ClenshawCurtisQuadratureRule,
)
from pyapprox.surrogates.affine.leja import (
    LejaSequence1D,
    ChristoffelWeighting,
)
from pyapprox.surrogates.sparsegrids import LejaLagrangeFactory

bkd = NumpyBkd()
```

--------------------------------

### Setup: Import PyApprox modules and initialize backend

Source: https://sandialabs.github.io/pyapprox/adaptive_sparse_grids.html

Imports required PyApprox and plotting modules, sets the backend to NumpyBkd, and seeds the random number generator for reproducibility. This setup is required before using the adaptive sparse grid fitter.

```python
import numpy as np
import matplotlib.pyplot as plt

from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox.surrogates.sparsegrids import (
    IsotropicSparseGridFitter,
    SingleFidelityAdaptiveSparseGridFitter,
    TensorProductSubspaceFactory,
)
from pyapprox.surrogates.sparsegrids.basis_factory import (
    ClenshawCurtisLagrangeFactory,
)
from pyapprox.surrogates.affine.indices import (
    ClenshawCurtisGrowthRule,
    MaxLevelCriteria,
    Max1DLevelsCriteria,
    CompositeCriteria,
)
from pyapprox.probability.univariate.uniform import UniformMarginal
from pyapprox.interface.functions.plot.plot2d_rectangular import (
    Plotter2DRectangularDomain,
)
from pyapprox.interface.functions.fromcallable.function import (
    FunctionFromCallable,
)
from pyapprox_tutorials.figures._quadrature import (
    plot_additive_function,
    plot_index_set,
    plot_points_compare,
    plot_adaptive_vs_iso,
    plot_4d_adaptive,
)

bkd = NumpyBkd()
np.random.seed(42)
```

--------------------------------

### Evaluate forward model and extract latent noise samples

Source: https://sandialabs.github.io/pyapprox/boed_kl_nonlinear_usage.html

Evaluates the observation model on the outer and inner parameter samples and extracts the latent noise samples from the outer loop for the reparameterization trick. Assumes outerloop_samples and innerloop_samples were produced by the OEDDataGenerator setup in the previous example.

```python
nparams = prior.nvars()

# Extract prior parameter samples from outer loop
theta_outer = outerloop_samples[:nparams, :]

# Evaluate the observation model at outer and inner parameter samples
outer_shapes = obs_model(theta_outer)     # shape (nobs, nouter)
inner_shapes = obs_model(innerloop_samples)  # shape (nobs, ninner)

# Latent noise samples for the reparameterization trick
latent_samples = outerloop_samples[nparams:, :]  # shape (nobs, nouter)
```

--------------------------------

### Set up MultiOutputEnsembleBenchmark model with three QoIs

Source: https://sandialabs.github.io/pyapprox/estimator_accuracy_mse.html

Sets up the lightweight analytical benchmark model with three QoIs used throughout the tutorial. The model is obtained from the MultiOutputEnsembleBenchmark with psd=True, and the highest-fidelity model is selected. The number of QoIs K is retrieved via model.nqoi().

```python
import numpy as np
import matplotlib.pyplot as plt
from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox_benchmarks.statest import (
    MultiOutputEnsembleBenchmark,
)

bkd = NumpyBkd()
bm = MultiOutputEnsembleBenchmark(bkd, psd=True)
model = bm.problem().models()[0]       # highest-fidelity model
variable = bm.problem().prior()
K = model.nqoi()              # 3 QoIs
```

--------------------------------

### Install PyApprox packages if missing

Source: https://sandialabs.github.io/pyapprox/notebooks/adaptive_gp_sampling.ipynb

Installs the required PyApprox packages from GitHub if they are not already available. Use this in Google Colab or when running locally without the packages installed.

```python
try:
    import pyapprox
    import pyapprox_benchmarks
    import pyapprox_tutorials
except ImportError:
    !pip install -q \
        "pyapprox[runtime-extras] @ git+https://github.com/sandialabs/pyapprox.git#subdirectory=packages/pyapprox" \
        "pyapprox-benchmarks @ git+https://github.com/sandialabs/pyapprox.git#subdirectory=packages/pyapprox-benchmarks" \
        "pyapprox-tutorials @ git+https://github.com/sandialabs/pyapprox.git#subdirectory=packages/pyapprox-tutorials"
```

--------------------------------

### Setup LotkaVolterraPredictionOEDProblem and extract obs_map and prior

Source: https://sandialabs.github.io/pyapprox/boed_kl_nonlinear_usage.html

Setup for the Lotka-Volterra OED problem: imports, random seed, and backend, then instantiate the problem and extract the observation model and prior.

```python
import numpy as np
import matplotlib.pyplot as plt
from pyapprox.util.backends.numpy import NumpyBkd
from pyapprox_benchmarks.expdesign.lotka_volterra import LotkaVolterraPredictionOEDProblem
from pyapprox.expdesign.data import OEDDataGenerator
from pyapprox.expdesign import (
    GaussianOEDInnerLoopLikelihood,
    create_kl_oed_objective_from_data,
    RelaxedKLOEDSolver,
    RelaxedOEDConfig,
)

np.random.seed(1)
bkd = NumpyBkd()
```

```python
oed = LotkaVolterraPredictionOEDProblem(bkd, final_time=50.0)
obs_model = oed.obs_map()   # maps theta -> (nobs_times * 2,)
prior     = oed.inference_problem().prior()  # IndependentJoint over theta
```
Tokens:
Raw