---
title: "Build v7.1 (`v7b`) — sea-turtle core habitat + the coverage floor"
params:
ver: v7b
ver_prev: v7
msens:
target_name: build_v7b
workflow_type: release
dependency: []
output: data/manifests/build_v7b.json
editor_options:
chunk_output_type: console
---
v7.1 is a **surgical patch of the public v7 release** — same `usa05` grid, same `mdl_seq`/`value`
schema, same 20 Program Areas — built the way v7 itself was built: clone the previous DuckDB, change
only what the fix touches, then re-execute v7's own scoring chunks verbatim. Two method changes,
nothing else:
1. **Sea-turtle core habitat.** AquaMaps suitability enters the turtle merge only where P ≥ 0.5
(AquaMaps' own core-habitat default), and a range cell with no retained suitability is
**dropped** — the floor-of-1 fill goes. A threshold alone removes none of that fill: for five of
the six turtles half to nine tenths of the merged US area *is* fill, with no AquaMaps value at
all. Critical habitat still overrides by `max`, so no designated cell is ever trimmed.
2. **Coverage floor κ = 5 %.** A component is reportable in a Program Area only if its scored cells
cover ≥ 5 % of the area. Below that it is *not reportable* — no row, excluded from the composite
exactly as a component with no cells is, never published as a zero.
Both are flags whose **control value is off**, so the same code path reproduces v7 exactly. That is
the first gate below, and it is what makes every other number here attributable.
Not in v7.1: the year-integrated productivity, the distribution-share / national score, occurrence
masks, per-taxon thresholds, any change to Formula 9.4. All of that is v10. **Every `primprod*` row
of v7b must equal v7's**, and a gate asserts it in every run.
## Design
```{mermaid}
%%| label: fig-design
%%| fig-cap: "v7 sdm.duckdb → copy-on-write clone → turtle rewrite under the SAME mdl_seq → v7's own scoring chunks → κ → v7b"
flowchart LR
v7[("v7 sdm.duckdb<br/>read-only")] -->|cp -c| cl[("{ver} sdm.duckdb")]
ts["turtle_src<br/>taxon x taxon_model x model_cell"] --> tr["msens::turtle_sql(<br/>suit_min, fill,<br/>half_even, ch_outer)"]
cl --> ts
tr -->|"DELETE + INSERT<br/>mdl_seq 54238-54243"| cl
cl --> sc["v7 scoring chunks, verbatim<br/>(old/calc_scores_v7.qmd)"]
sc --> k["msens::coverage_sql + coverage_floor<br/>kappa = 5%"]
k --> sc2["composite score<br/>SUM(value)/COUNT(value)"]
sc2 --> out[("v7b: model_cell, cell_metric,<br/>zone_metric, zone_taxon,<br/>release_method")]
```
## Setup
```{r}
#| label: setup
librarian::shelf(
DBI, dbplyr, dplyr, duckdb, fs, glue, here, jsonlite, logger, purrr, readr,
sf, stringr, tibble, tidyr,
quiet = T)
options(readr.show_col_types = F)
ver <- params$ver
ver_prev <- params$ver_prev
# paths come from msens::sdm_db_path(), never libs/paths.R (that file pins v9). The
# small per-release artifacts live under derived/{ver} (which IS big/{ver} on the
# laptop but NOT on the server), because backfill_versions.qmd reads layers_{ver}.csv
# from there.
is_server <- Sys.info()[["sysname"]] == "Linux"
dir_derived <- if (is_server) "/share/data/derived" else "~/_big/msens/derived"
sdm_db <- path.expand(msens::sdm_db_path(ver))
sdm_db_prev <- path.expand(msens::sdm_db_path(ver_prev))
dir_big_v <- path_dir(sdm_db)
dir_v <- path.expand(glue("{dir_derived}/{ver}"))
dir_v_prev <- path.expand(glue("{dir_derived}/{ver_prev}"))
dir_scratch <- path.expand(glue("{dir_big_v}/_controls"))
dir_create(c(dir_big_v, dir_v))
stopifnot(
"previous release database not found" = file_exists(sdm_db_prev),
"previous release Program Area gpkg not found" =
file_exists(glue("{dir_v_prev}/ply_programareas_2026_{ver_prev}.gpkg")),
"msens 0.42.0 or newer is required (turtle_sql knobs + coverage_sql)" =
packageVersion("msens") >= "0.42.0",
"grid registry does not know this version" = msens::grid_for_ver(ver) == "usa05")
# method flags -- CONTROL VALUE IS OFF for each, so all-off == v7 exactly ----
# TURTLE_SUIT_MIN: suitability strictly below this is absent BEFORE the multiply (0 = off)
# TURTLE_FILL : suitability assumed for an ER cell with none; NA drops the cell (1 = off)
# COVERAGE_FLOOR : kappa, the reportability floor as a fraction (0 = off)
env_num <- function(key, default) {
v <- trimws(Sys.getenv(key))
if (!nzchar(v)) return(default)
if (tolower(v) %in% c("na", "none", "")) return(NA_real_)
as.numeric(v)
}
TURTLE_SUIT_MIN <- env_num("TURTLE_SUIT_MIN", 50)
TURTLE_FILL <- env_num("TURTLE_FILL", NA_real_)
COVERAGE_FLOOR <- env_num("COVERAGE_FLOOR", 0.05)
do_build <- !file_exists(sdm_db) || nzchar(Sys.getenv("REDO_BUILD"))
# the turtle merge inputs, exactly v7's (old/merge_models_v7.qmd:361-371): suitability
# = am_0.05, spatial ER = the SWOT+DPS range, CH = the three designation datasets.
# rng_iucn / rng_fws are NOT turtle inputs.
DS_ER <- "rng_turtle_swot_dps"
DS_SUIT <- "am_0.05"
DS_CH <- c("ch_nmfs", "ch_fws", "ca_nmfs")
# size DuckDB to the machine, not to the laptop this was written on: duckdb_budget() /
# duckdb_tune() (40 % of available memory, <= 12 GB, threads <= cores; the why is in the file)
source(here("libs/duckdb_budget.R"))
f_gates <- file.path(dir_v, glue("tbl_gates_{ver}.csv"))
f_attr <- file.path(dir_v, glue("tbl_attribution_programarea_{ver}.csv"))
log_info(
"{ver} <- {ver_prev}: suit_min={TURTLE_SUIT_MIN}, fill={TURTLE_FILL}, ",
"kappa={COVERAGE_FLOOR}; build={do_build}; db={sdm_db}")
```
## The build
The whole patch is one function of `(flags, target dir)`, so the three control runs of
@sec-controls execute **the same code** on copy-on-write clones in a scratch directory. Nothing
here reads `libs/paths.R` and nothing writes to `{ver_prev}`: the previous release is opened
`read_only = TRUE` and never otherwise.
### v7's scoring chunks, read out of the archived notebook
v7 was never rendered from these notebooks — it is v6's database patched by `dev/build_v7.R`, which
**evals labelled chunks of the scoring notebook verbatim**. That technique is kept (it is the only
way to be bit-faithful to v7's arithmetic) but moved inside a committed, rendered notebook.
```{r}
#| label: chunk_reader
qmd_v7 <- here("old/calc_scores_v7.qmd")
qmd_lines <- readLines(qmd_v7)
get_chunk <- function(label) {
lab <- grep(glue("^#\\|\\s*label:\\s*{label}\\s*$"), qmd_lines)
stopifnot("chunk label not found exactly once in old/calc_scores_v7.qmd" = length(lab) == 1)
open <- max(grep("^```\\{", qmd_lines[seq_len(lab)]))
close <- lab + grep("^```\\s*$", qmd_lines[(lab + 1):length(qmd_lines)])[1]
code <- qmd_lines[(open + 1):(close - 1)]
paste(code[!grepl("^#\\|", code)], collapse = "\n")
}
# the two halves of v7's scoring, with kappa between them: the composite is
# SUM(value)/COUNT(value) over the component rows, so a not-reportable component must be
# gone BEFORE calc_programarea_score_from_programarea_metrics runs.
CHUNKS_PRE <- c(
"calc_cell_metric_redlist", "calc_ecoregion_minmax",
"extrisk_spcat_ecoregion_rescaled_cell", "cell_metrics_to_zone_metrics",
"apply_pctarea_to_programarea_components")
CHUNKS_POST <- c(
"cell_metrics_to_ecoregion_metrics", "cell_metric_score",
"calc_programarea_score_from_programarea_metrics", "metric_labels",
"add_full_studyarea_zone", "zone_taxon")
# `IUCN:TN` is a v1-v7 data artifact: an ESA "Threatened" code stamped with the IUCN
# authority, on 11 taxa (1 mammal -- the Arctic ringed seal -- and 10 birds). The msens that
# built v7 did not validate the code, so it fell through `score_iucn`'s `.default = 1L` and
# scored **1**; that is what `taxon.er_score` still stores for all 11. msens >= 0.19 ERRORS
# on it, which stops `calc_cell_metric_redlist` dead. Reproducing v7 means reproducing that
# 1 -- recoding those taxa is a data question, and it belongs to v10, not to a patch that
# must reproduce its predecessor exactly. The IUCN branch also ignores `is_mmpa` (the same
# documented quirk that made `IUCN:LC + MMPA` score 1); that too is v7's behavior and stays.
ER_LEGACY <- c("IUCN:TN" = 1L)
er_score_v7 <- function(extrisk_code, is_mmpa = FALSE, is_mbta = FALSE) {
n <- length(extrisk_code)
out <- rep(NA_integer_, n)
leg <- !is.na(extrisk_code) & extrisk_code %in% names(ER_LEGACY)
if (any(leg)) out[leg] <- unname(ER_LEGACY[extrisk_code[leg]])
if (any(!leg)) out[!leg] <- msens::compute_er_score(
extrisk_code[!leg],
is_mmpa = rep_len(is_mmpa, n)[!leg],
is_mbta = rep_len(is_mbta, n)[!leg])
out
}
# dev/build_v7.R:87-88 widened zone_taxon from (subregions, program areas) to include
# ecoregions. The archived notebook already carries that edit, so the substitution is a
# no-op -- asserted rather than assumed, because a silent no-op would drop 45k rows.
chunk_fix <- function(code)
gsub("msens::compute_er_score(", "er_score_v7(",
gsub("c(tbl_sr, tbl_pra)", "c(tbl_sr, tbl_pra, tbl_er)", code, fixed = TRUE), fixed = TRUE)
zt_fix <- chunk_fix
stopifnot("zone_taxon chunk does not cover ecoregions" =
grepl("c(tbl_sr, tbl_pra, tbl_er)", chunk_fix(get_chunk("zone_taxon")), fixed = TRUE),
"the extinction-risk chunk no longer calls compute_er_score" =
grepl("er_score_v7(", chunk_fix(get_chunk("calc_cell_metric_redlist")), fixed = TRUE))
# no primprod_* chunk is evaluated (v7.1 decision 4: the year-integrated productivity
# stays with v10), so those cell metrics remain v7's own rows.
stopifnot("a primprod chunk leaked into the run list" =
!any(grepl("primprod", c(CHUNKS_PRE, CHUNKS_POST))))
# er_score_v7 may differ from msens ONLY on the legacy code, and only by not erroring
.codes <- c("NMFS:EN", "NMFS:TN", "NMFS:LC", "FWS:EN", "FWS:TN", "FWS:LC",
"IUCN:CR", "IUCN:EN", "IUCN:VU", "IUCN:NT", "IUCN:LC", "IUCN:DD", NA)
stopifnot(
"er_score_v7 diverges from msens::compute_er_score on a valid code" =
identical(er_score_v7(.codes, is_mmpa = TRUE, is_mbta = TRUE),
msens::compute_er_score(.codes, is_mmpa = TRUE, is_mbta = TRUE)),
"the legacy IUCN:TN code must score 1, as it did for v1-v7" =
identical(er_score_v7(c("IUCN:TN", "IUCN:LC"), is_mmpa = TRUE), c(1L, 1L)))
```
### Clone, turtles, κ, `release_method`
```{r}
#| label: build_fns
# copy-on-write clone: APFS `cp -c` on the laptop, `cp --reflink=auto` on Linux. A 3.6 GB
# release database clones in milliseconds and only the pages the patch rewrites cost space.
clone_db <- function(src, dst) {
if (file_exists(dst)) file_delete(dst)
if (file_exists(glue("{dst}.wal"))) file_delete(glue("{dst}.wal"))
dir_create(path_dir(dst))
flag <- if (is_server) "--reflink=auto" else "-c"
rc <- system2("cp", c(flag, shQuote(src), shQuote(dst)))
stopifnot("copy-on-write clone failed" = identical(as.integer(rc), 0L))
log_info("cloned {src} -> {dst}")
}
# the turtle rewrite: same rule, same ids. Re-merging through the old loop would RENUMBER
# mdl_seq and break every published ?mdl_seq= deep link, so rows are replaced under the
# existing merged ids.
rewrite_turtles <- function(con, suit_min, fill) {
seqs <- dbGetQuery(con, "SELECT mdl_seq FROM taxon WHERE sp_cat = 'turtle' ORDER BY mdl_seq")$mdl_seq
stopifnot("expected v7's six merged turtle models at 54238-54243" =
identical(as.integer(seqs), 54238:54243))
seq_sql <- paste(seqs, collapse = ", ")
dbExecute(con, "CREATE OR REPLACE TEMP VIEW turtle_src AS
SELECT t.mdl_seq AS ms_merge_key, tm.ds_key, mc.cell_id, mc.value AS val
FROM taxon t
JOIN taxon_model tm ON tm.taxon_id = t.taxon_id AND tm.ds_key <> 'ms_merge'
JOIN model_cell mc ON mc.mdl_seq = tm.mdl_seq
WHERE t.sp_cat = 'turtle'")
# half_even = TRUE reproduces R's rounding, which is what built v1-v7 (DuckDB round()
# differs on 28,598 v7 cells where ER = 50); ch_outer = TRUE keeps the 31 CH cells that
# lie outside the ER footprint, as v7's union did.
sql <- msens::turtle_sql(
DS_ER, DS_SUIT, DS_CH, src = "turtle_src",
suit_min = suit_min, fill = fill, half_even = TRUE, ch_outer = TRUE)
# materialize BEFORE the delete: turtle_src reads model_cell, which is about to change
dbExecute(con, glue("CREATE OR REPLACE TEMP TABLE turtle_new AS {sql}"))
dbExecute(con, glue("DELETE FROM model_cell WHERE mdl_seq IN ({seq_sql})"))
n <- dbExecute(con, glue(
"INSERT INTO model_cell (mdl_seq, cell_id, value)
SELECT mdl_key::INTEGER, cell_id, val::INTEGER FROM turtle_new"))
rule <- glue(
"v7.1 turtle core habitat: suit_min={suit_min}, ",
"fill={if (is.na(fill)) 'none' else fill}, half_even=TRUE, ch_outer=TRUE")
# idempotent: strip any previous v7.1 clause before appending, so a re-run does not stack
dbExecute(con, glue(
"UPDATE model SET description =
regexp_replace(coalesce(description, ''), ' \\| v7\\.1 turtle.*$', '') || ' | {rule}'
WHERE mdl_seq IN ({seq_sql})"))
log_info("turtles rewritten: {format(n, big.mark = ',')} cells under mdl_seq {seq_sql}")
n
}
# kappa: a ZONE reporting rule. Cell-level rasters are untouched; only the Program-Area
# component row of a pair below the floor is deleted (its raw and _prepctareaweighting rows
# stay, so the diagnostic remains readable), and the composite then averages one fewer
# component -- exactly as it already does for an area with no turtles at all.
apply_kappa <- function(con, tbl_pra, component_keys, kappa, f_csv) {
if (is.na(kappa) || kappa <= 0) {
log_info("COVERAGE_FLOOR off -- no _coverage metrics written, no component dropped")
return(invisible(NULL))
}
d <- dbGetQuery(con, msens::coverage_sql(
tbl_pra, component_keys, weight = "pct_covered",
val_col = msens::sdm_val_col(con, "cell_metric"))) |>
msens::coverage_floor(kappa = kappa)
z_seqs <- dbGetQuery(con, glue(
"SELECT zone_seq FROM zone WHERE tbl = '{tbl_pra}'"))$zone_seq
z_sql <- paste(z_seqs, collapse = ", ")
for (k in component_keys) {
k_cov <- glue("{k}_coverage")
if (!nrow(dbGetQuery(con, glue("SELECT 1 FROM metric WHERE metric_key = '{k_cov}'"))))
dbExecute(con, glue(
"INSERT INTO metric (metric_key, description) VALUES
('{k_cov}', 'Percent of the zone covered by cells carrying {k} (coverage floor input)')"))
s_cov <- dbGetQuery(con, glue(
"SELECT metric_seq FROM metric WHERE metric_key = '{k_cov}'"))$metric_seq
dbExecute(con, glue(
"DELETE FROM zone_metric WHERE metric_seq = {s_cov} AND zone_seq IN ({z_sql})"))
dd <- d |> filter(metric_key == k)
if (nrow(dd))
dbExecute(con, glue(
"INSERT INTO zone_metric (zone_seq, metric_seq, value) VALUES {vals}",
vals = paste(glue("({dd$zone_seq}, {s_cov}, {dd$coverage * 100})"), collapse = ", ")))
}
d_drop <- d |> filter(!reportable)
for (i in seq_len(nrow(d_drop)))
dbExecute(con, glue(
"DELETE FROM zone_metric
WHERE zone_seq = {d_drop$zone_seq[i]}
AND metric_seq = (SELECT metric_seq FROM metric WHERE metric_key = '{d_drop$metric_key[i]}')"))
d_out <- d |>
left_join(dbGetQuery(con, glue(
"SELECT zone_seq, value AS programarea_key FROM zone WHERE tbl = '{tbl_pra}'")),
by = "zone_seq") |>
arrange(coverage) |>
select(programarea_key, metric_key, coverage, reportable)
write_csv(d_out |> filter(!reportable), f_csv)
log_info("kappa={kappa}: {nrow(d_drop)} of {nrow(d)} (zone x component) pairs not reportable")
invisible(d_out)
}
# the methods a release actually used, as data -- backfill_versions.qmd stages this table and
# manifest_build() carries it as `methods`, so the docs state the rules from the release itself
write_release_method <- function(con, suit_min, fill, kappa) {
d <- tibble(
method_key = c("turtle_suit_min", "turtle_fill", "coverage_floor_pct", "prev"),
value = c(
as.character(suit_min),
if (is.na(fill)) "none" else as.character(fill),
as.character(if (is.na(kappa)) 0 else kappa * 100),
ver_prev),
description = c(
"AquaMaps suitability strictly below this is absent before the turtle multiply (v7: 0)",
"Suitability assumed for a turtle range cell with none retained; 'none' drops the cell unless critical habitat covers it (v7: 1)",
"A Program-Area component is reportable only if its scored cells cover at least this percent of the area (v7: 0)",
"The release this one patches"))
dbWriteTable(con, "release_method", d, overwrite = TRUE)
d
}
```
### One function, four runs
```{r}
#| label: build_release
build_release <- function(db_path, out_dir, suit_min, fill, kappa, label) {
t0 <- Sys.time()
dir_create(out_dir)
clone_db(sdm_db_prev, db_path)
con <- dbConnect(duckdb(dbdir = db_path, read_only = FALSE))
# v7's zone / metric / model / dataset tables carry `DEFAULT CURRENT_DATE`, which DuckDB binds
# to the icu extension -- so EVERY write to them needs it, the zone-table rename included.
# A laptop autoloads icu from ~/.duckdb; the server's rstudio container has no copy for
# uid 1000 and halted here on its first UPDATE. Declared, as the other notebooks declare httpfs.
dbExecute(con, "INSTALL icu; LOAD icu;")
on.exit(dbDisconnect(con, shutdown = TRUE), add = TRUE)
duckdb_tune(con, glue("{out_dir}/duckdb_tmp"))
# the server's clone inherits a cell_model VIEW whose read_parquet path is absolute and
# points at {ver_prev}: left in place it would serve the PREVIOUS release's cells under
# this release's name. Dropped, then asserted gone.
dbExecute(con, "DROP VIEW IF EXISTS cell_model")
v_bad <- dbGetQuery(con, glue(
"SELECT view_name FROM duckdb_views() WHERE sql LIKE '%/{ver_prev}/%'"))
stopifnot("a view still references the previous release's path" = nrow(v_bad) == 0)
# zone tbl suffix _v7 -> _v7b (dev/build_v7.R:55-61)
for (old_tbl in dbGetQuery(con, glue(
"SELECT DISTINCT tbl FROM zone WHERE tbl LIKE '%\\_{ver_prev}' ESCAPE '\\'"))$tbl) {
new_tbl <- str_replace(old_tbl, glue("_{ver_prev}$"), glue("_{ver}"))
dbExecute(con, glue("UPDATE zone SET tbl = '{new_tbl}' WHERE tbl = '{old_tbl}'"))
dbExecute(con, glue("UPDATE zone_taxon SET zone_tbl = '{new_tbl}' WHERE zone_tbl = '{old_tbl}'"))
log_info("renamed zone tbl: {old_tbl} -> {new_tbl}")
}
n_turtle <- rewrite_turtles(con, suit_min, fill)
# v7's scoring chunks need these names; a child of globalenv() keeps the chunk code's own
# assignments out of the session while still resolving library functions
env <- new.env(parent = globalenv())
env$con_sdm <- con
env$tbl_pra <- glue("ply_programareas_2026_{ver}")
env$tbl_sr <- glue("ply_subregions_2026_{ver}")
env$tbl_er <- "ply_ecoregions_2025"
env$component_pct_pra_csv <- file.path(out_dir, glue("tbl_component_pct_programarea_{ver}.csv"))
env$zone_taxon_csv <- file.path(out_dir, glue("zone_taxon_{ver}.csv"))
# NEVER narrow sp_cats: the pct-area and composite chunks enumerate every component from it
env$sp_cats <- dbGetQuery(con,
"SELECT DISTINCT sp_cat FROM taxon WHERE is_ok ORDER BY sp_cat")$sp_cat
stopifnot("sp_cats lost a category" = length(env$sp_cats) == 7)
run <- function(labels)
for (lbl in labels) {
log_info("[{label}] chunk: {lbl}")
eval(parse(text = chunk_fix(get_chunk(lbl))), envir = env)
}
run(CHUNKS_PRE)
d_cov <- apply_kappa(
con, env$tbl_pra,
c(glue("extrisk_{env$sp_cats}_ecoregion_rescaled"), "primprod_ecoregion_rescaled"),
kappa, file.path(out_dir, glue("tbl_component_not_reportable_programarea_{ver}.csv")))
run(CHUNKS_POST)
write_release_method(con, suit_min, fill, kappa)
log_info("[{label}] built in {round(as.numeric(difftime(Sys.time(), t0, units = 'mins')), 1)} min")
invisible(list(n_turtle = n_turtle, coverage = d_cov))
}
```
## Build `v7b`
```{r}
#| label: build
if (!do_build) {
log_info("{sdm_db} exists -- rendering report + gates read-only (REDO_BUILD=1 to rebuild)")
} else {
build_release(sdm_db, dir_v, TURTLE_SUIT_MIN, TURTLE_FILL, COVERAGE_FLOOR, ver)
}
stopifnot("release database missing after build" = file_exists(sdm_db))
con <- dbConnect(duckdb(dbdir = sdm_db, read_only = TRUE))
duckdb_tune(con, glue("{dir_big_v}/duckdb_tmp"))
dbExecute(con, glue("ATTACH '{sdm_db_prev}' AS prev (READ_ONLY)"))
tbl_pra <- glue("ply_programareas_2026_{ver}")
tbl_sr <- glue("ply_subregions_2026_{ver}")
tbl_er <- "ply_ecoregions_2025"
sp_cats <- dbGetQuery(con, "SELECT DISTINCT sp_cat FROM taxon WHERE is_ok ORDER BY sp_cat")$sp_cat
MDL_TURTLE <- dbGetQuery(con, "SELECT mdl_seq FROM taxon WHERE sp_cat = 'turtle' ORDER BY mdl_seq")$mdl_seq
mdl_sql <- paste(MDL_TURTLE, collapse = ", ")
K_SCORE <- "score_extriskspcat_primprod_ecoregionrescaled_equalweights"
```
## Artifacts {#sec-artifacts}
`layers_{ver}.csv`, the pct-area diagnostic and `zone_taxon_{ver}.csv` are written by the build
itself; the Program-Area geopackage is **regenerated** (v7's own gpkg is md5-identical to v6's — it
was copied, so its metric columns are stale) from the previous release's geometry columns plus this
release's `zone_metric` pivoted wide, per `old/calc_scores_v7.qmd:2893-2940`.
```{r}
#| label: artifacts
f_lyrs <- file.path(dir_v, glue("layers_{ver}.csv"))
f_gpkg <- file.path(dir_v, glue("ply_programareas_2026_{ver}.gpkg"))
f_srpa <- file.path(dir_v, "subregion_programareas.csv")
if (do_build || !file_exists(f_lyrs))
bind_rows(
tibble(order = 1, category = "Overall",
layer = "score", lyr = K_SCORE),
tibble(order = 2, category = "Species, rescaled by Ecoregion",
layer = glue("{sp_cats}: ext. risk, ecorgn"),
lyr = glue("extrisk_{sp_cats}_ecoregion_rescaled")),
tibble(order = 3, category = "Primary Productivity, rescaled by Ecoregion",
layer = "prim prod, ecorgn", lyr = "primprod_ecoregion_rescaled"),
tibble(order = 4, category = "Species, raw Extinction Risk",
layer = glue("{sp_cats}: ext. risk"), lyr = glue("extrisk_{sp_cats}")),
tibble(order = 5, category = "Primary Productivity, raw Phytoplankton",
lyr = "primprod", layer = "prim prod, 2014-2023 avg (mg C/m^2/day)")) |>
write_csv(f_lyrs)
# version-independent lookup: carried forward, not recomputed
if (!file_exists(f_srpa))
file_copy(file.path(dir_v_prev, "subregion_programareas.csv"), f_srpa)
if (do_build || !file_exists(f_gpkg)) {
d_pra_metrics <- dbGetQuery(con, glue(
"SELECT z.value AS programarea_key, m.metric_key, round(zm.value, 1) AS value
FROM zone z JOIN zone_metric zm USING (zone_seq) JOIN metric m USING (metric_seq)
WHERE z.tbl = '{tbl_pra}' AND zm.value IS NOT NULL")) |>
pivot_wider(names_from = metric_key, values_from = value)
base_cols <- c("programarea_id", "region_key", "region_name", "planarea_key", "planarea_name",
"programarea_key", "programarea_name", "ctr_lon", "ctr_lat", "area_km2")
read_sf(file.path(dir_v_prev, glue("ply_programareas_2026_{ver_prev}.gpkg"))) |>
select(any_of(base_cols)) |>
left_join(d_pra_metrics, by = "programarea_key") |>
st_write(f_gpkg, delete_dsn = TRUE, quiet = TRUE)
}
msens::report_table(
tibble(file = basename(c(f_lyrs, f_gpkg, f_srpa,
file.path(dir_v, glue("tbl_component_pct_programarea_{ver}.csv")),
file.path(dir_v, glue("zone_taxon_{ver}.csv")))),
exists = file_exists(c(f_lyrs, f_gpkg, f_srpa,
file.path(dir_v, glue("tbl_component_pct_programarea_{ver}.csv")),
file.path(dir_v, glue("zone_taxon_{ver}.csv")))),
kb = round(as.numeric(file_size(c(f_lyrs, f_gpkg, f_srpa,
file.path(dir_v, glue("tbl_component_pct_programarea_{ver}.csv")),
file.path(dir_v, glue("zone_taxon_{ver}.csv"))))) / 1024)),
caption = glue("Artifacts in {dir_v}"))
```
## Comparing a build against v7 {#sec-compare}
One function, used for `v7b` and for each control. It compares the three surfaces that carry the
release — `cell_metric` (~22.0 M rows), `zone_metric` and `zone_taxon` — **by key, exactly**, and
reports the largest value movement per metric.
A note on "exactly". `model_cell.value` is INTEGER, so the six turtle models are compared
bit-for-bit and the only admissible answer is 0 diffs. The scored surfaces are DOUBLE and every
zone value is a `SUM()` over up to 662,075 cells: IEEE addition is not associative, and DuckDB's
parallel aggregate does not fix the reduction order, so re-running v7's own SQL on v7's own data
reproduces `zone_metric` to about 5 ULP (measured: max 8.9e-15 on values near 10) rather than
bit-identically. The gate is therefore **exact on keys** (no row gained, none lost, none renamed)
and **bounded on values at 1e-9 relative** — four orders of magnitude tighter than the smallest
change either method makes, and six orders looser than reassociation noise. `n_moved` below counts
values outside that bound; it is the number the gates assert.
```{r}
#| label: compare_fn
# rel difference guard: |a-b| / max(|b|, 1) -- absolute for small values, relative for large
REL_TOL <- 1e-9
compare_build <- function(db_path) {
cc <- dbConnect(duckdb(dbdir = db_path, read_only = TRUE))
on.exit(dbDisconnect(cc, shutdown = TRUE), add = TRUE)
duckdb_tune(cc, glue("{path_dir(db_path)}/duckdb_tmp"))
# ATTACH is per DATABASE INSTANCE, not per connection: comparing the release this notebook
# already holds open would otherwise fail with `database with name "prev" already exists`.
if (!nrow(dbGetQuery(cc, "SELECT 1 FROM duckdb_databases() WHERE database_name = 'prev'")))
dbExecute(cc, glue("ATTACH '{sdm_db_prev}' AS prev (READ_ONLY)"))
ms <- paste(MDL_TURTLE, collapse = ", ")
turtle <- dbGetQuery(cc, glue(
"SELECT
(SELECT count(*) FROM (SELECT mdl_seq, cell_id, value FROM model_cell WHERE mdl_seq IN ({ms})
EXCEPT ALL SELECT mdl_seq, cell_id, value FROM prev.model_cell WHERE mdl_seq IN ({ms}))) AS new_not_prev,
(SELECT count(*) FROM (SELECT mdl_seq, cell_id, value FROM prev.model_cell WHERE mdl_seq IN ({ms})
EXCEPT ALL SELECT mdl_seq, cell_id, value FROM model_cell WHERE mdl_seq IN ({ms}))) AS prev_not_new,
(SELECT count(*) FROM model_cell WHERE mdl_seq IN ({ms})) AS n_new,
(SELECT count(*) FROM prev.model_cell WHERE mdl_seq IN ({ms})) AS n_prev,
(SELECT count(*) FROM (SELECT mdl_seq, cell_id FROM model_cell WHERE mdl_seq IN ({ms})
EXCEPT SELECT mdl_seq, cell_id FROM prev.model_cell WHERE mdl_seq IN ({ms}))) AS n_gained"))
cell <- dbGetQuery(cc, glue(
"WITH a AS (SELECT m.metric_key AS k, cm.cell_id AS cid, cm.value AS v
FROM cell_metric cm JOIN metric m USING (metric_seq)),
b AS (SELECT m.metric_key AS k, cm.cell_id AS cid, cm.value AS v
FROM prev.cell_metric cm JOIN prev.metric m USING (metric_seq)),
j AS (SELECT coalesce(a.k, b.k) AS metric_key, a.v AS va, b.v AS vb
FROM a FULL JOIN b ON a.k = b.k AND a.cid = b.cid)
SELECT metric_key,
count(*) FILTER (WHERE vb IS NULL) AS n_only_new,
count(*) FILTER (WHERE va IS NULL) AS n_only_prev,
coalesce(max(abs(va - vb) / greatest(abs(vb), 1.0)), 0) AS max_rel,
count(*) FILTER (WHERE abs(va - vb) / greatest(abs(vb), 1.0) > {REL_TOL}) AS n_moved
FROM j GROUP BY 1 ORDER BY 1"))
zone <- dbGetQuery(cc, glue(
"WITH a AS (SELECT z.fld AS f, z.value AS zv, m.metric_key AS k, zm.value AS v
FROM zone_metric zm JOIN zone z USING (zone_seq) JOIN metric m USING (metric_seq)),
b AS (SELECT z.fld AS f, z.value AS zv, m.metric_key AS k, zm.value AS v
FROM prev.zone_metric zm JOIN prev.zone z USING (zone_seq)
JOIN prev.metric m USING (metric_seq)),
j AS (SELECT coalesce(a.k, b.k) AS metric_key, a.v AS va, b.v AS vb
FROM a FULL JOIN b ON a.f = b.f AND a.zv = b.zv AND a.k = b.k)
SELECT metric_key,
count(*) FILTER (WHERE vb IS NULL) AS n_only_new,
count(*) FILTER (WHERE va IS NULL) AS n_only_prev,
coalesce(max(abs(va - vb) / greatest(abs(vb), 1.0)), 0) AS max_rel,
count(*) FILTER (WHERE abs(va - vb) / greatest(abs(vb), 1.0) > {REL_TOL}) AS n_moved
FROM j GROUP BY 1 ORDER BY 1"))
# zone_taxon: key = (zone_fld, zone_value, mdl_seq); zone_tbl differs by the _v7b rename
ztx <- dbGetQuery(cc, glue(
"WITH a AS (SELECT zone_fld f, zone_value zv, mdl_seq, sp_cat, area_km2, avg_suit,
suit_rl, suit_rl_area, cat_suit_rl_area, pct_cat FROM zone_taxon),
b AS (SELECT zone_fld f, zone_value zv, mdl_seq, sp_cat, area_km2, avg_suit,
suit_rl, suit_rl_area, cat_suit_rl_area, pct_cat FROM prev.zone_taxon),
j AS (SELECT coalesce(a.sp_cat, b.sp_cat) AS sp_cat,
a.area_km2 aa, b.area_km2 ba, a.avg_suit sa, b.avg_suit sb,
a.suit_rl ra, b.suit_rl rb, a.pct_cat pa, b.pct_cat pb,
a.f AS af, b.f AS bf
FROM a FULL JOIN b ON a.f = b.f AND a.zv = b.zv AND a.mdl_seq = b.mdl_seq)
-- named `metric_key` so moved_keys() reads the real column: it looked for metric_key,
-- found nothing on a frame keyed by sp_cat, and the zone_taxon gate passed vacuously
SELECT sp_cat AS metric_key,
count(*) FILTER (WHERE bf IS NULL) AS n_only_new,
count(*) FILTER (WHERE af IS NULL) AS n_only_prev,
coalesce(greatest(
max(abs(aa - ba) / greatest(abs(ba), 1.0)),
max(abs(sa - sb) / greatest(abs(sb), 1.0)),
max(abs(ra - rb) / greatest(abs(rb), 1.0)),
max(abs(pa - pb) / greatest(abs(pb), 1.0))), 0) AS max_rel,
count(*) FILTER (WHERE
greatest(abs(aa - ba) / greatest(abs(ba), 1.0),
abs(sa - sb) / greatest(abs(sb), 1.0),
abs(ra - rb) / greatest(abs(rb), 1.0),
abs(pa - pb) / greatest(abs(pb), 1.0)) > {REL_TOL}) AS n_moved
FROM j GROUP BY 1 ORDER BY 1"))
pra <- dbGetQuery(cc, glue(
"SELECT z.value AS programarea_key, m.metric_key, zm.value
FROM zone_metric zm JOIN zone z USING (zone_seq) JOIN metric m USING (metric_seq)
WHERE z.fld = 'programarea_key'"))
list(turtle = turtle, cell = cell, zone = zone, ztx = ztx, pra = pra)
}
# the invariant that must hold in EVERY run, whatever the flags: decision 4 keeps the
# productivity fix with v10. `{key}_coverage` is a NEW zone-only key, not a primprod row --
# it says how much of the zone the component covers and never touches primprod's values.
moved_keys <- function(d) {
# a check that reads a column that is not there returns character(0) and passes forever
stopifnot("comparison frame has no metric_key column" = "metric_key" %in% names(d))
sort(d$metric_key[d$n_moved > 0 | d$n_only_new > 0 | d$n_only_prev > 0])
}
primprod_ok <- function(cmp) {
k <- c(moved_keys(cmp$cell), moved_keys(cmp$zone))
!any(grepl("^primprod", k) & !grepl("_coverage$", k))
}
```
## Control runs {#sec-controls}
Three clones in a scratch directory, each through `build_release()` with a different flag set, then
deleted. The point of the all-off run is that it is the *same code path* as `v7b`: if it reproduces
v7, every difference `v7b` shows is attributable to a flag and nothing else.
```{r}
#| label: controls
CONTROLS <- tribble(
~run, ~suit_min, ~fill, ~kappa,
"all_off", 0, 1, 0,
"turtles_only", TURTLE_SUIT_MIN, TURTLE_FILL, 0,
"kappa_only", 0, 1, COVERAGE_FLOOR)
do_controls <- do_build || !file_exists(f_gates)
if (!do_controls) {
log_info("control gate results cached at {f_gates} -- not rebuilding the clones")
# an EMPTY moved-key set is the passing answer, and read_csv() reads "" back as NA:
# nchar(NA) == 0 is NA, so every gate below would halt on a cached re-render of a good build
d_gates <- read_csv(f_gates) |>
mutate(across(ends_with("_moved"), \(x) coalesce(as.character(x), "")))
d_attr <- read_csv(f_attr)
} else {
cmps <- list()
pras <- list()
for (i in seq_len(nrow(CONTROLS))) {
r <- CONTROLS[i, ]
od <- file.path(dir_scratch, r$run)
build_release(file.path(od, "sdm.duckdb"), od, r$suit_min, r$fill, r$kappa, r$run)
cmps[[r$run]] <- compare_build(file.path(od, "sdm.duckdb"))
pras[[r$run]] <- cmps[[r$run]]$pra
}
cmps[[ver]] <- compare_build(sdm_db)
pras[[ver]] <- cmps[[ver]]$pra
# one tidy row per (run, check) so the gates and a cached re-render assert the same numbers
d_gates <- imap_dfr(cmps, function(cmp, run) tibble(
run = run,
turtle_n = cmp$turtle$n_new,
turtle_gained= cmp$turtle$n_gained,
turtle_diff = cmp$turtle$new_not_prev + cmp$turtle$prev_not_new,
cell_moved = paste(moved_keys(cmp$cell), collapse = " "),
zone_moved = paste(moved_keys(cmp$zone), collapse = " "),
ztx_moved = paste(moved_keys(cmp$ztx), collapse = " "),
cell_max_rel = max(cmp$cell$max_rel), zone_max_rel = max(cmp$zone$max_rel),
ztx_max_rel = max(cmp$ztx$max_rel),
primprod_ok = primprod_ok(cmp)))
write_csv(d_gates, f_gates)
d_attr <- imap_dfr(pras, \(d, run) d |> mutate(run = run)) |>
filter(metric_key %in% c("extrisk_turtle_ecoregion_rescaled", K_SCORE)) |>
mutate(what = if_else(metric_key == K_SCORE, "composite", "turtle")) |>
select(programarea_key, what, run, value) |>
bind_rows(
dbGetQuery(con, glue(
"SELECT z.value AS programarea_key, m.metric_key, zm.value
FROM prev.zone_metric zm JOIN prev.zone z USING (zone_seq)
JOIN prev.metric m USING (metric_seq)
WHERE z.fld = 'programarea_key'")) |>
filter(metric_key %in% c("extrisk_turtle_ecoregion_rescaled", K_SCORE)) |>
mutate(what = if_else(metric_key == K_SCORE, "composite", "turtle"), run = ver_prev) |>
select(programarea_key, what, run, value))
write_csv(d_attr, f_attr)
dir_delete(dir_scratch)
log_info("control clones removed from {dir_scratch}")
}
msens::report_table(
d_gates |> select(run, turtle_n, turtle_gained, turtle_diff, primprod_ok,
cell_max_rel, zone_max_rel, ztx_max_rel),
caption = "Control runs vs v7 (turtle_diff is an exact INTEGER multiset difference)")
```
## Gates that stop the line
```{r}
#| label: gate_all_off
g <- d_gates |> filter(run == "all_off")
stopifnot(
"all flags off: the six turtle models must be v7's 2,022,909 cells, bit-identical" =
g$turtle_n == 2022909 && g$turtle_diff == 0,
"all flags off: cell_metric must not gain, lose or move a (metric_key, cell_id) value" =
nchar(g$cell_moved) == 0,
"all flags off: zone_metric must not gain, lose or move a (zone, metric_key) value" =
nchar(g$zone_moved) == 0,
"all flags off: zone_taxon must not gain, lose or move a row" =
nchar(g$ztx_moved) == 0)
log_info("GATE all-off == v7: turtles {g$turtle_n} cells / {g$turtle_diff} diffs; ",
"max rel dev cell {signif(g$cell_max_rel, 3)}, zone {signif(g$zone_max_rel, 3)}, ",
"zone_taxon {signif(g$ztx_max_rel, 3)} (tol {REL_TOL})")
```
```{r}
#| label: gate_single_flags
# the turtles flag may move only the turtle component, the composite that averages it, and
# turtle zone_taxon rows. kappa may move only the ONE pair it drops on v7 (St. George Basin
# turtles, 1.5 % coverage), the composite that no longer averages it, and the new _coverage
# keys -- and no cell metric at all, because it is a zone reporting rule. Naming the dropped
# key rather than allowing any component is the point: if a second pair falls below the
# floor, that is news, not noise.
keys_of <- function(x) setdiff(strsplit(x, " ")[[1]], "")
ok_turtle <- function(keys) all(grepl("^extrisk_turtle|^score_", keys))
ok_kappa <- function(keys) all(
keys %in% c("extrisk_turtle_ecoregion_rescaled", K_SCORE) | grepl("_coverage$", keys))
gt <- d_gates |> filter(run == "turtles_only")
gk <- d_gates |> filter(run == "kappa_only")
stopifnot(
"turtles-only: only extrisk_turtle* and the score may move in cell_metric" =
ok_turtle(keys_of(gt$cell_moved)),
"turtles-only: only extrisk_turtle* and the score may move in zone_metric" =
ok_turtle(keys_of(gt$zone_moved)),
"turtles-only: only turtle rows may move in zone_taxon" =
all(keys_of(gt$ztx_moved) %in% "turtle"),
"turtles-only: the turtle surface must actually have changed" =
gt$turtle_diff > 0 && gt$turtle_n < 2022909,
"kappa-only: no cell metric may move -- kappa is a zone reporting rule" =
nchar(gk$cell_moved) == 0,
"kappa-only: zone_taxon must not move" =
nchar(gk$ztx_moved) == 0,
"kappa-only: only the dropped component, the score and the new _coverage keys may move" =
ok_kappa(keys_of(gk$zone_moved)),
"kappa-only: the turtle models must stay v7's" =
gk$turtle_diff == 0 && gk$turtle_n == 2022909)
log_info("GATE single flags: turtles-only moved cell[{gt$cell_moved}] zone[{gt$zone_moved}] ",
"ztx[{gt$ztx_moved}]; kappa-only moved zone[{gk$zone_moved}]")
```
```{r}
#| label: gate_primprod
# decision 4: the productivity fix is NOT in v7.1. Every primprod row, cell and zone, must
# equal v7's in EVERY run -- the control runs included.
stopifnot("a primprod row moved in some run -- the productivity fix is v10's, not v7.1's" =
all(d_gates$primprod_ok))
log_info("GATE primprod unchanged in all {nrow(d_gates)} runs")
```
```{r}
#| label: gate_turtle_surface
n_ch <- dbGetQuery(con, glue(
"WITH src AS (
SELECT t.mdl_seq AS mdl_seq, tm.ds_key, mc.cell_id, mc.value AS val
FROM prev.taxon t
JOIN prev.taxon_model tm ON tm.taxon_id = t.taxon_id
JOIN prev.model_cell mc ON mc.mdl_seq = tm.mdl_seq
WHERE t.sp_cat = 'turtle' AND tm.ds_key IN ({paste(sprintf(\"'%s'\", DS_CH), collapse = ', ')})),
ch AS (SELECT mdl_seq, cell_id FROM src GROUP BY 1, 2)
SELECT count(*) AS n_ch,
count(*) FILTER (WHERE n.value IS NULL) AS n_lost,
count(*) FILTER (WHERE n.value IS DISTINCT FROM o.value) AS n_changed
FROM ch
LEFT JOIN model_cell n ON n.mdl_seq = ch.mdl_seq AND n.cell_id = ch.cell_id
LEFT JOIN prev.model_cell o ON o.mdl_seq = ch.mdl_seq AND o.cell_id = ch.cell_id"))
d_typ <- dbGetQuery(con,
"SELECT (SELECT column_type FROM (DESCRIBE SELECT * FROM model_cell) WHERE column_name = 'value') AS mc_type,
(SELECT column_type FROM (DESCRIBE SELECT * FROM cell_metric) WHERE column_name = 'value') AS cm_type")
d_eco <- dbGetQuery(con, glue(
"SELECT z.value AS ecoregion, mc.mdl_seq, count(DISTINCT mc.cell_id) AS n
FROM model_cell mc
JOIN zone_cell zc ON zc.cell_id = mc.cell_id
JOIN zone z ON z.zone_seq = zc.zone_seq AND z.fld = 'ecoregion_key'
WHERE mc.mdl_seq IN ({mdl_sql})
GROUP BY 1, 2 ORDER BY 1, 2"))
d_goa <- d_eco |> filter(ecoregion == "GOA")
gv <- d_gates |> filter(run == ver)
stopifnot(
"critical habitat was trimmed -- every v7 CH cell must survive with its v7 value" =
n_ch$n_lost == 0 && n_ch$n_changed == 0 && n_ch$n_ch > 0,
"mdl_seq 54238-54243 must be unchanged" =
identical(as.integer(MDL_TURTLE), 54238:54243),
"model_cell.value must stay INTEGER and cell_metric.value DOUBLE (DuckDB hash() is type-sensitive)" =
d_typ$mc_type == "INTEGER" && d_typ$cm_type == "DOUBLE",
"no turtle cell may be GAINED" =
gv$turtle_gained == 0,
"no turtle cell may remain in the East Bering Sea, Chukchi/Beaufort Sea or High Arctic ecoregions" =
!any(d_eco$ecoregion %in% c("EBS", "CBS", "HAR")),
"Gulf of Alaska ecoregion must hold the leatherback alone, 9,608 cells" =
nrow(d_goa) == 1 && d_goa$mdl_seq == 54241 && d_goa$n == 9608)
log_info("GATE turtle surface: CH {n_ch$n_ch} cells, 0 lost / 0 changed; types {d_typ$mc_type}/",
"{d_typ$cm_type}; v7b turtle cells {gv$turtle_n} (gained {gv$turtle_gained}); GOA ",
"ecoregion mdl {d_goa$mdl_seq} x {d_goa$n}")
msens::report_table(d_eco, caption = "v7b turtle cells by ecoregion (EBS / CBS / HAR must be absent)")
```
## Report
### Attribution: what each flag moves, per Program Area
```{r}
#| label: tbl_attribution
lvl <- c(ver_prev, "all_off", "turtles_only", "kappa_only", ver)
d_w <- d_attr |>
filter(run %in% lvl) |>
mutate(run = factor(run, levels = lvl)) |>
pivot_wider(names_from = run, values_from = value) |>
arrange(what, programarea_key)
rank_of <- function(x) rank(-x, na.last = "keep", ties.method = "min")
d_comp <- d_w |> filter(what == "composite") |>
mutate(across(where(is.numeric), \(x) round(x, 1)))
d_turt <- d_w |> filter(what == "turtle") |>
mutate(across(where(is.numeric), \(x) round(x, 1)))
msens::report_table(
d_comp |> select(-what) |> mutate(rank_v7 = rank_of(.data[[ver_prev]]), rank_v7b = rank_of(.data[[ver]])),
caption = "Composite score per Program Area: v7 | +turtles | +kappa | v7.1, with ranks")
msens::report_table(d_turt |> select(-what),
caption = "Turtle component per Program Area (blank = not reportable / no cells)")
sp <- suppressWarnings(cor(d_comp[[ver_prev]], d_comp[[ver]], method = "spearman", use = "complete.obs"))
log_info("Spearman rank correlation of the composite, {ver_prev} vs {ver}: {round(sp, 4)}")
```
### Expected vs actual
The expectations are the emulation recorded in the plan (`§3.3`, `§5 P3`). They are **printed, not
asserted**: an emulation of v7's chain reproduces its published turtle scores to within 0.5, so a
Program Area further than about 1.0 from its expectation is a result to explain, not a tolerance to
widen.
```{r}
#| label: tbl_expected
d_exp <- tribble(
~programarea_key, ~turtle_v7_plan, ~turtle_exp, ~comp_v7_plan, ~comp_exp,
"ALA", 70.9, NA, 28.3, 22.3,
"GOA", 69.6, 11.6, 39.5, 32.2,
"KOD", 52.9, NA, 32.1, 29.1,
"COK", 48.8, NA, 52.6, 53.2,
"SHU", 46.4, NA, 29.7, 27.3,
"GEO", 0.9, NA, 33.9, 38.6,
"GAA", 54.4, 52.4, 33.1, 32.8,
"GAB", 26.4, 23.9, 25.7, 25.4,
"SOC", 26.8, 35.2, 29.7, 30.7,
"CEC", 18.6, 27.5, 21.9, 23.0,
"NOC", 14.4, 23.1, 20.7, 21.8)
d_chk <- d_exp |>
left_join(d_turt |> select(programarea_key, turtle_v7_act = all_of(ver_prev),
turtle_v7b = all_of(ver)), by = "programarea_key") |>
left_join(d_comp |> select(programarea_key, comp_v7_act = all_of(ver_prev),
comp_v7b = all_of(ver)), by = "programarea_key") |>
mutate(d_turtle = round(turtle_v7b - turtle_exp, 1),
d_comp = round(comp_v7b - comp_exp, 1),
flag = if_else(abs(coalesce(d_comp, 0)) > 1.0 |
abs(coalesce(d_turtle, 0)) > 1.0, "CHECK", "")) |>
select(programarea_key, turtle_v7 = turtle_v7_act, turtle_v7b, turtle_exp, d_turtle,
comp_v7 = comp_v7_act, comp_v7b, comp_exp, d_comp, flag)
if (any(abs(d_exp$comp_v7_plan - d_chk$comp_v7) > 0.2))
log_warn("the plan's v7 composite column disagrees with the release it was measured from: ",
"{paste(d_chk$programarea_key[abs(d_exp$comp_v7_plan - d_chk$comp_v7) > 0.2], collapse = ', ')}")
msens::report_table(d_chk, caption =
"Emulated expectation vs the built release (turtle blank = component absent)")
if (any(d_chk$flag == "CHECK"))
log_warn("Program Areas beyond ~1.0 of the emulated expectation: {paste(d_chk$programarea_key[d_chk$flag == 'CHECK'], collapse = ', ')}")
```
### Critical-habitat retention (§3.4)
```{r}
#| label: tbl_ch
msens::report_table(dbGetQuery(con, glue(
"WITH src AS (
SELECT t.mdl_seq, t.common_name, tm.ds_key, mc.cell_id, mc.value AS val
FROM prev.taxon t
JOIN prev.taxon_model tm ON tm.taxon_id = t.taxon_id
JOIN prev.model_cell mc ON mc.mdl_seq = tm.mdl_seq
WHERE t.sp_cat = 'turtle'),
ch AS (SELECT mdl_seq, common_name, cell_id FROM src
WHERE ds_key IN ({paste(sprintf(\"'%s'\", DS_CH), collapse = ', ')}) GROUP BY 1, 2, 3),
suit AS (SELECT mdl_seq, cell_id, max(val) AS p FROM src WHERE ds_key = '{DS_SUIT}' GROUP BY 1, 2)
SELECT ch.common_name AS species, count(*) AS ch_cells,
count(*) FILTER (WHERE suit.p IS NULL) AS no_aquamaps,
count(*) FILTER (WHERE suit.p < {TURTLE_SUIT_MIN}) AS below_threshold,
min(suit.p) AS min_p_in_ch,
count(*) FILTER (WHERE n.value IS NULL) AS lost_in_v7b
FROM ch LEFT JOIN suit USING (mdl_seq, cell_id)
LEFT JOIN model_cell n ON n.mdl_seq = ch.mdl_seq AND n.cell_id = ch.cell_id
GROUP BY 1 ORDER BY ch_cells DESC")),
caption = "Critical habitat: half of loggerhead and green CH has no AquaMaps value at all, so no threshold could cover it -- the override is what keeps it")
```
### Species presence per Program Area, before and after
```{r}
#| label: tbl_presence
msens::report_table(dbGetQuery(con, glue(
"WITH z AS (SELECT zone_seq, value AS pa FROM zone WHERE fld = 'programarea_key'),
n AS (SELECT z.pa, mc.mdl_seq, count(*) AS n
FROM model_cell mc JOIN zone_cell zc ON zc.cell_id = mc.cell_id
JOIN z ON z.zone_seq = zc.zone_seq WHERE mc.mdl_seq IN ({mdl_sql}) GROUP BY 1, 2),
o AS (SELECT z.pa, mc.mdl_seq, count(*) AS n
FROM prev.model_cell mc JOIN zone_cell zc ON zc.cell_id = mc.cell_id
JOIN z ON z.zone_seq = zc.zone_seq WHERE mc.mdl_seq IN ({mdl_sql}) GROUP BY 1, 2)
SELECT coalesce(n.pa, o.pa) AS programarea_key, t.common_name AS species,
coalesce(o.n, 0) AS cells_{ver_prev}, coalesce(n.n, 0) AS cells_{ver}
FROM o FULL JOIN n ON n.pa = o.pa AND n.mdl_seq = o.mdl_seq
JOIN taxon t ON t.mdl_seq = coalesce(n.mdl_seq, o.mdl_seq)
ORDER BY 1, 2")),
caption = glue("Turtle cells per Program Area, {ver_prev} -> {ver}"))
```
### The ten lowest coverages, and what κ removed
```{r}
#| label: tbl_coverage
component_keys <- c(glue("extrisk_{sp_cats}_ecoregion_rescaled"), "primprod_ecoregion_rescaled")
d_cov <- dbGetQuery(con, msens::coverage_sql(
tbl_pra, component_keys, weight = "pct_covered",
val_col = msens::sdm_val_col(con, "cell_metric"))) |>
msens::coverage_floor(kappa = if (is.na(COVERAGE_FLOOR)) 0 else COVERAGE_FLOOR) |>
left_join(dbGetQuery(con, glue(
"SELECT zone_seq, value AS programarea_key FROM zone WHERE tbl = '{tbl_pra}'")),
by = "zone_seq") |>
arrange(coverage) |>
select(programarea_key, metric_key, coverage, reportable)
msens::report_table(head(d_cov, 10) |> mutate(coverage = round(coverage * 100, 2)),
caption = glue("Ten lowest (Program Area x component) coverages in {ver}, percent"))
f_nr <- file.path(dir_v, glue("tbl_component_not_reportable_programarea_{ver}.csv"))
msens::report_table(
if (file_exists(f_nr)) read_csv(f_nr) else tibble(note = "kappa off -- nothing dropped"),
caption = "Pairs below the floor (removed from the composite, never published as a zero)")
```
### US-EEZ location quotient (§3.5) — diagnostic only
Tim's national reading: turtle weight per km², with **100 = the study-area average**. It is printed
here as a diagnostic, not published — the released metric is v10's distribution share.
```{r}
#| label: tbl_lq
lq_sql <- function(src) glue(
"WITH t AS (SELECT cell_id, sum(value) AS v FROM {src} WHERE mdl_seq IN ({mdl_sql}) GROUP BY 1),
nat AS (SELECT sum(t.v * c.area_km2) / sum(c.area_km2) AS m
FROM cell c LEFT JOIN t USING (cell_id))
SELECT z.value AS programarea_key,
sum(coalesce(t.v, 0) * c.area_km2 * zc.pct_covered / 100.0) /
sum(c.area_km2 * zc.pct_covered / 100.0) / (SELECT m FROM nat) * 100 AS lq,
sum(coalesce(t.v, 0) * c.area_km2 * zc.pct_covered / 100.0) /
(SELECT sum(coalesce(t.v, 0) * c.area_km2) FROM cell c LEFT JOIN t USING (cell_id)) * 100 AS pct_us
FROM zone z JOIN zone_cell zc USING (zone_seq) JOIN cell c ON c.cell_id = zc.cell_id
LEFT JOIN t ON t.cell_id = zc.cell_id
WHERE z.fld = 'programarea_key' GROUP BY 1 ORDER BY lq DESC")
d_lq <- dbGetQuery(con, lq_sql("prev.model_cell")) |>
rename(lq_v7 = lq, pct_us_v7 = pct_us) |>
full_join(dbGetQuery(con, lq_sql("model_cell")) |>
rename(lq_v7b = lq, pct_us_v7b = pct_us), by = "programarea_key") |>
mutate(across(where(is.numeric), \(x) round(x, 1))) |>
arrange(desc(lq_v7))
msens::report_table(d_lq, caption = "Turtle location quotient (100 = study-area average per km2) and share of US turtle weight")
ak <- c("ALA", "BOW", "CHU", "COK", "GEO", "GOA", "HOP", "KOD", "NAV", "NOR", "SHU", "MAT")
log_info("Alaska Program Areas hold {round(sum(d_lq$pct_us_v7[d_lq$programarea_key %in% ak], na.rm = TRUE), 2)}% ",
"of US turtle weight in {ver_prev}, {round(sum(d_lq$pct_us_v7b[d_lq$programarea_key %in% ak], na.rm = TRUE), 2)}% in {ver}")
```
### Sensitivity: other thresholds, and Tim's percentile rules {#sec-sensitivity}
v7's chain, emulated over the turtle component alone: rebuild the six merged surfaces, sum them per
cell (`extrisk_turtle`, since `er_score = 100` for a spatial-ER taxon), rescale on the ecoregional
min/max, take the `pct_covered`-weighted Program-Area mean, multiply by `pct_area`. Run against
v7's own surface it reproduces the published turtle column, which is what licenses the variants.
The percentile rows threshold the **merged** surface as Tim ran it (area-weighted percentile per
species), not the suitability before the merge — which is why they remove nothing for the five
species whose merged area is mostly floor-of-1 fill.
```{r}
#| label: emulate
# reads v7 through the `prev` attachment of the connection this notebook already holds: a
# second DuckDB instance on the same file would only invite the lock and alias trouble, and
# every input here (taxon_model, the raw model cells, zone_cell) is v7's, unchanged by v7.1.
emulate_turtle_pa <- function(cc, surface_sql) {
dbExecute(cc, glue("CREATE OR REPLACE TEMP TABLE e_surf AS {surface_sql}"))
dbGetQuery(cc, glue(
"WITH er AS (SELECT cell_id, sum(val) AS v FROM e_surf GROUP BY 1),
eco AS (SELECT zc.cell_id, zc.zone_seq, zc.pct_covered
FROM prev.zone_cell zc JOIN prev.zone z USING (zone_seq)
WHERE z.fld = 'ecoregion_key'),
tot AS (SELECT cell_id, sum(pct_covered) AS t FROM eco GROUP BY 1),
ecn AS (SELECT e.cell_id, e.zone_seq,
CASE WHEN count(*) OVER (PARTITION BY e.cell_id) = 1 THEN 100
ELSE e.pct_covered * 100.0 / t.t END AS npct
FROM eco e JOIN tot t USING (cell_id)),
mm AS (SELECT ecn.zone_seq, min(er.v) AS lo, max(er.v) AS hi
FROM ecn JOIN er USING (cell_id) GROUP BY 1),
rs AS (SELECT er.cell_id,
sum((er.v - mm.lo) / (mm.hi - mm.lo) * ecn.npct / 100.0) * 100 AS v
FROM er JOIN ecn USING (cell_id) JOIN mm ON mm.zone_seq = ecn.zone_seq
WHERE mm.hi > mm.lo GROUP BY 1),
pz AS (SELECT zone_seq, value AS pa FROM prev.zone WHERE fld = 'programarea_key'),
agg AS (SELECT pz.pa,
sum(rs.v * zc.pct_covered) / sum(zc.pct_covered) AS mean_v,
sum(zc.pct_covered) AS cov
FROM pz JOIN prev.zone_cell zc USING (zone_seq) JOIN rs USING (cell_id) GROUP BY 1),
den AS (SELECT pz.pa, sum(zc.pct_covered) AS tot
FROM pz JOIN prev.zone_cell zc USING (zone_seq) GROUP BY 1)
SELECT agg.pa AS programarea_key,
round(agg.mean_v * agg.cov / den.tot, 1) AS value,
round(agg.cov / den.tot * 100) AS pct_covered
FROM agg JOIN den USING (pa) ORDER BY 1"))
}
cc <- con
dbExecute(cc, "CREATE OR REPLACE TEMP VIEW turtle_src AS
SELECT t.mdl_seq AS ms_merge_key, tm.ds_key, mc.cell_id, mc.value AS val
FROM prev.taxon t
JOIN prev.taxon_model tm ON tm.taxon_id = t.taxon_id AND tm.ds_key <> 'ms_merge'
JOIN prev.model_cell mc ON mc.mdl_seq = tm.mdl_seq WHERE t.sp_cat = 'turtle'")
surf_v7 <- glue("SELECT mdl_seq AS mdl_key, cell_id, value AS val FROM prev.model_cell
WHERE mdl_seq IN ({mdl_sql})")
surf_thr <- function(sm, fl) msens::turtle_sql(
DS_ER, DS_SUIT, DS_CH, src = "turtle_src",
suit_min = sm, fill = fl, half_even = TRUE, ch_outer = TRUE)
# Tim's rule: an AREA-WEIGHTED percentile of the MERGED surface, per species, then keep >= it
surf_pct <- function(q) glue(
"WITH s AS (SELECT mc.mdl_seq, mc.cell_id, mc.value AS val, c.area_km2
FROM prev.model_cell mc JOIN prev.cell c USING (cell_id)
WHERE mc.mdl_seq IN ({mdl_sql})),
o AS (SELECT mdl_seq, val, area_km2,
sum(area_km2) OVER (PARTITION BY mdl_seq ORDER BY val
ROWS BETWEEN UNBOUNDED PRECEDING AND CURRENT ROW) AS cum,
sum(area_km2) OVER (PARTITION BY mdl_seq) AS tot FROM s),
q AS (SELECT mdl_seq, min(val) AS cut FROM o WHERE cum >= tot * {q} GROUP BY 1)
SELECT s.mdl_seq AS mdl_key, s.cell_id, s.val
FROM s JOIN q USING (mdl_seq) WHERE s.val >= q.cut")
variants <- list(
"v7 (as published)" = surf_v7,
"P >= 0.4, fill dropped" = surf_thr(40, NA),
"P >= 0.5, fill dropped" = surf_thr(50, NA),
"P >= 0.6, fill dropped" = surf_thr(60, NA),
"P >= 0.5, fill KEPT" = surf_thr(50, 1),
"Tim p15 on merged" = surf_pct(0.15),
"Tim p20 on merged" = surf_pct(0.20))
d_sens <- imap_dfr(variants, \(s, nm) emulate_turtle_pa(cc, s) |> mutate(rule = nm)) |>
mutate(cell = if_else(is.na(value), NA_character_, glue("{value} | {pct_covered}"))) |>
select(programarea_key, rule, cell) |>
pivot_wider(names_from = rule, values_from = cell)
msens::report_table(
d_sens |> filter(programarea_key %in% d_exp$programarea_key) |>
arrange(match(programarea_key, d_exp$programarea_key)),
caption = "Turtle component, score | % of area covered (emulated; blank = no turtle cells)")
msens::report_table(d_sens |> filter(!programarea_key %in% d_exp$programarea_key),
caption = "The remaining Program Areas")
```
## Manifest
Digests of the surfaces P4 must reproduce on the server: the six turtle models, `cell_metric`,
`zone_metric` and `zone_taxon`. `msens::hash_query()` is an order-independent reduction
(`count + bit_xor + sum`), so it does not care how DuckDB partitioned the scan.
It DOES care about the last bit of a DOUBLE, and @sec-compare measured that the zone sums are
not bit-reproducible (parallel `SUM()` reassociates; ~5 ULP). A digest of the raw doubles would
therefore differ between two correct builds, and the server-equals-laptop gate would be red for
no reason. So the INTEGER turtle models are digested as they are (exact), and every DOUBLE is
digested as `CAST(... AS FLOAT)`: ~7 significant digits at every magnitude (`area_km2` runs to
1e6, a score to 100), nine orders coarser than the noise and far finer than anything either
rule moves.
**The committed manifest is the reference.** A build with the same flags on any machine must
reproduce its digests, and this chunk stops if it does not -- that is P4's "server build equals
the laptop's" gate, run by the render itself rather than by eye. `BUILD_REFERENCE=1` is the
only way to move the reference, and says so in the log.
```{r}
#| label: manifest
h <- list(
turtle_models = msens::hash_query(con, glue(
"SELECT mdl_seq, cell_id, value FROM model_cell WHERE mdl_seq IN ({mdl_sql})")),
cell_metric = msens::hash_query(con,
"SELECT m.metric_key, cm.cell_id, CAST(cm.value AS FLOAT) AS value
FROM cell_metric cm JOIN metric m USING (metric_seq)"),
zone_metric = msens::hash_query(con,
"SELECT z.fld, z.value AS zone_value, m.metric_key, CAST(zm.value AS FLOAT) AS value
FROM zone_metric zm JOIN zone z USING (zone_seq) JOIN metric m USING (metric_seq)"),
zone_taxon = msens::hash_query(con,
"SELECT zone_fld, zone_value, mdl_seq, CAST(area_km2 AS FLOAT) AS area_km2,
CAST(avg_suit AS FLOAT) AS avg_suit, CAST(suit_rl AS FLOAT) AS suit_rl,
CAST(pct_cat AS FLOAT) AS pct_cat
FROM zone_taxon"))
# the reference gate: same flags => same digests as the committed manifest ----
f_mfst <- here(glue("data/manifests/build_{ver}.json"))
ref <- if (file_exists(f_mfst)) read_json(f_mfst) else NULL
fill_chr <- if (is.na(TURTLE_FILL)) "none" else as.character(TURTLE_FILL)
kappa <- if (is.na(COVERAGE_FLOOR)) 0 else COVERAGE_FLOOR
is_ref <- nzchar(Sys.getenv("BUILD_REFERENCE"))
same_flags <- !is.null(ref) &&
isTRUE(ref$turtle_suit_min == TURTLE_SUIT_MIN) && isTRUE(ref$turtle_fill == fill_chr) &&
isTRUE(ref$coverage_floor == kappa)
d_dig <- tibble(
surface = names(h),
digest = unlist(h),
reference = vapply(names(h), \(k) ref[[glue("digest_{k}")]] %||% NA_character_, ""))
if (is_ref || is.null(ref)) {
log_info("BUILD_REFERENCE / no manifest yet: these digests BECOME the reference")
} else if (!same_flags) {
log_info("flags differ from the committed manifest's -- a control render, reference not compared")
} else {
stopifnot("this build's digests differ from the committed reference (data/manifests); a correct build with the same flags reproduces them on any machine" =
all(d_dig$digest == d_dig$reference))
log_info("GATE reference digests reproduced: {paste(d_dig$surface, collapse = ', ')}")
}
gv <- d_gates |> filter(run == ver)
ga <- d_gates |> filter(run == "all_off")
# a control render (flags differ) never overwrites the reference it was not compared against
if (is_ref || is.null(ref) || same_flags) msens::write_manifest(
f_mfst,
target = glue("build_{ver}"),
content_hash = msens::hash_query(con, glue(
"SELECT '{h$turtle_models}' AS a, '{h$cell_metric}' AS b,
'{h$zone_metric}' AS c, '{h$zone_taxon}' AS d")),
stats = list(
ver = ver,
ver_prev = ver_prev,
turtle_suit_min = TURTLE_SUIT_MIN,
turtle_fill = fill_chr,
coverage_floor = kappa,
digest_turtle_models = h$turtle_models,
digest_cell_metric = h$cell_metric,
digest_zone_metric = h$zone_metric,
digest_zone_taxon = h$zone_taxon,
n_turtle_cells = gv$turtle_n,
n_turtle_cells_prev = ga$turtle_n,
n_turtle_gained = gv$turtle_gained,
allflagsoff_diffs = ga$turtle_diff,
n_ch_cells = n_ch$n_ch,
n_ch_lost = n_ch$n_lost,
goa_ecoregion_cells = d_goa$n,
n_not_reportable = sum(!d_cov$reportable),
spearman_composite = round(sp, 4)),
force = msens::force_target(glue("build_{ver}")))
msens::report_table(d_dig,
caption = "Content digests (DOUBLEs as FLOAT) vs the committed reference any same-flag build must reproduce")
dbDisconnect(con, shutdown = TRUE)
```