Getting Started with MCPost¶

# Import required libraries
import numpy as np
import matplotlib.pyplot as plt
from mcpost import gsa_pipeline, monte_carlo_integral, qmc_integral_auto

# Quieten sklearn's per-restart optimiser chatter. MCPost's own warnings --
# surrogate collapse, Sobol range violations, length scales on a bound -- are
# NOT suppressed: those are the ones you want to see.
import warnings
from sklearn.exceptions import ConvergenceWarning
warnings.filterwarnings("ignore", category=ConvergenceWarning)

# Set random seed for reproducibility
np.random.seed(42)

print("MCPost successfully imported!")

# Define a simple test function
def polynomial_function(X):
    """
    Simple polynomial: f(x1, x2, x3) = x1^2 + 2*x2 + 0.1*x3
    
    We expect x2 to be most influential, x1 moderately influential, 
    and x3 to have minimal influence.
    """
    x1, x2, x3 = X[:, 0], X[:, 1], X[:, 2]
    return x1**2 + 2*x2 + 0.1*x3

# Generate parameter samples
n_samples = 1000
X = np.random.uniform(-1, 1, (n_samples, 3))  # 3 parameters in [-1, 1]

# Evaluate function
y = polynomial_function(X)
Y = y.reshape(-1, 1)  # GSA expects 2D array

print(f"Generated {n_samples} samples")
print(f"Parameter ranges: {X.min(axis=0)} to {X.max(axis=0)}")
print(f"Output range: [{y.min():.3f}, {y.max():.3f}]")

# Run GSA analysis
param_names = ["x1", "x2", "x3"]
feature_names = ["polynomial"]

print("Running GSA analysis...")
results = gsa_pipeline(
    X, Y,
    param_names=param_names,
    feature_names=feature_names,
    scaler="minmax",
    enable_sobol=True,
    enable_gp=True,
    enable_perm=True,
    make_pdp=False,  # Skip PDPs for this simple example
    N_sobol=2048
)

# Display results
sensitivity_table = results["results"]["polynomial"]["table"]
print("\nSensitivity Analysis Results:")
print(sensitivity_table)

# Visualize the results
from mcpost import plot_sensitivity_metrics

fig, ax = plot_sensitivity_metrics(
    sensitivity_table,
    title="Polynomial Function - Sensitivity Analysis"
)
plt.show()

print("\nInterpretation:")
print("- x2 shows highest sensitivity (expected - linear coefficient of 2)")
print("- x1 shows moderate sensitivity (quadratic term)")
print("- x3 shows low sensitivity (small coefficient of 0.1)")

# Extract individual metrics for interpretation
table = sensitivity_table

print("GSA Metrics Explanation:")
print("=" * 25)

for param in param_names:
    print(f"\nParameter: {param}")
    print(f"  MI (Mutual Information): {table.loc[param, 'MI']:.3f}")
    print("    -> Statistical dependence, any shape (higher = more influential)")

    print(f"  dCor (Distance Correlation): {table.loc[param, 'dCor']:.3f}")
    print("    -> Nonlinear dependence, always in [0, 1]")

    print(f"  PermMean (Permutation Importance): {table.loc[param, 'PermMean']:.3f}")
    print("    -> Drop in model R^2 when this column is shuffled")

    print(f"  S1 (First-order Sobol): {table.loc[param, 'S1']:.3f}")
    print("    -> Variance explained by this parameter ALONE")

    print(f"  ST (Total-effect Sobol): {table.loc[param, 'ST']:.3f}")
    print("    -> Variance explained by it plus all its interactions")

# Sanity checks worth making every time. MCPost warns about these itself, but
# it costs nothing to assert them where you use the numbers.
s1, st = table["S1"], table["ST"]
in_range = bool(((s1 >= -0.05) & (s1 <= st + 0.05) & (st <= 1.05)).all())
print("\nSanity checks:")
print(f"  0 <= S1 <= ST <= 1 for every parameter : {in_range}")
print(f"  sum(S1) = {s1.sum():.3f}")
print("     Must not exceed 1. A value near 0 usually means the surrogate")
print("     carries no signal, NOT that the parameters are purely interacting.")

# This function is simple enough to have closed-form Sobol indices, so we can
# check the answer rather than just eyeballing the ranking. For
# f = x1^2 + 2*x2 + 0.1*x3 with x ~ U(-1, 1)^3 the terms are independent, so
#     V1 = Var(x^2)  = E[x^4] - E[x^2]^2 = 1/5 - 1/9
#     V2 = Var(2x)   = 4/3
#     V3 = Var(0.1x) = 0.01/3
# and S1_i = V_i / (V1 + V2 + V3). No interactions, so ST == S1 and sum(S1) = 1.
v1 = 1 / 5 - (1 / 3) ** 2
v2 = 4 * (1 / 3)
v3 = 0.01 * (1 / 3)
V = v1 + v2 + v3
exact = dict(zip(param_names, [v1 / V, v2 / V, v3 / V]))

print("\nAgainst the closed-form indices:")
print(f"{'Param':<7}{'S1 (GP)':>10}{'S1 (exact)':>12}{'abs err':>10}")
print("-" * 39)
for param in param_names:
    got = table.loc[param, "S1"]
    print(f"{param:<7}{got:>10.4f}{exact[param]:>12.4f}{abs(got - exact[param]):>10.4f}")
print("\nAgreeing with theory is a much stronger check than the ranking looking")
print("plausible. Whenever your test function has a known answer, compare to it.")

# Example: Integration over a custom parameter space
def custom_target(theta):
    """Custom target distribution (mixture of Gaussians)"""
    comp1 = 0.6 * np.exp(-0.5 * np.sum((theta - 1) ** 2, axis=1))
    comp2 = 0.4 * np.exp(-0.5 * np.sum((theta + 1) ** 2, axis=1))
    return (comp1 + comp2) / (2 * np.pi)

def complex_integrand(theta):
    """More complex integrand"""
    return np.exp(theta[:, 0]) * np.cos(theta[:, 1]) * theta[:, 0] ** 2

BOX = [(-3, 3), (-3, 3)]
result = qmc_integral_auto(
    N_samples=8192,
    N_params=2,
    data_func=complex_integrand,
    p_target=custom_target,
    bounds=BOX,
    method="sobol",   # or "halton"
    seed=0,           # omit for fresh scrambling on every call
)
qmc_estimate = float(np.ravel(result["integral"])[0])

# Cross-check against a dense tensor-product grid. This is a completely
# different quadrature rule, so agreement is evidence -- unlike comparing two
# QMC runs, which would share any bug the estimator has.
g = np.linspace(-3, 3, 1201)
gx, gy = np.meshgrid(g, g, indexing="ij")
pts = np.column_stack([gx.ravel(), gy.ravel()])
integrand_vals = complex_integrand(pts) * custom_target(pts)
h = g[1] - g[0]
grid_estimate = integrand_vals.reshape(gx.shape).sum() * h * h

print(f"QMC estimate            : {qmc_estimate:.6f}")
print(f"Grid quadrature (1201^2): {grid_estimate:.6f}")
print(f"Relative difference     : {abs(qmc_estimate - grid_estimate) / abs(grid_estimate):.2e}")
print(f"Effective sample size   : {result['effective_sample_size']:.0f}")

print("\nNote the integral is over the BOX, not all of R^2. exp(x)*x^2 grows fast")
print("enough that truncating at x = 3 matters here -- widen `bounds` and the")
print("answer moves. That is a property of the problem, not of the estimator.")

# Define integration problem: E[x^2] where x ~ N(0,1)
# Analytical solution: 1.0

def integrand(theta):
    """Function to integrate: f(x) = x^2"""
    return theta[:, 0]**2

def target_pdf(theta):
    """Standard normal PDF"""
    return np.exp(-0.5 * theta[:, 0]**2) / np.sqrt(2 * np.pi)

print("Integration Problem: E[X^2] where X ~ N(0,1)")
print("Analytical solution: 1.0")
print()

# Method 1: Standard Monte Carlo
n_samples = 5000
theta_samples = np.random.normal(0, 1, (n_samples, 1))
f_values = integrand(theta_samples)   # 1-D is fine: one integrand

mc_result = monte_carlo_integral(theta_samples, f_values, target_pdf)

print("Standard Monte Carlo:")
print(f"  Integral estimate    : {float(np.ravel(mc_result['integral'])[0]):.6f}")
print(f"  Uncertainty          : {float(np.ravel(mc_result['uncertainty'])[0]):.6f}")
print(f"  Effective sample size: {mc_result['effective_sample_size']:.0f}")
print(f"  Error vs analytic 1.0: {abs(float(np.ravel(mc_result['integral'])[0]) - 1.0):.6f}")

print("\nThe uncertainty is scaled by the EFFECTIVE sample size, not N. When the")
print("samples are drawn from something far from the target, ESS falls well below")
print("N and the error bar widens to match -- which is the whole point of quoting")
print("one. Compare ESS against N above to see how much of your sample is")
print("actually doing work.")

# Method 2: Quasi-Monte Carlo
#
# NOTE: before v0.2.0 this cell was marked "seems problematic". It was.
# qmc_integral_auto normalised its importance weights to sum to 1 -- turning the
# sum into a weighted mean -- and then multiplied by the box volume anyway, so
# every result came back exactly `volume` times too large. On the [-4, 4] box
# below that is a factor of 8: it printed 8.0 where the answer is 1.0.
qmc_result = qmc_integral_auto(
    N_samples=n_samples,
    N_params=1,
    data_func=integrand,
    p_target=target_pdf,
    bounds=[(-4, 4)],   # Integration bounds
    method="sobol",
    seed=0,
)

qmc_value = float(np.ravel(qmc_result["integral"])[0])
print("\nQuasi-Monte Carlo (Sobol):")
print(f"  Integral estimate    : {qmc_value:.6f}")
print(f"  Effective sample size: {qmc_result['effective_sample_size']:.0f}")
print(f"  Error vs analytic 1.0: {abs(qmc_value - 1.0):.6f}")

# Compare convergence rates
sample_sizes = [100, 500, 1000, 2000, 5000]
mc_errors = []
qmc_errors = []

print("\nConvergence Comparison (target: E[X^2] = 1 for X ~ N(0,1)):")
print(f"{'N_samples':<10} {'MC Error':<14} {'QMC Error':<14}")
print("-" * 40)

for n in sample_sizes:
    # Monte Carlo
    theta = np.random.normal(0, 1, (n, 1))
    mc_res = monte_carlo_integral(theta, integrand(theta), target_pdf)
    mc_error = abs(float(np.ravel(mc_res["integral"])[0]) - 1.0)
    mc_errors.append(mc_error)

    # Quasi-Monte Carlo, seeded so this table is reproducible
    qmc_res = qmc_integral_auto(
        N_samples=n, N_params=1, data_func=integrand,
        p_target=target_pdf, bounds=[(-4, 4)], seed=n,
    )
    qmc_error = abs(float(np.ravel(qmc_res["integral"])[0]) - 1.0)
    qmc_errors.append(qmc_error)

    print(f"{n:<10} {mc_error:<14.6f} {qmc_error:<14.6f}")

print("\nQMC converges faster than MC for smooth integrands: roughly O(1/N)")
print("against O(1/sqrt(N)). Both are compared against the exact value 1.0, not")
print("against each other -- two estimators agreeing tells you nothing if they")
print("share a bug.")

# Define a more complex model
def engineering_model(X):
    """
    Simplified engineering model: stress in a beam
    
    Parameters:
    - x1: Load [kN]
    - x2: Length [m] 
    - x3: Width [cm]
    - x4: Height [cm]
    
    Output: Maximum stress [MPa]
    """
    load = X[:, 0]      # Load in kN
    length = X[:, 1]    # Length in m
    width = X[:, 2]     # Width in cm
    height = X[:, 3]    # Height in cm
    
    # Moment of inertia (simplified)
    I = width * height**3 / 12
    
    # Maximum bending moment
    M_max = load * length**2 / 8  # Simply supported beam, uniform load
    
    # Maximum stress (bending)
    stress = M_max * (height/2) / I * 1000  # Convert to MPa
    
    return stress

# Parameter ranges (realistic engineering values)
param_ranges = {
    'Load': [10, 50],      # kN
    'Length': [2, 8],      # m
    'Width': [10, 30],     # cm
    'Height': [15, 40],    # cm
}

print("Engineering Model: Beam Stress Analysis")
print("Parameter ranges:")
for param, (low, high) in param_ranges.items():
    print(f"  {param}: [{low}, {high}]")

# Generate samples and run GSA
n_samples = 2000
X_eng = np.zeros((n_samples, 4))

# Sample parameters uniformly
for i, (param, (low, high)) in enumerate(param_ranges.items()):
    X_eng[:, i] = np.random.uniform(low, high, n_samples)

# Evaluate model
stress = engineering_model(X_eng)
Y_eng = stress.reshape(-1, 1)

print(f"Stress range: [{stress.min():.1f}, {stress.max():.1f}] MPa")

# Run GSA
param_names_eng = list(param_ranges.keys())
results_eng = gsa_pipeline(
    X_eng, Y_eng,
    param_names=param_names_eng,
    feature_names=["Stress_MPa"],
    scaler="minmax",
    enable_sobol=True,
    make_pdp=False
)

print("\nBeam Stress Sensitivity Analysis:")
eng_table = results_eng["results"]["Stress_MPa"]["table"]
eng_extras = results_eng["results"]["Stress_MPa"]["models"]
print(eng_table.round(4))

print(f"\nFit diagnostics:")
print(f"  surrogate collapsed     : {eng_extras['surrogate_collapsed']}")
print(f"  log marginal likelihood : {eng_extras['log_marginal_likelihood']:.1f}")
print(f"  length scales on a bound: {eng_extras['ard_at_bounds']}")
print(f"  Sobol out of range      : {eng_extras['sobol_out_of_range']}")
print(f"  sum(S1)                 : {eng_table['S1'].sum():.4f}")

print("\nPhysically, stress ~ M*(h/2)/I with I ~ w*h^3, so stress ~ 1/h^2:")
print("Height should dominate, then Length (squared in the moment), then Load")
print("(linear), then Width (linear in I). Check the table against that -- a GSA")
print("result you cannot explain from the model is a result worth doubting.")

# Compute expected stress using integration
def uniform_pdf(theta):
    """Uniform PDF over parameter ranges"""
    # Volume of parameter space
    volume = 1.0
    for low, high in param_ranges.values():
        volume *= (high - low)
    return np.ones(len(theta)) / volume

def stress_integrand(theta):
    """Stress function for integration"""
    return engineering_model(theta).reshape(-1, 1)

# Compute expected stress
bounds = list(param_ranges.values())
expected_stress = qmc_integral_auto(
    N_samples=4096,
    N_params=4,
    data_func=stress_integrand,
    p_target=uniform_pdf,
    bounds=bounds,
    seed=0,
)

exp_val = float(np.ravel(expected_stress['integral'])[0])
print(f"\nExpected stress: {exp_val:.2f} MPa")
print(f"Sample mean stress: {stress.mean():.2f} MPa")
print(f"Difference: {abs(exp_val - stress.mean()):.2f} MPa")
print("\nThese two should agree: p_target is uniform over the same box the QMC")
print("points fill, so the integral IS the mean. Before v0.2.0 this printed a")
print("value larger by the box volume -- here a factor of 40*6*20*25 = 1.2e6.")