#!/usr/bin/env python3
"""
DOES THE FILTER ENHANCE THE FIELD'S EXISTING LINES, OR ADD ITS OWN?
erosion_filter.inc:549 places gullies with phacelle_noise(p*freq, gully_slope,..):
the PATTERN is noise, merely ORIENTED downhill. So a priori the filter adds a
new gully system rather than deepening the drainage the diffusion field already
has. This measures whether that is true, and whether assumed_slope_mix changes it.
WHY THIS RE-TESTS A KNOB ALREADY "REFUTED". sweep_assumed_slope.py found mix
inert -- but it ran at the shipped operating point, where the erosion moves the
surface 0.227 m. A knob that modulates an effect which is not happening reads as
inert no matter what it does. That measurement proved nothing; this one runs at
10 m, where there is something to modulate.
THE TWO CORRELATIONS, both against the SOURCE field (never against a render, so
they cannot share a bug with what they judge):
corr(|ero|, slope) does the filter cut harder on steep ground
corr(|ero|, concavity) does it cut harder where the field is ALREADY a
valley. This is the "enhance existing lines" question.
Concavity = -laplacian, so positive = valley floor.
A filter that enhances has concavity correlation well above zero. A filter that
adds an independent pattern has it near zero however deep it cuts.
THE CONTROL is the correlation of the erosion-OFF field with itself-derived
concavity, which is meaningless by construction, plus a SHUFFLED |ero| whose
correlation must collapse to ~0. If the shuffle does not collapse, the statistic
is picking up something structural in the window rather than in the filter.
python follows_structure.py [--depth 10]
"""
import argparse, os
import numpy as np
from sweep_assumed_slope import render, MAP_PARAMS
from stitch_block import read_hblk
HERE = os.path.dirname(os.path.abspath(__file__))
HBLK = os.path.abspath(os.path.join(HERE, "..", "hyper-tps", "assets",
"terrain", "diffusion_gully_30km.hblk"))
MAP = dict(MAP_PARAMS)
MEAS_RMS, MEAS_STRENGTH = 0.227, MAP["en_strength"]
MIXES = [1.0, 0.75, 0.5, 0.25, 0.0]
def args_with(**over):
d = dict(MAP_PARAMS); d.update(over)
out = []
for k, v in d.items():
out += ["--set", k, str(v)]
return out
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--bin", default="/tmp/erosion_render")
ap.add_argument("--ox", type=float, default=50000.0)
ap.add_argument("--oz", type=float, default=34000.0)
ap.add_argument("--n", type=int, default=512)
ap.add_argument("--mpp", type=float, default=8.0)
ap.add_argument("--depth", type=float, default=10.0)
a = ap.parse_args()
strength = MEAS_STRENGTH * (a.depth / MEAS_RMS)
print(f"window {a.n*a.mpp/1000:.2f} km @ {a.mpp} m/px, depth {a.depth:g} m rms "
f"(strength {strength:.5g})\n")
# SOURCE field, sampled independently of any render.
arr, meta = read_hblk(HBLK)
ii = ((np.arange(a.n)*a.mpp + a.ox) - meta["origin_x"]) / meta["mpp"]
jj = ((np.arange(a.n)*a.mpp + a.oz) - meta["origin_z"]) / meta["mpp"]
gi = np.clip(np.round(ii).astype(int), 0, meta["side"]-1)
gj = np.clip(np.round(jj).astype(int), 0, meta["side"]-1)
hsrc = arr[np.ix_(gj, gi)].astype(np.float64)
sz, sx = np.gradient(hsrc, a.mpp)
slope = np.sqrt(sx*sx + sz*sz)
lap = (np.gradient(np.gradient(hsrc, a.mpp, axis=0), a.mpp, axis=0) +
np.gradient(np.gradient(hsrc, a.mpp, axis=1), a.mpp, axis=1))
concav = -lap # positive = already a valley
off = render(a.bin, a.ox, a.oz, a.n, a.mpp, "/tmp/fs_off.f32",
args_with(en_strength=0))
rng = np.random.default_rng(0)
print(f"{'mix':>6} {'rms m':>8} {'corr(|ero|,slope)':>19} {'corr(|ero|,concav)':>20} "
f"{'shuffled':>10}")
out = []
for m in MIXES:
h = render(a.bin, a.ox, a.oz, a.n, a.mpp, f"/tmp/fs_{m}.f32",
args_with(en_strength=strength, assumed_slope_mix=m))
ae = np.abs(h - off)
cs = np.corrcoef(ae.ravel(), slope.ravel())[0, 1]
cc = np.corrcoef(ae.ravel(), concav.ravel())[0, 1]
sh = np.corrcoef(rng.permutation(ae.ravel()), concav.ravel())[0, 1]
print(f"{m:>6.2f} {np.sqrt(((h-off)**2).mean()):>8.2f} {cs:>19.3f} "
f"{cc:>20.3f} {sh:>10.3f}")
out.append((m, cs, cc, sh))
cc0, cc1 = out[0][2], out[-1][2]
print(f"\nconcavity correlation: {cc0:+.3f} (mix 1.0) -> {cc1:+.3f} (mix 0.0)")
print(f"shuffle control : {max(abs(o[3]) for o in out):.3f} "
"(must be ~0 or the statistic is not measuring the filter)")
best = max(out, key=lambda o: abs(o[2]))
if abs(best[2]) < 0.15:
print("\nREADS AS: the filter ADDS an independent pattern. It does not\n"
" find the field's existing valleys at ANY mix setting.")
else:
print(f"\nREADS AS: at mix {best[0]:g} the filter DOES follow existing valleys "
f"(corr {best[2]:+.3f}).")
if __name__ == "__main__":
main()