#!/usr/bin/env python3 """ PROVENANCE: PROOF Computation 191 -- The Cauchy slab on a deterministic cloud: the slab-wise CAR limit needs equidistribution with a rate, not independence, and that is A3's kind of condition ========================================================================================== THE ITEM. comp 186 built the construction the slice-wise CAR limit was missing: no site is realised ON a measure-zero Cauchy slice, so the relations are posed on a SLAB of thickness h_n, whose density nu_n = lambda_n h_n carries the [mass]^3 the withdrawn delta^(4) reading lacked, and the slab sum converges to the slice integral when h_n -> 0 and lambda_n h_n -> oo. It measured that on the I.I.D. SURROGATE, where the second condition is the central limit theorem, and said so: "for the emergence-determined rho_n the equidistribution that part 3 leans on is exactly A3". So the item rested on a surrogate. This file removes the surrogate. It proves that the slab limit holds for EVERY sequence of clouds whose discrepancy vanishes faster than the thickness, with no probability anywhere, and measures that on three deterministic clouds, where it converges FASTER than on the surrogate, because independence was never the hypothesis. THE DECOMPOSITION, exact. Let P be any N points in the unit 4-box, t the time coordinate, F = the smearing integrand, and S_h(P) = (N h)^-1 sum_{a in P : |t_a - 1/2| <= h/2} F(p_a) the slab estimator. Write I(h) = integral of F over the slab. Then, as an identity, S_h(P) - integral_Sigma F = [ I(h)/h - integral_Sigma F ] + (1/h) [ N^-1 sum_slab F - I(h) ] = BIAS(h) + DISC(P, h) / h. The bias depends on h alone, not on the cloud: it is the slab average minus the slice value, O(h^2) for smooth F (here exactly 3h^2/4 + h^4/20). The second term is the cloud's equidistribution defect tested on F times the slab indicator. THEOREM (Koksma--Hlawka, applied). For a cloud of star discrepancy D_N and F of bounded Hardy--Krause variation, |N^-1 sum F 1_slab - I(h)| <= V(F 1_slab) D_N, and the indicator of a box multiplies the variation by at most a constant, so | S_h(P) - integral_Sigma F | <= c_F h^2 + C_F D_N / h. Hence the slab-wise CAR limit holds for ANY sequence of clouds with D_N -> 0, once h_n -> 0 and D_N / h_n -> 0, and the balanced thickness is h_n ~ D_N^{1/3} with total error D_N^{2/3}. No independence, no measure, no averaging over clouds: one cloud, one inequality. The i.i.d. surrogate is the special case D_N ~ N^{-1/2} (with probability one), for which this worst-case bound gives h ~ N^{-1/6} and error N^{-1/3}, weaker than the CLT's N^{-2/5}; that gap is exactly what independence was buying, a probabilistic improvement on a deterministic theorem that never needed it. A low-discrepancy cloud has D_N ~ (log N)^4 / N, hence h ~ N^{-1/3} and error N^{-2/3} up to logarithms. THE SAME THEOREM IN A3's TEST CLASS. The sharp slab tests the cloud on a discontinuous function. Replace the indicator by a smooth window w_h(t) of width h and unit mass; the estimator N^-1 sum F w_h tests the cloud on a C^1 function with ||F w_h||_{H^1} of order h^{-3/2}, so an H^1-tested weak discrepancy eps_N, the form in which the Mosco chain's rate condition A3 is stated, gives error <= c h^2 + eps_N h^{-3/2}, balanced at h ~ eps_N^{2/7}. Measured in part 4: the windowed estimator converges on every deterministic cloud, and its bias is exactly the window's second moment, 9h^2/7. WHAT THIS CLOSES, AND WHAT IT DOES NOT. The slab-wise CAR limit is a theorem given (i) a deterministic equidistribution rate for the realised cloud, tested on the smearing functions, and (ii) the time function tau in which the slab is cut. (i) is the kind of statement A3 makes for the emergence-determined rho_n and nothing more than it: the rate P3's maximum-entropy sublimation supplies, which equicoercivity A7 already rests on. (ii) is P3's. So the item no longer rests on a surrogate; it rests on the same two things the Mosco chain rests on, and W4 follows from it (comp 185). It is NOT closed here, because A3 is supplied by P3 and not derived, and because the slab needs the rate to VANISH, which is what "rate hypothesis" means, where coercivity needs only its leading coefficient inside the margin d_0^2/4. Part 5's tamper shows the count A1 is not the hypothesis, and part 6's shows what the time function must deliver. CHECKS. 1. The decomposition is an identity and the bias formula is exact: the slab integral I(h) = 2h + 3h^3/4 + h^5/20 against quadrature to 1e-8. 2. Three deterministic clouds (a jittered grid with one fixed stratification, the Halton sequence in bases 2, 3, 5, 7, and the Kronecker sequence of a badly approximable vector) and the i.i.d. surrogate beside them, at |D_n| = 2^12 .. 2^20 with h = |D_n|^{-1/3}. A single deterministic cloud's error oscillates in sign and size, so the RMS over eight DISJOINT BLOCKS of each sequence is reported, an envelope with no probability in it; the two low-discrepancy clouds' errors fall below the surrogate's at every size and their fitted rates approach the predicted -2/3. The jittered grid is the lesson in the other direction: it is i.i.d. below its cell scale N^{-1/4}, the slab is thinner than that, so it tracks the surrogate. Stratification at one scale is not equidistribution at every scale, and the theorem asks for D_N / h -> 0. 3. The discrepancy term alone, DISC/h, against its bound: at fixed |D_n| it grows as 1/h on the low-discrepancy clouds and as (|D_n| h)^{-1/2} on the surrogate; at fixed h it falls as |D_n|^{-1} and |D_n|^{-1/2} respectively. RMS over eight blocks or draws throughout. 4. The smooth window: bias 9h^2/7 exactly, and convergence on every cloud. 5. Tamper, A1 without equidistribution: a Halton cloud with the same count but a spatial void refilled elsewhere converges to the WRONG slice value, with a defect flat across a 16-fold count increase and equal to the analytic void bias. The count is not the condition. 6. Tamper, the time function: a cloud equidistributed in space but with time density 1 + eps cos(2 pi t) returns (1 - eps) times the slice value. tau must be the function in which the cloud is 4-volume equidistributed, which is what P3 supplies. A "pass": the identity and bias exact; the Halton and Kronecker errors below the surrogate's at every size, their rates at most -0.55 and -0.45, and the jittered grid still converging; the discrepancy-term exponents within 0.3 (deterministic) and 0.2 (surrogate) of the predictions; the windowed bias within 5% of 9h^2/7 and the windowed error falling monotonically on the two low-discrepancy clouds; the void defect flat within 20% across a 16-fold count increase and within 10% of the analytic value; the time tamper within 2% of 1 - eps. """ import math import time import numpy as np CHUNK = 1 << 17 # points generated at a time, for the browser's heap TARGET = 2.0 # the slice integral of F at t = 1/2 over the unit cube ALPHA = np.sqrt(np.array([2.0, 3.0, 5.0, 7.0])) # Kronecker direction def F(t, x): """comp 186's smearing integrand , even in t - 1/2 so the bias does not cancel.""" d = t - 0.5 f = 1.0 + 4.0 * d * d + np.sin(2 * np.pi * x[:, 0]) g = 2.0 + 1.0 * d * d + np.cos(2 * np.pi * x[:, 1]) return f * g def I_slab(h): """integral of F over the slab |t - 1/2| <= h/2, exact: the x-integrals kill every sin and cos, leaving int (1 + 4d^2)(2 + d^2) dd = 2h + 3h^3/4 + h^5/20.""" return 2 * h + 0.75 * h ** 3 + h ** 5 / 20 def bias(h): return I_slab(h) / h - TARGET # --- clouds, generated in chunks --------------------------------------------------- def radical_inverse(idx, base): x = np.zeros(idx.shape); f = 1.0 / base; i = idx.copy() while i.any(): x += f * (i % base); i //= base; f /= base return x _HALTON_CACHE = {} def halton(lo, hi): # The sequence is deterministic, and the same block is asked for by every slab width and # every window, so each block is generated once (since v32.162; the values are unchanged). key = (lo, hi) if key not in _HALTON_CACHE: if len(_HALTON_CACHE) > 64: _HALTON_CACHE.clear() idx = np.arange(lo + 1, hi + 1, dtype=np.int64) # the sequence from its start, skipping 0 _HALTON_CACHE[key] = np.column_stack([radical_inverse(idx, b) for b in (2, 3, 5, 7)]) return _HALTON_CACHE[key] def kronecker(lo, hi): n = np.arange(lo + 1, hi + 1, dtype=np.float64)[:, None] return np.mod(n * ALPHA[None, :] + 0.5, 1.0) def jittered(lo, hi, m, rng): """cell (i,j,k,l) of an m^4 grid holds one point, uniform inside the cell: the stratification is fixed, only the position inside each cell is drawn. Below the cell scale 1/m = N^{-1/4} this cloud is i.i.d. again, so a slab thinner than a cell sees the surrogate's fluctuation; it is included for exactly that lesson.""" idx = np.arange(lo, hi, dtype=np.int64) cells = np.column_stack([(idx // m ** k) % m for k in range(4)]).astype(np.float64) return (cells + rng.random((hi - lo, 4))) / m def cloud_chunks(kind, N, seed=0): """For the two sequences, `seed` selects the DISJOINT BLOCK [seed N, (seed+1) N) of the sequence: eight blocks are eight deterministic clouds, and their RMS error is an envelope of the oscillating single-cloud error, with no probability in it. For the surrogate and the jittered grid, `seed` seeds the draw.""" rng = np.random.default_rng(seed) m = int(round(N ** 0.25)) lo = 0 while lo < N: hi = min(lo + CHUNK, N) if kind == "iid": yield rng.random((hi - lo, 4)) elif kind == "halton": yield halton(seed * N + lo, seed * N + hi) elif kind == "kronecker": yield kronecker(seed * N + lo, seed * N + hi) elif kind == "jittered": yield jittered(lo, hi, m, rng) lo = hi def rms_error(kind, N, h, blocks=8, window=False): v = [slab_sum(kind, N, h, seed=b, window=window) / (1 if window else h) - TARGET for b in range(blocks)] return math.sqrt(float(np.mean(np.square(v)))) def slab_sum(kind, N, h, seed=0, window=False, warp=None): """N^-1 sum over the slab of F (sharp), or N^-1 sum of F w_h (smooth window).""" total = 0.0 for p in cloud_chunks(kind, N, seed): t, x = p[:, 0], p[:, 1:] if warp is not None: t = warp(t) d = t - 0.5 if window: m = np.abs(d) <= h if m.any(): u = d[m] / h total += float(np.sum(F(t[m], x[m]) * (15 / 16) / h * (1 - u * u) ** 2)) else: m = np.abs(d) <= h / 2 if m.any(): total += float(np.sum(F(t[m], x[m]))) return total / N def estimate(kind, N, h, seed=0): return slab_sum(kind, N, h, seed) / h def fit_slope(xs, ys): return float(np.polyfit(np.log(xs), np.log(np.abs(ys)), 1)[0]) # --- parts -------------------------------------------------------------------------- def part1(): print("\n (1) THE DECOMPOSITION IS AN IDENTITY AND THE BIAS IS EXACT") ok = True for h in (0.4, 0.1, 0.02): d = np.linspace(-h / 2, h / 2, 20001) y = (1 + 4 * d * d) * (2 + d * d) quad = float(np.sum(0.5 * (y[1:] + y[:-1]) * np.diff(d))) # trapezoid rule, written out for the browser's numpy ok &= abs(quad - I_slab(h)) < 1e-8 print(f" h = {h:5.2f}: I(h) quadrature {quad:.12f}, closed form {I_slab(h):.12f}, bias {bias(h):.6f}") # identity on one concrete cloud N, h = 1 << 16, 0.1 s = slab_sum("halton", N, h) lhs = s / h - TARGET; rhs = bias(h) + (s - I_slab(h)) / h ok &= abs(lhs - rhs) < 1e-12 print(f" identity on a Halton cloud, N = {N}, h = {h}: S - target = {lhs:.6e} = bias + disc/h = {rhs:.6e}") return ok def part2(): print("\n (2) EIGHT DISJOINT BLOCKS OF EACH SEQUENCE, RMS ERROR, h = N^(-1/3): against the surrogate") kinds = ("jittered", "halton", "kronecker", "iid") Ns = [2 ** k for k in range(12, 21)] err = {k: [] for k in kinds} print(f" {'N':>9} {'h':>8}" + "".join(f"{k:>12}" for k in kinds)) for N in Ns: h = N ** (-1 / 3) row = [] for k in kinds: e = rms_error(k, N, h); err[k].append(e); row.append(e) print(f" {N:9d} {h:8.4f}" + "".join(f"{e:12.2e}" for e in row)) slopes = {k: fit_slope(Ns, err[k]) for k in kinds} print(" fitted rate " + "".join(f"{slopes[k]:12.2f}" for k in kinds) + " (predicted -2/3 low-discrepancy)") below = all(abs(err[k][i]) < abs(err["iid"][i]) for k in ("halton", "kronecker") for i in range(len(Ns))) rates = slopes["halton"] <= -0.55 and slopes["kronecker"] <= -0.45 jit = abs(err["jittered"][-1]) < abs(err["jittered"][0]) print(f" Halton and Kronecker below the surrogate at every size: {below}; rates <= -0.55 and -0.45: {rates}") print(f" the jittered grid, i.i.d. below its cell scale N^(-1/4) > h, tracks the surrogate and still converges: {jit}") return below and rates and jit def part3(): print("\n (3) THE DISCREPANCY TERM ALONE, DISC/h, AGAINST ITS BOUND") ok = True def disc(k, N, h, seeds=8): # RMS over eight blocks of the sequence, or eight draws of the surrogate: one # cloud's discrepancy oscillates in h and N and carries no readable exponent return math.sqrt(np.mean([((slab_sum(k, N, h, seed=s) - I_slab(h)) / h) ** 2 for s in range(seeds)])) N = 1 << 18 hs = [0.4 * 2 ** (-k / 2) for k in range(8)] print(f" at N = {N}, versus h: " + "".join(f"{h:>9.3f}" for h in hs) + " slope predicted") pred = {"halton": -1.0, "kronecker": -1.0, "iid": -0.5} tol = {"halton": 0.3, "kronecker": 0.3, "iid": 0.2} for k in ("halton", "kronecker", "iid"): v = [disc(k, N, h) for h in hs] s = fit_slope(hs, v); ok &= abs(s - pred[k]) < tol[k] print(f" {k:>10}: " + "".join(f"{x:>9.1e}" for x in v) + f" {s:6.2f} {pred[k]:6.2f}") h = 0.1 Ns = [2 ** k for k in range(12, 21)] print(f" at h = {h}, versus N: " + "".join(f"{n:>9d}" for n in Ns) + " slope predicted") for k in ("halton", "kronecker", "iid"): v = [disc(k, n, h) for n in Ns] s = fit_slope(Ns, v); ok &= abs(s - pred[k]) < tol[k] print(f" {k:>10}: " + "".join(f"{x:>9.1e}" for x in v) + f" {s:6.2f} {pred[k]:6.2f}") print(" (the surrogate's slopes are the CLT's; the deterministic ones are the discrepancy's, 1/N up to logarithms)") return ok def part4(): print("\n (4) THE SMOOTH WINDOW, A3's TEST CLASS: bias 9h^2/7 exactly, convergence on every cloud") ok = True # the window's bias, from a huge Halton cloud so the discrepancy term is negligible print(f" {'h':>7} {'windowed - target':>18} {'9h^2/7 + 4h^4/21':>18} {'ratio':>7}") for h in (0.2, 0.1, 0.05): w = slab_sum("halton", 1 << 20, h, window=True) - TARGET pb = 9 * h * h / 7 + 4 * h ** 4 / 21 ok &= abs(w / pb - 1) < 0.05 print(f" {h:7.3f} {w:18.6f} {pb:18.6f} {w/pb:7.3f}") Ns = [4 ** k for k in range(6, 11)] print(f" {'N':>9} {'h':>8}" + "".join(f"{k:>12}" for k in ("jittered", "halton", "kronecker", "iid"))) errs = {k: [] for k in ("jittered", "halton", "kronecker", "iid")} for N in Ns: h = N ** (-2 / 7) row = [] for k in errs: e = rms_error(k, N, h, window=True); errs[k].append(e); row.append(e) print(f" {N:9d} {h:8.4f}" + "".join(f"{e:12.2e}" for e in row)) conv = all(abs(errs[k][-1]) < abs(errs[k][0]) for k in errs) mono = all(all(abs(errs[k][i + 1]) < abs(errs[k][i]) for i in range(len(Ns) - 1)) for k in ("halton", "kronecker")) print(f" every cloud's windowed error falls across the range: {conv}; monotonically on Halton and Kronecker: {mono}") return ok and conv and mono def void_refill(N, r): """Halton cloud with the points inside the spatial ball |x - 1/2| < r removed and replaced by later Halton points outside it: same count, a void, and a density 1/(1 - v) elsewhere. A1 holds, equidistribution does not.""" keep, lo = [], 0 have = 0 while have < N: p = halton(lo, lo + CHUNK); lo += CHUNK m = np.sum((p[:, 1:] - 0.5) ** 2, axis=1) >= r * r p = p[m]; keep.append(p); have += len(p) return np.vstack(keep)[:N] def part5(): print("\n (5) TAMPER, A1 WITHOUT EQUIDISTRIBUTION: a void refilled elsewhere, same count") r = 0.3; v = 4 / 3 * np.pi * r ** 3 # analytic limit: (1/(1 - v)) * int_{cube \ ball} F(1/2, x) d^3x, by 3D Halton quadrature q = halton(0, 1 << 20)[:, 1:] outside = np.sum((q - 0.5) ** 2, axis=1) >= r * r limit = float(np.mean(F(np.full(outside.sum(), 0.5), q[outside]))) * outside.mean() / (1 - v) print(f" ball radius {r}, void fraction v = {v:.4f}; analytic limit of the estimator = {limit:.5f} against the slice value {TARGET}") print(f" {'N':>9} {'h':>8} {'estimate':>10} {'defect':>10} {'analytic':>10}") defects = [] for N in (4 ** 8, 4 ** 9, 4 ** 10): h = N ** (-1 / 3) P = void_refill(N, r) t, x = P[:, 0], P[:, 1:] m = np.abs(t - 0.5) <= h / 2 est = float(np.sum(F(t[m], x[m]))) / (N * h) defects.append(est - TARGET) print(f" {N:9d} {h:8.4f} {est:10.5f} {est - TARGET:10.5f} {limit - TARGET:10.5f}") flat = max(defects) / min(defects) < 1.2 if min(defects) > 0 else False close = abs(defects[-1] / (limit - TARGET) - 1) < 0.10 print(f" defect flat across a 16-fold count increase: {flat}; final defect within 10% of the analytic void bias: {close}") return flat and close def part6(): print("\n (6) TAMPER, THE TIME FUNCTION: time density 1 + eps cos(2 pi t), space equidistributed") eps = 0.3 def warp(u): # invert the CDF T(t) = t + eps sin(2 pi t) / (2 pi) by Newton, so t has density 1 + eps cos(2 pi t) t = u.copy() for _ in range(30): t -= (t + eps * np.sin(2 * np.pi * t) / (2 * np.pi) - u) / (1 + eps * np.cos(2 * np.pi * t)) return t N = 1 << 20; h = N ** (-1 / 3) est = slab_sum("halton", N, h, warp=warp) / h print(f" N = {N}, h = {h:.4f}: estimate {est:.5f}, slice value {TARGET}, ratio {est / TARGET:.4f} against 1 - eps = {1 - eps}") ok = abs(est / TARGET - (1 - eps)) < 0.02 print(f" the slab reads the time marginal's density at the slice: {ok}") return ok def main(): t0 = time.time() print("=" * 100) print(" Computation 191 -- The Cauchy slab on a deterministic cloud: equidistribution, not independence") print("=" * 100) checks = [ ("the decomposition is an identity and the bias is exact", part1()), ("deterministic clouds beat the surrogate, Halton and Kronecker at the predicted rate", part2()), ("the discrepancy term scales as its bound says, in h and in N", part3()), ("the smooth window: bias 9h^2/7 and convergence on every cloud", part4()), ("A1 without equidistribution converges to the wrong value, flat, as the analytic void bias", part5()), ("a wrong time marginal is read off by the slab as 1 - eps", part6()), ] print("\n" + "=" * 100) for name, ok in checks: print(f" {'ok ' if ok else 'FAIL'} {name}") print("=" * 100) print(" THEOREM: |S_h - slice| <= c h^2 + C D_N / h for every cloud, so the slab-wise CAR limit holds") print(" whenever D_N -> 0 faster than h_n. Independence was never the hypothesis; a vanishing") print(" equidistribution rate in the time function is, which is what A3 and P3's tau are.") print(f" ({time.time() - t0:.0f} s)") print(" PASS" if all(ok for _, ok in checks) else " FAIL") main()