Rewording several parts of the findings to add clarity.

ergodic_sampling_test.py

@@ -4,28 +4,48 @@ from scipy.interpolate import griddata
 from scipy.stats import pearsonr, spearmanr
 import matplotlib.pyplot as plt
 
+# ---------------- 1. parameters ----------------
 rng = np.random.default_rng(0)
+
+# grid
 NX = NY = 64; N = NX*NY
+
+# ground truths
 BETA = 2.5                      # potential-field truth: power spectrum ~ k^-BETA (before attenuation)
 Z_SRC = 1.5                     # source depth (cells): upward-continuation attenuation exp(-2*pi*k*Z_SRC)
 N_MODES = 12                    # Fourier-sparse truth: number of modes (low-frequency biased)
+
+# sampling patterns
 STRIDE = 3                      # sparse (regular-grid) stride -> delta ~ 0.118
+NBR, NBA, RMAX, BLK = 45, 36, 90.0, 4   # BLK=4 ~ paper's patch size A/N_phi (2.9 px, non-integer)
+SA_ITERS = 4000
+
+# reconstruction solvers
 WAVELET, LEVELS, MODE = 'sym8', 2, 'periodization'   # 2 = pywt.dwt_max_level(64, 'sym8')
 SYM_SPINS = 8                   # cycle-spin shifts for the symlet reconstruction (our choice; the
                                 # paper never specifies cycle spinning or any solver parameters)
-NBR, NBA, RMAX, BLK = 45, 36, 90.0, 4   # BLK=4 ~ paper's patch size A/N_phi (2.9 px, non-integer)
-SA_ITERS = 4000
 POCS_ITERS, POCS_LMIN = 800, 1e-4
+
+# validation
 LOO_M = 60                      # leave-one-out calibration points
 ENS_K = 64                      # ensemble validation realizations (our construct, not the paper's)
+N_TRIALS = 25                   # Monte-Carlo trials per truth class
+RUN_SANITY_CHECK = False        # plt4 served its purpose; flip to True to re-run the check
+
+# ---------------- 2. ground truths (two tiers) ----------------
+# spectral grid shared by the truth generator and the Gaussian-process covariance
+KRAD = np.hypot(np.fft.fftfreq(NX)[:, None], np.fft.fftfreq(NY)[None, :])
+
+def power_spectrum(beta, z_src):
+    # power spectrum k^-beta * exp(-4*pi*k*z): squared amplitude of the attenuated generator
+    S = np.where(KRAD > 0, KRAD, 1.0)**(-beta)*np.exp(-4*np.pi*KRAD*z_src)
+    S[0, 0] = 0.0
+    return S
 
-# ---------------- ground truths (two tiers) ----------------
 def powerlaw_field(beta, z_src=Z_SRC):
     # random-phase power law with source-depth attenuation: without the exp(-2*pi*k*z) factor the
     # field is rough at pixel scale, which real (upward-continued) potential-field data never is
-    kx = np.fft.fftfreq(NX)[:, None]; ky = np.fft.fftfreq(NY)[None, :]
+    amp = np.sqrt(power_spectrum(beta, z_src))
-    k = np.hypot(kx, ky); k[0, 0] = 1.0
-    amp = k**(-beta/2)*np.exp(-2*np.pi*k*z_src); amp[0, 0] = 0.0
     f = np.real(np.fft.ifft2(amp*np.exp(2j*np.pi*rng.random((NX, NY)))))
     return (f - f.mean())/f.std()
 
@@ -49,10 +69,13 @@ def fourier_sparse_field(n_modes):
     f = np.fft.ifft2(C).real
     return (f - f.mean())/f.std()
 
-truth_pl = powerlaw_field(BETA)           # tier 2: hard, not l1-sparse in any basis
+def white_noise_field():
-truth_sp = fourier_sparse_field(N_MODES)  # tier 1: exactly sparse, the paper's Fig. A.1 setting
+    # incompressible control: no method can beat predicting zero at unsampled cells,
+    # whose rms is sqrt(1 - DELTA) ~ 0.94
+    f = rng.standard_normal((NX, NY))
+    return (f - f.mean())/f.std()
 
-# ---------------- sampling patterns ----------------
+# ---------------- 3. sampling patterns and ISA properties (Eq. 4 ingredients) ----------------
 mask_reg = np.zeros((NX, NY), bool); mask_reg[::STRIDE, ::STRIDE] = True
 NS = int(mask_reg.sum()); DELTA = NS/N
 
@@ -60,36 +83,38 @@ def random_mask():
     m = np.zeros(N, bool); m[rng.choice(N, NS, replace=False)] = True
     return m.reshape(NX, NY)
 
-# ---------------- ISA properties (Eq. 4 ingredients) ----------------
 def H(p):
+    # entropy of the absolute deviations (Eq. 4 applies it to histogram differences)
     p = np.abs(p); nz = p > 0
     return -np.sum(p[nz]*np.log(p[nz]))
 
+def pair_hists(dx, dy, weights=None):
+    hr, _ = np.histogram(np.hypot(dx, dy), NBR, (0, RMAX), weights=weights)
+    # angles folded to [0, pi): the all-pairs rose diagram is 180-degree symmetric
+    ha, _ = np.histogram(np.arctan2(dy, dx) % np.pi, NBA, (0, np.pi), weights=weights)
+    return hr/hr.max(), ha/ha.max()
+
 def dense_pair_hists(nx, ny):
     # all-pairs interval/angle histograms of the dense grid via closed-form pair counts
     DX, DY = np.meshgrid(np.arange(-(nx-1), nx), np.arange(0, ny), indexing='ij')
     keep = (DY > 0) | ((DY == 0) & (DX > 0))
     DX, DY = DX[keep], DY[keep]
     w = (nx - np.abs(DX))*(ny - np.abs(DY))
-    hr, _ = np.histogram(np.hypot(DX, DY), NBR, (0, RMAX), weights=w)
+    return pair_hists(DX, DY, w)
-    # angles folded to [0, pi): the all-pairs rose diagram is 180-degree symmetric
-    ha, _ = np.histogram(np.arctan2(DY, DX) % np.pi, NBA, (0, np.pi), weights=w)
-    return hr/hr.max(), ha/ha.max()
 
 def sparse_pair_hists(mask):
     c = np.argwhere(mask).astype(float)
     d = c[:, None, :] - c[None, :, :]
     iu = np.triu_indices(len(c), 1)
-    dx, dy = d[..., 0][iu], d[..., 1][iu]
+    return pair_hists(d[..., 0][iu], d[..., 1][iu])
-    hr, _ = np.histogram(np.hypot(dx, dy), NBR, (0, RMAX))
-    ha, _ = np.histogram(np.arctan2(dy, dx) % np.pi, NBA, (0, np.pi))
-    return hr/hr.max(), ha/ha.max()
 
 def block_density(mask):
     # per-block occupancy fraction gamma_phi; compared against delta*gamma_theta = DELTA (Eq. 2)
     return mask.reshape(NX//BLK, BLK, NY//BLK, BLK).sum((1, 3))/BLK**2
 
 def srf(mask):
+    # spatial response function of the pattern: mu = largest side lobe (the mutual coherence),
+    # xi = the full side-lobe vector
     P = np.abs(np.fft.fft2(mask.astype(float))); P /= P[0, 0]
     xi = P.ravel()[1:]
     return xi.max(), xi
@@ -117,12 +142,7 @@ def ergodic_mask(iters=SA_ITERS, T0=0.05):
             m.flat[i] = True; m.flat[j] = False
     return best[1], best[0]
 
-mask_erg, f_erg = ergodic_mask()
+# ---------------- 4. transforms and compressive-sensing solvers ----------------
-f_rand = np.mean([objective(random_mask()) for _ in range(20)])
-print(f"N_samples={NS} (delta={DELTA:.3f})  obj: sparse={objective(mask_reg):.3f}  "
-      f"random(mean)={f_rand:.3f}  ergodic={f_erg:.3f}")
-
-# ---------------- transforms ----------------
 _, SLICES = pywt.coeffs_to_array(pywt.wavedec2(np.zeros((NX, NY)), WAVELET, mode=MODE, level=LEVELS))
 def Wsym(x):  return pywt.coeffs_to_array(pywt.wavedec2(x, WAVELET, mode=MODE, level=LEVELS))[0]
 def Wtsym(c): return pywt.waverec2(pywt.array_to_coeffs(c, SLICES, output_format='wavedec2'), WAVELET, mode=MODE)
@@ -130,67 +150,67 @@ SYM_APP = np.zeros((NX, NY), bool); SYM_APP[SLICES[0]] = True    # approximation
 def Wfft(x):  return np.fft.fft2(x)
 def Wtfft(c): return np.fft.ifft2(c).real
 
-# ---------------- CS reconstruction: iterative thresholding + data reinsertion (POCS) ----------------
+def cs_reconstruct(y_img, mask, W=Wfft, Wt=Wtfft, protect=None,
-def cs_reconstruct(y_img, mask, W=Wfft, Wt=Wtfft, n_iter=POCS_ITERS, lam_min_frac=POCS_LMIN):
+                   n_iter=POCS_ITERS, lam_min_frac=POCS_LMIN):
+    # iterative thresholding + data reinsertion (POCS); coefficients flagged in `protect` are
+    # never thresholded and excluded from the threshold schedule
     x = y_img.copy()
-    lam_max = np.abs(W(y_img)).max()
+    c0 = W(y_img)
+    lam_max = np.abs(c0 if protect is None else np.where(protect, 0, c0)).max()
     for lam in np.geomspace(lam_max, lam_max*lam_min_frac, n_iter):
         c = W(x)
-        c *= np.maximum(1 - lam/np.maximum(np.abs(c), 1e-300), 0)   # soft threshold on magnitude
+        t = c*np.maximum(1 - lam/np.maximum(np.abs(c), 1e-300), 0)   # soft threshold on magnitude
-        x = Wt(c)
+        x = Wt(t if protect is None else np.where(protect, c, t))
         x[mask] = y_img[mask]
     return x
 
 # symlet variant (the paper's transform), tuned to reproduce the paper's accuracy in its regime:
 # never threshold the approximation subband, and average reconstructions over circular shifts
-# (cycle spinning) to remove the decimated transform's shift-variance
+# (cycle spinning) to remove the decimated transform's shift-variance. For the periodized
-def cs_reconstruct_sym(y_img, mask, n_iter=POCS_ITERS, lam_min_frac=POCS_LMIN, n_shift=SYM_SPINS):
+# 2-level transform, shifting data and mask together by any multiple of 4 gives an exactly
-    def single(y2, m2):
+# equivalent problem, so only the shift residues modulo 2^LEVELS matter; the fixed list below
-        x = y2.copy()
+# covers distinct residue pairs (randomness is not needed and keeps the solver deterministic).
-        lam_max = np.abs(np.where(SYM_APP, 0, Wsym(y2))).max()
+SYM_SHIFTS = [(0, 0), (1, 2), (2, 1), (3, 3), (0, 2), (1, 1), (2, 3), (3, 0)][:SYM_SPINS]
-        for lam in np.geomspace(lam_max, lam_max*lam_min_frac, n_iter):
+
-            c = Wsym(x)
+def cs_reconstruct_sym(y_img, mask, n_iter=POCS_ITERS, lam_min_frac=POCS_LMIN):
-            t = np.sign(c)*np.maximum(np.abs(c) - lam, 0)
-            x = Wtsym(np.where(SYM_APP, c, t))
-            x[m2] = y2[m2]
-        return x
     recs = []
-    for _ in range(n_shift):
+    for s in SYM_SHIFTS:
-        s = (int(rng.integers(NX)), int(rng.integers(NY)))
+        xs = cs_reconstruct(np.roll(y_img, s, (0, 1)), np.roll(mask, s, (0, 1)),
-        xs = single(np.roll(y_img, s, (0, 1)), np.roll(mask, s, (0, 1)))
+                            Wsym, Wtsym, SYM_APP, n_iter, lam_min_frac)
         recs.append(np.roll(xs, (-s[0], -s[1]), (0, 1)))
     x = np.mean(recs, 0)
     x[mask] = y_img[mask]
     return x
 
-# ---------------- kriging / GP: Bayes-optimal for the power-law field class ----------------
+# ---------------- 5. kriging / GP: Bayes-optimal for the power-law field class ----------------
-KXF = np.fft.fftfreq(NX)[:, None]; KYF = np.fft.fftfreq(NY)[None, :]
-KRAD = np.hypot(KXF, KYF)
 GRID = np.argwhere(np.ones((NX, NY), bool))
 
 def circulant_cov(beta, z_src):
-    # power spectrum k^-beta * exp(-4*pi*k*z): the amplitude-attenuated generator's power spectrum
+    c = np.fft.ifft2(power_spectrum(beta, z_src)).real
-    S = np.where(KRAD > 0, KRAD, 1.0)**(-beta)*np.exp(-4*np.pi*KRAD*z_src); S[0, 0] = 0.0
-    c = np.fft.ifft2(S).real
     return c/c[0, 0]
 
+def sample_cov(on, cov):
+    # sample-sample covariance with a small jitter for numerical stability
+    K = cov[(on[:, None, 0]-on[None, :, 0]) % NX, (on[:, None, 1]-on[None, :, 1]) % NY]
+    return K + 1e-6*np.eye(len(on))
+
 def gp_matrices(mask, beta, z_src):
     on = np.argwhere(mask)
     cov = circulant_cov(beta, z_src)
-    Kss = cov[(on[:, None, 0]-on[None, :, 0]) % NX, (on[:, None, 1]-on[None, :, 1]) % NY]
     Kxs = cov[(GRID[:, None, 0]-on[None, :, 0]) % NX, (GRID[:, None, 1]-on[None, :, 1]) % NY]
-    return Kss + 1e-6*np.eye(len(on)), Kxs
+    return sample_cov(on, cov), Kxs
+
+def gp_unit_variance(Kss, Kxs):
+    # posterior variance at every cell for a unit-variance field
+    return np.maximum(1.0 - np.einsum('ij,ji->i', Kxs, np.linalg.solve(Kss, Kxs.T)), 0)
 
 def gp_fit_params(y_img, mask):
     # truth-free ML over spectral slope and source depth, variance profiled out
     # (only the sample-sample covariance is needed here, not the full-grid cross-covariance)
-    yv = y_img[mask]; best = None
+    yv = y_img[mask]; on = np.argwhere(mask); best = None
-    on = np.argwhere(mask)
     for beta in np.arange(1.5, 4.01, 0.5):
         for z_src in (0.0, 0.5, 1.0, 1.5, 2.0, 3.0):
-            cov = circulant_cov(beta, z_src)
+            Kss = sample_cov(on, circulant_cov(beta, z_src))
-            Kss = cov[(on[:, None, 0]-on[None, :, 0]) % NX,
-                      (on[:, None, 1]-on[None, :, 1]) % NY] + 1e-6*np.eye(len(on))
             try:
                 L = np.linalg.cholesky(Kss)
             except np.linalg.LinAlgError:
@@ -210,10 +230,8 @@ def gp_predictor(mask, beta, z_src):
 def gp_reconstruct(y_img, mask, beta, z_src):
     yv = y_img[mask]
     Kss, Kxs = gp_matrices(mask, beta, z_src)
-    w = np.linalg.solve(Kss, yv)
+    xh = (Kxs @ np.linalg.solve(Kss, yv)).reshape(NX, NY)
-    xh = (Kxs @ w).reshape(NX, NY)
+    var = gp_unit_variance(Kss, Kxs)
-    A = np.linalg.solve(Kss, Kxs.T)
-    var = np.maximum(1.0 - np.einsum('ij,ji->i', Kxs, A), 0)
     s2 = yv @ np.linalg.solve(Kss, yv)/len(yv)          # profiled variance scale
     return xh, np.sqrt(var*s2).reshape(NX, NY)
 
@@ -224,8 +242,8 @@ def apriori_level_correction(beta, z_src, mask, n_top=400, m_draws=20000):
     # Y = sum_p s_p u_p (the realization variance), where u_p = 1 + cos(theta_p) are the random
     # mode powers. Both expectations involve only scalar mode powers, so the level bias of the
     # posterior sigma is computable a priori from the fitted spectrum and the mask alone.
-    Sfit = np.where(KRAD > 0, KRAD, 1.0)**(-beta)*np.exp(-4*np.pi*KRAD*z_src); Sfit[0, 0] = 0.0
+    Sn = power_spectrum(beta, z_src)
-    Sn = (Sfit/Sfit.sum()).ravel()
+    Sn = (Sn/Sn.sum()).ravel()
     idx = np.arange(N); ix, iy = np.divmod(idx, NY)
     neg = ((-ix) % NX)*NY + ((-iy) % NY)
     take = ((idx < neg) | (idx == neg)) & (Sn > 0)
@@ -243,7 +261,7 @@ def apriori_level_correction(beta, z_src, mask, n_top=400, m_draws=20000):
     g_tail = (1/Y).mean()                               # small modes are independent of Y
     Kss_c, Kxs_c = gp_matrices(mask, beta, z_src)
     Wc = Kxs_c @ np.linalg.inv(Kss_c)
-    sig2_u = np.maximum(1.0 - np.einsum('ij,ji->i', Kxs_c, np.linalg.solve(Kss_c, Kxs_c.T)), 0)
+    sig2_u = gp_unit_variance(Kss_c, Kxs_c)
     ph = np.exp(2j*np.pi*(np.outer(GRID[:, 0], kxp[:n_top])/NX
                           + np.outer(GRID[:, 1], kyp[:n_top])/NY))
     resp = Wc @ ph[mask.ravel()] - ph                   # error response of each dominant mode
@@ -252,15 +270,15 @@ def apriori_level_correction(beta, z_src, mask, n_top=400, m_draws=20000):
     sig2_c = t_top @ g_top + t_tail*g_tail
     return np.sqrt(sig2_c/np.maximum(sig2_u, 1e-30)).reshape(NX, NY)
 
-# ---------------- acquire + reconstruct ----------------
+# ---------------- 6. metrics, baselines, and validation helpers ----------------
 def leakage(x, truth): return np.corrcoef((x - truth).ravel(), truth.ravel())[0, 1]
 def rms(x, truth):     return np.sqrt(np.mean((x - truth)**2))
 
-def spline_interp(truth, mask):
+def spline_interp(img, mask):
     pts = np.argwhere(mask); X, Y = np.mgrid[0:NX, 0:NY]
-    xr = griddata(pts, truth[mask], (X, Y), method='cubic')
+    xr = griddata(pts, img[mask], (X, Y), method='cubic')
     nn = np.isnan(xr)
-    xr[nn] = griddata(pts, truth[mask], (X, Y), method='nearest')[nn]
+    xr[nn] = griddata(pts, img[mask], (X, Y), method='nearest')[nn]
     return xr
 
 def fourier_oracle(truth, k=NS):
@@ -268,55 +286,6 @@ def fourier_oracle(truth, k=NS):
     C[np.argsort(np.abs(C))[:-k]] = 0
     return np.fft.ifft2(C.reshape(NX, NY)).real
 
-# tier 1: Fourier-sparse truth: the CS claim of the paper (Fig. A.1 / Figs. 15, 22)
-print("\n--- tier 1: Fourier-sparse truth (paper's CS regime) ---")
-x_sp_spl = spline_interp(truth_sp, mask_reg)
-x_sp_erg_spl = spline_interp(truth_sp, mask_erg)
-x_sp_reg = cs_reconstruct(truth_sp*mask_reg, mask_reg)
-x_sp_erg = cs_reconstruct(truth_sp*mask_erg, mask_erg)
-x_sp_sym = cs_reconstruct_sym(truth_sp*mask_erg, mask_erg)   # paper's literal transform
-x_sp_reg_sym = cs_reconstruct_sym(truth_sp*mask_reg, mask_reg)
-# kriging with the power-law covariance is misspecified for a 12-mode signal; shown as the
-# symmetric counterpart to CS-on-the-random-field so every process appears on both truths
-GP_SP = gp_fit_params(truth_sp*mask_erg, mask_erg)
-x_sp_gp, sig_sp_raw = gp_reconstruct(truth_sp*mask_erg, mask_erg, *GP_SP)
-x_sp_reg_gp = gp_reconstruct(truth_sp*mask_reg, mask_reg, *GP_SP)[0]
-for name, x in [('sparse+spline', x_sp_spl), ('sparse+CS', x_sp_reg),
-                ('sparse+CS-symlet', x_sp_reg_sym), ('sparse+kriging', x_sp_reg_gp),
-                ('ergodic+spline', x_sp_erg_spl), ('ergodic+CS', x_sp_erg),
-                ('ergodic+CS-symlet', x_sp_sym), ('ergodic+kriging', x_sp_gp)]:
-    print(f"{name:18s} leakage={leakage(x, truth_sp):+.3f}  rms={rms(x, truth_sp):.3f}")
-
-# tier 2: potential-field truth: not l1-sparse; report against the oracle floor
-print(f"\n--- tier 2: potential-field truth (beta={BETA}, source depth {Z_SRC}) ---")
-print(f"oracle {NS}-term Fourier floor: rms={rms(fourier_oracle(truth_pl), truth_pl):.3f}")
-BETA_HAT, Z_HAT = gp_fit_params(truth_pl*mask_erg, mask_erg)
-print(f"GP ML estimates: beta_hat={BETA_HAT:.1f} z_hat={Z_HAT:.1f} "
-      f"(generator beta={BETA}, z={Z_SRC})")
-CORR_PL = apriori_level_correction(BETA_HAT, Z_HAT, mask_erg)
-print(f"a-priori sigma level correction for the normalized truth generator: "
-      f"x{CORR_PL[~mask_erg].mean():.4f} (spread {CORR_PL[~mask_erg].std():.4f})")
-y_erg = truth_pl*mask_erg
-x_pl_spl = spline_interp(truth_pl, mask_reg)
-x_pl_erg_spl = spline_interp(truth_pl, mask_erg)
-x_pl_reg = cs_reconstruct(truth_pl*mask_reg, mask_reg)
-x_pl_reg_sym = cs_reconstruct_sym(truth_pl*mask_reg, mask_reg)
-x_pl_reg_gp = gp_reconstruct(truth_pl*mask_reg, mask_reg, BETA_HAT, Z_HAT)[0]
-x_erg = cs_reconstruct(y_erg, mask_erg)                      # primary reconstruction (paper's method)
-x_erg_sym = cs_reconstruct_sym(y_erg, mask_erg)              # symlet-wavelet comparison
-x_gp, sig_gp_raw = gp_reconstruct(y_erg, mask_erg, BETA_HAT, Z_HAT)
-print("method            leakage  rms    -rms/std(truth)   (leakage ~ -rms/std when error is unrecovered signal)")
-for name, x in [('sparse+spline', x_pl_spl),
-                ('sparse+CS', x_pl_reg),
-                ('sparse+CS-symlet', x_pl_reg_sym),
-                ('sparse+kriging', x_pl_reg_gp),
-                ('ergodic+spline', x_pl_erg_spl),
-                ('ergodic+CS-Fourier', x_erg),
-                ('ergodic+CS-symlet', x_erg_sym),
-                ('ergodic+kriging', x_gp)]:
-    print(f"{name:18s} {leakage(x, truth_pl):+.3f}  {rms(x, truth_pl):.3f}   {-rms(x, truth_pl)/truth_pl.std():+.3f}")
-
-# ---------------- predicted sigma (truth-free: from y + mask only) ----------------
 # LOO calibration: reconstruct without one sample, compare prediction at that point to its value
 def loo_rms(y_img, mask, rec_fn, m=LOO_M):
     on = np.argwhere(mask)
@@ -337,8 +306,6 @@ def spearman_ceiling(s, ndraw=50):
     return np.mean([spearmanr(s.ravel()*np.abs(rng.standard_normal(s.size)), s.ravel())[0]
                     for _ in range(ndraw)])
 
-OFF = ~mask_erg      # evaluate z-scores at unsampled pixels only (err and sigma are both 0 at samples)
-
 def reliability(s, a, nb=10, reduce=np.mean):
     q = np.quantile(s, np.linspace(0, 1, nb+1)); q[-1] += 1e-12
     m, r = [], []
@@ -350,58 +317,51 @@ def reliability(s, a, nb=10, reduce=np.mean):
 
 def rms_reduce(a): return np.sqrt(np.mean(a**2))
 
-# ---------------- plots (no abbreviations in any figure text) ----------------
+# ---------------- 7. method registry and figure helpers ----------------
-C_ENS, C_ONE = '#1f77b4', '#ff7f0e'   # blue = ensemble comparison, orange = single realization
+# every pattern-method combination shown in tables and figures; the key prefix names the
+# sampling pattern, resolved through MASKS once the ergodic mask has been optimized
+METHODS = {  # key -> (figure column title, Monte-Carlo bar label)
+    'sparse+spline':      ('Sparse + spline interpolation\n(naive control)',
+                           'Sparse + spline interpolation'),
+    'sparse+CS-Fourier':  ('Sparse + compressive sensing\n(Fourier)',
+                           'Sparse + compressive sensing (Fourier)'),
+    'sparse+CS-symlet':   ('Sparse + compressive sensing\n(symlet wavelet, from paper)',
+                           'Sparse + compressive sensing (symlet wavelet)'),
+    'sparse+kriging':     ('Sparse + kriging\n(Gaussian process)',
+                           'Sparse + kriging (Gaussian process)'),
+    'ergodic+spline':     ('Ergodic + spline interpolation',
+                           'Ergodic + spline interpolation'),
+    'ergodic+CS-Fourier': ('Ergodic + compressive sensing\n(Fourier)',
+                           'Ergodic + compressive sensing (Fourier)'),
+    'ergodic+CS-symlet':  ('Ergodic + compressive sensing\n(symlet wavelet, from paper)',
+                           'Ergodic + compressive sensing (symlet wavelet)'),
+    'ergodic+kriging':    ('Ergodic + kriging\n(Gaussian process)',
+                           'Ergodic + kriging (Gaussian process)'),
+}
+MC_KEYS = [k for k in METHODS if k != 'sparse+CS-symlet']   # Monte-Carlo method set
 
-# --- figures 1a/1b: one figure per ground truth: reconstructions, errors, sampling locations ---
+def run_methods(truth, krig_sparse, krig_erg, keys=None):
-COLS = ['Ground truth',
+    # all reconstructions of one truth; the kriging callables are supplied by the caller
-        'Sparse + spline interpolation\n(naive control)',
+    # (gp_reconstruct closures in the tier sections, precomputed weights in the ensembles)
-        'Sparse + compressive sensing\n(Fourier)',
+    out = {}
-        'Sparse + compressive sensing\n(symlet wavelet, from paper)',
+    for k in (keys or list(METHODS)):
-        'Sparse + kriging\n(Gaussian process)',
+        pattern, method = k.split('+')
-        'Ergodic + spline interpolation',
+        m = MASKS[pattern]; y = truth*m
-        'Ergodic + compressive sensing\n(Fourier)',
+        if method == 'spline':       out[k] = spline_interp(y, m)
-        'Ergodic + compressive sensing\n(symlet wavelet, from paper)',
+        elif method == 'CS-Fourier': out[k] = cs_reconstruct(y, m)
-        'Ergodic + kriging\n(Gaussian process)']
+        elif method == 'CS-symlet':  out[k] = cs_reconstruct_sym(y, m)
-COL_MASKS = [np.ones_like(mask_reg), mask_reg, mask_reg, mask_reg, mask_reg,
+        else:                        out[k] = (krig_sparse if pattern == 'sparse' else krig_erg)(y, m)
-             mask_erg, mask_erg, mask_erg, mask_erg]
+    return out
-for fname, gt, imgs, note in [
+
-        ('plt1a', truth_sp,
+# plot colors (figure text carries no abbreviations, a standing convention of this project)
-         [truth_sp, x_sp_spl, x_sp_reg, x_sp_reg_sym, x_sp_reg_gp,
+C_ENS, C_ONE = '#1f77b4', '#ff7f0e'         # blue = ensemble comparison, orange = single realization
-          x_sp_erg_spl, x_sp_erg, x_sp_sym, x_sp_gp],
+C_LIGHT_BLUE, C_LIGHT_ORANGE = '#9ecae1', '#fdbe85'   # scatter points / pattern bar colors
-         f'Fourier-sparse truth ({N_MODES} modes); compressive sensing with the ergodic pattern '
+
-         'recovers the signal exactly'),
+def plot_mask_panel(a, mm, title):
-        ('plt1b', truth_pl,
-         [truth_pl, x_pl_spl, x_pl_reg, x_pl_reg_sym, x_pl_reg_gp,
-          x_pl_erg_spl, x_erg, x_erg_sym, x_gp],
-         f'potential-field truth (spectral exponent {BETA}, source depth {Z_SRC} cells); '
-         'kriging is the best reconstruction for this field class')]:
-    fig, axes = plt.subplots(3, 9, figsize=(36, 12.4), constrained_layout=True)
-    v = np.abs(gt).max()
-    for c, (a, x, t) in enumerate(zip(axes[0], imgs, COLS)):
-        im = a.imshow(x.T, origin='lower', vmin=-v, vmax=v, cmap='RdBu_r')
-        sub = '' if c == 0 else (f'\nleakage {leakage(x, gt):+.2f} · '
-                                 f'RMS error {rms(x, gt):.3f}')
-        a.set_title(t + sub, fontsize=9)
-    axes[0, 0].set_ylabel('Reconstruction\n\ngrid y (cells)')
-    fig.colorbar(im, ax=axes[0], shrink=0.9,
-                 label='field value\n(units of the truth standard deviation)')
-    errors = [gt - x for x in imgs[1:]]
-    ve = max(np.abs(e).max() for e in errors)
-    axes[1, 0].axis('off')
-    for a, e in zip(axes[1, 1:], errors):
-        im2 = a.imshow(e.T, origin='lower', vmin=-ve, vmax=ve, cmap='RdBu_r')
-    axes[1, 1].set_ylabel('Error: ground truth − reconstruction\n\ngrid y (cells)')
-    fig.colorbar(im2, ax=axes[1, 1:], shrink=0.9,
-                 label='error\n(units of the truth standard deviation)')
-    for a, mm in zip(axes[2], COL_MASKS):
     a.scatter(*np.argwhere(mm).T, s=2, c='k')
     a.set_aspect(1); a.set_xlim(-1, NX); a.set_ylim(-1, NY)
-        a.set_title(f'{int(mm.sum())} of {N} samples', fontsize=10)
+    a.set_title(title, fontsize=10)
     a.set_xlabel('grid x (cells)')
-    axes[2, 0].set_ylabel('Sampling locations\n\ngrid y (cells)')
-    fig.suptitle(f'Reconstruction from {DELTA:.1%} of the samples: {note}')
-    plt.savefig(f'./{fname}.png', dpi=120)
 
 # --- figures 2a1/2b1 (kriging) and 2a2/2b2 (symlet compressive sensing): predicted standard
 # --- deviation versus actual error, one figure per ground truth and reconstruction method ---
@@ -459,14 +419,12 @@ def sigma_study(fname, desc, make_truth, truth, x_rec, sig_raw, rec_fn, rec_name
     for r in (0, 1):
         axes[r, 0].set_ylabel('grid y (cells)')
 
-    a = axes[2, 0]   # the sampling pattern that produced everything above
+    # the sampling pattern that produced everything above
-    a.scatter(*np.argwhere(mask_erg).T, s=2, c='k')
+    plot_mask_panel(axes[2, 0], mask_erg, f'Ergodic sampling locations\n({NS} of {N} samples)')
-    a.set_aspect(1); a.set_xlim(-1, NX); a.set_ylim(-1, NY)
+    axes[2, 0].set_ylabel('grid y (cells)')
-    a.set_title(f'Ergodic sampling locations\n({NS} of {N} samples)', fontsize=10)
-    a.set_xlabel('grid x (cells)'); a.set_ylabel('grid y (cells)')
 
     a = axes[2, 1]   # scatter with decile-binned means: does the prediction rank the error?
-    a.scatter(sig[OFF], aerr[OFF], s=2, alpha=.2, color='#9ecae1')
+    a.scatter(sig[OFF], aerr[OFF], s=2, alpha=.2, color=C_LIGHT_BLUE)
     bm, bv = reliability(sig[OFF], aerr[OFF])
     a.plot(bm, bv, 'o-', color='k', lw=1.6)
     xs = np.linspace(0, sig[OFF].max(), 50)
@@ -477,7 +435,7 @@ def sigma_study(fname, desc, make_truth, truth, x_rec, sig_raw, rec_fn, rec_name
     a.set_title('Prediction versus error at unsampled cells\n(single-realization Pearson '
                 f'correlation is capped\nat about {pc:.2f} even for a perfect prediction)',
                 fontsize=10)
-    a.legend(handles=[plt.Line2D([], [], marker='o', ls='none', ms=4, color='#9ecae1',
+    a.legend(handles=[plt.Line2D([], [], marker='o', ls='none', ms=4, color=C_LIGHT_BLUE,
                                  label='unsampled cells'),
                       plt.Line2D([], [], marker='o', ls='-', ms=5, color='k',
                                  label='mean absolute error\nwithin prediction decile'),
@@ -499,104 +457,10 @@ def sigma_study(fname, desc, make_truth, truth, x_rec, sig_raw, rec_fn, rec_name
                 'the prediction)', fontsize=10)
     a.legend(loc='upper left', fontsize=8)
     fig.suptitle('Statistical analysis of reconstruction error\n'
-                 f'({desc}; ergodic sampling; {rec_name} reconstruction)\n'
+                 f'({desc}; ergodic sampling; {rec_name} reconstruction)', fontsize=11)
-                ,fontsize=11)
     plt.savefig(f'./{fname}.png', dpi=120)
 
-DESC_SP = f'Fourier-sparse truth ({N_MODES} modes)'
+# ---------------- 8. sanity check: is the Gaussian process sigma statistically valid? ----------------
-DESC_PL = f'potential-field truth (spectral exponent {BETA}, source depth {Z_SRC} cells)'
-REC_GP, REC_SYM = 'kriging (Gaussian process)', 'compressive sensing (symlet wavelet, from paper)'
-KRIG_FAST_SP = gp_predictor(mask_erg, *GP_SP)
-KRIG_FAST_PL = gp_predictor(mask_erg, BETA_HAT, Z_HAT)
-sigma_study('plt2a1', DESC_SP, lambda: fourier_sparse_field(N_MODES), truth_sp, x_sp_gp, sig_sp_raw,
-            lambda y, m: gp_reconstruct(y, m, *GP_SP)[0], REC_GP, ens_rec_fn=KRIG_FAST_SP)
-sigma_study('plt2a2', DESC_SP, lambda: fourier_sparse_field(N_MODES), truth_sp, x_sp_sym, sig_sp_raw,
-            cs_reconstruct_sym, REC_SYM)
-sigma_study('plt2b1', DESC_PL, lambda: powerlaw_field(BETA), truth_pl, x_gp, sig_gp_raw,
-            lambda y, m: gp_reconstruct(y, m, BETA_HAT, Z_HAT)[0], REC_GP, corr=CORR_PL,
-            ens_rec_fn=KRIG_FAST_PL)
-sigma_study('plt2b2', DESC_PL, lambda: powerlaw_field(BETA), truth_pl, x_erg_sym, sig_gp_raw,
-            cs_reconstruct_sym, REC_SYM)
-
-# ---------------- Monte-Carlo RMS table over truth classes ----------------
-N_TRIALS = 25
-
-def white_noise_field():
-    # incompressible control: no method can beat predicting zero at unsampled cells,
-    # whose rms is sqrt(1 - DELTA) ~ 0.94
-    f = rng.standard_normal((NX, NY))
-    return (f - f.mean())/f.std()
-
-def trial_methods(truth, krig_sparse, krig_erg):
-    return {
-        'sparse+spline': spline_interp(truth, mask_reg),
-        'sparse+CS-Fourier': cs_reconstruct(truth*mask_reg, mask_reg),
-        'sparse+kriging': krig_sparse(truth*mask_reg),
-        'ergodic+spline': spline_interp(truth, mask_erg),
-        'ergodic+CS-Fourier': cs_reconstruct(truth*mask_erg, mask_erg),
-        'ergodic+CS-symlet': cs_reconstruct_sym(truth*mask_erg, mask_erg),
-        'ergodic+kriging': krig_erg(truth*mask_erg),
-    }
-
-CLASSES = [('white noise (control)', white_noise_field),
-           ('Fourier-sparse', lambda: fourier_sparse_field(N_MODES)),
-           ('potential-field', lambda: powerlaw_field(BETA))]
-print(f"\n--- Monte-Carlo rms error, mean +/- std over {N_TRIALS} trials per truth class ---")
-results = {}
-for cname, gen in CLASSES:
-    t0 = gen()
-    gp_params = gp_fit_params(t0*mask_erg, mask_erg)   # truth-free model selection, once per class
-    krig_sparse = gp_predictor(mask_reg, *gp_params)   # weights precomputed once per class
-    krig_erg = gp_predictor(mask_erg, *gp_params)
-    rr = {}
-    for k in range(N_TRIALS):
-        truth = t0 if k == 0 else gen()
-        for m, x in trial_methods(truth, krig_sparse, krig_erg).items():
-            rr.setdefault(m, []).append(rms(x, truth))
-    results[cname] = {m: (np.mean(v), np.std(v)) for m, v in rr.items()}
-print('method'.ljust(20) + ''.join(c.rjust(24) for c, _ in CLASSES))
-for m in results[CLASSES[0][0]]:
-    print(m.ljust(20) + ''.join(f"{results[c][m][0]:.3f} +/- {results[c][m][1]:.3f}".rjust(24)
-                                for c, _ in CLASSES))
-
-# --- figure 3: the Monte-Carlo table as a chart with error bars ---
-DISPLAY = {'sparse+spline': 'Sparse + spline interpolation',
-           'sparse+CS-Fourier': 'Sparse + compressive sensing (Fourier)',
-           'sparse+kriging': 'Sparse + kriging (Gaussian process)',
-           'ergodic+spline': 'Ergodic + spline interpolation',
-           'ergodic+CS-Fourier': 'Ergodic + compressive sensing (Fourier)',
-           'ergodic+CS-symlet': 'Ergodic + compressive sensing (symlet wavelet)',
-           'ergodic+kriging': 'Ergodic + kriging (Gaussian process)'}
-methods = list(DISPLAY)
-ypos = np.arange(len(methods))[::-1]
-xmax = max(results[c][m][0] + results[c][m][1] for c, _ in CLASSES for m in methods)
-fig, axes = plt.subplots(1, 3, figsize=(15.5, 5.8), sharey=True, constrained_layout=True)
-for a, (cname, _) in zip(axes, CLASSES):
-    mus = np.array([results[cname][m][0] for m in methods])
-    sds = np.array([results[cname][m][1] for m in methods])
-    # bar color encodes the sampling pattern (blue = sparse grid, orange = ergodic); the method
-    # labels already carry this, so no legend is needed
-    a.barh(ypos, mus, xerr=sds, height=0.6,
-           color=['#9ecae1' if m.startswith('sparse') else '#fdbe85' for m in methods],
-           error_kw=dict(ecolor='k', lw=1, capsize=3))
-    for yp, mu, sd in zip(ypos, mus, sds):
-        a.text(mu + sd + 0.02*xmax, yp, f'{mu:.3f} ± {sd:.3f}', va='center', fontsize=8)
-    t = cname[0].upper() + cname[1:] + ' truth'*(not cname.endswith(')'))
-    if cname == CLASSES[0][0]:
-        t += f'\n(dashed line: best possible = predict zero everywhere, {np.sqrt(1 - DELTA):.2f})'
-    a.set_title(t, fontsize=10)
-    a.set_xlabel('RMS error\n(units of the truth standard deviation)')
-    a.set_xlim(0, xmax*1.32)
-    a.grid(axis='x', lw=0.4, alpha=0.4)
-    a.set_axisbelow(True)
-axes[0].set_yticks(ypos, [DISPLAY[m] for m in methods])
-axes[0].axvline(np.sqrt(1 - DELTA), ls='--', lw=1, color='k')
-fig.suptitle('Monte-Carlo RMS reconstruction error: '
-             f'mean ± standard deviation over {N_TRIALS} trials per truth class '
-             f'({DELTA:.1%} of samples, error bars = ± one standard deviation)')
-plt.savefig('./plt3.png', dpi=120)
-
-# ---------------- sanity check: is the Gaussian process sigma statistically valid? ----------------
 # With the mask fixed, kriging is a linear map, so a large ensemble is cheap: precompute the
 # weight matrix once and run SANITY_K fresh truths through it. Under the hypothesis that the
 # posterior standard deviation is correct, the empirical standard deviation over K trials should
@@ -604,15 +468,13 @@ plt.savefig('./plt3.png', dpi=120)
 # cells are spatially correlated within a single realization, so the honest global test averages
 # the normalized squared error within each trial and treats the SANITY_K per-trial means as the
 # independent observations.
-RUN_SANITY_CHECK = False   # plt4 served its purpose; flip to True to re-run the check
+def run_sanity_check():
-if RUN_SANITY_CHECK:
     from scipy.stats import norm as normal_dist
 
     SANITY_K = 512
     Kss_s, Kxs_s = gp_matrices(mask_erg, BETA_HAT, Z_HAT)
     W_krig = Kxs_s @ np.linalg.inv(Kss_s)
-    var_unit = np.maximum(1.0 - np.einsum('ij,ji->i', Kxs_s, np.linalg.solve(Kss_s, Kxs_s.T)), 0)
+    sig_pred_unit = np.sqrt(gp_unit_variance(Kss_s, Kxs_s)).reshape(NX, NY)   # unit-variance field
-    sig_pred_unit = np.sqrt(var_unit).reshape(NX, NY)      # model posterior sigma, unit-variance field
     sanity_errs = np.empty((SANITY_K, NX, NY))
     for k in range(SANITY_K):
         t = powerlaw_field(BETA)
@@ -639,8 +501,7 @@ if RUN_SANITY_CHECK:
     # pipeline generator uses fixed-amplitude random phases and normalizes each realization by its
     # sample standard deviation, which perturbs the ensemble covariance away from the model; the
     # exact-model truths isolate the question "is the posterior itself constructed correctly?"
-    SQ_S = np.sqrt(np.where(KRAD > 0, KRAD, 1.0)**(-BETA_HAT)*np.exp(-4*np.pi*KRAD*Z_HAT))
+    SQ_S = np.sqrt(power_spectrum(BETA_HAT, Z_HAT))
-    SQ_S[0, 0] = 0.0
     C0_S = np.fft.ifft2(SQ_S**2).real[0, 0]
     ctrl_errs = np.empty((SANITY_K, NX, NY))
     for k in range(SANITY_K):
@@ -673,7 +534,7 @@ if RUN_SANITY_CHECK:
     fig.colorbar(imR, ax=axes[0], shrink=0.9,
                  label='difference\n(units of the truth standard deviation)')
     a = axes[1]
-    a.scatter(sig_pred_unit[OFF], sig_emp_big[OFF], s=2, alpha=.25, color='#9ecae1')
+    a.scatter(sig_pred_unit[OFF], sig_emp_big[OFF], s=2, alpha=.25, color=C_LIGHT_BLUE)
     lim = max(sig_pred_unit[OFF].max(), sig_emp_big[OFF].max())*1.05
     a.plot([0, lim], [0, lim], 'k--', lw=1, label='equality')
     a.set_xlim(0, lim); a.set_ylim(0, lim); a.set_aspect(1)
@@ -686,10 +547,10 @@ if RUN_SANITY_CHECK:
     zn = np.sqrt(2*SANITY_K)*(ratio - 1)          # approximately standard normal if sigma is correct
     zn_ctrl = np.sqrt(2*SANITY_K)*(ratio_ctrl - 1)
     zn_corr = np.sqrt(2*SANITY_K)*(ratio_corr - 1)
-    a.hist(zn, bins=50, density=True, color='#9ecae1', label='pipeline truths')
+    a.hist(zn, bins=50, density=True, color=C_LIGHT_BLUE, label='pipeline truths')
     a.hist(zn_corr, bins=50, density=True, histtype='step', color='#2ca02c', lw=1.5,
            label='pipeline truths,\na-priori corrected sigma')
-    a.hist(zn_ctrl, bins=50, density=True, histtype='step', color='#ff7f0e', lw=1.5,
+    a.hist(zn_ctrl, bins=50, density=True, histtype='step', color=C_ONE, lw=1.5,
            label='exact-model truths (control)')
     xs = np.linspace(min(zn_corr.min(), -4.5), max(zn.max(), 4.5), 200)
     a.plot(xs, normal_dist.pdf(xs), 'k-', lw=1.2, label='expected from sampling\nnoise alone')
@@ -706,3 +567,150 @@ if RUN_SANITY_CHECK:
                  f'(z = {z_corr:+.2f}, p = {p_corr:.2f}); truths drawn exactly from the Gaussian '
                  f'model also pass (z = {z_ctrl:+.2f}, p = {p_ctrl:.2f}).', fontsize=11)
     plt.savefig('./plt4.png', dpi=120)
+
+# ================ EXPERIMENT (everything below acquires, reconstructs, and plots) ================
+
+# --- ground truths and sampling patterns ---
+truth_pl = powerlaw_field(BETA)           # tier 2: hard, not l1-sparse in any basis
+truth_sp = fourier_sparse_field(N_MODES)  # tier 1: exactly sparse, the paper's Fig. A.1 setting
+
+mask_erg, f_erg = ergodic_mask()
+MASKS = {'sparse': mask_reg, 'ergodic': mask_erg}
+OFF = ~mask_erg      # evaluate z-scores at unsampled pixels only (err and sigma are both 0 at samples)
+f_rand = np.mean([objective(random_mask()) for _ in range(20)])
+print(f"N_samples={NS} (delta={DELTA:.3f})  obj: sparse={objective(mask_reg):.3f}  "
+      f"random(mean)={f_rand:.3f}  ergodic={f_erg:.3f}")
+
+# --- tier 1: Fourier-sparse truth: the CS claim of the paper (Fig. A.1 / Figs. 15, 22) ---
+print("\n--- tier 1: Fourier-sparse truth (paper's CS regime) ---")
+# kriging with the power-law covariance is misspecified for a 12-mode signal; shown as the
+# symmetric counterpart to CS-on-the-random-field so every process appears on both truths
+GP_SP = gp_fit_params(truth_sp*mask_erg, mask_erg)
+krig_sp = lambda y, m: gp_reconstruct(y, m, *GP_SP)[0]
+recs_sp = run_methods(truth_sp, krig_sp, krig_sp)
+sig_sp_raw = gp_reconstruct(truth_sp*mask_erg, mask_erg, *GP_SP)[1]
+for name, x in recs_sp.items():
+    print(f"{name:18s} leakage={leakage(x, truth_sp):+.3f}  rms={rms(x, truth_sp):.3f}")
+
+# --- tier 2: potential-field truth: not l1-sparse; report against the oracle floor ---
+print(f"\n--- tier 2: potential-field truth (beta={BETA}, source depth {Z_SRC}) ---")
+print(f"oracle {NS}-term Fourier floor: rms={rms(fourier_oracle(truth_pl), truth_pl):.3f}")
+BETA_HAT, Z_HAT = gp_fit_params(truth_pl*mask_erg, mask_erg)
+print(f"GP ML estimates: beta_hat={BETA_HAT:.1f} z_hat={Z_HAT:.1f} "
+      f"(generator beta={BETA}, z={Z_SRC})")
+CORR_PL = apriori_level_correction(BETA_HAT, Z_HAT, mask_erg)
+print(f"a-priori sigma level correction for the normalized truth generator: "
+      f"x{CORR_PL[~mask_erg].mean():.4f} (spread {CORR_PL[~mask_erg].std():.4f})")
+krig_pl = lambda y, m: gp_reconstruct(y, m, BETA_HAT, Z_HAT)[0]
+recs_pl = run_methods(truth_pl, krig_pl, krig_pl)
+sig_gp_raw = gp_reconstruct(truth_pl*mask_erg, mask_erg, BETA_HAT, Z_HAT)[1]
+print("method             leakage  rms    -rms/std(truth)   (leakage ~ -rms/std when error is unrecovered signal)")
+for name, x in recs_pl.items():
+    print(f"{name:18s} {leakage(x, truth_pl):+.3f}  {rms(x, truth_pl):.3f}   {-rms(x, truth_pl)/truth_pl.std():+.3f}")
+
+# --- figures 1a/1b: one figure per ground truth: reconstructions, errors, sampling locations ---
+COL_TITLES = ['Ground truth'] + [METHODS[k][0] for k in METHODS]
+COL_MASKS = [np.ones_like(mask_reg)] + [MASKS[k.split('+')[0]] for k in METHODS]
+for fname, gt, recs, note in [
+        ('plt1a', truth_sp, recs_sp,
+         f'Fourier-sparse truth ({N_MODES} modes); compressive sensing with the ergodic pattern '
+         'recovers the signal exactly'),
+        ('plt1b', truth_pl, recs_pl,
+         f'potential-field truth (spectral exponent {BETA}, source depth {Z_SRC} cells); '
+         'kriging is the best reconstruction for this field class')]:
+    imgs = [gt] + [recs[k] for k in METHODS]
+    fig, axes = plt.subplots(3, 9, figsize=(36, 12.4), constrained_layout=True)
+    v = np.abs(gt).max()
+    for c, (a, x, t) in enumerate(zip(axes[0], imgs, COL_TITLES)):
+        im = a.imshow(x.T, origin='lower', vmin=-v, vmax=v, cmap='RdBu_r')
+        sub = '' if c == 0 else (f'\nleakage {leakage(x, gt):+.2f} · '
+                                 f'RMS error {rms(x, gt):.3f}')
+        a.set_title(t + sub, fontsize=9)
+    axes[0, 0].set_ylabel('Reconstruction\n\ngrid y (cells)')
+    fig.colorbar(im, ax=axes[0], shrink=0.9,
+                 label='field value\n(units of the truth standard deviation)')
+    errors = [gt - x for x in imgs[1:]]
+    ve = max(np.abs(e).max() for e in errors)
+    axes[1, 0].axis('off')
+    for a, e in zip(axes[1, 1:], errors):
+        im2 = a.imshow(e.T, origin='lower', vmin=-ve, vmax=ve, cmap='RdBu_r')
+    axes[1, 1].set_ylabel('Error: ground truth − reconstruction\n\ngrid y (cells)')
+    fig.colorbar(im2, ax=axes[1, 1:], shrink=0.9,
+                 label='error\n(units of the truth standard deviation)')
+    for a, mm in zip(axes[2], COL_MASKS):
+        plot_mask_panel(a, mm, f'{int(mm.sum())} of {N} samples')
+    axes[2, 0].set_ylabel('Sampling locations\n\ngrid y (cells)')
+    fig.suptitle(f'Reconstruction from {DELTA:.1%} of the samples: {note}')
+    plt.savefig(f'./{fname}.png', dpi=120)
+
+# --- figures 2a/2b: predicted standard deviation versus actual error ---
+DESC_SP = f'Fourier-sparse truth ({N_MODES} modes)'
+DESC_PL = f'potential-field truth (spectral exponent {BETA}, source depth {Z_SRC} cells)'
+REC_GP, REC_SYM = 'kriging (Gaussian process)', 'compressive sensing (symlet wavelet, from paper)'
+KRIG_FAST_SP = gp_predictor(mask_erg, *GP_SP)
+KRIG_FAST_PL = gp_predictor(mask_erg, BETA_HAT, Z_HAT)
+sigma_study('plt2a1', DESC_SP, lambda: fourier_sparse_field(N_MODES), truth_sp,
+            recs_sp['ergodic+kriging'], sig_sp_raw, krig_sp, REC_GP, ens_rec_fn=KRIG_FAST_SP)
+sigma_study('plt2a2', DESC_SP, lambda: fourier_sparse_field(N_MODES), truth_sp,
+            recs_sp['ergodic+CS-symlet'], sig_sp_raw, cs_reconstruct_sym, REC_SYM)
+sigma_study('plt2b1', DESC_PL, lambda: powerlaw_field(BETA), truth_pl,
+            recs_pl['ergodic+kriging'], sig_gp_raw, krig_pl, REC_GP, corr=CORR_PL,
+            ens_rec_fn=KRIG_FAST_PL)
+sigma_study('plt2b2', DESC_PL, lambda: powerlaw_field(BETA), truth_pl,
+            recs_pl['ergodic+CS-symlet'], sig_gp_raw, cs_reconstruct_sym, REC_SYM)
+
+# --- Monte-Carlo RMS table over truth classes ---
+CLASSES = [('white noise (control)', white_noise_field),
+           ('Fourier-sparse', lambda: fourier_sparse_field(N_MODES)),
+           ('potential-field', lambda: powerlaw_field(BETA))]
+print(f"\n--- Monte-Carlo rms error, mean +/- std over {N_TRIALS} trials per truth class ---")
+results = {}
+for cname, gen in CLASSES:
+    t0 = gen()
+    gp_params = gp_fit_params(t0*mask_erg, mask_erg)   # truth-free model selection, once per class
+    krig_sparse = gp_predictor(mask_reg, *gp_params)   # weights precomputed once per class
+    krig_erg = gp_predictor(mask_erg, *gp_params)
+    rr = {}
+    for k in range(N_TRIALS):
+        truth = t0 if k == 0 else gen()
+        for m, x in run_methods(truth, krig_sparse, krig_erg, MC_KEYS).items():
+            rr.setdefault(m, []).append(rms(x, truth))
+    results[cname] = {m: (np.mean(v), np.std(v)) for m, v in rr.items()}
+print('method'.ljust(20) + ''.join(c.rjust(24) for c, _ in CLASSES))
+for m in MC_KEYS:
+    print(m.ljust(20) + ''.join(f"{results[c][m][0]:.3f} +/- {results[c][m][1]:.3f}".rjust(24)
+                                for c, _ in CLASSES))
+
+# --- figure 3: the Monte-Carlo table as a chart with error bars ---
+ypos = np.arange(len(MC_KEYS))[::-1]
+xmax = max(results[c][m][0] + results[c][m][1] for c, _ in CLASSES for m in MC_KEYS)
+fig, axes = plt.subplots(1, 3, figsize=(15.5, 5.8), sharey=True, constrained_layout=True)
+for a, (cname, _) in zip(axes, CLASSES):
+    mus = np.array([results[cname][m][0] for m in MC_KEYS])
+    sds = np.array([results[cname][m][1] for m in MC_KEYS])
+    # bar color encodes the sampling pattern (blue = sparse grid, orange = ergodic); the method
+    # labels already carry this, so no legend is needed
+    a.barh(ypos, mus, xerr=sds, height=0.6,
+           color=[C_LIGHT_BLUE if m.startswith('sparse') else C_LIGHT_ORANGE for m in MC_KEYS],
+           error_kw=dict(ecolor='k', lw=1, capsize=3))
+    for yp, mu, sd in zip(ypos, mus, sds):
+        a.text(mu + sd + 0.02*xmax, yp, f'{mu:.3f} ± {sd:.3f}', va='center', fontsize=8)
+    t = cname[0].upper() + cname[1:]
+    if not cname.endswith(')'):
+        t += ' truth'
+    if cname == CLASSES[0][0]:
+        t += f'\n(dashed line: best possible = predict zero everywhere, {np.sqrt(1 - DELTA):.2f})'
+    a.set_title(t, fontsize=10)
+    a.set_xlabel('RMS error\n(units of the truth standard deviation)')
+    a.set_xlim(0, xmax*1.32)
+    a.grid(axis='x', lw=0.4, alpha=0.4)
+    a.set_axisbelow(True)
+axes[0].set_yticks(ypos, [METHODS[m][1] for m in MC_KEYS])
+axes[0].axvline(np.sqrt(1 - DELTA), ls='--', lw=1, color='k')
+fig.suptitle('Monte-Carlo RMS reconstruction error: '
+             f'mean ± standard deviation over {N_TRIALS} trials per truth class '
+             f'({DELTA:.1%} of samples, error bars = ± one standard deviation)')
+plt.savefig('./plt3.png', dpi=120)
+
+if RUN_SANITY_CHECK:
+    run_sanity_check()

findings.md

@@ -2,15 +2,15 @@
 
 Written 2026-09-12. Refer to [ergodic_sampling_test.py](ergodic_sampling_test.py) for analysis and figures.
 
-## What was tested
+## What the analysis tested
 
-The analysis uses a 64 by 64 grid and keeps only 11.8% of the cells (484 samples), placed either on a
+The analysis found here: https://code.karsttech.com/jeremy/ergodic-sampling, uses a 64 by 64 grid and
-sparse grid (every third cell) or in the "ergodic" irregular pattern from the Zhang and Li paper
+keeps only 11.8% of the cells (484 samples), placed either on a sparse grid (every third cell) or in
-(found by minimizing equation 4). From those samples we rebuild the full grid with several methods
+the "ergodic" irregular pattern from the Ergodic Sampling paper (found by minimizing equation 4). 
-and measure how much the reconstruction differs from the ground truth. The value of merit chosen is
+From those samples we rebuild the full grid with several methods and measure how much the reconstruction
-RMS error: the typical size of the difference between the reconstruction and the true
+differs from the ground truth. The value of merit plotted is RMS error: the typical size of the difference
-field at a cell, normalized such that 1.0 is 1 standard deviation (lower is better). Each RMS error 
+between the reconstruction and the true field at a cell, and normalized to be scale invariant.
-value shown is the average over 25 iterations with a freshly generated ground truth.
+Each RMS error value shown is the average over 25 iterations with a freshly generated ground truth.
 
 The reconstruction methods:
 
@@ -19,69 +19,77 @@ The reconstruction methods:
 - **Compressive sensing**: the paper's method, iterative thresholding that assumes the signal is
   built from a small number of components in some transform. It is run two ways: with a Fourier
   transform (components are sinusoidal waves) and with the symlet wavelet from the paper.
-- **Kriging**: an interpolation method from geostatistics that the paper does not use; it was added
+- **Kriging**: the interpolation method from geostatistics, not used in the paper; it was added
-  as a strong reference. It estimates each empty cell as a weighted average of the samples, where
+  as a strong reference and a frequently used tool for modeling gravity and magnetic data. It 
-  the weights come from a model of how similar the field tends to be at each separation distance.
+  estimates each empty cell as a weighted average of the samples, where the weights come from a 
-  That similarity model is itself estimated from the samples, so kriging adapts to the data: for a
+  model of how similar the field tends to be at each separation distance. The similarity model 
-  smooth field it interpolates broadly, and for structureless data it learns that the samples say
+  is itself estimated from the samples, so kriging adapts to the data: for a smooth field it 
-  nothing about their neighbors and backs off.
+  interpolates broadly, and for structureless data it learns that the samples saynothing about 
+  their neighbors and backs off.
 
 The three kinds of ground truth:
 
 - **White noise (control)**: every cell is an independent random number. There is no structure at
-  all, so nothing between the samples can actually be predicted. This serves as a good control
+  all, so nothing between the samples can be predicted. This serves as a good control because it 
-  because it emphasizes any tendency of the reconstruction to hallucinate structure based on its
+  emphasizes any tendency of the reconstruction to hallucinate structure based on its
   assumptions about the data.
 - **Fourier-sparse**: a sum of 12 sinusoidal waves, half of them varying too quickly for the
   sparse grid to follow. This is exactly the kind of signal compressive sensing is designed for,
   and matches the paper's own demonstration signal.
 - **Potential-field**: a smooth random map that mimics the paper's gravity survey data.
 
-## The numbers
+## Some summarized results
 
 RMS error, mean plus or minus spread over 25 trials (best per column in bold):
 
 | method | white noise (control) | Fourier-sparse | potential-field |
 |---|---|---|---|
-| sparse + spline | 1.221 ± 0.023 | 0.517 ± 0.066 | 0.069 ± 0.012 |
+| sparse + spline | 1.218 ± 0.018 | 0.514 ± 0.093 | 0.068 ± 0.012 |
-| sparse + compressive sensing (Fourier) | 1.076 ± 0.011 | 0.155 ± 0.155 | 0.866 ± 0.100 |
+| sparse + compressive sensing (Fourier) | 1.076 ± 0.014 | 0.111 ± 0.132 | 0.812 ± 0.152 |
-| sparse + kriging | 1.013 ± 0.009 | 0.485 ± 0.065 | **0.060 ± 0.010** |
+| sparse + kriging | 1.014 ± 0.005 | 0.480 ± 0.085 | **0.059 ± 0.010** |
-| ergodic + spline | 1.445 ± 0.049 | 0.674 ± 0.078 | 0.166 ± 0.039 |
+| ergodic + spline | 1.443 ± 0.044 | 0.693 ± 0.116 | 0.161 ± 0.030 |
-| ergodic + compressive sensing (Fourier) | 1.056 ± 0.010 | **0.000 ± 0.000** | 0.252 ± 0.042 |
+| ergodic + compressive sensing (Fourier) | 1.063 ± 0.009 | **0.000 ± 0.000** | 0.251 ± 0.044 |
-| ergodic + compressive sensing (symlet) | 1.304 ± 0.034 | 0.737 ± 0.107 | 0.145 ± 0.020 |
+| ergodic + compressive sensing (symlet) | 1.282 ± 0.030 | 0.716 ± 0.113 | 0.138 ± 0.026 |
-| ergodic + kriging | **1.002 ± 0.007** | 0.558 ± 0.074 | 0.094 ± 0.016 |
+| ergodic + kriging | **1.001 ± 0.005** | 0.569 ± 0.091 | 0.095 ± 0.016 |
 
 ![Figure 1](plt3.png)
 
 *Figure 1. Reconstruction error broken down by sampling pattern + reconstruction formulation. Bar
-color encodes the sampling pattern: blue for sparse grid, orange for ergodic. Bars are the mean 
+color shows the sampling pattern with blue for sparse grid and orange for ergodic. Bars are the mean 
 over 25 trials per truth class. The dashed line in the left panel marks the best score attainable
 on structureless / random data, which is predicting zero at every unsampled cell (0.94).*
 
-## What was found
+## Some probably unsurprising results
 
 **1. On truly random ground truth, ergodic + kriging came out best of everything tested.**
-When the field has no structure, the smartest possible move is to admit it: output zero at every
+When the field has no structure, the smartest possible move is to admit it by outputting the mean
-unsampled cell, which scores 0.94 on our scale. Ergodic + kriging lands essentially on that mark
+(zero) at every unsampled cell, which scores 0.94 on our scale. Ergodic + kriging lands essentially
-(1.002) because its similarity model, estimated from the samples, correctly concludes the samples
+on that mark (1.001) because its similarity model is estimated from the samples, and correctly 
-carry no information about their neighbors, so it barely interpolates at all. Every method with a
+concludes the samples carry no information about their neighbors, so it barely interpolates at all.
-built-in belief about structure does worse than the "give up" answer: spline draws smooth hills
+As expected, every method with a built-in belief about structure does worse than the 
-that do not exist (1.22 on the sparse grid, 1.45 on the ergodic pattern) and the symlet version
+guess-zero-everywhere answer: spline draws smooth hills that do not exist (1.22 on the sparse grid, 
-paints in wavelet texture (1.30). That gap, up to roughly 50% worse than guessing zero, is the
+1.44 on the ergodic pattern) and the symlet version paints in wavelet texture (1.28). The 
-cost of hallucinated structure made concrete.
+misprediction gap of up to ~50% worse than guessing zero is the cost of hallucinated 
+structure made concrete.
 
-**2. On the potential-field ground truth, the plain sparse grid is the best pattern (figure 3).** This
+**2. On the potential-field ground truth, the plain sparse grid is the best pattern (figure 3).** 
-field is smooth enough that every third cell is dense enough sampling, the situation classical
+This field is smooth enough that every third cell is dense enough sampling, the situation classical
-sampling theory covers. There the sparse grid wins simply on coverage: its farthest cell from
+sampling theory covers. Here the sparse grid wins simply on coverage. The ergodic pattern, being 
-any sample is 1.41 cells away, while the ergodic pattern, being irregular, leaves gaps up to 3
+irregular, leaves larger gaps which costs about 1.6 times the error margin. That cost is the
-cells. The ergodic pattern costs about 1.5 times the error here. That cost is the
 premium for insurance that pays out on the Fourier-sparse truth, where half the signal varies too
 fast for the sparse grid: there the sparse grid garbles the fast waves into false slow ones
-(this is aliasing, and is destructive; compare the sparse-grid and ergodic columns of
+(this is traditional destructive aliasing, compare the sparse-grid and ergodic columns of
-figure 2). THIS CASE IS MOSTLY THEORETICAL: If we knew that there was no small-scale variation
+figure 2). THIS CASE IS MOSTLY THEORETICAL: sparse sampling's win here does not carry over to
-to capture, the sparse grid would be the right choice outright; in practice that is rarely known
+real applications for two reasons. First, the rough scale of variation is always a known prior
-before sampling, which is exactly what motivates the signal-agnostic ergodic pattern. See
+in practice (it is how any sample spacing gets chosen at all), but acting on that prior does not
-the "Potential Future Work" section for some related thoughts.
+lead to a comfortably dense regular grid: any scale at where we knew there were no spatial
+frequencies smaller than what uniform spacing resolves would be better served by sensing the same 
+area with fewer, ergodically spaced samples. The regular grid's win in this column is really a
+statement that its sample budget was larger than the known field scale required. Second, for
+potential fields in particular, variation much finer than the survey scale reflects shallow
+sources rather than the deep signal of interest, which is why this ground truth is generated
+with a source-depth attenuation term. See the "Potential Future Work" section for some related
+thoughts.
 
 ![Figure 2](plt1a.png)
 
@@ -89,7 +97,7 @@ the "Potential Future Work" section for some related thoughts.
 combinations. Top row: the ground truth on the left followed by each reconstruction, labelled with 
 its leakage and RMS error. Middle row: the signed error of each reconstruction, on a single shared color
 scale. Bottom row: the sampling locations that produced the column above. Ergodic sampling with Fourier 
-compressive sensing on a fourier sparse ground truth recovers the field exactly, leaving a blank error
+compressive sensing on a Fourier-sparse ground truth recovers the field exactly, leaving a blank error
 panel. Every sparse sampling method captures only the slow part of the field and leaves the unobserved
 high frequency error.*
 
@@ -97,7 +105,7 @@ high frequency error.*
 
 *Figure 3. The same layout for one realization of the potential-field ground truth. Because this field is
 smooth on the scale of the sample spacing, every method except Fourier compressive sensing reproduces it
-closely; Fourier compressive sensing insists on a handful of dominant waves hallucinates high frequency
+closely; Fourier compressive sensing insists on a handful of dominant waves and hallucinates high frequency
 noise, which shows up as the striped error panel in its column. What error remains for the other methods
 gathers in the widest gaps between samples, which is why the more even coverage of the sparse grid wins 
 here: its farthest cell from a sample is 1.41 cells away, against 3 cells for the ergodic pattern.*
@@ -107,17 +115,17 @@ Each column of the table is won by the method whose built-in assumption mirrors
 kriging on the smooth random fields (its fitted similarity model actually recovers the
 generator's parameters exactly), Fourier compressive sensing on the sum-of-waves signal, and
 kriging again on noise because it alone can gracefully back off to predicting zero. Mismatches
-fail hard: Fourier compressive sensing assumes a few dominant waves and scores 0.866 on the
+are dramatically worse: Fourier compressive sensing assumes a few dominant waves and scores 0.812 
-smooth field, and the symlet version, which assumes the wrong kind of building block for a
+on the smooth field, and the symlet version, which assumes the wrong kind of building block for a
-sum-of-waves signal, reaches only 0.737 where the Fourier version is exact.
+sum-of-waves signal, reaches only 0.716 where the Fourier version is exact.
 
 ## Predicting the error without the ground truth
 
-Everything above scores a reconstruction against a truth we happen to know. In a real survey the truth
+Everything above scores a reconstruction against a truth we happen to know, but in a real survey the 
-is exactly what is missing, so the practical question is whether the reconstruction can say, by itself,
+the practical question is whether the reconstruction can say, by itself, where it is likely to be wrong. 
-where it is likely to be wrong. Kriging answers that for free: alongside an estimate at every cell it
+Kriging answers that for free: alongside an estimate at every cell it reports the standard deviation of 
-reports the standard deviation of that estimate, computed from the sample locations and from a
+that estimate, computed from the sample locations and from a similarity model fitted to the sample values 
-similarity model fitted to the sample values alone. No ground truth enters at any point.
+alone. No ground truth enters at any point.
 
 Figures 4 to 7 test that predicted error map for both ground truth classes and for both the kriging and
 the symlet compressive sensing reconstruction. Each figure has the same layout. The top row shows the
@@ -126,15 +134,22 @@ standard deviation beside two things it should resemble: the error actually meas
 truths drawn from the same generator and sampled with the same pattern, and the absolute error in the
 single realization above. The bottom row shows the sampling pattern, a cell-by-cell scatter of predicted
 against actual error, and a calibration curve in which cells are binned into deciles of the prediction.
+Fixed axes are chosen so that the results are directly comparable from figure to figure, not just subplots.
 
 Two cautions matter when reading these. First, one realization of an error is one draw of a random
 variable, so even a perfect prediction correlates with it only up to a ceiling, which is printed on each
-scatter panel; the comparison against the 64-truth ensemble is the decisive test, and there the
+scatter panel. The comparison against the 64-truth ensemble is a test showing the prediction tracks the 
-prediction tracks the actual error pattern with correlations of 0.91 to 0.99. Second, uniform ergodic
+actual error pattern with correlations of 0.90 to 0.99. Second, uniform ergodic coverage intentionally 
-coverage intentionally flattens the predicted map, since the whole point of the pattern is to leave no
+flattens the predicted map, since the whole point of the pattern is to leave no cell poorly covered, and 
-cell poorly covered, and that narrow range is what pushes the single-realization ceiling down. The
+that narrow range is what pushes the single-realization ceiling down. The calibration curves run close to 
-calibration curves run close to the one-to-one line in every case, slightly below it, meaning the
+the one-to-one line in every case, slightly below it, meaning the predicted magnitudes are roughly correct.
-predicted magnitudes are roughly correct and not merely correctly ranked.
+**One takeaway here** is not the magnitude of the error mismatch between the single realization and the 
+predicted error map, but it's structure. If the error map is dominated by high frequency variations on the 
+scale of a discovered anomaly like an ore deposit, then smaller variations would be expected to average 
+out over the span of the contiguous region of the anomaly and thus a more precise estimate of the size and 
+extent of the anomaly could be gathered than by a simply measuring the reconstruction. This is because we 
+can leverage the structure of the ergodic sampling points and their implied error structure to reign in 
+the uncertainty in our anomaly predictions.
 
 ![Figure 4](plt2a1.png)
 
@@ -170,7 +185,7 @@ worse than predicting zero. The paper's ergodic pattern is best understood in th
 deliberately designed without any knowledge of the signal, as insurance that keeps every option open.
 Alongside it, our error-prediction experiment (figures 4 to 7) shows that a map of the expected 
 error at every cell can be computed from nothing but the sample values and their locations, and 
-it tracks the actual error pattern closely (correlation 0.91 to 0.99 against the error measured over 64
+it tracks the actual error pattern closely (correlation 0.90 to 0.99 against the error measured over 64
 independent trials, with correctly calibrated magnitudes).
 
 ## Does kriging reconstruction strictly dominate symlet wavelet reconstruction?
@@ -182,7 +197,13 @@ violates compressive sensing's incoherence requirement. Two caveats: for signals
 wavelet-like anomalies (the paper's Figure 15 regime) symlet keeps a relative edge, though at
 11.8% sampling both methods there do worse than predicting zero; and for sparse-spectrum signals
 Fourier compressive sensing, not kriging, is the right choice (exact recovery kriging cannot
-match). Symlet was never the best method in any regime tested.
+match). Symlet was never the best with any of my ground-truth construction methods.
+
+The potential-field truths are generated from the same stationary covariance family that kriging
+fits, so it's no coincidence that kriging is the mathematically optimal reconstruction; its dominance
+is a statement about matched assumptions and generalizability, not about kriging being universally
+superior. This is tangential to the main purpose of the paper (chosen ergodic sample spacing),
+and I only mention it because kriging seems to be more applicable to magnetic and gravity sensing.
 
 ## Potential future work