Fixed a couple minor issues, improved performance and wording of findings section.

ergodic_sampling_test.py

@@ -9,14 +9,15 @@ NX = NY = 64; N = NX*NY
 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)
-STRIDE = 3                      # regular sparse stride -> delta ~ 0.118
+STRIDE = 3                      # sparse (regular-grid) stride -> delta ~ 0.118
 WAVELET, LEVELS, MODE = 'sym8', 2, 'periodization'   # 2 = pywt.dwt_max_level(64, 'sym8')
-SYM_SPINS = 16                  # cycle-spin shifts for the symlet reconstruction
+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
 LOO_M = 60                      # leave-one-out calibration points
-ENS_K = 128                     # ensemble validation realizations
+ENS_K = 64                      # ensemble validation realizations (our construct, not the paper's)
 
 # ---------------- ground truths (two tiers) ----------------
 def powerlaw_field(beta, z_src=Z_SRC):
@@ -30,7 +31,7 @@ def powerlaw_field(beta, z_src=Z_SRC):
 
 def fourier_sparse_field(n_modes):
     # visually structured sparse signal: amplitude ~ |k|^-1/2 so low frequencies dominate the look,
-    # with half the modes beyond the regular grid's Nyquist 1/(2*STRIDE) so that grid aliases
+    # with half the modes beyond the sparse grid's Nyquist 1/(2*STRIDE) so that grid aliases
     k_alias = 1/(2*STRIDE)
     C = np.zeros((NX, NY), complex); picked = set()
     m = 0
@@ -118,7 +119,7 @@ def ergodic_mask(iters=SA_ITERS, T0=0.05):
 
 mask_erg, f_erg = ergodic_mask()
 f_rand = np.mean([objective(random_mask()) for _ in range(20)])
-print(f"N_samples={NS} (delta={DELTA:.3f})  obj: regular={objective(mask_reg):.3f}  "
+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 ----------------
@@ -182,10 +183,14 @@ def gp_matrices(mask, beta, z_src):
 
 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
+    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):
-            Kss, _ = gp_matrices(mask, beta, z_src)
+            cov = 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:
@@ -195,6 +200,13 @@ def gp_fit_params(y_img, mask):
             if best is None or ll > best[0]: best = (ll, beta, z_src)
     return best[1], best[2]
 
+def gp_predictor(mask, beta, z_src):
+    # precomputed kriging weights: with the mask and model fixed, reconstruction is a single
+    # matrix multiplication per field (the second argument exists only for signature parity)
+    Kss, Kxs = gp_matrices(mask, beta, z_src)
+    Wm = Kxs @ np.linalg.inv(Kss)
+    return lambda y_img, m=None: (Wm @ y_img[mask]).reshape(NX, NY)
+
 def gp_reconstruct(y_img, mask, beta, z_src):
     yv = y_img[mask]
     Kss, Kxs = gp_matrices(mask, beta, z_src)
@@ -205,6 +217,41 @@ def gp_reconstruct(y_img, mask, beta, z_src):
     s2 = yv @ np.linalg.solve(Kss, yv)/len(yv)          # profiled variance scale
     return xh, np.sqrt(var*s2).reshape(NX, NY)
 
+def apriori_level_correction(beta, z_src, mask, n_top=400, m_draws=20000):
+    # The potential-field generator uses fixed-amplitude random phases and divides each
+    # realization by its sample standard deviation. For a linear reconstruction the resulting
+    # error variance is E[X/Y] with X = sum_p t_p u_p (per-mode error contributions) and
+    # 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 = (Sfit/Sfit.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)
+    s_pair = np.where(idx == neg, Sn, 2*Sn)[take]
+    kxp, kyp = ix[take], iy[take]
+    o = np.argsort(s_pair)[::-1]
+    s_pair, kxp, kyp = s_pair[o], kxp[o], kyp[o]
+    rs = np.random.default_rng(123)                     # local stream, leaves rng untouched
+    u_top = 1 + np.cos(2*np.pi*rs.random((m_draws, n_top)))
+    Y = u_top @ s_pair[:n_top]
+    for j0 in range(n_top, len(s_pair), 512):
+        blk = s_pair[j0:j0+512]
+        Y += (1 + np.cos(2*np.pi*rs.random((m_draws, len(blk))))) @ blk
+    g_top = (u_top/Y[:, None]).mean(0)                  # E[u_p / Y] for the dominant modes
+    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)
+    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
+    t_top = (np.abs(resp)**2)*s_pair[:n_top]/2
+    t_tail = np.maximum(sig2_u - t_top.sum(1), 0)
+    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 ----------------
 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))
@@ -228,15 +275,16 @@ 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 [('regular+spline', x_sp_spl), ('regular+CS', x_sp_reg),
+for name, x in [('sparse+spline', x_sp_spl), ('sparse+CS', x_sp_reg),
-                ('regular+kriging', x_sp_reg_gp), ('ergodic+spline', x_sp_erg_spl),
+                ('sparse+CS-symlet', x_sp_reg_sym), ('sparse+kriging', x_sp_reg_gp),
-                ('ergodic+CS', x_sp_erg), ('ergodic+CS-symlet', x_sp_sym),
+                ('ergodic+spline', x_sp_erg_spl), ('ergodic+CS', x_sp_erg),
-                ('ergodic+kriging', x_sp_gp)]:
+                ('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
@@ -245,18 +293,23 @@ print(f"oracle {NS}-term Fourier floor: rms={rms(fourier_oracle(truth_pl), truth
 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 [('regular+spline', x_pl_spl),
+for name, x in [('sparse+spline', x_pl_spl),
-                ('regular+CS', x_pl_reg),
+                ('sparse+CS', x_pl_reg),
-                ('regular+kriging', x_pl_reg_gp),
+                ('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),
@@ -299,44 +352,31 @@ def rms_reduce(a): return np.sqrt(np.mean(a**2))
 
 # ---------------- plots (no abbreviations in any figure text) ----------------
 C_ENS, C_ONE = '#1f77b4', '#ff7f0e'   # blue = ensemble comparison, orange = single realization
-OBJ_REG = objective(mask_reg)
 
-# --- figure 1: the three sampling patterns (paper Figures 3/14 analog) ---
+# --- figures 1a/1b: one figure per ground truth: reconstructions, errors, sampling locations ---
-fig, axes = plt.subplots(1, 3, figsize=(13.5, 4.8), constrained_layout=True)
-for a, m, t in zip(axes,
-                   [np.ones_like(mask_reg), mask_reg, mask_erg],
-                   [f'Dense reference grid\n{N} samples, objective = 0 by definition',
-                    f'Regular sparse subset\n{NS} samples, objective = {OBJ_REG:.1f}',
-                    f'Ergodic subset, optimized by equation 4\n{NS} samples, objective = {f_erg:.1f}']):
-    a.scatter(*np.argwhere(m).T, s=3, c='k')
-    a.set_aspect(1); a.set_xlim(-1, NX); a.set_ylim(-1, NY)
-    a.set_title(t); a.set_xlabel('grid x (cells)')
-axes[0].set_ylabel('grid y (cells)')
-fig.suptitle(f'Sampling patterns on the {NX}×{NY} grid: each sparse pattern keeps {DELTA:.1%} of the '
-             f'samples (mean objective of 20 purely random patterns: {f_rand:.1f}; lower is better)')
-plt.savefig('./plt1.png', dpi=120)
-
-# --- figures 2a/2b: one figure per ground truth: reconstructions (top) and errors (bottom) ---
 COLS = ['Ground truth',
-        'Regular + spline interpolation\n(naive control)',
+        'Sparse + spline interpolation\n(naive control)',
-        'Regular + compressive sensing\n(Fourier)',
+        'Sparse + compressive sensing\n(Fourier)',
-        'Regular + kriging\n(Gaussian process)',
+        'Sparse + compressive sensing\n(symlet wavelet, as in the paper)',
+        'Sparse + kriging\n(Gaussian process)',
         'Ergodic + spline interpolation',
         'Ergodic + compressive sensing\n(Fourier)',
         'Ergodic + compressive sensing\n(symlet wavelet, as in the paper)',
         'Ergodic + kriging\n(Gaussian process)']
-COL_MASKS = [np.ones_like(mask_reg), mask_reg, mask_reg, mask_reg,
+COL_MASKS = [np.ones_like(mask_reg), mask_reg, mask_reg, mask_reg, mask_reg,
              mask_erg, mask_erg, mask_erg, mask_erg]
 for fname, gt, imgs, note in [
-        ('plt2a', truth_sp,
+        ('plt1a', truth_sp,
-         [truth_sp, x_sp_spl, x_sp_reg, x_sp_reg_gp, x_sp_erg_spl, x_sp_erg, x_sp_sym, x_sp_gp],
+         [truth_sp, x_sp_spl, x_sp_reg, x_sp_reg_sym, x_sp_reg_gp,
+          x_sp_erg_spl, x_sp_erg, x_sp_sym, x_sp_gp],
          f'Fourier-sparse truth ({N_MODES} modes); compressive sensing with the ergodic pattern '
          'recovers the signal exactly'),
-        ('plt2b', truth_pl,
+        ('plt1b', truth_pl,
-         [truth_pl, x_pl_spl, x_pl_reg, x_pl_reg_gp, x_pl_erg_spl, x_erg, x_erg_sym, x_gp],
+         [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, 8, figsize=(32, 12.4), constrained_layout=True)
+    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')
@@ -363,21 +403,27 @@ for fname, gt, imgs, note in [
     fig.suptitle(f'Reconstruction from {DELTA:.1%} of the samples: {note}')
     plt.savefig(f'./{fname}.png', dpi=120)
 
-# --- figures 3a/3b: predicted standard deviation versus actual error, one per ground truth ---
+# --- figures 2a1/2b1 (kriging) and 2a2/2b2 (symlet compressive sensing): predicted standard
-def sigma_study(fname, desc, make_truth, truth, x_rec, sig_raw):
+# --- deviation versus actual error, one figure per ground truth and reconstruction method ---
+def sigma_study(fname, desc, make_truth, truth, x_rec, sig_raw, rec_fn, rec_name, corr=None,
+                ens_rec_fn=None):
     y_img = truth*mask_erg
-    loo = loo_rms(y_img, mask_erg, cs_reconstruct)
+    loo = loo_rms(y_img, mask_erg, rec_fn)
-    sig = recalibrate(sig_raw, loo)    # GP shape transfers to the l1 solver; LOO sets the level
+    if corr is None:
+        sig = recalibrate(sig_raw, loo); level_note = 'sigma level set by leave-one-out'
+    else:
+        sig = sig_raw*corr; level_note = 'sigma level corrected a priori'
     err = x_rec - truth; aerr = np.abs(err)
     print(f"\n--- sigma study ({fname}): {desc} ---")
-    print(f"LOO error rms (truth-free calibration level): {loo:.3f}   "
+    print(f"LOO error rms (truth-free check): {loo:.3f}   "
-          f"actual rms: {np.sqrt((err**2).mean()):.3f}")
+          f"actual rms: {np.sqrt((err**2).mean()):.3f}   ({level_note})")
     pr = pearsonr(aerr.ravel(), sig.ravel())[0]; sr = spearmanr(aerr.ravel(), sig.ravel())[0]
     pc, sc = pearson_ceiling(sig), spearman_ceiling(sig)
     print(f"|err| vs sigma_GP: pearson={pr:+.3f} (ceiling {pc:.3f}, ratio {pr/pc:+.2f})  "
           f"spearman={sr:+.3f} (ceiling {sc:.3f}, ratio {sr/sc:+.2f})  "
           f"z-std={np.std(err[OFF]/sig[OFF]):.2f}")
-    ens_errs = np.array([cs_reconstruct(t*mask_erg, mask_erg) - t
+    ens_rec = ens_rec_fn if ens_rec_fn is not None else rec_fn
+    ens_errs = np.array([ens_rec(t*mask_erg, mask_erg) - t
                          for t in (make_truth() for _ in range(ENS_K))])
     sig_emp = np.sqrt((ens_errs**2).mean(0))
     pr2 = pearsonr(sig_emp.ravel(), sig.ravel())[0]; sr2 = spearmanr(sig_emp.ravel(), sig.ravel())[0]
@@ -386,27 +432,26 @@ def sigma_study(fname, desc, make_truth, truth, x_rec, sig_raw):
           f"pooled z-std={np.std(ens_errs[:, OFF]/sig[OFF]):.2f}")
 
     fig, axes = plt.subplots(3, 3, figsize=(14.5, 13.8), constrained_layout=True)
-    vs = max(sig.max(), sig_emp.max(), aerr.max())
+    # fixed axes (3.0 truth standard deviations) so the plt2 figures are directly comparable;
-    v = np.abs(truth).max()
+    # some values clip by design
+    v = 3.0
     im0 = axes[0, 0].imshow(truth.T, origin='lower', vmin=-v, vmax=v, cmap='RdBu_r')
     axes[0, 0].set_title('Ground truth', fontsize=10)
     axes[0, 1].imshow(x_rec.T, origin='lower', vmin=-v, vmax=v, cmap='RdBu_r')
-    axes[0, 1].set_title('Reconstruction: ergodic + compressive sensing (Fourier)\n'
+    axes[0, 1].set_title(f'Reconstruction: ergodic + {rec_name}\n'
                          f'leakage {leakage(x_rec, truth):+.2f} · root-mean-square error '
                          f'{rms(x_rec, truth):.3f}', fontsize=10)
-    fig.colorbar(im0, ax=axes[0, :2], shrink=0.9,
+    axes[0, 2].imshow(err.T, origin='lower', vmin=-v, vmax=v, cmap='RdBu_r')
-                 label='field value\n(units of the truth standard deviation)')
-    imE = axes[0, 2].imshow(err.T, origin='lower', vmin=-vs, vmax=vs, cmap='RdBu_r')
     axes[0, 2].set_title('Signed error:\nreconstruction − truth', fontsize=10)
-    fig.colorbar(imE, ax=axes[0, 2], shrink=0.9,
+    fig.colorbar(im0, ax=axes[0], shrink=0.9,
-                 label='error\n(units of the truth standard deviation)')
+                 label='field value or error\n(units of the truth standard deviation)')
     for a, x, t in zip(axes[1],
                        [sig, sig_emp, aerr],
                        ['Predicted standard deviation:\nGaussian process posterior (truth-free)',
                         f'Empirical standard deviation\n({ENS_K} independent truths, '
                         'same sampling pattern)',
                         'Actual absolute error, this realization:\n|reconstruction − truth|']):
-        im1 = a.imshow(x.T, origin='lower', vmin=0, vmax=vs, cmap='magma')
+        im1 = a.imshow(x.T, origin='lower', vmin=0, vmax=v, cmap='turbo')
         a.set_title(t, fontsize=10)
         a.set_xlabel('grid x (cells)')
     fig.colorbar(im1, ax=axes[1], shrink=0.9,
@@ -427,6 +472,7 @@ def sigma_study(fname, desc, make_truth, truth, x_rec, sig_raw):
     xs = np.linspace(0, sig[OFF].max(), 50)
     a.plot(xs, np.sqrt(2/np.pi)*xs, '--', color='k', lw=1.2)
     a.set_xlabel('predicted standard deviation'); a.set_ylabel('absolute error, this realization')
+    a.set_ylim(0, v)
     a.ticklabel_format(style='sci', scilimits=(-2, 3))
     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)',
@@ -444,9 +490,8 @@ def sigma_study(fname, desc, make_truth, truth, x_rec, sig_raw):
                             (sig_emp, f'against the {ENS_K}-truth ensemble', C_ENS)]:
         mm, rr = reliability(sig, target, reduce=rms_reduce)
         a.plot(mm, rr, 'o-', color=cc, label=lab)
-    cal_lim = max(a.get_xlim()[1], a.get_ylim()[1])
+    a.plot([0, 1], [0, 1], 'k--', lw=.8, label='perfect calibration')
-    a.plot([0, cal_lim], [0, cal_lim], 'k--', lw=.8, label='perfect calibration')
+    a.set_xlim(0, 1.0); a.set_ylim(0, 1.0); a.set_aspect(1)
-    a.set_xlim(0, cal_lim); a.set_ylim(0, cal_lim); a.set_aspect(1)
     a.set_xlabel('mean predicted standard deviation within decile')
     a.set_ylabel('root-mean-square error or\nempirical standard deviation within decile')
     a.ticklabel_format(style='sci', scilimits=(-2, 3))
@@ -454,16 +499,26 @@ def sigma_study(fname, desc, make_truth, truth, x_rec, sig_raw):
                 'the prediction)', fontsize=10)
     a.legend(loc='upper left', fontsize=8)
     fig.suptitle('Truth-free predicted standard deviation versus actual reconstruction error\n'
-                 f'({desc}; ergodic sampling; compressive-sensing reconstruction)\n'
+                 f'({desc}; ergodic sampling; {rec_name} reconstruction)\n'
                  'Uniform ergodic coverage intentionally flattens the prediction; its narrow '
                  'range caps single-realization correlation;\nthe ensemble comparison is the '
                  'decisive test.', fontsize=11)
     plt.savefig(f'./{fname}.png', dpi=120)
 
-sigma_study('plt3a', f'Fourier-sparse truth ({N_MODES} modes)',
+DESC_SP = f'Fourier-sparse truth ({N_MODES} modes)'
-            lambda: fourier_sparse_field(N_MODES), truth_sp, x_sp_erg, sig_sp_raw)
+DESC_PL = f'potential-field truth (spectral exponent {BETA}, source depth {Z_SRC} cells)'
-sigma_study('plt3b', f'potential-field truth (spectral exponent {BETA}, source depth {Z_SRC} cells)',
+REC_GP, REC_SYM = 'kriging (Gaussian process)', 'compressive sensing (symlet wavelet, as in the paper)'
-            lambda: powerlaw_field(BETA), truth_pl, x_erg, sig_gp_raw)
+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
@@ -474,15 +529,15 @@ def white_noise_field():
     f = rng.standard_normal((NX, NY))
     return (f - f.mean())/f.std()
 
-def trial_methods(truth, gp_params):
+def trial_methods(truth, krig_sparse, krig_erg):
     return {
-        'regular+spline': spline_interp(truth, mask_reg),
+        'sparse+spline': spline_interp(truth, mask_reg),
-        'regular+CS-Fourier': cs_reconstruct(truth*mask_reg, mask_reg),
+        'sparse+CS-Fourier': cs_reconstruct(truth*mask_reg, mask_reg),
-        'regular+kriging': gp_reconstruct(truth*mask_reg, mask_reg, *gp_params)[0],
+        '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': gp_reconstruct(truth*mask_erg, mask_erg, *gp_params)[0],
+        'ergodic+kriging': krig_erg(truth*mask_erg),
     }
 
 CLASSES = [('white noise (control)', white_noise_field),
@@ -493,10 +548,12 @@ 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, gp_params).items():
+        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))
@@ -504,10 +561,10 @@ 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 4: the Monte-Carlo table as a chart with error bars ---
+# --- figure 3: the Monte-Carlo table as a chart with error bars ---
-DISPLAY = {'regular+spline': 'Regular + spline interpolation',
+DISPLAY = {'sparse+spline': 'Sparse + spline interpolation',
-           'regular+CS-Fourier': 'Regular + compressive sensing (Fourier)',
+           'sparse+CS-Fourier': 'Sparse + compressive sensing (Fourier)',
-           'regular+kriging': 'Regular + kriging (Gaussian process)',
+           '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)',
@@ -519,10 +576,10 @@ fig, axes = plt.subplots(1, 3, figsize=(15.5, 5.8), sharey=True, constrained_lay
 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 = regular, orange = ergodic); the method
+    # 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('regular') else '#fdbe85' for m in methods],
+           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)
@@ -539,4 +596,115 @@ axes[0].axvline(np.sqrt(1 - DELTA), ls='--', lw=1, color='k')
 fig.suptitle('Monte-Carlo root-mean-square 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
+# differ from it only by sampling noise (relative size about 1/sqrt(2K) per cell). The unsampled
+# 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
+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(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)
+        sanity_errs[k] = (W_krig @ t[mask_erg]).reshape(NX, NY) - t
+    sig_emp_big = np.sqrt((sanity_errs**2).mean(0))
+    resid = sig_emp_big - sig_pred_unit
+
+    # per-trial mean of error^2/sigma^2 over unsampled cells: one independent number per trial
+    m_k = (sanity_errs[:, OFF]**2/np.maximum(sig_pred_unit[OFF]**2, 1e-12)).mean(1)
+    z_glob = (m_k.mean() - 1)/(m_k.std(ddof=1)/np.sqrt(SANITY_K))
+    p_glob = 2*normal_dist.sf(abs(z_glob))
+    ratio = sig_emp_big[OFF]/sig_pred_unit[OFF]
+    print(f"\n--- sanity check: Gaussian process sigma vs {SANITY_K}-trial empirical sigma ---")
+    print(f"mean of per-trial normalized squared error: {m_k.mean():.4f} (1.0 if sigma is correct)")
+    print(f"global test over {SANITY_K} independent trials: z = {z_glob:+.2f}, "
+          f"two-sided p = {p_glob:.3f}")
+    print(f"per-cell ratio sigma_empirical/sigma_predicted at {OFF.sum()} unsampled cells: "
+          f"mean {ratio.mean():.4f}, spread {ratio.std():.4f} "
+          f"(pure sampling noise predicts spread ~ {1/np.sqrt(2*SANITY_K):.4f})")
+    corr_sanity = pearsonr(sig_pred_unit[OFF], sig_emp_big[OFF])[0]
+    print(f"correlation(sigma_predicted, sigma_empirical) = {corr_sanity:.4f}")
+
+    # control: truths drawn exactly from the fitted Gaussian model with a deterministic scale. The
+    # 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[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):
+        t = np.fft.ifft2(SQ_S*np.fft.fft2(rng.standard_normal((NX, NY)))).real/np.sqrt(C0_S)
+        ctrl_errs[k] = (W_krig @ t[mask_erg]).reshape(NX, NY) - t
+    m2_k = (ctrl_errs[:, OFF]**2/np.maximum(sig_pred_unit[OFF]**2, 1e-12)).mean(1)
+    z_ctrl = (m2_k.mean() - 1)/(m2_k.std(ddof=1)/np.sqrt(SANITY_K))
+    p_ctrl = 2*normal_dist.sf(abs(z_ctrl))
+    ratio_ctrl = np.sqrt((ctrl_errs**2).mean(0))[OFF]/sig_pred_unit[OFF]
+    print(f"control with exact-model truths (no per-realization normalization): "
+          f"mean ratio {m2_k.mean():.4f}, z = {z_ctrl:+.2f}, two-sided p = {p_ctrl:.3f}")
+
+    # same pipeline ensemble scored against the a-priori corrected sigma
+    sig_corr_unit = sig_pred_unit*CORR_PL
+    m3_k = (sanity_errs[:, OFF]**2/np.maximum(sig_corr_unit[OFF]**2, 1e-12)).mean(1)
+    z_corr = (m3_k.mean() - 1)/(m3_k.std(ddof=1)/np.sqrt(SANITY_K))
+    p_corr = 2*normal_dist.sf(abs(z_corr))
+    ratio_corr = np.sqrt((sanity_errs**2).mean(0))[OFF]/sig_corr_unit[OFF]
+    print(f"pipeline truths with the a-priori level correction: "
+          f"mean ratio {m3_k.mean():.4f}, z = {z_corr:+.2f}, two-sided p = {p_corr:.3f}")
+
+    fig, axes = plt.subplots(1, 3, figsize=(16.5, 4.9), constrained_layout=True)
+    vr = np.abs(resid).max()
+    imR = axes[0].imshow(np.where(OFF, resid, np.nan).T, origin='lower', vmin=-vr, vmax=vr,
+                         cmap='RdBu_r')
+    axes[0].set_title(f'Empirical ({SANITY_K} trials) − predicted standard deviation\n'
+                      '(sampled cells blanked; structure here would mean a construction error)',
+                      fontsize=10)
+    axes[0].set_xlabel('grid x (cells)'); axes[0].set_ylabel('grid y (cells)')
+    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')
+    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)
+    a.set_xlabel('predicted standard deviation')
+    a.set_ylabel(f'empirical standard deviation ({SANITY_K} trials)')
+    a.set_title('Per-cell agreement at unsampled cells\n'
+                f'(correlation {corr_sanity:.3f})', fontsize=10)
+    a.legend(loc='upper left', fontsize=8)
+    a = axes[2]
+    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_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,
+           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')
+    a.set_xlabel('normalized per-cell deviation of the ratio\nempirical/predicted standard deviation')
+    a.set_ylabel('probability density')
+    a.set_title(f'Deviations vs pure sampling noise\n(pipeline: z = {z_glob:+.2f}; corrected: '
+                f'z = {z_corr:+.2f}; control: z = {z_ctrl:+.2f})', fontsize=10)
+    a.legend(loc='upper right', fontsize=8)
+    fig.suptitle('Validity check of the Gaussian process standard deviation over '
+                 f'{SANITY_K} independent truths: shape agreement {corr_sanity:.3f}.\n'
+                 'The level offset for pipeline truths comes from their per-realization '
+                 'normalization, not from the posterior: it is computable a priori from the '
+                 'spectrum and mask,\nand correcting for it centers the test '
+                 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)

findings.md

@@ -4,8 +4,8 @@ Written 2026-09-12. Refer to [ergodic_sampling_test.py](ergodic_sampling_test.py
 
 ## What was tested
 
-We work on a 64 by 64 grid and keep only 11.8% of the cells (484 samples), placed either on a
+The analysis uses a 64 by 64 grid and keeps only 11.8% of the cells (484 samples), placed either on a
-regular 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 Zhang and Li 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
 root-mean-squared error: the typical size of the difference between the reconstruction and the true
@@ -30,10 +30,10 @@ 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
-  because it emphasizes any tendancy of the reconstruction to hallucinate structure based on it's
+  because it 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
-  regular grid to follow. This is exactly the kind of signal compressive sensing is designed for,
+  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.
 
@@ -43,75 +43,81 @@ Root-mean-square error, mean plus or minus spread over 25 trials (best per colum
 
 | method | white noise (control) | Fourier-sparse | potential-field |
 |---|---|---|---|
-| regular + spline | 1.220 ± 0.014 | 0.503 ± 0.074 | 0.067 ± 0.011 |
+| sparse + spline | 1.221 ± 0.023 | 0.517 ± 0.066 | 0.069 ± 0.012 |
-| regular + compressive sensing (Fourier) | 1.077 ± 0.010 | 0.117 ± 0.134 | 0.848 ± 0.140 |
+| sparse + compressive sensing (Fourier) | 1.076 ± 0.011 | 0.155 ± 0.155 | 0.866 ± 0.100 |
-| regular + kriging | 1.014 ± 0.006 | 0.475 ± 0.072 | **0.058 ± 0.010** |
+| sparse + kriging | 1.013 ± 0.009 | 0.485 ± 0.065 | **0.060 ± 0.010** |
-| ergodic + spline | 1.417 ± 0.035 | 0.689 ± 0.107 | 0.150 ± 0.023 |
+| ergodic + spline | 1.445 ± 0.049 | 0.674 ± 0.078 | 0.166 ± 0.039 |
-| ergodic + compressive sensing (Fourier) | 1.059 ± 0.008 | **0.000 ± 0.000** | 0.242 ± 0.046 |
+| ergodic + compressive sensing (Fourier) | 1.056 ± 0.010 | **0.000 ± 0.000** | 0.252 ± 0.042 |
-| ergodic + compressive sensing (symlet) | 1.286 ± 0.031 | 0.725 ± 0.124 | 0.136 ± 0.024 |
+| ergodic + compressive sensing (symlet) | 1.304 ± 0.034 | 0.737 ± 0.107 | 0.145 ± 0.020 |
-| ergodic + kriging | **1.000 ± 0.007** | 0.563 ± 0.086 | 0.092 ± 0.016 |
+| ergodic + kriging | **1.002 ± 0.007** | 0.558 ± 0.074 | 0.094 ± 0.016 |
 
 ## What was found
 
-**1. On truly random ground truth, ergodic + kriging came out best of everything we tested.**
+**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.000) because its similarity model, estimated from the samples, correctly concludes the samples
+(1.002) 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 regular grid, 1.42 on the ergodic pattern) and the symlet version
+that do not exist (1.22 on the sparse grid, 1.45 on the ergodic pattern) and the symlet version
-paints in wavelet texture (1.29). That gap, up to 50% worse than guessing zero, is the cost of
+paints in wavelet texture (1.30). That gap, up to roughly 50% worse than guessing zero, is the
-hallucinated structure made concrete.
+cost of hallucinated structure made concrete.
 
-**2. On the potential-field ground truth, the plain regular grid is the best pattern.** This
+**2. On the potential-field ground truth, the plain sparse grid is the best pattern.** This
 field is smooth enough that every third cell is dense enough sampling, the situation classical
-sampling theory covers. There the regular grid wins simply on coverage: its farthest cell from
+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
 premium for insurance that pays out on the Fourier-sparse truth, where half the signal varies too
-fast for the regular grid: there the regular grid garbles the fast waves into false slow ones
+fast for the sparse grid: there the sparse grid garbles the fast waves into false slow ones
 (this is aliasing, and no processing can undo it). THIS CASE IS MOSTLY THEORETICAL: If we knew
-that there were no small-scale variation to capture with our sparse sampling method
+that there were no small-scale variation to capture, the sparse grid would be the right choice
+outright; in practice that is rarely known before sampling, which is exactly what motivates the
+signal-agnostic ergodic pattern.
 
 **3. Matching the reconstruction method to how the ground truth was generated wins, as expected.**
 Each column of the table is won by the method whose built-in assumption mirrors the generator: 
 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.848 on the
+fail hard: Fourier compressive sensing assumes a few dominant waves and scores 0.866 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.725 where the Fourier version is exact.
+sum-of-waves signal, reaches only 0.737 where the Fourier version is exact.
 
 ## The bigger picture
 
 Choosing how to reconstruct in essence is a choice about how much you trust what you know about
-the underlying signall. A correct guess about signal structure rewards all the way up to perfect
+the underlying signal. A correct guess about signal structure rewards all the way up to perfect
 recovery with minimal samples when the assumption is exactly right. The same choices punish wrong
 assumptions: a method that expects structure will manufacture it out of nothing, and on our 
 structureless control every such method lost to simply predicting zero, with kriging only slightly
 worse than predicting zero. The paper's ergodic pattern is best understood in this light: it is 
 deliberately designed without any knowledge of the signal, as insurance that keeps every option open.
-Alongside it, our error-prediction experiment (plt3a and plt3b) shows that a map of the expected 
+Alongside it, our error-prediction experiment (the plt2 figures) 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.85 to 0.88 against the error measured 
+it tracks the actual error pattern closely (correlation 0.91 to 0.99 against the error measured over 64
-over 128 independent trials).
+independent trials, with correctly calibrated magnitudes).
 
 ## Potential future work
 
-Since we have shown that for cases were no small scale variation exists to capture, sparse sampling
+Since we have shown that for cases where no small scale variation exists to capture, sparse sampling
 wins due to optimal coverage; we know that in such a case Ergodic sampling would still be superior 
 if we knew at what distance scale local variation existed and further reduced our sample count.
 
 I propose the following:
 
-Develop a sampling procedure which optimizes for collection cost for several collection strategies
+Develop a sampling procedure which optimizes for collection cost in two collection scenarios
  - Where travel is the dominant cost
    - This applies to survey type collection where a vehicle carrying a sensor is used
    - Compute an optimal survey path which leverages Ergodic sensing methods but reframed as a
-     densely sampled path; chosing the path dynamically as sensed signal spatial scales are discovered.
+     densely sampled path; choosing the path dynamically as sensed signal spatial scales are discovered.
 
  - Where number of sample locations are the dominant cost
-   - This applies to sattelite-pointing type collection or ground-station collection
+   - This applies to satellite-pointing type collection or ground-station collection
    - Compute an optimal survey path which leverages Ergodic sensing methods, but is scale-adaptive.
      Rather than having coverage sparsity as a prior, this would discover the required coverage
      dynamically to attain a certain reconstruction confidence metric.
+
+Both scenarios should ideally attempt to iteratively estimate a confidence level based on how much future
+predictions match the current model (created from previous samples), and therefore determine when sampling
+should stop due to being sufficient for our confidence interval goals.