Winner models, code, notebooks, data and results

2b9ed19 verified 4 days ago

19.2 kB

"""
The adaptation layer: retrieval over whatever labelled target pixels we are
given, plus label-free neighbourhood features.
Three distinct channels, deliberately kept separate so their contributions can
be measured one at a time:
  A. NEIGHBOURHOOD FEATURES - no labels at all. The mean feature vector of the
     surrounding pixels. Buildings are built in blocks, so a pixel's
     neighbourhood says a great deal about its era, and this works at n=0.
  B. SPECTRAL RETRIEVAL - the k labelled pixels in this city whose features are
     closest to the query. Summarised as a class histogram, and/or rendered as
     in-context demonstrations.
  C. SPATIAL RETRIEVAL - the labelled pixels physically nearest the query.
     Weak when only 20 pixels are labelled, strong at 800.
Leakage discipline: when these fields are built for MADRID training rows, the
labelled pool is a RANDOM n-per-class sample drawn from spatial blocks that
exclude the query pixel. That reproduces exactly the situation the model meets
in Amsterdam, and a pixel can never retrieve itself or its own neighbours.
"""
from __future__ import annotations

import numpy as np
from sklearn.neighbors import NearestNeighbors

CLASSES = np.array([1, 2, 3, 4])


# --------------------------------------------------------------------------- #
# A. Label-free neighbourhood features
# --------------------------------------------------------------------------- #
def add_neighbourhood_features(df, coord_cols=("px_key", "py_key"),
                               feature_cols=None, k: int = 8,
                               max_radius: float = 10.0, prefix: str = "nbr_"):
    """
    For each pixel, average a handful of feature columns over its k nearest
    neighbours (excluding itself). Uses no labels, so it is legal in both
    cities and available at n=0.
    `max_radius` is in pixel units (1 = 30 m); neighbours further away are
    ignored and the pixel keeps its own value instead.
    """
    import pandas as pd

    if feature_cols is None:
        feature_cols = [c for c in ("sr_brightness_mean", "sr_NDVI_mean",
                                    "sr_SWIR1_mean", "sr_brightness_late_minus_early",
                                    "step_brightness_year", "step_brightness_mag",
                                    "sr_brightness_std") if c in df.columns]
    coords = df[list(coord_cols)].to_numpy(dtype=float)
    V = df[list(feature_cols)].to_numpy(dtype=float)

    nn = NearestNeighbors(n_neighbors=min(k + 1, len(df))).fit(coords)
    dist, idx = nn.kneighbors(coords)
    dist, idx = dist[:, 1:], idx[:, 1:]                 # drop self
    ok = dist <= max_radius
    w = ok.astype(float)
    denom = np.clip(w.sum(axis=1, keepdims=True), 1e-9, None)

    out = {}
    for j, c in enumerate(feature_cols):
        vals = V[idx, j]
        mean = (vals * w).sum(axis=1, keepdims=True) / denom
        mean = np.where(w.sum(axis=1, keepdims=True) > 0, mean, V[:, [j]])
        out[f"{prefix}{c}"] = mean.ravel()
        out[f"{prefix}{c}_diff"] = V[:, j] - mean.ravel()
    out[f"{prefix}n_within_radius"] = ok.sum(axis=1).astype(float)
    return pd.concat([df, pd.DataFrame(out, index=df.index)], axis=1)


# --------------------------------------------------------------------------- #
# B + C. Retrieval over a labelled pool
# --------------------------------------------------------------------------- #
class LabelledPool:
    """
    A set of labelled pixels in one city (in Amsterdam: the n-per-class sample
    the organisers hand us; in Madrid during training: a simulated sample).
    sklearn only - no FAISS - because the pool is at most a few thousand rows
    and the notebook must run with the internet switched off on Kaggle.
    """

    def __init__(self, X, y, coords, rows=None):
        self.X = np.asarray(X, dtype=np.float64)
        self.y = np.asarray(y, dtype=int)
        self.coords = np.asarray(coords, dtype=float)
        self.rows = rows                                  # for rendering shots
        self._nn_spec = NearestNeighbors(n_neighbors=min(16, len(self.X))).fit(self.X)
        self._nn_spat = NearestNeighbors(n_neighbors=min(16, len(self.X))).fit(self.coords)

    def __len__(self):
        return len(self.X)

    def spectral(self, Xq, k: int = 8):
        k = min(k, len(self.X))
        d, i = self._nn_spec.kneighbors(np.atleast_2d(Xq), n_neighbors=k)
        return d, i

    def spatial(self, coords_q, k: int = 4):
        k = min(k, len(self.coords))
        d, i = self._nn_spat.kneighbors(np.atleast_2d(coords_q), n_neighbors=k)
        return d, i

    def stratified_spectral(self, Xq, per_class: int = 1):
        """Closest `per_class` pool members of each class - a class-balanced
        set of demonstrations, so the prompt cannot be swamped by class 2."""
        picks = []
        for c in CLASSES:
            sel = np.flatnonzero(self.y == c)
            if len(sel) == 0:
                continue
            d = np.linalg.norm(self.X[sel] - np.asarray(Xq)[None, :], axis=1)
            order = np.argsort(d)[:per_class]
            picks += [(float(d[o]), int(sel[o])) for o in order]
        picks.sort()
        return picks


LETTERS = ["A", "B", "C", "D"]


def rag_block(pool: LabelledPool, Xq, coords_q, k_spec: int = 16, k_spat: int = 4,
              style: str = "compact", pixel_m: float = 30.0,
              max_spatial_px: float = 30.0, exclude: int | None = None) -> str:
    """
    Turn retrieval into text.
    exclude : index INTO THE POOL to leave out, for when the query pixel is
              itself a pool member. Without it the block states the answer: the
              identical row gives spectral distance 0.00 to its own class,
              makes itself "the single most similar known building", and shows
              up as "class X at 0 m". That is how the Amsterdam calibration
              shots were being rendered, so the per-class biases were fitted on
              a prompt distribution that never occurs at evaluation time. The
              Madrid training side avoids it structurally, by drawing each
              block's pool from the other blocks (`simulate_pools`,
              `exclude_block=`), but nothing stopped it here.
    The field that carries the information is the DISTANCE TO THE NEAREST
    LABELLED EXAMPLE OF EACH CLASS - not a class histogram. The organisers hand
    out n labelled pixels *per class*, so the pool is balanced by construction
    and a histogram of the k nearest neighbours carries almost no signal while
    actively inviting the model to copy the majority. Per-class distances remove
    that failure mode by construction instead of patching it afterwards.
    A second, independent channel reports the physically nearest labelled
    pixels in metres. The two are deliberately NOT fused into one tuned
    distance: there are no target labels to tune the weighting with, so the
    model is left to learn it.
    """
    if len(pool) == 0:
        return ""
    Xq = np.asarray(Xq, dtype=float).ravel()

    keep = np.ones(len(pool), dtype=bool)
    if exclude is not None:
        keep[int(exclude)] = False
        if not keep.any():
            return ""
    n_kept = int(keep.sum())
    lines = [f"Labelled buildings elsewhere in this city ({n_kept} of them, "
             f"{n_kept // 4} per age class)"]

    # --- channel 1: spectral distance to the nearest example of each class ---
    per_class = {}
    for c in CLASSES:
        sel = np.flatnonzero((pool.y == c) & keep)
        if len(sel) == 0:
            continue
        d = np.linalg.norm(pool.X[sel] - Xq[None, :], axis=1)
        per_class[int(c)] = float(d.min())
    if per_class:
        lines.append("- how far this pixel is from the most similar known "
                     "building of each class (smaller = more similar): "
                     + ", ".join(f"{LETTERS[c - 1]} {v:.2f}"
                                 for c, v in sorted(per_class.items())))
        best = min(per_class, key=per_class.get)
        lines.append(f"- the single most similar known building is class "
                     f"{LETTERS[best - 1]}")

    # --- a soft, distance-weighted vote over the k nearest, per-class normalised
    # Ask for one extra neighbour when a row is excluded, so dropping it still
    # leaves k_spec of them.
    d, i = pool.spectral(Xq, k_spec + (exclude is not None))
    d, i = d[0], i[0]
    if exclude is not None:
        m_keep = i != int(exclude)
        d, i = d[m_keep][:k_spec], i[m_keep][:k_spec]
    tau = float(np.median(d)) + 1e-9
    w = np.exp(-(d / tau) ** 2)
    score = np.zeros(4)
    for c in CLASSES:
        m = pool.y[i] == c
        n_c = max(int(((pool.y == c) & keep).sum()), 1)
        score[c - 1] = w[m].sum() / n_c
    if score.sum() > 0:
        score = score / score.sum()
        lines.append("- similarity-weighted vote of the "
                     f"{len(d)} most similar known buildings: "
                     + ", ".join(f"{LETTERS[k]} {100 * score[k]:.0f}%"
                                 for k in range(4)))

    # --- channel 2: physical distance ---
    ds, isp = pool.spatial(coords_q, k_spat + (exclude is not None))
    ds, isp = ds[0], isp[0]
    if exclude is not None:
        m_keep = isp != int(exclude)
        ds, isp = ds[m_keep][:k_spat], isp[m_keep][:k_spat]
    near = ds <= max_spatial_px
    if near.any():
        lines.append("- nearest known buildings on the ground: "
                     + ", ".join(f"class {LETTERS[int(pool.y[j]) - 1]} at {m * pixel_m:.0f} m"
                                 for j, m in zip(isp[near], ds[near])))
    else:
        lines.append(f"- no known building within "
                     f"{max_spatial_px * pixel_m:.0f} m of this pixel")
    return "\n".join(lines)


# --------------------------------------------------------------------------- #
# Metric space for retrieval
# --------------------------------------------------------------------------- #
def transfer_safe_columns(columns) -> list[str]:
    """
    Columns whose *level* does not encode which city we are in.
    Keep: spectral indices (a multiplicative gain cancels in a normalised
    ratio), every spread/std feature (an additive city offset cancels exactly in
    a standard deviation - the point most teams miss), year-on-year deltas,
    early-minus-late contrasts, slopes, chromaticity, and the change-event
    fields. Drop: raw band means and raw era means, which are precisely the
    "Madrid is brighter than Amsterdam" channel, and the data-availability
    flags, which encode that Amsterdam is cloudier.
    """
    keep, drop_exact = [], {"has_early_data", "has_late_data", "n_observed_years",
                            "first_observed_year", "last_observed_year"}
    # Anything that encodes WHICH INSTRUMENT took the picture. The satellite
    # selected for a given year differs between the two cities in 11 of the 25
    # years from 2001 on, and OLI reads about 0.015 darker than TM/ETM+, so
    # these are city markers dressed up as physics. They stay in the feature
    # table so a model can absorb the nuisance, and stay out of any distance.
    sensor_marks = ("sat_share", "sat_switches", "_oli_minus_tm",
                    "_oli_mean", "_tm_mean", "sat_code")

    # Spectral INDEX names, matched case-sensitively on purpose. "NBR" the
    # index and "nbr_" the neighbourhood prefix are not the same thing, and
    # matching "nbr" case-insensitively conflated them: every single one of the
    # 15 `nbr_*` neighbourhood columns was admitted to the metric, including
    # `nbr_sr_brightness_mean` (city AUC 0.966), `nbr_sr_SWIR1_mean` (0.957)
    # and `nbr_n_within_radius`, a pure urban-density marker. That is exactly
    # the "Madrid is brighter than Amsterdam" channel this function exists to
    # remove, let in through the one clause meant to keep a burn ratio.
    INDEX_NAMES = ("NDVI", "NDBI", "_UI_", "MNDWI", "BSI", "NBR", "_TCG",
                   "_TCW", "chrom")

    def judge(name: str) -> bool:
        """Is this column's LEVEL free of city identity?"""
        lc = name.lower()
        is_index = any(ix in name for ix in INDEX_NAMES) or "chrom" in lc
        # Offset- and gain-cancelling shapes. `_diff` is a pixel minus its own
        # neighbourhood, so any city-wide offset cancels in it.
        safe = ("_std" in lc or lc.startswith("d_") or "_minus_" in lc
                or "_slope" in lc or "chrom" in lc or lc.startswith("step_")
                or "_cv" in lc or lc.endswith("_diff") or is_index)
        # A raw level dressed up as a mean. A difference of two levels is fine.
        unsafe = ("_mean" in lc
                  and not ("_minus_" in lc or "_cv" in lc or lc.endswith("_diff"))
                  and not is_index)
        return safe and not unsafe

    for c in columns:
        if c in drop_exact or any(s in c for s in sensor_marks):
            continue
        # A neighbourhood aggregate inherits the transfer safety of whatever it
        # aggregates: the mean of a raw band over a 7x7 window is still a raw
        # band level, and the mean of NDVI is still a ratio.
        inner = c[4:] if c.lower().startswith("nbr_") else c
        if judge(inner):
            keep.append(c)
    return keep


class WhitenedMetric:
    """
    PCA-whitening fitted on the TARGET city's own unlabelled features.
    Plain Euclidean distance over these features triple-counts absolute
    brightness, because `Blue_mean`, `Blue_early_mean` and `Blue_late_mean` are
    near-duplicates. Whitening on the target city removes that redundancy and
    is the cheap stand-in for a learned metric.
    """

    def __init__(self, n_components: int = 24):
        self.n_components = n_components

    def fit(self, X):
        from sklearn.decomposition import PCA
        from sklearn.preprocessing import QuantileTransformer
        X = np.asarray(X, dtype=float)
        self.qt_ = QuantileTransformer(output_distribution="normal",
                                       n_quantiles=min(1000, len(X)),
                                       random_state=0).fit(X)
        k = min(self.n_components, X.shape[1], len(X))
        self.pca_ = PCA(n_components=k, whiten=True, random_state=0)
        self.pca_.fit(self.qt_.transform(X))
        return self

    def transform(self, X):
        return self.pca_.transform(self.qt_.transform(np.asarray(X, dtype=float)))

    def fit_transform(self, X):
        return self.fit(X).transform(X)


def window_stats(df, coord_cols=("px_key", "py_key"), feature_cols=None,
                 windows=(1, 3), prefix="win"):
    """
    Mean / std of a few features over square pixel windows around each pixel
    (radius 1 = 3x3, radius 3 = 7x7), plus the pixel-minus-window contrast.
    Label-free, so it works at n = 0, and it is a local contrast rather than an
    absolute level, which is exactly the kind of quantity that survives the jump
    to another city.
    """
    import pandas as pd
    from scipy.spatial import cKDTree

    if feature_cols is None:
        feature_cols = [c for c in ("sr_brightness_mean", "sr_NDVI_mean",
                                    "sr_SWIR1_mean", "step_brightness_year",
                                    "sr_brightness_late_minus_early")
                        if c in df.columns]
    coords = df[list(coord_cols)].to_numpy(float)
    V = df[list(feature_cols)].to_numpy(float)
    n, F = V.shape
    tree = cKDTree(coords)
    out = {}
    for r in windows:
        groups = tree.query_ball_point(coords, r=r * 1.4143, p=2)
        # One pass over pixels, all features at once. Looping per feature as
        # well would be F times as many Python iterations - 760k instead of
        # 152k on the real 76k-pixel city.
        mu = np.empty((n, F))
        sd = np.empty((n, F))
        cnt = np.empty(n)
        for i, g in enumerate(groups):
            block = V[g, :]
            mu[i] = block.mean(axis=0)
            sd[i] = block.std(axis=0, ddof=0)
            cnt[i] = len(g)
        w = 2 * r + 1
        for j, c in enumerate(feature_cols):
            out[f"{prefix}{w}_{c}_mean"] = mu[:, j]
            out[f"{prefix}{w}_{c}_std"] = sd[:, j]
            out[f"{prefix}{w}_{c}_contrast"] = V[:, j] - mu[:, j]
        out[f"{prefix}{w}_n"] = cnt
    return pd.concat([df, pd.DataFrame(out, index=df.index)], axis=1)


def shots_block(pool: LabelledPool, Xq, render_fn, per_class: int = 1,
                max_shots: int = 4) -> str:
    """
    In-context demonstrations: the nearest labelled pixel of each class,
    rendered with the same verbaliser and ending in its true answer.
    Class-balanced on purpose; majority-label bias is a real failure mode.
    """
    if len(pool) == 0 or pool.rows is None:
        return ""
    picks = pool.stratified_spectral(Xq, per_class=per_class)[:max_shots]
    parts = ["Worked examples from this city:"]
    for _, j in picks:
        parts.append(render_fn(pool.rows.iloc[j], int(pool.y[j])))
    return "\n\n".join(parts)


# --------------------------------------------------------------------------- #
# Simulating the adaptation regime inside the training city
# --------------------------------------------------------------------------- #
def spatial_blocks(coords, n_blocks: int = 20) -> np.ndarray:
    """
    Tile the city into roughly equal-area blocks by coordinate quantiles.
    Blocks, not random folds: neighbouring pixels share an era, so a random
    split would let a pixel retrieve its own back garden.
    """
    coords = np.asarray(coords, dtype=float)
    side = int(np.ceil(np.sqrt(n_blocks)))
    bx = np.digitize(coords[:, 0], np.quantile(coords[:, 0], np.linspace(0, 1, side + 1)[1:-1]))
    by = np.digitize(coords[:, 1], np.quantile(coords[:, 1], np.linspace(0, 1, side + 1)[1:-1]))
    return bx * side + by


def simulate_pools(y, blocks, n_per_class: int, rng, exclude_block,
                   allowed_blocks=None):
    """
    Draw an n-per-class labelled pool from blocks other than the query's own.
    exclude_block  : one block id, or any iterable of them.
    allowed_blocks : if given, the pool may ONLY come from these blocks. Used
                     to keep the spatial validation holdout out of the pool:
                     without it, `exclude_block=b` alone still let a training
                     block's retrieval text quote the labels of validation-block
                     pixels, which is the one thing the holdout exists to
                     prevent.
    """
    y = np.asarray(y)
    blocks = np.asarray(blocks)
    excl = ({int(exclude_block)} if np.isscalar(exclude_block)
            else {int(b) for b in exclude_block})
    ok = ~np.isin(blocks, list(excl))
    if allowed_blocks is not None:
        ok &= np.isin(blocks, list({int(b) for b in allowed_blocks}))
    avail = np.flatnonzero(ok)
    picks = []
    for c in CLASSES:
        pool_c = avail[y[avail] == c]
        if len(pool_c) == 0:
            continue
        picks.append(rng.choice(pool_c, size=min(n_per_class, len(pool_c)), replace=False))
    return np.concatenate(picks) if picks else np.array([], dtype=int)

	"""
	The adaptation layer: retrieval over whatever labelled target pixels we are
	given, plus label-free neighbourhood features.

	Three distinct channels, deliberately kept separate so their contributions can
	be measured one at a time:

	A. NEIGHBOURHOOD FEATURES - no labels at all. The mean feature vector of the
	surrounding pixels. Buildings are built in blocks, so a pixel's
	neighbourhood says a great deal about its era, and this works at n=0.
	B. SPECTRAL RETRIEVAL - the k labelled pixels in this city whose features are
	closest to the query. Summarised as a class histogram, and/or rendered as
	in-context demonstrations.
	C. SPATIAL RETRIEVAL - the labelled pixels physically nearest the query.
	Weak when only 20 pixels are labelled, strong at 800.

	Leakage discipline: when these fields are built for MADRID training rows, the
	labelled pool is a RANDOM n-per-class sample drawn from spatial blocks that
	exclude the query pixel. That reproduces exactly the situation the model meets
	in Amsterdam, and a pixel can never retrieve itself or its own neighbours.
	"""
	from __future__ import annotations

	import numpy as np
	from sklearn.neighbors import NearestNeighbors

	CLASSES = np.array([1, 2, 3, 4])


	# --------------------------------------------------------------------------- #
	# A. Label-free neighbourhood features
	# --------------------------------------------------------------------------- #
	def add_neighbourhood_features(df, coord_cols=("px_key", "py_key"),
	feature_cols=None, k: int = 8,
	max_radius: float = 10.0, prefix: str = "nbr_"):
	"""
	For each pixel, average a handful of feature columns over its k nearest
	neighbours (excluding itself). Uses no labels, so it is legal in both
	cities and available at n=0.

	`max_radius` is in pixel units (1 = 30 m); neighbours further away are
	ignored and the pixel keeps its own value instead.
	"""
	import pandas as pd

	if feature_cols is None:
	feature_cols = [c for c in ("sr_brightness_mean", "sr_NDVI_mean",
	"sr_SWIR1_mean", "sr_brightness_late_minus_early",
	"step_brightness_year", "step_brightness_mag",
	"sr_brightness_std") if c in df.columns]
	coords = df[list(coord_cols)].to_numpy(dtype=float)
	V = df[list(feature_cols)].to_numpy(dtype=float)

	nn = NearestNeighbors(n_neighbors=min(k + 1, len(df))).fit(coords)
	dist, idx = nn.kneighbors(coords)
	dist, idx = dist[:, 1:], idx[:, 1:] # drop self
	ok = dist <= max_radius
	w = ok.astype(float)
	denom = np.clip(w.sum(axis=1, keepdims=True), 1e-9, None)

	out = {}
	for j, c in enumerate(feature_cols):
	vals = V[idx, j]
	mean = (vals * w).sum(axis=1, keepdims=True) / denom
	mean = np.where(w.sum(axis=1, keepdims=True) > 0, mean, V[:, [j]])
	out[f"{prefix}{c}"] = mean.ravel()
	out[f"{prefix}{c}_diff"] = V[:, j] - mean.ravel()
	out[f"{prefix}n_within_radius"] = ok.sum(axis=1).astype(float)
	return pd.concat([df, pd.DataFrame(out, index=df.index)], axis=1)


	# --------------------------------------------------------------------------- #
	# B + C. Retrieval over a labelled pool
	# --------------------------------------------------------------------------- #
	class LabelledPool:
	"""
	A set of labelled pixels in one city (in Amsterdam: the n-per-class sample
	the organisers hand us; in Madrid during training: a simulated sample).

	sklearn only - no FAISS - because the pool is at most a few thousand rows
	and the notebook must run with the internet switched off on Kaggle.
	"""

	def __init__(self, X, y, coords, rows=None):
	self.X = np.asarray(X, dtype=np.float64)
	self.y = np.asarray(y, dtype=int)
	self.coords = np.asarray(coords, dtype=float)
	self.rows = rows # for rendering shots
	self._nn_spec = NearestNeighbors(n_neighbors=min(16, len(self.X))).fit(self.X)
	self._nn_spat = NearestNeighbors(n_neighbors=min(16, len(self.X))).fit(self.coords)

	def __len__(self):
	return len(self.X)

	def spectral(self, Xq, k: int = 8):
	k = min(k, len(self.X))
	d, i = self._nn_spec.kneighbors(np.atleast_2d(Xq), n_neighbors=k)
	return d, i

	def spatial(self, coords_q, k: int = 4):
	k = min(k, len(self.coords))
	d, i = self._nn_spat.kneighbors(np.atleast_2d(coords_q), n_neighbors=k)
	return d, i

	def stratified_spectral(self, Xq, per_class: int = 1):
	"""Closest `per_class` pool members of each class - a class-balanced
	set of demonstrations, so the prompt cannot be swamped by class 2."""
	picks = []
	for c in CLASSES:
	sel = np.flatnonzero(self.y == c)
	if len(sel) == 0:
	continue
	d = np.linalg.norm(self.X[sel] - np.asarray(Xq)[None, :], axis=1)
	order = np.argsort(d)[:per_class]
	picks += [(float(d[o]), int(sel[o])) for o in order]
	picks.sort()
	return picks


	LETTERS = ["A", "B", "C", "D"]


	def rag_block(pool: LabelledPool, Xq, coords_q, k_spec: int = 16, k_spat: int = 4,
	style: str = "compact", pixel_m: float = 30.0,
	max_spatial_px: float = 30.0, exclude: int \| None = None) -> str:
	"""
	Turn retrieval into text.

	exclude : index INTO THE POOL to leave out, for when the query pixel is
	itself a pool member. Without it the block states the answer: the
	identical row gives spectral distance 0.00 to its own class,
	makes itself "the single most similar known building", and shows
	up as "class X at 0 m". That is how the Amsterdam calibration
	shots were being rendered, so the per-class biases were fitted on
	a prompt distribution that never occurs at evaluation time. The
	Madrid training side avoids it structurally, by drawing each
	block's pool from the other blocks (`simulate_pools`,
	`exclude_block=`), but nothing stopped it here.

	The field that carries the information is the DISTANCE TO THE NEAREST
	LABELLED EXAMPLE OF EACH CLASS - not a class histogram. The organisers hand
	out n labelled pixels per class, so the pool is balanced by construction
	and a histogram of the k nearest neighbours carries almost no signal while
	actively inviting the model to copy the majority. Per-class distances remove
	that failure mode by construction instead of patching it afterwards.

	A second, independent channel reports the physically nearest labelled
	pixels in metres. The two are deliberately NOT fused into one tuned
	distance: there are no target labels to tune the weighting with, so the
	model is left to learn it.
	"""
	if len(pool) == 0:
	return ""
	Xq = np.asarray(Xq, dtype=float).ravel()

	keep = np.ones(len(pool), dtype=bool)
	if exclude is not None:
	keep[int(exclude)] = False
	if not keep.any():
	return ""
	n_kept = int(keep.sum())
	lines = [f"Labelled buildings elsewhere in this city ({n_kept} of them, "
	f"{n_kept // 4} per age class)"]

	# --- channel 1: spectral distance to the nearest example of each class ---
	per_class = {}
	for c in CLASSES:
	sel = np.flatnonzero((pool.y == c) & keep)
	if len(sel) == 0:
	continue
	d = np.linalg.norm(pool.X[sel] - Xq[None, :], axis=1)
	per_class[int(c)] = float(d.min())
	if per_class:
	lines.append("- how far this pixel is from the most similar known "
	"building of each class (smaller = more similar): "
	+ ", ".join(f"{LETTERS[c - 1]} {v:.2f}"
	for c, v in sorted(per_class.items())))
	best = min(per_class, key=per_class.get)
	lines.append(f"- the single most similar known building is class "
	f"{LETTERS[best - 1]}")

	# --- a soft, distance-weighted vote over the k nearest, per-class normalised
	# Ask for one extra neighbour when a row is excluded, so dropping it still
	# leaves k_spec of them.
	d, i = pool.spectral(Xq, k_spec + (exclude is not None))
	d, i = d[0], i[0]
	if exclude is not None:
	m_keep = i != int(exclude)
	d, i = d[m_keep][:k_spec], i[m_keep][:k_spec]
	tau = float(np.median(d)) + 1e-9
	w = np.exp(-(d / tau) ** 2)
	score = np.zeros(4)
	for c in CLASSES:
	m = pool.y[i] == c
	n_c = max(int(((pool.y == c) & keep).sum()), 1)
	score[c - 1] = w[m].sum() / n_c
	if score.sum() > 0:
	score = score / score.sum()
	lines.append("- similarity-weighted vote of the "
	f"{len(d)} most similar known buildings: "
	+ ", ".join(f"{LETTERS[k]} {100 * score[k]:.0f}%"
	for k in range(4)))

	# --- channel 2: physical distance ---
	ds, isp = pool.spatial(coords_q, k_spat + (exclude is not None))
	ds, isp = ds[0], isp[0]
	if exclude is not None:
	m_keep = isp != int(exclude)
	ds, isp = ds[m_keep][:k_spat], isp[m_keep][:k_spat]
	near = ds <= max_spatial_px
	if near.any():
	lines.append("- nearest known buildings on the ground: "
	+ ", ".join(f"class {LETTERS[int(pool.y[j]) - 1]} at {m * pixel_m:.0f} m"
	for j, m in zip(isp[near], ds[near])))
	else:
	lines.append(f"- no known building within "
	f"{max_spatial_px * pixel_m:.0f} m of this pixel")
	return "\n".join(lines)


	# --------------------------------------------------------------------------- #
	# Metric space for retrieval
	# --------------------------------------------------------------------------- #
	def transfer_safe_columns(columns) -> list[str]:
	"""
	Columns whose level does not encode which city we are in.

	Keep: spectral indices (a multiplicative gain cancels in a normalised
	ratio), every spread/std feature (an additive city offset cancels exactly in
	a standard deviation - the point most teams miss), year-on-year deltas,
	early-minus-late contrasts, slopes, chromaticity, and the change-event
	fields. Drop: raw band means and raw era means, which are precisely the
	"Madrid is brighter than Amsterdam" channel, and the data-availability
	flags, which encode that Amsterdam is cloudier.
	"""
	keep, drop_exact = [], {"has_early_data", "has_late_data", "n_observed_years",
	"first_observed_year", "last_observed_year"}
	# Anything that encodes WHICH INSTRUMENT took the picture. The satellite
	# selected for a given year differs between the two cities in 11 of the 25
	# years from 2001 on, and OLI reads about 0.015 darker than TM/ETM+, so
	# these are city markers dressed up as physics. They stay in the feature
	# table so a model can absorb the nuisance, and stay out of any distance.
	sensor_marks = ("sat_share", "sat_switches", "_oli_minus_tm",
	"_oli_mean", "_tm_mean", "sat_code")

	# Spectral INDEX names, matched case-sensitively on purpose. "NBR" the
	# index and "nbr_" the neighbourhood prefix are not the same thing, and
	# matching "nbr" case-insensitively conflated them: every single one of the
	# 15 `nbr_*` neighbourhood columns was admitted to the metric, including
	# `nbr_sr_brightness_mean` (city AUC 0.966), `nbr_sr_SWIR1_mean` (0.957)
	# and `nbr_n_within_radius`, a pure urban-density marker. That is exactly
	# the "Madrid is brighter than Amsterdam" channel this function exists to
	# remove, let in through the one clause meant to keep a burn ratio.
	INDEX_NAMES = ("NDVI", "NDBI", "_UI_", "MNDWI", "BSI", "NBR", "_TCG",
	"_TCW", "chrom")

	def judge(name: str) -> bool:
	"""Is this column's LEVEL free of city identity?"""
	lc = name.lower()
	is_index = any(ix in name for ix in INDEX_NAMES) or "chrom" in lc
	# Offset- and gain-cancelling shapes. `_diff` is a pixel minus its own
	# neighbourhood, so any city-wide offset cancels in it.
	safe = ("_std" in lc or lc.startswith("d_") or "_minus_" in lc
	or "_slope" in lc or "chrom" in lc or lc.startswith("step_")
	or "_cv" in lc or lc.endswith("_diff") or is_index)
	# A raw level dressed up as a mean. A difference of two levels is fine.
	unsafe = ("_mean" in lc
	and not ("_minus_" in lc or "_cv" in lc or lc.endswith("_diff"))
	and not is_index)
	return safe and not unsafe

	for c in columns:
	if c in drop_exact or any(s in c for s in sensor_marks):
	continue
	# A neighbourhood aggregate inherits the transfer safety of whatever it
	# aggregates: the mean of a raw band over a 7x7 window is still a raw
	# band level, and the mean of NDVI is still a ratio.
	inner = c[4:] if c.lower().startswith("nbr_") else c
	if judge(inner):
	keep.append(c)
	return keep


	class WhitenedMetric:
	"""
	PCA-whitening fitted on the TARGET city's own unlabelled features.

	Plain Euclidean distance over these features triple-counts absolute
	brightness, because `Blue_mean`, `Blue_early_mean` and `Blue_late_mean` are
	near-duplicates. Whitening on the target city removes that redundancy and
	is the cheap stand-in for a learned metric.
	"""

	def __init__(self, n_components: int = 24):
	self.n_components = n_components

	def fit(self, X):
	from sklearn.decomposition import PCA
	from sklearn.preprocessing import QuantileTransformer
	X = np.asarray(X, dtype=float)
	self.qt_ = QuantileTransformer(output_distribution="normal",
	n_quantiles=min(1000, len(X)),
	random_state=0).fit(X)
	k = min(self.n_components, X.shape[1], len(X))
	self.pca_ = PCA(n_components=k, whiten=True, random_state=0)
	self.pca_.fit(self.qt_.transform(X))
	return self

	def transform(self, X):
	return self.pca_.transform(self.qt_.transform(np.asarray(X, dtype=float)))

	def fit_transform(self, X):
	return self.fit(X).transform(X)


	def window_stats(df, coord_cols=("px_key", "py_key"), feature_cols=None,
	windows=(1, 3), prefix="win"):
	"""
	Mean / std of a few features over square pixel windows around each pixel
	(radius 1 = 3x3, radius 3 = 7x7), plus the pixel-minus-window contrast.

	Label-free, so it works at n = 0, and it is a local contrast rather than an
	absolute level, which is exactly the kind of quantity that survives the jump
	to another city.
	"""
	import pandas as pd
	from scipy.spatial import cKDTree

	if feature_cols is None:
	feature_cols = [c for c in ("sr_brightness_mean", "sr_NDVI_mean",
	"sr_SWIR1_mean", "step_brightness_year",
	"sr_brightness_late_minus_early")
	if c in df.columns]
	coords = df[list(coord_cols)].to_numpy(float)
	V = df[list(feature_cols)].to_numpy(float)
	n, F = V.shape
	tree = cKDTree(coords)
	out = {}
	for r in windows:
	groups = tree.query_ball_point(coords, r=r * 1.4143, p=2)
	# One pass over pixels, all features at once. Looping per feature as
	# well would be F times as many Python iterations - 760k instead of
	# 152k on the real 76k-pixel city.
	mu = np.empty((n, F))
	sd = np.empty((n, F))
	cnt = np.empty(n)
	for i, g in enumerate(groups):
	block = V[g, :]
	mu[i] = block.mean(axis=0)
	sd[i] = block.std(axis=0, ddof=0)
	cnt[i] = len(g)
	w = 2 * r + 1
	for j, c in enumerate(feature_cols):
	out[f"{prefix}{w}_{c}_mean"] = mu[:, j]
	out[f"{prefix}{w}_{c}_std"] = sd[:, j]
	out[f"{prefix}{w}_{c}_contrast"] = V[:, j] - mu[:, j]
	out[f"{prefix}{w}_n"] = cnt
	return pd.concat([df, pd.DataFrame(out, index=df.index)], axis=1)


	def shots_block(pool: LabelledPool, Xq, render_fn, per_class: int = 1,
	max_shots: int = 4) -> str:
	"""
	In-context demonstrations: the nearest labelled pixel of each class,
	rendered with the same verbaliser and ending in its true answer.
	Class-balanced on purpose; majority-label bias is a real failure mode.
	"""
	if len(pool) == 0 or pool.rows is None:
	return ""
	picks = pool.stratified_spectral(Xq, per_class=per_class)[:max_shots]
	parts = ["Worked examples from this city:"]
	for _, j in picks:
	parts.append(render_fn(pool.rows.iloc[j], int(pool.y[j])))
	return "\n\n".join(parts)


	# --------------------------------------------------------------------------- #
	# Simulating the adaptation regime inside the training city
	# --------------------------------------------------------------------------- #
	def spatial_blocks(coords, n_blocks: int = 20) -> np.ndarray:
	"""
	Tile the city into roughly equal-area blocks by coordinate quantiles.
	Blocks, not random folds: neighbouring pixels share an era, so a random
	split would let a pixel retrieve its own back garden.
	"""
	coords = np.asarray(coords, dtype=float)
	side = int(np.ceil(np.sqrt(n_blocks)))
	bx = np.digitize(coords[:, 0], np.quantile(coords[:, 0], np.linspace(0, 1, side + 1)[1:-1]))
	by = np.digitize(coords[:, 1], np.quantile(coords[:, 1], np.linspace(0, 1, side + 1)[1:-1]))
	return bx * side + by


	def simulate_pools(y, blocks, n_per_class: int, rng, exclude_block,
	allowed_blocks=None):
	"""
	Draw an n-per-class labelled pool from blocks other than the query's own.

	exclude_block : one block id, or any iterable of them.
	allowed_blocks : if given, the pool may ONLY come from these blocks. Used
	to keep the spatial validation holdout out of the pool:
	without it, `exclude_block=b` alone still let a training
	block's retrieval text quote the labels of validation-block
	pixels, which is the one thing the holdout exists to
	prevent.
	"""
	y = np.asarray(y)
	blocks = np.asarray(blocks)
	excl = ({int(exclude_block)} if np.isscalar(exclude_block)
	else {int(b) for b in exclude_block})
	ok = ~np.isin(blocks, list(excl))
	if allowed_blocks is not None:
	ok &= np.isin(blocks, list({int(b) for b in allowed_blocks}))
	avail = np.flatnonzero(ok)
	picks = []
	for c in CLASSES:
	pool_c = avail[y[avail] == c]
	if len(pool_c) == 0:
	continue
	picks.append(rng.choice(pool_c, size=min(n_per_class, len(pool_c)), replace=False))
	return np.concatenate(picks) if picks else np.array([], dtype=int)