Reorganized the code and added some comments to the findings document.

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

@@ -5,7 +5,7 @@ Written 2026-09-12. Refer to [ergodic_sampling_test.py](ergodic_sampling_test.py
 ## What was tested
 
 The analysis uses a 64 by 64 grid and keeps only 11.8% of the cells (484 samples), placed either on a
-sparse grid (every third cell) or in the "ergodic" irregular pattern from the Zhang and Li paper
+sparse grid (every third cell) or in the "ergodic" irregular pattern from the Ergodic Sampling paper
 (found by minimizing equation 4). From those samples we rebuild the full grid with several methods
 and measure how much the reconstruction differs from the ground truth. The value of merit chosen is
 RMS error: the typical size of the difference between the reconstruction and the true
@@ -24,7 +24,8 @@ The reconstruction methods:
   the weights come from a model of how similar the field tends to be at each separation distance.
   That similarity model is itself estimated from the samples, so kriging adapts to the data: for a
   smooth field it interpolates broadly, and for structureless data it learns that the samples say
-  nothing about their neighbors and backs off.
+  nothing about their neighbors and backs off. It is also the standard gridding tool for gravity
+  and magnetic data in industry practice.
 
 The three kinds of ground truth:
 
@@ -43,13 +44,13 @@ 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)
 
@@ -63,18 +64,18 @@ on structureless / random data, which is predicting zero at every unsampled cell
 **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
 unsampled cell, which scores 0.94 on our scale. Ergodic + kriging lands essentially on that mark
-(1.002) because its similarity model, estimated from the samples, correctly concludes the samples
+(1.001) because its similarity model, estimated from the samples, correctly concludes the samples
 carry no information about their neighbors, so it barely interpolates at all. Every method with a
 built-in belief about structure does worse than the "give up" answer: spline draws smooth hills
-that do not exist (1.22 on the sparse grid, 1.45 on the ergodic pattern) and the symlet version
+that do not exist (1.22 on the sparse grid, 1.44 on the ergodic pattern) and the symlet version
-paints in wavelet texture (1.30). That gap, up to roughly 50% worse than guessing zero, is the
+paints in wavelet texture (1.28). That gap, up to roughly 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
 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
 any sample is 1.41 cells away, while the ergodic pattern, being irregular, leaves gaps up to 3
-cells. The ergodic pattern costs about 1.5 times the error here. That cost is the
+cells. The ergodic pattern costs about 1.6 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
@@ -89,7 +90,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 +98,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,9 +108,9 @@ 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
+fail hard: Fourier compressive sensing assumes a few dominant waves and scores 0.812 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
 
@@ -130,7 +131,7 @@ against actual error, and a calibration curve in which cells are binned into dec
 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
-prediction tracks the actual error pattern with correlations of 0.91 to 0.99. Second, uniform ergodic
+prediction tracks the actual error pattern with correlations of 0.90 to 0.99. Second, uniform ergodic
 coverage intentionally flattens the predicted map, since the whole point of the pattern is to leave no
 cell poorly covered, and that narrow range is what pushes the single-realization ceiling down. The
 calibration curves run close to the one-to-one line in every case, slightly below it, meaning the
@@ -170,7 +171,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 +183,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