1import json
  2import time
  3from pathlib import Path
  4import numpy as np
  5from scipy.sparse import diags
  6from scipy.linalg import solve
  7
  8SEED = 2403
  9rng = np.random.default_rng(SEED)
 10OUT = Path("results.json")
 11
 12
 13def stiffness(n):
 14    """Dirichlet P1 stiffness on n interior nodes, with h-scaled nodal values."""
 15    h = 1.0 / (n + 1)
 16    K = (1.0 / h) * (2*np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1))
 17    x = np.arange(1, n+1) * h
 18    return x, K
 19
 20
 21def prolongation(nc, nf):
 22    """P1 injection from nc interior nodes to nf interior nodes, requiring nested grids."""
 23    assert (nf + 1) % (nc + 1) == 0
 24    ratio = (nf + 1) // (nc + 1)
 25    P = np.zeros((nf, nc))
 26    for j in range(1, nc+1):
 27        coarse_x = j * ratio
 28        for i in range(max(1, coarse_x-ratio), min(nf, coarse_x+ratio)+1):
 29            P[i-1, j-1] = max(0.0, 1.0 - abs(i-coarse_x)/ratio)
 30    return P
 31
 32
 33def make_data(n=128):
 34    x, K = stiffness(n)
 35    # Manufactured exact solution u*=sin(pi x); f=-u''=pi^2 sin(pi x).
 36    ustar = np.sin(np.pi*x)
 37    f = K @ ustar  # exact for the discrete problem, avoids load quadrature error
 38    return x, K, ustar, f
 39
 40
 41def riesz_score(u, K, f, P):
 42    A = P.T @ K @ P
 43    b = P.T @ (f - K @ u)
 44    z = solve(A, b, assume_a="pos")
 45    return float(np.sqrt(max(0.0, z @ A @ z)))
 46
 47
 48def energy_error(u, ustar, K):
 49    e = u-ustar
 50    return float(np.sqrt(max(0.0, e @ K @ e)))
 51
 52
 53def run():
 54    n = 127
 55    x, K, ustar, f = make_data(n)
 56    # Candidate archive: progressively improving solutions plus oscillatory and smooth errors.
 57    modes = [np.sin(np.pi*x), np.sin(2*np.pi*x), np.sin(7*np.pi*x), np.sin(15*np.pi*x)]
 58    archive = []
 59    for k in range(100):
 60        t = k/99.0
 61        # Main error decays, while a few checkpoints have low-frequency bias that is
 62        # hard for small auxiliary spaces to detect.
 63        amp = 0.95*(1-t) + 0.025
 64        e = amp*(0.55*modes[1] + 0.25*modes[2])
 65        e += (0.10 + 0.35*(1-t))*modes[3]
 66        e += 0.018*rng.normal(size=n)
 67        if k in (16, 48, 76):
 68            e += 0.16*modes[1]
 69        archive.append(ustar + e)
 70    archive = np.asarray(archive)
 71    oracle = np.array([energy_error(u, ustar, K) for u in archive])
 72
 73    levels = [3, 7, 15, 31, 63]
 74    # Full auxiliary space is the exact finite-dimensional Riesz representative.
 75    full_scores = np.array([riesz_score(u, K, f, np.eye(n)) for u in archive])
 76    full_identity_max_relerr = float(np.max(np.abs(full_scores-oracle) / np.maximum(oracle, 1e-15)))
 77    scores = {}
 78    ratios = {}
 79    for nc in levels:
 80        P = prolongation(nc, n)
 81        scores[nc] = np.array([riesz_score(u, K, f, P) for u in archive])
 82        ratios[nc] = scores[nc] / oracle
 83
 84    # Prediction 1: Galerkin score is a lower bound for every candidate.
 85    lower_bound = {str(m): float(np.max(scores[m] <= oracle + 2e-10)) for m in levels}
 86    # Prediction 2: nested auxiliary spaces are pointwise monotone.
 87    monotone = {str(a): {
 88        "min_increment": float(np.min(scores[b] - scores[a])),
 89        "max_monotonicity_violation": float(np.max(np.maximum(0.0, scores[a] - scores[b])))
 90    } for a,b in zip(levels[:-1], levels[1:])}
 91    # Prediction 3: refinement converges to the exact energy error; report relative gap.
 92    convergence = {str(m): float(np.mean(np.abs(1.0-ratios[m]))) for m in levels}
 93    max_gap = {str(m): float(np.max(1.0-ratios[m])) for m in levels}
 94
 95    # Ranking quality and selection comparison.
 96    oracle_rank = np.argsort(oracle)
 97    rank_corr = {}
 98    selected = {}
 99    for m in levels:
100        # Spearman correlation computed with rank arrays, no scipy.stats dependency needed.
101        r1 = np.empty(100, dtype=int); r1[np.argsort(oracle)] = np.arange(100)
102        r2 = np.empty(100, dtype=int); r2[np.argsort(scores[m])] = np.arange(100)
103        rank_corr[str(m)] = float(np.corrcoef(r1, r2)[0,1])
104        j = int(np.argmin(scores[m]))
105        selected[str(m)] = {"checkpoint": j, "score": float(scores[m][j]), "oracle_error": float(oracle[j])}
106
107    # A deliberately realistic raw training loss proxy: sampled pointwise residual,
108    # with noise, and a pointwise residual monitor on a separate coarse sample.
109    raw_loss = np.array([np.mean((f-K@u)**2) for u in archive])
110    raw_loss += rng.normal(0, 0.03*np.std(raw_loss), size=100)
111    raw_j = int(np.argmin(raw_loss))
112    # Measure the actual cheap post-processing wall time on this tiny problem.
113    t0 = time.perf_counter()
114    for nc in levels:
115        P = prolongation(nc, n)
116        _ = [riesz_score(u, K, f, P) for u in archive]
117    riesz_seconds = time.perf_counter() - t0
118    # Riesz postprocessing cost is measured as number of tiny solves.
119    result = {
120        "seed": SEED, "n_fine": n, "levels": levels,
121        "full_space_identity_max_relative_error": full_identity_max_relerr,
122        "riesz_postprocess_seconds": riesz_seconds,
123        "predictions": {
124            "lower_bound_fraction": lower_bound,
125            "nested_min_increment": monotone,
126            "mean_relative_gap_to_oracle": convergence,
127            "max_relative_gap_to_oracle": max_gap,
128            "interpretation": "All scores must be <= oracle; nested increments must be >= 0; gaps should decrease with level."
129        },
130        "selection": {
131            "raw_training_loss": {"checkpoint": raw_j, "oracle_error": float(oracle[raw_j])},
132            "riesz": selected,
133            "oracle_best": {"checkpoint": int(np.argmin(oracle)), "oracle_error": float(np.min(oracle))},
134            "spearman_rank_correlation": rank_corr
135        },
136        "oracle_error_range": [float(np.min(oracle)), float(np.max(oracle))],
137        "postprocess_solves": len(levels)*len(archive)
138    }
139    OUT.write_text(json.dumps(result, indent=2))
140    print(json.dumps(result, indent=2))
141
142if __name__ == "__main__":
143    run()