"""Setup: paths, palette, helpers. Run this first."""
from __future__ import annotations

import hashlib
import json
from pathlib import Path

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import rasterio
import yaml
from matplotlib.colors import LinearSegmentedColormap

REPO = Path.cwd().parent if Path.cwd().name == "notebooks" else Path.cwd()
OUT = REPO / "outputs" / "psf_disaggregation"
VAL = OUT / "validation"
RAS = OUT / "rasters" / "full_city"
CONFIG_PATH = REPO / "configs" / "psf_disaggregation.yaml"
SAVEFIG_PATH = REPO / "docs" / "figures"
SAVEFIG_PATH.mkdir(parents=True, exist_ok=True)

CITIES = ["usa_new_york", "india_delhi"]
CITY_LABEL = {"usa_new_york": "New York", "india_delhi": "Delhi"}

BROAD = RAS / "v2_broad_qa_reporting"
STRICT = RAS / "v1_resumable_tiled"
PRIMARY = "built_form_primary__no_water_prior__circular_mean_reference__broad_qa"
S2_STRICT = "s2_only_ablation__no_water_prior__circular_mean_reference"

# --- palette (validated categorical order; see dataviz reference palette) -----
BLUE, ORANGE, AQUA, YELLOW, MAGENTA, GREEN, VIOLET, RED = (
    "#2a78d6", "#eb6834", "#1baf7a", "#eda100",
    "#e87ba4", "#008300", "#4a3aa7", "#e34948",
)
INK, INK_2, MUTED = "#0b0b0b", "#52514e", "#8a8984"
GRID = "#e6e5e1"

# role assignments held fixed across every figure in this notebook
COLOR = {
    "structural": BLUE,
    "neighbors": ORANGE,
    "built_form_primary": BLUE,
    "s2_only_ablation": AQUA,
    "direct": MUTED,
    "uniform": INK_2,
    "improve": BLUE,   # diverging cool pole: error goes down
    "harm": RED,       # diverging warm pole: error goes up
}

# single-hue sequential ramp for radiance magnitude (never a rainbow)
BLUE_RAMP = LinearSegmentedColormap.from_list(
    "nocturne_blue",
    ["#fcfcfb", "#cde2fb", "#9ec5f4", "#5598e7", "#2a78d6", "#184f95", "#0d366b"],
)

mpl.rcParams.update({
    "figure.dpi": 120,
    "savefig.dpi": 120,
    "font.size": 9,
    "axes.edgecolor": GRID,
    "axes.labelcolor": INK_2,
    "axes.titlesize": 10,
    "axes.titlelocation": "left",
    "axes.titlecolor": INK,
    "axes.spines.top": False,
    "axes.spines.right": False,
    "xtick.color": INK_2,
    "ytick.color": INK_2,
    "grid.color": GRID,
    "grid.linewidth": 0.6,
    "legend.frameon": False,
})


def load_json(path: Path):
    with open(path) as fh:
        return json.load(fh)


def sha256_file(path: Path) -> str:
    h = hashlib.sha256()
    with open(path, "rb") as fh:
        for chunk in iter(lambda: fh.read(1 << 20), b""):
            h.update(chunk)
    return h.hexdigest()


def band(path: Path, index: int, window=None) -> np.ndarray:
    with rasterio.open(path) as src:
        return src.read(index, window=window, masked=False)


def band_names(path: Path) -> list[str]:
    with rasterio.open(path) as src:
        return list(src.descriptions)


CONFIG = yaml.safe_load(open(CONFIG_PATH))
print(f"repo   : {REPO}")
print(f"cities : {', '.join(CITY_LABEL.values())}")

repo   : /home/daynan/research/ntl-psf-disaggregation
cities : New York, Delhi

"""Artifact inventory. Every path this notebook reads, checked up front."""
ARTIFACTS = {
    "config": CONFIG_PATH,
    "selection_v2_broad": VAL / "gate2/primary_selection/v2_broad_qa_reporting/selection.json",
    "selection_v1_strict": VAL / "gate2/primary_selection/v1_post_gate2/selection.json",
    "evidence_v2": VAL / "gate2/closeout/v2_broad_qa_reporting/evidence_classification.json",
    "headline_v2": VAL / "gate2/closeout/v2_broad_qa_reporting/headline_metrics.csv",
    "headline_v1": VAL / "gate2/closeout/v1_bounded_method_closeout/headline_metrics.csv",
    "heldout_broad": VAL / "gate2/heldout/v3_physics_buffered_native_cell_broad_qa/summary.csv",
    "heldout_broad_gain": VAL / "gate2/heldout/v3_physics_buffered_native_cell_broad_qa/gain_summary.csv",
    "heldout_broad_retain": VAL / "gate2/heldout/v3_physics_buffered_native_cell_broad_qa/sample_retention.csv",
    "heldout_strict": VAL / "gate2/heldout/v2_physics_buffered_native_cell/summary.csv",
    "heldout_strict_gain": VAL / "gate2/heldout/v2_physics_buffered_native_cell/gain_summary.csv",
    "water_broad": VAL / "gate2/broad_reporting/v1_coverage_and_water_inland/water_inland_ratio.csv",
    "water_strict": VAL / "gate2/day4_allocated_direct_ratio/v1_existing_matched_inland_control/region_metrics.csv",
    "water_compare_strict": VAL / "gate2/day4_allocated_direct_ratio/v1_existing_matched_inland_control/water_vs_inland_comparison.csv",
    "shoreline_summary": VAL / "gate2/shoreline/v1_overture_mapped_water/summary.csv",
    "shoreline_strata": VAL / "gate2/shoreline/v1_overture_mapped_water/waterfront_infrastructure_strata.csv",
    "kernel_sensitivity": VAL / "gate2/sensitivity/v1_full_city_sensitivity/kernel_sensitivity.csv",
    "obs_conditions": VAL / "gate2/observation_conditions/v1_daily_coarse_operator/condition_summary.csv",
    "gate1_closeout": VAL / "gate1/full_city_machine_closeout.md",
    "gate0_builtform": VAL / "gate0/v3_qa_grid_halo_datatake_support/built_form_gate0_summary.csv",
    "gate0_s2": VAL / "gate0/v3_qa_grid_halo_datatake_support/s2_only_gate0_summary.csv",
    "gate0_regen": VAL / "gate0/v3_qa_grid_halo_datatake_support/gate0_figure_regeneration.json",
    "manifest_broad": BROAD / "artifact_manifest.json",
}
for city in CITIES:
    ARTIFACTS[f"primary_{city}"] = BROAD / city / PRIMARY / "products.tif"
    ARTIFACTS[f"direct_{city}"] = BROAD / city / "direct_upsample__broad_qa" / "products.tif"
    ARTIFACTS[f"uniform_{city}"] = BROAD / city / "uniform_normalized_convolution__broad_qa" / "products.tif"
    ARTIFACTS[f"trust_{city}"] = BROAD / city / PRIMARY / "trust_indicators.tif"
    ARTIFACTS[f"s2_{city}"] = STRICT / city / S2_STRICT / "products.tif"
    ARTIFACTS[f"structure_{city}"] = OUT / "rasters/day2_inputs" / city / "overture_structure_bundle.tif"

inventory = pd.DataFrame(
    [{"key": k, "exists": p.exists(),
      "size_mb": round(p.stat().st_size / 1e6, 1) if p.exists() else None}
     for k, p in ARTIFACTS.items()]
)
missing = inventory.loc[~inventory.exists, "key"].tolist()
print(f"{int(inventory.exists.sum())}/{len(inventory)} artifacts present")
if missing:
    print("MISSING:", missing)
inventory

35/35 artifacts present

"""What reaggregation error actually measures on the reporting-primary product."""
rows = []
for city in CITIES:
    m = load_json(BROAD / city / PRIMARY / "metrics.json")["metrics"]
    direct = band(ARTIFACTS[f"direct_{city}"], 1)
    mean_abs = float(np.nanmean(np.abs(direct)))
    rows.append({
        "city": CITY_LABEL[city],
        "consistency_MAE": m["operator_consistency_mae"],
        "consistency_RMSE": m["operator_consistency_rmse"],
        "consistency_bias": m["operator_consistency_bias"],
        "mean_abs_source_radiance": mean_abs,
        "MAE_pct_of_mean": 100 * m["operator_consistency_mae"] / mean_abs,
    })
    del direct

consistency = pd.DataFrame(rows)
display(consistency.round(3))

print(
    "\nReaggregating the 10 m field returns the coarse observation to within "
    f"{consistency.MAE_pct_of_mean.min():.1f}-{consistency.MAE_pct_of_mean.max():.1f}% "
    "of mean absolute source radiance."
    "\nThat is a real residual, not a rounding artifact. The registered term for "
    "this operator is\n'locally normalized', and the config sets "
    "maximum_error_for_conserving_language = null, i.e. no error threshold has\n"
    "been declared that would license conserving language."
)

/tmp/ipykernel_4640/2449412562.py:92: DeprecationWarning: Setting the shape on a NumPy array has been deprecated in NumPy 2.5.
As an alternative, you can create a new view using np.reshape (with copy=False if needed).
  return src.read(index, window=window, masked=False)

Reaggregating the 10 m field returns the coarse observation to within 10.4-12.2% of mean absolute source radiance.
That is a real residual, not a rounding artifact. The registered term for this operator is
'locally normalized', and the config sets maximum_error_for_conserving_language = null, i.e. no error threshold has
been declared that would license conserving language.

"""Band inventory of a reporting-primary product."""
p = ARTIFACTS["primary_usa_new_york"]
with rasterio.open(p) as src:
    meta = dict(crs=str(src.crs), shape=src.shape, res=src.res,
                dtype=src.dtypes[0], nodata=src.nodata)
print(json.dumps(meta, indent=2, default=str))

roles = {
    1: "the product", 2: "reaggregation residual (section 2)",
    3: "support", 4: "support", 5: "support", 6: "support",
    7: "validity", 8: "validity", 9: "validity",
    10: "invalid-input mask", 11: "invalid-input mask",
    12: "numerical guard", 13: "numerical guard",
    14: "numerical guard", 15: "numerical guard", 16: "numerical guard",
}
display(pd.DataFrame([{"band": i, "name": n, "role": roles[i]}
                      for i, n in enumerate(band_names(p), 1)]))

t = ARTIFACTS["trust_usa_new_york"]
print(f"\ncompanion file: {t.name}")
display(pd.DataFrame([{"band": i, "name": n} for i, n in enumerate(band_names(t), 1)]))
print("FINE-GAIN-001: the trust indicator is a SEPARATE file, computed from proxy "
      "and kernel only.\nIt never touches radiance. Bands 2-6 of products.tif are "
      "operator diagnostics, which is a different thing.")

{
  "crs": "EPSG:32618",
  "shape": [
    5000,
    5000
  ],
  "res": [
    10.0,
    10.0
  ],
  "dtype": "float32",
  "nodata": NaN
}

companion file: trust_indicators.tif

FINE-GAIN-001: the trust indicator is a SEPARATE file, computed from proxy and kernel only.
It never touches radiance. Bands 2-6 of products.tif are operator diagnostics, which is a different thing.

"""Cutout helper: geographic center -> window on the locked UTM grid."""
from pyproj import Transformer

def cutout(path: Path, index: int, lat: float, lon: float, half_km: float = 2.0):
    with rasterio.open(path) as src:
        tx = Transformer.from_crs("EPSG:4326", src.crs, always_xy=True)
        x, y = tx.transform(lon, lat)
        half = half_km * 1000.0
        row0, col0 = src.index(x - half, y + half)
        row1, col1 = src.index(x + half, y - half)
        row0, col0 = max(row0, 0), max(col0, 0)
        row1 = min(row1, src.height); col1 = min(col1, src.width)
        win = rasterio.windows.Window(col0, row0, col1 - col0, row1 - row0)
        return src.read(index, window=win, masked=False)


def structure_cutout(city: str, index: int, lat: float, lon: float, half_km: float = 2.0):
    """Structure bundle is 5200x5200 with a 100 px halo; same CRS and resolution."""
    return cutout(ARTIFACTS[f"structure_{city}"], index, lat, lon, half_km)


SITES = {
    "usa_new_york": {
        "dense":   ("Midtown Manhattan",        40.7580, -73.9855),
        "water":   ("Brooklyn / East River",    40.7033, -73.9903),
        "failure": ("Central Park",             40.7812, -73.9665),
    },
    "india_delhi": {
        "dense":   ("Connaught Place",          28.6315,  77.2167),
        "water":   ("Yamuna corridor",          28.6280,  77.2530),
        "failure": ("Central Ridge forest",     28.5930,  77.1830),
    },
}
print("cutout sites resolved against each city's locked UTM grid")

cutout sites resolved against each city's locked UTM grid

"""Triptych: coarse observation -> structural proxy -> allocation."""
def triptych(city: str, site_key: str, half_km: float = 2.0):
    label, lat, lon = SITES[city][site_key]
    d = cutout(ARTIFACTS[f"direct_{city}"], 1, lat, lon, half_km)
    rho = structure_cutout(city, 3, lat, lon, half_km)   # built_form_base_proxy
    a = cutout(ARTIFACTS[f"primary_{city}"], 1, lat, lon, half_km)

    vmax = float(np.nanpercentile(np.concatenate([d.ravel(), a.ravel()]), 99.0))
    fig, axes = plt.subplots(1, 3, figsize=(11, 3.9))
    for ax, arr, title, cmap, vm in [
        (axes[0], d, "VIIRS obs. radiance (corrected, upsample)", "viridis", vmax),
        (axes[1], rho, "structural proxy $\\rho$", "gray",
         float(np.nanpercentile(rho, 99.0))),
        (axes[2], a, "allocation after operator, 10 m", "viridis", vmax),
    ]:
        im = ax.imshow(arr, cmap=cmap, vmin=0, vmax=vm, interpolation="nearest")
        ax.set_title(title)
        ax.set_xticks([]); ax.set_yticks([])
        fig.colorbar(im, ax=ax, fraction=0.046, pad=0.03)
    fig.suptitle(f"{CITY_LABEL[city]} — {label}   ({half_km*2:.0f} x {half_km*2:.0f} km, "
                 "nighttime radiance allocated by structure)", x=0.5, y=1.02, fontsize=10, color=INK)
    fig.tight_layout()
    plt.savefig(SAVEFIG_PATH / f'tryptych_{city}_{site_key}.png', dpi=600, bbox_inches='tight')
    plt.close(fig)
    

triptych("usa_new_york", "dense")
triptych("usa_new_york", "water")

/tmp/ipykernel_4640/2247567552.py:14: DeprecationWarning: Setting the shape on a NumPy array has been deprecated in NumPy 2.5.
As an alternative, you can create a new view using np.reshape (with copy=False if needed).
  return src.read(index, window=win, masked=False)

triptych("india_delhi", "dense")
triptych("india_delhi", "water")

/tmp/ipykernel_4640/2247567552.py:14: DeprecationWarning: Setting the shape on a NumPy array has been deprecated in NumPy 2.5.
As an alternative, you can create a new view using np.reshape (with copy=False if needed).
  return src.read(index, window=win, masked=False)

"""Failure cutout: low proxy support, with the diagnostic bands that flag it."""
def failure_panel(city: str, half_km: float = 2.5):
    label, lat, lon = SITES[city]["failure"]
    panels = [
        (cutout(ARTIFACTS[f"direct_{city}"], 1, lat, lon, half_km),
         "direct upsample (observation)", "viridis", None),
        (cutout(ARTIFACTS[f"primary_{city}"], 1, lat, lon, half_km),
         "allocation", "viridis", None),
        (structure_cutout(city, 3, lat, lon, half_km),
         "proxy $\\rho$", "gray", None),
        (cutout(ARTIFACTS[f"primary_{city}"], 2, lat, lon, half_km),
         "band 2: reaggregation residual", "RdBu_r", "sym"),
        (cutout(ARTIFACTS[f"primary_{city}"], 1, lat, lon, half_km)
         - cutout(ARTIFACTS[f"direct_{city}"], 1, lat, lon, half_km),
         "allocation - direct: where radiance moved", "RdBu_r", "sym"),
        (cutout(ARTIFACTS[f"trust_{city}"], 1, lat, lon, half_km),
         "trust band: allocation gain", "RdBu_r", (0, 2)),
    ]
    fig, axes = plt.subplots(2, 3, figsize=(11, 6.6))
    for ax, (arr, title, cmap, scale) in zip(axes.ravel(), panels):
        if scale == "sym":
            lim = float(np.nanpercentile(np.abs(arr), 99))
            im = ax.imshow(arr, cmap=cmap, vmin=-lim, vmax=lim, interpolation="nearest")
        elif isinstance(scale, tuple):
            im = ax.imshow(arr, cmap=cmap, vmin=scale[0], vmax=scale[1], interpolation="nearest")
        else:
            im = ax.imshow(arr, cmap=cmap, vmin=0,
                           vmax=float(np.nanpercentile(arr, 99)), interpolation="nearest")
        ax.set_title(title); ax.set_xticks([]); ax.set_yticks([])
        fig.colorbar(im, ax=ax, fraction=0.046, pad=0.03)
    fig.suptitle(f"{CITY_LABEL[city]} — {label}: the operator's hardest honest case",
                 x=0.5, y=1.00, fontsize=10, color=INK)
    fig.tight_layout()
    plt.savefig(SAVEFIG_PATH / f'failure_panel_{city}.png', dpi=600, bbox_inches='tight')
    plt.close(fig)

failure_panel("usa_new_york")

/tmp/ipykernel_4640/2247567552.py:14: DeprecationWarning: Setting the shape on a NumPy array has been deprecated in NumPy 2.5.
As an alternative, you can create a new view using np.reshape (with copy=False if needed).
  return src.read(index, window=win, masked=False)

failure_panel("india_delhi")

/tmp/ipykernel_4640/2247567552.py:14: DeprecationWarning: Setting the shape on a NumPy array has been deprecated in NumPy 2.5.
As an alternative, you can create a new view using np.reshape (with copy=False if needed).
  return src.read(index, window=win, masked=False)

"""Gate 0 v3 summary tables."""
g0_bf = pd.read_csv(ARTIFACTS["gate0_builtform"])
g0_s2 = pd.read_csv(ARTIFACTS["gate0_s2"])
print("built-form primary:")
display(g0_bf)
print("S2-only ablation:")
display(g0_s2)

regen = load_json(ARTIFACTS["gate0_regen"])
print("figure regeneration record (deviation #2 from section 4):")
print(json.dumps({k: v for k, v in regen.items() if not isinstance(v, (list, dict))},
                 indent=2))

built-form primary:

S2-only ablation:

figure regeneration record (deviation #2 from section 4):
{
  "artifact_version": "v3_qa_grid_halo_datatake_support",
  "created_at_utc": "2026-07-30T15:53:47.609707+00:00",
  "defect_closed": "ARTIFACT-PROVENANCE-001",
  "schema_version": 1,
  "source_policy": "saved v3 samples and v3 summary tables only"
}

"""The regenerated v3 Gate 0 panels."""
from matplotlib.image import imread

g0_dir = VAL / "gate0/v3_qa_grid_halo_datatake_support"
panels = [
    ("usa_new_york_built_form_gate0.png", "New York: built form"),
    ("india_delhi_built_form_gate0.png", "Delhi: built form"),
    ("usa_new_york_s2_only_gate0.png", "New York: S2 only"),
    ("india_delhi_s2_only_gate0.png", "Delhi: S2 only"),
]
fig, axes = plt.subplots(2, 2, figsize=(12, 9))
for ax, (fn, title) in zip(axes.ravel(), panels):
    path = g0_dir / fn
    if path.exists():
        ax.imshow(imread(path))
        ax.set_title(title)
    ax.axis("off")
fig.tight_layout()
plt.savefig(SAVEFIG_PATH / "gate0-panels.png", dpi=300, bbox_inches="tight")
plt.close(fig)

"""Held-out design: target and neighbor retention across buffers (broad QA)."""
ret = pd.read_csv(ARTIFACTS["heldout_broad_retain"])
print("columns:", list(ret.columns))
display(ret.head(20))

columns: ['city_id', 'proxy', 'buffer_m', 'buffer_role', 'fold', 'target_count', 'retained_target_count', 'retained_target_fraction', 'median_nearest_neighbor_distance_m', 'median_distance_to_block_edge_m']

"""Confirm the design property: identical target counts across buffers."""
hb = pd.read_csv(ARTIFACTS["heldout_broad"])
overall = hb[(hb.fold.astype(str) == "all") & (hb.radiance_decile.astype(str) == "all")]
tgt = (overall[overall.method == "neighbors_only_prediction"]
       .pivot_table(index="city_id", columns="buffer_m", values="sample_count"))
display(tgt)
constant = tgt.nunique(axis=1).eq(1).all()
print(f"target set identical across all three buffers: {constant}")
print("=> differences across buffers measure distance dependence, not resampling.")

target set identical across all three buffers: True
=> differences across buffers measure distance dependence, not resampling.

"""Gain asymmetry, broad QA, primary 2550 m buffer."""
gb = pd.read_csv(ARTIFACTS["heldout_broad_gain"])
gp = gb[(gb.buffer_m == 2550)].copy()
gp["city"] = gp.city_id.map(CITY_LABEL)

tbl = gp.pivot_table(index=["city", "proxy"], columns="gain_stratum",
                     values="structural_minus_neighbors_absolute_error")
tbl = tbl[["<0.8", "0.8-1.25", ">=1.25"]]
display(tbl.round(3))
print("Negative = structural prediction beats neighbors-only. Units: radiance (nW/cm2/sr).")

Negative = structural prediction beats neighbors-only. Units: radiance (nW/cm2/sr).

"""Figure: gain asymmetry as small multiples.

Small multiples rather than grouped bars: the message is that all four
combinations share one shape, and color is already spent on polarity
(improve/harm), so series identity needs its own panel rather than a
fifth hue.
"""
strata = ["<0.8", "0.8-1.25", ">=1.25"]
combos = [(c, p) for c in ["usa_new_york", "india_delhi"]
          for p in ["built_form_primary", "s2_only_ablation"]]

lim = float(np.abs(gp.structural_minus_neighbors_absolute_error).max()) * 1.28
fig, axes = plt.subplots(1, 4, figsize=(11, 3.9), sharey=True)
xs = np.arange(len(strata))

for ax, (city, proxy) in zip(axes, combos):
    sub = gp[(gp.city_id == city) & (gp.proxy == proxy)].set_index("gain_stratum")
    vals = np.array([sub.loc[s, "structural_minus_neighbors_absolute_error"]
                     for s in strata])
    colors = [COLOR["improve"] if v < 0 else COLOR["harm"] for v in vals]
    ax.bar(xs, vals, 0.66, color=colors, edgecolor="white", linewidth=1.4)
    for x, v in zip(xs, vals):
        ax.annotate(f"{v:+.1f}", (x, v), ha="center",
                    va="top" if v < 0 else "bottom",
                    xytext=(0, -4 if v < 0 else 4), textcoords="offset points",
                    fontsize=8, color=INK_2)
    ax.axhline(0, color=INK, lw=1.0)
    ax.set_xticks(xs)
    ax.set_xticklabels(["< 0.8", "0.8-\n1.25", ">= 1.25"], fontsize=8)
    ax.set_title(f"{CITY_LABEL[city]}\n{'built form' if 'built' in proxy else 'S2 only'}",
                 fontsize=9)
    ax.set_ylim(-lim, lim)
    ax.grid(axis="y"); ax.set_axisbelow(True)
    ax.spines["left"].set_visible(False)
    ax.tick_params(axis="y", length=0)

axes[0].set_ylabel("structural MAE - neighbors MAE")
handles = [mpl.patches.Patch(color=COLOR["improve"], label="structural prediction better"),
           mpl.patches.Patch(color=COLOR["harm"], label="structural prediction worse")]
fig.legend(handles=handles, loc="lower center", ncol=2, fontsize=8,
           bbox_to_anchor=(0.5, -0.06))
fig.suptitle("Downweighting helps a lot. Upweighting hurts. Same shape in all four "
             "city x proxy combinations.\nx axis: allocation gain stratum "
             "(radiance-blind, computed from structure proxy and kernel only).",
             x=0.5, y=1.08, fontsize=9.5, color=INK, ha="center")
fig.tight_layout()
plt.savefig(SAVEFIG_PATH / f'gain_asymmetry.png', dpi=600, bbox_inches='tight')
plt.close(fig)

"""Noise simulation: can ratio noise alone reproduce the observed pattern?

Design mirrors the actual held-out estimator, at the coarse-cell scale where it
operates. Known truth L; a proxy rho = L * lognormal noise -- i.e. a proxy whose
ONLY defect is multiplicative noise, with no structural information missing.
Gain, strata, and the structural predictor are computed exactly as in the real
run. The kernel is held fixed; only rho is perturbed.

Second arm: a heavy right-tailed truth field, testing whether proxy tail shape
(built form is a sqrt-transformed building fraction plus saturating road density)
is a competing explanation.

Errors are reported as PERCENT CHANGE relative to the neighbors-only baseline in
the same stratum, so simulated and observed are on a common scale.
"""
from scipy import ndimage as _ndi

rng = np.random.default_rng(20260803)
N_CELLS = 140


def smooth_noise(n, scale, rng):
    a = _ndi.gaussian_filter(rng.normal(size=(n, n)), scale, mode="reflect")
    return a / a.std()


def neighbor_mean(c):
    ker = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]], dtype=float)
    return _ndi.convolve(c, ker / ker.sum(), mode="reflect")


def simulate(arm: str, sigma: float, n_rep: int = 8) -> pd.DataFrame:
    """Return the two summary statistics the real run also reports:
    overall percent change in MAE vs the neighbors-only baseline, and the
    dispersion of the gain field (share of cells outside the middle stratum).
    """
    rows = []
    for _ in range(n_rep):
        base = np.exp(1.1 * smooth_noise(N_CELLS, 4.0, rng))
        if arm == "heavy-tailed proxy":
            base = base ** 2.2
        L = 20.0 + 180.0 * base / base.mean()
        rho = L * np.exp(sigma * smooth_noise(N_CELLS, 0.8, rng))

        nb_L = neighbor_mean(L)
        gain = rho / np.maximum(neighbor_mean(rho), 1e-9)
        structural = nb_L * gain

        err_s, err_n = np.abs(structural - L), np.abs(nb_L - L)
        rows.append({
            "arm": arm, "sigma": sigma,
            "pct_delta": 100 * (err_s.mean() - err_n.mean()) / err_n.mean(),
            "dispersion": 100 * float(((gain < 0.8) | (gain >= 1.25)).mean()),
        })
    return (pd.DataFrame(rows)
            .groupby(["arm", "sigma"], as_index=False)[["pct_delta", "dispersion"]]
            .mean())


SIGMAS = [0.02, 0.05, 0.08, 0.12, 0.16, 0.20, 0.25, 0.30, 0.40, 0.50]
sim = pd.concat([simulate(arm, s) for s in SIGMAS
                 for arm in ("ratio noise only", "heavy-tailed proxy")],
                ignore_index=True)

# observed: same two summary statistics, from the broad-QA held-out run
obs_rows = []
for city in CITIES:
    for proxy in ("built_form_primary", "s2_only_ablation"):
        g = gb[(gb.buffer_m == 2550) & (gb.city_id == city) & (gb.proxy == proxy)]
        n = g.sample_count.sum()
        o = hb[(hb.buffer_m == 2550) & (hb.city_id == city) & (hb.proxy == proxy) &
               (hb.fold.astype(str) == "all") & (hb.radiance_decile.astype(str) == "all")]
        nb_mae = o[o.method == "neighbors_only_prediction"].mae.iloc[0]
        st_mae = o[o.method == "structural_prediction"].mae.iloc[0]
        obs_rows.append({
            "label": f"{CITY_LABEL[city]} {'built form' if 'built' in proxy else 'S2 only'}",
            "proxy": proxy,
            "pct_delta": 100 * (st_mae - nb_mae) / nb_mae,
            "dispersion": 100 * g[g.gain_stratum != "0.8-1.25"].sample_count.sum() / n,
        })
observed = pd.DataFrame(obs_rows)

display(sim.pivot_table(index="sigma", columns="arm",
                        values=["pct_delta", "dispersion"]).round(1))
print("\nOBSERVED:")
display(observed.round(1))

OBSERVED:

"""Figure: the noise-only trajectory cannot reach the observed operating point."""
fig, ax = plt.subplots(figsize=(9.2, 5.4))

# shade the quadrant no noise-only proxy can reach
ax.axhspan(-30, 0, xmin=0.0, xmax=1.0, color=BLUE, alpha=0.05, zorder=0)

# Both simulated arms are the SAME model with one assumption changed, so they
# share a hue and separate by line style + marker. This also keeps the on-screen
# palette to categorical slots 1-3 (blue / orange / aqua), the three that clear
# the all-pairs CVD floors -- avoiding the documented orange-yellow failure.
for arm, ls, marker in [("ratio noise only", "-", "o"),
                        ("heavy-tailed proxy", "--", "s")]:
    s = sim[sim.arm == arm].sort_values("sigma")
    ax.plot(s.dispersion, s.pct_delta, color=ORANGE, ls=ls, lw=2, marker=marker,
            ms=5, mec="white", mew=1.0, label=f"simulated: {arm}", zorder=2)

# annotate sigma on one curve only, to the right of the marker
s = sim[sim.arm == "heavy-tailed proxy"].sort_values("sigma")
for _, r in s.iterrows():
    if r.sigma in (0.20, 0.50):
        ax.annotate(f"$\\sigma$={r.sigma:g}", (r.dispersion, r.pct_delta),
                    xytext=(8, -2), textcoords="offset points",
                    fontsize=7.5, color=MUTED, va="center")

for _, r in observed.iterrows():
    ax.scatter([r.dispersion], [r.pct_delta], s=120,
               color=COLOR[r.proxy], edgecolor="white", linewidth=1.6, zorder=4)
    ax.annotate(r.label, (r.dispersion, r.pct_delta), xytext=(0, -13),
                textcoords="offset points", fontsize=8, color=INK_2,
                ha="center", va="top")

ax.axhline(0, color=INK, lw=1.0, zorder=1)
ax.set_xlabel("gain dispersion: cells outside the 0.8-1.25 stratum (%)")
ax.set_ylabel("overall change in MAE vs neighbors-only (%)")
ax.set_title("A noise-only proxy cannot be both widely dispersed and accurate")
ax.set_yscale("symlog", linthresh=25)
ax.set_yticks([-20, 0, 25, 100, 400])
ax.set_yticklabels(["-20  (better)", "0", "+25", "+100", "+400  (worse)"])
ax.set_xlim(-3, 74)
ax.set_ylim(-32, 900)
ax.annotate("observed proxies live in the shaded quadrant:\n"
            "wide gain dispersion AND lower error.\n"
            "No noise level reaches it.",
            (71, 220), fontsize=8.5, color=BLUE, ha="right", va="center",
            linespacing=1.5)
ax.grid(True); ax.set_axisbelow(True)
ax.legend(loc="upper left", fontsize=8)
fig.tight_layout()
plt.savefig(SAVEFIG_PATH / "noise-only-trajectory.png", dpi=300, bbox_inches="tight")
plt.close(fig)

"""Dose-response: structural advantage as a function of exclusion buffer."""
def advantage(frame: pd.DataFrame) -> pd.DataFrame:
    o = frame[(frame.fold.astype(str) == "all") &
              (frame.radiance_decile.astype(str) == "all")]
    p = o.pivot_table(index=["city_id", "proxy", "buffer_m"],
                      columns="method", values=["mae", "spearman"])
    out = pd.DataFrame({
        "delta_mae": p[("mae", "structural_prediction")] - p[("mae", "neighbors_only_prediction")],
        "delta_spearman": p[("spearman", "structural_prediction")] - p[("spearman", "neighbors_only_prediction")],
    }).reset_index()
    return out

adv_broad = advantage(hb)
adv_strict = advantage(pd.read_csv(ARTIFACTS["heldout_strict"]))
adv_broad["city"] = adv_broad.city_id.map(CITY_LABEL)
display(adv_broad.pivot_table(index=["city", "proxy"], columns="buffer_m",
                              values=["delta_mae", "delta_spearman"]).round(3))

"""Figure: dose-response. Two panels, one measure each -- never a dual axis."""
buffers = [1500.0, 2000.0, 2550.0]
series = [("usa_new_york", "built_form_primary"), ("usa_new_york", "s2_only_ablation"),
          ("india_delhi", "built_form_primary"), ("india_delhi", "s2_only_ablation")]

fig, axes = plt.subplots(1, 2, figsize=(10.5, 4.2))
for ax, metric, ylab, title in [
    (axes[0], "delta_mae", "structural MAE - neighbors MAE",
     "Error advantage grows with exclusion distance"),
    (axes[1], "delta_spearman", "structural rho - neighbors rho",
     "Rank advantage grows too"),
]:
    for city, proxy in series:
        sub = adv_broad[(adv_broad.city_id == city) & (adv_broad.proxy == proxy)]
        sub = sub.set_index("buffer_m").reindex(buffers)
        color = COLOR[proxy]
        ls = "-" if city == "usa_new_york" else "--"
        ax.plot(buffers, sub[metric], color=color, ls=ls, lw=2,
                marker="o", ms=6, mec="white", mew=1.2)
        ax.annotate(f"{CITY_LABEL[city]} {'built form' if 'built' in proxy else 'S2 only'}",
                    (buffers[-1], sub[metric].iloc[-1]), xytext=(6, 0),
                    textcoords="offset points", va="center", fontsize=7.5, color=INK_2)
    ax.axhline(0, color=MUTED, lw=1.0, ls=":")
    ax.set_xticks(buffers)
    ax.set_xticklabels([f"{int(b)} m" for b in buffers])
    ax.set_xlabel("held-out training-exclusion buffer")
    ax.set_ylabel(ylab)
    ax.set_title(title)
    ax.grid(axis="y"); ax.set_axisbelow(True)
    ax.set_xlim(1400, 3050)
fig.suptitle("Broad-QA held-out, identical target set at every buffer. "
             "Leakage would predict decay toward zero.",
             x=0.5, y=1.02, fontsize=9.5, color=INK_2)
fig.tight_layout()
plt.savefig(SAVEFIG_PATH / f'dose-response.png', dpi=600, bbox_inches='tight')
plt.close(fig)

"""The dose-response numbers, broad primary with the strict footnote."""
def dose_row(frame, city, proxy):
    s = frame[(frame.city_id == city) & (frame.proxy == proxy)].set_index("buffer_m")
    return {
        "city": CITY_LABEL[city],
        "proxy": "built form" if "built" in proxy else "S2 only",
        **{f"dMAE@{int(b)}": round(s.loc[b, "delta_mae"], 3) for b in buffers},
        **{f"drho@{int(b)}": round(s.loc[b, "delta_spearman"], 3) for b in buffers},
    }

rows = [dose_row(adv_broad, "usa_new_york", "built_form_primary"),
        dose_row(adv_broad, "india_delhi", "built_form_primary")]
display(pd.DataFrame(rows))

"""Degradation of each method between the 1500 m and 2550 m buffers."""
def degradation(frame):
    o = frame[(frame.fold.astype(str) == "all") & (frame.radiance_decile.astype(str) == "all")]
    p = o.pivot_table(index=["city_id", "proxy", "method"], columns="buffer_m", values="mae")
    d = (p[2550.0] - p[1500.0]).rename("degradation").reset_index()
    return d

deg = degradation(hb)
rows = []
for city in CITIES:
    for proxy in ["built_form_primary", "s2_only_ablation"]:
        n = deg[(deg.city_id == city) & (deg.proxy == proxy) &
                (deg.method == "neighbors_only_prediction")].degradation.iloc[0]
        s = deg[(deg.city_id == city) & (deg.proxy == proxy) &
                (deg.method == "structural_prediction")].degradation.iloc[0]
        rows.append({"city": CITY_LABEL[city],
                     "proxy": "built form" if "built" in proxy else "S2 only",
                     "neighbors_only_degradation": round(n, 3),
                     "structural_degradation": round(s, 3),
                     "skill_retained_pct": round(100 * (1 - s / n), 1)})
display(pd.DataFrame(rows))
print("Skill retained = 1 - (structural degradation / neighbors degradation).")
print("Higher means the structural prior substitutes better for withdrawn observations.")

Skill retained = 1 - (structural degradation / neighbors degradation).
Higher means the structural prior substitutes better for withdrawn observations.

"""Headline held-out table at the 2550 m primary buffer."""
prim = hb[(hb.buffer_m == 2550) & (hb.fold.astype(str) == "all") &
          (hb.radiance_decile.astype(str) == "all")]
wide = prim.pivot_table(index=["city_id", "proxy"], columns="method",
                        values=["mae", "rmse", "spearman"])

folds = hb[(hb.buffer_m == 2550) & (hb.fold.astype(str) != "all") &
           (hb.radiance_decile.astype(str) == "all")]
fp = folds.pivot_table(index=["city_id", "proxy", "fold"], columns="method", values="mae")
improved = ((fp["structural_prediction"] < fp["neighbors_only_prediction"])
            .groupby(level=[0, 1]).sum())

head = pd.DataFrame({
    "city": [CITY_LABEL[c] for c, _ in wide.index],
    "proxy": ["built form" if "built" in p else "S2 only" for _, p in wide.index],
    "neighbors MAE": wide[("mae", "neighbors_only_prediction")].values,
    "structural MAE": wide[("mae", "structural_prediction")].values,
    "neighbors RMSE": wide[("rmse", "neighbors_only_prediction")].values,
    "structural RMSE": wide[("rmse", "structural_prediction")].values,
    "structural rho": wide[("spearman", "structural_prediction")].values,
    "folds improved": [f"{improved.loc[i]}/5" for i in wide.index],
})
display(head.round(3))

"""Over-water allocation vs direct upsample. One metric basis, stated in section 8e."""
wb = pd.read_csv(ARTIFACTS["water_broad"])
ws = pd.read_csv(ARTIFACTS["water_strict"])

rows = []
for city in CITIES:
    r = wb[(wb.city_id == city) & (wb.region == "mapped_water")].iloc[0]
    i = wb[(wb.city_id == city) & (wb.region == "matched_inland_low_proxy")].iloc[0]
    rows.append({"city": CITY_LABEL[city], "proxy": "built form", "QA contract": "broad",
                 "water_reduction_pct": 100 * r.reduction_fraction_vs_direct,
                 "inland_control_reduction_pct": 100 * i.reduction_fraction_vs_direct})
    rs = ws[(ws.city_id == city) & (ws.proxy == "s2_only_ablation") &
            (ws.region == "mapped_water")].iloc[0]
    isx = ws[(ws.city_id == city) & (ws.proxy == "s2_only_ablation") &
             (ws.region == "matched_inland_low_proxy")].iloc[0]
    rows.append({"city": CITY_LABEL[city], "proxy": "S2 only", "QA contract": "strict",
                 "water_reduction_pct": 100 * rs.reduction_fraction_vs_direct,
                 "inland_control_reduction_pct": 100 * isx.reduction_fraction_vs_direct})

water = pd.DataFrame(rows)
water["water_minus_inland_pp"] = (water.water_reduction_pct
                                  - water.inland_control_reduction_pct)
display(water.round(1))

"""The uniform null: isolating structure from smoothing.

Shares are computed on the shoreline table's own basis (all valid pixels), which
is a DIFFERENT denominator from the reduction fractions above. Both the absolute
change and the share change are reported so the basis is never ambiguous.
"""
sl = pd.read_csv(ARTIFACTS["shoreline_summary"])
rows = []
for city in CITIES:
    d = sl[(sl.city_id == city) & (sl.configuration == "direct_upsample")].iloc[0]
    u = sl[(sl.city_id == city) &
           (sl.configuration == "uniform_normalized_convolution")].iloc[0]
    b = sl[(sl.city_id == city) &
           (sl.configuration ==
            "built_form_primary__no_water_prior__circular_mean_reference")].iloc[0]
    for name, r in [("uniform null", u), ("built form", b)]:
        rows.append({
            "city": CITY_LABEL[city], "method": name,
            "over_water_abs_change_pct":
                100 * (r.allocated_radiance_over_water
                       / d.allocated_radiance_over_water - 1),
            "over_water_share_change_pct":
                100 * (r.allocated_radiance_over_water_share
                       / d.allocated_radiance_over_water_share - 1),
        })
display(pd.DataFrame(rows).round(1))
print("Basis: shoreline summary.csv, all valid pixels. Strict contract "
      "(the full water factorial was\nnever re-run under broad). Not "
      "interchangeable with the reduction fractions in the previous cell.")

Basis: shoreline summary.csv, all valid pixels. Strict contract (the full water factorial was
never re-run under broad). Not interchangeable with the reduction fractions in the previous cell.

"""Figure: water is not special. The inland low-proxy control reduces more."""
fig, ax = plt.subplots(figsize=(8.5, 3.8))
sub = water[water.proxy == "built form"]
xs = np.arange(len(sub))
w = 0.36
ax.bar(xs - w / 2, sub.water_reduction_pct, w * 0.92, color=BLUE,
       edgecolor="white", linewidth=1.2, label="mapped water (Overture reference)")
ax.bar(xs + w / 2, sub.inland_control_reduction_pct, w * 0.92, color=ORANGE,
       edgecolor="white", linewidth=1.2, label="matched inland low-proxy control")
for x, (_, r) in zip(xs, sub.iterrows()):
    ax.annotate(f"{r.water_reduction_pct:.1f}%", (x - w / 2, r.water_reduction_pct),
                ha="center", va="bottom", xytext=(0, 3), textcoords="offset points",
                fontsize=8, color=INK_2)
    ax.annotate(f"{r.inland_control_reduction_pct:.1f}%",
                (x + w / 2, r.inland_control_reduction_pct),
                ha="center", va="bottom", xytext=(0, 3), textcoords="offset points",
                fontsize=8, color=INK_2)
ax.set_xticks(xs); ax.set_xticklabels(sub.city)
ax.set_ylabel("reduction vs direct upsample (%)")
ax.set_title("Reduction over water is not larger than over comparable dry low-proxy land")
ax.legend(loc="upper right")
ax.grid(axis="y"); ax.set_axisbelow(True)
fig.tight_layout()
plt.savefig(SAVEFIG_PATH / "water-control.png", dpi=300, bbox_inches="tight")
plt.close(fig)

"""Four independent diagnostics for Delhi S2-only."""
g0 = pd.read_csv(ARTIFACTS["gate0_s2"])
delhi_g0 = g0[g0.get("city_id", g0.columns[0]).astype(str).str.contains("delhi", case=False)]

heldout_delhi = head[(head.city == "Delhi") & (head.proxy == "S2 only")].iloc[0]
gain_delhi = gp[(gp.city_id == "india_delhi") & (gp.proxy == "s2_only_ablation") &
                (gp.gain_stratum == ">=1.25")].iloc[0]
water_delhi = water[(water.city == "Delhi") & (water.proxy == "S2 only")].iloc[0]

evidence = pd.DataFrame([
    {"#": 1, "diagnostic": "Gate 0 shape (strict)",
     "value": "mound-shaped; passes all positive-correlation thresholds",
     "verdict": "missed by the rule"},
    {"#": 2, "diagnostic": "Held-out MAE, 2550 m (broad)",
     "value": f"{heldout_delhi['structural MAE']:.3f} vs "
              f"{heldout_delhi['neighbors MAE']:.3f} baseline "
              f"({heldout_delhi['folds improved']} folds)",
     "verdict": "worse than interpolation"},
    {"#": 3, "diagnostic": "Gain asymmetry, >=1.25 stratum (broad)",
     "value": f"{gain_delhi.structural_minus_neighbors_absolute_error:+.3f}",
     "verdict": "steepest of all four combinations"},
    {"#": 4, "diagnostic": "Over-water allocation vs direct (strict)",
     "value": f"{water_delhi.water_reduction_pct:+.1f}% reduction "
              f"(i.e. {-water_delhi.water_reduction_pct:+.1f}% MORE radiance on water)",
     "verdict": "wrong direction"},
])
display(evidence)

"""The final evidence classification, verbatim."""
ev = load_json(ARTIFACTS["evidence_v2"])

print("CLASSIFICATION:", ev["classification"].upper().replace("_", " "))
print("\n" + "=" * 78)
print("EXPERIMENT CONTRACT")
print("=" * 78)
for k, v in ev["experiment_contract"].items():
    print(f"\n{k}:\n  {v}")

print("\n" + "=" * 78)
print("SUPPORTED CLAIMS")
print("=" * 78)
for c in ev["supported_claims"]:
    print(f"  [+] {c}")

print("\n" + "=" * 78)
print("UNSUPPORTED CLAIMS -- language that may not be used")
print("=" * 78)
for c in ev["unsupported_claims"]:
    print(f"  [-] {c}")

print("\n" + "=" * 78)
print("MIXED FINDINGS")
print("=" * 78)
for c in ev["mixed_findings"]:
    print(f"  [~] {c}")

print("\n" + "=" * 78)
print("INHERITED EVIDENCE")
print("=" * 78)
for c in ev["inherited_evidence"]:
    print(f"  [<] {c}")

CLASSIFICATION: COVERAGE COMPLETE BOUNDED METHOD RESULT

==============================================================================
EXPERIMENT CONTRACT
==============================================================================

fine_grid_limit:
  the 10 m products are structurally allocated observation fields, not independently observed or calibrated 10 m radiance

objective:
  test whether a static fine-grid structural proxy can allocate coarse VNP46A2 radiance more plausibly than declared null allocations while VNP46A2 remains the radiometric authority

primary_estimand:
  native-cell broad-QA held-out prediction error under the frozen 2550 m PSF-support exclusion

==============================================================================
SUPPORTED CLAIMS
==============================================================================
  [+] all v2 primary and baseline COGs have complete two-city spatial support
  [+] the New York S2-only broad-QA sensitivity has zero radiance-neighborhood gaps; remaining invalid pixels are unchanged S2 proxy-support gaps
  [+] on the frozen native-cell cohort, broad-QA built form improves overall MAE and RMSE in both cities at the 2550 m buffer
  [+] the radiance-blind gain asymmetry persists under broad QA in all four city-by-proxy comparisons
  [+] the broad water-versus-inland result remains consistent with general low-proxy reallocation rather than a water-specific effect

==============================================================================
UNSUPPORTED CLAIMS -- language that may not be used
==============================================================================
  [-] observed or calibrated 10 m nighttime radiance
  [-] pixel-level accuracy from gain strata
  [-] equivalent quality of MQF 0 and MQF 1
  [-] broad-QA daily observation-condition correction
  [-] validated water correction or recovered physical emissions

==============================================================================
MIXED FINDINGS
==============================================================================
  [~] Delhi built form improves broad-QA MAE in 4 of 5 folds; fold 3 worsens by only 0.044 radiance units
  [~] New York S2-only remains lower-error than built form, while Delhi S2-only remains negative overall
  [~] broad QA removes support holes but admits MQF 1 and lower cloud-mask quality
  [~] only the New York S2-only spatial raster was regenerated under broad QA; the Delhi S2-only spatial sensitivity remains strict-v1 and is not a like-for-like cross-city map comparison

==============================================================================
INHERITED EVIDENCE
==============================================================================
  [<] the daily observation-condition propagation audit remains strict-QA conservative sensitivity evidence and is not relabeled as broad-QA validation

"""Reproduction check against the private-runtime frozen configuration."""
FROZEN_CONFIG_SHA256 = "9a0338232881f7d6353259672d6a3c45f83b50d67dacbc00f3d8d986b72d718a"
config_bytes = CONFIG_PATH.read_bytes()
print(f"public redacted config sha256: {hashlib.sha256(config_bytes).hexdigest()}")
print(f"frozen runtime config sha256: {FROZEN_CONFIG_SHA256}")

manifests = {
    "full-city broad (v2)": BROAD / "artifact_manifest.json",
    "full-city strict (v1)": STRICT / "artifact_manifest.json",
    "held-out broad (v3)": VAL / "gate2/heldout/v3_physics_buffered_native_cell_broad_qa/manifest.json",
    "shoreline (v1)": VAL / "gate2/shoreline/v1_overture_mapped_water/manifest.json",
    "sensitivity (v1)": VAL / "gate2/sensitivity/v1_full_city_sensitivity/manifest.json",
    "observation conditions (v1)": VAL / "gate2/observation_conditions/v1_daily_coarse_operator/manifest.json",
}
rows = []
for name, path in manifests.items():
    if not path.exists():
        rows.append({"manifest": name, "configuration_sha256": "MISSING"}); continue
    blob = json.dumps(load_json(path))
    import re
    hits = sorted(set(re.findall(r"\b[0-9a-f]{64}\b", blob)))
    cfg = [h for h in hits if h == FROZEN_CONFIG_SHA256]
    rows.append({"manifest": name,
                 "configuration_sha256": (cfg[0][:16] + "...") if cfg else "not found"})
display(pd.DataFrame(rows))

public redacted config sha256: 67b984311cda71bd275afaff47ad517649e857e9bd1b5959d664945b3e7a1796
frozen runtime config sha256: 9a0338232881f7d6353259672d6a3c45f83b50d67dacbc00f3d8d986b72d718a

	city	consistency_MAE	consistency_RMSE	consistency_bias	mean_abs_source_radiance	MAE_pct_of_mean
0	New York	5.342	9.907	0.559	43.638	12.242
1	Delhi	3.885	6.379	0.572	37.222	10.437

	city_id	allocation_proxy	sample_count	input_sample_count	excluded_s2_coarse_support_fraction	minimum_s2_coarse_support_fraction	block_count	citywide_spearman	block_aggregate_spearman	block_mean_spearman	...	continuous_resampling	categorical_resampling	s2_composite_decision_id	s2_composite_method	s2_minimum_common_valid_observations	s2_coarse_support_decision_id	gate0_aggregation_halo_m	gate0_status	overture_release	built_form_temporal_semantics
0	usa_new_york	built_form_primary	15335	15335	0.0	0.9	121	0.811875	0.826067	0.826067	...	bilinear	nearest	S2-COMPOSITE-001	index_of_bandwise_temporal_medians_after_datat...	5	S2-SUPPORT-001	1000	current_corrected_rerun	2026-07-22.0	contemporary_structural_proxy_used_to_allocate...
1	india_delhi	built_form_primary	13280	13280	0.0	0.9	121	0.761388	0.900982	0.900982	...	bilinear	nearest	S2-COMPOSITE-001	index_of_bandwise_temporal_medians_after_datat...	5	S2-SUPPORT-001	1000	current_corrected_rerun	2026-07-22.0	contemporary_structural_proxy_used_to_allocate...

	city_id	allocation_proxy	sample_count	input_sample_count	excluded_s2_coarse_support_fraction	minimum_s2_coarse_support_fraction	block_count	citywide_spearman	block_aggregate_spearman	block_mean_spearman	...	vnp_quality_decision_id	resampling_decision_id	continuous_resampling	categorical_resampling	s2_composite_decision_id	s2_composite_method	s2_minimum_common_valid_observations	s2_coarse_support_decision_id	gate0_aggregation_halo_m	gate0_status
0	usa_new_york	s2_only_ablation	14963	15335	0.024258	0.9	121	0.783180	0.811195	0.811195	...	VNP-QA-001	RESAMPLE-001	bilinear	nearest	S2-COMPOSITE-001	index_of_bandwise_temporal_medians_after_datat...	5	S2-SUPPORT-001	1000	current_corrected_rerun
1	india_delhi	s2_only_ablation	13280	13280	0.000000	0.9	121	0.437784	0.575816	0.575816	...	VNP-QA-001	RESAMPLE-001	bilinear	nearest	S2-COMPOSITE-001	index_of_bandwise_temporal_medians_after_datat...	5	S2-SUPPORT-001	1000	current_corrected_rerun

		delta_mae			delta_spearman
	buffer_m	1500.0	2000.0	2550.0	1500.0	2000.0	2550.0
city	proxy
Delhi	built_form_primary	-1.071	-1.357	-1.347	0.057	0.066	0.069
Delhi	s2_only_ablation	0.322	0.228	0.161	0.017	0.018	0.014
New York	built_form_primary	-1.618	-1.992	-2.617	0.077	0.097	0.109
New York	s2_only_ablation	-1.988	-2.456	-2.929	0.061	0.076	0.088

	city	proxy	dMAE@1500	dMAE@2000	dMAE@2550	drho@1500	drho@2000	drho@2550
0	New York	built form	-1.618	-1.992	-2.617	0.077	0.097	0.109
1	Delhi	built form	-1.071	-1.357	-1.347	0.057	0.066	0.069

Measuring lights at night¶

0. Reporting contract¶

1. Objective¶

2. The operator¶

Properties of the kernel¶

3. The Fork et al. analogy, and its two breaks¶

4. Data, construction, and how the choices were governed¶

Governance¶

5. The product¶

A failure cutout¶

6. Gate 0: does the prior track the observation at all?¶

Delhi S2-Only is a confident miss (bad)¶

7. Held-out validation design¶

8. Results¶

8a. The gain asymmetry: lead with the mechanism¶

Is this just noise?¶

8b. Dose–response across buffers¶

8c. The substitution mechanism¶

8d. The headline held-out result¶

8e. Water: an out-of-model check¶

8f. Delhi S2-only: a worked failure¶

8g. A couple notes¶

9. Data and experiment contracts¶

The scope/limits of this experiment¶

10. Future work¶

11. Data provenance, licensing, and reproduction¶

Provenance¶

Reproduction¶

	key	exists	size_mb
0	config	True	0.0
1	selection_v2_broad	True	0.0
2	selection_v1_strict	True	0.0
3	evidence_v2	True	0.0
4	headline_v2	True	0.0
5	headline_v1	True	0.0
6	heldout_broad	True	0.1
7	heldout_broad_gain	True	0.0
8	heldout_broad_retain	True	0.0
9	heldout_strict	True	0.1
10	heldout_strict_gain	True	0.0
11	water_broad	True	0.0
12	water_strict	True	0.0
13	water_compare_strict	True	0.0
14	shoreline_summary	True	0.0
15	shoreline_strata	True	0.1
16	kernel_sensitivity	True	0.0
17	obs_conditions	True	0.0
18	gate1_closeout	True	0.0
19	gate0_builtform	True	0.0
20	gate0_s2	True	0.0
21	gate0_regen	True	0.0
22	manifest_broad	True	0.0
23	primary_usa_new_york	True	275.4
24	direct_usa_new_york	True	13.8
25	uniform_usa_new_york	True	272.4
26	trust_usa_new_york	True	105.6
27	s2_usa_new_york	True	276.8
28	structure_usa_new_york	True	41.5
29	primary_india_delhi	True	277.6
30	direct_india_delhi	True	13.2
31	uniform_india_delhi	True	273.5
32	trust_india_delhi	True	127.6
33	s2_india_delhi	True	277.6
34	structure_india_delhi	True	35.6

	band	name	role
0	1	allocation	the product
1	2	operator_consistency_error	reaggregation residual (section 2)
2	3	radiance_support_fraction	support
3	4	proxy_support_fraction	support
4	5	allocation_support_fraction	support
5	6	geometric_support_fraction	support
6	7	valid_output_mask	validity
7	8	operator_consistency_valid_mask	validity
8	9	boundary_mask	validity
9	10	invalid_radiance_neighborhood_mask	invalid-input mask
10	11	invalid_proxy_neighborhood_mask	invalid-input mask
11	12	denominator_floor_mask	numerical guard
12	13	denominator_instability_mask	numerical guard
13	14	negative_radiance_input_mask	numerical guard
14	15	negative_preclip_output_mask	numerical guard
15	16	insufficient_proxy_support_mask	numerical guard

buffer_m	1500.0	2000.0	2550.0
city_id
india_delhi	13280.0	13280.0	13280.0
usa_new_york	15335.0	15335.0	15335.0

	dispersion		pct_delta
arm	heavy-tailed proxy	ratio noise only	heavy-tailed proxy	ratio noise only
sigma
0.02	0.9	0.0	-88.3	-67.1
0.05	0.9	0.0	-73.4	-19.6
0.08	1.0	0.0	-49.5	27.4
0.12	1.4	0.0	-31.9	107.3
0.16	2.2	0.3	-5.2	167.0
0.20	3.9	1.5	11.7	244.3
0.25	7.4	5.0	42.1	305.5
0.30	12.7	10.3	65.0	383.0
0.40	24.0	22.2	144.0	574.9
0.50	34.9	33.6	171.2	741.3

	label	proxy	pct_delta	dispersion
0	New York built form	built_form_primary	-14.0	54.1
1	New York S2 only	s2_only_ablation	-15.7	47.3
2	Delhi built form	built_form_primary	-10.0	63.3
3	Delhi S2 only	s2_only_ablation	1.2	50.0

	city	proxy	neighbors_only_degradation	structural_degradation	skill_retained_pct
0	New York	built form	1.671	0.672	59.8
1	New York	S2 only	1.671	0.730	56.3
2	Delhi	built form	0.810	0.534	34.1
3	Delhi	S2 only	0.810	0.649	19.8

	city	proxy	neighbors MAE	structural MAE	neighbors RMSE	structural RMSE	structural rho	folds improved
0	Delhi	built form	13.493	12.147	19.615	18.660	0.857	4/5
1	Delhi	S2 only	13.493	13.655	19.615	20.256	0.802	2/5
2	New York	built form	18.717	16.100	30.433	28.199	0.879	5/5
3	New York	S2 only	18.717	15.787	30.433	27.539	0.857	5/5

	city	proxy	QA contract	water_reduction_pct	inland_control_reduction_pct	water_minus_inland_pp
0	New York	built form	broad	39.6	81.0	-41.4
1	New York	S2 only	strict	6.8	68.6	-61.7
2	Delhi	built form	broad	52.6	55.0	-2.4
3	Delhi	S2 only	strict	-19.0	70.5	-89.5

	city	method	over_water_abs_change_pct	over_water_share_change_pct
0	New York	uniform null	3.3	4.3
1	New York	built form	-42.0	-40.7
2	Delhi	uniform null	4.8	4.8
3	Delhi	built form	-52.8	-52.1

	#	diagnostic	value	verdict
0	1	Gate 0 shape (strict)	mound-shaped; passes all positive-correlation ...	missed by the rule
1	2	Held-out MAE, 2550 m (broad)	13.655 vs 13.493 baseline (2/5 folds)	worse than interpolation
2	3	Gain asymmetry, >=1.25 stratum (broad)	+3.407	steepest of all four combinations
3	4	Over-water allocation vs direct (strict)	-19.0% reduction (i.e. +19.0% MORE radiance on...	wrong direction

	band	name
0	1	fine_grid_allocation_gain
1	2	validated_gain_stratum

Input	Source	Version / pin	Notes
Nighttime radiance	NASA Black Marble VNP46A2 (`DNB_BRDF_Corrected_NTL` + QA bands)	2024-01-11 to 2024-04-20	NASA-led mission data are CC0 unless marked otherwise; cite VNP46A2.002. Accessed via Earth Engine. Collection 2. The Oct 2024 NASA C2 guide and the EE catalog disagree about MQF 1, which is why both QA contracts are retained
Daytime structure	Sentinel-2 SR Harmonized + Cloud Score+	same window	Copernicus Sentinel free/full/open terms; modified outputs require the Copernicus source notice. Cloud Score+ is CC BY 4.0. B11/B12 and derived indices retain 20 m effective resolution though sampled at 10 m
Buildings and roads	Overture Maps	`2026-07-22.0`	Buildings and Transportation themes are ODbL and require Overture/OSM attribution plus applicable upstream notices. A contemporary proxy used to allocate 2024 radiance
Water reference (validation)	Overture areal water	`2026-06-17.0`	Base theme is ODbL. Deliberately separate from the proxy's internal prior. Primarily OSM-derived; line-only waterways unbuffered, pools/wastewater excluded
Persistent-water prior (internal)	JRC Global Surface Water occurrence	1984–2021, 30 m	A temporal frequency, not a 10 m spatial fraction and not a 2024 observation. Inactive in the reporting primary (`no_water_prior`)

	manifest	configuration_sha256
0	full-city broad (v2)	9a0338232881f7d6...
1	full-city strict (v1)	9a0338232881f7d6...
2	held-out broad (v3)	9a0338232881f7d6...
3	shoreline (v1)	9a0338232881f7d6...
4	sensitivity (v1)	9a0338232881f7d6...
5	observation conditions (v1)	9a0338232881f7d6...