1import json, time
 2import numpy as np
 3
 4
 5def fourier_weights(mask, bandwidth, clip=False):
 6    """Periodic Fourier projection of a sampled domain indicator."""
 7    n0, n1 = mask.shape
 8    coeff = np.fft.fft2(mask)
 9    k0 = np.fft.fftfreq(n0) * n0
10    k1 = np.fft.fftfreq(n1) * n1
11    keep = (np.abs(k0[:, None]) <= bandwidth[0]) & (np.abs(k1[None, :]) <= bandwidth[1])
12    weights = np.fft.ifft2(coeff * keep).real
13    if clip:
14        weights = np.maximum(weights, 0.0)
15    return weights
16
17
18def domains(x, y):
19    X, Y = np.meshgrid(x, y, indexing="ij")
20    ellipse = ((X - .48) / .34) ** 2 + ((Y - .52) / .27) ** 2 <= 1
21    hole = ellipse & (((X - .48) / .105) ** 2 + ((Y - .52) / .085) ** 2 > 1)
22    # A cusp-like cardioid, plus a nonconvex boundary.
23    xx, yy = X - .5, Y - .5
24    r = np.sqrt(xx * xx + yy * yy)
25    theta = np.arctan2(yy, xx)
26    cusp = (r <= .27 * (1 + .55 * np.cos(theta)))
27    polygon = (Y > .18 + .18 * X) & (Y < .84 - .26 * X) & (X > .16) & (X < .84)
28    return {"ellipse": ellipse, "ellipse_hole": hole, "cusp": cusp, "polygon": polygon}
29
30
31def fields(x, y):
32    X, Y = np.meshgrid(x, y, indexing="ij")
33    return np.stack([
34        0.4 + 0.7 * X + 0.3 * Y,
35        np.sin(2*np.pi*(1.3*X + .7*Y)) + .2*np.cos(2*np.pi*Y),
36        np.exp(-35*((X-.3)**2 + (Y-.7)**2)),
37        np.sin(2*np.pi*8*X) * np.cos(2*np.pi*5*Y),
38    ])
39
40
41def normalized_pool(f, w):
42    return np.sum(f * w) / np.sum(w)
43
44
45def run():
46    np.random.seed(7)
47    # Reference is much finer than the tested grid; all fields are evaluated analytically.
48    nr = 768
49    xr = (np.arange(nr) + .5) / nr
50    ref_domains = domains(xr, xr)
51    ref_fields = fields(xr, xr)
52    references = {name: np.array([normalized_pool(f, m.astype(float)) for f in ref_fields])
53                  for name, m in ref_domains.items()}
54
55    rows = []
56    for n in (24, 48, 96):
57        x = (np.arange(n) + .5) / n
58        ms = domains(x, x)
59        fs = fields(x, x)
60        for name, m in ms.items():
61            m = m.astype(float)
62            ref = references[name]
63            base = np.array([normalized_pool(f, m) for f in fs])
64            for frac in (.25, .5):
65                w = fourier_weights(m, (int(frac*n), int(frac*n)), clip=False)
66                fw = np.array([normalized_pool(f, w) for f in fs])
67                wc = fourier_weights(m, (int(frac*n), int(frac*n)), clip=True)
68                fc = np.array([normalized_pool(f, wc) for f in fs])
69                rows.append({"n": n, "domain": name, "bandwidth_fraction": frac,
70                             "masked_mean_mae": float(np.mean(np.abs(base-ref))),
71                             "fourier_mae": float(np.mean(np.abs(fw-ref))),
72                             "fourier_clipped_mae": float(np.mean(np.abs(fc-ref))),
73                             "weight_min": float(w.min()), "weight_max": float(w.max()),
74                             "mask_fraction": float(m.mean())})
75
76    # Claimed reusable-weight property: pooling multiple channels is one matrix-weighted sum,
77    # and geometry preprocessing is amortized over channels.
78    n = 64; x = (np.arange(n)+.5)/n; m = domains(x,x)["cusp"].astype(float)
79    t0 = time.perf_counter(); w = fourier_weights(m, (n//2,n//2)); prep = time.perf_counter()-t0
80    channels = fields(x,x)
81    t0 = time.perf_counter()
82    for _ in range(300):
83        _ = np.sum(channels*w, axis=(1,2))/np.sum(w)
84    pool_ms = (time.perf_counter()-t0)*1000/300
85    # Fourier projection's direct indicator approximation error (L2), for the representative cusp.
86    projection_l2 = float(np.sqrt(np.mean((w-m)**2)))
87
88    summary = {
89        "rows": rows,
90        "representative": {"preprocess_ms": prep*1000, "pool_4_channels_ms": pool_ms,
91                            "projection_l2_n64_Fhalf": projection_l2},
92        "mean_ratio_fourier_over_mask": float(np.mean([r["fourier_mae"]/max(r["masked_mean_mae"],1e-12) for r in rows])),
93        "mean_ratio_clipped_over_mask": float(np.mean([r["fourier_clipped_mae"]/max(r["masked_mean_mae"],1e-12) for r in rows]))
94    }
95    with open("results.json", "w") as f: json.dump(summary, f, indent=2)
96    print(json.dumps(summary, indent=2))
97
98if __name__ == "__main__": run()