diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index d9511d1f..05c5b767 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -642,6 +642,39 @@ jobs: fi echo "No pack artifacts tracked" + # apps/atlas-outcomes: the outcome-real place-decision pack (stdlib only). + # Lint is a HARD gate (the app started at 0 findings). The tests run on a + # synthetic fixture, with a floor so a green run that collected nothing fails. + # The rebuild step proves `curate` reproduces the committed PACK_MANIFEST.json + # and sample byte for byte from the committed sources/ snapshots, offline. + atlas-outcomes: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - name: Setup Python 3.11 + uses: actions/setup-python@v5 + with: + python-version: "3.11" + - name: Ruff lint — apps/atlas-outcomes (HARD gate, at 0) + run: | + pip install ruff==0.15.22 + ruff check apps/atlas-outcomes + - name: Unit tests (floor 14) + run: | + set -euo pipefail + python3 -m unittest discover -s apps/atlas-outcomes/tests -v 2>&1 | tee /tmp/atlas.txt + ran=$(grep -oE 'Ran [0-9]+ tests?' /tmp/atlas.txt | tail -1 | grep -oE '[0-9]+' || echo 0) + if [ "${ran}" -lt 14 ]; then + echo "FAIL: only ${ran} tests ran (floor 14)." + exit 1 + fi + grep -qE '^OK$' /tmp/atlas.txt + - name: Offline rebuild matches the committed manifest + run: | + set -euo pipefail + python3 apps/atlas-outcomes/run.py curate > /dev/null + git diff --exit-code -- apps/atlas-outcomes/PACK_MANIFEST.json apps/atlas-outcomes/sample + bluehenre-checks: # The Dottie site (apps/bluehenre) — zero-dep, bare-node. Enforces the # provenance guarantees: pure parsers, the Hub exporter parser, and that the diff --git a/HANDOFF.md b/HANDOFF.md index 0dfab5da..743474f7 100644 --- a/HANDOFF.md +++ b/HANDOFF.md @@ -18,8 +18,8 @@ before writing "current" anywhere in this file. ## 📌 Session continuation — 2026-09-23 (supersedes every block below) -**Stamped 2026-09-23 at HEAD `07a2146`** (branch `claude/decision-plane`, the -merge of `origin/main` after PR #58's squash into this branch). Measured with +**Stamped 2026-09-23 at HEAD `8266088`** (`main`, the squash merge of PR #60; +the earlier stamp named a pre-squash branch commit that is not on `main`). Measured with `git rev-parse --short HEAD` + `scripts/check_handoff_fresh.py --check`. **What landed: decision plane, phase 1** (`docs/ARCHITECTURE.md` is now the diff --git a/apps/atlas-outcomes/BASELINE.json b/apps/atlas-outcomes/BASELINE.json new file mode 100644 index 00000000..dab9519b --- /dev/null +++ b/apps/atlas-outcomes/BASELINE.json @@ -0,0 +1,139 @@ +{ + "features": [ + "pct_date", + "pct_all", + "log_flow_over_p90", + "above_p90_today", + "chg1", + "chg3", + "chg7", + "up_present", + "up_pct_date", + "up_chg1", + "up_chg3", + "warn_active", + "warn_gauge_30d", + "warn_wfo_7d", + "warn_wfo_1d", + "doy_sin", + "doy_cos" + ], + "holdout": { + "days": [ + "2025-09-22", + "2026-09-21" + ], + "rows": 20921 + }, + "questions": { + "flow_change": { + "climatology": { + "accuracy": 0.4348, + "brier_multiclass": 0.69199, + "mae_expected_level": 0.8194 + }, + "holdout_levels": { + "falls 5% to 20%": 5976, + "falls more than 20%": 1531, + "rises 5% to 20%": 2716, + "rises more than 20%": 2334, + "within 5% either way": 8364 + }, + "logistic": { + "accuracy": 0.4653, + "brier_multiclass": 0.66652, + "mae_expected_level": 0.7431 + }, + "persistence": { + "accuracy": 0.5003, + "brier_multiclass": 0.99947, + "mae_expected_level": 0.7828 + }, + "type": "score" + }, + "high_next": { + "climatology": { + "accuracy": 0.9522, + "auc": 0.7303, + "brier": 0.04984 + }, + "holdout_positives": 999, + "holdout_rate": 0.0478, + "logistic": { + "accuracy": 0.9799, + "auc": 0.9818, + "brier": 0.01568 + }, + "logistic_weights": { + "above_p90_today": 0.6524, + "chg1": 0.4213, + "chg3": -0.2192, + "chg7": -0.0909, + "doy_cos": 0.1636, + "doy_sin": 0.0821, + "intercept": -6.2591, + "log_flow_over_p90": 3.2662, + "pct_all": 0.683, + "pct_date": -0.2581, + "up_chg1": 0.2979, + "up_chg3": 0.0537, + "up_pct_date": 0.5385, + "up_present": -0.4769, + "warn_active": 0.1392, + "warn_gauge_30d": 0.0272, + "warn_wfo_1d": 0.4312, + "warn_wfo_7d": -0.0074 + }, + "persistence": { + "accuracy": 0.9772, + "auc": 0.8744, + "brier": 0.0228 + }, + "type": "noul" + }, + "warn_next": { + "climatology": { + "accuracy": 0.9909, + "auc": 0.6968, + "brier": 0.00894 + }, + "holdout_positives": 190, + "holdout_rate": 0.0091, + "logistic": { + "accuracy": 0.9909, + "auc": 0.8313, + "brier": 0.00887 + }, + "logistic_weights": { + "above_p90_today": 0.0041, + "chg1": 0.0008, + "chg3": -0.0016, + "chg7": -0.0235, + "doy_cos": -0.1225, + "doy_sin": -0.0773, + "intercept": -5.19, + "log_flow_over_p90": -0.1188, + "pct_all": 0.4422, + "pct_date": 0.0294, + "up_chg1": 0.0263, + "up_chg3": 0.0548, + "up_pct_date": 0.0673, + "up_present": -0.0432, + "warn_active": 0.0264, + "warn_gauge_30d": 0.1802, + "warn_wfo_1d": 0.3904, + "warn_wfo_7d": 0.1468 + }, + "persistence": { + "accuracy": 0.9833, + "auc": 0.5326, + "brier": 0.01673 + }, + "type": "noul" + } + }, + "train": { + "rows": 79491, + "weighted_rows": 226290.0 + } +} diff --git a/apps/atlas-outcomes/PACK_MANIFEST.json b/apps/atlas-outcomes/PACK_MANIFEST.json new file mode 100644 index 00000000..b91ffbc7 --- /dev/null +++ b/apps/atlas-outcomes/PACK_MANIFEST.json @@ -0,0 +1,175 @@ +{ + "balance": { + "holdout": { + "flow_change": { + "falls 5% to 20%": 5976, + "falls more than 20%": 1531, + "rises 5% to 20%": 2716, + "rises more than 20%": 2334, + "within 5% either way": 8364 + }, + "high_next": { + "natural_rate": 0.0478, + "positive_rate": 0.0478, + "positives": 999 + }, + "rows": 20921, + "warn_next": { + "natural_rate": 0.0091, + "positive_rate": 0.0091, + "positives": 190 + }, + "weighted_rows": 20921.0 + }, + "train": { + "flow_change": { + "falls 5% to 20%": 23807, + "falls more than 20%": 9639, + "rises 5% to 20%": 9284, + "rises more than 20%": 10584, + "within 5% either way": 26177 + }, + "high_next": { + "natural_rate": 0.1033, + "positive_rate": 0.294, + "positives": 23368 + }, + "rows": 79491, + "warn_next": { + "natural_rate": 0.0094, + "positive_rate": 0.0268, + "positives": 2131 + }, + "weighted_rows": 226290.0 + } + }, + "consent": { + "capture_training": true, + "champion": false, + "public_sources": true + }, + "downsampling": { + "applies_to": "train only", + "easy_keep": 0.1, + "easy_negatives_before": { + "holdout": 16641, + "train": 163136 + }, + "rule": "both noul labels false, pct_date < 90, pct_all < 90, no warning over the gauge in 30 days" + }, + "files": { + "holdout.jsonl": { + "bytes": 37903287, + "sha256": "1ff022452ca5ca36a595de6b855d0b225f707843adbcae6c2047966e30a5ef4a" + }, + "provenance.jsonl": { + "bytes": 43704421, + "sha256": "fb0004042da132f828aec6639965d57b29b962ca860c8738ff5e9ed9d154e68b" + }, + "train.jsonl": { + "bytes": 144326335, + "sha256": "bbd6a2dc5c03ae5340df26819a8b186fc8e67d3b6c362c7705c5d0942b79a0ef" + } + }, + "gauges": { + "harvest_end": "2026-09-22", + "in_pack": 60, + "offices": [ + "BGM", + "BOX", + "BTV", + "CRP", + "CTP", + "DVN", + "EAX", + "EWX", + "FFC", + "FWD", + "HGX", + "ILM", + "ILN", + "ILX", + "LMK", + "LSX", + "LWX", + "MHX", + "MTR", + "PAH", + "PBZ", + "PHI", + "RAH", + "RLX", + "SEW" + ], + "states": [ + "CA", + "CT", + "GA", + "IA", + "IL", + "IN", + "KY", + "MA", + "MD", + "MO", + "NC", + "NJ", + "OH", + "PA", + "TX", + "VA", + "VT", + "WA", + "WV" + ] + }, + "label_sources": { + "flow_change": "usgs-nwis-dv", + "high_next": "usgs-nwis-dv", + "warn_next": "iem-vtec-sbw" + }, + "pack": "atlas-outcomes-1", + "provenance": "outcome-real", + "questions": { + "flow_change": "score (5 levels)", + "high_next": "noul", + "warn_next": "noul" + }, + "rejected": { + "easy negative downsampled (train)": 146825, + "embargo (label window touches the holdout)": 58 + }, + "rows": { + "candidates": 247295, + "holdout": 20921, + "holdout_frac": 0.2084, + "total": 100412, + "train": 79491 + }, + "rules": [ + "labels are recorded futures (USGS daily values on t+1, IEM VTEC warning polygons issued after T): provenance tier outcome-real", + "features use only data dated <= t and warnings issued <= T (end of local standard day t); percentiles from days strictly before t", + "holdout is the latest 365 days; a one-day embargo keeps every train label window out of it", + "easy negatives are downsampled in train only; sample_weight in provenance.jsonl restores natural rates", + "may train candidates; eval on the time-split holdout; consent.champion=false: nothing auto-promotes" + ], + "schema": "jev-decision-schema-1.0.0", + "seed": 20260923, + "sources": { + "flows.jsonl.gz": "6c0e7305d69e1afab7756e4bf6fc3a30c8019c187884eefc4da7029d1a152f3e", + "gauges.jsonl.gz": "5de5a28331aac5ed67a3c728e2e1bc9f3aee7c4c143c52aa8277323f390d26ae", + "warnings.jsonl.gz": "ca5d4e7756e2cc78922f36e7df0bb651012d3bba06e692915b12e92b04ce8b93" + }, + "split": { + "embargo_days": 1, + "holdout_days": [ + "2025-09-22", + "2026-09-21" + ], + "kind": "time", + "train_days": [ + "2015-01-01", + "2025-09-20" + ] + } +} diff --git a/apps/atlas-outcomes/README.md b/apps/atlas-outcomes/README.md new file mode 100644 index 00000000..d2e3a70d --- /dev/null +++ b/apps/atlas-outcomes/README.md @@ -0,0 +1,208 @@ +# atlas-outcomes: System One place decisions labelled by what happened next + +A `jev-decision-schema-1.0.0` pack for place decisions on the eye.jcamd.com +Atlas. Each record is one USGS streamgage at the end of day t: its flow and how +that flow ranks against the gauge's own history, the upstream gauge, the +construct stack it sits in (county, HUC-8, NWS office, state, as the Atlas +strata rail names them), which layers are active there (the streamflow +condition class, a flood-type warning in effect), and recent NWS warning +activity. The three questions are answered by **recorded futures**: the next +day's flow and the NWS warnings actually issued afterwards. No teacher and no +template writes a label. + +Provenance tier **`outcome-real`** (factory/datasets.json): the rows may train +System One candidates and are evaluated on a time-split holdout. The pack is +`consent.champion=false`. Nothing here trains, promotes or serves. + +## At a glance (pack `atlas-outcomes-1`, harvest end 2026-09-22) + +| | | +|---|---| +| Gauges | 60 in 19 states, 25 NWS offices; 25 with a paired upstream gauge | +| Sources | 626,718 gauge-days of daily flow (1995-10 -> 2026-09-22); 37,383 warning polygons (FL 16,467, FA 7,112, FF 13,804), 2015 -> 2026-09-23 | +| Candidate gauge-days | 247,295 (2015-01-01 -> 2026-09-21) | +| Train | 79,491 rows, 2015-01-01 -> 2025-09-20 (226,290 after weights) | +| Holdout | 20,921 rows, 2025-09-22 -> 2026-09-21, natural distribution | + +Class balance (raw rate in the pack / natural rate after `sample_weight`): + +| Question | Train | Holdout | +|---|---|---| +| `high_next` | 23,368 positives, 29.4% / 10.3% | 999 positives, 4.8% | +| `warn_next` | 2,131 positives, 2.7% / 0.94% | 190 positives, 0.91% | +| `flow_change` (falls >20 / 5-20 / within 5 / rises 5-20 / >20) | 9,639 / 23,807 / 26,177 / 9,284 / 10,584 | 1,531 / 5,976 / 8,364 / 2,716 / 2,334 | + +The holdout year was drier than the train decade (4.8% vs 10.3% of days +above p90), which is part of what a time split is for. + +Baselines on the holdout (`BASELINE.json`): + +| Question | Metric | Climatology | Persistence | Logistic (stdlib) | +|---|---|---|---|---| +| `high_next` | AUC / accuracy / Brier | 0.730 / 0.952 / 0.0498 | 0.874 / 0.977 / 0.0228 | **0.982** / 0.980 / **0.0157** | +| `warn_next` | AUC / accuracy / Brier | 0.697 / 0.991 / 0.0089 | 0.533 / 0.983 / 0.0167 | **0.831** / 0.991 / **0.0089** | +| `flow_change` | accuracy / multi-class Brier / MAE of E[level] | 0.435 / 0.692 / 0.819 | **0.500** / 0.999 / 0.783 | 0.465 / **0.667** / **0.743** | + + +## Stages + +```bash +python3 apps/atlas-outcomes/run.py harvest [--end YYYY-MM-DD] [--refresh] # network: USGS + IEM -> sources/*.jsonl.gz +python3 apps/atlas-outcomes/run.py features # every candidate gauge-day -> data/features/ (inspection only) +python3 apps/atlas-outcomes/run.py curate # offline: sources/ -> data/packs/atlas-outcomes-1/ + PACK_MANIFEST.json + sample/ +python3 apps/atlas-outcomes/run.py baseline # offline: CPU baselines on the holdout -> BASELINE.json +python3 apps/atlas-outcomes/run.py baseline --score preds.jsonl # a candidate's holdout answers on the same metrics +python3 apps/atlas-outcomes/run.py train-command # prints the GPU-host jev-v0 command; runs nothing +``` + +Stdlib only; the network goes through `curl`, throttled per host (1 s gap, +retries on 429/5xx). Raw responses are cached gzipped under `data/raw/` +(gitignored, like all of `data/`), so a rerun of `harvest` costs nothing. +`sources/` holds the committed, compact snapshots (gzip with a fixed mtime), and +`curate` rebuilds the pack from them offline, byte for byte (CI checks this). + +| Source | Endpoint | What it gives | +|---|---|---| +| USGS site service | `waterservices.usgs.gov/nwis/site` | name, lat/lon, county, HUC, drainage area, time zone | +| IEM UGC table | `mesonet.agron.iastate.edu/api/1/nws/ugcs.json` | county UGC -> NWS office (WFO) and county name | +| USGS daily values | `waterservices.usgs.gov/nwis/dv` (00060, statistic 00003) | mean daily discharge, 1995-10-01 -> `--end` | +| IEM VTEC polygons | `mesonet.agron.iastate.edu/api/1/vtec/sbw_interval.geojson` | flood (FL), areal flood (FA) and flash flood (FF) warning polygons at issuance, every gauge office, 2015 -> `--end` | + +`sources/gauges.jsonl.gz` (one row per gauge), `sources/flows.jsonl.gz` (per +gauge a start date and a dense daily cfs array, null where the service +reported nothing), `sources/warnings.jsonl.gz` (one row per warning polygon, +rings rounded to 1e-4 degrees) and `sources/harvest_stats.json`. + +The 60 gauges (`atlas_outcomes/gauges.py`) span 19 states and the NWS offices +that cover them, with upstream -> downstream pairs on the Potomac, +Susquehanna, Delaware, Ohio, Mississippi, Missouri, Iowa/Cedar, Trinity, +Brazos, Guadalupe, Neuse, Cape Fear, Russian, Snoqualmie/Snohomish and Skagit. + +## Features (what the state holds) + +The decision time of a row is **T = the end of local standard day t**, in +UTC: the moment day t's mean flow exists. Every feature is dated at or before +t, or is a warning issued at or before T (`atlas_outcomes/features.py`): + +| Feature | Definition | +|---|---| +| flow today | mean daily cfs on t | +| percentile for the date | rank of today's flow among the gauge's own flows within +-7 days of the same day of year, from days at least 8 days before t (so an event does not rank against itself); at least five years of such days | +| percentile of all prior days | rank among every day on record strictly before t | +| change vs 1/3/7 days ago | percent change | +| upstream | the paired upstream gauge's percentile for the date and 1/3-day change on day t | +| construct stack | County, HUC-8, NWS office, State (smallest first, the Atlas strata rail's kind names) | +| layers | `water`: the USGS WaterWatch class of today's flow (much below / below / normal / above / much above normal, the eye.jcamd.com streamflow condition); `nws-alert`: whether a flood-type warning polygon covering the gauge is in effect at T; `upstream water` where paired | +| warnings | polygons over the gauge issued in the last 30 days; polygons issued by the gauge's office in the last 7 days and 24 hours | + +## Questions and labels (recorded futures) + +| Question | Type | Label | +|---|---|---| +| `high_next` | noul | mean flow on t+1 above the threshold stated in the question: the gauge's 90th percentile over every day before t, to three significant figures | +| `warn_next` | noul | an NWS FL, FA or FF warning whose polygon contains the gauge is issued in (T, T + 24 h] | +| `flow_change` | score, 5 levels | t -> t+1 change: falls > 20%, falls 5-20%, within 5%, rises 5-20%, rises > 20% | + +The state never carries its own answer: `render_state` is a function of the +features alone (tests flip every label and the next-day flow and check the +state does not move). + +## Curation + +- **Time split.** Holdout = the latest 365 decision days; the day before it is + an embargo (its label window reaches into the holdout); train = everything + earlier. No gauge-day is in both, and every train label window closes before + the holdout starts. +- **Downsampled easy negatives (train only).** A row with both noul labels + false, flow under the 90th percentile for the date and of all days (not + "much above normal"), and no + warning over the gauge in 30 days is kept with probability 0.10 (a seeded + hash of site and date) and carries `sample_weight` 10 in `provenance.jsonl`, + so weighted counts give back the natural rates. The holdout keeps the + natural distribution so its metrics mean what they say. +- **Validate** every record with the frozen `apps/jev-v0/decision_io.validate_record`; + decontaminate with the dottie-os bench regex; drop duplicates and both sides + of any conflicting duplicate. Rejects are counted by reason in the manifest. +- **Shuffle** train with the seed: the jev-v0 trainer walks its file in order. + +`data/packs/atlas-outcomes-1/`: `train.jsonl`, `holdout.jsonl` (strict jev +records: `schema, id, state, questions, labels` only), `provenance.jsonl` (per +id: site, date, decision time, office, split, easy-negative flag, +sample_weight, provenance tier, label sources, consent) and `MANIFEST.json` +(copied here as `PACK_MANIFEST.json`: counts, class balance raw and weighted, +rejects, split dates, sha256 of every file and source snapshot). +`sample/pack_sample.jsonl` has a few rows per (split, label) combination. + +## Baselines + +`baseline` fits on train (weighted) and scores the holdout (unweighted): + +- **climatology**: the gauge's rate (or level distribution) for the calendar + month in train, shrunk towards the gauge's overall rate; +- **persistence**: today carried forward (flow above the threshold today; a + warning in effect over the gauge at T; today's 1-day change level); +- **logistic**: stdlib ridge logistic regression by IRLS on 17 numeric + features (one-vs-rest for `flow_change`). + +Noul: AUC, accuracy at 0.5, Brier. Score: accuracy of the most likely level, +multi-class Brier, mean absolute error of the expected level. The committed +`BASELINE.json` is the bar a System One candidate has to clear on the same +holdout. + +## Training (GPU host, not run here) + +`python3 apps/atlas-outcomes/run.py train-command` prints the exact commands. +Today it prints: + +```bash +# GPU host, from the repo root. Not run here. +# pack atlas-outcomes-1: 79491 train rows x 3 questions = 238473 steps (one pass) +pip install -r apps/jev-v0/requirements-jev-v0.txt # plus the host CUDA torch wheel +python3 apps/atlas-outcomes/run.py curate # rebuilds apps/atlas-outcomes/data/packs/atlas-outcomes-1/ offline from sources/ +python apps/jev-v0/train_pointer_lora.py --dry-run --fixtures apps/atlas-outcomes/data/packs/atlas-outcomes-1/train.jsonl +python apps/jev-v0/train_pointer_lora.py --go --fixtures apps/atlas-outcomes/data/packs/atlas-outcomes-1/train.jsonl --out apps/jev-v0/runs/atlas-outcomes-1 --steps 238473 +python apps/jev-v0/serve_decide.py --checkpoint apps/jev-v0/runs/atlas-outcomes-1 --port 8771 +``` + +`--steps` defaults to one pass over every question of every train record +(the trainer takes one question per step). Evaluate by serving the checkpoint +(`serve_decide.py --checkpoint`), answering `holdout.jsonl`, and scoring the +answers with `run.py baseline --score`. Promotion is the usual human stamp; +nothing here promotes. + +The jev-v0 trainer has no per-example weights: it sees the pack's raw rates +(29% / 2.7% positives), not the natural ones. Compare a candidate on the +holdout's AUC and Brier, where the baselines are, not on train loss. + +## Judgment calls and limits + +- **Revised values.** USGS daily values are fetched as they stand today + (approved or provisional), not as they were published at t. Revisions are + usually small, but they are information the decision maker at t did not + have exactly. +- **Standard time.** T uses the gauge's local standard time all year; daily + values follow local time, so during daylight time T is an hour late for + the warning windows. +- **Polygons, not counties.** "Covers the gauge" is point-in-polygon on the + warning's polygon at issuance. Extensions and polygon updates (CON/EXT) are + not harvested, so `warnings in effect at T` uses the issuance expiry. Only + the offices that contain a gauge are harvested; a neighbouring office's + polygon over a border gauge would be missed. +- **Connecticut** replaced its counties with planning regions in 2022; the NWS + UGC still uses the legacy county, so `01184000` maps to Hartford (`CTC003`). +- **Holdout at the natural rate.** Positives are rare in the holdout; read + AUC and Brier together, not accuracy alone. + +## Tests + +```bash +python3 -m unittest discover -s apps/atlas-outcomes/tests -v +``` + +Offline, on a synthetic two-gauge fixture built in the test (never written to +`sources/`): features ignore every flow and warning after t; p90 and +percentiles are strictly before t; the warning label window (T, T + 24 h] and +the polygon test; next-day labels match the stated threshold; the time split +has an embargo and no straddle; downsampling touches train only and sets the +weights; every record validates with jev-v0 and has exactly the five keys; +parsers and the baseline maths. diff --git a/apps/atlas-outcomes/atlas_outcomes/__init__.py b/apps/atlas-outcomes/atlas_outcomes/__init__.py new file mode 100644 index 00000000..86c4ba46 --- /dev/null +++ b/apps/atlas-outcomes/atlas_outcomes/__init__.py @@ -0,0 +1,12 @@ +"""atlas-outcomes: System One place decisions labelled by recorded futures. + +Stdlib only. Stages: harvest (USGS daily flow + IEM VTEC flood-warning +polygons) -> features (per gauge-day, only what was known at the end of day t) +-> curate (strict jev records, time-split holdout, manifest) -> baseline +(CPU logistic regression vs persistence and climatology). Labels are what +happened next (provenance tier outcome-real); see README.md. +""" + +SCHEMA_ID = "jev-decision-schema-1.0.0" +PACK_VERSION = "atlas-outcomes-1" +PROVENANCE_TIER = "outcome-real" diff --git a/apps/atlas-outcomes/atlas_outcomes/baseline.py b/apps/atlas-outcomes/atlas_outcomes/baseline.py new file mode 100644 index 00000000..9a965dd5 --- /dev/null +++ b/apps/atlas-outcomes/atlas_outcomes/baseline.py @@ -0,0 +1,307 @@ +"""Stage 4: CPU baselines on the time-split holdout. Stdlib only. + +Three predictors per question, fit on train (weighted by sample_weight, so the +downsampled easy negatives count at their natural rate) and scored on the +holdout (natural distribution, unweighted): + + climatology the gauge's rate for the calendar month in train, shrunk to the gauge rate + persistence today's state carried forward: high_next <- flow above the threshold + today; warn_next <- a warning in effect over the gauge at T; + flow_change <- today's 1-day change bucket + logistic ridge logistic regression (IRLS) on the numeric features + (flow_change: one-vs-rest, renormalised) + +Noul questions: AUC, accuracy at 0.5, Brier. Score question: accuracy of the +most likely level, multi-class Brier, and mean absolute error of the expected +level. The report goes to data/baseline/ and is copied to BASELINE.json here. +These are the numbers a System One candidate has to beat on the same holdout. +""" + +from __future__ import annotations + +import argparse +import json +import math +from collections import defaultdict +from operator import mul +from pathlib import Path +from typing import Any + +from . import PACK_VERSION +from .common import APP_ROOT, DATA, read_jsonl +from .curate import SEED, assemble +from .features import CHANGE_LEVELS, Corpus, bucket_of_change + +RIDGE = 1.0 +SHRINK = 20.0 +FEATURES = [ + "pct_date", "pct_all", "log_flow_over_p90", "above_p90_today", "chg1", "chg3", "chg7", + "up_present", "up_pct_date", "up_chg1", "up_chg3", + "warn_active", "warn_gauge_30d", "warn_wfo_7d", "warn_wfo_1d", "doy_sin", "doy_cos", +] + + +def _lchg(c: float | None) -> float: + if c is None: + return 0.0 + return max(-3.0, min(3.0, math.log(max(0.01, 1.0 + c / 100.0)))) + + +def vector(f: dict[str, Any]) -> list[float]: + ang = 2.0 * math.pi * f["doy"] / 365.0 + up = f.get("up_site") is not None + return [ + f["pct_date"] / 100.0, + (50.0 if f["pct_all"] is None else f["pct_all"]) / 100.0, + max(-6.0, min(6.0, math.log((f["flow"] + 1.0) / (f["p90"] + 1.0)))), + 1.0 if f["flow"] > f["p90"] else 0.0, + _lchg(f["chg1"]), _lchg(f["chg3"]), _lchg(f["chg7"]), + 1.0 if up else 0.0, + (f.get("up_pct_date") if up and f.get("up_pct_date") is not None else 50.0) / 100.0, + _lchg(f.get("up_chg1")), _lchg(f.get("up_chg3")), + 1.0 if f["warn_active"] else 0.0, + math.log1p(f["warn_gauge_30d"]), math.log1p(f["warn_wfo_7d"]), math.log1p(f["warn_wfo_1d"]), + math.sin(ang), math.cos(ang), + ] + + +def _sigmoid(z: float) -> float: + if z >= 0: + return 1.0 / (1.0 + math.exp(-z)) + e = math.exp(z) + return e / (1.0 + e) + + +def _solve(a: list[list[float]], b: list[float]) -> list[float]: + n = len(b) + m = [[*row, b[i]] for i, row in enumerate(a)] + for c in range(n): + p = max(range(c, n), key=lambda r: abs(m[r][c])) + m[c], m[p] = m[p], m[c] + if abs(m[c][c]) < 1e-12: + continue + for r in range(n): + if r != c and m[r][c]: + f = m[r][c] / m[c][c] + m[r] = [x - f * y for x, y in zip(m[r], m[c], strict=True)] + return [m[i][n] / m[i][i] if abs(m[i][i]) >= 1e-12 else 0.0 for i in range(n)] + + +class Standardizer: + def __init__(self, xs: list[list[float]], w: list[float]): + tw = sum(w) + d = len(xs[0]) + self.mu = [sum(wi * x[j] for wi, x in zip(w, xs, strict=True)) / tw for j in range(d)] + self.sd = [] + for j in range(d): + var = sum(wi * (x[j] - self.mu[j]) ** 2 for wi, x in zip(w, xs, strict=True)) / tw + self.sd.append(math.sqrt(var) or 1.0) + + def cols(self, xs: list[list[float]]) -> list[list[float]]: + """Column-major, standardized, with an intercept column first.""" + out = [[1.0] * len(xs)] + for j in range(len(self.mu)): + mu, sd = self.mu[j], self.sd[j] + out.append([(x[j] - mu) / sd for x in xs]) + return out + + +def fit_logistic(cols: list[list[float]], y: list[float], w: list[float], ridge: float = RIDGE, iters: int = 30) -> list[float]: + """Weighted ridge logistic regression by IRLS (Newton). The intercept is not penalised.""" + d, n = len(cols), len(y) + beta = [0.0] * d + for _ in range(iters): + eta = [0.0] * n + for j in range(d): + if beta[j]: + bj = beta[j] + eta = [e + bj * x for e, x in zip(eta, cols[j], strict=True)] + p = [_sigmoid(e) for e in eta] + r = [wi * (yi - pi) for wi, yi, pi in zip(w, y, p, strict=True)] + ww = [wi * pi * (1.0 - pi) for wi, pi in zip(w, p, strict=True)] + grad = [sum(map(mul, cols[j], r)) - (ridge * beta[j] if j else 0.0) for j in range(d)] + hess = [[0.0] * d for _ in range(d)] + for i in range(d): + wc = list(map(mul, ww, cols[i])) + for j in range(i, d): + hess[i][j] = hess[j][i] = sum(map(mul, wc, cols[j])) + if i: + hess[i][i] += ridge + step = _solve(hess, grad) + beta = [b + s for b, s in zip(beta, step, strict=True)] + if max(abs(s) for s in step) < 1e-7: + break + return beta + + +def predict(beta: list[float], cols: list[list[float]]) -> list[float]: + n = len(cols[0]) + eta = [0.0] * n + for j, bj in enumerate(beta): + eta = [e + bj * x for e, x in zip(eta, cols[j], strict=True)] + return [_sigmoid(e) for e in eta] + + +def auc(scores: list[float], y: list[int]) -> float | None: + pos = sum(y) + neg = len(y) - pos + if not pos or not neg: + return None + order = sorted(range(len(scores)), key=lambda i: scores[i]) + ranks = [0.0] * len(scores) + i = 0 + while i < len(order): + j = i + while j + 1 < len(order) and scores[order[j + 1]] == scores[order[i]]: + j += 1 + for k in range(i, j + 1): + ranks[order[k]] = (i + j) / 2.0 + 1.0 + i = j + 1 + u = sum(r for r, yi in zip(ranks, y, strict=True) if yi) - pos * (pos + 1) / 2.0 + return u / (pos * neg) + + +def noul_metrics(p: list[float], y: list[int]) -> dict[str, Any]: + a = auc(p, y) + return { + "auc": None if a is None else round(a, 4), + "accuracy": round(sum((pi >= 0.5) == bool(yi) for pi, yi in zip(p, y, strict=True)) / len(y), 4), + "brier": round(sum((pi - yi) ** 2 for pi, yi in zip(p, y, strict=True)) / len(y), 5), + } + + +def score_metrics(dists: list[list[float]], y: list[int]) -> dict[str, Any]: + k = len(CHANGE_LEVELS) + acc = sum(max(range(k), key=lambda c: d[c]) == yi for d, yi in zip(dists, y, strict=True)) / len(y) + brier = sum(sum((d[c] - (1.0 if c == yi else 0.0)) ** 2 for c in range(k)) for d, yi in zip(dists, y, strict=True)) / len(y) + mae = sum(abs(sum(c * d[c] for c in range(k)) - yi) for d, yi in zip(dists, y, strict=True)) / len(y) + return {"accuracy": round(acc, 4), "brier_multiclass": round(brier, 5), "mae_expected_level": round(mae, 4)} + + +def climatology(train: list[dict[str, Any]], key: str, classes: int | None = None) -> Any: + """(site, month) -> rate (noul) or distribution (score), shrunk to the site, then to the whole train.""" + k = classes or 2 + + def target(v: dict[str, Any]) -> int: + lab = v["row"]["labels"][key] + return int(lab) if classes is None else lab + + glob = [0.0] * k + site: dict[str, list[float]] = defaultdict(lambda: [0.0] * k) + cell: dict[tuple[str, int], list[float]] = defaultdict(lambda: [0.0] * k) + for v in train: + w, c = v["meta"]["sample_weight"], target(v) + glob[c] += w + site[v["meta"]["site"]][c] += w + cell[(v["meta"]["site"], v["row"]["features"]["month"])][c] += w + gt = sum(glob) + gd = [x / gt for x in glob] + + def shrink(counts: list[float], prior: list[float]) -> list[float]: + t = sum(counts) + return [(counts[c] + SHRINK * prior[c]) / (t + SHRINK) for c in range(k)] + + def dist(s: str, month: int) -> list[float]: + sd = shrink(site[s], gd) if s in site else gd + return shrink(cell[(s, month)], sd) if (s, month) in cell else sd + + return dist + + +def run(seed: int = SEED) -> dict[str, Any]: + kept, info = assemble(Corpus.load(), seed) + train = [v for v in kept if v["meta"]["split"] == "train"] + hold = [v for v in kept if v["meta"]["split"] == "holdout"] + xtr = [vector(v["row"]["features"]) for v in train] + xho = [vector(v["row"]["features"]) for v in hold] + w = [v["meta"]["sample_weight"] for v in train] + st = Standardizer(xtr, w) + ctr, cho = st.cols(xtr), st.cols(xho) + report: dict[str, Any] = { + "holdout": {"rows": len(hold), "days": [info["first_holdout_day"], info["last_day"]]}, + "train": {"rows": len(train), "weighted_rows": round(sum(w), 1)}, + "features": FEATURES, + "questions": {}, + } + for key in ("high_next", "warn_next"): + ytr = [1.0 if v["row"]["labels"][key] else 0.0 for v in train] + yho = [1 if v["row"]["labels"][key] else 0 for v in hold] + clim = climatology(train, key) + p_clim = [clim(v["meta"]["site"], v["row"]["features"]["month"])[1] for v in hold] + if key == "high_next": + p_pers = [1.0 if v["row"]["features"]["flow"] > v["row"]["features"]["p90"] else 0.0 for v in hold] + else: + p_pers = [1.0 if v["row"]["features"]["warn_active"] else 0.0 for v in hold] + beta = fit_logistic(ctr, ytr, w) + p_lr = predict(beta, cho) + report["questions"][key] = { + "type": "noul", + "holdout_positives": sum(yho), + "holdout_rate": round(sum(yho) / len(yho), 4), + "climatology": noul_metrics(p_clim, yho), + "persistence": noul_metrics(p_pers, yho), + "logistic": noul_metrics(p_lr, yho), + "logistic_weights": dict(zip(["intercept", *FEATURES], [round(b, 4) for b in beta], strict=True)), + } + k = len(CHANGE_LEVELS) + yho = [v["row"]["labels"]["change_bucket"] for v in hold] + clim = climatology(train, "change_bucket", classes=k) + d_clim = [clim(v["meta"]["site"], v["row"]["features"]["month"]) for v in hold] + d_pers = [] + for v in hold: + f = v["row"]["features"] + b = 2 if f["chg1"] is None else bucket_of_change(f["chg1"]) + d_pers.append([1.0 if c == b else 0.0 for c in range(k)]) + per_class = [] + for c in range(k): + yc = [1.0 if v["row"]["labels"]["change_bucket"] == c else 0.0 for v in train] + per_class.append(predict(fit_logistic(ctr, yc, w), cho)) + d_lr = [] + for i in range(len(hold)): + s = sum(per_class[c][i] for c in range(k)) + d_lr.append([per_class[c][i] / s for c in range(k)]) + report["questions"]["flow_change"] = { + "type": "score", + "holdout_levels": {CHANGE_LEVELS[c]: yho.count(c) for c in range(k)}, + "climatology": score_metrics(d_clim, yho), + "persistence": score_metrics(d_pers, yho), + "logistic": score_metrics(d_lr, yho), + } + out = DATA / "baseline" + out.mkdir(parents=True, exist_ok=True) + text = json.dumps(report, indent=2, sort_keys=True) + "\n" + (out / "baseline_report.json").write_text(text, encoding="utf-8") + (APP_ROOT / "BASELINE.json").write_text(text, encoding="utf-8") + return report + + +def score_predictions(pred_path: Path, pack: str = PACK_VERSION) -> dict[str, Any]: + """A candidate's holdout answers, one JSON line {"id", "answers": {qid: /decide answer}}, on the baseline metrics.""" + hold = {r["id"]: r for r in read_jsonl(DATA / "packs" / pack / "holdout.jsonl")} + preds = {p["id"]: p["answers"] for p in read_jsonl(pred_path) if p.get("id") in hold} + if not preds: + raise SystemExit(f"no prediction ids match {pack}/holdout.jsonl") + ids = sorted(preds) + out: dict[str, Any] = {"pack": pack, "holdout_rows": len(hold), "scored_rows": len(ids)} + for key in ("high_next", "warn_next"): + y = [int(hold[i]["labels"][key]["noul"] >= 0.5) for i in ids] + out[key] = noul_metrics([float(preds[i][key]["noul"]) for i in ids], y) + k = len(CHANGE_LEVELS) + y = [round(hold[i]["labels"]["flow_change"]["score"]) for i in ids] + dists = [[float(preds[i]["flow_change"]["probabilities"][str(c)]) for c in range(k)] for i in ids] + out["flow_change"] = score_metrics(dists, y) + return out + + +def main(argv: list[str] | None = None) -> int: + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--seed", type=int, default=SEED) + ap.add_argument("--score", type=Path, default=None, help="score a candidate's holdout predictions (JSONL) instead of fitting baselines") + args = ap.parse_args(argv) + if args.score: + print(json.dumps(score_predictions(args.score), indent=2)) + return 0 + r = run(args.seed) + for q, m in r["questions"].items(): + print(q, json.dumps({k: m[k] for k in ("climatology", "persistence", "logistic")})) + return 0 diff --git a/apps/atlas-outcomes/atlas_outcomes/common.py b/apps/atlas-outcomes/atlas_outcomes/common.py new file mode 100644 index 00000000..567f3139 --- /dev/null +++ b/apps/atlas-outcomes/atlas_outcomes/common.py @@ -0,0 +1,137 @@ +"""Shared helpers: polite cached HTTP, JSONL (optionally gzip), stable ids. + +Same conventions as apps/arxiviq-factory/arxiviq_factory/common.py (curl for +the network, gzip with mtime=0 so identical content gives identical bytes). +""" + +from __future__ import annotations + +import gzip +import hashlib +import json +import shutil +import subprocess +import time +import urllib.error +import urllib.request +from pathlib import Path +from typing import TYPE_CHECKING, Any + +if TYPE_CHECKING: + from collections.abc import Iterable + +APP_ROOT = Path(__file__).resolve().parent.parent +SOURCES = APP_ROOT / "sources" +# data/ is gitignored repo-wide: raw responses and curated packs are regenerable. +DATA = APP_ROOT / "data" +RAW = DATA / "raw" +USER_AGENT = "atlas-outcomes/0.1 (+https://github.com/jcdavis131/dottie)" + +_last: dict[str, float] = {} + + +class HttpError(Exception): + def __init__(self, code: int, url: str): + super().__init__(f"HTTP {code} from {url[:160]}") + self.code = code + + +def _fetch(url: str, timeout: float) -> bytes: + if not shutil.which("curl"): + req = urllib.request.Request(url, headers={"User-Agent": USER_AGENT}) # noqa: S310 - https URLs built here + try: + with urllib.request.urlopen(req, timeout=timeout) as r: # noqa: S310 - https URLs built here + return r.read() + except urllib.error.HTTPError as e: + raise HttpError(e.code, url) from e + cmd = ["curl", "-sS", "-L", "--compressed", "--max-time", str(int(timeout)), "-A", USER_AGENT, + "-w", "\n%{http_code}", url] + proc = subprocess.run(cmd, capture_output=True, timeout=timeout + 15, check=False) + if proc.returncode != 0: + raise urllib.error.URLError(proc.stderr.decode("utf-8", "ignore")[:200]) + body, _, code = proc.stdout.rpartition(b"\n") + status = int(code or b"0") + if status >= 400: + raise HttpError(status, url) + return body + + +def polite_get(url: str, *, host_gap_s: float = 1.0, timeout: float = 120.0, tries: int = 4) -> bytes: + """GET with a per-host minimum gap and retries on 429/5xx and transport errors.""" + host = url.split("/")[2] + wait = 3.0 + for attempt in range(1, tries + 1): + gap = host_gap_s - (time.monotonic() - _last.get(host, 0.0)) + if gap > 0: + time.sleep(gap) + _last[host] = time.monotonic() + try: + return _fetch(url, timeout) + except HttpError as e: + if e.code in (429, 500, 502, 503, 504) and attempt < tries: + time.sleep(wait) + wait *= 2 + continue + raise + except (urllib.error.URLError, TimeoutError, subprocess.TimeoutExpired): + if attempt < tries: + time.sleep(wait) + wait *= 2 + continue + raise + raise RuntimeError("unreachable") + + +def cached_get(url: str, cache_name: str, *, refresh: bool = False, host_gap_s: float = 1.0) -> bytes: + """polite_get through a gzip cache under data/raw (gitignored), so a rerun costs no requests.""" + path = RAW / f"{cache_name}.gz" + if path.exists() and not refresh: + return gzip.decompress(path.read_bytes()) + body = polite_get(url, host_gap_s=host_gap_s) + path.parent.mkdir(parents=True, exist_ok=True) + tmp = path.with_suffix(".tmp") + tmp.write_bytes(gzip.compress(body, mtime=0)) + tmp.replace(path) + return body + + +def read_jsonl(path: Path) -> list[dict[str, Any]]: + if not path.exists(): + return [] + opener = gzip.open if path.suffix == ".gz" else open + with opener(path, "rt", encoding="utf-8") as f: + return [json.loads(line) for line in f if line.strip()] + + +def write_jsonl(path: Path, rows: Iterable[dict[str, Any]]) -> int: + path.parent.mkdir(parents=True, exist_ok=True) + n = 0 + if path.suffix == ".gz": + with path.open("wb") as raw, gzip.GzipFile(fileobj=raw, mode="wb", mtime=0) as gz: + for row in rows: + gz.write((json.dumps(row, ensure_ascii=False, sort_keys=True, separators=(",", ":")) + "\n").encode("utf-8")) + n += 1 + return n + with path.open("w", encoding="utf-8") as f: + for row in rows: + f.write(json.dumps(row, ensure_ascii=False, sort_keys=True) + "\n") + n += 1 + return n + + +def sha256_file(path: Path) -> str: + h = hashlib.sha256() + with path.open("rb") as f: + for chunk in iter(lambda: f.read(1 << 20), b""): + h.update(chunk) + return h.hexdigest() + + +def stable_unit(*parts: str) -> float: + """A deterministic number in [0, 1) from the parts (seeded sampling without RNG state).""" + h = hashlib.sha1("|".join(parts).encode("utf-8"), usedforsecurity=False).hexdigest() + return int(h[:12], 16) / float(1 << 48) + + +def stable_id(*parts: str, n: int = 12) -> str: + return hashlib.sha1("|".join(parts).encode("utf-8"), usedforsecurity=False).hexdigest()[:n] diff --git a/apps/atlas-outcomes/atlas_outcomes/curate.py b/apps/atlas-outcomes/atlas_outcomes/curate.py new file mode 100644 index 00000000..3ac90f56 --- /dev/null +++ b/apps/atlas-outcomes/atlas_outcomes/curate.py @@ -0,0 +1,243 @@ +"""Stage 3: curate the gauge-days into a training-ready pack. + + 1. build every (gauge, day) row from the committed snapshots (features.py) + 2. time split: the holdout is the latest HOLDOUT_DAYS days of decision + dates; the day before it is an embargo (its label window touches the + holdout), everything earlier is train. No gauge-day is in both. + 3. downsample the easy negative mass in TRAIN only: a row with both noul + labels false, flow below the 90th percentile for the date and of all + days, and no warning over the gauge in the last 30 days is kept with + probability EASY_KEEP (a seeded hash of site and date). Kept easy rows + carry sample_weight 1 / EASY_KEEP in provenance.jsonl, so weighted counts + recover the natural rates. The holdout keeps the natural distribution. + 4. render strict jev records, decontaminate with the dottie-os bench regex, + validate with the frozen jev-v0 validator (rejects counted, never repaired), + drop duplicates and both sides of any conflicting duplicate + 5. write data/packs//{train,holdout,provenance}.jsonl and MANIFEST.json + (copied to PACK_MANIFEST.json here) and sample/pack_sample.jsonl. + +Train rows are shuffled with the seed (the jev-v0 trainer walks its file in +order). Nothing here trains or promotes: consent.champion is false and the +labels are recorded futures (provenance tier outcome-real). +""" + +from __future__ import annotations + +import argparse +import hashlib +import importlib.util +import json +import random +import re +import sys +from collections import Counter, defaultdict +from datetime import date, timedelta +from typing import Any + +from . import PACK_VERSION, PROVENANCE_TIER, SCHEMA_ID +from .common import APP_ROOT, DATA, SOURCES, sha256_file, stable_unit, write_jsonl +from .features import CHANGE_LEVELS, ROW_START, Corpus, build_rows +from .records import to_record + +REPO = APP_ROOT.parent.parent +JEV_IO = REPO / "apps" / "jev-v0" / "decision_io.py" +HOLDOUT_DAYS = 365 +EMBARGO_DAYS = 1 +EASY_KEEP = 0.10 +# Easy: flow not 'much above normal' (under the 90th percentile for the date and of all days). +EASY_BELOW_PCT = 90 +SEED = 20260923 +CONSENT = {"capture_training": True, "public_sources": True, "champion": False} +LABEL_SOURCES = {"high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw", "flow_change": "usgs-nwis-dv"} + +# Same bench-contamination guard as apps/dottie-os and apps/arxiviq-factory. +_DECONTAM = re.compile(r"(arxiviq|openjev|jevbench|nanojev|typesafe\s+teacher|jev-v0\s+champion)", re.I) + + +def load_validator() -> Any: + spec = importlib.util.spec_from_file_location("jev_decision_io", JEV_IO) + if spec is None or spec.loader is None: + raise SystemExit(f"cannot load the jev-v0 validator at {JEV_IO}") + mod = importlib.util.module_from_spec(spec) + sys.modules["jev_decision_io"] = mod + spec.loader.exec_module(mod) + if mod.SCHEMA_ID != SCHEMA_ID: + raise SystemExit(f"jev-v0 schema is {mod.SCHEMA_ID}, this factory writes {SCHEMA_ID}") + return mod + + +def split_dates(last_day: date, holdout_days: int = HOLDOUT_DAYS, embargo_days: int = EMBARGO_DAYS) -> tuple[date, date]: + """(first holdout day, last train day).""" + first_hold = last_day - timedelta(days=holdout_days - 1) + return first_hold, first_hold - timedelta(days=embargo_days + 1) + + +def assign_split(d: str, first_hold: date, last_train: date) -> str | None: + day = date.fromisoformat(d) + if day >= first_hold: + return "holdout" + if day <= last_train: + return "train" + return None + + +def is_easy_negative(row: dict[str, Any]) -> bool: + f, lab = row["features"], row["labels"] + return (not lab["high_next"] and not lab["warn_next"] and f["pct_date"] < EASY_BELOW_PCT + and (f["pct_all"] is None or f["pct_all"] < EASY_BELOW_PCT) and not f["warn_active"] and f["warn_gauge_30d"] == 0) + + +def _key(rec: dict[str, Any], with_label: bool) -> str: + body = {"state": rec["state"], "questions": rec["questions"]} + if with_label: + body["labels"] = rec["labels"] + return hashlib.sha1(json.dumps(body, sort_keys=True).encode(), usedforsecurity=False).hexdigest() + + +def assemble(corpus: Corpus, seed: int = SEED, easy_keep: float = EASY_KEEP, *, row_start: date = ROW_START, + holdout_days: int = HOLDOUT_DAYS) -> tuple[list[dict[str, Any]], dict[str, Any]]: + """Rows kept for the pack, each {"row", "record", "meta"}, plus curation counts.""" + io = load_validator() + gauges = {g["site"]: g for g in corpus.gauges} + rows = build_rows(corpus, row_start) + last_day = date.fromisoformat(max(r["date"] for r in rows)) + first_hold, last_train = split_dates(last_day, holdout_days) + rejects: Counter = Counter() + info: dict[str, Any] = {"candidate_rows": len(rows), "first_holdout_day": first_hold.isoformat(), + "last_train_day": last_train.isoformat(), "last_day": last_day.isoformat(), + "easy_negatives": Counter()} + valid: list[dict[str, Any]] = [] + for r in rows: + split = assign_split(r["date"], first_hold, last_train) + if split is None: + rejects["embargo (label window touches the holdout)"] += 1 + continue + easy = is_easy_negative(r) + weight = 1.0 + if easy: + info["easy_negatives"][split] += 1 + if split == "train": + if stable_unit(str(seed), r["site"], r["date"]) >= easy_keep: + rejects["easy negative downsampled (train)"] += 1 + continue + weight = round(1.0 / easy_keep, 4) + rec = to_record(gauges[r["site"]], r) + if _DECONTAM.search(json.dumps(rec["state"]) + json.dumps(rec["questions"])): + rejects["decontam"] += 1 + continue + try: + rec = io.validate_record(rec) + except io.SchemaError as err: + rejects[f"schema: {str(err)[:60]}"] += 1 + continue + g = gauges[r["site"]] + meta = { + "id": rec["id"], "site": r["site"], "date": r["date"], "decision_utc": r["decision_utc"], + "wfo": g["wfo"], "state": g["state"], "split": split, "easy_negative": easy, "sample_weight": weight, + "provenance": PROVENANCE_TIER, "label_sources": LABEL_SOURCES, + } + valid.append({"row": r, "record": rec, "meta": meta}) + + by_q: dict[str, set[str]] = defaultdict(set) + for v in valid: + by_q[_key(v["record"], False)].add(_key(v["record"], True)) + conflicted = {k for k, labs in by_q.items() if len(labs) > 1} + seen: set[str] = set() + kept: list[dict[str, Any]] = [] + for v in valid: + if _key(v["record"], False) in conflicted: + rejects["conflicting labels"] += 1 + continue + k = _key(v["record"], True) + if k in seen: + rejects["duplicate"] += 1 + continue + seen.add(k) + kept.append(v) + info["rejected"] = dict(rejects) + info["easy_negatives"] = dict(info["easy_negatives"]) + return kept, info + + +def balance(items: list[dict[str, Any]]) -> dict[str, Any]: + """Class balance per label: raw counts, and weighted rates (the natural distribution before downsampling).""" + if not items: + return {} + w = [v["meta"]["sample_weight"] for v in items] + tw = sum(w) + out: dict[str, Any] = {"rows": len(items), "weighted_rows": round(tw, 1)} + for lab in ("high_next", "warn_next"): + pos = [v["row"]["labels"][lab] for v in items] + out[lab] = {"positives": sum(pos), "positive_rate": round(sum(pos) / len(items), 4), + "natural_rate": round(sum(wi for wi, p in zip(w, pos, strict=True) if p) / tw, 4)} + buckets = Counter(v["row"]["labels"]["change_bucket"] for v in items) + out["flow_change"] = {CHANGE_LEVELS[k]: buckets.get(k, 0) for k in range(len(CHANGE_LEVELS))} + return out + + +def curate(seed: int = SEED, version: str = PACK_VERSION) -> dict[str, Any]: + corpus = Corpus.load() + kept, info = assemble(corpus, seed) + train = [v for v in kept if v["meta"]["split"] == "train"] + hold = [v for v in kept if v["meta"]["split"] == "holdout"] + random.Random(seed).shuffle(train) + if not train or not hold: + raise SystemExit("split left train or holdout empty") + + out = DATA / "packs" / version + write_jsonl(out / "train.jsonl", (v["record"] for v in train)) + write_jsonl(out / "holdout.jsonl", (v["record"] for v in hold)) + write_jsonl(out / "provenance.jsonl", ({**v["meta"], "consent": CONSENT} for v in kept)) + stats = json.loads((SOURCES / "harvest_stats.json").read_text(encoding="utf-8")) + manifest = { + "pack": version, + "schema": SCHEMA_ID, + "seed": seed, + "provenance": PROVENANCE_TIER, + "consent": CONSENT, + "rows": {"total": len(kept), "train": len(train), "holdout": len(hold), + "holdout_frac": round(len(hold) / len(kept), 4), "candidates": info["candidate_rows"]}, + "split": {"kind": "time", "train_days": [min(v["meta"]["date"] for v in train), info["last_train_day"]], + "embargo_days": EMBARGO_DAYS, "holdout_days": [info["first_holdout_day"], info["last_day"]]}, + "downsampling": {"easy_keep": EASY_KEEP, "applies_to": "train only", "easy_negatives_before": info["easy_negatives"], + "rule": f"both noul labels false, pct_date < {EASY_BELOW_PCT}, pct_all < {EASY_BELOW_PCT}, no warning over the gauge in 30 days"}, + "balance": {"train": balance(train), "holdout": balance(hold)}, + "gauges": {"in_pack": len({v["meta"]["site"] for v in kept}), "offices": sorted({v["meta"]["wfo"] for v in kept}), + "states": sorted({v["meta"]["state"] for v in kept}), "harvest_end": stats["end"]}, + "questions": {"high_next": "noul", "warn_next": "noul", "flow_change": "score (5 levels)"}, + "label_sources": LABEL_SOURCES, + "rejected": info["rejected"], + "files": {n: {"sha256": sha256_file(out / n), "bytes": (out / n).stat().st_size} for n in ("train.jsonl", "holdout.jsonl", "provenance.jsonl")}, + "sources": {p.name: sha256_file(p) for p in sorted(SOURCES.glob("*.jsonl.gz"))}, + "rules": [ + "labels are recorded futures (USGS daily values on t+1, IEM VTEC warning polygons issued after T): provenance tier outcome-real", + "features use only data dated <= t and warnings issued <= T (end of local standard day t); percentiles from days strictly before t", + "holdout is the latest 365 days; a one-day embargo keeps every train label window out of it", + "easy negatives are downsampled in train only; sample_weight in provenance.jsonl restores natural rates", + "may train candidates; eval on the time-split holdout; consent.champion=false: nothing auto-promotes", + ], + } + (out / "MANIFEST.json").write_text(json.dumps(manifest, indent=2, sort_keys=True) + "\n", encoding="utf-8") + (APP_ROOT / "PACK_MANIFEST.json").write_text(json.dumps(manifest, indent=2, sort_keys=True) + "\n", encoding="utf-8") + + # Reviewable sample: one of each (split, high_next, warn_next) combination that exists, then a few easy rows. + sample: list[dict[str, Any]] = [] + per: Counter = Counter() + for v in sorted(kept, key=lambda v: (v["meta"]["date"], v["meta"]["site"])): + lab = v["row"]["labels"] + k = (v["meta"]["split"], lab["high_next"], lab["warn_next"]) + if per[k] < 2: + per[k] += 1 + sample.append({"record": v["record"], "provenance": v["meta"]}) + write_jsonl(APP_ROOT / "sample" / "pack_sample.jsonl", sample) + return manifest + + +def main(argv: list[str] | None = None) -> int: + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--seed", type=int, default=SEED) + ap.add_argument("--version", default=PACK_VERSION) + args = ap.parse_args(argv) + m = curate(args.seed, args.version) + print(json.dumps({k: m[k] for k in ("rows", "split", "balance", "rejected")}, indent=2)) + return 0 diff --git a/apps/atlas-outcomes/atlas_outcomes/features.py b/apps/atlas-outcomes/atlas_outcomes/features.py new file mode 100644 index 00000000..42b00383 --- /dev/null +++ b/apps/atlas-outcomes/atlas_outcomes/features.py @@ -0,0 +1,330 @@ +"""Stage 2: features per (gauge, day t), from what was known at the end of day t. + +The decision time of a row is T = the end of local (standard) day t, in UTC: +the moment day t's mean flow exists. Everything in a row's features comes +from data dated at or before t, or from warnings issued at or before T: + + flow the day's mean flow (cfs) + pct_date percentile of that flow among the gauge's own flows on the + same day of year (+- 7 days) in prior years. Only days at + least DOY_LAG_DAYS before t enter, so the current event does + not rank against itself. Day-of-year aware, like USGS + WaterWatch and the eye.jcamd.com streamflow condition. + pct_all, p90 percentile of the flow among, and the 90th percentile of, + every prior day on record (strictly before t) + chg1/3/7 percent change against t-1, t-3, t-7 + up_* the upstream gauge's pct_date, chg1 and chg3 on day t, where a pair exists + warn_active an NWS flood-type warning polygon covering the gauge in effect at T + warn_gauge_30d warnings covering the gauge issued in (T - 30 d, T] + warn_wfo_7d/1d warnings (any polygon) issued by the gauge's NWS office in (T - 7 d, T] / (T - 24 h, T] + season, month calendar + +Labels are recorded futures, never inputs: + + high_next mean flow on t+1 > the threshold shown in the question (p90, 3 significant figures) + warn_next a flood-type warning polygon covering the gauge is issued in (T, T + 24 h] + change_bucket the change from t to t+1 in five buckets (CHANGE_LEVELS) + +`python3 run.py features` writes every candidate row (before curation) to +data/features/features.jsonl.gz with a stats file, for inspection. +""" + +from __future__ import annotations + +import argparse +import bisect +import json +import math +from collections import Counter +from dataclasses import dataclass, field +from datetime import UTC, date, datetime, timedelta +from typing import Any + +from .common import DATA, SOURCES, read_jsonl, write_jsonl +from .gauges import TZ_HOURS +from .geo import polygons_bbox, polygons_contain + +ROW_START = date(2015, 1, 1) +DOY_WINDOW = 7 +DOY_LAG_DAYS = 8 +# At least five prior years of same-season days before a row is emitted. +MIN_DOY_HISTORY = 5 * (2 * DOY_WINDOW + 1) +MIN_ALL_HISTORY = 5 * 365 +DAY_S = 86400 +ACTIVE_LOOKBACK_S = 20 * DAY_S +CHANGE_EDGES = (-20.0, -5.0, 5.0, 20.0) +CHANGE_LEVELS = [ + "falls more than 20%", + "falls 5% to 20%", + "within 5% either way", + "rises 5% to 20%", + "rises more than 20%", +] +SEASONS = {12: "winter", 1: "winter", 2: "winter", 3: "spring", 4: "spring", 5: "spring", + 6: "summer", 7: "summer", 8: "summer", 9: "fall", 10: "fall", 11: "fall"} + + +def doy_index(d: date) -> int: + """0..364; Feb 29 shares Feb 28's bin so every calendar day keeps its bin across years.""" + y = d.timetuple().tm_yday + leap = d.year % 4 == 0 and (d.year % 100 != 0 or d.year % 400 == 0) + if leap and y > 59: + y -= 1 + return y - 1 + + +def sig3(x: float) -> float: + """Round to three significant figures (the threshold as the question states it).""" + if x <= 0: + return 0.0 + digits = 2 - math.floor(math.log10(x)) + return float(round(x, digits)) + + +def pct_change(now: float | None, then: float | None) -> float | None: + if now is None or then is None: + return None + if then == 0: + return 0.0 if now == 0 else None + return round((now / then - 1.0) * 100.0, 1) + + +def change_bucket(today: float, tomorrow: float) -> int: + if today == 0: + return 2 if tomorrow == 0 else 4 + return bucket_of_change((tomorrow / today - 1.0) * 100.0) + + +def bucket_of_change(c: float) -> int: + """Percent change -> level index: < -20, [-20, -5), [-5, 5], (5, 20], > 20.""" + c = round(c, 6) # 95 / 100 - 1 is not exactly -5% in floating point + for i, edge in enumerate(CHANGE_EDGES): + if (c < edge) if i < 2 else (c <= edge): + return i + return 4 + + +def flow_class(pct: float | None) -> str: + """The USGS WaterWatch classes the eye.jcamd.com streamflow condition uses.""" + if pct is None: + return "not rated" + if pct < 10: + return "much below normal" + if pct < 25: + return "below normal" + if pct <= 75: + return "normal" + if pct <= 90: + return "above normal" + return "much above normal" + + +def _rank_pct(sorted_vals: list[float], x: float) -> tuple[int, int]: + lo = bisect.bisect_left(sorted_vals, x) + hi = bisect.bisect_right(sorted_vals, x) + return lo, hi - lo + + +@dataclass +class SiteHistory: + """Per-day percentiles for one gauge, each computed only from days before it.""" + + start: date + cfs: list[float | None] + pct_date: list[float | None] = field(default_factory=list) + pct_all: list[float | None] = field(default_factory=list) + p90: list[float | None] = field(default_factory=list) + + def index(self, d: date) -> int | None: + i = (d - self.start).days + return i if 0 <= i < len(self.cfs) else None + + def value(self, d: date) -> float | None: + i = self.index(d) + return None if i is None else self.cfs[i] + + +def site_history(start: date, cfs: list[float | None]) -> SiteHistory: + """Walk the record once. Day i's statistics see days < i (all-days) and days <= i - DOY_LAG_DAYS (day of year).""" + h = SiteHistory(start, cfs) + all_sorted: list[float] = [] + bins: list[list[float]] = [[] for _ in range(365)] + for i, x in enumerate(cfs): + d = start + timedelta(days=i) + lag = i - DOY_LAG_DAYS + if lag >= 0 and cfs[lag] is not None: + bisect.insort(bins[doy_index(start + timedelta(days=lag))], cfs[lag]) + pd = pa = q = None + if x is not None: + k = doy_index(d) + less = eq = n = 0 + for off in range(-DOY_WINDOW, DOY_WINDOW + 1): + b = bins[(k + off) % 365] + lo, e = _rank_pct(b, x) + less += lo + eq += e + n += len(b) + if n >= MIN_DOY_HISTORY: + pd = round(100.0 * (less + 0.5 * eq) / n, 1) + if len(all_sorted) >= MIN_ALL_HISTORY: + lo, e = _rank_pct(all_sorted, x) + pa = round(100.0 * (lo + 0.5 * e) / len(all_sorted), 1) + q = all_sorted[round(0.9 * (len(all_sorted) - 1))] + h.pct_date.append(pd) + h.pct_all.append(pa) + h.p90.append(q) + if x is not None: + bisect.insort(all_sorted, x) + return h + + +def iso_ts(s: str) -> int: + return int(datetime.strptime(s[:19], "%Y-%m-%dT%H:%M:%S").replace(tzinfo=UTC).timestamp()) + + +def decision_ts(d: date, tz: str) -> int: + """T: the end of local standard day d, as UTC epoch seconds.""" + midnight = datetime(d.year, d.month, d.day, tzinfo=UTC) + timedelta(days=1, hours=TZ_HOURS.get(tz, 5)) + return int(midnight.timestamp()) + + +@dataclass +class WarningIndex: + """Warnings covering each gauge (issue, expire) and issue times per office, all sorted.""" + + covering: dict[str, list[tuple[int, int]]] + by_wfo: dict[str, list[int]] + last_ts: int + + @classmethod + def build(cls, gauges: list[dict[str, Any]], warnings: list[dict[str, Any]], last_ts: int) -> WarningIndex: + covering: dict[str, list[tuple[int, int]]] = {g["site"]: [] for g in gauges} + by_wfo: dict[str, list[int]] = {} + for w in warnings: + t0, t1 = iso_ts(w["issue"]), iso_ts(w["expire"]) + by_wfo.setdefault(w["wfo"], []).append(t0) + x0, y0, x1, y1 = polygons_bbox(w["polygons"]) + for g in gauges: + if x0 <= g["lon"] <= x1 and y0 <= g["lat"] <= y1 and polygons_contain(w["polygons"], g["lon"], g["lat"]): + covering[g["site"]].append((t0, t1)) + for v in covering.values(): + v.sort() + for v in by_wfo.values(): + v.sort() + return cls(covering, by_wfo, last_ts) + + def features(self, site: str, wfo: str, t: int) -> dict[str, Any]: + cov = self.covering.get(site, []) + hi = bisect.bisect_right(cov, (t, 1 << 62)) + lo30 = bisect.bisect_right(cov, (t - 30 * DAY_S, 1 << 62)) + lo_act = bisect.bisect_right(cov, (t - ACTIVE_LOOKBACK_S, 1 << 62)) + issues = self.by_wfo.get(wfo, []) + end = bisect.bisect_right(issues, t) + return { + "warn_active": any(e > t for _, e in cov[lo_act:hi]), + "warn_gauge_30d": hi - lo30, + "warn_wfo_7d": end - bisect.bisect_right(issues, t - 7 * DAY_S), + "warn_wfo_1d": end - bisect.bisect_right(issues, t - DAY_S), + } + + def issued_after(self, site: str, t: int, horizon_s: int = DAY_S) -> bool | None: + """A covering warning issued in (t, t + horizon]; None when that window runs past the harvest.""" + if t + horizon_s > self.last_ts: + return None + cov = self.covering.get(site, []) + i = bisect.bisect_right(cov, (t, 1 << 62)) + return i < len(cov) and cov[i][0] <= t + horizon_s + + +@dataclass +class Corpus: + gauges: list[dict[str, Any]] + flows: list[dict[str, Any]] + warnings: list[dict[str, Any]] + warnings_through: str # last day (UTC, exclusive) the warning harvest covers + + @classmethod + def load(cls) -> Corpus: + stats = json.loads((SOURCES / "harvest_stats.json").read_text(encoding="utf-8")) + end = date.fromisoformat(stats["end"]) + timedelta(days=2) + return cls(read_jsonl(SOURCES / "gauges.jsonl.gz"), read_jsonl(SOURCES / "flows.jsonl.gz"), + read_jsonl(SOURCES / "warnings.jsonl.gz"), end.isoformat()) + + +def build_rows(corpus: Corpus, row_start: date = ROW_START) -> list[dict[str, Any]]: + """Every (gauge, day) with a full feature set and all three labels observed.""" + hist = {f["site"]: site_history(date.fromisoformat(f["start"]), f["cfs"]) for f in corpus.flows} + gauges = [g for g in corpus.gauges if g["site"] in hist] + last_ts = iso_ts(corpus.warnings_through + "T00:00:00") + widx = WarningIndex.build(gauges, corpus.warnings, last_ts) + names = {g["site"]: g for g in gauges} + rows: list[dict[str, Any]] = [] + for g in gauges: + h = hist[g["site"]] + up = hist.get(g["upstream"] or "") + for i in range(len(h.cfs) - 1): + d = h.start + timedelta(days=i) + if d < row_start: + continue + x, nxt = h.cfs[i], h.cfs[i + 1] + if x is None or nxt is None or h.pct_date[i] is None or h.p90[i] is None: + continue + t = decision_ts(d, g["tz"]) + warn_next = widx.issued_after(g["site"], t) + if warn_next is None: + continue + threshold = sig3(h.p90[i]) + feats: dict[str, Any] = { + "flow": x, + "pct_date": h.pct_date[i], + "pct_all": h.pct_all[i], + "p90": threshold, + "chg1": pct_change(x, h.cfs[i - 1] if i >= 1 else None), + "chg3": pct_change(x, h.cfs[i - 3] if i >= 3 else None), + "chg7": pct_change(x, h.cfs[i - 7] if i >= 7 else None), + "month": d.month, + "doy": doy_index(d), + "season": SEASONS[d.month], + **widx.features(g["site"], g["wfo"], t), + } + if up is not None: + j = up.index(d) + ux = up.value(d) + feats["up_site"] = g["upstream"] + feats["up_name"] = names.get(g["upstream"], {}).get("name", g["upstream"]) + feats["up_pct_date"] = up.pct_date[j] if j is not None else None + feats["up_chg1"] = pct_change(ux, up.value(d - timedelta(days=1))) + feats["up_chg3"] = pct_change(ux, up.value(d - timedelta(days=3))) + rows.append({ + "site": g["site"], + "date": d.isoformat(), + "decision_utc": datetime.fromtimestamp(t, UTC).strftime("%Y-%m-%dT%H:%MZ"), + "features": feats, + "labels": { + "high_next": nxt > threshold, + "warn_next": warn_next, + "change_bucket": change_bucket(x, nxt), + }, + "next_flow": nxt, + }) + return rows + + +def main(argv: list[str] | None = None) -> int: + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.parse_args(argv) + rows = build_rows(Corpus.load()) + out = DATA / "features" + write_jsonl(out / "features.jsonl.gz", rows) + stats = { + "rows": len(rows), + "gauges": len({r["site"] for r in rows}), + "first_day": min(r["date"] for r in rows), + "last_day": max(r["date"] for r in rows), + "high_next_rate": round(sum(r["labels"]["high_next"] for r in rows) / len(rows), 4), + "warn_next_rate": round(sum(r["labels"]["warn_next"] for r in rows) / len(rows), 4), + "change_buckets": dict(sorted(Counter(r["labels"]["change_bucket"] for r in rows).items())), + } + (out / "features_stats.json").write_text(json.dumps(stats, indent=2) + "\n", encoding="utf-8") + print(json.dumps(stats, indent=2)) + return 0 diff --git a/apps/atlas-outcomes/atlas_outcomes/gauges.py b/apps/atlas-outcomes/atlas_outcomes/gauges.py new file mode 100644 index 00000000..df38f9f0 --- /dev/null +++ b/apps/atlas-outcomes/atlas_outcomes/gauges.py @@ -0,0 +1,99 @@ +"""The gauges: 60 USGS streamgages across 19 states and ~20 NWS offices. + +Chosen for long daily-discharge records on rivers the NWS forecasts, with +upstream -> downstream pairs on the same river where one exists (the +upstream gauge's state at t is a feature of the downstream gauge's row). +Every site id's daily discharge was checked against the USGS site service on 2026-09-23; the +names, coordinates, county, HUC and drainage area come from that service at +harvest, not from this file. +""" + +from __future__ import annotations + +# site -> the upstream gauge on the same main stem (None when unpaired). +GAUGES: dict[str, str | None] = { + # Potomac (LWX) + "01610000": None, # Potomac at Paw Paw, WV + "01613000": "01610000", # Potomac at Hancock, MD + "01636500": None, # Shenandoah at Millville, WV + "01638500": "01613000", # Potomac at Point of Rocks, MD + "01643000": None, # Monocacy at Jug Bridge, MD + "01646500": "01638500", # Potomac at Little Falls, DC + "01668000": None, # Rappahannock near Fredericksburg, VA + # Susquehanna (BGM, CTP, LWX) + "01531500": None, # Towanda, PA + "01536500": "01531500", # Wilkes-Barre, PA + "01570500": "01536500", # Harrisburg, PA + "01578310": "01570500", # Conowingo, MD + # Delaware / Schuylkill / Choptank (BGM, PHI) + "01438500": None, # Delaware at Montague, NJ + "01446500": "01438500", # Delaware at Belvidere, NJ + "01463500": "01446500", # Delaware at Trenton, NJ + "01474500": None, # Schuylkill at Philadelphia, PA + "01491000": None, # Choptank near Greensboro, MD + # New England (BOX, BTV) + "01100000": None, # Merrimack at Lowell, MA + "01184000": None, # Connecticut at Thompsonville, CT + "04290500": None, # Winooski near Essex Junction, VT + # Carolinas / Georgia (RAH, MHX, ILM, FFC) + "02087500": None, # Neuse near Clayton, NC + "02089500": "02087500", # Neuse at Kinston, NC + "02102500": None, # Cape Fear at Lillington, NC + "02105769": "02102500", # Cape Fear at Lock 1, NC + "02336000": None, # Chattahoochee at Atlanta, GA + # Ohio River and tributaries (PBZ, ILN, LMK, PAH) + "03049500": None, # Allegheny at Natrona, PA + "03085000": None, # Monongahela at Braddock, PA + "03086000": "03049500", # Ohio at Sewickley, PA + "03150000": None, # Muskingum at McConnelsville, OH + "03234000": None, # Paint Creek near Bourneville, OH + "03234500": None, # Scioto at Higby, OH + "03274000": None, # Great Miami at Hamilton, OH + # Ohio at Cincinnati (03255000) publishes stage only; Markland Dam is the next discharge gauge up. + "03277200": None, # Ohio at Markland Dam near Warsaw, KY + "03294500": "03277200", # Ohio at Louisville, KY + "03303280": "03294500", # Ohio at Cannelton, IN + "03611500": "03303280", # Ohio at Metropolis, IL + # Upper Mississippi, Iowa, Illinois, Missouri (DVN, LSX, EAX) + "05420500": None, # Mississippi at Clinton, IA + "05474500": "05420500", # Mississippi at Keokuk, IA + "07010000": "05474500", # Mississippi at St. Louis, MO + "07022000": "07010000", # Mississippi at Thebes, IL + "05446500": None, # Rock River near Joslin, IL + "05464500": None, # Cedar at Cedar Rapids, IA + "05465500": "05464500", # Iowa River at Wapello, IA (below the Cedar confluence) + "05586100": None, # Illinois at Valley City, IL + "06893000": None, # Missouri at Kansas City, MO + "06934500": "06893000", # Missouri at Hermann, MO + # Texas (FWD, HGX, EWX) + "08057000": None, # Trinity at Dallas + "08066500": "08057000", # Trinity at Romayor + "08114000": None, # Brazos at Richmond + "08116650": "08114000", # Brazos near Rosharon + "08158000": None, # Colorado at Austin + "08167500": None, # Guadalupe near Spring Branch + "08176500": "08167500", # Guadalupe at Victoria + "08178000": None, # San Antonio River at San Antonio + # Pacific coast (MTR, SEW) + "11463500": None, # Russian River at Geyserville, CA + "11467000": "11463500", # Russian River at Hacienda Bridge, CA + "12027500": None, # Chehalis near Grand Mound, WA + "12144500": None, # Snoqualmie near Snoqualmie, WA + "12150800": "12144500", # Snohomish near Monroe, WA + "12194000": None, # Skagit near Concrete, WA + "12200500": "12194000", # Skagit near Mount Vernon, WA +} + +# The site service's county, where the NWS county UGC differs. Connecticut +# replaced its counties with planning regions in 2022 (FIPS 09110 = Capitol); +# NWS county UGCs still name the legacy counties (Enfield is in Hartford, CTC003). +COUNTY_UGC_OVERRIDE = {"01184000": "CTC003"} + +STATE_USPS = { + "06": "CA", "09": "CT", "13": "GA", "17": "IL", "18": "IN", "19": "IA", "21": "KY", + "24": "MD", "25": "MA", "29": "MO", "34": "NJ", "37": "NC", "39": "OH", "42": "PA", + "48": "TX", "50": "VT", "51": "VA", "53": "WA", "54": "WV", +} + +# Local standard time offset (hours behind UTC) by the site service's tz_cd. +TZ_HOURS = {"EST": 5, "CST": 6, "MST": 7, "PST": 8} diff --git a/apps/atlas-outcomes/atlas_outcomes/geo.py b/apps/atlas-outcomes/atlas_outcomes/geo.py new file mode 100644 index 00000000..71f40fa4 --- /dev/null +++ b/apps/atlas-outcomes/atlas_outcomes/geo.py @@ -0,0 +1,32 @@ +"""Point-in-polygon for warning polygons (lon/lat, even-odd rule, holes honoured).""" + +from __future__ import annotations + +Ring = list[list[float]] +Polygon = list[Ring] + + +def ring_contains(ring: Ring, x: float, y: float) -> bool: + inside = False + n = len(ring) + j = n - 1 + for i in range(n): + xi, yi = ring[i][0], ring[i][1] + xj, yj = ring[j][0], ring[j][1] + if (yi > y) != (yj > y) and x < (xj - xi) * (y - yi) / (yj - yi) + xi: + inside = not inside + j = i + return inside + + +def polygons_bbox(polys: list[Polygon]) -> tuple[float, float, float, float]: + xs = [p[0] for poly in polys for p in poly[0]] + ys = [p[1] for poly in polys for p in poly[0]] + return min(xs), min(ys), max(xs), max(ys) + + +def polygons_contain(polys: list[Polygon], lon: float, lat: float) -> bool: + for poly in polys: + if poly and ring_contains(poly[0], lon, lat) and not any(ring_contains(h, lon, lat) for h in poly[1:]): + return True + return False diff --git a/apps/atlas-outcomes/atlas_outcomes/harvest.py b/apps/atlas-outcomes/atlas_outcomes/harvest.py new file mode 100644 index 00000000..cb16c011 --- /dev/null +++ b/apps/atlas-outcomes/atlas_outcomes/harvest.py @@ -0,0 +1,227 @@ +"""Stage 1: harvest gauges, daily flow and NWS flood-warning polygons. + + USGS site service site metadata (name, lat/lon, county, HUC, drainage area, time zone) + IEM /api/1/nws/ugcs county UGC -> NWS office (WFO) and county name, per state + USGS NWIS daily values mean daily discharge (00060), 1995-10-01 -> --end + IEM /api/1/vtec/sbw_interval + storm-based (polygon) warnings issued in each gauge WFO, + flood (FL), areal flood (FA) and flash flood (FF), 2015 -> --end + +Raw responses are cached gzipped under data/raw (gitignored), so a rerun costs +nothing; the normalized, compact snapshots go to sources/ and are committed: + + sources/gauges.jsonl.gz one row per gauge + sources/flows.jsonl.gz one row per gauge: start date and a dense cfs array (null = no value) + sources/warnings.jsonl.gz one row per warning polygon at issuance (NEW), rings rounded to 1e-4 deg + sources/harvest_stats.json counts, the window, and the endpoints + +A gauge or office that fails contributes nothing and is listed in the stats; +nothing is filled in. +""" + +from __future__ import annotations + +import argparse +import json +from datetime import UTC, date, datetime, timedelta +from typing import Any + +from .common import SOURCES, HttpError, cached_get, stable_id, write_jsonl +from .gauges import COUNTY_UGC_OVERRIDE, GAUGES, STATE_USPS + +FLOW_START = "1995-10-01" +WARN_START_YEAR = 2015 +DEFAULT_END = "2026-09-22" +SITE_URL = "https://waterservices.usgs.gov/nwis/site/?format=rdb&siteOutput=expanded&sites={sites}" +DV_URL = "https://waterservices.usgs.gov/nwis/dv/?format=json&sites={site}¶meterCd=00060&statCd=00003&startDT={start}&endDT={end}" +UGC_URL = "https://mesonet.agron.iastate.edu/api/1/nws/ugcs.json?state={state}" +SBW_URL = "https://mesonet.agron.iastate.edu/api/1/vtec/sbw_interval.geojson?begints={b}T00:00Z&endts={e}T00:00Z&wfo={wfo}&only_new=true&{ph}" +# The endpoint takes at most two phenomena per call. +PH_GROUPS = ("ph=FL&ph=FA", "ph=FF") +FLOOD_PHENOMENA = {"FL", "FA", "FF"} + + +def parse_rdb(text: str) -> list[dict[str, str]]: + lines = [ln for ln in text.splitlines() if ln and not ln.startswith("#")] + if len(lines) < 2: + return [] + head = lines[0].split("\t") + return [dict(zip(head, ln.split("\t"), strict=False)) for ln in lines[2:]] + + +def parse_dv(payload: dict[str, Any]) -> dict[str, float | None]: + """date -> mean daily cfs from a NWIS dv JSON answer; the longest method series wins; no-data sentinels are None.""" + series = (payload.get("value") or {}).get("timeSeries") or [] + best: list[dict[str, Any]] = [] + nodata = None + for ts in series: + nodata = (ts.get("variable") or {}).get("noDataValue", nodata) + for block in ts.get("values") or []: + vals = block.get("value") or [] + if len(vals) > len(best): + best = vals + out: dict[str, float | None] = {} + for v in best: + day = str(v.get("dateTime", ""))[:10] + try: + x = float(v.get("value")) + except (TypeError, ValueError): + x = None + if x is not None and (x < 0 or (nodata is not None and x == float(nodata))): + x = None + if day: + out[day] = x + return out + + +def dense(values: dict[str, float | None]) -> tuple[str, list[float | None]]: + """(start date, contiguous daily array) with None for days the service did not report.""" + if not values: + return "", [] + days = sorted(values) + d0, d1 = date.fromisoformat(days[0]), date.fromisoformat(days[-1]) + out: list[float | None] = [] + d = d0 + while d <= d1: + x = values.get(d.isoformat()) + out.append(None if x is None else (round(x, 2) if x != int(x) else int(x))) + d += timedelta(days=1) + return d0.isoformat(), out + + +def _round_ring(ring: list[list[float]]) -> list[list[float]]: + return [[round(float(x), 4), round(float(y), 4)] for x, y, *_ in ring] + + +def normalize_geometry(geom: dict[str, Any] | None) -> list[list[list[list[float]]]]: + """GeoJSON Polygon/MultiPolygon -> a list of polygons, each [outer ring, holes...].""" + if not geom: + return [] + if geom.get("type") == "Polygon": + return [[_round_ring(r) for r in geom["coordinates"]]] + if geom.get("type") == "MultiPolygon": + return [[_round_ring(r) for r in poly] for poly in geom["coordinates"]] + return [] + + +def normalize_warning(feature: dict[str, Any]) -> dict[str, Any] | None: + p = feature.get("properties") or {} + if p.get("phenomena") not in FLOOD_PHENOMENA or p.get("significance") != "W": + return None + polys = normalize_geometry(feature.get("geometry")) + if not polys or not p.get("utc_issue") or not p.get("utc_expire"): + return None + return { + "wfo": p.get("wfo"), + "phenomena": p["phenomena"], + "significance": "W", + "eventid": int(p.get("eventid") or 0), + "year": int(p.get("year") or str(p["utc_issue"])[:4]), + "issue": p["utc_issue"], + "expire": p["utc_expire"], + "product_id": p.get("product_id"), + "ugcs": sorted(u.strip() for u in str(p.get("ugclist") or "").split(",") if u.strip()), + "polygons": polys, + } + + +def harvest(end: str = DEFAULT_END, refresh: bool = False) -> dict[str, Any]: + stats: dict[str, Any] = {"end": end, "flow_start": FLOW_START, "warn_start_year": WARN_START_YEAR, "failed": [], + "endpoints": [SITE_URL, DV_URL, UGC_URL, SBW_URL]} + sites = sorted(GAUGES) + meta = {r["site_no"]: r for r in parse_rdb(cached_get(SITE_URL.format(sites=",".join(sites)), f"site/expanded_{stable_id(*sites)}", refresh=refresh).decode("utf-8"))} + ugc: dict[str, dict[str, str]] = {} + for fips, usps in sorted(STATE_USPS.items()): + try: + rows = json.loads(cached_get(UGC_URL.format(state=usps), f"ugc/{usps}", refresh=refresh, host_gap_s=1.0))["data"] + except (HttpError, OSError, ValueError, KeyError) as e: + stats["failed"].append({"ugc": usps, "error": str(e)[:120]}) + continue + for r in rows: + ugc[r["ugc"]] = {"wfo": r["wfo"], "name": r["name"], "state": fips} + + gauges: list[dict[str, Any]] = [] + flows: list[dict[str, Any]] = [] + for site in sites: + m = meta.get(site) + if not m: + stats["failed"].append({"site": site, "error": "not in the site service"}) + continue + usps = STATE_USPS.get(m["state_cd"], "") + county_ugc = COUNTY_UGC_OVERRIDE.get(site, f"{usps}C{m['county_cd']}") + u = ugc.get(county_ugc) + if not u: + stats["failed"].append({"site": site, "error": f"no UGC {county_ugc}"}) + continue + try: + payload = json.loads(cached_get(DV_URL.format(site=site, start=FLOW_START, end=end), f"dv/{site}_{FLOW_START}_{end}", refresh=refresh, host_gap_s=1.0)) + except (HttpError, OSError, ValueError) as e: + stats["failed"].append({"site": site, "error": str(e)[:120]}) + continue + start, cfs = dense(parse_dv(payload)) + if not cfs: + stats["failed"].append({"site": site, "error": "no daily values"}) + continue + drain = m.get("drain_area_va", "").strip() + gauges.append({ + "site": site, + "name": m["station_nm"].strip(), + "lat": round(float(m["dec_lat_va"]), 5), + "lon": round(float(m["dec_long_va"]), 5), + "state": usps, + "county_ugc": county_ugc, + "county": u["name"], + "wfo": u["wfo"], + "huc8": m["huc_cd"][:8], + "drainage_sq_mi": float(drain) if drain else None, + "tz": m["tz_cd"], + "upstream": GAUGES[site], + }) + flows.append({"site": site, "start": start, "cfs": cfs}) + + wfos = sorted({g["wfo"] for g in gauges}) + end_d = date.fromisoformat(end) + warnings: dict[tuple, dict[str, Any]] = {} + for wfo in wfos: + for year in range(WARN_START_YEAR, end_d.year + 1): + b = f"{year}-01-01" + e = min(date(year + 1, 1, 1), end_d + timedelta(days=2)).isoformat() + for ph in PH_GROUPS: + tag = ph.replace("ph=", "").replace("&", "") + try: + body = cached_get(SBW_URL.format(b=b, e=e, wfo=wfo, ph=ph), f"sbw/{wfo}_{year}_{tag}_{e}", refresh=refresh, host_gap_s=1.0) + feats = json.loads(body).get("features") or [] + except (HttpError, OSError, ValueError) as err: + stats["failed"].append({"wfo": wfo, "year": year, "ph": tag, "error": str(err)[:120]}) + continue + for f in feats: + w = normalize_warning(f) + if w: + warnings[(w["wfo"], w["phenomena"], w["year"], w["eventid"], w["product_id"])] = w + rows = sorted(warnings.values(), key=lambda w: (w["issue"], w["wfo"], w["phenomena"], w["eventid"], w["product_id"] or "")) + + SOURCES.mkdir(parents=True, exist_ok=True) + write_jsonl(SOURCES / "gauges.jsonl.gz", gauges) + write_jsonl(SOURCES / "flows.jsonl.gz", flows) + write_jsonl(SOURCES / "warnings.jsonl.gz", rows) + stats.update({ + "gauges": len(gauges), + "paired": sum(1 for g in gauges if g["upstream"]), + "wfos": wfos, + "flow_days": sum(sum(1 for x in f["cfs"] if x is not None) for f in flows), + "warnings": len(rows), + "warnings_by_phenomena": {p: sum(1 for w in rows if w["phenomena"] == p) for p in sorted(FLOOD_PHENOMENA)}, + "harvested_at": datetime.now(UTC).strftime("%Y-%m-%dT%H:%M:%SZ"), + }) + (SOURCES / "harvest_stats.json").write_text(json.dumps(stats, indent=2, sort_keys=True) + "\n", encoding="utf-8") + return stats + + +def main(argv: list[str] | None = None) -> int: + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--end", default=DEFAULT_END, help="last day of flow and warnings to fetch (YYYY-MM-DD)") + ap.add_argument("--refresh", action="store_true", help="ignore the data/raw cache") + args = ap.parse_args(argv) + s = harvest(args.end, args.refresh) + print(json.dumps({k: s[k] for k in ("gauges", "paired", "wfos", "flow_days", "warnings", "warnings_by_phenomena", "failed")}, indent=2)) + return 0 if s["gauges"] else 1 diff --git a/apps/atlas-outcomes/atlas_outcomes/records.py b/apps/atlas-outcomes/atlas_outcomes/records.py new file mode 100644 index 00000000..86d5ff62 --- /dev/null +++ b/apps/atlas-outcomes/atlas_outcomes/records.py @@ -0,0 +1,128 @@ +"""Render a feature row as a strict jev-decision-schema-1.0.0 record. + +The state is what a place decision on eye.jcamd.com sees at the end of day t: +the gauge, its construct stack (county, HUC-8, NWS office, state, named as the +Atlas strata rail names them) with the layers active there (the streamflow +condition class, whether a flood-type warning polygon is in effect over the +gauge), the flow and its recent change, the upstream gauge, and recent +warning activity. It never carries t+1's flow, the next day's warnings or any +label; tests/test_outcomes.py checks that. + +Three questions per record, all answered by what happened next: + + high_next noul tomorrow's mean flow above the stated p90 threshold + warn_next noul a flood-type warning polygon over the gauge issued in the next 24 h + flow_change score tomorrow's change in five ordered levels +""" + +from __future__ import annotations + +from typing import Any + +from . import SCHEMA_ID +from .common import stable_id +from .features import CHANGE_LEVELS, flow_class + + +def _cfs(x: float) -> str: + return f"{x:,.0f} cfs" if x >= 10 else f"{x:g} cfs" + + +def _pct(x: float | None) -> str: + return "not reported" if x is None else f"{x:+.1f}%" + + +def construct_stack(g: dict[str, Any], f: dict[str, Any]) -> dict[str, Any]: + """Smallest construct first, like the Atlas strata rail; the layers active at the gauge.""" + return { + "constructs": [ + {"kind": "County", "id": g["county_ugc"], "name": f"{g['county']}, {g['state']}"}, + {"kind": "HUC-8", "id": g["huc8"]}, + {"kind": "NWS office", "id": g["wfo"]}, + {"kind": "State", "id": g["state"]}, + ], + "layers": { + "water": f"streamflow {flow_class(f['pct_date'])} for the date", + "nws-alert": "flood-type warning in effect over the gauge" if f["warn_active"] else "no flood-type warning over the gauge", + **({"upstream water": f"streamflow {flow_class(f.get('up_pct_date'))} for the date"} if f.get("up_site") else {}), + }, + } + + +def render_state(g: dict[str, Any], row: dict[str, Any]) -> dict[str, Any]: + f = row["features"] + state: dict[str, Any] = { + "place": { + "gauge": f"USGS {g['site']}", + "name": g["name"], + "drainage_area": "not published" if g.get("drainage_sq_mi") is None else f"{g['drainage_sq_mi']:,.0f} sq mi", + }, + "as_of": f"end of day {row['date']} local time ({f['season']})", + **construct_stack(g, f), + "flow": { + "today": _cfs(f["flow"]), + "percentile_for_the_date": f["pct_date"], + "percentile_of_all_prior_days": f["pct_all"], + "change_vs_1_day_ago": _pct(f["chg1"]), + "change_vs_3_days_ago": _pct(f["chg3"]), + "change_vs_7_days_ago": _pct(f["chg7"]), + }, + "warnings": { + "flood_type_warnings_over_gauge_last_30_days": f["warn_gauge_30d"], + "flood_type_warnings_by_office_last_7_days": f["warn_wfo_7d"], + "flood_type_warnings_by_office_last_24_hours": f["warn_wfo_1d"], + }, + } + if f.get("up_site"): + state["upstream"] = { + "gauge": f"USGS {f['up_site']}", + "name": f["up_name"], + "percentile_for_the_date": f["up_pct_date"], + "change_vs_1_day_ago": _pct(f["up_chg1"]), + "change_vs_3_days_ago": _pct(f["up_chg3"]), + } + else: + state["upstream"] = "no paired upstream gauge" + return state + + +def questions(row: dict[str, Any]) -> dict[str, Any]: + p90 = row["features"]["p90"] + return { + "high_next": { + "type": "noul", + "instructions": f"Will tomorrow's mean daily flow at this gauge be above {_cfs(p90)}, the 90th percentile of every prior day on its record?", + }, + "warn_next": { + "type": "noul", + "instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", + }, + "flow_change": { + "type": "score", + "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", + "criteria": list(CHANGE_LEVELS), + }, + } + + +def labels(row: dict[str, Any]) -> dict[str, Any]: + lab = row["labels"] + return { + "high_next": {"type": "noul", "noul": 1.0 if lab["high_next"] else 0.0}, + "warn_next": {"type": "noul", "noul": 1.0 if lab["warn_next"] else 0.0}, + "flow_change": {"type": "score", "score": float(lab["change_bucket"])}, + } + + +def record_id(row: dict[str, Any]) -> str: + return f"atlas-{row['site']}-{row['date']}-{stable_id(row['site'], row['date'], n=6)}" + + +def to_record(g: dict[str, Any], row: dict[str, Any]) -> dict[str, Any]: + return { + "schema": SCHEMA_ID, + "id": record_id(row), + "state": render_state(g, row), + "questions": questions(row), + "labels": labels(row), + } diff --git a/apps/atlas-outcomes/atlas_outcomes/train_command.py b/apps/atlas-outcomes/atlas_outcomes/train_command.py new file mode 100644 index 00000000..82539505 --- /dev/null +++ b/apps/atlas-outcomes/atlas_outcomes/train_command.py @@ -0,0 +1,64 @@ +"""Stage 5: print the GPU-host System One training command for this pack. Runs nothing. + +The trainer is apps/jev-v0/train_pointer_lora.py (the same argv +dottie_loop.router_training.train_command builds for router packs): `--go` +trains LoRA + pointer heads on `--fixtures`, one question per step, walking the +file in order (curate shuffles train.jsonl for that reason). The default step +count is one pass over every question of every train record. +""" + +from __future__ import annotations + +import argparse +import json +import shlex + +from . import PACK_VERSION +from .common import APP_ROOT, DATA + +REPO = APP_ROOT.parent.parent + + +def command(pack: str = PACK_VERSION, steps: int | None = None, base_model: str | None = None) -> dict[str, object]: + pack_dir = DATA / "packs" / pack + manifest = json.loads((pack_dir / "MANIFEST.json").read_text(encoding="utf-8")) + questions = len(manifest["questions"]) + n_steps = steps or manifest["rows"]["train"] * questions + rel = pack_dir.relative_to(REPO).as_posix() + out = f"apps/jev-v0/runs/{pack}" + go = ["python", "apps/jev-v0/train_pointer_lora.py", "--go", "--fixtures", f"{rel}/train.jsonl", + "--out", out, "--steps", str(n_steps)] + if base_model: + go += ["--base-model", base_model] + return { + "pack": pack, + "train_rows": manifest["rows"]["train"], + "questions_per_row": questions, + "steps": n_steps, + "setup": "pip install -r apps/jev-v0/requirements-jev-v0.txt # plus the host CUDA torch wheel", + "rebuild": f"python3 apps/atlas-outcomes/run.py curate # rebuilds {rel}/ offline from sources/", + "dry_run": shlex.join(["python", "apps/jev-v0/train_pointer_lora.py", "--dry-run", "--fixtures", f"{rel}/train.jsonl"]), + "train": shlex.join(go), + "serve": shlex.join(["python", "apps/jev-v0/serve_decide.py", "--checkpoint", out, "--port", "8771"]), + "eval": f"score {rel}/holdout.jsonl with the served checkpoint and compare against apps/atlas-outcomes/BASELINE.json " + "(python3 apps/atlas-outcomes/run.py baseline --score ); nothing promotes without a human stamp", + } + + +def main(argv: list[str] | None = None) -> int: + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--pack", default=PACK_VERSION) + ap.add_argument("--steps", type=int, default=None) + ap.add_argument("--base-model", default=None) + ap.add_argument("--json", action="store_true") + args = ap.parse_args(argv) + c = command(args.pack, args.steps, args.base_model) + if args.json: + print(json.dumps(c, indent=2)) + return 0 + print("# GPU host, from the repo root. Not run here.") + print(f"# pack {c['pack']}: {c['train_rows']} train rows x {c['questions_per_row']} questions = {c['steps']} steps (one pass)") + for k in ("setup", "rebuild", "dry_run", "train", "serve"): + print(c[k]) + print(f"# eval: {c['eval']}") + return 0 diff --git a/apps/atlas-outcomes/run.py b/apps/atlas-outcomes/run.py new file mode 100644 index 00000000..5f37e69e --- /dev/null +++ b/apps/atlas-outcomes/run.py @@ -0,0 +1,25 @@ +"""Entry point: `python3 apps/atlas-outcomes/run.py [args]` from anywhere. + +Stages, in order: harvest, features, curate, baseline, train-command. +""" + +from __future__ import annotations + +import importlib +import sys +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parent)) + +STAGES = { + "harvest": "atlas_outcomes.harvest", + "features": "atlas_outcomes.features", + "curate": "atlas_outcomes.curate", + "baseline": "atlas_outcomes.baseline", + "train-command": "atlas_outcomes.train_command", +} + +if __name__ == "__main__": + if len(sys.argv) < 2 or sys.argv[1] not in STAGES: + raise SystemExit(f"usage: run.py {{{'|'.join(STAGES)}}} [args]") + raise SystemExit(importlib.import_module(STAGES[sys.argv[1]]).main(sys.argv[2:])) diff --git a/apps/atlas-outcomes/sample/pack_sample.jsonl b/apps/atlas-outcomes/sample/pack_sample.jsonl new file mode 100644 index 00000000..8ddc98a9 --- /dev/null +++ b/apps/atlas-outcomes/sample/pack_sample.jsonl @@ -0,0 +1,16 @@ +{"provenance": {"date": "2015-01-01", "decision_utc": "2015-01-02T05:00Z", "easy_negative": true, "id": "atlas-01610000-2015-01-01-ac21a6", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 10.0, "site": "01610000", "split": "train", "state": "MD", "wfo": "LWX"}, "record": {"id": "atlas-01610000-2015-01-01-ac21a6", "labels": {"flow_change": {"score": 2.0, "type": "score"}, "high_next": {"noul": 0.0, "type": "noul"}, "warn_next": {"noul": 0.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 8,340 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2015-01-01 local time (winter)", "constructs": [{"id": "MDC001", "kind": "County", "name": "Allegany, MD"}, {"id": "02070003", "kind": "HUC-8"}, {"id": "LWX", "kind": "NWS office"}, {"id": "MD", "kind": "State"}], "flow": {"change_vs_1_day_ago": "-13.2%", "change_vs_3_days_ago": "-23.7%", "change_vs_7_days_ago": "-34.9%", "percentile_for_the_date": 24.2, "percentile_of_all_prior_days": 41.5, "today": "1,510 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "water": "streamflow below normal for the date"}, "place": {"drainage_area": "3,129 sq mi", "gauge": "USGS 01610000", "name": "POTOMAC RIVER AT PAW PAW, WV"}, "upstream": "no paired upstream gauge", "warnings": {"flood_type_warnings_by_office_last_24_hours": 0, "flood_type_warnings_by_office_last_7_days": 0, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2015-01-01", "decision_utc": "2015-01-02T05:00Z", "easy_negative": false, "id": "atlas-02087500-2015-01-01-90bd42", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "02087500", "split": "train", "state": "NC", "wfo": "RAH"}, "record": {"id": "atlas-02087500-2015-01-01-90bd42", "labels": {"flow_change": {"score": 2.0, "type": "score"}, "high_next": {"noul": 1.0, "type": "noul"}, "warn_next": {"noul": 0.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 2,750 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2015-01-01 local time (winter)", "constructs": [{"id": "NCC101", "kind": "County", "name": "Johnston, NC"}, {"id": "03020201", "kind": "HUC-8"}, {"id": "RAH", "kind": "NWS office"}, {"id": "NC", "kind": "State"}], "flow": {"change_vs_1_day_ago": "+23.8%", "change_vs_3_days_ago": "+126.8%", "change_vs_7_days_ago": "-23.4%", "percentile_for_the_date": 93.5, "percentile_of_all_prior_days": 94.2, "today": "3,380 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "water": "streamflow much above normal for the date"}, "place": {"drainage_area": "1,150 sq mi", "gauge": "USGS 02087500", "name": "NEUSE RIVER NEAR CLAYTON, NC"}, "upstream": "no paired upstream gauge", "warnings": {"flood_type_warnings_by_office_last_24_hours": 0, "flood_type_warnings_by_office_last_7_days": 0, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2015-01-01", "decision_utc": "2015-01-02T05:00Z", "easy_negative": false, "id": "atlas-02089500-2015-01-01-859f8e", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "02089500", "split": "train", "state": "NC", "wfo": "MHX"}, "record": {"id": "atlas-02089500-2015-01-01-859f8e", "labels": {"flow_change": {"score": 2.0, "type": "score"}, "high_next": {"noul": 1.0, "type": "noul"}, "warn_next": {"noul": 0.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 6,080 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2015-01-01 local time (winter)", "constructs": [{"id": "NCC107", "kind": "County", "name": "Lenoir, NC"}, {"id": "03020202", "kind": "HUC-8"}, {"id": "MHX", "kind": "NWS office"}, {"id": "NC", "kind": "State"}], "flow": {"change_vs_1_day_ago": "+7.8%", "change_vs_3_days_ago": "+40.0%", "change_vs_7_days_ago": "+164.9%", "percentile_for_the_date": 100.0, "percentile_of_all_prior_days": 98.2, "today": "11,100 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "upstream water": "streamflow much above normal for the date", "water": "streamflow much above normal for the date"}, "place": {"drainage_area": "2,692 sq mi", "gauge": "USGS 02089500", "name": "NEUSE RIVER AT KINSTON, NC"}, "upstream": {"change_vs_1_day_ago": "+23.8%", "change_vs_3_days_ago": "+126.8%", "gauge": "USGS 02087500", "name": "NEUSE RIVER NEAR CLAYTON, NC", "percentile_for_the_date": 93.5}, "warnings": {"flood_type_warnings_by_office_last_24_hours": 0, "flood_type_warnings_by_office_last_7_days": 0, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2015-01-01", "decision_utc": "2015-01-02T05:00Z", "easy_negative": true, "id": "atlas-03234500-2015-01-01-469b48", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 10.0, "site": "03234500", "split": "train", "state": "OH", "wfo": "ILN"}, "record": {"id": "atlas-03234500-2015-01-01-469b48", "labels": {"flow_change": {"score": 1.0, "type": "score"}, "high_next": {"noul": 0.0, "type": "noul"}, "warn_next": {"noul": 0.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 15,300 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2015-01-01 local time (winter)", "constructs": [{"id": "OHC141", "kind": "County", "name": "Ross, OH"}, {"id": "05060002", "kind": "HUC-8"}, {"id": "ILN", "kind": "NWS office"}, {"id": "OH", "kind": "State"}], "flow": {"change_vs_1_day_ago": "-19.4%", "change_vs_3_days_ago": "-17.1%", "change_vs_7_days_ago": "+27.0%", "percentile_for_the_date": 13.0, "percentile_of_all_prior_days": 32.6, "today": "1,790 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "water": "streamflow below normal for the date"}, "place": {"drainage_area": "5,131 sq mi", "gauge": "USGS 03234500", "name": "Scioto River at Higby OH"}, "upstream": "no paired upstream gauge", "warnings": {"flood_type_warnings_by_office_last_24_hours": 0, "flood_type_warnings_by_office_last_7_days": 0, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2015-01-04", "decision_utc": "2015-01-05T08:00Z", "easy_negative": false, "id": "atlas-12027500-2015-01-04-da1f5b", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "12027500", "split": "train", "state": "WA", "wfo": "SEW"}, "record": {"id": "atlas-12027500-2015-01-04-da1f5b", "labels": {"flow_change": {"score": 4.0, "type": "score"}, "high_next": {"noul": 1.0, "type": "noul"}, "warn_next": {"noul": 1.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 7,660 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2015-01-04 local time (winter)", "constructs": [{"id": "WAC067", "kind": "County", "name": "Thurston, WA"}, {"id": "17100103", "kind": "HUC-8"}, {"id": "SEW", "kind": "NWS office"}, {"id": "WA", "kind": "State"}], "flow": {"change_vs_1_day_ago": "-3.1%", "change_vs_3_days_ago": "-18.7%", "change_vs_7_days_ago": "-48.8%", "percentile_for_the_date": 21.1, "percentile_of_all_prior_days": 66.2, "today": "2,790 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "water": "streamflow below normal for the date"}, "place": {"drainage_area": "895 sq mi", "gauge": "USGS 12027500", "name": "CHEHALIS RIVER NEAR GRAND MOUND, WA"}, "upstream": "no paired upstream gauge", "warnings": {"flood_type_warnings_by_office_last_24_hours": 2, "flood_type_warnings_by_office_last_7_days": 2, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2015-01-04", "decision_utc": "2015-01-05T08:00Z", "easy_negative": false, "id": "atlas-12144500-2015-01-04-b6a1b4", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "12144500", "split": "train", "state": "WA", "wfo": "SEW"}, "record": {"id": "atlas-12144500-2015-01-04-b6a1b4", "labels": {"flow_change": {"score": 4.0, "type": "score"}, "high_next": {"noul": 1.0, "type": "noul"}, "warn_next": {"noul": 1.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 5,360 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2015-01-04 local time (winter)", "constructs": [{"id": "WAC033", "kind": "County", "name": "King, WA"}, {"id": "17110010", "kind": "HUC-8"}, {"id": "SEW", "kind": "NWS office"}, {"id": "WA", "kind": "State"}], "flow": {"change_vs_1_day_ago": "+1.1%", "change_vs_3_days_ago": "-8.5%", "change_vs_7_days_ago": "-55.9%", "percentile_for_the_date": 33.0, "percentile_of_all_prior_days": 43.5, "today": "1,840 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "water": "streamflow normal for the date"}, "place": {"drainage_area": "375 sq mi", "gauge": "USGS 12144500", "name": "SNOQUALMIE RIVER NEAR SNOQUALMIE, WA"}, "upstream": "no paired upstream gauge", "warnings": {"flood_type_warnings_by_office_last_24_hours": 2, "flood_type_warnings_by_office_last_7_days": 2, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2015-01-17", "decision_utc": "2015-01-18T05:00Z", "easy_negative": false, "id": "atlas-01463500-2015-01-17-02bde2", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "01463500", "split": "train", "state": "NJ", "wfo": "PHI"}, "record": {"id": "atlas-01463500-2015-01-17-02bde2", "labels": {"flow_change": {"score": 1.0, "type": "score"}, "high_next": {"noul": 0.0, "type": "noul"}, "warn_next": {"noul": 1.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 27,300 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2015-01-17 local time (winter)", "constructs": [{"id": "NJC021", "kind": "County", "name": "Mercer, NJ"}, {"id": "02040105", "kind": "HUC-8"}, {"id": "PHI", "kind": "NWS office"}, {"id": "NJ", "kind": "State"}], "flow": {"change_vs_1_day_ago": "-17.3%", "change_vs_3_days_ago": "-24.4%", "change_vs_7_days_ago": "-27.1%", "percentile_for_the_date": 16.7, "percentile_of_all_prior_days": 29.6, "today": "6,200 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "upstream water": "streamflow below normal for the date", "water": "streamflow below normal for the date"}, "place": {"drainage_area": "6,780 sq mi", "gauge": "USGS 01463500", "name": "Delaware River at Trenton NJ"}, "upstream": {"change_vs_1_day_ago": "-14.1%", "change_vs_3_days_ago": "-19.9%", "gauge": "USGS 01446500", "name": "Delaware River at Belvidere NJ", "percentile_for_the_date": 17.2}, "warnings": {"flood_type_warnings_by_office_last_24_hours": 0, "flood_type_warnings_by_office_last_7_days": 0, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2015-01-17", "decision_utc": "2015-01-18T05:00Z", "easy_negative": false, "id": "atlas-01474500-2015-01-17-dc2f7b", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "01474500", "split": "train", "state": "PA", "wfo": "PHI"}, "record": {"id": "atlas-01474500-2015-01-17-dc2f7b", "labels": {"flow_change": {"score": 4.0, "type": "score"}, "high_next": {"noul": 0.0, "type": "noul"}, "warn_next": {"noul": 1.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 6,620 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2015-01-17 local time (winter)", "constructs": [{"id": "PAC101", "kind": "County", "name": "Philadelphia, PA"}, {"id": "02040203", "kind": "HUC-8"}, {"id": "PHI", "kind": "NWS office"}, {"id": "PA", "kind": "State"}], "flow": {"change_vs_1_day_ago": "-12.6%", "change_vs_3_days_ago": "-36.5%", "change_vs_7_days_ago": "-12.6%", "percentile_for_the_date": 18.8, "percentile_of_all_prior_days": 35.4, "today": "1,600 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "water": "streamflow below normal for the date"}, "place": {"drainage_area": "1,893 sq mi", "gauge": "USGS 01474500", "name": "Schuylkill River at Philadelphia, PA"}, "upstream": "no paired upstream gauge", "warnings": {"flood_type_warnings_by_office_last_24_hours": 0, "flood_type_warnings_by_office_last_7_days": 0, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2025-09-22", "decision_utc": "2025-09-23T05:00Z", "easy_negative": true, "id": "atlas-01100000-2025-09-22-df0d4a", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "01100000", "split": "holdout", "state": "MA", "wfo": "BOX"}, "record": {"id": "atlas-01100000-2025-09-22-df0d4a", "labels": {"flow_change": {"score": 1.0, "type": "score"}, "high_next": {"noul": 0.0, "type": "noul"}, "warn_next": {"noul": 0.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 19,100 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2025-09-22 local time (fall)", "constructs": [{"id": "MAC017", "kind": "County", "name": "Middlesex, MA"}, {"id": "01070006", "kind": "HUC-8"}, {"id": "BOX", "kind": "NWS office"}, {"id": "MA", "kind": "State"}], "flow": {"change_vs_1_day_ago": "-7.2%", "change_vs_3_days_ago": "-14.6%", "change_vs_7_days_ago": "-8.1%", "percentile_for_the_date": 0.5, "percentile_of_all_prior_days": 0.1, "today": "753 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "water": "streamflow much below normal for the date"}, "place": {"drainage_area": "4,635 sq mi", "gauge": "USGS 01100000", "name": "MERRIMACK RIVER BL CONCORD RIVER AT LOWELL, MA"}, "upstream": "no paired upstream gauge", "warnings": {"flood_type_warnings_by_office_last_24_hours": 0, "flood_type_warnings_by_office_last_7_days": 0, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2025-09-22", "decision_utc": "2025-09-23T05:00Z", "easy_negative": true, "id": "atlas-01184000-2025-09-22-4b01ca", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "01184000", "split": "holdout", "state": "CT", "wfo": "BOX"}, "record": {"id": "atlas-01184000-2025-09-22-4b01ca", "labels": {"flow_change": {"score": 1.0, "type": "score"}, "high_next": {"noul": 0.0, "type": "noul"}, "warn_next": {"noul": 0.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 39,900 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2025-09-22 local time (fall)", "constructs": [{"id": "CTC003", "kind": "County", "name": "Hartford, CT"}, {"id": "01080205", "kind": "HUC-8"}, {"id": "BOX", "kind": "NWS office"}, {"id": "CT", "kind": "State"}], "flow": {"change_vs_1_day_ago": "+38.5%", "change_vs_3_days_ago": "+23.6%", "change_vs_7_days_ago": "+59.3%", "percentile_for_the_date": 27.7, "percentile_of_all_prior_days": 6.5, "today": "4,030 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "water": "streamflow normal for the date"}, "place": {"drainage_area": "9,660 sq mi", "gauge": "USGS 01184000", "name": "CONNECTICUT RIVER AT THOMPSONVILLE, CT"}, "upstream": "no paired upstream gauge", "warnings": {"flood_type_warnings_by_office_last_24_hours": 0, "flood_type_warnings_by_office_last_7_days": 0, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2025-10-06", "decision_utc": "2025-10-07T06:00Z", "easy_negative": false, "id": "atlas-03303280-2025-10-06-490c42", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "03303280", "split": "holdout", "state": "IN", "wfo": "LMK"}, "record": {"id": "atlas-03303280-2025-10-06-490c42", "labels": {"flow_change": {"score": 4.0, "type": "score"}, "high_next": {"noul": 0.0, "type": "noul"}, "warn_next": {"noul": 1.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 302,000 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2025-10-06 local time (fall)", "constructs": [{"id": "INC123", "kind": "County", "name": "Perry, IN"}, {"id": "05140201", "kind": "HUC-8"}, {"id": "LMK", "kind": "NWS office"}, {"id": "IN", "kind": "State"}], "flow": {"change_vs_1_day_ago": "-12.5%", "change_vs_3_days_ago": "+36.7%", "change_vs_7_days_ago": "+166.8%", "percentile_for_the_date": 15.2, "percentile_of_all_prior_days": 5.1, "today": "17,500 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "upstream water": "streamflow below normal for the date", "water": "streamflow below normal for the date"}, "place": {"drainage_area": "97,000 sq mi", "gauge": "USGS 03303280", "name": "OHIO RIVER AT CANNELTON DAM AT CANNELTON, IN"}, "upstream": {"change_vs_1_day_ago": "-11.2%", "change_vs_3_days_ago": "+2.9%", "gauge": "USGS 03294500", "name": "OHIO RIVER AT LOUISVILLE, KY", "percentile_for_the_date": 18.4}, "warnings": {"flood_type_warnings_by_office_last_24_hours": 0, "flood_type_warnings_by_office_last_7_days": 0, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2025-10-24", "decision_utc": "2025-10-25T06:00Z", "easy_negative": false, "id": "atlas-08057000-2025-10-24-505292", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "08057000", "split": "holdout", "state": "TX", "wfo": "FWD"}, "record": {"id": "atlas-08057000-2025-10-24-505292", "labels": {"flow_change": {"score": 4.0, "type": "score"}, "high_next": {"noul": 1.0, "type": "noul"}, "warn_next": {"noul": 0.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 7,030 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2025-10-24 local time (fall)", "constructs": [{"id": "TXC113", "kind": "County", "name": "Dallas, TX"}, {"id": "12030105", "kind": "HUC-8"}, {"id": "FWD", "kind": "NWS office"}, {"id": "TX", "kind": "State"}], "flow": {"change_vs_1_day_ago": "+422.0%", "change_vs_3_days_ago": "+394.4%", "change_vs_7_days_ago": "+456.2%", "percentile_for_the_date": 80.9, "percentile_of_all_prior_days": 70.3, "today": "1,780 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "water": "streamflow above normal for the date"}, "place": {"drainage_area": "6,106 sq mi", "gauge": "USGS 08057000", "name": "Trinity Rv at Dallas, TX"}, "upstream": "no paired upstream gauge", "warnings": {"flood_type_warnings_by_office_last_24_hours": 3, "flood_type_warnings_by_office_last_7_days": 3, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2025-10-24", "decision_utc": "2025-10-25T06:00Z", "easy_negative": false, "id": "atlas-08178000-2025-10-24-b085ed", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "08178000", "split": "holdout", "state": "TX", "wfo": "EWX"}, "record": {"id": "atlas-08178000-2025-10-24-b085ed", "labels": {"flow_change": {"score": 4.0, "type": "score"}, "high_next": {"noul": 1.0, "type": "noul"}, "warn_next": {"noul": 0.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 38 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2025-10-24 local time (fall)", "constructs": [{"id": "TXC029", "kind": "County", "name": "Bexar, TX"}, {"id": "12100301", "kind": "HUC-8"}, {"id": "EWX", "kind": "NWS office"}, {"id": "TX", "kind": "State"}], "flow": {"change_vs_1_day_ago": "+32.7%", "change_vs_3_days_ago": "+75.7%", "change_vs_7_days_ago": "+229.9%", "percentile_for_the_date": 25.3, "percentile_of_all_prior_days": 22.2, "today": "20 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "water": "streamflow normal for the date"}, "place": {"drainage_area": "42 sq mi", "gauge": "USGS 08178000", "name": "San Antonio Rv at San Antonio, TX"}, "upstream": "no paired upstream gauge", "warnings": {"flood_type_warnings_by_office_last_24_hours": 0, "flood_type_warnings_by_office_last_7_days": 0, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2025-11-19", "decision_utc": "2025-11-20T06:00Z", "easy_negative": false, "id": "atlas-08057000-2025-11-19-4724e9", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "08057000", "split": "holdout", "state": "TX", "wfo": "FWD"}, "record": {"id": "atlas-08057000-2025-11-19-4724e9", "labels": {"flow_change": {"score": 4.0, "type": "score"}, "high_next": {"noul": 1.0, "type": "noul"}, "warn_next": {"noul": 1.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 7,030 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2025-11-19 local time (fall)", "constructs": [{"id": "TXC113", "kind": "County", "name": "Dallas, TX"}, {"id": "12030105", "kind": "HUC-8"}, {"id": "FWD", "kind": "NWS office"}, {"id": "TX", "kind": "State"}], "flow": {"change_vs_1_day_ago": "+18.6%", "change_vs_3_days_ago": "+14.1%", "change_vs_7_days_ago": "+15.7%", "percentile_for_the_date": 23.2, "percentile_of_all_prior_days": 14.3, "today": "420 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "water": "streamflow below normal for the date"}, "place": {"drainage_area": "6,106 sq mi", "gauge": "USGS 08057000", "name": "Trinity Rv at Dallas, TX"}, "upstream": "no paired upstream gauge", "warnings": {"flood_type_warnings_by_office_last_24_hours": 0, "flood_type_warnings_by_office_last_7_days": 0, "flood_type_warnings_over_gauge_last_30_days": 0}}}} +{"provenance": {"date": "2025-11-23", "decision_utc": "2025-11-24T06:00Z", "easy_negative": false, "id": "atlas-08057000-2025-11-23-0f88c1", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "08057000", "split": "holdout", "state": "TX", "wfo": "FWD"}, "record": {"id": "atlas-08057000-2025-11-23-0f88c1", "labels": {"flow_change": {"score": 4.0, "type": "score"}, "high_next": {"noul": 0.0, "type": "noul"}, "warn_next": {"noul": 1.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 7,030 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2025-11-23 local time (fall)", "constructs": [{"id": "TXC113", "kind": "County", "name": "Dallas, TX"}, {"id": "12030105", "kind": "HUC-8"}, {"id": "FWD", "kind": "NWS office"}, {"id": "TX", "kind": "State"}], "flow": {"change_vs_1_day_ago": "-66.4%", "change_vs_3_days_ago": "-80.1%", "change_vs_7_days_ago": "+288.6%", "percentile_for_the_date": 79.3, "percentile_of_all_prior_days": 66.9, "today": "1,430 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "water": "streamflow above normal for the date"}, "place": {"drainage_area": "6,106 sq mi", "gauge": "USGS 08057000", "name": "Trinity Rv at Dallas, TX"}, "upstream": "no paired upstream gauge", "warnings": {"flood_type_warnings_by_office_last_24_hours": 0, "flood_type_warnings_by_office_last_7_days": 6, "flood_type_warnings_over_gauge_last_30_days": 2}}}} +{"provenance": {"date": "2025-12-07", "decision_utc": "2025-12-08T08:00Z", "easy_negative": false, "id": "atlas-12144500-2025-12-07-ce98aa", "label_sources": {"flow_change": "usgs-nwis-dv", "high_next": "usgs-nwis-dv", "warn_next": "iem-vtec-sbw"}, "provenance": "outcome-real", "sample_weight": 1.0, "site": "12144500", "split": "holdout", "state": "WA", "wfo": "SEW"}, "record": {"id": "atlas-12144500-2025-12-07-ce98aa", "labels": {"flow_change": {"score": 4.0, "type": "score"}, "high_next": {"noul": 1.0, "type": "noul"}, "warn_next": {"noul": 1.0, "type": "noul"}}, "questions": {"flow_change": {"criteria": ["falls more than 20%", "falls 5% to 20%", "within 5% either way", "rises 5% to 20%", "rises more than 20%"], "instructions": "How will tomorrow's mean daily flow at this gauge compare with today's?", "type": "score"}, "high_next": {"instructions": "Will tomorrow's mean daily flow at this gauge be above 5,280 cfs, the 90th percentile of every prior day on its record?", "type": "noul"}, "warn_next": {"instructions": "Will the National Weather Service issue a flood, areal flood or flash flood warning whose polygon covers this gauge within the next 24 hours?", "type": "noul"}}, "schema": "jev-decision-schema-1.0.0", "state": {"as_of": "end of day 2025-12-07 local time (winter)", "constructs": [{"id": "WAC033", "kind": "County", "name": "King, WA"}, {"id": "17110010", "kind": "HUC-8"}, {"id": "SEW", "kind": "NWS office"}, {"id": "WA", "kind": "State"}], "flow": {"change_vs_1_day_ago": "+7.5%", "change_vs_3_days_ago": "+261.9%", "change_vs_7_days_ago": "+182.6%", "percentile_for_the_date": 85.7, "percentile_of_all_prior_days": 88.9, "today": "5,030 cfs"}, "layers": {"nws-alert": "no flood-type warning over the gauge", "water": "streamflow above normal for the date"}, "place": {"drainage_area": "375 sq mi", "gauge": "USGS 12144500", "name": "SNOQUALMIE RIVER NEAR SNOQUALMIE, WA"}, "upstream": "no paired upstream gauge", "warnings": {"flood_type_warnings_by_office_last_24_hours": 0, "flood_type_warnings_by_office_last_7_days": 0, "flood_type_warnings_over_gauge_last_30_days": 0}}}} diff --git a/apps/atlas-outcomes/sources/flows.jsonl.gz b/apps/atlas-outcomes/sources/flows.jsonl.gz new file mode 100644 index 00000000..e0b75e5f Binary files /dev/null and b/apps/atlas-outcomes/sources/flows.jsonl.gz differ diff --git a/apps/atlas-outcomes/sources/gauges.jsonl.gz b/apps/atlas-outcomes/sources/gauges.jsonl.gz new file mode 100644 index 00000000..e26f2ddc Binary files /dev/null and b/apps/atlas-outcomes/sources/gauges.jsonl.gz differ diff --git a/apps/atlas-outcomes/sources/harvest_stats.json b/apps/atlas-outcomes/sources/harvest_stats.json new file mode 100644 index 00000000..a312ed94 --- /dev/null +++ b/apps/atlas-outcomes/sources/harvest_stats.json @@ -0,0 +1,49 @@ +{ + "end": "2026-09-22", + "endpoints": [ + "https://waterservices.usgs.gov/nwis/site/?format=rdb&siteOutput=expanded&sites={sites}", + "https://waterservices.usgs.gov/nwis/dv/?format=json&sites={site}¶meterCd=00060&statCd=00003&startDT={start}&endDT={end}", + "https://mesonet.agron.iastate.edu/api/1/nws/ugcs.json?state={state}", + "https://mesonet.agron.iastate.edu/api/1/vtec/sbw_interval.geojson?begints={b}T00:00Z&endts={e}T00:00Z&wfo={wfo}&only_new=true&{ph}" + ], + "failed": [], + "flow_days": 649341, + "flow_start": "1995-10-01", + "gauges": 60, + "harvested_at": "2026-09-23T23:20:12Z", + "paired": 25, + "warn_start_year": 2015, + "warnings": 37383, + "warnings_by_phenomena": { + "FA": 7112, + "FF": 13804, + "FL": 16467 + }, + "wfos": [ + "BGM", + "BOX", + "BTV", + "CRP", + "CTP", + "DVN", + "EAX", + "EWX", + "FFC", + "FWD", + "HGX", + "ILM", + "ILN", + "ILX", + "LMK", + "LSX", + "LWX", + "MHX", + "MTR", + "PAH", + "PBZ", + "PHI", + "RAH", + "RLX", + "SEW" + ] +} diff --git a/apps/atlas-outcomes/sources/warnings.jsonl.gz b/apps/atlas-outcomes/sources/warnings.jsonl.gz new file mode 100644 index 00000000..34ef339b Binary files /dev/null and b/apps/atlas-outcomes/sources/warnings.jsonl.gz differ diff --git a/apps/atlas-outcomes/tests/test_outcomes.py b/apps/atlas-outcomes/tests/test_outcomes.py new file mode 100644 index 00000000..f00ac9b3 --- /dev/null +++ b/apps/atlas-outcomes/tests/test_outcomes.py @@ -0,0 +1,284 @@ +"""Offline tests for atlas-outcomes: no network, a small synthetic two-gauge corpus. + + python3 -m unittest discover -s apps/atlas-outcomes/tests -v + +The fixture's flows and warnings are made up here to exercise the code; they +never reach sources/ or a pack. +""" + +from __future__ import annotations + +import copy +import json +import math +import random +import sys +import unittest +from datetime import UTC, date, datetime, timedelta +from pathlib import Path + +APP = Path(__file__).resolve().parent.parent +sys.path.insert(0, str(APP)) + +from atlas_outcomes import SCHEMA_ID +from atlas_outcomes.baseline import auc, fit_logistic, noul_metrics, predict +from atlas_outcomes.curate import ( + assemble, + assign_split, + is_easy_negative, + load_validator, + split_dates, +) +from atlas_outcomes.features import ( + DAY_S, + Corpus, + build_rows, + change_bucket, + decision_ts, + doy_index, + flow_class, + site_history, +) +from atlas_outcomes.geo import polygons_contain +from atlas_outcomes.harvest import dense, normalize_warning, parse_dv +from atlas_outcomes.records import render_state, to_record + +START = date(2008, 1, 1) +END = date(2015, 6, 30) +ROW_START = date(2015, 1, 1) + + +def _iso(ts: int) -> str: + return datetime.fromtimestamp(ts, UTC).strftime("%Y-%m-%dT%H:%M:%SZ") + + +def synthetic_flow(seed: int, n: int) -> list[float | None]: + rng = random.Random(seed) + out: list[float | None] = [] + for i in range(n): + season = 1.0 + 0.6 * math.sin(2 * math.pi * i / 365.25) + out.append(round(1000 * season * math.exp(rng.gauss(0, 0.35)), 1)) + out[100] = None # a gap the service did not report + return out + + +def square(lon: float, lat: float, r: float = 0.05) -> list: + return [[[[lon - r, lat - r], [lon + r, lat - r], [lon + r, lat + r], [lon - r, lat + r], [lon - r, lat - r]]]] + + +GAUGE_A = {"site": "00000001", "name": "TEST RIVER AT UPPER", "lat": 39.0, "lon": -77.0, "state": "MD", + "county_ugc": "MDC001", "county": "Allegany", "wfo": "LWX", "huc8": "02070002", "drainage_sq_mi": 500.0, + "tz": "EST", "upstream": None} +GAUGE_B = {**GAUGE_A, "site": "00000002", "name": "TEST RIVER AT LOWER", "lat": 38.5, "lon": -76.5, "upstream": "00000001"} + + +def fixture_corpus() -> Corpus: + n = (END - START).days + 1 + flows = [{"site": g["site"], "start": START.isoformat(), "cfs": synthetic_flow(k, n)} for k, g in enumerate((GAUGE_A, GAUGE_B))] + t = decision_ts(date(2015, 3, 10), "EST") + warnings = [ + # over gauge B, issued one hour after the end of 2015-03-10, in effect 12 h + {"wfo": "LWX", "phenomena": "FL", "significance": "W", "eventid": 1, "year": 2015, "issue": _iso(t + 3600), + "expire": _iso(t + 13 * 3600), "product_id": "p1", "ugcs": ["MDC001"], "polygons": square(-76.5, 38.5)}, + # same office, far from both gauges: office activity, never coverage + {"wfo": "LWX", "phenomena": "FA", "significance": "W", "eventid": 2, "year": 2015, "issue": _iso(t - 2 * DAY_S), + "expire": _iso(t - DAY_S), "product_id": "p2", "ugcs": ["VAC001"], "polygons": square(-80.0, 37.0)}, + ] + return Corpus([GAUGE_A, GAUGE_B], flows, warnings, (END + timedelta(days=2)).isoformat()) + + +class TestNoLeak(unittest.TestCase): + def test_history_ignores_the_future(self): + n = 9 * 365 + cfs = synthetic_flow(3, n) + h = site_history(START, cfs) + k = 8 * 365 + future = cfs[: k + 1] + [x * 50 if x is not None else None for x in cfs[k + 1 :]] + h2 = site_history(START, future) + self.assertEqual(h.pct_date[: k + 1], h2.pct_date[: k + 1]) + self.assertEqual(h.pct_all[: k + 1], h2.pct_all[: k + 1]) + self.assertEqual(h.p90[: k + 1], h2.p90[: k + 1]) + self.assertIsNotNone(h.pct_date[k]) + + def test_p90_is_strictly_before_t(self): + n = 6 * 365 + cfs = [100.0] * n + k = n - 10 + cfs[k] = 1e6 # a record flood on day k cannot raise its own threshold + h = site_history(START, cfs) + self.assertEqual(h.p90[k], 100.0) + self.assertEqual(h.pct_all[k], 100.0) + + def test_rows_ignore_future_flows_and_warnings(self): + corpus = fixture_corpus() + rows = {(r["site"], r["date"]): r for r in build_rows(corpus, ROW_START)} + cut = date(2015, 3, 5) + later = copy.deepcopy(corpus) + for f in later.flows: + i0 = (cut - START).days + 1 + f["cfs"] = f["cfs"][:i0] + [None if x is None else x * 7 for x in f["cfs"][i0:]] + t_cut = decision_ts(cut, "EST") + later.warnings.append({"wfo": "LWX", "phenomena": "FF", "significance": "W", "eventid": 9, "year": 2015, + "issue": _iso(t_cut + 60), "expire": _iso(t_cut + DAY_S), "product_id": "p9", + "ugcs": [], "polygons": square(-77.0, 39.0)}) + rows2 = {(r["site"], r["date"]): r for r in build_rows(later, ROW_START)} + checked = 0 + for key, r in rows.items(): + if date.fromisoformat(key[1]) <= cut: + self.assertEqual(r["features"], rows2[key]["features"], key) + checked += 1 + self.assertGreater(checked, 50) + # ...while the label of the cut day itself does see the future + self.assertTrue(rows2[("00000001", cut.isoformat())]["labels"]["warn_next"]) + self.assertFalse(rows[("00000001", cut.isoformat())]["labels"]["warn_next"]) + + def test_state_never_carries_the_answer(self): + corpus = fixture_corpus() + rows = build_rows(corpus, ROW_START) + for r in rows[:40]: + flipped = copy.deepcopy(r) + flipped["labels"] = {"high_next": not r["labels"]["high_next"], "warn_next": not r["labels"]["warn_next"], + "change_bucket": (r["labels"]["change_bucket"] + 2) % 5} + flipped["next_flow"] = r["next_flow"] * 3 + 1 + self.assertEqual(render_state(GAUGE_A, r), render_state(GAUGE_A, flipped)) + text = json.dumps(render_state(GAUGE_A, r)) + for word in ("high_next", "warn_next", "next_flow", "tomorrow"): + self.assertNotIn(word, text) + + +class TestLabels(unittest.TestCase): + def setUp(self): + self.rows = {(r["site"], r["date"]): r for r in build_rows(fixture_corpus(), ROW_START)} + + def test_warning_label_window(self): + b = "00000002" + self.assertTrue(self.rows[(b, "2015-03-10")]["labels"]["warn_next"]) # issued T + 1 h + self.assertFalse(self.rows[(b, "2015-03-09")]["labels"]["warn_next"]) # T(03-09) + 24 h = T(03-10) < issue + self.assertFalse(self.rows[(b, "2015-03-11")]["labels"]["warn_next"]) + self.assertFalse(self.rows[("00000001", "2015-03-10")]["labels"]["warn_next"]) # polygon is not over gauge A + f = self.rows[(b, "2015-03-11")]["features"] + self.assertEqual(f["warn_gauge_30d"], 1) + self.assertFalse(f["warn_active"]) # expired 12 h after issue, before T(03-11) + self.assertEqual(self.rows[(b, "2015-03-10")]["features"]["warn_wfo_7d"], 1) # the far FA warning only + self.assertEqual(self.rows[(b, "2015-03-11")]["features"]["warn_wfo_7d"], 2) + + def test_high_flow_label_is_next_day_over_stated_threshold(self): + corpus = fixture_corpus() + rows = self.rows + for (site, d), r in rows.items(): + f = next(x for x in corpus.flows if x["site"] == site) + nxt = f["cfs"][(date.fromisoformat(d) - START).days + 1] + self.assertEqual(r["next_flow"], nxt) + self.assertEqual(r["labels"]["high_next"], nxt > r["features"]["p90"]) + self.assertEqual(r["labels"]["change_bucket"], change_bucket(r["features"]["flow"], nxt)) + self.assertTrue(any(r["labels"]["high_next"] for r in rows.values())) + + def test_change_buckets_and_classes(self): + self.assertEqual([change_bucket(100, x) for x in (70, 80, 90, 95, 100, 105, 110, 120, 130)], [0, 1, 1, 2, 2, 2, 3, 3, 4]) + self.assertEqual(change_bucket(0, 0), 2) + self.assertEqual(change_bucket(0, 5), 4) + self.assertEqual([flow_class(p) for p in (5, 20, 50, 80, 95, None)], + ["much below normal", "below normal", "normal", "above normal", "much above normal", "not rated"]) + self.assertEqual(doy_index(date(2016, 3, 1)), doy_index(date(2015, 3, 1))) + + def test_upstream_features_come_from_the_upstream_gauge(self): + r = self.rows[("00000002", "2015-02-01")] + self.assertEqual(r["features"]["up_site"], "00000001") + self.assertEqual(r["features"]["up_pct_date"], self.rows[("00000001", "2015-02-01")]["features"]["pct_date"]) + self.assertNotIn("up_site", self.rows[("00000001", "2015-02-01")]["features"]) + + +class TestSplitAndRecords(unittest.TestCase): + def test_time_split_has_no_straddle(self): + first_hold, last_train = split_dates(date(2026, 9, 21)) + self.assertEqual(first_hold, date(2025, 9, 22)) + self.assertEqual(last_train, date(2025, 9, 20)) + self.assertIsNone(assign_split("2025-09-21", first_hold, last_train)) # embargo + self.assertEqual(assign_split("2025-09-22", first_hold, last_train), "holdout") + self.assertEqual(assign_split("2025-09-20", first_hold, last_train), "train") + + def test_assemble_validates_splits_and_downsamples_train_only(self): + kept, info = assemble(fixture_corpus(), easy_keep=0.25, row_start=ROW_START, holdout_days=60) + io = load_validator() + by_day: dict[str, set[str]] = {} + for v in kept: + self.assertEqual(io.validate_record(v["record"]), v["record"]) + self.assertEqual(set(v["record"]), {"schema", "id", "state", "questions", "labels"}) + self.assertEqual(v["record"]["schema"], SCHEMA_ID) + by_day.setdefault(v["meta"]["date"], set()).add(v["meta"]["split"]) + if v["meta"]["split"] == "holdout" or not v["meta"]["easy_negative"]: + self.assertEqual(v["meta"]["sample_weight"], 1.0) + else: + self.assertEqual(v["meta"]["sample_weight"], 4.0) + self.assertEqual(v["meta"]["provenance"], "outcome-real") + self.assertTrue(all(len(s) == 1 for s in by_day.values())) + last_train = max(d for d, s in by_day.items() if "train" in s) + first_hold = min(d for d, s in by_day.items() if "holdout" in s) + self.assertGreaterEqual((date.fromisoformat(first_hold) - date.fromisoformat(last_train)).days, 2) + self.assertIn("easy negative downsampled (train)", info["rejected"]) + hold_easy = sum(1 for v in kept if v["meta"]["split"] == "holdout" and v["meta"]["easy_negative"]) + self.assertEqual(hold_easy, info["easy_negatives"].get("holdout", 0)) + + def test_easy_negative_rule(self): + row = {"features": {"pct_date": 40.0, "pct_all": 50.0, "warn_active": False, "warn_gauge_30d": 0}, + "labels": {"high_next": False, "warn_next": False}} + self.assertTrue(is_easy_negative(row)) + for patch in ({"labels": {"high_next": True, "warn_next": False}}, + {"features": {**row["features"], "pct_date": 95.0}}, + {"features": {**row["features"], "warn_gauge_30d": 1}}): + self.assertFalse(is_easy_negative({**row, **patch})) + + def test_record_shape(self): + r = build_rows(fixture_corpus(), ROW_START)[0] + rec = load_validator().validate_record(to_record(GAUGE_A, r)) + self.assertEqual(rec["questions"]["high_next"]["type"], "noul") + self.assertEqual(rec["questions"]["flow_change"]["type"], "score") + self.assertEqual(len(rec["questions"]["flow_change"]["criteria"]), 5) + self.assertIn(f"{r['features']['p90']:,.0f} cfs", rec["questions"]["high_next"]["instructions"]) + self.assertEqual(rec["state"]["constructs"][0]["kind"], "County") + + +class TestParsersAndBaseline(unittest.TestCase): + def test_parse_dv_and_dense(self): + payload = {"value": {"timeSeries": [{"variable": {"noDataValue": -999999.0}, "values": [{"value": [ + {"value": "10", "dateTime": "2020-01-01T00:00:00.000"}, + {"value": "-999999", "dateTime": "2020-01-02T00:00:00.000"}, + {"value": "12.5", "dateTime": "2020-01-04T00:00:00.000"}]}]}]}} + vals = parse_dv(payload) + self.assertEqual(vals, {"2020-01-01": 10.0, "2020-01-02": None, "2020-01-04": 12.5}) + self.assertEqual(dense(vals), ("2020-01-01", [10, None, None, 12.5])) + + def test_normalize_warning_keeps_flood_warnings_only(self): + feat = {"properties": {"phenomena": "FL", "significance": "W", "eventid": 3, "year": 2024, "wfo": "LWX", + "utc_issue": "2024-01-09T15:20:00Z", "utc_expire": "2024-01-10T03:50:00Z", + "ugclist": "MDC021, VAC107", "product_id": "x"}, + "geometry": {"type": "Polygon", "coordinates": [[[-77.123456, 39.1], [-77.0, 39.1], [-77.0, 39.2], [-77.123456, 39.1]]]}} + w = normalize_warning(feat) + self.assertEqual(w["ugcs"], ["MDC021", "VAC107"]) + self.assertEqual(w["polygons"][0][0][0], [-77.1235, 39.1]) + self.assertIsNone(normalize_warning({**feat, "properties": {**feat["properties"], "significance": "A"}})) + self.assertIsNone(normalize_warning({**feat, "properties": {**feat["properties"], "phenomena": "SV"}})) + + def test_point_in_polygon_with_hole(self): + outer = [[0, 0], [10, 0], [10, 10], [0, 10], [0, 0]] + hole = [[4, 4], [6, 4], [6, 6], [4, 6], [4, 4]] + polys = [[outer, hole]] + self.assertTrue(polygons_contain(polys, 2, 2)) + self.assertFalse(polygons_contain(polys, 5, 5)) + self.assertFalse(polygons_contain(polys, 11, 5)) + + def test_auc_and_logistic(self): + self.assertEqual(auc([0.1, 0.4, 0.35, 0.8], [0, 0, 1, 1]), 0.75) + self.assertEqual(auc([0.5, 0.5], [0, 1]), 0.5) + rng = random.Random(0) + xs = [rng.gauss(0, 1) for _ in range(400)] + ys = [1.0 if x + rng.gauss(0, 0.5) > 0 else 0.0 for x in xs] + cols = [[1.0] * 400, xs] + beta = fit_logistic(cols, ys, [1.0] * 400) + self.assertGreater(beta[1], 1.0) + m = noul_metrics(predict(beta, cols), [int(y) for y in ys]) + self.assertGreater(m["auc"], 0.9) + + +if __name__ == "__main__": + unittest.main() diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index 780905df..e1dd5668 100644 --- a/docs/ARCHITECTURE.md +++ b/docs/ARCHITECTURE.md @@ -129,6 +129,7 @@ plan. All three run in the threadpool, never on the event loop. | Code graph hits | personal-graphify `graph.json` | opt-in (`DOTTIE_CONTEXT_GRAPH=`) | off by default | | Goal features | computed | **yes** | lexical features, hash; always | | dottie_loop `MemoryStore`, scout brain `MEMORY.md`, ava-skills memory | their own files | **no** | not decision inputs; any future use comes through a `ContextProvider` adapter, not a new store | +| Place state (USGS gauge, Atlas construct stack, NWS flood warnings) | `apps/atlas-outcomes/sources/` snapshots | **no** | an `outcome-real` System One decision pack for place decisions; a future place `ContextProvider` would serve the same state | | harness-api analytics / corpus files | `apps/dottie-harness-api/lib` | **no** | dashboard only | Every provider runs under a per-provider timeout (default 250 ms) and fails @@ -157,6 +158,18 @@ source, id; the digest hashes exactly those items). 6. **Serve** from dottie-os `/decide`. The served checkpoint's identity hash in `/health` is what the stamp is checked against. +**Place decisions (`apps/atlas-outcomes`).** A second source of real labels, +outside the router: System One records about a place on the eye.jcamd.com +Atlas (a USGS gauge, its construct stack and the NWS warnings over it) at the +end of day t, labelled by what happened next (flow on t+1 above the gauge's +p90, an NWS flood-type warning polygon over it within 24 h, the next-day +change level). Provenance tier `outcome-real`: the rows may train System One +candidates and are scored on a time-split holdout (the latest 365 days) +against the persistence and climatology baselines in +`apps/atlas-outcomes/BASELINE.json`; the same stamp rule applies. The same +state shape is what a future place `ContextProvider` would hand the decision +plane; none is wired today. + **Today:** the runner's executors are deterministic stubs except the MCP operator, so almost every outcome is `executor: stub` and no pack can be built from real labels yet. That is the honest state; Phase 2 adds real executors. @@ -203,7 +216,9 @@ tested with a fake model); its GPU latency is not measured here. says `gate_passed: true` and a human stamped its exact bytes. Today nothing is stamped: the heuristic is the authority and every learned or System One answer is advisory, logged and displayed. -- Synthetic, test or stub-executor rows never train a champion. +- Synthetic, teacher, test or stub-executor rows never train a champion. + Rows whose labels are recorded futures (provenance `outcome-real`) may + train candidates; they are evaluated on a time-split holdout. - No new "confidence" in code or copy. System One reports `shape_concentration`; external probabilities are `backend_confidence`. The router's existing `confidence` key (a keyword-score ratio) keeps its name. diff --git a/docs/FACTORY.md b/docs/FACTORY.md index b7635bbe..c5927662 100644 --- a/docs/FACTORY.md +++ b/docs/FACTORY.md @@ -113,7 +113,7 @@ runs, not something this repo installs. ```json {"datasets": [{ "id": "arxiviq-papers", "repo": "arxiviq", "path": "site/public/data/papers.json", - "provenance": "real", // real | honest-synthetic | placeholder | unknown + "provenance": "real", // real | outcome-real | honest-synthetic | placeholder | unknown "source": "arXiv API via scripts/fetch_topics.py", "refresh": "python scripts/fetch_topics.py", "cadence_days": 7, "fresh_key": "json:generated_at", // or "mtime" @@ -128,6 +128,11 @@ runs, not something this repo installs. | `data refresh ID` | runs the refresh command in the owning repo, then re-checks | | `data restore ID` | copies the first existing `restore_from` source into `path` and writes `.manifest.json` (source, sha256, size, when); refuses to overwrite a present file without `--force` | +`outcome-real` marks a decision pack whose labels are recorded futures (what +happened after the state was observed, e.g. `apps/atlas-outcomes`). It may +train candidates, is evaluated on a time-split holdout, and never promotes +itself. + Freshness is declared, not inferred: no `cadence_days` means the dataset is static and can only be missing, never stale. diff --git a/factory/config.py b/factory/config.py index bec3b44e..0c0ce7d0 100644 --- a/factory/config.py +++ b/factory/config.py @@ -25,7 +25,10 @@ DATASETS_PATH = HERE / "datasets.json" RUNS_DIR = HERE / "runs" -PROVENANCE = {"real", "honest-synthetic", "placeholder", "unknown"} +# outcome-real: labels are recorded futures (what happened after the state was +# observed). Such rows may train candidates and are evaluated on a time-split +# holdout; teacher/synthetic rows never train a champion (docs/ARCHITECTURE.md). +PROVENANCE = {"real", "outcome-real", "honest-synthetic", "placeholder", "unknown"} ROLES = {"center", "game", "site", "service", "library", "archived"} GATE_OPS = {">=", "<=", ">", "<"} diff --git a/factory/datasets.json b/factory/datasets.json index 60938bff..aadf9120 100644 --- a/factory/datasets.json +++ b/factory/datasets.json @@ -317,6 +317,45 @@ "required": false, "restore_from": [], "consumers": [] + }, + { + "id": "atlas-outcomes-flows", + "repo": "dottie", + "path": "apps/atlas-outcomes/sources/flows.jsonl.gz", + "provenance": "real", + "source": "USGS NWIS daily values (mean daily discharge, 00060) for 60 gauges, 1995-10-01 to the harvest end; gauge metadata in sources/gauges.jsonl.gz (USGS site service + IEM county UGC -> NWS office)", + "refresh": "{python} apps/atlas-outcomes/run.py harvest --end ", + "cadence_days": null, + "fresh_key": "mtime", + "required": false, + "restore_from": [], + "consumers": [] + }, + { + "id": "atlas-outcomes-warnings", + "repo": "dottie", + "path": "apps/atlas-outcomes/sources/warnings.jsonl.gz", + "provenance": "real", + "source": "IEM VTEC storm-based warning polygons (FL, FA, FF warnings at issuance) for every office with a gauge, 2015 to the harvest end", + "refresh": "{python} apps/atlas-outcomes/run.py harvest --end ", + "cadence_days": null, + "fresh_key": "mtime", + "required": false, + "restore_from": [], + "consumers": [] + }, + { + "id": "atlas-outcomes-pack", + "repo": "dottie", + "path": "apps/atlas-outcomes/data/packs/atlas-outcomes-1/train.jsonl", + "provenance": "outcome-real", + "source": "jev-decision-schema-1.0.0 place-decision records (gauge state and Atlas construct stack at the end of day t) labelled by recorded futures: flow on t+1 above p90, an NWS flood-type warning polygon over the gauge within 24 h, next-day change level; time-split holdout (latest 365 days); consent.champion=false", + "refresh": "{python} apps/atlas-outcomes/run.py curate", + "cadence_days": null, + "fresh_key": "mtime", + "required": false, + "restore_from": [], + "consumers": [] } ] }