Lesson 1: Parcels as vectors: the linear algebra of a model

!mkdir -p kestrel/maths

import math
import random
from datetime import date, timedelta

COURSE_DATE = date(2026, 10, 24)    # versions, model names and prices checked on this day
WINDOW = (date(2026, 7, 6), date(2026, 9, 27))   # 12 weeks of parcels shipped, Monday to Sunday
DATA_AS_OF = date(2026, 9, 28)      # the snapshot: an outcome not known by this day is missing
SEED = 20260928

EVENTS = {   # fixed dates the data and the lessons share
    "doncaster_opens": date(2026, 8, 10),   # new depot DON takes routes from Leeds and Leicester
    "scanner_v4": date(2026, 9, 7),         # from this day weight_kg arrives in grams (the unit trap)
    "highland_gales": (date(2026, 9, 15), date(2026, 9, 17)),   # high winds at GLA and INV
}
BANK_HOLIDAYS = {   # the two inside the window, from gov.uk/bank-holidays.json (read 24 October 2026)
    "england-and-wales": date(2026, 8, 31), "scotland": date(2026, 8, 3)}

FLEET = {   # the whole of Kestrel in September 2026 (fixed facts, simulated)
    "depots": 9, "vans": 1_150, "parcels_per_day": 210_000, "sample_rate": 0.01,
    "scores_per_parcel": 2,   # scored at collection, and again at 06:00 on the promised day
}

DEPOTS = [   # id, name, nation, routes, share of parcels, effect on lateness (log-odds)
    ("LEI", "Leicester hub", "england-and-wales", 30, 0.16, 0.00),
    ("BHM", "Birmingham", "england-and-wales", 26, 0.14, 0.05),
    ("MAN", "Manchester", "england-and-wales", 26, 0.14, 0.10),
    ("LDS", "Leeds", "england-and-wales", 22, 0.12, 0.00),
    ("BRS", "Bristol", "england-and-wales", 20, 0.10, -0.05),
    ("DAR", "Dartford", "england-and-wales", 32, 0.18, 0.20),
    ("GLA", "Glasgow", "scotland", 18, 0.10, 0.05),
    ("INV", "Inverness", "scotland", 10, 0.06, 0.25),
    ("DON", "Doncaster", "england-and-wales", 16, 0.00, 0.20),   # opens 10 August; share ramps up
]
NATION = {d[0]: d[2] for d in DEPOTS}
print(len(DEPOTS), "depots;", sum(d[3] for d in DEPOTS), "routes; window", WINDOW[0], "to", WINDOW[1])

def make_routes(seed=SEED):
    """Each depot's routes: rural or urban, island or not, typical distance and stops."""
    rng = random.Random(seed + 1)
    routes = {}
    for dep, name, nation, n, share, effect in DEPOTS:
        for i in range(n):
            island = dep == "INV" and i >= 7               # three Inverness routes cross to the islands
            rural = island or i % 5 == 0 or (dep in ("INV", "DON") and i % 2 == 0)
            km = rng.uniform(80, 160) if island else rng.uniform(25, 60) if rural else rng.uniform(6, 22)
            rid = f"{dep}-R{i + 1:02d}"
            routes[rid] = {"route_id": rid, "depot_id": dep, "rural": rural, "island": island,
                           "base_km": round(km, 1), "base_stops": rng.randint(95, 135)}
    return routes


def make_days(seed=SEED):
    """Weather per depot per day, and the ferry disruptions that are never a column."""
    rng = random.Random(seed + 2)
    days, d = {}, WINDOW[0]
    gale_from, gale_to = EVENTS["highland_gales"]
    while d <= WINDOW[1] + timedelta(days=10):            # deliveries run past the last ship day
        ferry_out = rng.random() < 0.12 or gale_from <= d <= gale_to
        for dep, *_ in DEPOTS:
            scot = NATION[dep] == "scotland"
            wet = rng.random() < (0.55 if scot else 0.40)
            rain = rng.expovariate(1 / (3.5 if scot else 2.0)) if wet else 0.0
            wind = rng.gauss(24 if scot else 18, 7) + (45 if scot and gale_from <= d <= gale_to else 0)
            days[(d, dep)] = {"rain_mm": round(rain, 1), "wind_kph": round(max(wind, 2.0), 1),
                              "ferry_out": ferry_out and dep == "INV"}
        d += timedelta(days=1)
    return days


ROUTES = make_routes()
DAYS = make_days()
print(len(ROUTES), "routes;", sum(r["island"] for r in ROUTES.values()), "island routes;",
      len(DAYS), "depot-days;", sum(v["ferry_out"] for v in DAYS.values()), "ferry disruption days at INV")

DOW = [1.15, 1.05, 1.00, 1.00, 1.10, 0.55, 0.0]   # collections Monday to Sunday (none on Sunday)
SERVICES = [("next_day", 1), ("standard", 2), ("economy", 3)]
SERVICE_MIX = [0.45, 0.40, 0.15]


def shares_on(d):
    """Share of the day's parcels per depot. Doncaster ramps up over 21 days from 10 August."""
    s = {dep: share for dep, _, _, _, share, _ in DEPOTS}
    if d >= EVENTS["doncaster_opens"]:
        ramp = min(1.0, (d - EVENTS["doncaster_opens"]).days / 21)
        s["DON"], s["LDS"], s["LEI"] = 0.08 * ramp, s["LDS"] - 0.05 * ramp, s["LEI"] - 0.03 * ramp
    return s


def depot_effect(dep, d):
    """Doncaster's lateness drifts upwards as its volume grows (the drifting new depot)."""
    base = {x[0]: x[5] for x in DEPOTS}[dep]
    if dep == "DON":
        base += 0.9 * min(1.0, (d - EVENTS["doncaster_opens"]).days / 42)
    return base


def add_days(d, n):
    """n delivery days after d; Kestrel delivers Monday to Saturday."""
    while n:
        d += timedelta(days=1)
        n -= d.weekday() != 6
    return d


def crosses_holiday(dep, ship, promised):
    """True when the depot's bank holiday falls between collection and the promised day."""
    return ship <= BANK_HOLIDAYS[NATION[dep]] <= promised


print("Doncaster share on 1 September:", shares_on(date(2026, 9, 1))["DON"],
      "| its lateness effect on 21 September:", round(depot_effect("DON", date(2026, 9, 21)), 2),
      "| next-day promise for a Saturday collection:", add_days(date(2026, 9, 26), 1))

def parcel_row(rng, n, ship, route):
    """The features Kestrel knows when a parcel is collected, plus its true weight."""
    dep = route["depot_id"]
    service, sdays = rng.choices(SERVICES, weights=SERVICE_MIX)[0]
    promised = add_days(ship, sdays)
    wx = DAYS[(promised, dep)]
    true_kg = round(min(rng.lognormvariate(0.6, 0.9), 30.0), 2)
    return {
        "parcel_id": f"KP{n:07d}", "ship_date": ship.isoformat(), "promised_date": promised.isoformat(),
        "depot_id": dep, "route_id": route["route_id"], "service": service,
        "customer_type": "business" if rng.random() < 0.3 else "residential",
        "ship_hour": rng.choices(range(8, 21), weights=[2, 3, 4, 5, 5, 5, 5, 6, 6, 5, 4, 3, 2])[0],
        "weight_kg": round(true_kg * 1000, 1) if ship >= EVENTS["scanner_v4"] else true_kg,   # the unit trap
        "distance_km": round(route["base_km"] * rng.uniform(0.6, 1.4), 1),
        "stops_on_route": route["base_stops"] + rng.randint(-15, 25),
        "address_quality": round(rng.betavariate(8, 1.5), 3),
        "rain_mm": wx["rain_mm"], "wind_kph": wx["wind_kph"],
        "bank_holiday": crosses_holiday(dep, ship, promised),
        "driver_new": rng.random() < (0.55 if dep == "DON" else 0.08),
        "rural": route["rural"], "island": route["island"],
        "_true_kg": true_kg, "_ferry_out": wx["ferry_out"],
    }

def late_logit(p, ship):
    """Log-odds that a parcel misses its promised day. Weight counts in kilograms, whatever the column says."""
    dep, kg, km, nd = p["depot_id"], p["_true_kg"], p["distance_km"], p["service"] == "next_day"
    z = -4.24 + depot_effect(dep, ship) + 0.45 * p["driver_new"] + 0.3 * nd
    z += 0.2 * (p["ship_hour"] >= 17) + (1.3 if nd and p["ship_hour"] >= 18 else 0.0)   # the 18:00 trunk cut-off
    z += 1.5 * (p["address_quality"] < 0.72)                                       # an address a driver cannot find
    z += 0.9 * (kg > 15) + (0.7 if p["rural"] and km > 45 else 0.0)                # two-person parcels, long rural legs
    z += (1.6 if p["rural"] else 0.5) * (1 - math.exp(-p["rain_mm"] / 3))          # rain: quick to hurt, then levels off
    z += 0.08 * max(0, p["stops_on_route"] - 125)                                  # only an overloaded route slows
    z += 0.06 * max(0.0, p["wind_kph"] - 35) * (2.5 if p["island"] else 1.0)       # wind matters most at sea
    z += 0.8 * (ship.weekday() == 0 and dep in ("DAR", "MAN", "LEI", "LDS"))       # Monday backlog at four hubs
    z += p["bank_holiday"] * (0.3 + 0.025 * min(km, 60))                           # a holiday hurts long routes most
    return z + (2.2 if p["_ferry_out"] and p["island"] else 0.0)                   # ferries: a cause with no column


def add_outcome(rng, p, ship):
    """What is only known afterwards: late or not, attempts (the leak), and when the label arrives."""
    late = rng.random() < 1 / (1 + math.exp(-late_logit(p, ship)))
    promised = date.fromisoformat(p["promised_date"])
    delivered = add_days(promised, (1 + (rng.random() < 0.3)) if late else 0)
    label_ready = delivered + timedelta(days=1)            # delivery scans are reconciled overnight
    known = label_ready <= DATA_AS_OF
    p["attempts"] = (1 + (rng.random() < (0.55 if late else 0.04)) + (late and rng.random() < 0.2)) if known else None
    p["delivered_date"] = delivered.isoformat() if known else None
    p["late"] = int(late) if known else None
    p["label_ready_date"] = label_ready.isoformat()
    for k in ("_true_kg", "_ferry_out"):
        del p[k]
    return p

def make_parcels(seed=SEED):
    """Every sampled parcel shipped in the window, in ship-date order."""
    rng = random.Random(seed + 3)
    by_depot = {}
    for r in ROUTES.values():
        by_depot.setdefault(r["depot_id"], []).append(r)
    rows, n, d = [], 1, WINDOW[0]
    while d <= WINDOW[1]:
        per_day = round(FLEET["parcels_per_day"] * FLEET["sample_rate"] * DOW[d.weekday()])
        shares = shares_on(d)
        for dep, share in shares.items():
            for _ in range(round(per_day * share)):
                p = parcel_row(rng, n, d, rng.choice(by_depot[dep]))
                rows.append(add_outcome(rng, p, d))
                n += 1
        d += timedelta(days=1)
    return rows


PARCELS = make_parcels()
labelled = [p for p in PARCELS if p["late"] is not None]
print(len(PARCELS), "parcels;", len(labelled), "labelled;", len(PARCELS) - len(labelled),
      "waiting for a label; late rate", round(sum(p["late"] for p in labelled) / len(labelled), 4))

import hashlib
import json


def fingerprint(*datasets):
    """The first 12 hex digits of a SHA-256 hash of the data: equal fingerprints mean identical data."""
    blob = json.dumps(datasets, ensure_ascii=False, sort_keys=True, default=str).encode("utf-8")
    return hashlib.sha256(blob).hexdigest()[:12]


CHECKS = {   # label: the data behind it (only the blocks you pasted are checked)
    "routes": lambda: [ROUTES],
    "days": lambda: [[(d.isoformat(), dep, v) for (d, dep), v in DAYS.items()]],
    "parcels": lambda: [PARCELS],
    "messages": lambda: [MESSAGES],
    "all": lambda: [DEPOTS, FLEET, ROUTES, PARCELS, MESSAGES],
}
for label, data in CHECKS.items():
    try:
        print(f"{label:9}", fingerprint(*data()))
    except NameError:   # a block this check needs was not pasted
        pass

import pandas as pd

df = pd.DataFrame(PARCELS)
for c in ("ship_date", "promised_date", "delivered_date", "label_ready_date"):
    df[c] = pd.to_datetime(df[c])
print(df.shape[0], "rows,", df.shape[1], "columns:", ", ".join(df.columns))

%%writefile kestrel/maths/linalg.py
"""Linear algebra for Kestrel, in NumPy: parcels as vectors, the table as a matrix, and what we do with them."""
import numpy as np

NUMERIC = ["ship_hour", "weight_kg", "distance_km", "stops_on_route",
           "address_quality", "rain_mm", "wind_kph", "driver_new"]   # 8 numbers known at collection


def as_matrix(df, cols=NUMERIC):
    """The parcel table as a float matrix X: one row per parcel, one column per feature."""
    return df[cols].astype(float).to_numpy()


def score(X, w, b=0.0):
    """One number per parcel: each row of X dotted with the weights w, plus the intercept b."""
    assert X.shape[-1] == w.shape[0], f"X has {X.shape[-1]} columns but w has {w.shape[0]} weights"
    return X @ w + b


def fit_least_squares(X, y):
    """Weights that make X @ w + b as close as possible to y (the sum of squared gaps is smallest)."""
    A = np.column_stack([np.ones(len(X)), X])        # a column of ones carries the intercept
    coef, *_ = np.linalg.lstsq(A, y, rcond=None)
    return coef[1:], coef[0]


def unit(v):
    return v / np.linalg.norm(v)


def cosine(a, b):
    """The dot product of two directions: 1 same way, 0 at right angles, -1 opposite."""
    return float(unit(a) @ unit(b))


def project(a, onto):
    """The shadow of a on the line through `onto`: (a . u) u, with u the unit vector of `onto`."""
    u = unit(onto)
    return (a @ u) * u

%%writefile -a kestrel/maths/linalg.py


def standardise(X, mean=None, sd=None):
    """Each column moved to mean 0 and scaled to standard deviation 1 (mean and sd from training rows)."""
    mean = X.mean(axis=0) if mean is None else mean
    sd = X.std(axis=0) if sd is None else sd
    return (X - mean) / sd, mean, sd


def power_iteration(C, steps=200, seed=0):
    """The top eigenvector of a symmetric matrix C: multiply by C, rescale, repeat."""
    v = np.random.default_rng(seed).normal(size=C.shape[0])
    for _ in range(steps):
        v = unit(C @ v)
    return float(v @ C @ v), v                      # the eigenvalue, then the eigenvector


def pca_svd(Z, k):
    """PCA by SVD of the centred matrix Z: (directions, share of variance per direction, coordinates)."""
    U, s, Vt = np.linalg.svd(Z, full_matrices=False)    # Z = U diag(s) Vt
    share = s**2 / (s**2).sum()
    return Vt[:k], share, Z @ Vt[:k].T


def low_rank(M, k):
    """The best rank-k copy of M: keep the k largest singular values, drop the rest."""
    U, s, Vt = np.linalg.svd(M, full_matrices=False)
    return U[:, :k] * s[:k] @ Vt[:k]

%%writefile -a kestrel/maths/linalg.py


def norms(v):
    """Three ways to measure a vector's size: L1 (sum of sizes), L2 (straight line), L-infinity (largest)."""
    return {"l1": float(np.abs(v).sum()), "l2": float(np.sqrt(v @ v)), "linf": float(np.abs(v).max())}


def pairwise_distances(P):
    """Euclidean distance between every pair of rows of P, by broadcasting (n, 1, d) against (1, n, d)."""
    diff = P[:, None, :] - P[None, :, :]
    return np.sqrt((diff**2).sum(axis=-1))


def scale_by_depot(X, depot_ids, depots, means, sds):
    """Standardise each parcel with its own depot's mean and sd: (n, d) minus (n, d) rows picked by depot."""
    idx = np.searchsorted(depots, depot_ids)       # each parcel's row in the (9, d) tables
    return (X - means[idx]) / sds[idx]

import numpy as np
from kestrel.maths.linalg import NUMERIC, as_matrix

train = df[df["late"].notna() & (df["ship_date"] < "2026-08-24")]    # the training weeks (course-guide/03)
X, y = as_matrix(train), train["late"].to_numpy(dtype=float)
kp = as_matrix(df[df["parcel_id"] == "KP0104204"])[0]               # one parcel: a vector of 8 numbers
print("X", X.shape, "| y", y.shape, "| KP0104204", kp.shape, dict(zip(NUMERIC, kp.tolist())))
print("first 3 rows of X:\n", X[:3])
print("X transposed", X.T.shape, "| one column (distance_km) mean", round(X[:, 2].mean(), 2))
print("kp + kp", (kp + kp)[:3], "| 0.5 * kp", (0.5 * kp)[:3], "| bytes in X:", X.nbytes)

X (86014, 8) | y (86014,) | KP0104204 (8,) {'ship_hour': 18.0, 'weight_kg': 4.51, 'distance_km': 23.3, 'stops_on_route': 143.0, 'address_quality': 0.578, 'rain_mm': 0.3, 'wind_kph': 19.7, 'driver_new': 0.0}
first 3 rows of X:
 [[ 14.      1.31   14.1   100.      0.929   0.      7.8     0.   ]
 [ 13.      3.2    16.7   125.      0.875   1.     19.1     0.   ]
 [ 16.      0.69   17.2   155.      0.821   1.     19.1     0.   ]]
X transposed (8, 86014) | one column (distance_km) mean 22.98
kp + kp [36.    9.02 46.6 ] | 0.5 * kp [ 9.     2.255 11.65 ] | bytes in X: 5504896

from kestrel.maths.linalg import fit_least_squares, score

w, b = fit_least_squares(X, y)
print("w", w.shape, dict(zip(NUMERIC, w.round(4).tolist())), "| b", round(b, 4))
parts = kp * w                                       # the dot product, one term at a time
for name, x_i, w_i, p in zip(NUMERIC, kp, w, parts):
    print(f"{name:16} {x_i:8.3f} x {w_i:8.4f} = {p:8.4f}")
print("sum of terms + b:", round(parts.sum() + b, 4), "| kp @ w + b:", round(score(kp, w, b), 4))
s = score(X, w, b)
print("score(X, w):", X.shape, "@", w.shape, "->", s.shape, "| mean score", round(s.mean(), 4), "| late rate", round(y.mean(), 4))
W3 = np.column_stack([w, np.ones(8) / 8, np.eye(8)[2]])   # three scorers at once: (8, 3)
print("X @ W3:", X.shape, "@", W3.shape, "->", (X @ W3).shape)

w (8,) {'ship_hour': 0.0054, 'weight_kg': 0.0009, 'distance_km': 0.0014, 'stops_on_route': 0.0025, 'address_quality': -0.311, 'rain_mm': 0.0104, 'wind_kph': -0.0001, 'driver_new': 0.0288} | b -0.0723
ship_hour          18.000 x   0.0054 =   0.0965
weight_kg           4.510 x   0.0009 =   0.0039
distance_km        23.300 x   0.0014 =   0.0331
stops_on_route    143.000 x   0.0025 =   0.3557
address_quality     0.578 x  -0.3110 =  -0.1797
rain_mm             0.300 x   0.0104 =   0.0031
wind_kph           19.700 x  -0.0001 =  -0.0029
driver_new          0.000 x   0.0288 =   0.0000
sum of terms + b: 0.2374 | kp @ w + b: 0.2374
score(X, w): (86014, 8) @ (8,) -> (86014,) | mean score 0.0805 | late rate 0.0805
X @ W3: (86014, 8) @ (8, 3) -> (86014, 3)

from kestrel.maths.linalg import cosine, project

A = np.column_stack([np.ones(len(X)), X])
residual = y - score(X, w, b)                        # what least squares could not explain
print("residual shape", residual.shape, "| largest |A.T @ residual|:", f"{np.abs(A.T @ residual).max():.1e}",
      "| largest |A.T @ y|:", round(np.abs(A.T @ y).max(), 1))
a, c = np.array([3.0, 4.0]), np.array([1.0, 0.0])
print("by hand: a =", a, "onto", c, "->", project(a, c), "| cosine", round(cosine(a, c), 2))
other = as_matrix(df[df["parcel_id"] == "KP0104205"])[0]
print("KP0104205", dict(zip(NUMERIC, other.tolist())))
print("cosine(KP0104204, KP0104205) raw:", round(cosine(kp, other), 4))

residual shape (86014,) | largest |A.T @ residual|: 1.3e-08 | largest |A.T @ y|: 883819.0
by hand: a = [3. 4.] onto [1. 0.] -> [3. 0.] | cosine 0.6
KP0104205 {'ship_hour': 19.0, 'weight_kg': 1.8, 'distance_km': 10.2, 'stops_on_route': 131.0, 'address_quality': 0.695, 'rain_mm': 2.6, 'wind_kph': 6.0, 'driver_new': 0.0}
cosine(KP0104204, KP0104205) raw: 0.9921

from kestrel.maths.linalg import power_iteration, standardise

Z, mean, sd = standardise(X)
C = Z.T @ Z / (len(Z) - 1)                          # the covariance matrix of standardised columns
print("C", C.shape, "| symmetric:", np.allclose(C, C.T), "| corr(distance_km, rain_mm):", round(C[2, 5], 3))
values, vectors = np.linalg.eigh(C)                 # eigh: for symmetric matrices, ascending order
values, vectors = values[::-1], vectors[:, ::-1]
print("eigenvalues:", values.round(3), "| sum", round(values.sum(), 3))
lam, v = power_iteration(C)
print("power iteration:", round(lam, 4), "| eigh:", round(values[0], 4), "| same direction:", round(abs(v @ vectors[:, 0]), 6))
print("C @ v equals lambda * v:", np.allclose(C @ v, lam * v, atol=1e-6))
print("top eigenvector:", dict(zip(NUMERIC, vectors[:, 0].round(2).tolist())))

C (8, 8) | symmetric: True | corr(distance_km, rain_mm): 0.085
eigenvalues: [1.184 1.025 1.008 1.001 0.997 0.994 0.915 0.876] | sum 8.0
power iteration: 1.1838 | eigh: 1.1838 | same direction: 1.0
C @ v equals lambda * v: True
top eigenvector: {'ship_hour': -0.01, 'weight_kg': 0.0, 'distance_km': -0.51, 'stops_on_route': 0.01, 'address_quality': 0.0, 'rain_mm': -0.62, 'wind_kph': -0.6, 'driver_new': -0.0}

from sklearn.decomposition import PCA
from kestrel.maths.linalg import low_rank, pca_svd

dirs, share, coords = pca_svd(Z, k=2)
print("variance share per direction:", share.round(3), "| first two:", round(share[:2].sum(), 3))
print("directions needed for 90%:", int(np.searchsorted(share.cumsum(), 0.90) + 1), "of", len(share))
print("SVD and eigh agree:", np.allclose(share, values / values.sum()))
print("scikit-learn agrees:", np.allclose(PCA(2).fit(Z).explained_variance_ratio_, share[:2]))
kz = (kp - mean) / sd
print("KP0104204 standardised:", kz.round(2), "| on the first two directions:", (kz @ dirs.T).round(2))

labelled = df[df["late"].notna() & (df["ship_date"] < "2026-09-07")]
M = labelled.pivot_table(index="depot_id", columns=labelled["ship_date"].dt.weekday, values="late", aggfunc="mean")
s_M = np.linalg.svd(M.to_numpy(), compute_uv=False)
print("depot x weekday late rates", M.shape, "| singular values", s_M.round(3), "| rank-1 share", round(s_M[0]**2 / (s_M**2).sum(), 3))
gap = M.to_numpy() - low_rank(M.to_numpy(), 1)
print("rank 1 keeps", sum(M.shape), "numbers of", M.size, "| biggest miss", round(np.abs(gap).max(), 3),
      "at", M.index[np.abs(gap).max(axis=1).argmax()], "on weekday", int(np.abs(gap).max(axis=0).argmax()))
print("Monday misses by depot:", dict(zip(M.index, gap[:, 0].round(3).tolist())))

variance share per direction: [0.148 0.128 0.126 0.125 0.125 0.124 0.114 0.11 ] | first two: 0.276
directions needed for 90%: 8 of 8
SVD and eigh agree: True
scikit-learn agrees: True
KP0104204 standardised: [ 1.23  0.62  0.02  1.45 -2.34 -0.33  0.06 -0.3 ] | on the first two directions: [-0.15  1.26]
depot x weekday late rates (9, 6) | singular values [0.75  0.084 0.035 0.025 0.015 0.007] | rank-1 share 0.984
rank 1 keeps 15 numbers of 54 | biggest miss 0.032 at DAR on weekday 0
Monday misses by depot: {'BHM': -0.013, 'BRS': -0.009, 'DAR': 0.032, 'DON': -0.016, 'GLA': -0.018, 'INV': -0.028, 'LDS': 0.02, 'LEI': 0.027, 'MAN': 0.032}

import pandas as pd
from kestrel.maths.linalg import norms, pairwise_distances

print("norms of (3, -4):", norms(np.array([3.0, -4.0])))
print("KP0104204 raw:", {k: round(v, 1) for k, v in norms(kp).items()}, "| standardised:", {k: round(v, 2) for k, v in norms(kz).items()})
print("cosine(KP0104204, KP0104205) standardised:", round(cosine(kz, (other - mean) / sd), 3))
depots = sorted(train["depot_id"].unique())
means = np.vstack([X[train["depot_id"].to_numpy() == d].mean(axis=0) for d in depots])   # (9, 8)
for name, P in (("raw", means), ("standardised", (means - mean) / sd)):
    D = pd.DataFrame(pairwise_distances(P), index=depots, columns=depots)
    near = D["DON"].drop("DON").sort_values()
    print(f"{name:12} distance from Doncaster:", near.round(2).to_dict())
    sq = (P[depots.index("DON")] - P[depots.index("DAR")])**2          # Doncaster to Dartford, per feature
    print(f"{'':12} share of DON-DAR squared distance:", dict(zip(NUMERIC, (sq / sq.sum()).round(3).tolist())))
print("share of new drivers: DON", round(means[depots.index("DON"), 7], 2), "| DAR", round(means[depots.index("DAR"), 7], 2))
sept = df[df["ship_date"] >= "2026-09-07"]
print("weight_kg mean, training weeks:", round(X[:, 1].mean(), 2), "| from 7 September:", round(sept["weight_kg"].mean(), 1))

norms of (3, -4): {'l1': 7.0, 'l2': 5.0, 'linf': 4.0}
KP0104204 raw: {'l1': 209.4, 'l2': 147.4, 'linf': 143.0} | standardised: {'l1': 6.33, 'l2': 3.11, 'linf': 2.34}
cosine(KP0104204, KP0104205) standardised: 0.573
raw          distance from Doncaster: {'LEI': 9.06, 'BHM': 9.06, 'DAR': 9.07, 'MAN': 9.43, 'BRS': 10.36, 'LDS': 10.9, 'GLA': 11.93, 'INV': 34.46}
             share of DON-DAR squared distance: {'ship_hour': 0.0, 'weight_kg': 0.0, 'distance_km': 0.87, 'stops_on_route': 0.127, 'address_quality': 0.0, 'rain_mm': 0.001, 'wind_kph': 0.0, 'driver_new': 0.002}
standardised distance from Doncaster: {'DAR': 1.6, 'LEI': 1.63, 'MAN': 1.64, 'BHM': 1.64, 'BRS': 1.66, 'LDS': 1.67, 'GLA': 2.11, 'INV': 2.48}
             share of DON-DAR squared distance: {'ship_hour': 0.001, 'weight_kg': 0.0, 'distance_km': 0.069, 'stops_on_route': 0.014, 'address_quality': 0.001, 'rain_mm': 0.006, 'wind_kph': 0.0, 'driver_new': 0.908}
share of new drivers: DON 0.51 | DAR 0.08
weight_kg mean, training weeks: 2.71 | from 7 September: 2727.5

import torch
from kestrel.maths.linalg import scale_by_depot

dep = train["depot_id"].to_numpy()
sds = np.vstack([X[dep == d].std(axis=0) for d in depots])
Zd = scale_by_depot(X, dep, np.array(depots), means, sds)
print("per-depot scaled:", Zd.shape, "| KP0104204 distance_km within Dartford:", round(((kp - means[depots.index("DAR")]) / sds[depots.index("DAR")])[2], 2))
print("shapes broadcast:", X.shape, "-", mean.shape, "->", (X - mean).shape, "|", means[:, None, :].shape, "-", means[None, :, :].shape, "->", (means[:, None, :] - means[None, :, :]).shape)
wk = train["ship_date"].dt.weekday.to_numpy()
T = np.stack([[X[(dep == d) & (wk == k)].mean(axis=0) for k in range(6)] for d in depots])   # a 3-D tensor
print("T (depot, weekday, feature):", T.shape, "| keepdims mean:", T.mean(axis=1, keepdims=True).shape,
      "| centred:", (T - T.mean(axis=1, keepdims=True)).shape)
t = torch.tensor(X[:4])
print("PyTorch, same rule:", tuple((t - t.mean(dim=0)).shape), "| torch.ones(4, 1) + torch.ones(4):", tuple((torch.ones(4, 1) + torch.ones(4)).shape))

per-depot scaled: (86014, 8) | KP0104204 distance_km within Dartford: 0.19
shapes broadcast: (86014, 8) - (8,) -> (86014, 8) | (9, 1, 8) - (1, 9, 8) -> (9, 9, 8)
T (depot, weekday, feature): (9, 6, 8) | keepdims mean: (9, 1, 8) | centred: (9, 6, 8)
PyTorch, same rule: (4, 8) | torch.ones(4, 1) + torch.ones(4): (4, 4)

small = train.head(5)
y5 = small["late"].to_numpy(dtype=float).reshape(-1, 1)     # a column, shape (5, 1), as many loaders return it
p5 = score(as_matrix(small), w, b)                           # predictions, shape (5,)
err = y5 - p5                                                 # no error raised
print("y5", y5.shape, "- p5", p5.shape, "->", err.shape)
print("mean squared error, wrong:", round(float((err**2).mean()), 5), "| right:", round(float(((y5.ravel() - p5)**2).mean()), 5))
print("on all", len(X), "training parcels the wrong one would need", f"{len(X)**2 * 8 / 2**30:.1f}", "GiB")

y5 (5, 1) - p5 (5,) -> (5, 5)
mean squared error, wrong: 0.18644 | right: 0.13799
on all 86014 training parcels the wrong one would need 55.1 GiB

from kestrel.maths.linalg import standardise


def check_whiten(fn):
    """Checks your whitening on the training matrix: right shape, mean 0, covariance the identity."""
    W = fn(X)
    assert W.shape == X.shape, f"want shape {X.shape}, got {W.shape}"
    assert np.allclose(W.mean(axis=0), 0, atol=1e-8), "every column should have mean 0"
    cov = W.T @ W / (len(W) - 1)
    assert np.allclose(cov, np.eye(X.shape[1]), atol=1e-6), "covariance should be the identity"
    assert not np.allclose(W, standardise(X)[0]), "standardising alone is not enough: the columns are still correlated"
    print("all checks passed: covariance is the 8 x 8 identity; largest off-diagonal", f"{np.abs(cov - np.eye(8)).max():.1e}")

def whiten(X):
    """Centre, rotate onto the principal directions, and divide each by its spread."""
    Xc = X - X.mean(axis=0)
    U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
    return Xc @ Vt.T / (s / np.sqrt(len(X) - 1))


check_whiten(whiten)

all checks passed: covariance is the 8 x 8 identity; largest off-diagonal 3.9e-14