"""exercise_4b.py Exercise IV-B. The GeoAI Feedback Reflex Assessment. Builds a predictive allocation loop on real tract geography with an event rate set identical everywhere, so any concentration the model produces is an artifact and nothing else. Then adds a reserved random fraction and sweeps it. python exercise_4b.py """ from __future__ import annotations import sys from pathlib import Path import matplotlib import numpy as np import pandas as pd matplotlib.use("Agg") import matplotlib.pyplot as plt # noqa: E402 sys.path.insert(0, str(Path(__file__).resolve().parent)) import common as C # noqa: E402 from mapping import load_geometry # noqa: E402 from matplotlib.collections import PolyCollection # noqa: E402 GEO = C.SNAP / "ms_tracts.geojson" COUNTY = "28049" # Hinds County ITERS, BUDGET = 20, 100 TRUE_RATE = 0.05 # identical in every tract, by construction EXPLORE = [0.0, 0.10, 0.20, 0.30] TOLERANCE = 0.20 def run_loop(n, explore, seed): """One allocation loop. Returns effort history and cumulative effort.""" rng = np.random.default_rng(seed) score = rng.uniform(0.9, 1.1, n) # small initial noise only seen = np.zeros(n) # inspections to date hits = np.zeros(n) # events observed to date cum = np.zeros(n) share = [] for _ in range(ITERS): n_rand = int(round(BUDGET * explore)) n_score = BUDGET - n_rand alloc = np.zeros(n) p = score / score.sum() alloc += rng.multinomial(n_score, p) if n_rand: alloc += rng.multinomial(n_rand, np.full(n, 1.0 / n)) # Events occur at the same rate everywhere; we only see them where we # looked, which is the entire mechanism. obs = rng.binomial(alloc.astype(int), TRUE_RATE) seen += alloc hits += obs cum += alloc with np.errstate(divide="ignore", invalid="ignore"): rate = np.where(seen > 0, hits / np.maximum(seen, 1), 0.0) score = np.maximum(rate, 1e-6) + 1e-3 # refit on observed rate top = np.argsort(cum)[-max(1, n // 10):] share.append(float(alloc[top].sum() / BUDGET)) return np.array(share), cum, seen, hits def main(): C.ensure_dirs() m = C.load_acs() tr = m[(m.level == "tract") & (m.code.str.startswith(COUNTY))].copy() tr = tr.reset_index(drop=True) n = len(tr) C.banner("exercise IV-B: the GeoAI feedback reflex assessment", "Census tract geography, synthetic uniform event rate", "acs2023_mississippi_asi.csv", {"county": COUNTY, "tracts": n, "true rate": TRUE_RATE}) out = {"tracts": n, "county": COUNTY, "true_rate": TRUE_RATE, "iterations": ITERS, "budget": BUDGET} print("\n[1] THE UNCORRECTED LOOP") share0, cum0, seen0, hits0 = run_loop(n, 0.0, C.SEED) print(" the true rate is identical in every tract, so any lasting " "concentration comes from the loop itself") print(f" {'iteration':>10}{'top decile share':>19}") for it in (1, 5, 10, 20): print(f" {it:>10d}{share0[it-1]*100:>18.1f}%") out["uncorrected"] = {f"iter_{it}": round(float(share0[it-1]) * 100, 1) for it in (1, 5, 10, 20)} print("\n[2] FEEDBACK LAUNDERING") # Apparent agreement: correlation between cumulative inspection effort # and the observed event rate, the indicator Table 12.10.5.1 lists. The # program chose where to look, so the agreement cannot separate validity # from self confirmation. with np.errstate(divide="ignore", invalid="ignore"): obs_rate = np.where(seen0 > 0, hits0 / np.maximum(seen0, 1), np.nan) r_self, _ = C.corr(cum0, np.nan_to_num(obs_rate)) print(f" inspections concentrated in {int((cum0 > 0).sum())} of {n} " f"tracts") print(f" tracts never inspected at all: {int((cum0 == 0).sum())}") print(f" apparent agreement between effort and observed rate " f"r = {r_self:+.3f}") print(" The true rate is identical everywhere, so this agreement " "measures self consistency.") out["laundering"] = {"tracts_inspected": int((cum0 > 0).sum()), "tracts_never_inspected": int((cum0 == 0).sum()), "effort_vs_observed_r": round(r_self, 3)} print("\n[3] EXPLORATION SWEEP") print(f" {'explore':>9}{'top decile share at 20':>25}" f"{'within tolerance':>19} (mean of 25 seeds)") # A single run is noisy enough to make the sweep non-monotonic, which # would teach the wrong lesson. Each fraction is the mean of 25 seeds. N_REP = 25 curves, out["sweep"] = {}, [] for e in EXPLORE: reps = np.array([run_loop(n, e, C.SEED + k)[0] for k in range(N_REP)]) s = reps.mean(axis=0) curves[e] = s ok = "yes" if s[-1] <= TOLERANCE else "no" print(f" {e*100:>8.0f}%{s[-1]*100:>24.1f}%{ok:>19}") out["sweep"].append({"explore_pct": e * 100, "share_at_20": round(float(s[-1]) * 100, 1), "replicates": N_REP, "within_tolerance": ok}) passing = [e for e in EXPLORE if curves[e][-1] <= TOLERANCE and e > 0] out["smallest_sufficient_explore_pct"] = (min(passing) * 100 if passing else None) print(f" smallest exploration fraction holding the share at or below " f"{TOLERANCE:.0%}: " f"{out['smallest_sufficient_explore_pct']}%") print("\n[4] INTERVENTION AWARE VALIDATION") # Validate only on the uniformly sampled part, which is the only sample # the allocation policy did not choose. s_r, c_r, seen_r, hits_r = run_loop(n, 0.30, C.SEED) with np.errstate(divide="ignore", invalid="ignore"): rate_r = np.where(seen_r > 0, hits_r / np.maximum(seen_r, 1), np.nan) est = float(np.nansum(hits_r) / max(np.nansum(seen_r), 1)) print(f" pooled rate estimated from the biased loop " f"{float(np.nansum(hits0)/max(np.nansum(seen0),1)):.4f}") print(f" pooled rate with 30 percent random sampling {est:.4f}") print(f" the true rate is {TRUE_RATE:.4f}") out["validation"] = { "biased_estimate": round(float(np.nansum(hits0) / max(np.nansum(seen0), 1)), 4), "explored_estimate": round(est, 4), "true": TRUE_RATE} # ---- figures ---------------------------------------------------------- fig, ax = plt.subplots(figsize=(8.5, 5.5)) for e in EXPLORE: ax.plot(range(1, ITERS + 1), curves[e] * 100, marker="o", ms=3, label=f"{e*100:.0f}% random") ax.axhline(10, ls="--", color="#b2182b", label="ten percent reference") ax.set_xlabel("iteration") ax.set_ylabel("share of budget in the top decile (%)") ax.set_title("A loop finding what it was sent to look for", loc="left", fontweight="bold") ax.text(0, 1.005, "Event rate identical in every tract by construction", transform=ax.transAxes, fontsize=9, va="bottom", color="#444") ax.legend(frameon=False, fontsize=9) ax.text(0, -0.13, f"Synthetic simulation on {n} Hinds County tracts, " f"mean of {N_REP} seeds from {C.SEED}. EPSG:5070.", transform=ax.transAxes, fontsize=7, color="#666") fig.savefig(C.FIG / "exercise_4b_concentration.png", dpi=150, bbox_inches="tight") plt.close(fig) geoms = load_geometry(GEO) polys, vals = [], [] for code, v in zip(tr.code, cum0): for ring in geoms.get(str(code), []): polys.append(ring) vals.append(v) fig, ax = plt.subplots(figsize=(8, 7)) pc = PolyCollection(polys, array=np.array(vals), cmap="magma", edgecolors="white", linewidths=0.3) ax.add_collection(pc) ax.autoscale_view(), ax.set_aspect("equal"), ax.axis("off") ax.set_title("Cumulative inspection effort at iteration 20", loc="left", fontweight="bold") ax.text(0, 1.005, "The underlying event rate is the same in every tract", transform=ax.transAxes, fontsize=9, va="bottom", color="#444") fig.colorbar(pc, ax=ax, shrink=0.7, label="inspections") ax.text(0, -0.03, f"Synthetic simulation on {n} Hinds County tracts, " f"seed {C.SEED}. EPSG:5070.", transform=ax.transAxes, fontsize=7, va="top", color="#666") fig.savefig(C.FIG / "exercise_4b_effort_map.png", dpi=140, bbox_inches="tight") plt.close(fig) tbl = tr[["code"]].copy() tbl["cumulative_inspections"] = cum0 tbl["events_observed"] = hits0 tbl["observed_rate"] = obs_rate tbl["true_rate"] = TRUE_RATE C.write_handoff("exercise_4b", tbl, { "true_rate": TRUE_RATE, "iterations": ITERS, "budget": BUDGET, "seed": C.SEED, "initial_score": "uniform noise between 0.9 and 1.1", "allocation": "multinomial, proportional to score", "refit": "observed events per inspection, floored at 0.001 so a " "tract with no observed event keeps a small share", "tolerance": TOLERANCE, "color_ramp": "magma, unjustified", }, [ "The initial score was random noise. Real systems seed from historical " "records, which are themselves the output of an earlier loop. Does " "that make the problem better or worse?", "No smoothing beyond that floor was applied when refitting. A tract " "inspected once with one event reads as a rate of 1.0. What would a " "prior do here?", "The exploration fraction costs inspections that go where the model " "says not to look. Who pays that cost and who benefits?", "The map is persuasive and describes nothing real. What caption keeps " "a reader from drawing the obvious wrong conclusion?", ]) C.save_checkvalues("exercise_4b", out) print(f"\nwrote {C.FIG / 'exercise_4b_concentration.png'}") print(f"wrote {C.FIG / 'exercise_4b_effort_map.png'}") if __name__ == "__main__": main()