#!/usr/bin/env python3
"""
SPECTRUM PROBE — where does genuine detail stop in terrain-diffusion output?
WHY THIS EXISTS
hyper-games #46 §3. The band split between a learned coarse generator and our
erosion filter is a decision about a wavelength, and it is only sound if we know
the finest wavelength at which the model puts real energy.
WHAT IT MEASURES
The radially-averaged power spectrum of an elevation tile, P(k), fitted to a
power law P(k) ~ k^-beta over a mid band, plus the residual of the top octave
against that fit. A generator that produces genuine detail to its own Nyquist
holds the power law all the way up. One that upsamples falls off a cliff.
THE INSTRUMENT IS A LAYER TOO (#45 section 5)
Every number below is worthless if the spectrum code is wrong, and a broken
ruler reports confidently. So --selftest runs FIRST and must pass:
- white noise -> beta ~ 0 (flat spectrum), the null case
- a single sinusoid -> one sharp peak at the known bin, and nowhere else
If either fails, no terrain number is printed at all.
USAGE
python spectrum_probe.py --selftest
python spectrum_probe.py --url http://127.0.0.1:8817 --size 512 --scale 1
python spectrum_probe.py --url http://127.0.0.1:8817 --compare 1 2
"""
import argparse
import sys
import numpy as np
# ---------------------------------------------------------------------------
# The instrument.
# ---------------------------------------------------------------------------
def radial_power_spectrum(field):
"""
Radially-averaged power spectrum of a 2-D real field.
Returns (k, power) where k is in cycles per field-width (so k=1 is one
cycle across the whole tile) and power is the mean |FFT|^2 in that
annulus. Bin k covers wavelength lambda = N / k pixels.
A Hann window is applied in both axes: the tile is a crop of an unbounded
field, so its edges do not match, and without a window that discontinuity
leaks broadband energy across every bin -- which would look exactly like
the "genuine high-frequency detail" this probe exists to detect.
"""
a = np.asarray(field, dtype=np.float64)
if a.ndim != 2 or a.shape[0] != a.shape[1]:
raise ValueError(f"expected a square 2-D field, got {a.shape}")
n = a.shape[0]
# Remove mean and best-fit plane. A tilt is a k=0/k=1 feature and would
# otherwise dominate every low bin.
yy, xx = np.mgrid[0:n, 0:n].astype(np.float64)
A = np.column_stack([np.ones(n * n), xx.ravel(), yy.ravel()])
coef, *_ = np.linalg.lstsq(A, a.ravel(), rcond=None)
a = a - (A @ coef).reshape(n, n)
w = np.hanning(n)
a = a * w[:, None] * w[None, :]
f = np.fft.fftshift(np.fft.fft2(a))
power = (np.abs(f) ** 2) / (n * n)
cy = cx = n // 2
yy, xx = np.mgrid[0:n, 0:n]
r = np.sqrt((yy - cy) ** 2 + (xx - cx) ** 2)
rint = r.astype(np.int64)
nbins = n // 2
tot = np.bincount(rint.ravel(), weights=power.ravel(), minlength=nbins + 1)
cnt = np.bincount(rint.ravel(), minlength=nbins + 1)
with np.errstate(invalid="ignore", divide="ignore"):
radial = tot / np.maximum(cnt, 1)
k = np.arange(1, nbins) # drop k=0 (DC, removed anyway)
return k, radial[1:nbins]
def fit_power_law(k, p, kmin, kmax):
"""Least-squares fit of log10(P) = c - beta*log10(k) over [kmin, kmax]."""
m = (k >= kmin) & (k <= kmax) & (p > 0)
if m.sum() < 4:
raise ValueError(f"only {m.sum()} usable bins in [{kmin},{kmax}]")
lk, lp = np.log10(k[m]), np.log10(p[m])
slope, intercept = np.polyfit(lk, lp, 1)
pred = intercept + slope * lk
ss_res = float(np.sum((lp - pred) ** 2))
ss_tot = float(np.sum((lp - lp.mean()) ** 2))
r2 = 1.0 - ss_res / ss_tot if ss_tot > 0 else float("nan")
return -slope, intercept, r2
def top_octave_residual(k, p, kmin, kmax):
"""
How far the top octave sits below the power law fitted BELOW it, in dB.
Fit over [kmin, kmax], extrapolate into (kmax, nyquist], and report the
median dB shortfall. ~0 dB means the power law continues -- genuine
detail. A large negative number means the energy is not there.
"""
beta, intercept, r2 = fit_power_law(k, p, kmin, kmax)
m = (k > kmax) & (p > 0)
if m.sum() < 3:
return beta, r2, float("nan"), 0
pred_log = intercept - beta * np.log10(k[m])
resid_db = 10.0 * (np.log10(p[m]) - pred_log)
return beta, r2, float(np.median(resid_db)), int(m.sum())
# ---------------------------------------------------------------------------
# Self-test. Nothing downstream runs unless this passes.
# ---------------------------------------------------------------------------
def selftest(n=512, verbose=True):
ok = True
# 1. White noise -> flat spectrum.
rng = np.random.default_rng(20260831)
k, p = radial_power_spectrum(rng.standard_normal((n, n)))
beta, _, r2 = fit_power_law(k, p, 8, n // 2 - 4)
passed = abs(beta) < 0.10
ok &= passed
if verbose:
print(f" [{'PASS' if passed else 'FAIL'}] white noise "
f"beta = {beta:+.4f} (expect 0.000 +/- 0.100)")
# 2. Single sinusoid at a known wavelength -> one peak, right bin.
lam = 32.0 # pixels per cycle
expect = int(round(n / lam)) # -> bin 16
xx = np.arange(n)
field = np.sin(2 * np.pi * xx / lam)[None, :] * np.ones((n, 1))
k, p = radial_power_spectrum(field)
peak = int(k[np.argmax(p)])
passed = peak == expect
ok &= passed
if verbose:
print(f" [{'PASS' if passed else 'FAIL'}] sinusoid lam={lam:.0f}px "
f"peak at k = {peak} (expect {expect})")
# 3a. Leakage check, on the UN-AVERAGED 2-D spectrum. A 1-D sinusoid puts
# all its energy in two bins; everything else is windowing leakage and
# must sit far below. This is the honest test of the FFT + window.
ff = np.fft.fftshift(np.fft.fft2(
(field - field.mean()) * np.hanning(n)[:, None] * np.hanning(n)[None, :]))
pw = np.abs(ff) ** 2
cy = cx = n // 2
yy2, xx2 = np.mgrid[0:n, 0:n]
near = (np.abs(xx2 - cx) <= expect + 2) & (np.abs(xx2 - cx) >= expect - 2) & \
(np.abs(yy2 - cy) <= 2)
far = ~near
ratio_db = 10 * np.log10(pw[near].max() / pw[far].max())
passed = ratio_db > 30
ok &= passed
if verbose:
print(f" [{'PASS' if passed else 'FAIL'}] leakage (2-D) "
f"peak is {ratio_db:.1f} dB above off-axis (expect > 30.0)")
# 3b. Purity through the RADIAL path, with an isotropic signal. Radial
# averaging divides an annulus's energy by its pixel count, so it
# dilutes a directional peak by ~10*log10(2*pi*r/2) -- about 17 dB at
# r=16 -- by construction, not by fault. The instrument is built for
# isotropic fields, so the check uses one: energy on a ring at k=16
# and nowhere else.
kx = np.fft.fftfreq(n) * n
KX, KY = np.meshgrid(kx, kx)
KR = np.sqrt(KX ** 2 + KY ** 2)
ring = (np.abs(KR - expect) < 0.5).astype(np.float64)
ph = rng.uniform(0, 2 * np.pi, (n, n))
field_iso = np.real(np.fft.ifft2(ring * np.exp(1j * ph)))
k, p = radial_power_spectrum(field_iso)
peak_iso = int(k[np.argmax(p)])
mask = np.abs(k - expect) > 2
ratio_db = 10 * np.log10(p.max() / p[mask].max())
passed = (peak_iso == expect) and ratio_db > 30
ok &= passed
if verbose:
print(f" [{'PASS' if passed else 'FAIL'}] ring k={expect} isotropic "
f"peak at k = {peak_iso}, {ratio_db:.1f} dB above next "
f"(expect {expect}, > 30.0)")
# 4. Known power law in -> same power law out.
for target in (1.0, 2.0):
f = np.fft.fftfreq(n) * n
kx, ky = np.meshgrid(f, f)
kr = np.sqrt(kx ** 2 + ky ** 2)
kr[0, 0] = 1.0
amp = kr ** (-target / 2.0)
ph = rng.uniform(0, 2 * np.pi, (n, n))
spec = amp * np.exp(1j * ph)
field = np.real(np.fft.ifft2(spec))
k, p = radial_power_spectrum(field)
beta, _, r2 = fit_power_law(k, p, 8, n // 4)
passed = abs(beta - target) < 0.20
ok &= passed
if verbose:
print(f" [{'PASS' if passed else 'FAIL'}] synthetic k^-{target:.0f} "
f"beta = {beta:.3f} (expect {target:.3f} +/- 0.200), R2 = {r2:.3f}")
return ok
# ---------------------------------------------------------------------------
# The terrain fetch.
# ---------------------------------------------------------------------------
def fetch_tile(url, i1, j1, size, scale, timeout=1800):
import urllib.request
q = (f"{url}/terrain?i1={i1}&j1={j1}&i2={i1 + size}&j2={j1 + size}"
f"&scale={scale}")
with urllib.request.urlopen(q, timeout=timeout) as r:
h = int(r.headers["X-Height"])
w = int(r.headers["X-Width"])
body = r.read()
nelev = h * w * 2
elev = np.frombuffer(body[:nelev], dtype="<i2").reshape(h, w).astype(np.float64)
climate = None
if len(body) > nelev:
climate = np.frombuffer(body[nelev:], dtype="<f4").reshape(h, w, 4)
return elev, climate
def describe(name, elev, px_metres, kmin, kmax):
k, p = radial_power_spectrum(elev)
n = elev.shape[0]
beta, r2, resid_db, nbins = top_octave_residual(k, p, kmin, kmax)
print(f"\n {name}")
print(f" tile {n} x {n} px @ {px_metres:g} m/px "
f"= {n * px_metres / 1000:.1f} km across")
print(f" elevation min {elev.min():.0f} m max {elev.max():.0f} m "
f"mean {elev.mean():.0f} m sd {elev.std():.0f} m")
print(f" fit band k = {kmin}..{kmax} "
f"(lambda {n * px_metres / kmax:.0f} m .. {n * px_metres / kmin / 1000:.1f} km)")
print(f" power law P(k) ~ k^-{beta:.3f} R2 = {r2:.4f}")
print(f" above k={kmax} {resid_db:+.1f} dB vs the fit, "
f"median over {nbins} bins")
print(f" Nyquist lambda = {2 * px_metres:g} m")
return k, p, beta, resid_db
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--url", default="http://127.0.0.1:8817")
ap.add_argument("--size", type=int, default=512)
ap.add_argument("--i1", type=int, default=0)
ap.add_argument("--j1", type=int, default=0)
ap.add_argument("--native-metres", type=float, default=30.0,
help="metres per pixel at scale=1 for the loaded model")
ap.add_argument("--compare", type=int, nargs="*", default=[1],
help="scale factors to fetch and compare")
ap.add_argument("--selftest", action="store_true",
help="run the instrument checks and exit")
ap.add_argument("--csv", default=None, help="write spectra to this CSV")
args = ap.parse_args()
print("INSTRUMENT SELF-TEST")
if not selftest():
print("\n FAILED. The ruler is broken; no terrain number is meaningful.")
return 2
print(" all pass\n")
if args.selftest:
return 0
n = args.size
kmin, kmax = 6, n // 8 # fit below the top two octaves, extrapolate up
print(f"TERRAIN {args.url} region ({args.i1},{args.j1}) size {n}")
results = {}
for s in args.compare:
# Same WORLD region at every scale: coordinates are in target pixels.
# Same WORLD footprint at every scale: n native px wide. At scale s
# that is n*s target px starting at i1*s. Do NOT crop to n -- cropping
# would compare a smaller region, not a finer sampling of the same one.
elev, climate = fetch_tile(args.url, args.i1 * s, args.j1 * s, n * s, s)
elev_c = elev
px = args.native_metres / s
# kmin/kmax stay FIXED across scales. n*px is the same physical width
# for every scale, so a given k is the same physical wavelength -- the
# fit band is therefore identical and the residual above it asks
# exactly the question: does a finer scale add energy up there?
k, p, beta, resid = describe(
f"scale={s} ({px:g} m/px)", elev_c, px, kmin, kmax)
results[s] = (k, p, beta, resid)
if climate is not None:
c = climate
print(f" climate temp {c[...,0].mean():.1f} C "
f"t_season {c[...,1].mean():.0f} "
f"precip {c[...,2].mean():.0f} mm/yr "
f"p_cv {c[...,3].mean():.0f}%")
if args.csv:
import csv
with open(args.csv, "w", newline="") as fh:
wr = csv.writer(fh)
wr.writerow(["scale", "k", "power"])
for s, (k, p, _, _) in results.items():
for kk, pp in zip(k, p):
wr.writerow([s, int(kk), f"{pp:.10e}"])
print(f"\n spectra written to {args.csv}")
return 0
if __name__ == "__main__":
sys.exit(main())