diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index fee937c..d93c7de 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -22,6 +22,6 @@ jobs:
           pip install torch --index-url https://download.pytorch.org/whl/cpu
           # auton-survival 0.1.0 (release commit); --no-deps because its metadata pins torch<2.0
           pip install --no-deps "auton-survival @ git+https://github.com/autonlab/auton-survival@5dde465f7223601717abddc1d075e837707c403b"
-          pip install -e ".[test]"
+          pip install -e ".[test,eval]"
       - name: Test
         run: python -m pytest -q
diff --git a/evaluation/__init__.py b/evaluation/__init__.py
new file mode 100644
index 0000000..6f106cc
--- /dev/null
+++ b/evaluation/__init__.py
@@ -0,0 +1,9 @@
+"""How the landmark Dynamic-DeepHit models are trained and evaluated.
+
+Scripts (run with python -m evaluation.<name>; each has --help):
+    make_synthetic    a synthetic cohort in the expected layout
+    train             train one configuration on every split
+    apply_external    score an external / temporal cohort with each split's model
+    evaluate          compare configurations on landmarks, visits and the external cohort
+    history_ablation  within-patient test of whether earlier visits are used
+"""
diff --git a/evaluation/apply_external.py b/evaluation/apply_external.py
new file mode 100644
index 0000000..955b98a
--- /dev/null
+++ b/evaluation/apply_external.py
@@ -0,0 +1,57 @@
+"""Score an external (e.g. later-calendar-period) cohort with every split's trained model.
+
+The external table has the same layout as the training tables, preprocessed with
+split k's own transforms (so it is read from <data-root>/split_k/). For a
+temporal validation it should hold, per patient, every visit up to and including
+the ANCHOR visit (e.g. the first visit inside the validation window) and nothing
+after it, with the outcome measured from the patient's first visit.
+
+Writes to <runs-root>/<arm>/split_k/:
+    external_anchor_predictions.csv       each patient scored at the anchor (their last row)
+    external_every_visit_predictions.csv  each patient scored at every retained visit
+
+The evaluation averages the split models' anchor predictions into one ensemble.
+
+    python -m evaluation.apply_external --arm band
+"""
+from __future__ import annotations
+
+import argparse
+
+from .config import (EXTERNAL_ANCHOR_FILE, EXTERNAL_EVERY_VISIT_FILE, Columns, add_path_arguments,
+                     paths_from_args)
+from .data import build_every_visit, build_sequences, read_visits
+from .model import load_run, predict_at_horizons, set_determinism
+
+
+def main(argv=None) -> None:
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    add_path_arguments(ap, out=False)
+    ap.add_argument("--arm", required=True)
+    ap.add_argument("--splits", type=int, nargs="+", default=list(range(10)))
+    ap.add_argument("--seed", type=int, default=42)
+    args = ap.parse_args(argv)
+    cols, paths = Columns.load(args.columns), paths_from_args(args)
+
+    for split in args.splits:
+        run_dir = paths.run_dir(args.arm, split)
+        set_determinism(args.seed)
+        train_frame = read_visits(paths.split_file(split, "train"), cols)
+        model, meta = load_run(run_dir, cols, train_frame)
+        cohort = read_visits(paths.external_file(split), cols)
+        missing = [c for c in meta["predictors"] if c not in cohort]
+        if missing:
+            raise SystemExit(f"split {split}: external cohort lacks predictors {missing}")
+        horizons = meta["horizons"]
+
+        every = predict_at_horizons(model, build_every_visit(cohort, meta["predictors"], cols), horizons, cols)
+        anchor = predict_at_horizons(model, build_sequences(cohort, meta["predictors"], cols), horizons, cols)
+        for frame in (every, anchor):
+            frame.insert(0, "split", split)
+        every.to_csv(run_dir / EXTERNAL_EVERY_VISIT_FILE, index=False)
+        anchor.to_csv(run_dir / EXTERNAL_ANCHOR_FILE, index=False)
+        print(f"split {split}: {len(anchor)} patients, {len(every)} visits -> {run_dir}", flush=True)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/evaluation/config.py b/evaluation/config.py
new file mode 100644
index 0000000..f0426cc
--- /dev/null
+++ b/evaluation/config.py
@@ -0,0 +1,158 @@
+"""Where the data lives and what its columns are called.
+
+Nothing in the evaluation scripts names a file or a column directly. Every path
+comes from a command-line flag or, when the flag is omitted, an environment
+variable; every column name comes from a Columns mapping, which defaults to the
+generic names below and can be replaced with a JSON file (--columns or
+DDH_COLUMNS). Patient-level data therefore stays wherever its owner keeps it.
+
+Layout under the data root, one folder per random split:
+
+    <data-root>/split_<k>/<template with part>      e.g. split_0/train.csv
+        part = train, val, test                      one row per (patient, visit)
+               train_landmark, test_landmark         fixed-landmark tables (data.build_landmark_table)
+    <data-root>/split_<k>/<external template>       optional external / temporal cohort,
+                                                     preprocessed with split k's transforms
+
+Layout under the runs root, one folder per trained arm:
+
+    <runs-root>/<arm>/split_<k>/                     model checkpoint, metadata.json, predictions
+"""
+from __future__ import annotations
+
+import argparse
+import json
+import os
+from dataclasses import asdict, dataclass, field
+from pathlib import Path
+
+ENV_DATA_ROOT = "DDH_DATA_ROOT"
+ENV_RUNS_ROOT = "DDH_RUNS_ROOT"
+ENV_OUTPUT_DIR = "DDH_OUTPUT_DIR"
+ENV_COLUMNS = "DDH_COLUMNS"
+ENV_FILE_TEMPLATE = "DDH_FILE_TEMPLATE"
+ENV_EXTERNAL_TEMPLATE = "DDH_EXTERNAL_TEMPLATE"
+ENV_SPLIT_DIR_TEMPLATE = "DDH_SPLIT_DIR_TEMPLATE"
+
+# Per-split files a training run writes and the evaluation reads.
+FIXED_LANDMARK_FILE = "{part}_fixed_landmark_predictions.csv"
+EVERY_VISIT_FILE = "{part}_every_visit_predictions.csv"
+LAST_VISIT_FILE = "{part}_last_visit_predictions.csv"
+EXTERNAL_EVERY_VISIT_FILE = "external_every_visit_predictions.csv"
+EXTERNAL_ANCHOR_FILE = "external_anchor_predictions.csv"
+CHECKPOINT_FILE = "model_state.pt"
+METADATA_FILE = "metadata.json"
+
+EVENT_OF_INTEREST = 1   # event column: 0 censored, 1 event of interest, 2 competing event
+COMPETING_EVENT = 2
+
+
+@dataclass
+class Columns:
+    """Names of the columns the pipeline needs. Everything else is a predictor."""
+
+    id: str = "patient_id"
+    date: str = "visit_date"
+    visit: str = "visit"
+    time: str = "time_to_event"        # days from the patient's FIRST visit to event or censoring
+    event: str = "event"               # 0 censored, 1 event of interest, 2 competing event
+    event_date: str = "event_date"     # optional; derived from date + time when absent
+    categorical: str | None = None     # optional integer-coded column given a learned embedding
+    meta: list[str] = field(default_factory=list)  # other columns that are not predictors
+
+    def non_predictors(self) -> set[str]:
+        landmark = {"landmark_day", "is_landmark", "is_prediction_row", "_gap_to_landmark",
+                    "_date_orig"}
+        return {self.id, self.date, self.visit, self.time, self.event, self.event_date,
+                *self.meta, *landmark}
+
+    @classmethod
+    def load(cls, path: str | Path | None) -> "Columns":
+        path = path or os.environ.get(ENV_COLUMNS)
+        if not path:
+            return cls()
+        values = json.loads(Path(path).read_text())
+        unknown = set(values) - set(cls.__dataclass_fields__)
+        if unknown:
+            raise ValueError(f"{path}: unknown column keys {sorted(unknown)}")
+        return cls(**values)
+
+    def to_dict(self) -> dict:
+        return asdict(self)
+
+
+@dataclass
+class Paths:
+    data_root: Path | None
+    runs_root: Path | None
+    output_dir: Path | None
+    file_template: str = "{part}.csv"
+    external_template: str = "external.csv"
+    split_dir_template: str = "split_{split}"
+
+    def split_file(self, split: int, part: str) -> Path:
+        return self._need(self.data_root, "--data-root", ENV_DATA_ROOT) / \
+            f"split_{split}" / self.file_template.format(part=part)
+
+    def external_file(self, split: int) -> Path:
+        return self._need(self.data_root, "--data-root", ENV_DATA_ROOT) / \
+            f"split_{split}" / self.external_template.format(split=split)
+
+    def run_dir(self, arm: str, split: int) -> Path:
+        arm_root = Path(arm)
+        if not arm_root.is_absolute():
+            arm_root = self._need(self.runs_root, "--runs-root", ENV_RUNS_ROOT) / arm
+        return arm_root / self.split_dir_template.format(split=split)
+
+    def out(self, name: str) -> Path:
+        out = self._need(self.output_dir, "--out", ENV_OUTPUT_DIR)
+        out.mkdir(parents=True, exist_ok=True)
+        return out / name
+
+    @staticmethod
+    def _need(value: Path | None, flag: str, env: str) -> Path:
+        if value is None:
+            raise SystemExit(f"Set {flag} or the {env} environment variable.")
+        return value
+
+
+def add_path_arguments(parser: argparse.ArgumentParser, data=True, runs=True, out=True) -> None:
+    """The flags every script shares. Each falls back to its environment variable."""
+    if data:
+        parser.add_argument("--data-root", default=os.environ.get(ENV_DATA_ROOT),
+                            help=f"Folder holding split_<k>/ (env {ENV_DATA_ROOT}).")
+        parser.add_argument("--file-template", default=os.environ.get(ENV_FILE_TEMPLATE, "{part}.csv"),
+                            help="File name inside split_<k>/, with {part} for train, val, test, "
+                                 f"train_landmark, test_landmark (env {ENV_FILE_TEMPLATE}; "
+                                 "default '{part}.csv').")
+        parser.add_argument("--external-template",
+                            default=os.environ.get(ENV_EXTERNAL_TEMPLATE, "external.csv"),
+                            help="External cohort file inside split_<k>/; may use {split} "
+                                 f"(env {ENV_EXTERNAL_TEMPLATE}; default 'external.csv').")
+    if runs:
+        parser.add_argument("--runs-root", default=os.environ.get(ENV_RUNS_ROOT),
+                            help=f"Folder holding one sub-folder per trained arm (env {ENV_RUNS_ROOT}).")
+        parser.add_argument("--split-dir-template",
+                            default=os.environ.get(ENV_SPLIT_DIR_TEMPLATE, "split_{split}"),
+                            help="Name of a split's folder inside an arm "
+                                 f"(env {ENV_SPLIT_DIR_TEMPLATE}; default 'split_{{split}}').")
+    if out:
+        parser.add_argument("--out", default=os.environ.get(ENV_OUTPUT_DIR),
+                            help=f"Folder for result tables (env {ENV_OUTPUT_DIR}).")
+    parser.add_argument("--columns", default=None,
+                        help=f"JSON file mapping Columns fields to your column names (env {ENV_COLUMNS}).")
+
+
+def paths_from_args(args: argparse.Namespace) -> Paths:
+    def p(name):
+        value = getattr(args, name, None)
+        return Path(value).expanduser().resolve() if value else None
+
+    return Paths(
+        data_root=p("data_root"),
+        runs_root=p("runs_root"),
+        output_dir=p("out"),
+        file_template=getattr(args, "file_template", "{part}.csv"),
+        external_template=getattr(args, "external_template", "external.csv"),
+        split_dir_template=getattr(args, "split_dir_template", "split_{split}"),
+    )
diff --git a/evaluation/data.py b/evaluation/data.py
new file mode 100644
index 0000000..0d10679
--- /dev/null
+++ b/evaluation/data.py
@@ -0,0 +1,268 @@
+"""Turn a per-visit table into the ragged sequences DynamicDeepHit reads.
+
+The input is one row per (patient, visit), already preprocessed (imputed, scaled),
+with the outcome repeated on every row: time from the patient's FIRST visit to the
+event or censoring, and the event code. Three views are built from it:
+
+  sequences       one sample per patient, their full observed history; the landmark
+                  s is the day of their last visit (training, and "score at last visit")
+  every visit     one sample per (patient, visit), history truncated at that visit;
+                  s is that visit's day ("score at every visit")
+  fixed landmark  one sample per (patient, landmark day) for patients still at risk,
+                  history up to the landmark plus a forward-filled snapshot row dated
+                  on the landmark ("score everyone at day 0, 90, 180, ...")
+
+Day 0 is each patient's first visit throughout.
+"""
+from __future__ import annotations
+
+from dataclasses import dataclass
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+
+from .config import Columns
+
+
+@dataclass
+class SequenceDataset:
+    x: np.ndarray              # (N,) object array of (visits_i, features) arrays
+    t: np.ndarray              # (N,) object array of per-visit outcome times
+    e: np.ndarray              # (N,) object array of per-visit event codes
+    sample_index: pd.DataFrame  # one row per sample: who, which visit, outcome
+    last_obs_days: np.ndarray  # landmark s per sample, days from the first visit
+    visit_days: np.ndarray | None = None  # (N,) object array of each visit's day
+
+
+def ragged(arrays: list[np.ndarray]) -> np.ndarray:
+    """1-D object array of per-sample arrays.
+
+    np.array(list, dtype=object) stacks the arrays into a dense block when every
+    sample has the same number of visits (always, for a first-visit-only cohort),
+    which the model's padding code cannot read.
+    """
+    out = np.empty(len(arrays), dtype=object)
+    for i, a in enumerate(arrays):
+        out[i] = a
+    return out
+
+
+def read_visits(path: str | Path, cols: Columns) -> pd.DataFrame:
+    """Read one per-visit table, check the required columns, and sort it."""
+    df = pd.read_csv(path, low_memory=False)
+    missing = [c for c in (cols.id, cols.date, cols.visit, cols.time, cols.event) if c not in df]
+    if missing:
+        raise ValueError(f"{path} is missing required column(s) {missing}; "
+                         "map your names with --columns.")
+    df[cols.date] = pd.to_datetime(df[cols.date])
+    return df.sort_values([cols.id, cols.date, cols.visit]).reset_index(drop=True)
+
+
+def predictor_columns(df: pd.DataFrame, cols: Columns) -> list[str]:
+    """Every column that is not pipeline metadata, in file order."""
+    skip = cols.non_predictors()
+    return [c for c in df.columns if c not in skip]
+
+
+def _visit_days(patient: pd.DataFrame, cols: Columns) -> np.ndarray:
+    dates = patient[cols.date]
+    return (dates - dates.iloc[0]).dt.days.to_numpy(dtype=float)
+
+
+def build_sequences(df: pd.DataFrame, features: list[str], cols: Columns) -> SequenceDataset:
+    """One sample per patient, using their complete visit sequence."""
+    x, t, e, last_obs, visit_days, rows = [], [], [], [], [], []
+    for pid, patient in df.sort_values([cols.id, cols.date, cols.visit]).groupby(cols.id, sort=False):
+        patient = patient.reset_index(drop=True)
+        days = _visit_days(patient, cols)
+        x.append(patient[features].to_numpy(dtype=float))
+        t.append(patient[cols.time].to_numpy(dtype=float))
+        e.append(patient[cols.event].to_numpy(dtype=float))
+        visit_days.append(days)
+        last_obs.append(days[-1])
+        last = patient.iloc[-1]
+        rows.append({cols.id: pid, cols.visit: last[cols.visit], cols.date: last[cols.date],
+                     cols.time: float(last[cols.time]), cols.event: int(last[cols.event]),
+                     "num_visits": len(patient)})
+    return SequenceDataset(x=ragged(x), t=ragged(t), e=ragged(e),
+                           sample_index=pd.DataFrame(rows),
+                           last_obs_days=np.array(last_obs, dtype=float),
+                           visit_days=ragged(visit_days))
+
+
+def build_every_visit(df: pd.DataFrame, features: list[str], cols: Columns) -> SequenceDataset:
+    """One sample per (patient, visit), with the input truncated at that visit.
+
+    The outcome is the patient-level one (time from the first visit) on every row;
+    entry_time is the visit's day and time_ref the time remaining from it.
+    """
+    x, t, e, last_obs, rows = [], [], [], [], []
+    for pid, patient in df.sort_values([cols.id, cols.date, cols.visit]).groupby(cols.id, sort=False):
+        patient = patient.reset_index(drop=True)
+        full_x = patient[features].to_numpy(dtype=float)
+        days = _visit_days(patient, cols)
+        terminal_t = float(patient.iloc[-1][cols.time])
+        terminal_e = int(patient.iloc[-1][cols.event])
+        for n in range(len(patient)):
+            x.append(full_x[: n + 1])
+            t.append(np.array([terminal_t]))
+            e.append(np.array([terminal_e]))
+            last_obs.append(days[n])
+            rows.append({cols.id: pid, cols.visit: patient.iloc[n][cols.visit],
+                         cols.date: patient.iloc[n][cols.date],
+                         cols.time: terminal_t, cols.event: terminal_e,
+                         "entry_time": int(days[n]), "time_ref": terminal_t - days[n]})
+    return SequenceDataset(x=ragged(x), t=ragged(t), e=ragged(e),
+                           sample_index=pd.DataFrame(rows),
+                           last_obs_days=np.array(last_obs, dtype=float))
+
+
+def build_fixed_landmarks(
+    table: pd.DataFrame,
+    features: list[str],
+    cols: Columns,
+    landmark_days: list[int] | None = None,
+    first_visit_only: bool = False,
+) -> SequenceDataset:
+    """One sample per (patient, landmark day) from a landmark table.
+
+    The table (see build_landmark_table) holds each patient's real visits
+    (is_landmark = 0; the first visit has landmark_day = 0) and one synthetic
+    snapshot per later landmark day (is_landmark = 1), forward-filled from the most
+    recent real visit. The input at landmark L is every real visit on or before
+    day L, followed by the snapshot for L. s = L, so the CIF is conditioned on
+    surviving to the landmark.
+
+    first_visit_only gives every landmark the first visit alone as input (for a
+    model trained on first visits only); the at-risk set and s are unchanged, so
+    the rows pair one-for-one with a dynamic model's.
+    """
+    for c in ("landmark_day", "is_landmark"):
+        if c not in table:
+            raise ValueError(f"landmark table has no '{c}' column; build it with build_landmark_table")
+    table = table.copy()
+    table[cols.date] = pd.to_datetime(table[cols.date])
+    table = table.sort_values([cols.id, cols.date]).reset_index(drop=True)
+    wanted = None if landmark_days is None else {float(d) for d in landmark_days}
+
+    x, t, e, last_obs, rows = [], [], [], [], []
+    for pid, patient in table.groupby(cols.id, sort=False):
+        patient = patient.reset_index(drop=True)
+        real = patient[patient["is_landmark"] == 0].sort_values(cols.date).reset_index(drop=True)
+        first_date = real.iloc[0][cols.date]
+        terminal_t = float(patient[cols.time].iloc[0])
+        terminal_e = int(patient[cols.event].iloc[0])
+        days = sorted(patient.loc[patient["landmark_day"].notna(), "landmark_day"].unique())
+        if wanted is not None:
+            days = [d for d in days if float(d) in wanted]
+        for day in days:
+            day = float(day)
+            landmark_date = first_date + pd.Timedelta(days=int(day))
+            real_before = real[real[cols.date] <= landmark_date]
+            snapshot = patient[(patient["landmark_day"] == day) & (patient["is_landmark"] == 1)]
+            if first_visit_only:
+                seq = real.iloc[[0]].copy()
+            elif snapshot.empty:          # day 0: the first visit is the prediction point
+                seq = real_before.copy()
+            else:
+                seq = pd.concat([real_before, snapshot], ignore_index=True)
+            seq = seq.sort_values(cols.date).reset_index(drop=True)
+            last_real = seq.iloc[-1] if first_visit_only else real_before.iloc[-1]
+
+            x.append(seq[features].to_numpy(dtype=float))
+            t.append(np.array([terminal_t]))
+            e.append(np.array([terminal_e]))
+            last_obs.append(day)
+            rows.append({cols.id: pid, cols.visit: last_real[cols.visit],
+                         cols.time: terminal_t, cols.event: terminal_e,
+                         "landmark_day": int(day),
+                         "landmark_date": landmark_date.date().isoformat(),
+                         "num_visits_used": len(seq)})
+    return SequenceDataset(x=ragged(x), t=ragged(t), e=ragged(e),
+                           sample_index=pd.DataFrame(rows),
+                           last_obs_days=np.array(last_obs, dtype=float))
+
+
+def build_landmark_table(
+    df: pd.DataFrame,
+    cols: Columns,
+    landmark_days: list[int],
+    reanchor_cols: list[str] | None = None,
+    gap_col: str | None = None,
+    no_ffill: list[str] | None = None,
+) -> pd.DataFrame:
+    """Add one forward-filled snapshot row per patient and landmark day.
+
+    For every patient and every L in landmark_days a row dated first_visit + L is
+    merged into their visit history. Snapshot rows (is_landmark = 1; landmark day 0
+    coincides with the first visit, which is real) take every covariate from the
+    most recent real visit, except:
+
+    - reanchor_cols, "days since X" columns, which are advanced by the gap between
+      that visit and the landmark instead of being copied;
+    - gap_col, if given, which becomes the number of days since the last real visit.
+
+    Rows on or after the patient's event/censoring date are dropped, so a patient
+    appears at a landmark only while still at risk. Run this on the per-visit table
+    BEFORE scaling, so the day arithmetic is in days.
+    """
+    reanchor_cols = list(reanchor_cols or [])
+    df = df.copy()
+    df[cols.date] = pd.to_datetime(df[cols.date], errors="coerce")
+    df = df.sort_values([cols.id, cols.date]).reset_index(drop=True)
+    df["_date_orig"] = df[cols.date]
+    if cols.event_date in df:
+        event_date = pd.to_datetime(df[cols.event_date])
+    else:
+        first = df.groupby(cols.id)[cols.date].transform("min")
+        event_date = first + pd.to_timedelta(df[cols.time], unit="D")
+    df["_event_date"] = event_date
+
+    first_dates = (df.groupby(cols.id, as_index=False)[cols.date].min()
+                     .rename(columns={cols.date: "_first_date"}))
+    frames = []
+    for day in landmark_days:
+        tmp = first_dates.copy()
+        tmp["landmark_day"] = day
+        tmp[cols.date] = tmp["_first_date"] + pd.to_timedelta(day, unit="D")
+        frames.append(tmp)
+    base = pd.concat(frames, ignore_index=True)
+    base["is_prediction_row"] = 1
+
+    keep_as_is = {cols.id, cols.date, cols.visit, "is_prediction_row", "is_landmark",
+                  "landmark_day", *reanchor_cols, *(no_ffill or [])}
+    if gap_col:
+        keep_as_is.add(gap_col)
+    ffill_cols = [c for c in df.columns if c not in keep_as_is]
+
+    lm = base.merge(df, on=[cols.id, cols.date], how="outer")
+    lm["is_prediction_row"] = lm["is_prediction_row"].fillna(0).astype(int)
+    lm = lm.sort_values([cols.id, cols.date, "landmark_day"]).reset_index(drop=True)
+    lm["is_landmark"] = lm[cols.visit].isna().astype(int)   # no real visit on that date
+    lm[ffill_cols] = lm.groupby(cols.id, sort=False)[ffill_cols].ffill()
+
+    lm = lm.loc[lm[cols.date] < lm["_event_date"]].copy()
+    lm["_gap_to_landmark"] = (lm[cols.date] - lm["_date_orig"]).dt.days.astype(float)
+    snap = lm["is_landmark"].eq(1)
+    if gap_col:
+        lm.loc[snap, gap_col] = lm.loc[snap, "_gap_to_landmark"]
+    for c in reanchor_cols:
+        at_visit = lm.groupby([cols.id, "_date_orig"])[c].transform("first")
+        lm[c] = at_visit + lm["_gap_to_landmark"]
+    return lm.drop(columns=["_first_date", "_event_date"])
+
+
+def choose_time_grid(train: SequenceDataset, val: SequenceDataset | None,
+                     horizons: list[int], num_quantiles: int = 100) -> list[float]:
+    """Discretisation grid: quantiles of train (+ val) outcome times, plus the horizons.
+
+    Including the validation times keeps validation outcomes inside the grid.
+    """
+    times = np.array([seq[-1] for seq in train.t], dtype=float)
+    if val is not None:
+        times = np.concatenate([times, np.array([seq[-1] for seq in val.t], dtype=float)])
+    grid = set(float(v) for v in np.quantile(times, np.linspace(0.0, 1.0, num_quantiles + 1)[1:]))
+    grid.update(float(h) for h in horizons)
+    grid.add(float(times.max()))
+    return [0.0] + sorted(v for v in grid if v > 0)
diff --git a/evaluation/evaluate.py b/evaluation/evaluate.py
new file mode 100644
index 0000000..8768244
--- /dev/null
+++ b/evaluation/evaluate.py
@@ -0,0 +1,215 @@
+"""Score trained arms against each other on three populations, raw (no recalibration).
+
+  landmarks  internal test split, every patient still at risk at each fixed landmark
+             day; metrics per split, then mean and SD over splits, and each arm's
+             per-split difference from the reference arm (Wilcoxon signed-rank,
+             Benjamini-Hochberg within each horizon and metric)
+  visits     internal test split, every visit, grouped by days since the first visit;
+             metrics per split, then mean and SD over splits
+  external   external cohort at each patient's anchor, the split models averaged into
+             one ensemble; strata by history; patient bootstrap for each arm and for
+             each arm's paired difference from the reference
+
+Arms are compared on identical rows: a row any arm cannot score (NaN risk) is
+dropped for every arm. A comparator (for example an existing clinical score) is
+evaluated as one more arm: write its predicted risks into the same files an arm
+writes (risk_<h>d columns) and pass its folder with --arm; rows it does not
+cover then drop out for everyone, keeping the comparison paired.
+
+Metric definitions: evaluation/metrics.py.
+
+    python -m evaluation.evaluate --arm band --arm "1/V_i=patient_weighting" \\
+        --reference band --stage landmarks visits external
+"""
+from __future__ import annotations
+
+import argparse
+
+import numpy as np
+import pandas as pd
+
+from .config import (EVENT_OF_INTEREST, EVERY_VISIT_FILE, EXTERNAL_ANCHOR_FILE,
+                     FIXED_LANDMARK_FILE, Columns, add_path_arguments, paths_from_args)
+from .metrics import benjamini_hochberg, bootstrap, evaluate, percentile_ci, wilcoxon_paired
+
+VISIT_BANDS = [("first visit", 0, 0), ("days 1-90", 1, 90), ("days 91-180", 91, 180),
+               ("days 181-365", 181, 365), (">365 days", 366, 10**6), ("all visits", 0, 10**6)]
+
+
+def parse_arms(specs: list[str]) -> dict[str, str]:
+    """'Label=folder' or 'folder' (label = folder name), in the order given."""
+    arms = {}
+    for spec in specs:
+        label, _, folder = spec.rpartition("=")
+        arms[label or folder] = folder
+    return arms
+
+
+def read_predictions(paths, folder: str, split: int, name: str) -> pd.DataFrame:
+    df = pd.read_csv(paths.run_dir(folder, split) / name)
+    return df[df.time_ref > 0].reset_index(drop=True)
+
+
+def n_events(t, e, h):
+    return int(((t <= h) & (e == EVENT_OF_INTEREST)).sum())
+
+
+# ------------------------------------------------------------------ landmarks
+def stage_landmarks(arms, paths, cols, splits, horizons, landmarks):
+    rows = []
+    for split in splits:
+        frames = {a: read_predictions(paths, f, split, FIXED_LANDMARK_FILE.format(part="test"))
+                  for a, f in arms.items()}
+        ref = next(iter(frames.values()))
+        for h in horizons:
+            for lm in landmarks:
+                base = (ref[ref.landmark_day == lm].drop_duplicates(cols.id).sort_values(cols.id))
+                keys = base[cols.id].to_numpy()
+                P = {a: f[f.landmark_day == lm].drop_duplicates(cols.id).set_index(cols.id)
+                        .reindex(keys)[f"risk_{h}d"].to_numpy() for a, f in frames.items()}
+                t, e = base.time_ref.to_numpy(float), base[cols.event].to_numpy()
+                ok = np.all(np.isfinite(np.vstack(list(P.values()))), axis=0)
+                for a, r in evaluate(t[ok], e[ok], {a: p[ok] for a, p in P.items()}, h).items():
+                    rows.append(dict(split=split, horizon=h, landmark=lm, n=int(ok.sum()),
+                                     events=n_events(t[ok], e[ok], h), arm=a, **r))
+        print(f"landmarks: split {split}", flush=True)
+    return pd.DataFrame(rows)
+
+
+def summarise_splits(per_split: pd.DataFrame, keys: list[str], reference: str) -> pd.DataFrame:
+    """Mean and SD over splits; per-split paired differences from the reference arm."""
+    metrics = [m for m in ("auc", "ipa", "op", "slope", "citl", "top_op", "obs", "pred")
+               if m in per_split]
+    summary = per_split.groupby(keys + ["arm"], sort=False)[metrics + ["n", "events"]].agg(["mean", "std"])
+    summary.columns = [f"{m}_{s}" if s == "std" else m for m, s in summary.columns]
+    summary = summary.reset_index()
+    tests = []
+    for key, cell in per_split.groupby(keys, sort=False):
+        key = key if isinstance(key, tuple) else (key,)
+        r = cell[cell.arm == reference].set_index("split")
+        for arm, g in cell[cell.arm != reference].groupby("arm", sort=False):
+            g = g.set_index("split")
+            for m in ("auc", "ipa"):
+                d = (g[m] - r[m]).dropna().to_numpy()
+                tests.append(dict(zip(keys, key), arm=arm, metric=m, delta=d.mean() if len(d) else np.nan,
+                                  delta_sd=d.std(ddof=1) if len(d) > 1 else np.nan,
+                                  splits_up=int((d > 0).sum()), p=wilcoxon_paired(d)))
+    tests = pd.DataFrame(tests)
+    if not tests.empty:
+        family = [k for k in keys if k == "horizon"] + ["metric"]
+        tests["q_bh"] = tests.groupby(family)["p"].transform(lambda p: benjamini_hochberg(p.to_numpy()))
+        wide = tests.pivot_table(index=keys + ["arm"], columns="metric",
+                                 values=["delta", "delta_sd", "splits_up", "p", "q_bh"], sort=False)
+        wide.columns = [f"{m}_{s}_vs_ref" for s, m in wide.columns]
+        summary = summary.merge(wide.reset_index(), on=keys + ["arm"], how="left")
+    return summary
+
+
+# ------------------------------------------------------------------ visits
+def stage_visits(arms, paths, cols, splits, horizons):
+    rows, key = [], [cols.id, cols.visit]
+    for split in splits:
+        frames = {a: read_predictions(paths, f, split, EVERY_VISIT_FILE.format(part="test"))
+                  for a, f in arms.items()}
+        base = next(iter(frames.values()))
+        day = base.entry_time.to_numpy(float)
+        t, e = base.time_ref.to_numpy(float), base[cols.event].to_numpy()
+        for h in horizons:
+            P = {a: base[key].merge(f[key + [f"risk_{h}d"]], on=key, how="left")[f"risk_{h}d"].to_numpy()
+                 for a, f in frames.items()}
+            for label, lo, hi in VISIT_BANDS:
+                ok = (day >= lo) & (day <= hi) & np.all(np.isfinite(np.vstack(list(P.values()))), axis=0)
+                for a, r in evaluate(t[ok], e[ok], {a: p[ok] for a, p in P.items()}, h).items():
+                    rows.append(dict(split=split, horizon=h, stratum=label, n=int(ok.sum()),
+                                     events=n_events(t[ok], e[ok], h), arm=a, **r))
+        print(f"visits: split {split}", flush=True)
+    return pd.DataFrame(rows)
+
+
+# ------------------------------------------------------------------ external
+def stage_external(arms, paths, cols, splits, horizons, reference, n_boot, seed):
+    first = next(iter(arms.values()))
+    ref = read_predictions(paths, first, splits[0], EXTERNAL_ANCHOR_FILE).sort_values(cols.id)
+    ids = ref[cols.id].to_numpy()
+    t, e = ref.time_ref.to_numpy(float), ref[cols.event].to_numpy()
+    day, n_visits = ref.entry_time.to_numpy(float), ref.num_visits.to_numpy()
+    ensemble = {}
+    for h in horizons:
+        for a, folder in arms.items():
+            ensemble[(h, a)] = np.vstack([
+                read_predictions(paths, folder, s, EXTERNAL_ANCHOR_FILE).set_index(cols.id)
+                .reindex(ids)[f"risk_{h}d"].to_numpy() for s in splits]).mean(0)
+    covered = np.all(np.isfinite(np.vstack(list(ensemble.values()))), axis=0)
+    strata = {"all": covered, "first visit": covered & (n_visits == 1),
+              "history <=365 d": covered & (n_visits > 1) & (day <= 365),
+              ">365 d": covered & (day > 365)}
+    rows = []
+    for h in horizons:
+        for name, mask in strata.items():
+            preds = {a: ensemble[(h, a)][mask] for a in arms}
+            tm, em = t[mask], e[mask]
+            point = evaluate(tm, em, preds, h)
+            boots = bootstrap(tm, em, preds, h, n_boot=n_boot, seed=seed)
+            for a, r in point.items():
+                row = dict(horizon=h, stratum=name, n=int(mask.sum()), events=n_events(tm, em, h),
+                           arm=a, **r)
+                for m in ("auc", "ipa"):
+                    row[f"{m}_lo"], row[f"{m}_hi"], _ = percentile_ci(boots[a][m])
+                    if a != reference:
+                        lo, hi, p = percentile_ci(boots[a][m] - boots[reference][m])
+                        row.update({f"{m}_delta_vs_ref": r[m] - point[reference][m],
+                                    f"{m}_delta_lo": lo, f"{m}_delta_hi": hi, f"{m}_delta_p": p})
+                rows.append(row)
+            print(f"external: horizon {h}, {name}", flush=True)
+    return pd.DataFrame(rows)
+
+
+def main(argv=None) -> None:
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    add_path_arguments(ap, data=False)
+    ap.add_argument("--arm", action="append", required=True,
+                    help="'Label=folder' (folder under --runs-root, or absolute); repeatable.")
+    ap.add_argument("--reference", default=None, help="Arm label every other arm is compared with "
+                                                      "(default: the first --arm).")
+    ap.add_argument("--stage", nargs="+", choices=["landmarks", "visits", "external"],
+                    default=["landmarks", "visits", "external"])
+    ap.add_argument("--splits", type=int, nargs="+", default=list(range(10)))
+    ap.add_argument("--horizons", type=int, nargs="+", default=[90, 180])
+    ap.add_argument("--landmarks", type=int, nargs="+", default=[0, 90, 180, 270])
+    ap.add_argument("--n-boot", type=int, default=1000)
+    ap.add_argument("--seed", type=int, default=7)
+    args = ap.parse_args(argv)
+
+    cols, paths = Columns.load(args.columns), paths_from_args(args)
+    arms = parse_arms(args.arm)
+    reference = args.reference or next(iter(arms))
+    if reference not in arms:
+        raise SystemExit(f"--reference {reference!r} is not one of {list(arms)}")
+    # The reference arm defines the rows and their order.
+    arms = {reference: arms[reference], **{a: f for a, f in arms.items() if a != reference}}
+
+    pd.set_option("display.width", 200)
+    if "landmarks" in args.stage:
+        per_split = stage_landmarks(arms, paths, cols, args.splits, args.horizons, args.landmarks)
+        per_split.to_csv(paths.out("landmarks_per_split.csv"), index=False)
+        summary = summarise_splits(per_split, ["horizon", "landmark"], reference)
+        summary.to_csv(paths.out("landmarks_summary.csv"), index=False)
+        print(summary.pivot_table(index=["horizon", "landmark"], columns="arm",
+                                  values=["auc", "op"], sort=False).round(3).to_string())
+    if "visits" in args.stage:
+        per_split = stage_visits(arms, paths, cols, args.splits, args.horizons)
+        per_split.to_csv(paths.out("visits_per_split.csv"), index=False)
+        summary = summarise_splits(per_split, ["horizon", "stratum"], reference)
+        summary.to_csv(paths.out("visits_summary.csv"), index=False)
+        print(summary.pivot_table(index=["horizon", "stratum"], columns="arm",
+                                  values=["auc", "op"], sort=False).round(3).to_string())
+    if "external" in args.stage:
+        ext = stage_external(arms, paths, cols, args.splits, args.horizons, reference,
+                             args.n_boot, args.seed)
+        ext.to_csv(paths.out("external.csv"), index=False)
+        print(ext.pivot_table(index=["horizon", "stratum"], columns="arm",
+                              values=["auc", "op"], sort=False).round(3).to_string())
+
+
+if __name__ == "__main__":
+    main()
diff --git a/evaluation/history_ablation.py b/evaluation/history_ablation.py
new file mode 100644
index 0000000..1a81e47
--- /dev/null
+++ b/evaluation/history_ablation.py
@@ -0,0 +1,136 @@
+"""Does the model actually use a patient's earlier visits? Within-patient ablation.
+
+Comparing patients with long records against patients with short ones cannot answer
+this: they are different patients. This scores the SAME patients with the SAME
+trained models on a changed record:
+
+  full        the patient's record up to the anchor, landmark s = anchor day
+              (must equal external_anchor_predictions.csv; checked when present)
+  masked      the anchor visit alone, s still the anchor day. What the encoder
+              extracts from earlier visits, holding the time frame fixed. The anchor
+              row keeps any history it already summarises (forward-filled labs, gaps).
+  reanchored  the anchor visit alone, s = 0, and --gap-col reset to the value it has
+              on first visits: the patient as if seen for the first time.
+              Only with --gap-col.
+
+Patients with at least --min-visits visits at the anchor are used. Outcomes come
+from `full` and are shared, so every contrast is paired. Each variant is scored as
+the ensemble mean of the split models, with a patient bootstrap of each contrast
+(full - masked, full - reanchored, masked - reanchored).
+
+Writes ablation_predictions.csv (patient level: keep it with the data) and
+ablation_summary.csv / ablation_contrasts.csv (aggregates).
+
+    python -m evaluation.history_ablation --arm band --gap-col days_since_last_obs
+"""
+from __future__ import annotations
+
+import argparse
+
+import numpy as np
+import pandas as pd
+from scipy.stats import spearmanr
+
+from .config import EXTERNAL_ANCHOR_FILE, Columns, add_path_arguments, paths_from_args
+from .data import build_sequences, read_visits
+from .metrics import bootstrap, evaluate, percentile_ci
+from .model import load_run, predict_at_horizons, set_determinism
+
+PAIRS = [("full", "masked"), ("full", "reanchored"), ("masked", "reanchored")]
+
+
+def score(model, frame, predictors, horizons, cols, s=None):
+    ds = build_sequences(frame, predictors, cols)
+    if s is not None:
+        ds.last_obs_days = ds.sample_index[cols.id].map(s).to_numpy(dtype=float)
+    return predict_at_horizons(model, ds, horizons, cols).set_index(cols.id)
+
+
+def main(argv=None) -> None:
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    add_path_arguments(ap)
+    ap.add_argument("--arm", required=True)
+    ap.add_argument("--splits", type=int, nargs="+", default=list(range(10)))
+    ap.add_argument("--min-visits", type=int, default=2)
+    ap.add_argument("--gap-col", default=None, help="Days-since-last-visit predictor, for 'reanchored'.")
+    ap.add_argument("--horizons", type=int, nargs="+", default=[90, 180])
+    ap.add_argument("--n-boot", type=int, default=1000)
+    ap.add_argument("--seed", type=int, default=19)
+    args = ap.parse_args(argv)
+    cols, paths = Columns.load(args.columns), paths_from_args(args)
+    variants = ["full", "masked"] + (["reanchored"] if args.gap_col else [])
+
+    rows = []
+    for split in args.splits:
+        run_dir = paths.run_dir(args.arm, split)
+        set_determinism(42)
+        model, meta = load_run(run_dir, cols, read_visits(paths.split_file(split, "train"), cols))
+        predictors, horizons = meta["predictors"], sorted(set(meta["horizons"]) | set(args.horizons))
+        cohort = read_visits(paths.external_file(split), cols)
+        n_visits = cohort.groupby(cols.id).size()
+        hist = cohort[cohort[cols.id].isin(n_visits[n_visits >= args.min_visits].index)]
+        anchor = hist.groupby(cols.id, sort=False).tail(1)
+
+        out = {"full": score(model, hist, predictors, horizons, cols)}
+        out["masked"] = score(model, anchor, predictors, horizons, cols, s=out["full"]["entry_time"])
+        if args.gap_col:
+            first_gap = cohort.groupby(cols.id, sort=False).head(1)[args.gap_col]
+            if first_gap.round(6).nunique() != 1:
+                raise SystemExit(f"--gap-col {args.gap_col} is not constant on first visits")
+            fresh = anchor.copy()
+            fresh[args.gap_col] = float(first_gap.iloc[0])
+            out["reanchored"] = score(model, fresh, predictors, horizons, cols)
+
+        saved = run_dir / EXTERNAL_ANCHOR_FILE
+        if saved.exists():
+            s = pd.read_csv(saved).set_index(cols.id).reindex(out["full"].index)
+            drift = max(float(np.nanmax(np.abs(out["full"][f"risk_{h}d"] - s[f"risk_{h}d"])))
+                        for h in meta["horizons"])
+            if drift > 1e-6:
+                raise SystemExit(f"split {split}: 'full' differs from {saved.name} by {drift:.2e}")
+        for v in variants:
+            frame = out[v].reindex(out["full"].index)
+            for pid, r in frame.iterrows():
+                rows.append({"split": split, "variant": v, cols.id: pid,
+                             "n_visits": int(n_visits[pid]),
+                             "time_ref": float(out["full"].loc[pid, "time_ref"]),
+                             "event": int(out["full"].loc[pid, cols.event]),
+                             **{f"risk_{h}d": float(r[f"risk_{h}d"]) for h in args.horizons}})
+        print(f"split {split}: {len(out['full'])} patients with >= {args.min_visits} visits", flush=True)
+
+    preds = pd.DataFrame(rows)
+    preds.to_csv(paths.out("ablation_predictions.csv"), index=False)
+
+    ref = preds[(preds.variant == "full") & (preds.split == args.splits[0])].sort_values(cols.id)
+    ids, t, e = ref[cols.id].to_numpy(), ref.time_ref.to_numpy(float), ref.event.to_numpy()
+    summary, contrasts = [], []
+    for h in args.horizons:
+        ens = {v: preds[preds.variant == v].pivot(index="split", columns=cols.id, values=f"risk_{h}d")
+               [ids].to_numpy().mean(0) for v in variants}
+        point = evaluate(t, e, ens, h)
+        boots = bootstrap(t, e, ens, h, n_boot=args.n_boot, seed=args.seed, metrics=("auc", "ipa", "op"))
+        for v, r in point.items():
+            row = dict(horizon=h, variant=v, n=len(t), **r)
+            for m in ("auc", "ipa", "op"):
+                row[f"{m}_lo"], row[f"{m}_hi"], _ = percentile_ci(boots[v][m])
+            summary.append(row)
+        for a, b in PAIRS:
+            if a not in ens or b not in ens:
+                continue
+            row = dict(horizon=h, contrast=f"{a} - {b}")
+            for m in ("auc", "ipa", "op"):
+                lo, hi, p = percentile_ci(boots[a][m] - boots[b][m])
+                row.update({f"d_{m}": point[a][m] - point[b][m], f"d_{m}_lo": lo, f"d_{m}_hi": hi,
+                            f"d_{m}_p": p})
+            if (a, b) == ("full", "masked"):
+                deciles = [pd.qcut(ens[v], 10, labels=False, duplicates="drop") for v in (a, b)]
+                row.update(spearman=float(spearmanr(ens[a], ens[b]).correlation),
+                           share_changing_decile=float(np.mean(deciles[0] != deciles[1])))
+            contrasts.append(row)
+    pd.DataFrame(summary).to_csv(paths.out("ablation_summary.csv"), index=False)
+    pd.DataFrame(contrasts).to_csv(paths.out("ablation_contrasts.csv"), index=False)
+    print(pd.DataFrame(contrasts).round(3).to_string(index=False))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/evaluation/make_synthetic.py b/evaluation/make_synthetic.py
new file mode 100644
index 0000000..d1958f7
--- /dev/null
+++ b/evaluation/make_synthetic.py
@@ -0,0 +1,120 @@
+"""Write a synthetic cohort in the layout the evaluation scripts read.
+
+No real patient is involved. Each simulated patient has a categorical type, an
+age, and a lab value measured at every visit; the event of interest (1) and a
+competing event (2) arrive at exponential times, and follow-up is censored
+administratively. Two properties of real follow-up data are built in:
+
+- visits stop at the event, so patients with an early event have few visits
+  (the informative visit pattern that makes 1/V_i landmark weighting biased);
+- the lab rises over the weeks before an event of interest, so later visits
+  carry information the first visit does not (what a dynamic model can use).
+
+    python -m evaluation.make_synthetic --out /tmp/ddh_synthetic
+
+writes <out>/split_<k>/{train,val,test,train_landmark,test_landmark,external}.csv
+and <out>/columns.json, which tells the other scripts the column names.
+"""
+from __future__ import annotations
+
+import argparse
+import json
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+
+from .config import Columns
+from .data import build_landmark_table
+
+CONTINUOUS = ["age", "lab", "days_since_last_obs"]
+COLUMNS = Columns(id="patient_id", date="visit_date", visit="visit", time="time_to_event",
+                  event="event", event_date="event_date", categorical="cancer_type")
+LANDMARK_DAYS = [0, 30, 60, 90, 120, 150, 180, 270]
+
+
+def simulate(n: int, rng: np.random.Generator, start: str, span_days: int, first_id: int) -> pd.DataFrame:
+    type_effect = np.array([-0.6, -0.2, 0.0, 0.4, 0.9])
+    rows = []
+    for i in range(n):
+        ctype = int(rng.integers(0, len(type_effect)))
+        age = rng.normal(62, 12)
+        lab0 = rng.normal(0, 1)
+        frailty = rng.normal(0, 0.4)
+        t1 = rng.exponential(1 / (0.0003 * np.exp(type_effect[ctype] + 0.5 * lab0 + frailty)))
+        t2 = rng.exponential(1 / (0.0003 * np.exp(0.03 * (age - 62))))
+        cens = rng.uniform(150, 900)
+        T = min(t1, t2, cens)
+        event = 1 if T == t1 else 2 if T == t2 else 0
+        first = pd.Timestamp(start) + pd.Timedelta(days=int(rng.integers(0, span_days)))
+
+        day, visit, last_day = 0.0, 1, 0.0
+        while day < T:
+            rise = 1.5 * np.exp(-(t1 - day) / 45) if event == 1 else 0.0
+            rows.append({"patient_id": first_id + i, "visit": visit,
+                         "visit_date": first + pd.Timedelta(days=int(day)),
+                         "time_to_event": float(np.ceil(T)), "event": event,
+                         "event_date": first + pd.Timedelta(days=int(np.ceil(T))),
+                         "cancer_type": ctype, "age": age,
+                         "lab": lab0 + rise + rng.normal(0, 0.3),
+                         "days_since_last_obs": day - last_day})
+            last_day = day
+            day += 7 + rng.exponential(28)
+            visit += 1
+    return pd.DataFrame(rows)
+
+
+def truncate_at_anchor(df: pd.DataFrame, rng: np.random.Generator) -> pd.DataFrame:
+    """Keep each patient's visits up to a random anchor visit (a temporal-validation cohort)."""
+    keep = []
+    for _, g in df.groupby("patient_id", sort=False):
+        keep.append(g.iloc[: int(rng.integers(1, len(g) + 1))])
+    return pd.concat(keep, ignore_index=True)
+
+
+def scale(frames: dict[str, pd.DataFrame], reference: pd.DataFrame) -> dict[str, pd.DataFrame]:
+    mean, sd = reference[CONTINUOUS].mean(), reference[CONTINUOUS].std().replace(0, 1)
+    out = {}
+    for name, f in frames.items():
+        f = f.copy()
+        f[CONTINUOUS] = (f[CONTINUOUS] - mean) / sd
+        out[name] = f
+    return out
+
+
+def main(argv=None) -> None:
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--out", required=True, help="Data root to create.")
+    ap.add_argument("--splits", type=int, default=3)
+    ap.add_argument("--patients", type=int, default=1500)
+    ap.add_argument("--external-patients", type=int, default=500)
+    ap.add_argument("--seed", type=int, default=0)
+    args = ap.parse_args(argv)
+    out = Path(args.out)
+    rng = np.random.default_rng(args.seed)
+
+    cohort = simulate(args.patients, rng, "2018-01-01", 4 * 365, first_id=1)
+    external = truncate_at_anchor(
+        simulate(args.external_patients, rng, "2023-01-01", 365, first_id=10**6), rng)
+    ids = cohort.patient_id.unique()
+    for k in range(args.splits):
+        perm = np.random.default_rng(args.seed + 100 + k).permutation(ids)
+        n_tr, n_va = int(0.7 * len(perm)), int(0.15 * len(perm))
+        parts = {"train": perm[:n_tr], "val": perm[n_tr:n_tr + n_va], "test": perm[n_tr + n_va:]}
+        raw = {name: cohort[cohort.patient_id.isin(p)] for name, p in parts.items()}
+        for name in ("train", "test"):
+            raw[f"{name}_landmark"] = build_landmark_table(
+                raw[name], COLUMNS, LANDMARK_DAYS, gap_col="days_since_last_obs")
+        raw["external"] = external
+        split_dir = out / f"split_{k}"
+        split_dir.mkdir(parents=True, exist_ok=True)
+        for name, frame in scale(raw, raw["train"]).items():
+            frame.to_csv(split_dir / f"{name}.csv", index=False)
+    (out / "columns.json").write_text(json.dumps(COLUMNS.to_dict(), indent=2))
+    events = cohort.groupby("patient_id").event.first().value_counts().to_dict()
+    print(f"{args.patients} patients ({events}), {len(cohort)} visits, "
+          f"{args.splits} splits, {external.patient_id.nunique()} external patients -> {out}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/evaluation/metrics.py b/evaluation/metrics.py
new file mode 100644
index 0000000..707f239
--- /dev/null
+++ b/evaluation/metrics.py
@@ -0,0 +1,209 @@
+"""Discrimination and calibration of a risk by horizon h, scored from a landmark.
+
+Every function takes, per row (one patient at one landmark):
+
+    t   time from the landmark to the event or censoring (days)
+    e   0 censored, 1 event of interest, 2 competing event
+    p   predicted risk of the event of interest by h
+
+and treats the competing event as a competing event, never as censoring.
+
+observed         Aalen-Johansen cumulative incidence of cause 1 at h.
+O:E              observed / mean(p). 1 is calibrated; below 1 the model over-predicts.
+AUC              IPCW cumulative/dynamic AUC with a competing event (Blanche et al.,
+                 Stat Med 2013). Cases had the event of interest by h, each weighted
+                 1/G(T-); controls were still event-free after h, weighted 1/G(h).
+                 Rows censored, or with a competing event, before h are in neither
+                 group; ties in p count 1/2. G is the reverse Kaplan-Meier censoring
+                 distribution, with a day's events leaving its risk set before that
+                 day's censorings (the prodlim / timeROC convention). Identical to R
+                 timeROC's AUC_1 except that timeROC counts only events strictly
+                 before h as cases.
+IPCW weights     Graf et al.: 1/G(T-) for a row with any event by h, 1/G(h) for a row
+                 still event-free at h, 0 for a row censored before h. G is the
+                 Kaplan-Meier estimate of the censoring distribution (here a day's
+                 events stay in its risk set, as in the original analysis code; this
+                 differs from the AUC's G only on days with both).
+Brier, IPA       IPCW Brier score of the status at h (y = 1 for the event of interest by
+                 h, y = 0 otherwise, including a competing event by h). IPA = 1 - Brier
+                 / Brier of a null model that predicts `observed` for every row.
+slope, CITL      IPCW-weighted logistic regression of y on logit(p): its slope, and the
+                 intercept with logit(p) as an offset (calibration-in-the-large).
+top-decile O:E   O:E among the 10% of rows with the highest p.
+
+These are the definitions the project's reports were built with; the functions are
+NumPy ports of that code so the numbers reproduce exactly.
+"""
+from __future__ import annotations
+
+import numpy as np
+from scipy import stats
+
+from .config import EVENT_OF_INTEREST as EOI
+
+EPS = 1e-6
+MIN_CASES = 3
+
+
+def logit(p: np.ndarray) -> np.ndarray:
+    p = np.clip(p, EPS, 1 - EPS)
+    return np.log(p / (1 - p))
+
+
+def km_steps(t: np.ndarray, event: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+    """Kaplan-Meier survival of `event` (boolean) at each distinct time."""
+    order = np.argsort(t, kind="mergesort")
+    ts, ev = t[order], event[order].astype(float)
+    uniq, first = np.unique(ts, return_index=True)
+    at_risk = len(ts) - first
+    return uniq, np.cumprod(1 - np.add.reduceat(ev, first) / at_risk)
+
+
+def step_at(uniq: np.ndarray, surv: np.ndarray, x: np.ndarray) -> np.ndarray:
+    pos = np.searchsorted(uniq, x, side="right") - 1
+    return np.where(pos >= 0, surv[np.clip(pos, 0, None)], 1.0)
+
+
+def ipcw_weights(t: np.ndarray, e: np.ndarray, h: float) -> tuple[np.ndarray, np.ndarray]:
+    """Graf weights and the binary status y at h (see module docstring)."""
+    uniq, g = km_steps(t, e == 0)
+    g_h = max(float(step_at(uniq, g, np.array([h], float))[0]), EPS)
+    g_t = np.maximum(step_at(uniq, g, np.maximum(t - EPS, 0)), EPS)
+    w = np.where((t <= h) & (e != 0), 1 / g_t, 0.0) + np.where(t > h, 1 / g_h, 0.0)
+    return w, ((t <= h) & (e == EOI)).astype(float)
+
+
+def aalen_johansen(t: np.ndarray, e: np.ndarray, h: float) -> float:
+    """Cumulative incidence of the event of interest at h, competing events as such."""
+    order = np.argsort(t, kind="mergesort")
+    ts, es = t[order], e[order]
+    uniq, first = np.unique(ts, return_index=True)
+    at_risk = len(ts) - first
+    d_eoi = np.add.reduceat((es == EOI).astype(float), first)
+    d_any = np.add.reduceat((es != 0).astype(float), first)
+    surv = np.cumprod(1 - d_any / at_risk)
+    cif = np.cumsum(np.concatenate([[1.0], surv[:-1]]) * d_eoi / at_risk)
+    pos = np.searchsorted(uniq, h, side="right") - 1
+    return float(cif[pos]) if pos >= 0 else 0.0
+
+
+def censoring_survival(t: np.ndarray, e: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+    """Reverse Kaplan-Meier G(t) = P(C > t), a day's events leaving the risk set first."""
+    order = np.argsort(t, kind="mergesort")
+    ts, es = t[order], e[order]
+    times, first = np.unique(ts, return_index=True)
+    censored = np.add.reduceat((es == 0).astype(float), first)
+    at_risk = (len(ts) - first) - np.add.reduceat((es != 0).astype(float), first)
+    with np.errstate(divide="ignore", invalid="ignore"):
+        step = np.where(at_risk > 0, 1 - censored / at_risk, 1.0)
+    return times, np.cumprod(step)
+
+
+def auc(t: np.ndarray, e: np.ndarray, p: np.ndarray, h: float) -> float:
+    """IPCW cumulative/dynamic AUC at h with a competing event (see module docstring)."""
+    case, ctrl = (t <= h) & (e == EOI), t > h
+    if case.sum() < MIN_CASES or ctrl.sum() == 0:
+        return np.nan
+    times, g = censoring_survival(t, e)
+    w_case = 1.0 / np.maximum(step_at(times, g, np.maximum(t[case] - EPS, 0)), EPS)
+    # Every control carries the same weight 1/G(h), which cancels from the ratio:
+    # each case contributes the share of controls it outranks (ties 1/2), weighted.
+    pk = np.sort(p[ctrl])
+    below = np.searchsorted(pk, p[case], side="left")
+    ties = np.searchsorted(pk, p[case], side="right") - below
+    return float(np.sum(w_case * (below + 0.5 * ties) / len(pk)) / np.sum(w_case))
+
+
+def slope_citl(w: np.ndarray, y: np.ndarray, p: np.ndarray) -> tuple[float, float]:
+    import statsmodels.api as sm
+    lp, k = logit(p), w > 0
+    if k.sum() < 20 or y[k].sum() < 5:
+        return np.nan, np.nan
+    binomial = sm.families.Binomial()
+    fit_citl = sm.GLM(y[k], np.ones((int(k.sum()), 1)), family=binomial, freq_weights=w[k],
+                      offset=lp[k]).fit()
+    citl = float(np.asarray(fit_citl.params)[0])
+    if np.ptp(lp[k]) == 0:          # constant predictions: no slope to estimate
+        return np.nan, citl
+    fit_slope = sm.GLM(y[k], sm.add_constant(lp[k]), family=binomial, freq_weights=w[k]).fit()
+    return float(np.asarray(fit_slope.params)[1]), citl
+
+
+def top_decile_oe(t: np.ndarray, e: np.ndarray, p: np.ndarray, h: float) -> float:
+    top = p >= np.quantile(p, 0.9)
+    return aalen_johansen(t[top], e[top], h) / float(p[top].mean())
+
+
+def evaluate(t: np.ndarray, e: np.ndarray, preds: dict[str, np.ndarray], h: float,
+             with_slope: bool = True, with_top_decile: bool = True) -> dict[str, dict[str, float]]:
+    """All metrics for every predictor in `preds`, on the same rows.
+
+    The IPCW weights and the observed incidence depend only on (t, e), so they are
+    computed once and shared, which keeps predictors strictly paired.
+    """
+    w, y = ipcw_weights(t, e, h)
+    observed = aalen_johansen(t, e, h)
+    brier_null = float(np.mean(w * (y - observed) ** 2))
+    out = {}
+    for name, p in preds.items():
+        brier = float(np.mean(w * (y - p) ** 2))
+        r = dict(auc=auc(t, e, p, h),
+                 ipa=1 - brier / brier_null if brier_null > 0 else np.nan,
+                 brier=brier, op=observed / float(p.mean()), obs=observed, pred=float(p.mean()))
+        if with_slope:
+            r["slope"], r["citl"] = slope_citl(w, y, p)
+        if with_top_decile:
+            r["top_op"] = top_decile_oe(t, e, p, h)
+        out[name] = r
+    return out
+
+
+def bootstrap(t: np.ndarray, e: np.ndarray, preds: dict[str, np.ndarray], h: float,
+              n_boot: int = 1000, seed: int = 7,
+              metrics: tuple[str, ...] = ("auc", "ipa")) -> dict[str, dict[str, np.ndarray]]:
+    """Bootstrap over patients, the same resampled rows for every predictor.
+
+    Returns {name: {metric: one value per kept resample}}. Because all predictors
+    share each resample, values[a] - values[b] is a paired difference. Resamples
+    with fewer than MIN_CASES cases or no controls are skipped (the AUC is
+    undefined there), for every predictor alike.
+    """
+    rng = np.random.default_rng(seed)
+    n = len(t)
+    values = {name: {m: [] for m in metrics} for name in preds}
+    for _ in range(n_boot):
+        k = rng.choice(np.arange(n), n, replace=True)
+        tk, ek = t[k], e[k]
+        if ((tk <= h) & (ek == EOI)).sum() < MIN_CASES or (tk > h).sum() == 0:
+            continue
+        r = evaluate(tk, ek, {name: p[k] for name, p in preds.items()}, h,
+                     with_slope=False, with_top_decile=False)
+        for name in preds:
+            for m in metrics:
+                values[name][m].append(r[name][m])
+    return {name: {m: np.asarray(v, float) for m, v in d.items()} for name, d in values.items()}
+
+
+def percentile_ci(v: np.ndarray) -> tuple[float, float, float]:
+    """95% percentile interval and a two-sided bootstrap p-value for a zero difference."""
+    v = v[np.isfinite(v)]
+    if not len(v):
+        return np.nan, np.nan, np.nan
+    lo, hi = np.percentile(v, [2.5, 97.5])
+    return float(lo), float(hi), float(min(1.0, 2 * min((v <= 0).mean(), (v >= 0).mean())))
+
+
+def wilcoxon_paired(diffs: np.ndarray) -> float:
+    """Wilcoxon signed-rank p-value for per-split paired differences."""
+    d = diffs[np.isfinite(diffs)]
+    if len(d) == 0 or np.abs(d).sum() == 0:
+        return 1.0
+    return float(stats.wilcoxon(d).pvalue)
+
+
+def benjamini_hochberg(p: np.ndarray) -> np.ndarray:
+    p = np.asarray(p, float)
+    order = np.argsort(p)
+    q = np.empty_like(p)
+    q[order] = np.minimum.accumulate((p[order] * len(p) / np.arange(1, len(p) + 1))[::-1])[::-1]
+    return np.clip(q, 0, 1)
diff --git a/evaluation/model.py b/evaluation/model.py
new file mode 100644
index 0000000..77e354c
--- /dev/null
+++ b/evaluation/model.py
@@ -0,0 +1,299 @@
+"""Build, train, save, load and score a DynamicDeepHit model for one split.
+
+Predictions are always made from a landmark s: the raw model output F_k(t) is a
+CIF from day 0, and the quantity a clinician needs at day s is
+
+    F_k(s + h | T > s) = (F_k(s + h) - F_k(s)) / S(s),   S(s) = 1 - sum_k F_k(s)
+
+for each horizon h. condition_cif_on_landmark and predict_at_horizons implement
+exactly that; every evaluation table is built from their output.
+"""
+from __future__ import annotations
+
+import hashlib
+import json
+import os
+import platform
+import random
+import subprocess
+import sys
+from dataclasses import asdict, dataclass, field
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+import pandas as pd
+import torch
+
+from ddh import DynamicDeepHit
+
+from .config import (CHECKPOINT_FILE, COMPETING_EVENT, EVENT_OF_INTEREST, METADATA_FILE,
+                     Columns)
+from .data import SequenceDataset, choose_time_grid
+
+
+@dataclass
+class Architecture:
+    """Network shape. Stored with every run so a checkpoint can be rebuilt exactly."""
+
+    hidden_rnn: int = 32
+    layers_rnn: int = 1
+    typ: str = "GRU"
+    long_layers: list[int] = field(default_factory=lambda: [64])
+    att_layers: list[int] = field(default_factory=lambda: [32])
+    cs_layers: list[int] = field(default_factory=lambda: [64, 32])
+    dropout: float = 0.2
+    embedding_dim: int = 16
+    hazard: bool = False
+    proportional: bool = False
+
+
+@dataclass
+class Objective:
+    """Training objective. Defaults reproduce the library's original loss."""
+
+    alpha: float = 0.01                  # ranking term
+    beta: float = 0.75                   # survival NLL; 1 - alpha - beta goes to reconstruction
+    sigma: float = 0.5                   # ranking sharpness
+    all_landmark_nll: bool = False       # NLL at every (patient, visit) landmark
+    landmark_row_weighting: str = "patient"   # patient (1/V_i) | visit | band
+    landmark_band_edges: list[float] | None = None
+    landmark_decay_days: float | None = None
+    landmark_augment_prob: float = 0.0
+    gamma: float = 0.0                   # horizon-specific likelihood term
+    horizon_days: list[int] | None = None
+    horizon_ipcw_clip: float = 10.0
+    first_visit_only: bool = False
+
+
+def set_determinism(seed: int) -> None:
+    os.environ["PYTHONHASHSEED"] = str(seed)
+    random.seed(seed)
+    np.random.seed(seed)
+    torch.manual_seed(seed)
+    if torch.cuda.is_available():
+        torch.cuda.manual_seed_all(seed)
+    torch.use_deterministic_algorithms(True, warn_only=True)
+    if hasattr(torch.backends, "cudnn"):
+        torch.backends.cudnn.deterministic = True
+        torch.backends.cudnn.benchmark = False
+
+
+def categorical_embedding(train: pd.DataFrame, features: list[str], cols: Columns) -> tuple[int | None, int | None]:
+    """(vocabulary size, column index) for the embedded column, or (None, None)."""
+    if cols.categorical is None or cols.categorical not in features:
+        return None, None
+    return int(train[cols.categorical].max()) + 1, features.index(cols.categorical)
+
+
+def build_model(split_time: list[float], arch: Architecture, objective: Objective,
+                vocab_size: int | None, categorical_index: int | None) -> DynamicDeepHit:
+    embed = vocab_size is not None
+    return DynamicDeepHit(
+        split=split_time, layers_rnn=arch.layers_rnn, hidden_rnn=arch.hidden_rnn, typ=arch.typ,
+        long_param={"layers": arch.long_layers, "dropout": arch.dropout},
+        att_param={"layers": arch.att_layers, "dropout": arch.dropout},
+        cs_param={"layers": arch.cs_layers, "dropout": arch.dropout},
+        alpha=objective.alpha, beta=objective.beta, sigma=objective.sigma, cuda=False,
+        cancer_vocab_size=vocab_size,
+        cancer_embedding_dim=arch.embedding_dim if embed else None,
+        cancer_id_col=categorical_index,
+        hazard=arch.hazard, proportional=arch.proportional,
+    )
+
+
+def train(train: SequenceDataset, val: SequenceDataset, horizons: list[int],
+          arch: Architecture, objective: Objective, vocab_size: int | None,
+          categorical_index: int | None, seed: int = 42, iters: int = 64,
+          batch_size: int = 256, learning_rate: float = 1e-3, patience_limit: int = 10,
+          grad_clip_norm: float = 1.0) -> DynamicDeepHit:
+    """Fit on one split. The likelihood is conditioned on each patient's landmark s."""
+    set_determinism(seed)
+    split_time = choose_time_grid(train, val, horizons)
+    model = build_model(split_time, arch, objective, vocab_size, categorical_index)
+    needs_days = objective.all_landmark_nll or objective.landmark_augment_prob > 0
+    model.fit(
+        train.x, train.t, train.e,
+        val_data=(val.x, val.t, val.e),
+        iters=iters, learning_rate=learning_rate, batch_size=batch_size, random_state=seed,
+        s=train.last_obs_days, s_val=val.last_obs_days,
+        patience_limit=patience_limit, grad_clip_norm=grad_clip_norm,
+        visit_days=train.visit_days if needs_days else None,
+        visit_days_val=val.visit_days if objective.all_landmark_nll else None,
+        landmark_augment_prob=objective.landmark_augment_prob,
+        all_landmark_nll=objective.all_landmark_nll,
+        landmark_decay_days=objective.landmark_decay_days,
+        landmark_row_weighting=objective.landmark_row_weighting,
+        landmark_band_edges=objective.landmark_band_edges,
+        gamma=objective.gamma,
+        horizon_days=list(objective.horizon_days) if objective.horizon_days else None,
+        horizon_ipcw_clip=objective.horizon_ipcw_clip,
+        event_of_interest=EVENT_OF_INTEREST - 1,
+    )
+    return model
+
+
+# ---------------------------------------------------------------- prediction
+def predict_cif_curves(model: DynamicDeepHit, ds: SequenceDataset, risk: int,
+                       bs: int = 256) -> tuple[np.ndarray, np.ndarray]:
+    """Unconditional CIF of cause `risk` on the model's time grid, one row per sample."""
+    grid = np.asarray(model.split_time[1:], dtype=float)
+    curves = []
+    with torch.no_grad():
+        for start in range(0, len(ds.x), bs):
+            batch = ds.x[start:start + bs]
+            if len(batch) == 0:
+                continue
+            _, outcomes = model.torch_model(model._prepocess_test_data(batch))
+            curves.append(torch.cumsum(outcomes[risk - 1], dim=1).cpu().numpy())
+    if not curves:
+        return np.empty((0, len(grid))), grid
+    return np.concatenate(curves, axis=0), grid
+
+
+def condition_cif_on_landmark(cif_1: np.ndarray, cif_2: np.ndarray, grid: np.ndarray,
+                              s: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
+    """F_k(t | T > s) = (F_k(t) - F_k(s)) / S(s) for both causes.
+
+    Samples with s before the first grid point keep the raw CIF (exact at s = 0;
+    the model has no resolution below the first point), and so do samples with s
+    at or past the last grid point, where the conditional CIF is undefined.
+    """
+    s_idx = np.clip(np.searchsorted(grid, s, side="right") - 1, 0, len(grid) - 1)
+    rows = np.arange(len(cif_1))
+    f1_s, f2_s = cif_1[rows, s_idx], cif_2[rows, s_idx]
+    surv_s = np.clip(1.0 - f1_s - f2_s, 1e-8, 1.0)
+    c1 = np.clip((cif_1 - f1_s[:, None]) / surv_s[:, None], 0.0, 1.0)
+    c2 = np.clip((cif_2 - f2_s[:, None]) / surv_s[:, None], 0.0, 1.0)
+    raw = (s < grid[0]) | (s >= grid[-1])
+    c1[raw], c2[raw] = cif_1[raw], cif_2[raw]
+    return c1, c2
+
+
+def predict_at_horizons(model: DynamicDeepHit, ds: SequenceDataset,
+                        horizons: list[int], cols: Columns) -> pd.DataFrame:
+    """Risk of each cause by s + h, conditioned on surviving to s, for every sample.
+
+    Adds risk_{h}d (event of interest), competing_risk_{h}d and survival_{h}d to
+    the sample index, with entry_time = s and time_ref = time remaining after s.
+    Samples whose s is at or past the model's grid get NaN risks and
+    beyond_model_horizon = True, so they cannot slip into a metric.
+    """
+    cif_1, grid = predict_cif_curves(model, ds, EVENT_OF_INTEREST)
+    cif_2, _ = predict_cif_curves(model, ds, COMPETING_EVENT)
+    cif_1, cif_2 = condition_cif_on_landmark(cif_1, cif_2, grid, ds.last_obs_days)
+    at = ds.last_obs_days[:, None] + np.asarray(horizons, dtype=float)[None, :]
+    r1 = np.vstack([np.interp(at[i], grid, cif_1[i], left=0.0, right=float(cif_1[i, -1]))
+                    for i in range(len(cif_1))]) if len(cif_1) else np.empty((0, len(horizons)))
+    r2 = np.vstack([np.interp(at[i], grid, cif_2[i], left=0.0, right=float(cif_2[i, -1]))
+                    for i in range(len(cif_2))]) if len(cif_2) else np.empty((0, len(horizons)))
+
+    out = ds.sample_index.copy()
+    out["entry_time"] = ds.last_obs_days.astype(int)
+    out["time_ref"] = (out[cols.time] - ds.last_obs_days).clip(lower=0)
+    for j, h in enumerate(horizons):
+        out[f"risk_{h}d"] = r1[:, j]
+        out[f"competing_risk_{h}d"] = r2[:, j]
+        out[f"survival_{h}d"] = np.clip(1.0 - r1[:, j] - r2[:, j], 0.0, 1.0)
+    beyond = ds.last_obs_days >= grid[-1]
+    out["beyond_model_horizon"] = beyond
+    if beyond.any():
+        risk_cols = [f"{p}_{h}d" for p in ("risk", "competing_risk", "survival") for h in horizons]
+        out.loc[beyond, risk_cols] = np.nan
+    return out
+
+
+# ---------------------------------------------------------------- persistence
+def file_fingerprint(paths: dict[str, Path]) -> dict[str, Any]:
+    """Content hash of each input file: two runs on different data cannot be compared silently."""
+    out = {}
+    for name, path in paths.items():
+        if not Path(path).exists():
+            continue
+        digest = hashlib.sha256()
+        with open(path, "rb") as fh:
+            for chunk in iter(lambda: fh.read(1 << 20), b""):
+                digest.update(chunk)
+        out[name] = {"file": Path(path).name, "sha256": digest.hexdigest()[:16],
+                     "bytes": Path(path).stat().st_size}
+    return out
+
+
+def environment() -> dict[str, Any]:
+    from importlib import metadata as im
+    versions = {}
+    for pkg in ("torch", "numpy", "pandas", "scikit-learn", "statsmodels", "scipy"):
+        try:
+            versions[pkg] = im.version(pkg)
+        except Exception:
+            versions[pkg] = None
+    repo = Path(__file__).resolve().parents[1]
+    try:
+        commit = subprocess.run(["git", "rev-parse", "HEAD"], cwd=repo, capture_output=True,
+                                text=True, timeout=10).stdout.strip() or None
+    except Exception:
+        commit = None
+    return {"python": sys.version.split()[0], "platform": platform.platform(),
+            "packages": versions, "code_commit": commit}
+
+
+def save_run(run_dir: Path, model: DynamicDeepHit, *, predictors: list[str], horizons: list[int],
+             arch: Architecture, objective: Objective, cols: Columns, vocab_size: int | None,
+             categorical_index: int | None, inputs: dict[str, Path], seed: int,
+             predictions: dict[str, pd.DataFrame]) -> None:
+    run_dir.mkdir(parents=True, exist_ok=True)
+    torch.save(model.torch_model.state_dict(), run_dir / CHECKPOINT_FILE)
+    for name, frame in predictions.items():
+        frame.to_csv(run_dir / name, index=False)
+    meta = {
+        "predictors": predictors, "horizons": horizons,
+        "split_time": [float(v) for v in model.split_time],
+        "architecture": asdict(arch), "objective": asdict(objective), "columns": cols.to_dict(),
+        "categorical_vocab_size": vocab_size, "categorical_index": categorical_index,
+        "seed": seed, "inputs": file_fingerprint(inputs), "environment": environment(),
+    }
+    (run_dir / METADATA_FILE).write_text(json.dumps(meta, indent=2, default=_json_default))
+
+
+def load_run(run_dir: Path, cols: Columns, train_frame: pd.DataFrame | None = None) -> tuple[DynamicDeepHit, dict]:
+    """Rebuild the architecture recorded in a run's metadata and load its weights.
+
+    Runs written by older pipelines carry the objective under
+    validations.training_objective and no architecture block; those fall back to
+    the defaults above, and the categorical vocabulary is recomputed from the
+    split's training table.
+    """
+    meta = json.loads((run_dir / METADATA_FILE).read_text())
+    predictors = meta["predictors"]
+    legacy = meta.get("validations", {}).get("training_objective", {})
+    arch = Architecture(**meta["architecture"]) if "architecture" in meta else Architecture(
+        hidden_rnn=int(legacy.get("hidden_rnn") or 32), hazard=bool(legacy.get("hazard")),
+        proportional=bool(legacy.get("proportional")))
+    objective = Objective(**meta["objective"]) if "objective" in meta else Objective()
+    if "categorical_vocab_size" in meta:
+        vocab, index = meta["categorical_vocab_size"], meta["categorical_index"]
+    elif train_frame is not None:
+        vocab, index = categorical_embedding(train_frame, predictors, cols)
+    else:
+        raise ValueError(f"{run_dir}: metadata has no categorical vocabulary; pass the split's train table")
+    model = build_model(meta["split_time"], arch, objective, vocab, index)
+    model.torch_model = model._gen_torch_model(len(predictors), optimizer="Adam", risks=2)
+    model.torch_model.load_state_dict(torch.load(run_dir / _checkpoint_name(run_dir), map_location="cpu"))
+    model.torch_model.eval()
+    model.fitted = True
+    return model, meta
+
+
+def _checkpoint_name(run_dir: Path) -> str:
+    # Older runs named the checkpoint ddh_model_state.pt.
+    return CHECKPOINT_FILE if (run_dir / CHECKPOINT_FILE).exists() else "ddh_model_state.pt"
+
+
+def _json_default(value: Any) -> Any:
+    if isinstance(value, np.integer):
+        return int(value)
+    if isinstance(value, np.floating):
+        return float(value)
+    if isinstance(value, (pd.Timestamp, Path)):
+        return str(value)
+    raise TypeError(f"Object of type {type(value).__name__} is not JSON serializable")
diff --git a/evaluation/train.py b/evaluation/train.py
new file mode 100644
index 0000000..3162a8e
--- /dev/null
+++ b/evaluation/train.py
@@ -0,0 +1,149 @@
+"""Train one arm (one training configuration) on each split and write its predictions.
+
+For every split k this reads <data-root>/split_k/{train,val,test} (and the
+fixed-landmark tables when they exist), trains on train with early stopping on
+val, and writes to <runs-root>/<arm>/split_k/:
+
+    model_state.pt, metadata.json             checkpoint and everything needed to rebuild it
+    {train,val,test}_last_visit_predictions   each patient scored at their last visit
+    {train,val,test}_every_visit_predictions  each patient scored at every visit
+    {train,test}_fixed_landmark_predictions   each at-risk patient scored at each landmark day
+
+The arm the project reports (hazard head, all-landmark likelihood, band weighting):
+
+    python -m evaluation.train --arm band --all-landmark-nll --alpha 0 --hazard \\
+        --landmark-row-weighting band
+"""
+from __future__ import annotations
+
+import argparse
+import json
+import traceback
+
+import pandas as pd
+
+from .config import (EVERY_VISIT_FILE, FIXED_LANDMARK_FILE, LAST_VISIT_FILE, Columns,
+                     add_path_arguments, paths_from_args)
+from .data import (build_every_visit, build_fixed_landmarks, build_sequences, predictor_columns,
+                   read_visits)
+from .model import Architecture, Objective, categorical_embedding, predict_at_horizons, save_run, train
+
+
+def parse_args(argv=None) -> argparse.Namespace:
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    add_path_arguments(ap, out=False)
+    ap.add_argument("--arm", required=True, help="Name of this configuration; its runs go to <runs-root>/<arm>.")
+    ap.add_argument("--splits", type=int, nargs="+", default=list(range(10)))
+    ap.add_argument("--iters", type=int, default=64, help="Maximum epochs (early stopping on val).")
+    ap.add_argument("--batch-size", type=int, default=256)
+    ap.add_argument("--learning-rate", type=float, default=1e-3)
+    ap.add_argument("--seed", type=int, default=42)
+    ap.add_argument("--horizons", type=int, nargs="+", default=[30, 90, 180],
+                    help="Days after the landmark at which risk is reported.")
+    ap.add_argument("--landmark-days", type=int, nargs="+", default=[0, 30, 60, 90, 120, 150, 180, 270],
+                    help="Fixed landmark days scored from the landmark tables.")
+
+    obj = ap.add_argument_group("training objective (defaults reproduce the original loss)")
+    obj.add_argument("--alpha", type=float, default=0.01, help="Ranking term weight; 0 removes it.")
+    obj.add_argument("--beta", type=float, default=0.75, help="Survival likelihood weight.")
+    obj.add_argument("--sigma", type=float, default=0.5, help="Ranking sharpness.")
+    obj.add_argument("--all-landmark-nll", action="store_true",
+                     help="Score the conditional likelihood at every (patient, visit) landmark.")
+    obj.add_argument("--landmark-row-weighting", choices=["patient", "visit", "band"], default="patient",
+                     help="With --all-landmark-nll: patient = 1/V_i per row; visit = equal rows; "
+                          "band = each patient once per landmark-day band.")
+    obj.add_argument("--landmark-band-edges", type=float, nargs="+", default=[0, 30, 90, 180, 365, 730],
+                     help="Band edges in days for --landmark-row-weighting band.")
+    obj.add_argument("--landmark-decay-days", type=float, default=None,
+                     help="With --all-landmark-nll and patient weighting: weight landmarks by exp(-day/T).")
+    obj.add_argument("--landmark-augment-prob", type=float, default=0.0,
+                     help="Per-epoch probability of truncating a patient to a random earlier visit.")
+    obj.add_argument("--gamma", type=float, default=0.0, help="Horizon-specific likelihood weight.")
+    obj.add_argument("--horizon-days", type=int, nargs="+", default=[90, 180],
+                     help="Horizons the --gamma term scores.")
+    obj.add_argument("--horizon-ipcw-clip", type=float, default=10.0)
+    obj.add_argument("--first-visit-only", action="store_true",
+                     help="Static model: keep each patient's first visit only, in every part.")
+
+    arch = ap.add_argument_group("architecture")
+    arch.add_argument("--hidden-rnn", type=int, default=32)
+    arch.add_argument("--hazard", action="store_true", help="Discrete-time cause-specific hazard head.")
+    arch.add_argument("--proportional", action="store_true", help="With --hazard: proportional factorisation.")
+    args = ap.parse_args(argv)
+
+    if args.landmark_row_weighting != "patient" and not args.all_landmark_nll:
+        ap.error("--landmark-row-weighting needs --all-landmark-nll")
+    if args.all_landmark_nll and args.landmark_augment_prob > 0:
+        ap.error("--all-landmark-nll already covers every landmark; drop --landmark-augment-prob")
+    if args.gamma and not args.all_landmark_nll:
+        ap.error("--gamma needs --all-landmark-nll")
+    return args
+
+
+def run_split(split: int, args: argparse.Namespace, cols: Columns) -> dict:
+    paths = paths_from_args(args)
+    parts = {part: paths.split_file(split, part) for part in ("train", "val", "test")}
+    frames = {part: read_visits(path, cols) for part, path in parts.items()}
+    if args.first_visit_only:
+        frames = {part: f.groupby(cols.id, sort=False).head(1).reset_index(drop=True)
+                  for part, f in frames.items()}
+
+    predictors = predictor_columns(frames["train"], cols)
+    vocab, cat_index = categorical_embedding(frames["train"], predictors, cols)
+    sequences = {part: build_sequences(f, predictors, cols) for part, f in frames.items()}
+
+    arch = Architecture(hidden_rnn=args.hidden_rnn, hazard=args.hazard, proportional=args.proportional)
+    objective = Objective(
+        alpha=args.alpha, beta=args.beta, sigma=args.sigma, all_landmark_nll=args.all_landmark_nll,
+        landmark_row_weighting=args.landmark_row_weighting,
+        landmark_band_edges=args.landmark_band_edges if args.landmark_row_weighting == "band" else None,
+        landmark_decay_days=args.landmark_decay_days, landmark_augment_prob=args.landmark_augment_prob,
+        gamma=args.gamma, horizon_days=args.horizon_days, horizon_ipcw_clip=args.horizon_ipcw_clip,
+        first_visit_only=args.first_visit_only)
+    model = train(sequences["train"], sequences["val"], args.horizons, arch, objective, vocab, cat_index,
+                  seed=args.seed, iters=args.iters, batch_size=args.batch_size,
+                  learning_rate=args.learning_rate)
+
+    predictions = {}
+    for part, ds in sequences.items():
+        predictions[LAST_VISIT_FILE.format(part=part)] = predict_at_horizons(model, ds, args.horizons, cols)
+    if not args.first_visit_only:
+        for part, f in frames.items():
+            ds = build_every_visit(f, predictors, cols)
+            predictions[EVERY_VISIT_FILE.format(part=part)] = predict_at_horizons(model, ds, args.horizons, cols)
+    inputs = dict(parts)
+    for part in ("train", "test"):
+        table_path = paths.split_file(split, f"{part}_landmark")
+        if not table_path.exists():
+            continue
+        inputs[f"{part}_landmark"] = table_path
+        ds = build_fixed_landmarks(pd.read_csv(table_path, low_memory=False), predictors, cols,
+                                   args.landmark_days, first_visit_only=args.first_visit_only)
+        predictions[FIXED_LANDMARK_FILE.format(part=part)] = predict_at_horizons(model, ds, args.horizons, cols)
+
+    run_dir = paths.run_dir(args.arm, split)
+    save_run(run_dir, model, predictors=predictors, horizons=args.horizons, arch=arch,
+             objective=objective, cols=cols, vocab_size=vocab, categorical_index=cat_index,
+             inputs=inputs, seed=args.seed, predictions=predictions)
+    print(f"split {split}: {len(predictors)} predictors, "
+          f"{len(sequences['train'].x)} train patients -> {run_dir}", flush=True)
+    return {"split": split, "patients": {p: len(s.x) for p, s in sequences.items()}}
+
+
+def main(argv=None) -> None:
+    args = parse_args(argv)
+    cols = Columns.load(args.columns)
+    done, failed = [], []
+    for split in args.splits:
+        try:
+            done.append(run_split(split, args, cols))
+        except Exception as exc:
+            print(f"[error] split {split}: {exc}\n{traceback.format_exc()}", flush=True)
+            failed.append(split)
+    print(json.dumps({"completed": done, "failed": failed}, indent=2))
+    if failed:
+        raise SystemExit(f"failed splits: {failed}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/pyproject.toml b/pyproject.toml
index 809769d..3828313 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -24,9 +24,11 @@ dependencies = [
 
 [project.optional-dependencies]
 test = ["pytest"]
+# evaluation/ (training and evaluation scripts)
+eval = ["statsmodels", "scipy"]
 
 [tool.setuptools]
-packages = ["ddh"]
+packages = ["ddh", "evaluation"]
 
 [tool.pytest.ini_options]
 testpaths = ["tests"]
diff --git a/tests/test_evaluation_metrics.py b/tests/test_evaluation_metrics.py
new file mode 100644
index 0000000..3123138
--- /dev/null
+++ b/tests/test_evaluation_metrics.py
@@ -0,0 +1,145 @@
+"""Metric and data-building helpers of the evaluation package, on hand-checkable inputs."""
+import numpy as np
+import pandas as pd
+import pytest
+
+pytest.importorskip("statsmodels")
+
+from evaluation.config import Columns
+from evaluation.data import build_fixed_landmarks, build_landmark_table
+from evaluation.metrics import (aalen_johansen, auc, benjamini_hochberg, bootstrap, evaluate,
+                                ipcw_weights)
+from evaluation.model import condition_cif_on_landmark
+
+
+def cohort(seed=0, n=400):
+    rng = np.random.default_rng(seed)
+    t1, t2, c = rng.exponential(300, n), rng.exponential(500, n), rng.uniform(50, 600, n)
+    t = np.ceil(np.minimum.reduce([t1, t2, c]))
+    e = np.where(t1 <= np.minimum(t2, c), 1, np.where(t2 <= c, 2, 0))
+    return t, e
+
+
+def test_aalen_johansen_is_one_minus_kaplan_meier_without_competing_events():
+    t, e = cohort()
+    e = np.where(e == 2, 0, e)            # competing events become censoring
+    order = np.argsort(t)
+    surv, h = 1.0, 180
+    for time in np.unique(t[order]):
+        if time > h:
+            break
+        at_risk, d = (t >= time).sum(), ((t == time) & (e == 1)).sum()
+        surv *= 1 - d / at_risk
+    assert aalen_johansen(t, e, h) == pytest.approx(1 - surv, abs=1e-12)
+
+
+def test_aalen_johansen_matches_scikit_survival():
+    sksurv = pytest.importorskip("sksurv.nonparametric")
+    t, e = cohort(1)
+    times, cif = sksurv.cumulative_incidence_competing_risks(e, t)[:2]
+    for h in (30, 90, 180, 365):
+        idx = np.searchsorted(times, h, side="right") - 1
+        assert aalen_johansen(t, e, h) == pytest.approx(float(cif[1][idx]), abs=1e-12)
+
+
+def test_auc_counts_cases_by_h_against_controls_beyond_h():
+    t = np.array([10, 20, 30, 200, 300, 400, 50, 60.])
+    e = np.array([1, 1, 1, 0, 0, 2, 0, 2])      # last two: censored / competing before h
+    p = np.array([.9, .8, .7, .1, .2, .3, .99, .99])
+    assert auc(t, e, p, 100) == 1.0            # rows 6-7 must not enter either group
+    assert auc(t, e, 1 - p, 100) == 0.0
+    assert auc(t, e, np.full(8, .5), 100) == 0.5
+    assert np.isnan(auc(t, e, p, 15))          # fewer than 3 cases
+
+
+def test_auc_weights_cases_by_inverse_probability_of_censoring():
+    # Reverse KM of censoring: one of five at risk is censored at day 20, so G = 0.8
+    # from then on; the cases at 30 and 40 weigh 1/0.8 = 1.25, the case at 10 weighs 1.
+    t = np.array([10, 20, 30, 40, 200, 300.])
+    e = np.array([1, 0, 1, 1, 0, 0])
+    p = np.array([.9, .4, .2, .6, .5, .3])     # case shares of controls outranked: 1, 0, 1
+    assert auc(t, e, p, 100) == pytest.approx((1 * 1 + 1.25 * 0 + 1.25 * 1) / 3.5)
+
+
+def test_auc_matches_timeroc():
+    # R timeROC 0.4, weighting = "marginal", AUC_1 for cause 1, on this exact cohort
+    # (continuous times, so no event falls on day h, where the two definitions differ).
+    expected = {(0, 90): 0.670443158985, (0, 180): 0.686828729985,
+                (1, 90): 0.650093708109, (1, 180): 0.727709021032}
+    for (seed, h), value in expected.items():
+        rng = np.random.default_rng(seed)
+        x = rng.normal(size=400)
+        t1, t2 = rng.exponential(300 * np.exp(-0.7 * x)), rng.exponential(500, 400)
+        c = rng.uniform(30, 600, 400) * np.exp(-0.3 * x)
+        t = np.minimum.reduce([t1, t2, c])
+        e = np.where(t1 <= np.minimum(t2, c), 1, np.where(t2 <= c, 2, 0))
+        p = np.round(1 / (1 + np.exp(-(x + rng.normal(0, .5, 400)))), 2)
+        assert auc(t, e, p, h) == pytest.approx(value, abs=1e-11)
+
+
+def test_ipcw_weights_without_censoring_are_one():
+    t = np.array([10, 20, 200, 300.])
+    e = np.array([1, 2, 1, 2])
+    w, y = ipcw_weights(t, e, 100)
+    assert np.allclose(w, 1.0)
+    assert y.tolist() == [1, 0, 0, 0]
+
+
+def test_ipa_is_zero_for_the_null_model_and_oe_is_observed_over_predicted():
+    t, e = cohort(2)
+    obs = aalen_johansen(t, e, 180)
+    r = evaluate(t, e, {"null": np.full(len(t), obs), "half": np.full(len(t), obs / 2)}, 180)
+    assert r["null"]["ipa"] == pytest.approx(0.0, abs=1e-12)
+    assert r["null"]["op"] == pytest.approx(1.0)
+    assert r["half"]["op"] == pytest.approx(2.0)
+
+
+def test_bootstrap_is_paired():
+    t, e = cohort(3)
+    p = np.random.default_rng(0).uniform(size=len(t))
+    b = bootstrap(t, e, {"a": p, "b": p.copy()}, 180, n_boot=50)
+    assert np.array_equal(b["a"]["auc"], b["b"]["auc"])
+
+
+def test_benjamini_hochberg():
+    q = benjamini_hochberg(np.array([0.01, 0.04, 0.03, 0.20]))
+    assert np.allclose(q, [0.04, 0.16 / 3, 0.16 / 3, 0.20])
+
+
+def test_conditioning_on_the_landmark():
+    grid = np.array([10., 20., 30., 40.])
+    f1 = np.array([[.1, .2, .3, .4]] * 3)
+    f2 = np.array([[.0, .1, .1, .2]] * 3)
+    c1, _ = condition_cif_on_landmark(f1, f2, grid, np.array([0., 20., 45.]))
+    assert np.allclose(c1[0], f1[0])                       # s = 0: unchanged
+    assert np.allclose(c1[1], np.clip((f1[1] - .2) / (1 - .2 - .1), 0, 1))
+    assert np.allclose(c1[2], f1[2])                       # past the grid: raw CIF kept
+
+
+def test_landmark_table_and_fixed_landmark_sequences():
+    cols = Columns(categorical=None)
+    visits = pd.DataFrame({
+        "patient_id": [1, 1, 1, 2],
+        "visit": [1, 2, 3, 1],
+        "visit_date": pd.to_datetime(["2020-01-01", "2020-01-21", "2020-03-01", "2020-01-01"]),
+        "time_to_event": [100.0, 100.0, 100.0, 40.0],
+        "event": [1, 1, 1, 0],
+        "lab": [1.0, 2.0, 3.0, 5.0],
+        "days_since_x": [5.0, 25.0, 65.0, 0.0],
+        "gap": [0.0, 20.0, 40.0, 0.0],
+    })
+    table = build_landmark_table(visits, cols, [0, 30, 90], reanchor_cols=["days_since_x"], gap_col="gap")
+    snap = table[table.is_landmark == 1].set_index(["patient_id", "landmark_day"])
+    # patient 2 is censored at day 40, so it has no day-90 snapshot
+    assert sorted(snap.index) == [(1, 30), (1, 90), (2, 30)]
+    assert snap.loc[(1, 30), "lab"] == 2.0                 # forward-filled from day 20
+    assert snap.loc[(1, 30), "gap"] == 10.0                # days since the last real visit
+    assert snap.loc[(1, 30), "days_since_x"] == 35.0       # 25 at day 20, advanced 10 days
+    assert snap.loc[(1, 90), "lab"] == 3.0
+
+    ds = build_fixed_landmarks(table, ["lab", "days_since_x", "gap"], cols)
+    idx = ds.sample_index.set_index(["patient_id", "landmark_day"])
+    lengths = dict(zip(idx.index, [len(x) for x in ds.x]))
+    assert lengths == {(1, 0): 1, (1, 30): 3, (1, 90): 4, (2, 0): 1, (2, 30): 2}
+    assert dict(zip(idx.index, ds.last_obs_days)) == {(1, 0): 0, (1, 30): 30, (1, 90): 90,
+                                                      (2, 0): 0, (2, 30): 30}
diff --git a/tests/test_evaluation_pipeline.py b/tests/test_evaluation_pipeline.py
new file mode 100644
index 0000000..c75d42a
--- /dev/null
+++ b/tests/test_evaluation_pipeline.py
@@ -0,0 +1,38 @@
+"""Every evaluation script, end to end, on a tiny synthetic cohort (no patient data)."""
+import pandas as pd
+import pytest
+
+pytest.importorskip("statsmodels")
+
+from evaluation import apply_external, evaluate, history_ablation, make_synthetic, train
+
+
+def test_synthetic_end_to_end(tmp_path, monkeypatch):
+    data, runs, out = tmp_path / "data", tmp_path / "runs", tmp_path / "out"
+    make_synthetic.main(["--out", str(data), "--splits", "2", "--patients", "160",
+                         "--external-patients", "60", "--seed", "1"])
+    monkeypatch.setenv("DDH_DATA_ROOT", str(data))
+    monkeypatch.setenv("DDH_RUNS_ROOT", str(runs))
+    monkeypatch.setenv("DDH_OUTPUT_DIR", str(out))
+    monkeypatch.setenv("DDH_COLUMNS", str(data / "columns.json"))
+
+    common = ["--splits", "0", "1", "--iters", "2", "--batch-size", "64", "--all-landmark-nll",
+              "--alpha", "0", "--hazard"]
+    train.main(["--arm", "patient", *common])
+    train.main(["--arm", "band", *common, "--landmark-row-weighting", "band"])
+    for arm in ("patient", "band"):
+        apply_external.main(["--arm", arm, "--splits", "0", "1"])
+    evaluate.main(["--arm", "band", "--arm", "1/V_i=patient", "--splits", "0", "1",
+                   "--horizons", "90", "--n-boot", "20"])
+    history_ablation.main(["--arm", "band", "--splits", "0", "1", "--horizons", "90",
+                           "--gap-col", "days_since_last_obs", "--n-boot", "20"])
+
+    landmarks = pd.read_csv(out / "landmarks_summary.csv")
+    assert set(landmarks.arm) == {"band", "1/V_i"}
+    assert landmarks["pred"].between(0, 1).all()
+    external = pd.read_csv(out / "external.csv")
+    assert {"auc_delta_vs_ref", "auc_delta_lo", "auc_delta_hi"} <= set(external.columns)
+    assert external.loc[external.arm == "band", "auc_delta_vs_ref"].isna().all()
+    assert set(pd.read_csv(out / "ablation_summary.csv").variant) == {"full", "masked", "reanchored"}
+    for split in (0, 1):
+        assert (runs / "band" / f"split_{split}" / "metadata.json").exists()