# CLAUDE.md — Orientation for a future agent session

> **NGÔN NGỮ / LANGUAGE (user request 2026-09-28):** Trao đổi với người dùng **bằng tiếng Việt**
> trong các phiên làm việc trên dự án này (giải thích, tóm tắt, câu hỏi). Mã nguồn, tên biến,
> comment kỹ thuật và log giữ nguyên tiếng Anh như hiện tại.

Operational subseasonal-to-seasonal (S2S) **extreme-heat-day (EHD)** forecast for Vietnam.
Two consumers: a duty forecaster running it on a schedule (outlooks to 45 days) and a
peer-reviewed paper. Method must survive review **and** be runnable by someone who did not
build it.

> Global rule (from the top-level repo CLAUDE.md): only write inside this working directory
> or its subdirectories. All data are read-only symlinks into `/work/projects/G09/...`; never
> modify, move, or delete them.

## Repository map

```
QG_1/
  config/config.yaml     # single source of truth: paths, grid, index, lead blocks, CV
  data/                  # read-only symlinks (see docs/DATA_INVENTORY.md)
  src/ehd/               # importable package (IO, climatology, index, calibration, verify)
  scripts/               # one runnable stage per file, NN_name.py, config-driven
  tests/                 # leakage-rule test + index-computation test (pytest)
  docs/                  # DATA_INVENTORY.md, operational guide, results report
  results/               # diagnostics, figures, tables, scores (evidence pack)
  outputs/ -> /work/.../QG1   # large NetCDF products (gridded P and binary fields)
  logs/                  # run logs
  MEMORY.md              # one lesson per entry; decisions, confirmations, dead ends
```

## Environment

`conda activate ehd` (Python 3.11; xarray, netCDF4, sklearn, xgboost, lightgbm, torch-CPU,
cartopy, shapely, pyproj). Pinned spec: `environment.yml`. Missing libs are either added
there or worked around (CRPS computed directly; region zonal stats via the gridded `reg`
mask instead of geopandas).

## How to run each stage (config-driven, seeded, logged) — ACTUAL scripts

```
conda activate ehd     # /home/tanpv/.conda/envs/ehd, Python 3.11 (pinned: environment.yml)
python scripts/01_climatology_threshold_diagnostic.py   # Tx90p, EHD index, min-vs-max diag
python scripts/02_obs_reconciliation.py                 # gridded vs station, per region
python scripts/03_build_forecast_targets.py             # 6h mx2t6 -> daily Tx per init (all yrs)
python scripts/04_calibrate_singlefold.py --holdout 2022  # dev smoke-test (one fold)
python scripts/06_run_loyo.py [--folds YYYY ...]        # FULL 20-fold LOYO -> outputs/loyo_pred/
python scripts/08_verify.py                             # pooled CV scores, bootstrap, figures
python scripts/09_regional_products.py                  # per-cell FDR skill grid, maps, partition
```
Recalibration + per-region×block decision thresholds (was planned as 05/07) are folded INTO
`06_run_loyo.py`. `06` writes one NPZ of out-of-sample predictions per held-out year; `08`
and `09` consume those (predictions computed once). `06` is ~40 min; run it detached and let
`pred_<year>.npz` files accumulate (its stdout buffers under nohup).

## CURRENT STATUS — resume here (updated 2026-09-13)

Pipeline is COMPLETE end-to-end: CV verification, regional products, season/severity
stratification, the operational per-init producer, the independent validator, and all written
deliverables. All 20 LOYO folds done (`outputs/loyo_pred/pred_*.npz`), 20 forecast years
extracted (`outputs/forecast_daily_tx/`), 10/10 tests pass. **All tasks 1–10 done.**

DONE (tasks 1–10):
- Data inventory + reconciliation (`docs/DATA_INVENTORY.md`, `docs/OBS_RECONCILIATION.md`).
- Climatology/Tx90p/EHD index + threshold diagnostic; obs reconciliation (gridded=target,
  station=validation); forecast daily-Tx extraction; calibration families; isotonic
  recalibration; per-region×block decision thresholds; 20-fold LOYO; pooled scores + year
  bootstrap; per-cell FDR field significance; regional partition table + maps + skill NetCDF.
- Task 8 tail — season (DJF/MAM/JJA/SON) + base-rate-severity stratification of CV scores
  (`results/tables/cv_scores_by_season_block.csv`, `cv_scores_by_severity_block.csv`).
- Task 9 tail — operational producer `scripts/10_operational_product.py`: PRODUCTION model
  refit on ALL years (`outputs/production_model/`) + `produce` writes one CF-1.8 NetCDF per
  init `outputs/operational/ehd_product_<YYYYMMDD>.nc` with prob_ehd + ehd_flag (block,lat,lon).
- Task 10 — `docs/RESULTS_REPORT.md` (paper-ready + evidence pack), `docs/OPERATIONAL_GUIDE.md`
  (Vietnamese primary + English summary), and independent validator
  `scripts/11_independent_validator.py` → `results/tables/validator_headline_rederivation.csv`
  (re-derives headline block scores from raw NPZ in plain numpy; **PASS**).
- Headline result recorded in MEMORY.md L-10/L-11.

REMAINING: none blocking. Optional follow-ups only — e.g. run `produce` over the full set of
inits/years for a complete operational archive (currently a 3-init demo for 2022 exists), and
have a human review the two written docs.

Real-time note for the operational producer: the reforecast is CONTROL-ONLY (1 member); a
real-time run must feed the same single-member features. If a real-time ensemble is used,
reduce it to the control/mean before applying the model, and check the operational model
version matches the reforecast (document any mismatch). See MEMORY.md L-01.

Key numbers (do not re-derive from memory — read the CSVs): headline in
`results/tables/cv_scores_by_block.csv` + `cv_bootstrap_ci.csv`; regional in
`cv_scores_by_region_block.csv` + `regional_partition.csv`.

## REQUESTED ENHANCEMENTS (v2) — user requests 2026-09-14, IMPLEMENTED + RUN (2026-09-14)

**STATUS: DONE.** v2 (severity predictand + expanded zoo + two-level cascade) is implemented
and run end-to-end on SLURM (`submit_loyo_v2.slurm`): 20-fold LOYO →
`outputs/loyo_pred_v2/pred_v2_*.npz`, then stage-13 pooled 2003–2022 out-of-sample binary
verification → `results/tables/v2_cv_scores_by_model_block.csv` (126 rows = 7 families × 3
cutoffs × 6 blocks), `v2_best_model_by_block.csv`, `results/figures/v2_binary_skill_by_lead.png`.
Headline unchanged from v1: skilful to ~W3 (best-by-SEDI W1 hgb 0.678 → W6 hgb_sev 0.153); the
LV1 easy-negative cascade screen is ~neutral for skill (cutoff sweep moves scores in the 3rd
decimal). See MEMORY.md L-12.

**v2 Phase A DONE (MEMORY L-13):** EHD_ens ensemble (ens_mean + leakage-safe LOYO ens_stack;
`scripts/14_ensemble_v2.py`), per-cell binary skill maps per lead week
(`scripts/15_skill_maps_v2.py` → `outputs/v2_skill_maps.nc` + `results/figures/v2_skillmap_*.png`),
and a self-contained HTML report (`scripts/16_report_v2_html.py` → `results/v2_report.html`).
ens_stack best at short lead, ens_mean best at long lead.
**v2 Phase B SUBMITTED (MEMORY L-13):** CNN (fully-convolutional) + deep torch ANN added
(`src/ehd/gridmodels.py`, `scripts/12b_run_cnn_ann.py`, `submit_cnn_ann.slurm`, job 1401269);
on completion the SLURM script re-runs stages 13→14→15→16 to fold them into EHD_ens + maps + HTML.
**v2 Phase B DONE (2026-09-15):** job 1401269 finished 02:59; 20/20 cnn+ann_torch folds in
`outputs/loyo_pred_v2_cnn/pred_v2cnn_*.npz`; stages 13→16 re-ran → 9-member ensemble
(`v2_scores_with_ensemble.csv`), maps (`v2_skill_maps.nc`), `results/v2_report.html`. Headline:
`ens_stack` beats best single member at every lead (W1 SEDI 0.678→0.715, BSS 0.351→0.391), meaningful
through W3, fades by W4–W6 (confirms v1). **v2 now written into `docs/RESULTS_REPORT.md` §10 +
evidence pack §11.**
**Remaining v2 (optional):** none blocking — full operational `produce` archive + human doc review.

## PROMPT §6 (vote ensemble + regional event table) — DONE 2026-09-15

User added `## 6. Additional requirements` to `S2S_EHD_Vietnam_Opus_Prompt.md`. Implemented in
`scripts/17_vote_ensemble_regional.py` (reads existing leakage-safe LOYO npz only; ~1m40 on login
node). Vote pool = **23 members** (7 sklearn families × 3 cascade cutoffs + cnn_c1 + ann_torch_c1),
each votes yEHDE with its stored per-region×block threshold → `P(yEHDE)=m/23`; deterministic
`EHDE=1 if P>P0` with per-fold Peirce-optimal P0 (leakage-safe); regional §6.4 `P=Mj/Nj` per
(year,init,region) with per-week P_thres (+ fixed-0.5 ref). Outputs `results/tables/v6_*.csv`,
`results/figures/v6_*.png`, `outputs/v6_vote_maps.nc`; folded into `results/v2_report.html`
(§"Vote ensemble & regions") and `docs/RESULTS_REPORT.md` §11 + evidence §12. See MEMORY
[[qg1-vote-ensemble-regional]]. Headline: vote ≈ best single member (< ens_stack at short lead),
P0 ≈ majority vote; CNN best single member at long lead; region×week much more skilful than
per-cell (W1 Peirce 0.65/SEDI 0.80). Caveat: per-fold regional P_thres jumpy (W5→0.10) — prefer
the fixed-0.5 columns operationally.

The original request text is preserved below for provenance. These extend the
calibration/verification and do **not** change the v1 hard constraints (six fixed blocks,
`EHD = min(Tx90p,35)`, leakage rule).

1. **Predictand = severity, not Tx.** Use **`Sev = Tx − Threshold`** (Threshold = `min(Tx90p,35°C)`)
   as the predictand instead of raw Tx. By construction `Sev ≥ 0 ⇔ the day is an EHD`, so severity
   regression maps directly to P(EHD) and also carries a "how extreme" magnitude.
2. **Expand the model zoo.** Add more AI model families to forecast EHD/EHDE from `Sev` (v1 has
   only logistic + HGB). Cf. the sibling HW report zoo (8 tier-1 regression + 6 tier-2 classifiers,
   see `../GHICHU_BaoCao_3PhuongAn.md`).
3. **Two-level (cascade) classification.**
   - **LV1 — screen easy negatives:** train a model to forecast **nEHDE** (no EHD/EHDE event).
     Samples predicted **nearly certain nEHDE** (probability ≈ 1) are **removed** from the pool.
   - **LV2 — model the remainder:** build AI models on the surviving (harder) samples to forecast
     **yEHDE** (yes EHD/EHDE event).
   - Rationale: strip the abundant, trivially-classified no-event samples so LV2 learns the hard
     discrimination without being swamped by easy negatives.
4. **Whole-dataset verification 2003–2022** with the binary contingency-table skill scores:
   **POD, POFD, FAR, EDI, SEDI** (+ the rest of the binary inventory: CSI, ETS/HSS, Peirce,
   frequency bias). MUST use out-of-sample (LOYO-pooled) predictions to keep the leakage rule —
   "whole dataset" means the 20 pooled held-out years, not an in-sample refit.

RESOLVED (user 2026-09-14):
- **EHDE = ≥1 EHD in block** (v1 "primary target") — "try with 1 EHD first". Same event as v1,
  so v2 vs v1 stays directly comparable.
- **LV1 "nearly certain nEHDE" cutoff = tunable, sweep starting 0.90 and increasing**
  (test 0.90 → 0.95 → 0.99…). Report skill vs cutoff; keep the sweep out-of-sample (leakage rule).
- **Scope = QG_1** (may reuse ideas/code from `../HW*`). User delegated method design to the agent.

EXPERT DESIGN DECISIONS (agent, 2026-09-14) — how v2 is being built:
- **Predictand `sev`** = block-max observed severity `max_{d∈block}(Tx_obs_d − Thr_d)`,
  Thr = `min(Tx90p,35)`. `y = 1[sev ≥ 0]` = yEHDE (identical to v1 primary). Regression families
  learn `sev`; predicted `sev` → P(yEHDE) via isotonic map fit on TRAIN (pred_sev, y). Classifier
  families learn `y` directly. Forecast severity added as a feature (`f_sevmax`).
- **Model zoo** (all single-member MOS, LOYO-nested): classifiers logistic, hgb, rf, **mlp (ANN)**;
  Sev-regressors ridge_sev, hgb_sev, **mlp_sev (ANN)**; **cnn** (fully-convolutional per-cell,
  torch-CPU, block-aggregated forecast maps as channels) added as a follow-on family.
- **Two-level cascade:** LV1 = hgb classifier on all TRAIN → screen samples with
  P(yEHDE) ≤ 1−cutoff as certain-nEHDE (final prob = LV1 value); LV2 = zoo models trained on the
  surviving TRAIN samples, predict the survivors. Final = screened ? LV1 : LV2.
- **New code (v1 left intact):** `src/ehd/cascade.py`, extended `calibration.py`,
  `scripts/12_run_loyo_v2.py` (→ `outputs/loyo_pred_v2/`), `scripts/13_verify_v2.py`
  (pooled 2003–2022 out-of-sample POD/POFD/FAR/EDI/SEDI/CSI/ETS/HSS/Peirce/bias per block×model×cutoff).

## Hard constraints (do not change without re-running verification)

- **Six fixed, non-overlapping lead blocks** W1 0–10, W2 11–17, W3 18–24, W4 25–31,
  W5 32–38, W6 39–45 days since init. **No running/sliding windows** anywhere in
  aggregation or verification.
- **EHD index** (user-fixed): a day is an EHD iff `Tx >= min(Tx90p, 35°C)`. Implement
  `min()` exactly. The `max(Tx90p,35)` variant is produced **only** as a diagnostic, never
  substituted into the product.
- **Tx90p** = calendar-day 90th percentile of observed Tx over the baseline period, using a
  centered multi-day window (`percentile_window_days`). This climatology-construction
  window is **not** a forecast running window and must be described separately in the paper.
- **Targets:** primary = P(≥1 EHD in block); secondary = P(heatwave, ≥3 consecutive EHDs),
  with block-boundary straddling handled explicitly and documented.
- **7 climate regions** from the gridded `reg` mask (ids 1–7); shapefile only for the map.
- **Output:** gridded 1/0 binary EHD field per block per init (CF-compliant NetCDF) with the
  underlying probability field retained alongside.

## LEAKAGE RULE (load-bearing — check at every handoff)

Cross-validation is **leave-one-year-out** over the 20 reforecast years, with hyperparameter
selection **nested** inside each left-out fold. For the held-out year *y*, none of the
following may use data from *y*:
1. the observed climatology and the **Tx90p threshold**,
2. the calibration model parameters (and their hyperparameters),
3. the probability recalibration (isotonic/logistic), and
4. the decision probability threshold (fit per block per region).

A threshold or climatology fit on data used for verification is a **failed deliverable**.
Each stage writes which years it was fit on into its output metadata; `06_verify.py`
asserts the held-out year is absent from every fit-year set before scoring.

## Data realities that shape the method (full detail in docs/DATA_INVENTORY.md)

- **Control-only, no ensemble.** Calibration is single-member MOS / distributional
  regression: map the deterministic forecast (+ covariates: calendar harmonics, lead, grid
  location, elevation, and the other control S2S fields) to a predictive distribution, then
  to P(EHD). P(EHD) is **never** a raw ensemble exceedance fraction.
- **Forecast Tx already on the 0.25° obs grid.** No downscaling stage; upstream
  interpolation is a documented limitation.
- **Forecast Tx is 6-hourly `mx2t6`** → daily Tx = daily max of four 6-hourly values.
- **Obs calendar is 365_day** (no 29 Feb); drop 29 Feb from forecast-derived Tx when
  matching to obs/climatology. Match on (month, day).
- **Real-time vs reforecast:** operational runs may have a different member count/model
  version than the reforecast the calibration was trained on — handled explicitly in the
  operational stage and flagged in the report.

## Scoring (always vs BOTH references)

Score calibrated forecasts against **climatology** and the **raw uncalibrated forecast**.
Continuous: MAE/RMSE/bias/CRPS/CRPSS, rank histogram (n/a single-member — note it),
spread–error. Probabilistic: Brier/BSS (+decomposition), reliability, ROC/AUC, sharpness.
Binary: POD, POFD, FAR, EDI, SEDI, CSI, ETS, HSS, Peirce, frequency bias. Stratify by lead
block, region, season, severity. Significance: block bootstrap resampling **whole years**;
field significance across cells controlled with **FDR**.