""" 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)