Region-scoped observation pooling as bundle data: lnk_species_pooling() (#290)

CLAUDE.md

@@ -60,6 +60,31 @@ To compare two code versions, prepare once and re-run classify/connect on the
 same schema (`data-raw/logs/habitat_thresholds_282/reclassify.R`). Issue drafted,
 awaiting body review.
 
+## Status (2026-09-26, late) — observation pooling is data (#290); #284 step 5 unparked
+
+**Which observation species count as which model species is a bundle tracker now, not code.**
+- `configs/default/species_pooling.csv` holds one dated, sourced row per decision, scoped to a region, sub-region or WSG.
+- `inst/extdata/wsg_regions.csv` is package-level geography: the region is the first segment of each group's *outlet* wscode.
+- `species_groups.csv` lets either side of a row be a group.
+- `lnk_species_pooling()` resolves it and knows no species.
+- **DV → BT pools in the Fraser, Mackenzie and Columbia/Kootenay, and in the Skeena only above Hazelton** (the sub-region `Skeena above Hazelton`: BULK, MORR, KISP, BABL, BABR, SUST, MSKE, USKE). Everything unlisted is not pooled.
+  - **Why the split:** interior DV is a legacy name (the DV share of char records is ~90 % before 1990 and under 15 % after 2000). Skeena DV is not: 86 % or more in every decade, and 85 % of Skeena DV streams were re-sampled and still recorded DV. The Skeena split was the operator's call, 2026-09-27.
+  - **The optional `obs_year_max` column** limits a row to records dated in or before a year; undated records do not qualify. It is empty in the seed.
+- The state of knowledge lives in `research/species_pooling.md`, with the scenarios from `data-raw/species_pooling_evidence.R`.
+- **Do not re-hard-code a species pair anywhere**; add a row.
+- **`default_tuned`'s BT `rear_gradient_max` of 0.1349 hangs by a thread.**
+  - Under the split, the pooled P95 is 0.1309, which clears the 0.13 line that would give 0.1249 by 0.0009.
+  - With no Skeena pooling, or with interior pooling limited to pre-1995 records, it gives 0.1249.
+  - #284 step 5 scoring should test both values.
+
+**Facts not worth re-deriving:**
+- **The region rule alone moved no #284 or #283 number** (all 55 + 59 WSGs are in pooled regions). **The Hazelton split does move them:** 852 DV records in LSKE, KLUM, LKEL and ZYMO drop out. No #284 verdict changes.
+- **Where the new rule bites:** 29 BT-present WSGs in the Nass, Stikine, Taku, Yukon and the coast, which pooled before and do not now.
+- **"No tracker" means no pooling in the resolver**, but the validator's list default (global DV→BT) is kept for back-compat. The driver resolves pooling **once**, from the first bundle declaring a tracker, and applies it to both bundles.
+- **Two presence tables must agree, WSG by WSG**: the pooling bundle's and each scored bundle's. Both drivers stop if they do not, and the comparison is keyed by WSG, never by row position (review found the row-position version).
+- **The #284 obs producer is not bit-reproducible.** A double-precision `sum(length_metre)` gives 1-ulp differences in `candidates.csv` between identical runs. This is pre-existing: #293.
+- **Hand-kept `consumed_by` file:line refs drift on every edit above them.** The two new dictionaries are now test-checked with a whole-token match.
+
 ## Status (2026-09-26) — v0.51.1: first calibrated threshold in `default_tuned` (#284; step 5 open)
 
 **`research/habitat_thresholds.md` is the staging table for tuning CH/BT numbers.** One

NAMESPACE

@@ -53,6 +53,7 @@ export(lnk_rollup_wsg)
 export(lnk_rules_build)
 export(lnk_score)
 export(lnk_source)
+export(lnk_species_pooling)
 export(lnk_stamp)
 export(lnk_stamp_finish)
 export(lnk_thresholds)

R/lnk_habitat_validate.R

@@ -135,10 +135,20 @@
 #'   matched on that key). Species and WSG codes are compared upper-cased
 #'   and trimmed. Optional: `match_type`, `source`, `is_spawn`, `is_rear`
 #'   (logical, 0/1 or t/true/yes), `activity_code`, `activity`,
-#'   `life_stage`.
-#' @param species_obs Named list mapping a model species to the observation
-#'   species codes that count as it. Species not named map to themselves.
-#'   Default pools DV records with BT (`list(BT = c("BT", "DV"))`).
+#'   `life_stage`, and `observation_date`, which is required when
+#'   `species_obs` carries a year limit (`obs_year_max`).
+#' @param species_obs Which observation species count as each model species.
+#'   Either a named list mapping a model species to observation species
+#'   codes, applied in every WSG (default `list(BT = c("BT", "DV"))`, which
+#'   pools DV records with BT), or a per-WSG data frame with
+#'   `watershed_group_code`, `species_code` and `obs_species`, such as
+#'   [lnk_species_pooling()] returns, optionally with `obs_year_max` (a
+#'   pooled record counts only if dated in or before that year; an undated
+#'   one does not). In the list form a species not named
+#'   maps to itself; in the data frame form a species always counts as
+#'   itself, and a WSG and species pair it does not list maps to itself only.
+#'   Either way a species is scored only where
+#'   `loaded$wsg_species_presence` marks it present.
 #' @param match_types Character vector of `match_type` classes (first
 #'   letter) to keep, or `NULL` for no match-type filter. Default
 #'   `c("A", "B")`.
@@ -216,8 +226,6 @@ lnk_habitat_validate <- function(conn, aoi, cfg, loaded, species, schema,
     is.character(schema), length(schema) == 1L, !is.na(schema),
     grepl("^[a-z_][a-z0-9_]*$", schema),
     is.list(species_obs),
-    length(species_obs) == 0L || !is.null(names(species_obs)),
-    all(vapply(species_obs, is.character, logical(1))),
     is.null(match_types) ||
       (is.character(match_types) && length(match_types) >= 1L &&
          all(grepl("^[A-Z]$", match_types))),
@@ -240,8 +248,7 @@ lnk_habitat_validate <- function(conn, aoi, cfg, loaded, species, schema,
   }
   aoi <- unique(aoi)
   species <- unique(toupper(species))
-  species_obs <- lapply(species_obs, toupper)
-  names(species_obs) <- toupper(names(species_obs))
+  species_obs <- .lnk_hv_species_obs(species_obs)
 
   .lnk_hv_check_schema(conn, schema, aoi, species)
   logged <- .lnk_hv_check_log(conn, schema, aoi, cfg)
@@ -355,26 +362,79 @@ lnk_habitat_validate <- function(conn, aoi, cfg, loaded, species, schema,
 #' the model species is present.
 #' @noRd
 .lnk_hv_spec <- function(presence, aoi, species, species_obs) {
+  per_wsg <- is.data.frame(species_obs)
   rows <- lapply(aoi, function(w) {
     row <- presence[presence$watershed_group_code == w, , drop = FALSE]
     if (nrow(row) == 0L) return(NULL)
     present <- intersect(species, .lnk_wsg_species_present(row[1, ]))
     if (length(present) == 0L) return(NULL)
     do.call(rbind, lapply(present, function(sp) {
-      data.frame(watershed_group_code = w, species_code = sp,
-                 obs_species = unique(species_obs[[sp]] %||% sp),
-                 stringsAsFactors = FALSE)
+      if (per_wsg) {
+        k <- species_obs$watershed_group_code == w &
+          species_obs$species_code == sp
+        obs <- c(sp, species_obs$obs_species[k])
+        yr <- c(NA_integer_, species_obs$obs_year_max[k])
+      } else {
+        obs <- species_obs[[sp]] %||% sp
+        yr <- rep(NA_integer_, length(obs))
+      }
+      d <- unique(data.frame(watershed_group_code = w, species_code = sp,
+                             obs_species = obs, obs_year_max = yr,
+                             stringsAsFactors = FALSE))
+      if (anyDuplicated(d$obs_species)) {
+        stop("species_obs gives ", w, " ", sp, " two year limits for one ",
+             "observation species", call. = FALSE)
+      }
+      d
     }))
   })
   out <- do.call(rbind, rows)
   if (is.null(out)) {
     out <- data.frame(watershed_group_code = character(0),
                       species_code = character(0),
-                      obs_species = character(0), stringsAsFactors = FALSE)
+                      obs_species = character(0),
+                      obs_year_max = integer(0), stringsAsFactors = FALSE)
   }
   out
 }
 
+#' Normalise `species_obs`: a named list (applied in every WSG) or a per-WSG
+#' data frame (`watershed_group_code`, `species_code`, `obs_species`).
+#' Checked as a data frame first, because a data frame is also a list and
+#' would otherwise pass the list checks and silently map nothing.
+#' @noRd
+.lnk_hv_species_obs <- function(species_obs) {
+  up <- function(x) toupper(trimws(as.character(x)))
+  if (is.data.frame(species_obs)) {
+    cols <- c("watershed_group_code", "species_code", "obs_species")
+    miss <- setdiff(cols, names(species_obs))
+    if (length(miss) > 0L) {
+      stop("species_obs data frame lacks columns: ",
+           paste(miss, collapse = ", "), call. = FALSE)
+    }
+    out <- data.frame(lapply(species_obs[cols], up), stringsAsFactors = FALSE)
+    if (anyNA(out) || !all(vapply(out, function(x) all(nzchar(x)), TRUE))) {
+      stop("species_obs data frame has empty or NA codes", call. = FALSE)
+    }
+    # Optional record-date limit from lnk_species_pooling(); NA is no limit.
+    out$obs_year_max <- if ("obs_year_max" %in% names(species_obs)) {
+      .lnk_sp_year(species_obs[["obs_year_max"]])
+    } else {
+      rep(NA_integer_, nrow(out))
+    }
+    return(unique(out))
+  }
+  if (!(length(species_obs) == 0L || !is.null(names(species_obs))) ||
+      !all(vapply(species_obs, is.character, logical(1)))) {
+    stop("species_obs must be a named list of character vectors or a data ",
+         "frame of watershed_group_code, species_code, obs_species",
+         call. = FALSE)
+  }
+  out <- lapply(species_obs, up)
+  names(out) <- up(names(out))
+  out
+}
+
 #' Spawn / rear stage from bcfishobs activity and life stage.
 #'
 #' Spawning from activity; rearing from activity or a juvenile life stage.
@@ -424,7 +484,8 @@ lnk_habitat_validate <- function(conn, aoi, cfg, loaded, species, schema,
   cols_req <- c("species_code", "watershed_group_code", "blue_line_key",
                 "downstream_route_measure")
   cols_opt <- c("observation_key", "match_type", "source", "is_spawn",
-                "is_rear", "activity_code", "activity", "life_stage")
+                "is_rear", "activity_code", "activity", "life_stage",
+                "observation_date")
   key_generated <- FALSE
   if (is.data.frame(observations)) {
     d <- as.data.frame(observations)
@@ -506,6 +567,10 @@ lnk_habitat_validate <- function(conn, aoi, cfg, loaded, species, schema,
   DBI::dbWriteTable(conn, "lnk_vd_excl",
                     data.frame(observation_key = as.character(keys)),
                     temporary = TRUE, overwrite = TRUE)
+  if (any(!is.na(spec$obs_year_max)) && !has("observation_date")) {
+    stop(s$src, " has no observation_date column, which the year limits in ",
+         "species_obs (obs_year_max) need", call. = FALSE)
+  }
   DBI::dbWriteTable(conn, "lnk_vd_spec", spec,
                     temporary = TRUE, overwrite = TRUE)
   uhc <- loaded$user_habitat_classification
@@ -541,13 +606,16 @@ lnk_habitat_validate <- function(conn, aoi, cfg, loaded, species, schema,
          JOIN pg_temp.lnk_vd_spec sp
            ON sp.watershed_group_code = upper(trim(o.watershed_group_code::text))
           AND sp.obs_species = upper(trim(o.species_code::text))
+          -- a year-limited pooling admits only records dated within it
+          AND (sp.obs_year_max IS NULL
+               OR extract(year FROM %10$s) <= sp.obs_year_max)
         WHERE TRUE %8$s %9$s) o
       WHERE NOT EXISTS (SELECT 1 FROM pg_temp.lnk_vd_excl e
                          WHERE e.observation_key = o.observation_key)",
     s$src, opt("match_type", "text"), opt("activity_code", "text"),
     opt("activity", "text"), opt("life_stage", "text"),
     opt_flag("is_spawn"), opt_flag("is_rear"),
-    where_match, where_source),
+    where_match, where_source, opt("observation_date", "date")),
     params = if (is.null(match_types)) NULL else
       list(paste0("{", paste(match_types, collapse = ","), "}")))
 

R/lnk_log.R

@@ -98,6 +98,20 @@
     digest::digest(file = p, algo = "sha256")
   }, character(1), USE.NAMES = FALSE)
 
+  # A pooling tracker (#290) is scoped by the package's region lookup, which
+  # sits outside the bundle: a region reassignment changes what pools, so it
+  # enters the hash too, by a fixed name.
+  if (!is.null(cfg$files$species_pooling)) {
+    regions <- .lnk_wsg_regions_path()
+    paths <- c(paths, regions)
+    rel <- c(rel, "link:wsg_regions.csv")
+    digests <- c(digests, if (file.exists(regions)) {
+      digest::digest(file = regions, algo = "sha256")
+    } else {
+      "MISSING"
+    })
+  }
+
   if (length(fallback) == 1L && nzchar(fallback)) {
     paths <- c(paths, fallback)
     rel <- c(rel, "fresh:parameters_habitat_thresholds.csv")

R/lnk_presence.R

@@ -91,8 +91,7 @@ lnk_presence <- function(
          call. = FALSE)
   }
 
-  species_cols <- setdiff(names(wsg_species_presence),
-                          c("watershed_group_code", "notes"))
+  species_cols <- .lnk_presence_species_cols(wsg_species_presence)
 
   raw_present <- vapply(species_cols, function(sp) {
     .lnk_presence_truthy(row[[sp]])