From 23f76e8e5f0d83b08e5f42fbb9c26666003ea40b Mon Sep 17 00:00:00 2001
From: dadudu13 <dulin.wang@hotmail.com>
Date: Mon, 5 Oct 2026 14:49:22 -0500
Subject: [PATCH 1/5] Add the evaluation package: configuration, sequences,
 model I/O, metrics

The training and evaluation code behind the project's results lived in a
private analysis repository, mixed with cohort-specific paths and column
names. This ports the reusable part into evaluation/, with every path and
column name a setting:

- config.py: Columns (names of the id, date, visit, time, event and
  optional categorical columns, loadable from JSON) and Paths (data root,
  runs root, output dir, file-name templates), each defaulting to an
  environment variable.
- data.py: the three views a per-visit table is scored under -- one sample
  per patient (last visit), one per visit, one per fixed landmark day -- and
  build_landmark_table, which adds forward-filled snapshot rows for
  patients still at risk on each landmark day.
- model.py: build / train / save / load a split's model with its full
  architecture and objective in metadata.json, and predictions conditioned
  on the landmark, F_k(s + h | T > s).
- metrics.py: Aalen-Johansen incidence, Graf IPCW weights, the IPCW
  cumulative/dynamic AUC with a competing event (cases by h weighted by
  1/G(T-), controls event-free after h), Brier/IPA, calibration slope and
  intercept, top-decile O:E, patient bootstrap, Wilcoxon, BH.

The functions are ports, not rewrites: on the original cohort they reproduce
that repository's saved predictions to within 1e-16 and its metric tables to
within floating-point rounding. tests/test_evaluation_metrics.py checks the
metrics on hand-computed cases, Aalen-Johansen against scikit-survival, and
the AUC against R timeROC (identical apart from timeROC excluding events
that fall exactly on day h).

The slope fit returns NaN instead of failing when every prediction is equal
(statsmodels drops a constant regressor), which a coarse comparator score
can produce.
---
 evaluation/__init__.py           |   9 +
 evaluation/config.py             | 158 ++++++++++++++++
 evaluation/data.py               | 268 +++++++++++++++++++++++++++
 evaluation/metrics.py            | 209 +++++++++++++++++++++
 evaluation/model.py              | 299 +++++++++++++++++++++++++++++++
 pyproject.toml                   |   4 +-
 tests/test_evaluation_metrics.py | 145 +++++++++++++++
 7 files changed, 1091 insertions(+), 1 deletion(-)
 create mode 100644 evaluation/__init__.py
 create mode 100644 evaluation/config.py
 create mode 100644 evaluation/data.py
 create mode 100644 evaluation/metrics.py
 create mode 100644 evaluation/model.py
 create mode 100644 tests/test_evaluation_metrics.py

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/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/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/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}

From 85ba6365a2e77e6aad3428a8b29780ee23a4c55e Mon Sep 17 00:00:00 2001
From: dadudu13 <dulin.wang@hotmail.com>
Date: Mon, 5 Oct 2026 14:49:22 -0500
Subject: [PATCH 2/5] Add scripts to train an arm per split and to score an
 external cohort

evaluation.train trains one configuration on every split and writes, per
split, the checkpoint, metadata and three prediction files: each patient
at their last visit, at every visit, and at every fixed landmark day they
are still at risk. Its flags mirror the library's: --all-landmark-nll,
--hazard/--proportional, --landmark-row-weighting, --gamma and the rest,
with defaults that reproduce the original loss.

evaluation.apply_external loads each split's model and scores an external
cohort (for a temporal validation: every visit up to each patient's anchor)
at the anchor and at every retained visit.

Retraining the reported configuration on one split of the original cohort
with this script gives identical weights and predictions to the original
run.
---
 evaluation/apply_external.py |  57 ++++++++++++++
 evaluation/train.py          | 149 +++++++++++++++++++++++++++++++++++
 2 files changed, 206 insertions(+)
 create mode 100644 evaluation/apply_external.py
 create mode 100644 evaluation/train.py

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/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()

From 4e4f73fe6cdcddaf934c80db8cceed60a22f59bc Mon Sep 17 00:00:00 2001
From: dadudu13 <dulin.wang@hotmail.com>
Date: Mon, 5 Oct 2026 14:49:35 -0500
Subject: [PATCH 3/5] Add the evaluation script: fixed landmarks, every visit,
 external cohort

evaluation.evaluate scores several arms against a reference arm, raw (no
recalibration), on three populations:

- landmarks: internal test patients still at risk at days 0/90/180/270;
  per-split metrics, mean and SD, and each arm's per-split difference from
  the reference (Wilcoxon signed-rank, BH within horizon and metric);
- visits: internal test patients at every visit, by days since the first
  visit;
- external: the external cohort at the anchor, split models averaged into
  one ensemble, overall and by history; patient bootstrap for each arm and
  for its paired difference from the reference.

Arms are scored on identical rows: a row any arm cannot score is dropped for
all of them. That is also how an existing clinical score is compared: its
risks are written in the same file layout and passed as one more arm.

On the original cohort this reproduces the project's raw-calibration tables
row for row, including the bootstrap intervals and p-values.
---
 evaluation/evaluate.py | 215 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 215 insertions(+)
 create mode 100644 evaluation/evaluate.py

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

From fc9d86346c7d52d8f6182bc6a2a9c54608aaafc0 Mon Sep 17 00:00:00 2001
From: dadudu13 <dulin.wang@hotmail.com>
Date: Mon, 5 Oct 2026 14:49:35 -0500
Subject: [PATCH 4/5] Add a within-patient ablation of the visit history

Whether the model uses earlier visits cannot be read off a comparison of
patients with long and short records: they are different patients.
evaluation.history_ablation re-scores the SAME patients with the SAME
models on a changed record:

  full        the record up to the anchor (must equal the saved predictions)
  masked      the anchor visit alone, landmark kept at the anchor day
  reanchored  the anchor visit alone, landmark 0, gap column reset as for a
              first visit (with --gap-col)

and reports each variant's metrics and the paired bootstrap contrasts,
plus how far the masked record moves the ranking (Spearman, share of
patients changing risk decile).

On the original cohort its predictions for all three variants are identical
to the project's own ablation.
---
 evaluation/history_ablation.py | 136 +++++++++++++++++++++++++++++++++
 1 file changed, 136 insertions(+)
 create mode 100644 evaluation/history_ablation.py

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

From 96668a8e51bd49fabc950365c83ffcf4f346b319 Mon Sep 17 00:00:00 2001
From: dadudu13 <dulin.wang@hotmail.com>
Date: Mon, 5 Oct 2026 14:50:06 -0500
Subject: [PATCH 5/5] Add a synthetic cohort and run every evaluation script in
 CI

evaluation.make_synthetic writes a cohort in the expected layout, with no
real patient involved, so the workflow can be run, read and tested without
access to the data it was built for. It keeps the two properties that
matter here: visits stop at the event, so patients with an early event have
few visits (what biases 1/V_i weighting), and a lab value rises before the
event of interest (what a dynamic model can use).

tests/test_evaluation_pipeline.py runs make_synthetic -> train (two arms) ->
apply_external -> evaluate -> history_ablation on a tiny cohort, and CI now
installs the eval extra so it runs on every push.
---
 .github/workflows/tests.yml       |   2 +-
 evaluation/make_synthetic.py      | 120 ++++++++++++++++++++++++++++++
 tests/test_evaluation_pipeline.py |  38 ++++++++++
 3 files changed, 159 insertions(+), 1 deletion(-)
 create mode 100644 evaluation/make_synthetic.py
 create mode 100644 tests/test_evaluation_pipeline.py

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/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/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()