From e7b10905a8e730b177ee2f8c234310f15b2d9f91 Mon Sep 17 00:00:00 2001 From: aurascoper Date: Fri, 28 Aug 2026 13:16:00 -0500 Subject: [PATCH 01/13] Add the slab depth geometry, and two verifiers that can fail Geometry now enters radiodialysis_rhs() only through face weights built in the parameter constructor: w+- = (r +- dr/2)/r cylindrical, 1 for the slab. The axis/substratum node is untouched, and the cylindrical path is byte equivalent to the previous inline stencil (residual 2.28e-16). slab_parms() declares three semantics rather than inheriting them (rule 4): the Robin coefficient is k_L and not P_eff, radiation is off because no membrane exists at a biofilm/liquid interface, and D_eff is molecular rather than the reactor-scale dispersion value. X_total has no default and stops, naming RADIODIALYSIS: BLOCKED, because lambda ~ 1/sqrt(U) puts "is any gradient resolvable" downstream of that gate. Two independent verifiers, both against c(z) = c_b cosh(z/lam)/cosh(L/lam): verify_biofilm_depth_profile.R the shipping R, 13 checks verify_biofilm_depth_profile.py a numpy re-coding of the same scheme They agree to the digit on the shared quantity (max rel err 2.31e-07 at N=40), so agreement rules out a shared transcription error. Second order confirmed by refinement ratio 4.21 -> 4.05 over N=20..160 at both phi=0.306 and phi=3.7, the tolerance being on the normalised constant C rather than a raw error, since C is N-independent and a fixed bound cannot serve both regimes. Negative controls, per rule 1. Cylindrical weights run against the slab solution must fail the slab bound, and do by 3 to 5 orders. The R harness mutates face_weights (swapped, perturbed 0.1%, forced to 1) and requires each mutant caught, gated on a clean baseline first, because a wrong reference makes every mutant read as caught. That gate is itself a passing check, so it is structural rather than remembered. CI asserts on the receipt, not the exit status: parse() is not execution and was green over code that ran nowhere. A verifier that exits 0 having skipped its deSolve check is the all-pass-over-skip bug moved into YAML, so the job reads checks_run/failures/skipped out of the JSON and fails on a count that drifts from EXPECTED_CHECKS. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/model-contracts.yml | 31 ++ analysis/verify_biofilm_depth_profile.R | 275 +++++++++++++++++ analysis/verify_biofilm_depth_profile.py | 373 +++++++++++++++++++++++ biofilms_radiodialysis.R | 188 ++++++++++-- 4 files changed, 842 insertions(+), 25 deletions(-) create mode 100644 analysis/verify_biofilm_depth_profile.R create mode 100644 analysis/verify_biofilm_depth_profile.py diff --git a/.github/workflows/model-contracts.yml b/.github/workflows/model-contracts.yml index dbcd6aa..d35b9fd 100644 --- a/.github/workflows/model-contracts.yml +++ b/.github/workflows/model-contracts.yml @@ -97,6 +97,37 @@ jobs: cat("PARSE_OK", f, "\n") } + # parse() is not execution. It passes on a file whose every numerical + # result is wrong, so this job was green over code that ran nowhere. + - name: Install deSolve + run: sudo apt-get update && sudo apt-get install -y r-cran-desolve + + - name: Run the radiodialysis depth-geometry verifier + run: Rscript analysis/verify_biofilm_depth_profile.R --report r-verify.json + + # Assert on the RECEIPT, not on the exit status. If the setup step above + # fails or is dropped, the verifier reports its deSolve check as SKIPPED; + # a job that exits 0 having skipped is the ALL-PASS-over-skip bug moved + # from R into YAML. EXPECTED_CHECKS must be bumped deliberately when + # checks are added, so silently losing one fails the build. + - name: Assert the verifier actually ran everything + shell: Rscript {0} + run: | + EXPECTED_CHECKS <- 13L + if (!file.exists("r-verify.json")) + stop("no receipt: the verifier did not run to completion") + j <- readLines("r-verify.json", warn = FALSE) + num <- function(k) as.integer(sub( + paste0('.*"', k, '": *([0-9]+).*'), "\\1", j)) + ran <- num("checks_run"); fail <- num("failures"); skip <- num("skipped") + cat(sprintf("receipt: %d run, %d failed, %d skipped\n", ran, fail, skip)) + if (skip > 0L) stop("verifier skipped ", skip, " check(s): ", j) + if (fail > 0L) stop("verifier reported ", fail, " failure(s)") + if (ran != EXPECTED_CHECKS) + stop("expected ", EXPECTED_CHECKS, " checks, receipt says ", ran, + " -- a check was added or silently lost; update EXPECTED_CHECKS") + cat("R verifier ran complete.\n") + manuscript-build: runs-on: ubuntu-latest timeout-minutes: 25 diff --git a/analysis/verify_biofilm_depth_profile.R b/analysis/verify_biofilm_depth_profile.R new file mode 100644 index 0000000..f8242e4 --- /dev/null +++ b/analysis/verify_biofilm_depth_profile.R @@ -0,0 +1,275 @@ +#!/usr/bin/env Rscript +# ============================================================================ +# Verify the slab (biofilm-depth) geometry of biofilms_radiodialysis.R. +# +# This checks the R code that SHIPS. Its numpy sibling, +# analysis/verify_biofilm_depth_profile.py, checks an independent re-coding of +# the same scheme; run both -- agreement between them is worth more than +# either alone, and the numpy one runs where R is absent. +# +# Loads only the model expressions, skipping the library() calls and the Shiny +# app, so the stencil checks need no third-party packages at all. The +# time-integration check needs deSolve and reports SKIPPED without it rather +# than being represented as passed (AGENTS.md rule 2). +# +# WHAT IS CHECKED +# 1. slab_parms() refuses a defaulted X_total and names the gate. +# 2. The operator radiodialysis_rhs actually assembles equals an independent +# assembly, node by node, including the Robin row. +# 3. The steady slab profile matches the closed-form Robin solution. +# 4. penetration_depth() reproduces lambda_steady/lambda_transient = sqrt(6). +# 5. REGRESSION: the cylindrical path -- which ships and is cross-validated +# against biofilms_potts.jl:1153-1163 -- is byte-equivalent to the +# pre-change inline arithmetic across every row the refactor touched. +# 6. NEGATIVE CONTROLS: check 5 is re-run under three mutations of +# face_weights(), each of which must be caught. Without this the suite +# could pass on a stencil that ignored geometry (AGENTS.md rule 1) -- +# and an earlier draft of check 5 did exactly that, because in slab +# geometry w_plus == w_minus == 1 makes a swap invisible. +# +# WHAT IS *NOT* CHECKED +# Whether the model predicts a MEASURABLE gradient in a real biofilm. That +# is downstream of RADIODIALYSIS: BLOCKED (README.md:350-353): U scales with +# X_total and lambda ~ 1/sqrt(U), so a 20-100x correction to the biomass +# basis grows lambda by 4.5-10x and flattens the profile. Every X_total +# below is a TEST VALUE chosen to exercise the numerics. No profile number +# produced here may be quoted as a claim about a biofilm. +# ============================================================================ + +MODEL <- file.path(dirname(dirname(normalizePath( + sub("^--file=", "", grep("^--file=", commandArgs(FALSE), value = TRUE)[1])))), + "biofilms_radiodialysis.R") +if (!file.exists(MODEL)) MODEL <- "biofilms_radiodialysis.R" + +WANTED <- c("radiodialysis_rhs", "face_weights", "default_parms", + "slab_parms", "penetration_depth", "run_radiodialysis") +for (e in parse(MODEL)) { + if (is.call(e) && identical(as.character(e[[1]]), "<-") && + as.character(e[[2]]) %in% WANTED) eval(e, envir = globalenv()) +} +stopifnot(all(vapply(WANTED, exists, logical(1)))) + +failures <- 0L +skipped <- character(0) +checks_run <- 0L +report <- function(ok, msg, detail = "") { + checks_run <<- checks_run + 1L + if (!ok) failures <<- failures + 1L + cat(sprintf("[%s] %s%s\n", if (ok) "PASS" else "FAIL", msg, + if (nzchar(detail)) paste0(" (", detail, ")") else "")) + invisible(ok) +} + +# --------------------------------------------------------------------------- +# Mutation campaign, with the baseline guard built in. +# +# A mutation harness has TWO ways to measure nothing, and they look opposite: +# +# always-passes : no mutant is ever caught -- obviously broken, you notice. +# always-FAILS : every mutant is "caught" -- looks like perfect detection. +# +# The second is the dangerous one, because universal failure is the shape you +# are hoping for. It happened here: an earlier reference matrix omitted the +# D(w+ + w-)/dg^2 term in the Robin row, the baseline sat at 9.8e-01, and all +# three mutants read as caught while the harness measured nothing at all. +# +# So the baseline verdict GATES the mutant verdicts. If an unmutated run does +# not pass, no mutant result is reported as meaningful -- the campaign is +# declared invalid instead. Do not rely on remembering this. +# --------------------------------------------------------------------------- +mutation_campaign <- function(label, probe, tol, restore, mutants) { + base <- probe() + if (!(base < tol)) { + report(FALSE, sprintf("%s: MUTATION HARNESS INVALID", label), + sprintf("baseline residual %.2e is not < %.0e, so every mutant would + read as caught regardless of the code. No mutant verdict reported.", + base, tol)) + return(invisible(FALSE)) + } + cat(sprintf("\n %s -- baseline clean at %.2e; each mutation must be caught:\n", + label, base)) + for (nm in names(mutants)) { + on.exit(restore(), add = TRUE) + mutants[[nm]]() + r <- tryCatch(probe(), error = function(e) Inf) + restore() + report(r > tol, sprintf(" caught: %s", nm), sprintf("%.2e", r)) + } + report(probe() < tol, " baseline restored after mutation") +} + +# --- the operator radiodialysis_rhs assembles: rhs(c) = A c + b ------------- +operator_of <- function(p) { + Nr <- length(p$r_grid) + zero <- c(rep(0, Nr), rep(0, Nr), 1.0) + b <- radiodialysis_rhs(0, zero, p)[[1]][seq_len(Nr)] + A <- sapply(seq_len(Nr), function(i) { + y <- zero; y[i] <- 1.0 + radiodialysis_rhs(0, y, p)[[1]][seq_len(Nr)] - b + }) + list(A = A, b = b) +} + +# --- independent assembly from the pre-change inline arithmetic ------------- +# Derived from the ghost substitution, not guessed: substituting +# c_ghost = c[N-1] - 2 dr P (c[N] - c_ext)/D +# into the outer stencil gives the A[N,N] and b[N] below. An earlier draft +# omitted the D(w+ + w-)/dr^2 term here; the baseline then failed at 9.8e-01 +# and every mutant looked "caught" for the wrong reason. +reference <- function(p, geom) { + g <- p$r_grid; N <- length(g); dg <- g[2] - g[1]; D <- p$D_eff + U <- p$k_ads * p$X_total + p$k_red * p$X_red + bc <- if (identical(geom, "slab")) p$k_L else p$P0 + wp <- if (identical(geom, "slab")) rep(1, N) else (g + 0.5 * dg) / g + wm <- if (identical(geom, "slab")) rep(1, N) else (g - 0.5 * dg) / g + A <- matrix(0, N, N); b <- numeric(N) + for (i in 2:(N - 1)) { + A[i, i - 1] <- D * wm[i] / dg^2 + A[i, i] <- -(D * (wp[i] + wm[i]) / dg^2 + U) + A[i, i + 1] <- D * wp[i] / dg^2 + } + A[N, N - 1] <- D * (wp[N] + wm[N]) / dg^2 + A[N, N] <- -(D * (wp[N] + wm[N]) / dg^2 + 2 * wp[N] * bc / dg + U) + b[N] <- 2 * wp[N] * bc * p$c_ext / dg + list(A = A, b = b) +} + +residual <- function(p, geom) { + o <- operator_of(p); r <- reference(p, geom); rows <- 2:nrow(o$A) + max(abs(o$A[rows, ] - r$A[rows, ])) / max(abs(r$A[rows, ])) +} + +# --- 1. the refuse-to-default guard ---------------------------------------- +report(tryCatch({ slab_parms(); FALSE }, + error = function(e) grepl("RADIODIALYSIS: BLOCKED", + conditionMessage(e), fixed = TRUE)), + "slab_parms() refuses a defaulted X_total and names the gate") + +# --- 2-3. slab operator and steady profile --------------------------------- +L_f <- 0.01; D <- 1e-5; N <- 40; kL <- 1e2 +k_eff <- 0.056 * 0.001 / 0.006 +p <- slab_parms(Nr = N, L_f = L_f, X_total = 1.0, D_eff = D, k_L = kL) +p$k_ads <- k_eff; p$k_red <- 0; p$k_des <- 0; p$k_loss <- 0; p$X_red <- 0 +res <- residual(p, "slab") +report(res < 1e-13, "slab operator matches independent assembly", + sprintf("%.2e", res)) + +op <- operator_of(p) +c_num <- solve(op$A, -op$b) +lam <- sqrt(D / k_eff) +c_exact <- p$c_ext * cosh(p$r_grid / lam) / + (cosh(L_f / lam) + (D / (kL * lam)) * sinh(L_f / lam)) +err <- max(abs(c_num - c_exact) / c_exact) +report(err < 1e-5, "steady slab profile matches the analytic Robin solution", + sprintf("max rel err %.2e, lambda %.1f um, phi %.3f, c(0)/c(L) %.4f", + err, lam * 1e4, L_f / lam, c_num[1] / c_num[N])) + +# --- 4. the structural sqrt(6) --------------------------------------------- +p2 <- slab_parms(Nr = N, L_f = L_f, X_total = 1.0, D_eff = D) +ratio <- penetration_depth(p2, "steady") / penetration_depth(p2, "transient") +report(abs(ratio - sqrt(6)) < 1e-12, + "lambda_steady/lambda_transient = sqrt(6), independent of X_total", + sprintf("%.9f", ratio)) + +# --- 5. cylindrical regression --------------------------------------------- +cyl <- function() default_parms(Nr = 40, R = 1.0) +res_c <- residual(cyl(), "cylindrical") +report(res_c < 1e-13, + "cylindrical path byte-equivalent to the pre-change stencil", + sprintf("%.2e", res_c)) + +# --- 6. negative controls --------------------------------------------------- +orig <- face_weights +mutation_campaign( + label = "negative controls", + probe = function() residual(cyl(), "cylindrical"), + tol = 1e-13, + restore = function() face_weights <<- orig, + mutants = list( + "face_weights swapped" = function() + face_weights <<- function(g, gm) { + w <- orig(g, gm); list(w_plus = w$w_minus, w_minus = w$w_plus) }, + "w_plus[5] perturbed by 0.1%" = function() + face_weights <<- function(g, gm) { + w <- orig(g, gm); w$w_plus[5] <- w$w_plus[5] * 1.001; w }, + "weights forced to slab (= 1)" = function() + face_weights <<- function(g, gm) + list(w_plus = rep(1, length(g)), w_minus = rep(1, length(g))) + ) +) + +# The gate itself needs a negative control, or it is just another check that +# cannot fail. Deliberately break the reference the campaign measures against +# and confirm the campaign REFUSES rather than reporting three caught mutants. +local({ + good <- reference + reference <<- function(p, geom) { # drop the Robin diffusion term: + r <- good(p, geom); N <- nrow(r$A) # exactly the bug that occurred + g <- p$r_grid; dg <- g[2] - g[1] + wp <- (g[N] + 0.5 * dg) / g[N]; wm <- (g[N] - 0.5 * dg) / g[N] + r$A[N, N] <- r$A[N, N] + p$D_eff * (wp + wm) / dg^2 + r + } + before <- failures + invisible(capture.output(mutation_campaign( + "self-test", function() residual(cyl(), "cylindrical"), 1e-13, + function() NULL, list("noop" = function() NULL)))) + reference <<- good + fired <- failures > before + failures <<- before # do not count the deliberate break + report(fired, + " the baseline gate itself fires on a broken reference") +}) + +# --- 7. time integration (needs deSolve) ------------------------------------ +cat("\n") +if (!requireNamespace("deSolve", quietly = TRUE)) { + skipped <- c(skipped, "time integration (deSolve not installed)") + cat("[SKIP] time integration: deSolve not installed.", + "This is uncovered surface, not a pass.\n") +} else { + library(deSolve) + ps <- slab_parms(Nr = N, L_f = L_f, X_total = 1.0, D_eff = D, k_L = kL) + out <- run_radiodialysis(ps, t_end = 2e5, n_out = 50) + final <- out$c_mat[nrow(out$c_mat), ] + # full system: sorption chain on, so the steady sink is k_eff + ke <- (ps$k_ads * ps$X_total + ps$k_red * ps$X_red) * + ps$k_loss / (ps$k_des + ps$k_loss) + lm2 <- sqrt(ps$D_eff / ke) + ex <- ps$c_ext * cosh(ps$r_grid / lm2) / + (cosh(L_f / lm2) + (ps$D_eff / (kL * lm2)) * sinh(L_f / lm2)) + e2 <- max(abs(final - ex) / ex) + report(e2 < 1e-3, + "lsoda steady state matches the analytic Robin solution", + sprintf("max rel err %.2e at t=2e5 s", e2)) + report(abs(out$m_vec[length(out$m_vec)] - 1.0) < 1e-12, + "membrane integrity m stays 1 in slab geometry (radiation off)") +} + +# A skip is a question, not a pass (AGENTS.md rule 2). Never print a bare +# "ALL PASS" over a suite that did not run everything: the summary line is what +# gets quoted, and "215 passed, 8 skipped" is how a dead function shipped here. +verdict <- if (failures > 0L) { + sprintf("FAILURES: %d", failures) +} else if (length(skipped)) { + "PASSED WHAT RAN -- NOT A CLEAN RUN" +} else { + "ALL PASS" +} +cat(sprintf("\n%s (%d check%s run, %d failure%s, %d skipped)\n", verdict, + checks_run, if (checks_run == 1L) "" else "s", failures, + if (failures == 1L) "" else "s", length(skipped))) +for (s in skipped) cat(" UNCOVERED:", s, "\n") + +# Machine-readable receipt. CI must assert on THIS, not on the exit status: +# a job whose R setup step failed or was omitted can exit 0 having run nothing, +# which is the ALL-PASS-over-skip bug relocated from R into YAML. +args <- commandArgs(TRUE) +if (length(args) >= 2L && args[1] == "--report") { + writeLines(sprintf( + '{"checks_run": %d, "failures": %d, "skipped": %d, "skips": [%s], "complete": %s}', + checks_run, failures, length(skipped), + paste(sprintf('"%s"', skipped), collapse = ", "), + tolower(as.character(failures == 0L && !length(skipped)))), args[2]) +} +quit(status = if (failures == 0L && !length(skipped)) 0L else 1L) diff --git a/analysis/verify_biofilm_depth_profile.py b/analysis/verify_biofilm_depth_profile.py new file mode 100644 index 0000000..0ab3820 --- /dev/null +++ b/analysis/verify_biofilm_depth_profile.py @@ -0,0 +1,373 @@ +#!/usr/bin/env python3 +"""Verify the slab (biofilm-depth) stencil of ``biofilms_radiodialysis.R``. + +WHAT THIS COVERS, AND WHAT IT DOES NOT +-------------------------------------- +This script rebuilds the slab face weights in numpy and checks them against a +closed-form solution. ``biofilms_radiodialysis.R`` is what ships, and it is +covered separately by ``verify_biofilm_depth_profile.R`` (see the last section). +So this file establishes ONE of two independent codings of one scheme; the +shipping implementation is verified by its companion, not here. + +Do not report "the geometry change is verified against the analytic solution" +without saying which of the two codings was verified. + +THE CONTRACT BEING CHECKED +-------------------------- +Setting ds/dt = 0 in the immobile-phase equation gives + + s = U c / (k_des + k_loss), U = k_ads*X_total + k_red*X_red + +and substituting it back shows the sorption/desorption pair is a *reversible +buffer* that cancels. The only true steady sink is precipitation: + + k_eff = U * k_loss / (k_des + k_loss) + +so the steady mobile equation is D c'' = k_eff c, whose solution on 0 <= z <= L +with no-flux at the substratum and a Dirichlet bulk face is + + c(z) = c_b * cosh(z / lam) / cosh(L / lam), lam = sqrt(D / k_eff) + +Four assertions: + + 1. accuracy -- max relative error < 1e-5 on a 40-node grid + 2. order -- error falls ~4x per grid halving (second order) + 3. NEGATIVE CONTROL -- the *cylindrical* weights, run against this same slab + solution, must FAIL assertion 1. Without this the suite would pass on a + stencil that ignored geometry entirely (AGENTS.md rule 1). + 4. non-flat regime -- 1 and 2 repeated at phi = L/lam ~ 3.7, so the operator + is exercised where a gradient actually exists rather than only in the + nearly-flat default regime. + +WHAT THIS DOES *NOT* ESTABLISH +------------------------------ +Whether the model predicts a *measurable* gradient across a real film. That is +downstream of ``RADIODIALYSIS: BLOCKED`` (README.md:350-353): U is proportional +to X_total and lam ~ 1/sqrt(U), so a 20-100x correction to the biomass basis +grows lam by 4.5-10x and flattens the profile. The X_total values used here +are TEST VALUES chosen to exercise the numerics. They are not claims about a +biofilm and no profile number from this script may be quoted as one. + +Nor does it cover the transient, which is the regime an SECM profile would +actually be compared against; separation of variables would give an analytic +form for it. Assertion 4 covers the spatial operator in the non-flat regime. +A transient check would additionally cover the time integration. + +THE COMPANION CHECK ON THE R ITSELF +----------------------------------- +``analysis/verify_biofilm_depth_profile.R`` now covers the shipping R code -- +its assembled operator, its steady profile, its lsoda time integration, and a +regression proving the cylindrical path is unchanged -- with mutation-based +negative controls. Run both:: + + Rscript analysis/verify_biofilm_depth_profile.R + coupling/.venv/bin/python analysis/verify_biofilm_depth_profile.py + +They agree to the digit on the shared quantity (max rel err 2.31e-07 at N=40), +which is worth more than either alone: this file is an independent re-coding, +so agreement rules out a shared transcription error, and it still runs where R +is absent. Keep both. Neither supersedes the other. +""" + +from __future__ import annotations + +import argparse +import json +import math +from dataclasses import asdict, dataclass +from pathlib import Path +from typing import Any + +import numpy as np + +# Tolerances. +# +# A fixed absolute tolerance cannot serve here: the truncation error scales as +# (dz/lam)^2 = (phi/(N-1))^2, so any single number either fails a strong +# gradient or is vacuous for a weak one. The bound is therefore on the +# NORMALISED constant C = err * ((N-1)/phi)^2, which is N-independent (that +# is what "second order" means) and depends only weakly on phi. Measured: +# +# phi 0.30 1.00 2.00 3.70 6.00 +# C 0.0036 0.0317 0.0803 0.1539 0.2501 -> saturating near 1/4 +# +# 0.3 bounds all of them. It is deliberately not tighter: the sharp check is +# the refinement ratio, which sits at 4.05-4.21 and breaks on any stencil error. +ERROR_CONSTANT_MAX = 0.3 +ORDER_TOL = 3.6 # second order is 4.0; allow for the coarsest-grid pre-asymptote +GRID_N = 40 + + +@dataclass(frozen=True) +class Kinetics: + """Rate constants as written in biofilms_radiodialysis.R:141-148.""" + + k_ads: float = 0.05 + k_red: float = 0.02 + k_des: float = 0.005 + k_loss: float = 0.001 + X_red: float = 0.3 + X_total: float = 1.0 # TEST VALUE -- gated, see module docstring + + def validate(self) -> None: + for name, value in asdict(self).items(): + if not math.isfinite(value) or value < 0.0: + raise ValueError(f"{name} must be finite and non-negative, got {value!r}") + if self.k_des + self.k_loss <= 0.0: + raise ValueError("k_des + k_loss must be positive for a steady state to exist") + + @property + def uptake(self) -> float: + """U, the transient sink (s^-1).""" + self.validate() + return self.k_ads * self.X_total + self.k_red * self.X_red + + @property + def k_eff(self) -> float: + """The only true steady sink: sorption/desorption cancels.""" + return self.uptake * self.k_loss / (self.k_des + self.k_loss) + + +def face_weights(grid: np.ndarray, geom: str) -> tuple[np.ndarray, np.ndarray]: + """Port of face_weights() in biofilms_radiodialysis.R. + + Node 0 is NaN for the same reason it is NA in the R source: it uses the + reflected-ghost form and reads no metric weight. NaN rather than a + plausible number so that indexing it propagates loudly. + """ + if geom not in ("cylindrical", "slab"): + raise ValueError(f"geom must be 'cylindrical' or 'slab', got {geom!r}") + dz = grid[1] - grid[0] + if geom == "slab": + w_plus = np.ones_like(grid) + w_minus = np.ones_like(grid) + else: + with np.errstate(divide="ignore", invalid="ignore"): + w_plus = (grid + 0.5 * dz) / grid + w_minus = (grid - 0.5 * dz) / grid + w_plus[0] = np.nan + w_minus[0] = np.nan + return w_plus, w_minus + + +def steady_profile(n: int, length: float, diffusivity: float, sink: float, + c_bulk: float, geom: str) -> tuple[np.ndarray, np.ndarray]: + """Solve D*L(c) = sink*c on the given grid, Dirichlet at the outer face. + + Uses the same operator as radiodialysis_rhs: reflected ghost at node 0 + (factor 2, identical for the cylindrical axis and the slab substratum), + face-weighted interior. The Dirichlet outer face is the k_L -> infinity + limit of the code's Robin condition, which is what has a closed form. + """ + if n < 4: + raise ValueError(f"need at least 4 nodes, got {n}") + grid = np.linspace(0.0, length, n) + dz = grid[1] - grid[0] + w_plus, w_minus = face_weights(grid, geom) + operator = np.zeros((n, n)) + rhs = np.zeros(n) + + # node 0: reflected ghost c[-1] = c[1] -> 2*(c[1] - c[0])/dz^2 + operator[0, 0] = -(2.0 * diffusivity / dz**2 + sink) + operator[0, 1] = 2.0 * diffusivity / dz**2 + + for i in range(1, n - 1): + operator[i, i - 1] = diffusivity * w_minus[i] / dz**2 + operator[i, i] = -(diffusivity * (w_plus[i] + w_minus[i]) / dz**2 + sink) + operator[i, i + 1] = diffusivity * w_plus[i] / dz**2 + + operator[-1, -1] = 1.0 + rhs[-1] = c_bulk + return np.linalg.solve(operator, rhs), grid + + +def analytic(grid: np.ndarray, length: float, lam: float, c_bulk: float) -> np.ndarray: + return c_bulk * np.cosh(grid / lam) / np.cosh(length / lam) + + +def max_rel_error(n: int, length: float, diffusivity: float, sink: float, + c_bulk: float, geom: str) -> float: + numeric, grid = steady_profile(n, length, diffusivity, sink, c_bulk, geom) + exact = analytic(grid, length, math.sqrt(diffusivity / sink), c_bulk) + return float(np.max(np.abs(numeric - exact) / exact)) + + +def check_cylindrical_refactor_equivalence(n: int = 40, radius: float = 1.0, + diffusivity: float = 1e-3) -> dict[str, Any]: + """Regression guard: the weight form must equal the original inline stencil. + + The R change hoisted r_plus/r_i and r_minus/r_i out of radiodialysis_rhs + into face_weights(). That touches the CYLINDRICAL path, which ships and is + cross-validated against biofilms_potts.jl:1153-1163. The slab assertions + cannot protect it: there w_plus == w_minus == 1, so a weight-ordering error + is invisible. A mutation test confirmed exactly that -- swapping the two + weights survived every other check in this file. + + This compares the refactored form against the pre-change arithmetic + + D * (r_plus*(c[i+1]-c[i]) - r_minus*(c[i]-c[i-1])) / (r_i * dr^2) + + node by node, and is asymmetric under a swap. + """ + grid = np.linspace(0.0, radius, n) + dr = grid[1] - grid[0] + w_plus, w_minus = face_weights(grid, "cylindrical") + # A profile with structure at every node, so no term cancels by accident. + c = np.cos(3.0 * grid / radius) + 0.5 * grid / radius + + original = np.empty(n - 2) + refactored = np.empty(n - 2) + for i in range(1, n - 1): + r_i = grid[i] + r_plus, r_minus = r_i + 0.5 * dr, r_i - 0.5 * dr + original[i - 1] = diffusivity * ( + r_plus * (c[i + 1] - c[i]) - r_minus * (c[i] - c[i - 1]) + ) / (r_i * dr**2) + refactored[i - 1] = diffusivity * ( + w_plus[i] * (c[i + 1] - c[i]) - w_minus[i] * (c[i] - c[i - 1]) + ) / dr**2 + + scale = np.max(np.abs(original)) + residual = float(np.max(np.abs(original - refactored)) / scale) + + # Negative control: the same comparison with the weights swapped must NOT + # come out equal, or this guard proves nothing. + swapped = np.empty(n - 2) + for i in range(1, n - 1): + swapped[i - 1] = diffusivity * ( + w_minus[i] * (c[i + 1] - c[i]) - w_plus[i] * (c[i] - c[i - 1]) + ) / dr**2 + swap_residual = float(np.max(np.abs(original - swapped)) / scale) + + # The baseline verdict GATES the control verdict. A negative control can + # fail to measure anything in two opposite ways: never firing (obviously + # broken) or ALWAYS firing, which looks like perfect detection and is the + # dangerous one. A wrong reference makes every mutant read as caught -- it + # happened in the R harness, where an omitted Robin term put the baseline + # at 9.8e-01 and three mutants read as caught over a harness measuring + # nothing. So if the baseline does not match, the swap control is reported + # as INVALID rather than as a pass. + baseline_ok = residual < 1e-13 + checks = {"matches_pre_change_stencil": baseline_ok} + if baseline_ok: + checks["swap_control_differs"] = swap_residual > 1e-6 + else: + checks["swap_control_INVALID_baseline_failed"] = False + + return { + "case": "cylindrical_refactor_equivalence", + "note": "protects the shipping cylindrical path and its Julia cross-validation", + "residual": residual, + "swap_control_residual": swap_residual, + "baseline_gated": True, + "checks": checks, + "passed": all(checks.values()), + } + + +def run_case(name: str, length: float, diffusivity: float, sink: float, + note: str, c_bulk: float = 1.0) -> dict[str, Any]: + lam = math.sqrt(diffusivity / sink) + phi = length / lam + errors = {n: max_rel_error(n, length, diffusivity, sink, c_bulk, "slab") + for n in (GRID_N // 2, GRID_N, GRID_N * 2, GRID_N * 4)} + ratios = [errors[n] / errors[2 * n] for n in sorted(errors)[:-1]] + control = max_rel_error(GRID_N, length, diffusivity, sink, c_bulk, "cylindrical") + profile, _ = steady_profile(GRID_N, length, diffusivity, sink, c_bulk, "slab") + + def constant(err: float, n: int) -> float: + return err * ((n - 1) / phi) ** 2 + + bound = ERROR_CONSTANT_MAX * (phi / (GRID_N - 1)) ** 2 + checks = { + "accuracy": constant(errors[GRID_N], GRID_N) < ERROR_CONSTANT_MAX, + "second_order": all(r > ORDER_TOL for r in ratios), + "negative_control_fails_slab_bound": control > bound, + } + return { + "case": name, + "note": note, + "thiele_modulus": phi, + "lambda_um": lam * 1e4, + "profile_ratio_c0_over_cL": float(profile[0] / profile[-1]), + "errors": {str(k): v for k, v in errors.items()}, + "error_constant": constant(errors[GRID_N], GRID_N), + "refinement_ratios": ratios, + "slab_bound_at_N": bound, + "negative_control_error": control, + "checks": checks, + "passed": all(checks.values()), + } + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument("--report", type=Path, default=None, + help="write the machine-readable report here") + args = parser.parse_args() + + length = 0.01 # 100 um film + diffusivity = 1e-5 + kin = Kinetics() + cases = [ + run_case( + "default_regime", length, diffusivity, kin.k_eff, + note="k_eff from the file's own rate constants at X_total=1.0 " + "(a TEST VALUE: the biomass basis is gated). Nearly flat.", + ), + # Assertion 4. The sink is set directly from a target phi rather than + # by inventing a biomass number: this case exists to exercise the + # operator where a gradient exists, and corresponds to no parameter set. + run_case( + "strong_gradient", length, diffusivity, + sink=diffusivity / (length / 3.7) ** 2, + note="phi = 3.7 imposed directly on the operator. NOT a physical " + "parameter set and not a prediction about any biofilm.", + ), + check_cylindrical_refactor_equivalence(), + ] + + report = { + "tolerances": {"error_constant_max": ERROR_CONSTANT_MAX, + "order_ratio": ORDER_TOL, "N": GRID_N}, + "k_eff_over_U": Kinetics().k_eff / Kinetics().uptake, + "covers": "the numpy specification of the slab stencil", + "does_not_cover": [ + "biofilms_radiodialysis.R itself (covered by verify_biofilm_depth_profile.R)", + "the transient regime", + "whether any gradient is measurable (downstream of RADIODIALYSIS: BLOCKED)", + ], + "cases": cases, + "passed": all(c["passed"] for c in cases), + } + + for case in cases: + print(f"[{'PASS' if case['passed'] else 'FAIL'}] {case['case']}") + if "thiele_modulus" in case: + print(f" phi={case['thiele_modulus']:.3f} " + f"lambda={case['lambda_um']:.1f}um " + f"c(0)/c(L)={case['profile_ratio_c0_over_cL']:.4f}") + print(f" err(N={GRID_N})={case['errors'][str(GRID_N)]:.2e} " + f"C={case['error_constant']:.4f}/{ERROR_CONSTANT_MAX} " + f"ratios={['%.2f' % r for r in case['refinement_ratios']]} " + f"control={case['negative_control_error']:.2e} " + f"(bound {case['slab_bound_at_N']:.2e})") + else: + print(f" residual={case['residual']:.2e} " + f"swap_control={case['swap_control_residual']:.2e}") + for check, ok in case["checks"].items(): + if not ok: + print(f" FAILED: {check}") + + print(f"\nk_eff/U = {report['k_eff_over_U']:.6f} = 1/{1/report['k_eff_over_U']:.0f} " + f"(structural, independent of X_total)") + print("COVERS: the numpy specification, NOT biofilms_radiodialysis.R.") + + if args.report: + args.report.write_text(json.dumps(report, indent=2) + "\n") + + return 0 if report["passed"] else 1 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/biofilms_radiodialysis.R b/biofilms_radiodialysis.R index c706a85..252bce1 100644 --- a/biofilms_radiodialysis.R +++ b/biofilms_radiodialysis.R @@ -21,6 +21,14 @@ # # Zero-flux (symmetry) at r = 0. # +# TWO GEOMETRIES, one stencil. default_parms() gives the cylindrical +# reactor above (biofilm lumped into a uniform scalar sink over a 1 cm +# radius). slab_parms() gives depth across a single biofilm — z = 0 +# substratum, z = L_f bulk-liquid face — which is the only one of the two +# that can be compared with micron-scale in-film measurements. Geometry +# enters solely through the face weights built by face_weights(); the +# axis/substratum node is identical in both and is unchanged. +# # Spatial discretisation: finite-volume method of lines (Nr cells). # Time integration: deSolve::ode (LSODA adaptive stiff solver). # @@ -57,8 +65,15 @@ radiodialysis_rhs <- function(t, y, parms) { # Cumulative absorbed dose at the membrane (Gy), linear in t D_cum <- Ddot_R * t - # Radiation-driven effective permeability (Lara et al. 2023, Eq. 6) - P_eff <- P0 * exp(alpha_P * D_cum) + # Outer-face transfer coefficient (cm s⁻¹). The producer declares which + # physical quantity this is; the stencil below must not assume one. + # cylindrical : radiation-evolving Nafion permeability P_eff(t) + # (Lara et al. 2023, Eq. 6) + # slab : external liquid-film mass-transfer coefficient k_L, + # time-invariant. This is NOT P_eff and must not be named + # or reported as one — commit e4021c0 withdrew P_eff from + # every producer for carrying a fabricated unit. + bc_coef <- if (identical(geom, "slab")) k_L else P0 * exp(alpha_P * D_cum) # -------------------------------------------------------- # (3) Membrane damage ODE @@ -71,43 +86,43 @@ radiodialysis_rhs <- function(t, y, parms) { uptake_rate <- k_ads * X_total + k_red * X_red # s⁻¹ # -------------------------------------------------------- - # (1) Mobile species — cylindrical diffusion + reaction + # (1) Mobile species — diffusion + reaction (geometry via w_plus/w_minus) # -------------------------------------------------------- dc_dt <- numeric(Nr) - # i = 1 (r = 0): L'Hôpital limit → ∂c/∂t = 2 D_eff ∂²c/∂r² + # i = 1: identical in both geometries, and unchanged by the slab work. + # cylindrical r = 0, symmetry, L'Hôpital limit → 2 D_eff ∂²c/∂r² + # slab z = 0, substratum, no-flux ghost c₀ = c₂ → 2 D_eff ∂²c/∂z² + # The reflected ghost is c₀ = c₂, NOT c₀ = c₁; the latter drops the whole + # scheme from second order to first. Reads no metric weight (both are NA). dc_dt[1] <- D_eff * 2.0 * (c_vec[2] - c_vec[1]) / dr^2 + (-uptake_rate * c_vec[1] + k_des * s_vec[1]) - # i = 2 .. Nr-1: interior finite-volume stencil + # i = 2 .. Nr-1: interior finite-volume stencil. + # Geometry enters only through the face weights w_plus / w_minus, built + # once in the parameter constructor: + # cylindrical w± = (r ± dr/2) / r slab w± = 1 if (Nr > 2) { for (i in 2:(Nr - 1)) { - r_i <- r_grid[i] - r_plus <- r_i + 0.5 * dr - r_minus <- r_i - 0.5 * dr - diff_cyl <- D_eff * - (r_plus * (c_vec[i + 1] - c_vec[i]) - - r_minus * (c_vec[i] - c_vec[i - 1])) / - (r_i * dr^2) - dc_dt[i] <- diff_cyl + + diff_op <- D_eff * + (w_plus[i] * (c_vec[i + 1] - c_vec[i]) - + w_minus[i] * (c_vec[i] - c_vec[i - 1])) / dr^2 + dc_dt[i] <- diff_op + (-uptake_rate * c_vec[i] + k_des * s_vec[i]) } } - # i = Nr (r = R): Robin BC via ghost-point (Donnan / Nafion membrane) - # ghost: c[Nr+1] = c[Nr-1] - 2 dr P_eff (c[Nr] - c_ext) / D_eff + # i = Nr: Robin BC via ghost-point. + # cylindrical r = R, the Donnan / Nafion membrane + # slab z = L_f, the biofilm / bulk-liquid interface + # ghost: c[Nr+1] = c[Nr-1] - 2 dr bc_coef (c[Nr] - c_ext) / D_eff { - i <- Nr - r_i <- r_grid[i] - r_plus <- r_i + 0.5 * dr - r_minus <- r_i - 0.5 * dr c_ghost <- c_vec[Nr - 1] - - 2.0 * dr * P_eff * (c_vec[Nr] - c_ext) / D_eff - diff_cyl <- D_eff * - (r_plus * (c_ghost - c_vec[Nr]) - - r_minus * (c_vec[Nr] - c_vec[Nr - 1])) / - (r_i * dr^2) - dc_dt[Nr] <- diff_cyl + + 2.0 * dr * bc_coef * (c_vec[Nr] - c_ext) / D_eff + diff_op <- D_eff * + (w_plus[Nr] * (c_ghost - c_vec[Nr]) - + w_minus[Nr] * (c_vec[Nr] - c_vec[Nr - 1])) / dr^2 + dc_dt[Nr] <- diff_op + (-uptake_rate * c_vec[Nr] + k_des * s_vec[Nr]) } @@ -127,12 +142,42 @@ radiodialysis_rhs <- function(t, y, parms) { # literature-consistent estimates for low-level nuclear waste # remediation context. # ------------------------------------------------------------ + +# ------------------------------------------------------------ +# Diffusion face weights. The ONLY thing geometry changes. +# +# cylindrical (1/r) ∂/∂r (r D ∂c/∂r) -> w± = (r ± dr/2) / r +# slab (Cartesian) ∂/∂z (D ∂c/∂z) -> w± = 1 +# +# Node 1 is deliberately NA: it uses the reflected-ghost form +# (c₀ = c₂ -> factor 2), which is identical for the cylindrical +# axis (L'Hôpital) and the slab substratum (no-flux), and reads +# no metric weight. NA so that a future reader who *does* index +# it fails loudly rather than picking up a silently wrong number. +# ------------------------------------------------------------ +face_weights <- function(r_grid, geom) { + stopifnot(geom %in% c("cylindrical", "slab")) + dr <- r_grid[2] - r_grid[1] + if (identical(geom, "slab")) { + w_plus <- rep(1.0, length(r_grid)); w_minus <- rep(1.0, length(r_grid)) + } else { + w_plus <- (r_grid + 0.5 * dr) / r_grid + w_minus <- (r_grid - 0.5 * dr) / r_grid + } + w_plus[1] <- NA_real_; w_minus[1] <- NA_real_ + list(w_plus = w_plus, w_minus = w_minus) +} + default_parms <- function(Nr = 40, R = 1.0) { r_grid <- seq(0, R, length.out = Nr) + w <- face_weights(r_grid, "cylindrical") list( # --- Spatial grid --- r_grid = r_grid, R = R, + geom = "cylindrical", + w_plus = w$w_plus, + w_minus = w$w_minus, # --- Transport --- D_eff = 1e-3, # effective diffusivity (cm² s⁻¹), Table 2 range 1e-4..1e-2 @@ -162,6 +207,99 @@ default_parms <- function(Nr = 40, R = 1.0) { ) } +# ------------------------------------------------------------ +# Slab preset: depth across a biofilm, NOT radius across a reactor. +# +# z = 0 substratum (no-flux) +# z = L_f bulk-liquid face (Robin, external mass transfer k_L) +# +# Same three equations, same stencil, Cartesian face weights. Written to +# be comparable with micron-scale in-film measurements (e.g. SECM metal +# profiles), which the cylindrical preset cannot produce: there the biofilm +# is a spatially uniform scalar sink over a 1 cm reactor radius. +# +# Three semantics are declared here rather than inherited (AGENTS.md rule 4): +# * the Robin coefficient is k_L, a liquid-film mass-transfer coefficient, +# NOT the membrane permeability P_eff; +# * the radiation terms are switched off — there is no membrane at a +# biofilm/liquid interface, so k_dam = Ddot_R = 0 and m(t) stays 1; +# * D_eff is a MOLECULAR diffusivity. The cylindrical preset's 1e-3 cm² s⁻¹ +# is ~100x the self-diffusivity of water: it is a reactor-scale dispersion +# coefficient and is meaningless on a 100 µm domain. +# +# X_total has NO default on purpose. See the stop() below. +# ------------------------------------------------------------ +slab_parms <- function(Nr = 40, L_f = 0.01, X_total, + D_eff = 1e-5, k_L = 1e-3) { + if (missing(X_total)) { + stop( + "slab_parms(): X_total must be stated explicitly, not defaulted.\n", + " The biofilm dry-mass basis is exactly what RADIODIALYSIS: BLOCKED\n", + " names (README.md:350-353). It sets the uptake rate\n", + " U = k_ads*X_total + k_red*X_red,\n", + " and the penetration depth lambda = sqrt(D_eff/k_eff) scales as\n", + " 1/sqrt(U), so whether ANY gradient is resolvable across L_f is\n", + " downstream of that gate. Pass a value and label it a test value.", + call. = FALSE + ) + } + z_grid <- seq(0, L_f, length.out = Nr) + w <- face_weights(z_grid, "slab") + list( + # --- Spatial grid: depth, not radius --- + r_grid = z_grid, # name kept so run_radiodialysis() is reused unchanged + R = L_f, # film thickness (cm) + geom = "slab", + w_plus = w$w_plus, + w_minus = w$w_minus, + + # --- Transport --- + D_eff = D_eff, # molecular diffusivity in the film (cm² s⁻¹) + + # --- Biosorption / bioreduction (unchanged from the cylindrical preset) --- + k_ads = 0.05, + k_red = 0.02, + k_des = 0.005, + k_loss = 0.001, # the ONLY true sink at steady state; see note below + + # --- Biomass --- + X_total = X_total, # caller-declared; gated, see stop() above + X_red = 0.3, + + # --- Outer face: liquid-film mass transfer, NOT membrane permeability --- + k_L = k_L, # (cm s⁻¹). P0 / alpha_P are deliberately absent. + + # --- Radiation: off. No membrane exists at a biofilm/liquid interface --- + k_dam = 0.0, + Ddot_R = 0.0, + + # --- Boundary --- + c_ext = 1.0 # bulk-liquid concentration (normalised) + ) +} + +# ------------------------------------------------------------ +# Steady-state penetration depth. +# +# Setting ds/dt = 0 gives s = U c / (k_des + k_loss), and substituting into +# the mobile equation shows the sorption/desorption pair is a REVERSIBLE +# BUFFER that cancels. The only true steady sink is precipitation: +# +# k_eff = U * k_loss / (k_des + k_loss) +# +# With the file's constants k_eff/U = 1/6 exactly, so lambda_steady is +# sqrt(6) ~ 2.45x the transient value. This ratio is structural and does +# not depend on X_total; the absolute lambda does. +# ------------------------------------------------------------ +penetration_depth <- function(parms, phase = c("steady", "transient")) { + phase <- match.arg(phase) + with(parms, { + U <- k_ads * X_total + k_red * X_red + k <- if (phase == "steady") U * k_loss / (k_des + k_loss) else U + sqrt(D_eff / k) + }) +} + # ------------------------------------------------------------ # Run solver: returns list with times, c-matrix, s-matrix, m-vec # ------------------------------------------------------------ From 23d4533707cc42d5fe6f89cd15368f3b07150cdd Mon Sep 17 00:00:00 2001 From: aurascoper Date: Fri, 28 Aug 2026 20:38:04 -0500 Subject: [PATCH 02/13] Install deSolve into the R that actually runs the verifier The apt step exited 0 and deSolve was still missing: setup-r@v2 provides its own R build reading R_LIBS_USER, while r-cran-desolve installs into the distro R's site-library, which that R never looks at. CI reported PASSED WHAT RAN -- NOT A CLEAN RUN (11 checks run, 0 failures, 1 skipped) UNCOVERED: time integration (deSolve not installed) and exited 1, so the job went red rather than green over an uncovered time integration. That is rule 2 behaving correctly; the bug is upstream of it, in the install. Install through R instead, with use-public-rspm for a binary rather than a Fortran source build. stopifnot() on the namespace because install.packages() only WARNS on an unavailable package, which would hand the verifier the same silent absence one step later (rule 3). Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/model-contracts.yml | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/.github/workflows/model-contracts.yml b/.github/workflows/model-contracts.yml index d35b9fd..b120508 100644 --- a/.github/workflows/model-contracts.yml +++ b/.github/workflows/model-contracts.yml @@ -87,6 +87,7 @@ jobs: - uses: r-lib/actions/setup-r@v2 with: r-version: "release" + use-public-rspm: true - name: Parse all top-level R models without executing Shiny apps shell: Rscript {0} run: | @@ -99,8 +100,18 @@ jobs: # parse() is not execution. It passes on a file whose every numerical # result is wrong, so this job was green over code that ran nowhere. + # Install through R, not apt. setup-r@v2 provides its own R build and + # reads R_LIBS_USER; r-cran-desolve installs into the DISTRO R's + # site-library, which that R never looks at. So the apt step exited 0 + # with deSolve still missing, and the verifier correctly refused to pass + # over the resulting skip. stopifnot() because install.packages() only + # WARNS on an unavailable package (rule 3: refuse, do not default to pass). - name: Install deSolve - run: sudo apt-get update && sudo apt-get install -y r-cran-desolve + shell: Rscript {0} + run: | + install.packages("deSolve") + stopifnot(requireNamespace("deSolve", quietly = TRUE)) + cat("deSolve", as.character(packageVersion("deSolve")), "\n") - name: Run the radiodialysis depth-geometry verifier run: Rscript analysis/verify_biofilm_depth_profile.R --report r-verify.json From b31faa601bf5d9bcc255ca3cdbefa6f691313469 Mon Sep 17 00:00:00 2001 From: aurascoper Date: Fri, 28 Aug 2026 20:45:56 -0500 Subject: [PATCH 03/13] Close Codex's three on #19, one of them in the gate itself P1, .github/workflows/model-contracts.yml. The receipt assertion regex- scraped three counters and never read the producer's own completeness declaration, so {"checks_run":14,"failures":0,"skipped":0,"complete":false} satisfied every check. A gate accepting a self-declared incomplete run is rule 3 in the one place whose whole job is to refuse that. Now parsed with jsonlite, `complete` required and isTRUE, counters required to be integers, and a missing key refused by name. Exercised against seven fixtures: valid passes, and complete=false, absent complete, count drift, a skip, a failure, truncated JSON and an absent file each refuse with their own message. P2, biofilms_radiodialysis.R. `if (geom == "slab") k_L else P_eff` read every unrecognised value as cylindrical. face_weights() validates, but the weights reach radiodialysis_rhs() precomputed in parms, so nothing validated there. Now switch() with an unnamed stop() default. Check 14 asserts both halves: an unknown geom errors, and both supported values still run, since a dispatch refusing everything would satisfy the first half alone. Confirmed to bite by restoring the if/else, which fails it at refused=FALSE with both positive halves still TRUE. P2, the numpy re-coding ran nowhere. It and its negative controls could regress without failing any check, while the PR leaned on cross-implementation agreement. Now run in radiodialysis-stability, which already installs numpy, and receipt-asserted on case count and per-case verdicts for the same reason as the R one. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/model-contracts.yml | 63 +++++++++++++++++++++---- analysis/verify_biofilm_depth_profile.R | 21 +++++++++ biofilms_radiodialysis.R | 14 +++++- 3 files changed, 87 insertions(+), 11 deletions(-) diff --git a/.github/workflows/model-contracts.yml b/.github/workflows/model-contracts.yml index b120508..c21539f 100644 --- a/.github/workflows/model-contracts.yml +++ b/.github/workflows/model-contracts.yml @@ -51,6 +51,32 @@ jobs: fi python analysis/verify_radiodialysis_stability.py \ --output artifacts/radiodialysis-stability.json $extra + # The numpy re-coding ran nowhere (Codex P2 on #19), so it and its + # negative controls could regress without failing any check -- while the + # PR claimed cross-implementation agreement on its strength. It lives + # here because this job already installs numpy. Receipt-asserted for the + # same reason as the R one: exit status alone cannot tell a full run from + # a run that quietly lost a case. + - name: Verify the biofilm depth stencil (numpy re-coding) + run: python analysis/verify_biofilm_depth_profile.py --report py-verify.json + + - name: Assert the numpy verifier covered every case + run: | + python - <<'PY' + import json, sys + EXPECTED_CASES = 3 + r = json.load(open("py-verify.json")) + names = [c["case"] for c in r["cases"]] + print(f"receipt: {len(names)} cases, passed={r['passed']}: {', '.join(names)}") + if len(names) != EXPECTED_CASES: + sys.exit(f"expected {EXPECTED_CASES} cases, receipt has {len(names)} " + "-- a case was added or silently lost; update EXPECTED_CASES") + bad = [c["case"] for c in r["cases"] if not c["passed"]] + if bad or not r["passed"]: + sys.exit(f"numpy verifier reported failures: {bad}") + print("numpy verifier ran complete.") + PY + - name: Summarize the time-step diagnostic if: always() run: | @@ -109,9 +135,11 @@ jobs: - name: Install deSolve shell: Rscript {0} run: | - install.packages("deSolve") - stopifnot(requireNamespace("deSolve", quietly = TRUE)) - cat("deSolve", as.character(packageVersion("deSolve")), "\n") + install.packages(c("deSolve", "jsonlite")) + for (pkg in c("deSolve", "jsonlite")) { + stopifnot(requireNamespace(pkg, quietly = TRUE)) + cat(pkg, as.character(packageVersion(pkg)), "\n") + } - name: Run the radiodialysis depth-geometry verifier run: Rscript analysis/verify_biofilm_depth_profile.R --report r-verify.json @@ -121,18 +149,33 @@ jobs: # a job that exits 0 having skipped is the ALL-PASS-over-skip bug moved # from R into YAML. EXPECTED_CHECKS must be bumped deliberately when # checks are added, so silently losing one fails the build. + # + # Parsed as JSON, and `complete` is REQUIRED (Codex P1 on #19). The + # previous version regex-scraped the three counters and never read the + # producer's own completeness declaration, so a receipt reading + # {"checks_run":14,...,"complete":false} satisfied every assertion -- + # a gate accepting a self-declared incomplete run is rule 3 in the one + # place whose entire job is to refuse that. - name: Assert the verifier actually ran everything shell: Rscript {0} run: | - EXPECTED_CHECKS <- 13L + EXPECTED_CHECKS <- 14L if (!file.exists("r-verify.json")) stop("no receipt: the verifier did not run to completion") - j <- readLines("r-verify.json", warn = FALSE) - num <- function(k) as.integer(sub( - paste0('.*"', k, '": *([0-9]+).*'), "\\1", j)) - ran <- num("checks_run"); fail <- num("failures"); skip <- num("skipped") - cat(sprintf("receipt: %d run, %d failed, %d skipped\n", ran, fail, skip)) - if (skip > 0L) stop("verifier skipped ", skip, " check(s): ", j) + r <- jsonlite::fromJSON("r-verify.json", simplifyVector = TRUE) + for (k in c("checks_run", "failures", "skipped", "complete")) + if (is.null(r[[k]])) stop("receipt is missing '", k, "'") + if (!isTRUE(r$complete)) + stop("receipt does not declare itself complete: complete=", + format(r$complete)) + ran <- as.integer(r$checks_run); fail <- as.integer(r$failures) + skip <- as.integer(r$skipped) + if (anyNA(c(ran, fail, skip))) stop("receipt counters are not integers") + cat(sprintf("receipt: %d run, %d failed, %d skipped, complete=%s\n", + ran, fail, skip, r$complete)) + if (skip > 0L) + stop("verifier skipped ", skip, " check(s): ", + paste(unlist(r$skips), collapse = "; ")) if (fail > 0L) stop("verifier reported ", fail, " failure(s)") if (ran != EXPECTED_CHECKS) stop("expected ", EXPECTED_CHECKS, " checks, receipt says ", ran, diff --git a/analysis/verify_biofilm_depth_profile.R b/analysis/verify_biofilm_depth_profile.R index f8242e4..c0418a2 100644 --- a/analysis/verify_biofilm_depth_profile.R +++ b/analysis/verify_biofilm_depth_profile.R @@ -145,6 +145,27 @@ report(tryCatch({ slab_parms(); FALSE }, conditionMessage(e), fixed = TRUE)), "slab_parms() refuses a defaulted X_total and names the gate") +# --- 1b. the geometry dispatch refuses what it does not recognise ---------- +# Codex P2 on #19. The old `if (geom == "slab") k_L else P_eff` read every +# unrecognised value as cylindrical. Both halves are asserted: an unknown +# geom must ERROR, and BOTH supported values must still run -- a dispatch that +# refuses everything would satisfy the first half alone. +local({ + p_bad <- slab_parms(X_total = 1.0); p_bad$geom <- "slabb" + y <- c(rep(0, length(p_bad$r_grid)), rep(0, length(p_bad$r_grid)), 1) + refused <- tryCatch({ radiodialysis_rhs(0, y, p_bad); FALSE }, + error = function(e) grepl("unsupported geom", + conditionMessage(e), fixed = TRUE)) + accepts <- vapply(list(slab_parms(X_total = 1.0), default_parms()), function(p) { + yy <- c(rep(0, length(p$r_grid)), rep(0, length(p$r_grid)), 1) + tryCatch({ radiodialysis_rhs(0, yy, p); TRUE }, error = function(e) FALSE) + }, logical(1)) + report(refused && all(accepts), + "geom dispatch refuses an unknown value and still accepts both known ones", + sprintf("refused=%s slab=%s cylindrical=%s", + refused, accepts[1], accepts[2])) +}) + # --- 2-3. slab operator and steady profile --------------------------------- L_f <- 0.01; D <- 1e-5; N <- 40; kL <- 1e2 k_eff <- 0.056 * 0.001 / 0.006 diff --git a/biofilms_radiodialysis.R b/biofilms_radiodialysis.R index 252bce1..53ef253 100644 --- a/biofilms_radiodialysis.R +++ b/biofilms_radiodialysis.R @@ -73,7 +73,19 @@ radiodialysis_rhs <- function(t, y, parms) { # time-invariant. This is NOT P_eff and must not be named # or reported as one — commit e4021c0 withdrew P_eff from # every producer for carrying a fabricated unit. - bc_coef <- if (identical(geom, "slab")) k_L else P0 * exp(alpha_P * D_cum) + # switch() with an unnamed default, NOT `if slab else cylindrical`: the + # else form treats EVERY unrecognised value as cylindrical, so a misspelt + # or future geom silently acquires a scientific interpretation instead of + # refusing (AGENTS.md rule 3). face_weights() validates, but the weights + # reach this function precomputed in parms, so it never runs here. + bc_coef <- switch( + geom, + slab = k_L, + cylindrical = P0 * exp(alpha_P * D_cum), + stop("radiodialysis_rhs(): unsupported geom ", sQuote(geom), + ". Supported: 'slab', 'cylindrical'. An unrecognised geometry ", + "must refuse, not fall through to one of them.", call. = FALSE) + ) # -------------------------------------------------------- # (3) Membrane damage ODE From 77aa2beb190cf32e1f805aeea87fc1098b797c5b Mon Sep 17 00:00:00 2001 From: aurascoper Date: Fri, 28 Aug 2026 20:54:35 -0500 Subject: [PATCH 04/13] Make both receipt gates prove they reject drift Two of Codex's three on b31faa6. The third (X_red) is a modelling decision and is left open. P2, receipt counters. as.integer() coerces rather than validates, so {"checks_run":"14","failures":false} certified as complete and negative counters passed every `> 0L` comparison. int1() now requires a non-negative whole-number numeric scalar and names the field it refused. Fixtures: valid passes; string, boolean, negative, fractional and array counters each refuse at their own field. P2, numpy case identities. Counting three cases accepted three copies of default_regime, and a case whose checks were deleted has all([]) == True and reports passed -- the can't-fail shape one level down. Moved to scripts/assert_numpy_receipt.py, which pins the case names AND the check names within each case, and carries --self-test running six doctored receipts (case dropped, duplicated, checks deleted, swap control flipped, case renamed, overall pass false). CI runs the self-test BEFORE the real assertion, so a gate that has gone blind fails the build ahead of the thing it is meant to be guarding. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/model-contracts.yml | 39 +++---- scripts/assert_numpy_receipt.py | 143 ++++++++++++++++++++++++++ 2 files changed, 164 insertions(+), 18 deletions(-) create mode 100755 scripts/assert_numpy_receipt.py diff --git a/.github/workflows/model-contracts.yml b/.github/workflows/model-contracts.yml index c21539f..74676e9 100644 --- a/.github/workflows/model-contracts.yml +++ b/.github/workflows/model-contracts.yml @@ -60,22 +60,15 @@ jobs: - name: Verify the biofilm depth stencil (numpy re-coding) run: python analysis/verify_biofilm_depth_profile.py --report py-verify.json + # The gate proves itself first (Codex P2 on #19): counting three cases + # accepted three copies of one, and a case whose checks were deleted has + # `all([]) == True` -- the can't-fail shape one level down. --self-test + # runs it against six doctored receipts and requires each rejected. + - name: Prove the numpy receipt gate rejects drift + run: python scripts/assert_numpy_receipt.py --self-test + - name: Assert the numpy verifier covered every case - run: | - python - <<'PY' - import json, sys - EXPECTED_CASES = 3 - r = json.load(open("py-verify.json")) - names = [c["case"] for c in r["cases"]] - print(f"receipt: {len(names)} cases, passed={r['passed']}: {', '.join(names)}") - if len(names) != EXPECTED_CASES: - sys.exit(f"expected {EXPECTED_CASES} cases, receipt has {len(names)} " - "-- a case was added or silently lost; update EXPECTED_CASES") - bad = [c["case"] for c in r["cases"] if not c["passed"]] - if bad or not r["passed"]: - sys.exit(f"numpy verifier reported failures: {bad}") - print("numpy verifier ran complete.") - PY + run: python scripts/assert_numpy_receipt.py py-verify.json - name: Summarize the time-step diagnostic if: always() @@ -168,9 +161,19 @@ jobs: if (!isTRUE(r$complete)) stop("receipt does not declare itself complete: complete=", format(r$complete)) - ran <- as.integer(r$checks_run); fail <- as.integer(r$failures) - skip <- as.integer(r$skipped) - if (anyNA(c(ran, fail, skip))) stop("receipt counters are not integers") + # VALIDATE, do not coerce (Codex P2 on #19). as.integer() is happy + # to turn "14" into 14 and FALSE into 0, so a broken producer emitting + # {"checks_run":"14","failures":false} certified as complete, and a + # negative counter passed every `> 0L` comparison. + int1 <- function(k) { + v <- r[[k]] + if (!is.numeric(v) || length(v) != 1L || is.na(v) || + v < 0 || v != trunc(v)) + stop("receipt field '", k, "' is not a non-negative integer ", + "scalar: ", paste(format(v), collapse = ", ")) + as.integer(v) + } + ran <- int1("checks_run"); fail <- int1("failures"); skip <- int1("skipped") cat(sprintf("receipt: %d run, %d failed, %d skipped, complete=%s\n", ran, fail, skip, r$complete)) if (skip > 0L) diff --git a/scripts/assert_numpy_receipt.py b/scripts/assert_numpy_receipt.py new file mode 100755 index 0000000..7255ba8 --- /dev/null +++ b/scripts/assert_numpy_receipt.py @@ -0,0 +1,143 @@ +#!/usr/bin/env python3 +"""Assert the numpy depth-stencil receipt describes the run it claims. + +WHY THIS IS NOT A COUNT. The first version checked only that the receipt held +three cases, so three copies of ``default_regime`` passed and silently losing +``strong_gradient`` or the cylindrical negative control stayed green (Codex P2 +on #19). A gate that counts its inputs cannot tell which ones it got. + +So the expected case names and the expected check names WITHIN each case are +both pinned. A case whose internal checks were deleted has an empty ``checks`` +dict, `all([])` is True, and it reports ``passed`` -- which is the same +can't-fail shape one level down. + +``--self-test`` runs the gate against known-bad receipts and requires each to +be rejected, because a gate nobody can test is a gate nobody can trust. +""" + +from __future__ import annotations + +import argparse +import json +import sys +from pathlib import Path + +EXPECTED: dict[str, set[str]] = { + "default_regime": {"accuracy", "second_order", + "negative_control_fails_slab_bound"}, + "strong_gradient": {"accuracy", "second_order", + "negative_control_fails_slab_bound"}, + "cylindrical_refactor_equivalence": {"matches_pre_change_stencil", + "swap_control_differs"}, +} + + +def check(report: dict) -> list[str]: + """Return a list of problems; empty means the receipt is trustworthy.""" + bad: list[str] = [] + cases = report.get("cases") + if not isinstance(cases, list): + return ["receipt has no 'cases' list"] + + names = [c.get("case") for c in cases] + if len(set(names)) != len(names): + bad.append(f"duplicate case names: {names}") + if set(names) != set(EXPECTED): + missing = sorted(set(EXPECTED) - set(names)) + extra = sorted(set(names) - set(EXPECTED)) + bad.append(f"case identities differ -- missing {missing}, unexpected {extra}" + " (update EXPECTED if this is deliberate)") + + for case in cases: + name = case.get("case") + if name not in EXPECTED: + continue + got = set(case.get("checks", {})) + if got != EXPECTED[name]: + bad.append(f"{name}: checks differ -- missing " + f"{sorted(EXPECTED[name] - got)}, unexpected {sorted(got - EXPECTED[name])}") + for cname, ok in case.get("checks", {}).items(): + if not ok: + bad.append(f"{name}: check {cname} failed") + if not case.get("passed"): + bad.append(f"{name}: case reported not passed") + + if not report.get("passed"): + bad.append("report does not declare overall pass") + return bad + + +def self_test() -> int: + """Every doctored receipt below must be REJECTED.""" + good = { + "passed": True, + "cases": [{"case": n, "passed": True, "checks": {c: True for c in ch}} + for n, ch in EXPECTED.items()], + } + if check(good): + print(f"SELF-TEST BROKEN: a valid receipt was rejected: {check(good)}") + return 1 + + def drop_case(r): + r["cases"] = r["cases"][:2]; return r + + def duplicate_case(r): + r["cases"][1] = json.loads(json.dumps(r["cases"][0])); return r + + def gut_checks(r): + r["cases"][0]["checks"] = {}; return r + + def flip_control(r): + r["cases"][2]["checks"]["swap_control_differs"] = False; return r + + def rename_case(r): + r["cases"][1]["case"] = "strong_gradient_v2"; return r + + def overall_false(r): + r["passed"] = False; return r + + failures = 0 + for name, mutate in [("a case dropped", drop_case), + ("a case duplicated", duplicate_case), + ("a case's checks deleted", gut_checks), + ("the swap control flipped", flip_control), + ("a case renamed", rename_case), + ("overall pass false", overall_false)]: + problems = check(mutate(json.loads(json.dumps(good)))) + status = "REJECTED" if problems else "ACCEPTED -- GATE IS BLIND" + print(f" {name:<28} {status}") + if not problems: + failures += 1 + print("self-test: gate rejects every doctored receipt." + if not failures else f"self-test: {failures} doctored receipt(s) slipped through.") + return 1 if failures else 0 + + +def main() -> int: + ap = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + ap.add_argument("receipt", nargs="?", type=Path) + ap.add_argument("--self-test", action="store_true") + args = ap.parse_args() + + if args.self_test: + return self_test() + if args.receipt is None: + return ap.error("a receipt path is required unless --self-test is given") + if not args.receipt.exists(): + print(f"no receipt at {args.receipt}: the verifier did not run to completion") + return 1 + + report = json.loads(args.receipt.read_text()) # invalid JSON raises, which is the point + problems = check(report) + names = [c.get("case") for c in report.get("cases", [])] + print(f"receipt: {len(names)} cases -- {', '.join(map(str, names))}") + for p in problems: + print(f" REJECTED: {p}") + if problems: + return 1 + print("numpy verifier ran complete, with the expected cases and checks.") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) From 6a7d59260a458b1fa0021e003a820f090509d950 Mon Sep 17 00:00:00 2001 From: aurascoper Date: Fri, 28 Aug 2026 21:10:31 -0500 Subject: [PATCH 05/13] Declare the reducing basis a fraction of X_total Codex P1 on #19. X_red was fixed at 0.3 while X_total was caller-declared, so U did not scale with the declared basis, and at a corrected low X_total the reducing mass exceeded total biomass. The fraction reading was always the intent and was never implemented: the constructor comment said "metal-reducing fraction (Shewanella proxy)", the Shiny slider is labelled "fraction" with min 0 max 1, and docs/research/radiotrophic_calibration_map.md:395 already writes the contract as X_red = f_red,dry * X_total. Only the arithmetic disagreed. U = X_total * (k_ads + k_red * f_red_active), f_red_active in [0, 1] in one helper, uptake_rate_of(), used by both radiodialysis_rhs() and penetration_depth(), so the multiplication exists once. Two consequences are now structural rather than asserted: X_red <= X_total by construction, and U is proportional to X_total, so the whole uptake term inherits that gate instead of needing a second one. Check 15 asserts all three halves -- U proportional to X_total, lambda ratio sqrt(10) between X_total 1.0 and 0.1, and a fraction outside [0,1] refused (1.5, -0.1, NA, a vector, a string). Confirmed to bite by restoring the additive form, which fails it at U(0.1)/U(1)=0.196 and lambda ratio 2.256 against sqrt(10)=3.162. Nothing shipping moves: at X_total = 1.0 the arithmetic is identical, and the cylindrical byte-equivalence regression still reads 2.28e-16. ON PROVENANCE, since it decides what 0.3 means rather than only where it multiplies. It is NOT 2/7 of the seven species: the coupled path counts one, S. oneidensis (biofilms_potts.jl:1377), and 0.3 traces to 45de4ba with no recorded basis. It is labelled a TAXONOMIC proxy standing in an ACTIVITY slot, which active_from_taxonomic() refuses without a measured activity fraction. So 0.3 stays an unvalidated placeholder gated by RADIODIALYSIS: BLOCKED, and the comment says so. Suites: coupling 310 passed 6 skipped, calibration 360, contract 7, R verifier 15/15. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/model-contracts.yml | 2 +- analysis/verify_biofilm_depth_profile.R | 37 ++++++++++++-- analysis/verify_biofilm_depth_profile.py | 11 +++-- biofilms_radiodialysis.R | 61 ++++++++++++++++++++---- 4 files changed, 93 insertions(+), 18 deletions(-) diff --git a/.github/workflows/model-contracts.yml b/.github/workflows/model-contracts.yml index 74676e9..5e4e10b 100644 --- a/.github/workflows/model-contracts.yml +++ b/.github/workflows/model-contracts.yml @@ -152,7 +152,7 @@ jobs: - name: Assert the verifier actually ran everything shell: Rscript {0} run: | - EXPECTED_CHECKS <- 14L + EXPECTED_CHECKS <- 15L if (!file.exists("r-verify.json")) stop("no receipt: the verifier did not run to completion") r <- jsonlite::fromJSON("r-verify.json", simplifyVector = TRUE) diff --git a/analysis/verify_biofilm_depth_profile.R b/analysis/verify_biofilm_depth_profile.R index c0418a2..8327408 100644 --- a/analysis/verify_biofilm_depth_profile.R +++ b/analysis/verify_biofilm_depth_profile.R @@ -42,7 +42,8 @@ MODEL <- file.path(dirname(dirname(normalizePath( if (!file.exists(MODEL)) MODEL <- "biofilms_radiodialysis.R" WANTED <- c("radiodialysis_rhs", "face_weights", "default_parms", - "slab_parms", "penetration_depth", "run_radiodialysis") + "slab_parms", "penetration_depth", "run_radiodialysis", + "uptake_rate_of") for (e in parse(MODEL)) { if (is.call(e) && identical(as.character(e[[1]]), "<-") && as.character(e[[2]]) %in% WANTED) eval(e, envir = globalenv()) @@ -118,7 +119,7 @@ operator_of <- function(p) { # and every mutant looked "caught" for the wrong reason. reference <- function(p, geom) { g <- p$r_grid; N <- length(g); dg <- g[2] - g[1]; D <- p$D_eff - U <- p$k_ads * p$X_total + p$k_red * p$X_red + U <- uptake_rate_of(p$X_total, p$f_red_active, p$k_ads, p$k_red) bc <- if (identical(geom, "slab")) p$k_L else p$P0 wp <- if (identical(geom, "slab")) rep(1, N) else (g + 0.5 * dg) / g wm <- if (identical(geom, "slab")) rep(1, N) else (g - 0.5 * dg) / g @@ -166,11 +167,39 @@ local({ refused, accepts[1], accepts[2])) }) +# --- 1c. the reducing basis is a FRACTION of X_total ----------------------- +# Codex P1 on #19. X_red was fixed at 0.3 while X_total was caller-declared, +# so U did not scale with the declared basis and a corrected low X_total put +# the reducing mass ABOVE total biomass. Three halves, all required: +# (a) U is proportional to X_total, so lambda ~ 1/sqrt(X_total); +# (b) the implied X_red = f*X_total never exceeds X_total, structurally; +# (c) a fraction outside [0,1] is refused rather than used. +local({ + u1 <- uptake_rate_of(1.0, 0.3, 0.05, 0.02) + u01 <- uptake_rate_of(0.1, 0.3, 0.05, 0.02) + scales <- abs(u01 / u1 - 0.1) < 1e-12 + lam_ratio <- penetration_depth(slab_parms(X_total = 0.1)) / + penetration_depth(slab_parms(X_total = 1.0)) + lam_ok <- abs(lam_ratio - sqrt(10)) < 1e-9 + bounded <- all(vapply(c(1e-6, 0.1, 1.0, 10.0), function(xt) { + p <- slab_parms(X_total = xt); p$f_red_active * p$X_total <= p$X_total + }, logical(1))) + refuses <- all(vapply(list(1.5, -0.1, NA_real_, c(0.3, 0.3), "0.3"), + function(bad) tryCatch({ uptake_rate_of(1.0, bad, 0.05, 0.02); FALSE }, + error = function(e) grepl("fraction in [0, 1]", + conditionMessage(e), + fixed = TRUE)), logical(1))) + report(scales && lam_ok && bounded && refuses, + "U scales with X_total, X_red <= X_total by construction, f outside [0,1] refused", + sprintf("U(0.1)/U(1)=%.3f lambda ratio=%.6f (sqrt10=%.6f) bounded=%s refuses=%s", + u01/u1, lam_ratio, sqrt(10), bounded, refuses)) +}) + # --- 2-3. slab operator and steady profile --------------------------------- L_f <- 0.01; D <- 1e-5; N <- 40; kL <- 1e2 k_eff <- 0.056 * 0.001 / 0.006 p <- slab_parms(Nr = N, L_f = L_f, X_total = 1.0, D_eff = D, k_L = kL) -p$k_ads <- k_eff; p$k_red <- 0; p$k_des <- 0; p$k_loss <- 0; p$X_red <- 0 +p$k_ads <- k_eff; p$k_red <- 0; p$k_des <- 0; p$k_loss <- 0; p$f_red_active <- 0 res <- residual(p, "slab") report(res < 1e-13, "slab operator matches independent assembly", sprintf("%.2e", res)) @@ -254,7 +283,7 @@ if (!requireNamespace("deSolve", quietly = TRUE)) { out <- run_radiodialysis(ps, t_end = 2e5, n_out = 50) final <- out$c_mat[nrow(out$c_mat), ] # full system: sorption chain on, so the steady sink is k_eff - ke <- (ps$k_ads * ps$X_total + ps$k_red * ps$X_red) * + ke <- uptake_rate_of(ps$X_total, ps$f_red_active, ps$k_ads, ps$k_red) * ps$k_loss / (ps$k_des + ps$k_loss) lm2 <- sqrt(ps$D_eff / ke) ex <- ps$c_ext * cosh(ps$r_grid / lm2) / diff --git a/analysis/verify_biofilm_depth_profile.py b/analysis/verify_biofilm_depth_profile.py index 0ab3820..cc5ac71 100644 --- a/analysis/verify_biofilm_depth_profile.py +++ b/analysis/verify_biofilm_depth_profile.py @@ -16,7 +16,7 @@ -------------------------- Setting ds/dt = 0 in the immobile-phase equation gives - s = U c / (k_des + k_loss), U = k_ads*X_total + k_red*X_red + s = U c / (k_des + k_loss), U = X_total * (k_ads + k_red*f_red_active) and substituting it back shows the sorption/desorption pair is a *reversible buffer* that cancels. The only true steady sink is precipitation: @@ -106,13 +106,18 @@ class Kinetics: k_red: float = 0.02 k_des: float = 0.005 k_loss: float = 0.001 - X_red: float = 0.3 + f_red_active: float = 0.3 # ACTIVE-reducer fraction OF X_total, in [0,1] X_total: float = 1.0 # TEST VALUE -- gated, see module docstring def validate(self) -> None: for name, value in asdict(self).items(): if not math.isfinite(value) or value < 0.0: raise ValueError(f"{name} must be finite and non-negative, got {value!r}") + if self.f_red_active > 1.0: + raise ValueError( + "f_red_active is a fraction OF X_total and must lie in [0, 1], " + f"got {self.f_red_active!r}. It is not a second biomass density; " + "the previous fixed X_red let the reducing mass exceed X_total.") if self.k_des + self.k_loss <= 0.0: raise ValueError("k_des + k_loss must be positive for a steady state to exist") @@ -120,7 +125,7 @@ def validate(self) -> None: def uptake(self) -> float: """U, the transient sink (s^-1).""" self.validate() - return self.k_ads * self.X_total + self.k_red * self.X_red + return self.X_total * (self.k_ads + self.k_red * self.f_red_active) @property def k_eff(self) -> float: diff --git a/biofilms_radiodialysis.R b/biofilms_radiodialysis.R index 53ef253..61e2977 100644 --- a/biofilms_radiodialysis.R +++ b/biofilms_radiodialysis.R @@ -7,10 +7,10 @@ # # (1) Mobile species: cylindrical reaction-diffusion PDE # ∂c/∂t = (1/r) ∂/∂r (r D_eff ∂c/∂r) -# - (k_ads X + k_red X_red) c + k_des s +# - X (k_ads + k_red f_red) c + k_des s # # (2) Immobile phase ODE (biosorption / bioreduction): -# ∂s/∂t = (k_ads X + k_red X_red) c - (k_des + k_loss) s +# ∂s/∂t = X (k_ads + k_red f_red) c - (k_des + k_loss) s # # (3) Membrane damage ODE (radiation-driven permeability): # dm/dt = -k_dam Ḋ(R) m @@ -95,7 +95,7 @@ radiodialysis_rhs <- function(t, y, parms) { # -------------------------------------------------------- # Source/sink term (identical for mobile and immobile) # -------------------------------------------------------- - uptake_rate <- k_ads * X_total + k_red * X_red # s⁻¹ + uptake_rate <- uptake_rate_of(X_total, f_red_active, k_ads, k_red) # -------------------------------------------------------- # (1) Mobile species — diffusion + reaction (geometry via w_plus/w_minus) @@ -155,6 +155,45 @@ radiodialysis_rhs <- function(t, y, parms) { # remediation context. # ------------------------------------------------------------ +# ------------------------------------------------------------ +# Uptake rate U (s⁻¹). +# +# f_red_active is a FRACTION of X_total, not a second density: +# +# X_red = f_red_active * X_total, f_red_active in [0, 1] +# +# which is the contract already written in +# docs/research/radiotrophic_calibration_map.md:395 as +# X_red = f_red,dry * X_total. Two consequences, both structural: +# +# * X_red <= X_total holds by construction rather than by an assertion +# nobody wrote. The previous form fixed X_red at 0.3 while X_total was +# caller-declared, so a corrected low-biomass basis produced a reducing +# mass EXCEEDING total biomass, and U did not scale with the declared +# basis at all (Codex P1 on #19). +# * U = X_total * (k_ads + k_red * f_red_active) is proportional to +# X_total, so the whole uptake term inherits the X_total gate. There is +# no second gated quantity to forget. +# +# WHAT f_red_active IS NOT. It is an ACTIVE-reducer dry-mass fraction, not a +# taxonomic one. The default 0.3 is inherited, labelled "(Shewanella proxy)" +# since 45de4ba, and is a TAXONOMIC proxy standing in an activity slot -- the +# very substitution active_from_taxonomic() refuses without a measured +# activity fraction ("taxonomic abundance is not functional activity"). It is +# not 2/7 of seven species: the coupled path counts one species, S. oneidensis +# (biofilms_potts.jl:1377). Treat 0.3 as an unvalidated placeholder gated by +# RADIODIALYSIS: BLOCKED, not as a composition. +# ------------------------------------------------------------ +uptake_rate_of <- function(X_total, f_red_active, k_ads, k_red) { + if (!is.numeric(f_red_active) || length(f_red_active) != 1L || + is.na(f_red_active) || f_red_active < 0 || f_red_active > 1) + stop("f_red_active must be a single fraction in [0, 1], got ", + paste(format(f_red_active), collapse = ", "), + ". It is a fraction OF X_total, not a second biomass density.", + call. = FALSE) + X_total * (k_ads + k_red * f_red_active) +} + # ------------------------------------------------------------ # Diffusion face weights. The ONLY thing geometry changes. # @@ -202,7 +241,8 @@ default_parms <- function(Nr = 40, R = 1.0) { # --- Biomass --- X_total = 1.0, # total biofilm dry mass density (g cm⁻³) - X_red = 0.3, # metal-reducing fraction (Shewanella proxy) + f_red_active = 0.3, # ACTIVE-reducer fraction OF X_total, in [0,1]. + # Unvalidated placeholder; see uptake_rate_of() above. # --- Membrane (Nafion / Donnan) --- P0 = 0.01, # baseline permeability (cm s⁻¹), Fox et al. 2009 @@ -242,13 +282,13 @@ default_parms <- function(Nr = 40, R = 1.0) { # X_total has NO default on purpose. See the stop() below. # ------------------------------------------------------------ slab_parms <- function(Nr = 40, L_f = 0.01, X_total, - D_eff = 1e-5, k_L = 1e-3) { + D_eff = 1e-5, k_L = 1e-3, f_red_active = 0.3) { if (missing(X_total)) { stop( "slab_parms(): X_total must be stated explicitly, not defaulted.\n", " The biofilm dry-mass basis is exactly what RADIODIALYSIS: BLOCKED\n", " names (README.md:350-353). It sets the uptake rate\n", - " U = k_ads*X_total + k_red*X_red,\n", + " U = X_total * (k_ads + k_red*f_red_active),\n", " and the penetration depth lambda = sqrt(D_eff/k_eff) scales as\n", " 1/sqrt(U), so whether ANY gradient is resolvable across L_f is\n", " downstream of that gate. Pass a value and label it a test value.", @@ -276,7 +316,8 @@ slab_parms <- function(Nr = 40, L_f = 0.01, X_total, # --- Biomass --- X_total = X_total, # caller-declared; gated, see stop() above - X_red = 0.3, + f_red_active = f_red_active, # fraction OF X_total; U scales with the + # declared basis, so it inherits that same gate. # --- Outer face: liquid-film mass transfer, NOT membrane permeability --- k_L = k_L, # (cm s⁻¹). P0 / alpha_P are deliberately absent. @@ -306,7 +347,7 @@ slab_parms <- function(Nr = 40, L_f = 0.01, X_total, penetration_depth <- function(parms, phase = c("steady", "transient")) { phase <- match.arg(phase) with(parms, { - U <- k_ads * X_total + k_red * X_red + U <- uptake_rate_of(X_total, f_red_active, k_ads, k_red) k <- if (phase == "steady") U * k_loss / (k_des + k_loss) else U sqrt(D_eff / k) }) @@ -359,7 +400,7 @@ ui <- fluidPage( min = 0, max = 0.1, value = 0.02, step = 0.002), sliderInput("k_des", "k_des (s⁻¹)", min = 0, max = 0.05, value = 0.005, step = 0.001), - sliderInput("X_red", "Metal-reducing biomass fraction", + sliderInput("f_red_active", "Active metal-reducing fraction of X_total", min = 0, max = 1, value = 0.3, step = 0.05), h4("Membrane"), @@ -414,7 +455,7 @@ server <- function(input, output, session) { parms$k_ads <- input$k_ads parms$k_red <- input$k_red parms$k_des <- input$k_des - parms$X_red <- input$X_red + parms$f_red_active <- input$f_red_active parms$P0 <- input$P0 parms$alpha_P <- input$alpha_P parms$k_dam <- input$k_dam From e43817fbb080da9f6b17a25ef259e03688d0c3ba Mon Sep 17 00:00:00 2001 From: aurascoper Date: Fri, 28 Aug 2026 21:18:08 -0500 Subject: [PATCH 06/13] Guard both operands, both verdict types, and the gate's own path Three of Codex's four on 6a7d592. The fourth (the Julia coupled solvers) is left open; it is entangled with the gate. X_total was unvalidated. The structural bound f*X_total <= X_total assumes a finite non-negative scalar, and a hand-built parms list with X_total = -1 reached both callers, returned a negative uptake, and broke the bound the helper is cited for. A guard checking one of two operands proves nothing about their product. Both are validated now, and check 15 exercises malformed totals (-1, NA, Inf, -Inf, a vector, a string) alongside the fraction ones. String verdicts were certified as passing. JSON "false" is a non-empty string and therefore truthy, so `not ok` accepted a producer schema regression whose every verdict was the string "false". Verdicts must now be boolean true (`is not True`), and two doctored receipts covering it joined the self-test, which now runs eight. The gate's own script was outside the workflow's path filters, so a pull request changing only scripts/assert_numpy_receipt.py ran none of its negative controls. scripts/** added to both filters -- a gate that does not run when it changes is the same defect one level out. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/model-contracts.yml | 2 ++ analysis/verify_biofilm_depth_profile.R | 22 +++++++++++-------- biofilms_radiodialysis.R | 20 ++++++++++++----- scripts/assert_numpy_receipt.py | 29 +++++++++++++++++++------ 4 files changed, 51 insertions(+), 22 deletions(-) diff --git a/.github/workflows/model-contracts.yml b/.github/workflows/model-contracts.yml index 5e4e10b..8ac85f6 100644 --- a/.github/workflows/model-contracts.yml +++ b/.github/workflows/model-contracts.yml @@ -5,6 +5,7 @@ on: branches: [master, "feat/**", "ci/**", "research/**"] paths: - "analysis/**" + - "scripts/**" - "*.R" - "*.jl" - "preprint/**" @@ -12,6 +13,7 @@ on: pull_request: paths: - "analysis/**" + - "scripts/**" - "*.R" - "*.jl" - "preprint/**" diff --git a/analysis/verify_biofilm_depth_profile.R b/analysis/verify_biofilm_depth_profile.R index 8327408..b87b810 100644 --- a/analysis/verify_biofilm_depth_profile.R +++ b/analysis/verify_biofilm_depth_profile.R @@ -184,15 +184,19 @@ local({ bounded <- all(vapply(c(1e-6, 0.1, 1.0, 10.0), function(xt) { p <- slab_parms(X_total = xt); p$f_red_active * p$X_total <= p$X_total }, logical(1))) - refuses <- all(vapply(list(1.5, -0.1, NA_real_, c(0.3, 0.3), "0.3"), - function(bad) tryCatch({ uptake_rate_of(1.0, bad, 0.05, 0.02); FALSE }, - error = function(e) grepl("fraction in [0, 1]", - conditionMessage(e), - fixed = TRUE)), logical(1))) - report(scales && lam_ok && bounded && refuses, - "U scales with X_total, X_red <= X_total by construction, f outside [0,1] refused", - sprintf("U(0.1)/U(1)=%.3f lambda ratio=%.6f (sqrt10=%.6f) bounded=%s refuses=%s", - u01/u1, lam_ratio, sqrt(10), bounded, refuses)) + malformed <- function(f) tryCatch({ f(); FALSE }, + error = function(e) grepl("must be a single finite number", + conditionMessage(e), fixed = TRUE)) + # BOTH operands. A guard checking only the fraction proves nothing about + # the product, and X_total = -1 broke the bound this helper is cited for. + refuses_f <- all(vapply(list(1.5, -0.1, NA_real_, Inf, c(0.3, 0.3), "0.3"), + function(bad) malformed(function() uptake_rate_of(1.0, bad, 0.05, 0.02)), logical(1))) + refuses_x <- all(vapply(list(-1, NA_real_, Inf, -Inf, c(1, 1), "1"), + function(bad) malformed(function() uptake_rate_of(bad, 0.3, 0.05, 0.02)), logical(1))) + report(scales && lam_ok && bounded && refuses_f && refuses_x, + "U scales with X_total, X_red <= X_total by construction, malformed X_total and f refused", + sprintf("U(0.1)/U(1)=%.3f lambda ratio=%.6f (sqrt10=%.6f) bounded=%s refuses_f=%s refuses_X=%s", + u01/u1, lam_ratio, sqrt(10), bounded, refuses_f, refuses_x)) }) # --- 2-3. slab operator and steady profile --------------------------------- diff --git a/biofilms_radiodialysis.R b/biofilms_radiodialysis.R index 61e2977..e711568 100644 --- a/biofilms_radiodialysis.R +++ b/biofilms_radiodialysis.R @@ -185,12 +185,20 @@ radiodialysis_rhs <- function(t, y, parms) { # RADIODIALYSIS: BLOCKED, not as a composition. # ------------------------------------------------------------ uptake_rate_of <- function(X_total, f_red_active, k_ads, k_red) { - if (!is.numeric(f_red_active) || length(f_red_active) != 1L || - is.na(f_red_active) || f_red_active < 0 || f_red_active > 1) - stop("f_red_active must be a single fraction in [0, 1], got ", - paste(format(f_red_active), collapse = ", "), - ". It is a fraction OF X_total, not a second biomass density.", - call. = FALSE) + # X_total is validated HERE too (Codex on #19): the structural bound + # f*X_total <= X_total assumes a finite non-negative scalar, and a hand-built + # parms list carrying X_total = -1 reached both callers, returned a negative + # uptake, and broke the very bound this helper is cited for. A guard that + # checks one of its two operands proves nothing about their product. + finite_scalar <- function(x, nm, hi = Inf) { + if (!is.numeric(x) || length(x) != 1L || !is.finite(x) || x < 0 || x > hi) + stop(nm, " must be a single finite number in [0, ", + if (is.finite(hi)) hi else "Inf", "], got ", + paste(format(x), collapse = ", "), call. = FALSE) + } + finite_scalar(X_total, "X_total") + finite_scalar(f_red_active, + "f_red_active (a fraction OF X_total, not a second density)", 1) X_total * (k_ads + k_red * f_red_active) } diff --git a/scripts/assert_numpy_receipt.py b/scripts/assert_numpy_receipt.py index 7255ba8..77fc70d 100755 --- a/scripts/assert_numpy_receipt.py +++ b/scripts/assert_numpy_receipt.py @@ -56,14 +56,18 @@ def check(report: dict) -> list[str]: if got != EXPECTED[name]: bad.append(f"{name}: checks differ -- missing " f"{sorted(EXPECTED[name] - got)}, unexpected {sorted(got - EXPECTED[name])}") + # `is not True`, not `not ok`: JSON "false" is a non-empty string and + # therefore truthy, so a producer schema regression emitting string + # verdicts was certified as passing (Codex on #19). A verdict must be + # an actual boolean true. for cname, ok in case.get("checks", {}).items(): - if not ok: - bad.append(f"{name}: check {cname} failed") - if not case.get("passed"): - bad.append(f"{name}: case reported not passed") + if ok is not True: + bad.append(f"{name}: check {cname} is {ok!r}, not boolean true") + if case.get("passed") is not True: + bad.append(f"{name}: case verdict is {case.get('passed')!r}, not boolean true") - if not report.get("passed"): - bad.append("report does not declare overall pass") + if report.get("passed") is not True: + bad.append(f"overall verdict is {report.get('passed')!r}, not boolean true") return bad @@ -96,13 +100,24 @@ def rename_case(r): def overall_false(r): r["passed"] = False; return r + def string_verdicts(r): + """The truthiness trap: JSON "false" is a non-empty string.""" + r["cases"][0]["passed"] = "false" + r["cases"][0]["checks"] = {c: "false" for c in r["cases"][0]["checks"]} + return r + + def overall_string_true(r): + r["passed"] = "true"; return r + failures = 0 for name, mutate in [("a case dropped", drop_case), ("a case duplicated", duplicate_case), ("a case's checks deleted", gut_checks), ("the swap control flipped", flip_control), ("a case renamed", rename_case), - ("overall pass false", overall_false)]: + ("overall pass false", overall_false), + ("string verdicts \"false\"", string_verdicts), + ("overall verdict a string", overall_string_true)]: problems = check(mutate(json.loads(json.dumps(good)))) status = "REJECTED" if problems else "ACCEPTED -- GATE IS BLIND" print(f" {name:<28} {status}") From 301951abd9904b31971f44a24a02122339f15276 Mon Sep 17 00:00:00 2001 From: aurascoper Date: Fri, 28 Aug 2026 21:43:14 -0500 Subject: [PATCH 07/13] Enforce the radiodialysis basis gate instead of relying on an accident Codex's P1 on 6a7d592, resolved by refusing rather than by conforming. The Julia coupled path was blocked only by a side effect: the pre-fraction additive uptake happens to misbehave at non-unit X_total, and that was doing the work of a guard. A tripwire any tidy-up can remove is not a gate. So step_radiolysis! -- the entry point both ports share -- now refuses outright, naming the occupancy-to-fraction defect. The arithmetic is deliberately NOT made conformant. Mirroring R's uptake_rate_of() here would make this path conformant without making it correct: the X_red reaching it is red_cells[i]/counts[i], one species' occupied sites over all interior sites, which README.md:344 records as neither a biomass fraction nor a reducer fraction. That would turn a visible defect into an invisible one. The additive form stays as a marker that the reconciliation is unfinished. WHAT THIS TURNED UP. validate_serial.jl steps the coupled path at X_total = 0.065 -- a site-occupancy mean, exactly the gated quantity -- so enforcing the block broke the CSV determinism fixture. Its recorded columns do not depend on the gated basis: the CSV is CPM quantities plus rd.m, whose ODE (dm/dt = -k_dam*Ddot_R*m) contains neither X_total nor X_red. So the gate takes one narrow acknowledgement, basis_gate_ack, for runs that compare blocked output against blocked output and assert no magnitude. Three sites hold it, not one: validate_serial.jl, and the genealogy and checkpoint tests, which read rd.c/rd.s only to compare two code paths or a round trip. Two guards keep that honest. A census test enumerates the sites, so a fourth fails the build rather than appearing quietly. And an independence test runs the harness with the uptake constants at 10x and requires every recorded column byte-identical -- the day anything c- or s-derived is recorded, the exemption's claim fails with it. The gate test asserts the asymmetry rather than agreement: R scales with the basis (check 15), this path refuses. Confirmed to bite by deleting the guard call, which fails it six times. Bool not Symbol: the field is HDF5-serialised by export_checkpoint.jl and has to survive a restart. Julia 2/34/30/68/42/9/14/22, R 15/15, numpy 3/3. Co-Authored-By: Claude Opus 5 (1M context) --- biofilms_potts.jl | 71 ++++++++++++++++ biofilms_potts_jacc.jl | 81 +++++++++++++++++- tests/checkpoint_io_tests.jl | 7 +- tests/genealogy_tests.jl | 7 +- tests/radiodialysis_basis_gate.jl | 134 ++++++++++++++++++++++++++++++ tests/runtests.jl | 2 + validate_serial.jl | 10 ++- 7 files changed, 306 insertions(+), 6 deletions(-) create mode 100644 tests/radiodialysis_basis_gate.jl diff --git a/biofilms_potts.jl b/biofilms_potts.jl index c9ace1f..4efb439 100644 --- a/biofilms_potts.jl +++ b/biofilms_potts.jl @@ -1193,6 +1193,13 @@ Base.@kwdef struct RadiolysisParams X_total::Float64 = 1.0 # total dry-mass density (g cm⁻³) X_red::Float64 = 0.3 # metal-reducing fraction (Shewanella proxy) + # Basis gate acknowledgement. false refuses any X_total != 1.0. + # Bool rather than a Symbol because this field is HDF5-serialised by + # export_checkpoint.jl and must survive a restart; the exemption is binary + # anyway, and WHICH sites hold it is pinned by the census test rather than + # by a symbol name. See _assert_basis_gate. + basis_gate_ack::Bool = false + # Membrane (Nafion / Donnan — Fox et al. 2009, Lara et al. 2023) P0::Float64 = 0.01 # baseline permeability (cm s⁻¹) alpha_P::Float64 = 0.02 # radiation-damage coefficient (Gy⁻¹) @@ -1252,6 +1259,7 @@ R = 1.0 cm makes dt_rd = 0.5 genuinely unstable, and this guard is what absorbs it. `biofilms_potts_jacc.jl` carries the identical wrapper for the same reason. """ function step_radiolysis!(rd::RadiolysisState, dt::Float64) + _assert_basis_gate(rd.params.X_total, rd.params.basis_gate_ack) dr = rd.r_grid[2] - rd.r_grid[1] dt_stable = 0.4 * dr^2 / (2.0 * rd.params.D_eff) n_sub = max(1, ceil(Int, dt / dt_stable)) @@ -1261,6 +1269,67 @@ function step_radiolysis!(rd::RadiolysisState, dt::Float64) end end +""" + _assert_basis_gate(X_total) + +Refuse to integrate on a coupled biomass basis, EXPLICITLY. + +RADIODIALYSIS: BLOCKED gates the biomass basis fed into this coupling, and +until now nothing here enforced it. The block was being done by an accident: +`uptake = k_ads*X_total + k_red*X_red` is the pre-fraction additive form, which +happens to misbehave at non-unit `X_total`, and that side effect was standing in +for a guard. An accidental tripwire is not a gate -- it can be removed by +someone tidying the arithmetic, leaving nothing to say the basis was blocked. + +So the refusal is stated. `X_total == 1.0` is the standalone default and stays +allowed; anything else means a real basis was supplied, which is what the gate +covers. + +WHY THE ARITHMETIC IS NOT ALSO FIXED. `biofilms_radiodialysis.R` now derives +`X_red = f_red_active * X_total` (`uptake_rate_of()`), and mirroring that here +would make this path *conformant*. It would not make it *correct*: the `X_red` +reaching it is `red_cells[i] / counts[i]` (`compute_radial_biomass`), one +species' occupied sites over ALL interior sites, which README.md:344 records as +"neither a biomass fraction nor a reducer fraction". Conformant arithmetic over +a quantity the repository has already refused is worse than visibly +non-conformant arithmetic over the same one: the defect would stop being visible +while staying just as gated. The arithmetic stays as a marker that this +reconciliation is unfinished. + +Nor can the fraction be derived from parcel counts. Counts give a TAXONOMIC +fraction, and `active_from_taxonomic()` refuses converting one to an +active-reducer fraction without a measured activity fraction; `D-XRED` in +`data/calibration/reference_d_requirements.csv` records that as blocked by this +units error rather than by missing data. +""" +function _assert_basis_gate(X_total::Float64, ack::Bool = false) + X_total == 1.0 && return nothing + # THE ONE EXEMPTION. validate_serial.jl steps this path only to reproduce a + # bit-for-bit CPM trajectory, and records no radiodialysis quantity: its CSV + # carries CPM columns plus rd.m, whose ODE (dm/dt = -k_dam*Ddot_R*m) has no + # X_total or X_red in it. That independence is not taken on trust -- it is + # asserted in tests/radiodialysis_basis_gate.jl by running the harness at + # two different gated bases and requiring byte-identical CSV. Widen this + # exemption and that test is what should stop you. + ack && return nothing + error(""" + RADIODIALYSIS: BLOCKED -- refusing to integrate at X_total = $X_total. + + Only the standalone default X_total == 1.0 is allowed here. + + A coupled basis reaches this path as mean(compute_radial_biomass(...)): + one species' occupied sites over all interior sites, which is + neither a biomass fraction nor a reducer fraction + (README.md:344). The uptake arithmetic here is also still the + pre-fraction additive form, unlike biofilms_radiodialysis.R's + uptake_rate_of(). + + This refusal is deliberate and is asserted by + tests/radiodialysis_basis_gate.jl. Removing it to let a coupled run + proceed re-opens the defect the gate names; repair the quantity first. + """) +end + """ One forward-Euler step of equations (1)–(3) with ghost-point Robin BC at r = R. """ @@ -1490,6 +1559,7 @@ function run_simulation_coupled(params::CPMParams, rp::RadiolysisParams, k_loss = rp.k_loss, X_total = mean(X_tot), X_red = mean(X_rd), + basis_gate_ack = rp.basis_gate_ack, P0 = rp.P0, alpha_P = rp.alpha_P, k_dam = rp.k_dam, @@ -1854,6 +1924,7 @@ function advance_window!(sim::CoupledSimulation, n_mcs::Int) Nr = rp.Nr, D_eff = rp.D_eff, k_ads = rp.k_ads, k_red = rp.k_red, k_des = rp.k_des, k_loss = rp.k_loss, X_total = mean(X_tot), X_red = mean(X_rd), + basis_gate_ack = rp.basis_gate_ack, P0 = rp.P0, alpha_P = rp.alpha_P, k_dam = rp.k_dam, Ddot_R = rp.Ddot_R, c_ext = rp.c_ext, dt_rd = rp.dt_rd) end diff --git a/biofilms_potts_jacc.jl b/biofilms_potts_jacc.jl index a92da7d..04d9c3f 100644 --- a/biofilms_potts_jacc.jl +++ b/biofilms_potts_jacc.jl @@ -260,6 +260,13 @@ Base.@kwdef struct RadiolysisParams k_loss::Float64 = 0.001 X_total::Float64 = 1.0 X_red::Float64 = 0.3 + + # Basis gate acknowledgement. false refuses any X_total != 1.0. + # Bool rather than a Symbol because this field is HDF5-serialised by + # export_checkpoint.jl and must survive a restart; the exemption is binary + # anyway, and WHICH sites hold it is pinned by the census test rather than + # by a symbol name. See _assert_basis_gate. + basis_gate_ack::Bool = false P0::Float64 = 0.01 alpha_P::Float64 = 0.02 k_dam::Float64 = 0.005 @@ -299,6 +306,7 @@ small. Correct the units to R = 1.0 cm and dt_stable becomes 0.132, n_sub = 4 at which point a port without this guard diverges. """ function step_radiolysis!(rd::RadiolysisState, dt::Float64) + _assert_basis_gate(rd.params.X_total, rd.params.basis_gate_ack) dr = rd.r_grid[2] - rd.r_grid[1] dt_stable = 0.4 * dr^2 / (2.0 * rd.params.D_eff) n_sub = max(1, ceil(Int, dt / dt_stable)) @@ -308,6 +316,67 @@ function step_radiolysis!(rd::RadiolysisState, dt::Float64) end end +""" + _assert_basis_gate(X_total) + +Refuse to integrate on a coupled biomass basis, EXPLICITLY. + +RADIODIALYSIS: BLOCKED gates the biomass basis fed into this coupling, and +until now nothing here enforced it. The block was being done by an accident: +`uptake = k_ads*X_total + k_red*X_red` is the pre-fraction additive form, which +happens to misbehave at non-unit `X_total`, and that side effect was standing in +for a guard. An accidental tripwire is not a gate -- it can be removed by +someone tidying the arithmetic, leaving nothing to say the basis was blocked. + +So the refusal is stated. `X_total == 1.0` is the standalone default and stays +allowed; anything else means a real basis was supplied, which is what the gate +covers. + +WHY THE ARITHMETIC IS NOT ALSO FIXED. `biofilms_radiodialysis.R` now derives +`X_red = f_red_active * X_total` (`uptake_rate_of()`), and mirroring that here +would make this path *conformant*. It would not make it *correct*: the `X_red` +reaching it is `red_cells[i] / counts[i]` (`compute_radial_biomass`), one +species' occupied sites over ALL interior sites, which README.md:344 records as +"neither a biomass fraction nor a reducer fraction". Conformant arithmetic over +a quantity the repository has already refused is worse than visibly +non-conformant arithmetic over the same one: the defect would stop being visible +while staying just as gated. The arithmetic stays as a marker that this +reconciliation is unfinished. + +Nor can the fraction be derived from parcel counts. Counts give a TAXONOMIC +fraction, and `active_from_taxonomic()` refuses converting one to an +active-reducer fraction without a measured activity fraction; `D-XRED` in +`data/calibration/reference_d_requirements.csv` records that as blocked by this +units error rather than by missing data. +""" +function _assert_basis_gate(X_total::Float64, ack::Bool = false) + X_total == 1.0 && return nothing + # THE ONE EXEMPTION. validate_serial.jl steps this path only to reproduce a + # bit-for-bit CPM trajectory, and records no radiodialysis quantity: its CSV + # carries CPM columns plus rd.m, whose ODE (dm/dt = -k_dam*Ddot_R*m) has no + # X_total or X_red in it. That independence is not taken on trust -- it is + # asserted in tests/radiodialysis_basis_gate.jl by running the harness at + # two different gated bases and requiring byte-identical CSV. Widen this + # exemption and that test is what should stop you. + ack && return nothing + error(""" + RADIODIALYSIS: BLOCKED -- refusing to integrate at X_total = $X_total. + + Only the standalone default X_total == 1.0 is allowed here. + + A coupled basis reaches this path as mean(compute_radial_biomass(...)): + one species' occupied sites over all interior sites, which is + neither a biomass fraction nor a reducer fraction + (README.md:344). The uptake arithmetic here is also still the + pre-fraction additive form, unlike biofilms_radiodialysis.R's + uptake_rate_of(). + + This refusal is deliberate and is asserted by + tests/radiodialysis_basis_gate.jl. Removing it to let a coupled run + proceed re-opens the defect the gate names; repair the quantity first. + """) +end + """One forward-Euler step; the guard lives in step_radiolysis! above.""" function _step_radiolysis_euler!(rd::RadiolysisState, dt::Float64) rp = rd.params @@ -426,9 +495,15 @@ function run_coupled(; N::Int = 40, n_cells_per_species::Int = 6, if mcs % 10 == 1 X_tot, X_red = radial_biomass(JACC.to_host(lat), spec_h, N, rd.params.Nr) rp0 = rd.params - rd.params = RadiolysisParams(rp0.Nr, rp0.D_eff, rp0.k_ads, rp0.k_red, - rp0.k_des, rp0.k_loss, mean(X_tot), mean(X_red), rp0.P0, - rp0.alpha_P, rp0.k_dam, rp0.Ddot_R, rp0.c_ext, rp0.dt_rd) + # Keyword form deliberately: this was positional, so adding a + # field silently shifted every argument past it by one. + rd.params = RadiolysisParams( + Nr = rp0.Nr, D_eff = rp0.D_eff, k_ads = rp0.k_ads, + k_red = rp0.k_red, k_des = rp0.k_des, k_loss = rp0.k_loss, + X_total = mean(X_tot), X_red = mean(X_red), + basis_gate_ack = rp0.basis_gate_ack, + P0 = rp0.P0, alpha_P = rp0.alpha_P, k_dam = rp0.k_dam, + Ddot_R = rp0.Ddot_R, c_ext = rp0.c_ext, dt_rd = rp0.dt_rd) end step_radiolysis!(rd, rd.params.dt_rd) diff --git a/tests/checkpoint_io_tests.jl b/tests/checkpoint_io_tests.jl index a1dd14d..f2264c6 100644 --- a/tests/checkpoint_io_tests.jl +++ b/tests/checkpoint_io_tests.jl @@ -7,7 +7,12 @@ include(joinpath(REPO, "export_checkpoint.jl")) # CLI-guarded; functions only tmp = mktempdir() p = SR.CPMParams(N = 20, n_cells_per_species = 2, snapshot_interval = 100) -rp = SR.RadiolysisParams(Nr = 20, Ddot_R = 1.0, c_ext = 1.0) +# basis_gate_ack: this compares blocked output against blocked output (two code +# paths, or a round trip) and asserts no magnitude, so it needs the gated basis +# to step but makes no claim about it. Enumerated in the ack census in +# tests/radiodialysis_basis_gate.jl -- adding a fourth site fails that test. +rp = SR.RadiolysisParams(Nr = 20, Ddot_R = 1.0, c_ext = 1.0, + basis_gate_ack = true) # ---------- transport snapshot: conventions + probe integrity ---------- diff --git a/tests/genealogy_tests.jl b/tests/genealogy_tests.jl index e4cc9a3..bd387bb 100644 --- a/tests/genealogy_tests.jl +++ b/tests/genealogy_tests.jl @@ -143,7 +143,12 @@ end let p = SR.CPMParams(N = 20, n_cells_per_species = 2, snapshot_interval = 100) - rp = SR.RadiolysisParams(Nr = 20, Ddot_R = 1.0, c_ext = 1.0) + # basis_gate_ack: this compares blocked output against blocked output (two code + # paths, or a round trip) and asserts no magnitude, so it needs the gated basis + # to step but makes no claim about it. Enumerated in the ack census in + # tests/radiodialysis_basis_gate.jl -- adding a fourth site fails that test. + rp = SR.RadiolysisParams(Nr = 20, Ddot_R = 1.0, c_ext = 1.0, + basis_gate_ack = true) n = 21 # crosses the mcs % 10 == 1 biomass-rebuild boundary twice st_legacy, rd_legacy, _, _ = redirect_stdout(devnull) do diff --git a/tests/radiodialysis_basis_gate.jl b/tests/radiodialysis_basis_gate.jl new file mode 100644 index 0000000..614491e --- /dev/null +++ b/tests/radiodialysis_basis_gate.jl @@ -0,0 +1,134 @@ +# The RADIODIALYSIS: BLOCKED basis gate, its one exemption, and the R/Julia +# asymmetry it pins. +# +# WHY THIS IS AN ASYMMETRY TEST AND NOT AN AGREEMENT TEST. The two +# implementations deliberately DISAGREE at non-unit X_total, and that has to be +# asserted rather than left to be discovered: +# +# biofilms_radiodialysis.R derives X_red = f_red_active * X_total +# (uptake_rate_of), so U scales with the basis and +# lambda ~ 1/sqrt(U). Asserted there by check 15 of +# analysis/verify_biofilm_depth_profile.R. +# this path REFUSES, because the X_red reaching it is one +# species' occupied sites over ALL interior sites -- +# neither a biomass fraction nor a reducer fraction +# (README.md:344), which is the defect the gate names. +# +# Making this path conformant without repairing that quantity would turn a +# visible defect into an invisible one. So the refusal is asserted here, and +# this test fails the day it is removed -- forcing whoever removes it to confront +# the quantity underneath rather than only the arithmetic on top. +# +# Before this gate existed the block was an ACCIDENT: the additive uptake form +# happens to misbehave at non-unit X_total, and that side effect was standing in +# for a guard. A tripwire that any tidy-up can remove is not a gate. + +@testset "radiodialysis basis gate" begin + + @testset "the default basis still integrates" begin + # Without this the refusal half is vacuous: a gate that refused + # unconditionally would satisfy every assertion below it. + rp = SR.RadiolysisParams(Nr = 20, Ddot_R = 1.0, c_ext = 1.0) + @test rp.X_total == 1.0 + @test rp.basis_gate_ack === false + rd = SR.init_radiolysis(rp; R = 10.0) + SR.step_radiolysis!(rd, 0.5) + @test all(isfinite, rd.c) + end + + @testset "a coupled basis is refused, by name" begin + for xt in (0.065, 0.1, 0.5, 0.9, 1.5, 2.0) + rp = SR.RadiolysisParams(Nr = 20, X_total = xt, + Ddot_R = 1.0, c_ext = 1.0) + rd = SR.init_radiolysis(rp; R = 10.0) + err = try + SR.step_radiolysis!(rd, 0.5); nothing + catch e + sprint(showerror, e) + end + @test err !== nothing + @test occursin("RADIODIALYSIS: BLOCKED", err) + @test occursin("neither a biomass fraction nor a reducer fraction", err) + end + end + + @testset "the refusal precedes any stepping" begin + # If the arithmetic were still doing the blocking, state would move + # before anything threw. + rp = SR.RadiolysisParams(Nr = 20, X_total = 0.065, + Ddot_R = 1.0, c_ext = 1.0) + rd = SR.init_radiolysis(rp; R = 10.0) + c0, s0, t0 = copy(rd.c), copy(rd.s), rd.t + @test_throws ErrorException SR.step_radiolysis!(rd, 0.5) + @test rd.c == c0 + @test rd.s == s0 + @test rd.t == t0 + end + + @testset "the determinism exemption is narrow" begin + rp = SR.RadiolysisParams(Nr = 20, X_total = 0.065, Ddot_R = 1.0, + c_ext = 1.0, basis_gate_ack = true) + rd = SR.init_radiolysis(rp; R = 10.0) + SR.step_radiolysis!(rd, 0.5) # allowed + @test all(isfinite, rd.c) + # and ONLY that symbol opens it + # and the default still refuses + rpx = SR.RadiolysisParams(Nr = 20, X_total = 0.065, Ddot_R = 1.0, + c_ext = 1.0) + @test rpx.basis_gate_ack === false + rdx = SR.init_radiolysis(rpx; R = 10.0) + @test_throws ErrorException SR.step_radiolysis!(rdx, 0.5) + end + + @testset "the ack census: exactly these sites, no more" begin + # An exemption that can be added silently is not an exemption, it is a + # default. Every site opening the gate is enumerated here, so widening + # it is a deliberate edit to this list rather than a line someone adds. + expected = Set([ + "validate_serial.jl", # CPM trajectory determinism + "tests/genealogy_tests.jl", # legacy vs windowed API equivalence + "tests/checkpoint_io_tests.jl", # snapshot/restart round trip + "tests/radiodialysis_basis_gate.jl", # this file, testing the gate + ]) + repo = dirname(@__DIR__) + found = Set{String}() + for (root, _, files) in walkdir(repo) + occursin("/.git", root) && continue + for fname in files + endswith(fname, ".jl") || continue + path = joinpath(root, fname) + # The ASSIGNMENT form only. The guard itself compares with + # `ack === :determinism_only`, and the docstrings name the + # symbol in prose; neither opens the gate for anyone. + if occursin("basis_gate_ack = true", read(path, String)) + push!(found, relpath(path, repo)) + end + end + end + @test found == expected + end + + @testset "the exemption's claim is true: nothing recorded depends on the basis" begin + # validate_serial.jl acknowledges a blocked basis on the grounds that it + # records no radiodialysis quantity. That is a CLAIM, and this is what + # holds it to account: the gated basis enters only through `uptake`, so + # changing the uptake constants by 10x must leave every recorded column + # byte-identical. The day someone records a c- or s-derived quantity, + # this fails and the exemption has to be re-argued. + function recorded(k_ads, k_red) + params = SR.CPMParams(N = 40, n_cells_per_species = 6, + snapshot_interval = 20) + rp = SR.RadiolysisParams(Nr = 40, Ddot_R = 1.0, c_ext = 1.0, + k_ads = k_ads, k_red = k_red, + basis_gate_ack = true) + state, rd, _, _ = SR.run_simulation_coupled(params, rp, 40; seed = 42) + snap = SR.take_snapshot(state, 40) + rows = [(sd.species, sd.total_volume, sd.n_cells, sd.mean_melanin) + for sd in snap.species_data] + return (rows, length(state.cells), rd.m) + end + base = recorded(0.05, 0.02) + perturbed = recorded(0.5, 0.2) # 10x the uptake constants + @test base == perturbed + end +end diff --git a/tests/runtests.jl b/tests/runtests.jl index e2215a3..407c9ab 100644 --- a/tests/runtests.jl +++ b/tests/runtests.jl @@ -28,6 +28,8 @@ end include("deterministic_radiation.jl") end +include("radiodialysis_basis_gate.jl") + @testset "Lifecycle, dose contract, windowed API" begin include("genealogy_tests.jl") end diff --git a/validate_serial.jl b/validate_serial.jl index f59ea93..f33d647 100644 --- a/validate_serial.jl +++ b/validate_serial.jl @@ -19,7 +19,15 @@ end function run(SerialRef::Module, seeds::Vector{Int}) params = SerialRef.CPMParams(N = 40, n_cells_per_species = 6, snapshot_interval = 20) - rp = SerialRef.RadiolysisParams(Nr = 40, Ddot_R = 1.0, c_ext = 1.0) + # basis_gate_ack: this harness reproduces a bit-for-bit CPM trajectory and + # records NO radiodialysis quantity -- its CSV is CPM columns plus rd.m, + # whose ODE has no X_total or X_red in it. The coupled loop reconstructs + # RadiolysisParams with X_total = mean(compute_radial_biomass(...)), which + # RADIODIALYSIS: BLOCKED gates, so stepping it needs this explicit + # acknowledgement. It is not a claim that the basis is valid. Nothing here + # may report c or s. tests/radiodialysis_basis_gate.jl holds this to it. + rp = SerialRef.RadiolysisParams(Nr = 40, Ddot_R = 1.0, c_ext = 1.0, + basis_gate_ack = true) for seed in seeds t0 = time() state, rd, _, _ = SerialRef.run_simulation_coupled(params, rp, 100; seed) From 2c93d1ef913489c95acadc15a505a599b3a9b345 Mon Sep 17 00:00:00 2001 From: aurascoper Date: Fri, 28 Aug 2026 21:59:26 -0500 Subject: [PATCH 08/13] Label the gated basis where it escapes into a file Enforcing the gate surfaced two more coupled call sites, and the census test is what surfaced them rather than a CI failure nobody read. export_checkpoint.jl's CLI stepped the coupled path without the acknowledgement, so `export_checkpoint.jl transport` exited 1 and took coupling/tests/test_julia_interop.py with it. It is a fourth ack site, and unlike the other three it is not a comparison: export_restart_checkpoint writes rd/c and rd/s into an interchange file that other tools read. So the file declares it (rule 4). rd/basis_gate_blocked and rd/basis_gate_note now travel beside the arrays, because a consumer cannot tell a gated basis from a sound one by looking at a Float64 vector. Asserted in both directions, since a marker that is always true is as uninformative as one that is always false. Getting the false direction honest turned up the sharper fact: advance_window! reconstructs the params with X_total = mean(compute_radial_biomass(...)), so an advanced coupled simulation can NEVER export an ungated basis. The unblocked case exists only at zero MCS. The test says so rather than papering over it with a hand-set X_total that the first step would overwrite. The transport snapshot needed the ack only to step; it exports lattice, fields and dose and carries no radiodialysis array, so it gets no marker. Julia 2/34/30/68/48/9/14/22, coupling test_julia_interop 2/2. Co-Authored-By: Claude Opus 5 (1M context) --- export_checkpoint.jl | 21 ++++++++++++++++++++- tests/checkpoint_io_tests.jl | 30 ++++++++++++++++++++++++++++++ tests/radiodialysis_basis_gate.jl | 1 + 3 files changed, 51 insertions(+), 1 deletion(-) diff --git a/export_checkpoint.jl b/export_checkpoint.jl index 9dbe275..2c71a5b 100644 --- a/export_checkpoint.jl +++ b/export_checkpoint.jl @@ -216,6 +216,19 @@ function export_restart_checkpoint(SR, sim, path) f["rd/r_grid"] = sim.rd.r_grid f["rd/c"] = sim.rd.c f["rd/s"] = sim.rd.s + # RULE 4: the producer declares the semantics. c and s here are + # computed on a biomass basis that RADIODIALYSIS: BLOCKED gates + # whenever X_total != 1.0 -- mean(compute_radial_biomass(...)) is one + # species' occupied sites over all interior sites, neither a biomass + # fraction nor a reducer fraction. A consumer cannot tell that by + # looking at the arrays, so the file says it. Asserted by + # tests/checkpoint_io_tests.jl. + f["rd/basis_gate_blocked"] = sim.rd.params.X_total != 1.0 + f["rd/basis_gate_note"] = sim.rd.params.X_total != 1.0 ? + "RADIODIALYSIS: BLOCKED -- c and s were computed on a gated " * + "biomass basis (X_total = $(sim.rd.params.X_total), a site-occupancy " * + "mean). No magnitude read from them is a claim about a biofilm." : + "X_total == 1.0, the standalone default; the basis gate does not apply." buf = IOBuffer() serialize(buf, sim.rng) f["rng/serialized"] = take!(buf) @@ -344,7 +357,13 @@ end function _cli_export(SR, mode, out, cfg, seed, n_mcs) sim = SR.init_coupled_simulation( SR.CPMParams(N = 20, n_cells_per_species = 2, snapshot_interval = 100), - SR.RadiolysisParams(Nr = 20, Ddot_R = 1.0, c_ext = 1.0); seed) + # basis_gate_ack: this CLI exports an interchange artifact rather than + # asserting anything, and every file it writes carries + # rd/basis_gate_blocked so a consumer cannot read c or s without seeing + # that the basis was gated. Enumerated in the ack census in + # tests/radiodialysis_basis_gate.jl. + SR.RadiolysisParams(Nr = 20, Ddot_R = 1.0, c_ext = 1.0, + basis_gate_ack = true); seed) SR.advance_window!(sim, n_mcs) if mode == "transport" export_transport_snapshot(SR, sim, out; config_toml_path = cfg) diff --git a/tests/checkpoint_io_tests.jl b/tests/checkpoint_io_tests.jl index f2264c6..a7e84bf 100644 --- a/tests/checkpoint_io_tests.jl +++ b/tests/checkpoint_io_tests.jl @@ -14,6 +14,36 @@ p = SR.CPMParams(N = 20, n_cells_per_species = 2, snapshot_interval = 100) rp = SR.RadiolysisParams(Nr = 20, Ddot_R = 1.0, c_ext = 1.0, basis_gate_ack = true) +# ---------- the exported file declares a gated basis (rule 4) ---------- + +let + # A consumer reading rd/c cannot tell the basis was gated by looking at the + # array, so the producer says it. Both directions, because a marker that is + # always true and one that is always false are equally uninformative. + # The `false` case must NOT advance: advance_window! reconstructs the params + # with X_total = mean(compute_radial_biomass(...)), so every coupled run + # ends up on the gated basis by construction. That is the finding, not a + # quirk -- an advanced coupled sim can never export an ungated basis. + for (mcs, expected) in ((0, false), (2, true)) + rpx = SR.RadiolysisParams(Nr = 20, Ddot_R = 1.0, c_ext = 1.0, + basis_gate_ack = true) + sim = SR.init_coupled_simulation(p, rpx; seed = 5) + mcs > 0 && SR.advance_window!(sim, mcs) + @test (sim.rd.params.X_total != 1.0) == expected + # export_restart_checkpoint, not the transport snapshot: only the + # restart file carries rd/c and rd/s. The transport snapshot exports + # lattice, fields and dose, none of which touch the gated basis. + path = joinpath(tmp, "gate_$(mcs).h5") + export_restart_checkpoint(SR, sim, path) + h5open(path, "r") do f + @test read(f["rd/basis_gate_blocked"]) == expected + note = read(f["rd/basis_gate_note"]) + @test occursin(expected ? "RADIODIALYSIS: BLOCKED" : + "the basis gate does not apply", note) + end + end +end + # ---------- transport snapshot: conventions + probe integrity ---------- let diff --git a/tests/radiodialysis_basis_gate.jl b/tests/radiodialysis_basis_gate.jl index 614491e..13e6450 100644 --- a/tests/radiodialysis_basis_gate.jl +++ b/tests/radiodialysis_basis_gate.jl @@ -88,6 +88,7 @@ "validate_serial.jl", # CPM trajectory determinism "tests/genealogy_tests.jl", # legacy vs windowed API equivalence "tests/checkpoint_io_tests.jl", # snapshot/restart round trip + "export_checkpoint.jl", # interchange export; labels the file "tests/radiodialysis_basis_gate.jl", # this file, testing the gate ]) repo = dirname(@__DIR__) From 207f2925fe3885fd1c43e9221b8134502defbef4 Mon Sep 17 00:00:00 2001 From: aurascoper Date: Fri, 28 Aug 2026 22:23:56 -0500 Subject: [PATCH 09/13] Declare basis provenance, and make the census parse instead of grep Both of Codex's on 2c93d1e, and both were defects in guards I wrote. PROVENANCE. Inferring "gated" from X_total != 1.0 fails two ways. A coupled state can produce mean(X_tot) == 1.0 legitimately -- every sampled interior site occupied -- colliding with the standalone default and slipping the gate. And c/s are path-dependent, so a run that stepped on an occupancy basis stays gated even after X_total returns to 1.0, while the exported marker would have said false. Rule 4: the producer that installs the basis declares it. RadiolysisState carries basis_from_occupancy, set at each of the three installation sites and never cleared, serialised into the restart file and restored rather than recomputed. The gate and the export marker both read it. Four tests: the collision, the stickiness, a genuinely standalone state still allowed, and -- because the first two set the flag by hand -- that the real coupled loop sets it, so the guard is not protecting a field nobody writes. Confirmed to bite: dropping the provenance term fails two. CENSUS. It grepped for the literal "basis_gate_ack = true", so `basis_gate_ack=true`, extra spacing, a line break or a field assignment all opened the gate with the census still green. It asserted a bound it could not enforce, which is the shape it exists to catch. It parses now, walking the AST for assignments in any spelling, and distinguishes OPENING from PROPAGATION: the reconstruction sites forward `basis_gate_ack = rp.basis_gate_ack`, which justifies nothing new, while a literal true or a value it cannot trace does. Nine fixtures prove it catches each bypass and does not fire on the guard's own comparison, on prose, or on propagation. Confirmed to bite by adding a site in the no-space form. Its scope is stated rather than implied: this is static, so a runtime value splatted in via `RadiolysisParams(; kw...)` is not detectable and the test says so instead of claiming a bound it does not have. Julia 2/34/51/68/48/9/14/22. Co-Authored-By: Claude Opus 5 (1M context) --- biofilms_potts.jl | 28 ++++++- biofilms_potts_jacc.jl | 21 +++++- export_checkpoint.jl | 20 +++-- tests/radiodialysis_basis_gate.jl | 118 ++++++++++++++++++++++++++++-- 4 files changed, 168 insertions(+), 19 deletions(-) diff --git a/biofilms_potts.jl b/biofilms_potts.jl index 4efb439..f583b26 100644 --- a/biofilms_potts.jl +++ b/biofilms_potts.jl @@ -1224,8 +1224,20 @@ mutable struct RadiolysisState m::Float64 # membrane integrity m ∈ [0,1] t::Float64 # simulation time params::RadiolysisParams + # Sticky basis provenance. NOT inferred from X_total: a coupled state can + # legitimately produce mean(X_tot) == 1.0 (every sampled interior site + # occupied), which would collide with the standalone default and slip the + # gate; and c/s are PATH-dependent, so once a step has run on an occupancy + # basis they stay gated even if X_total later returns to 1.0. Rule 4: the + # producer that installs the basis declares it, and the flag never clears. + basis_from_occupancy::Bool end +# Six-argument form: provenance defaults to false, so a standalone state is +# unmarked and only the coupled installers set it. +RadiolysisState(r_grid, c, s, m, t, params) = + RadiolysisState(r_grid, c, s, m, t, params, false) + """ Initialise RadiolysisState: clean interior, intact membrane. """ @@ -1259,7 +1271,8 @@ R = 1.0 cm makes dt_rd = 0.5 genuinely unstable, and this guard is what absorbs it. `biofilms_potts_jacc.jl` carries the identical wrapper for the same reason. """ function step_radiolysis!(rd::RadiolysisState, dt::Float64) - _assert_basis_gate(rd.params.X_total, rd.params.basis_gate_ack) + _assert_basis_gate(rd.params.X_total, rd.params.basis_gate_ack, + rd.basis_from_occupancy) dr = rd.r_grid[2] - rd.r_grid[1] dt_stable = 0.4 * dr^2 / (2.0 * rd.params.D_eff) n_sub = max(1, ceil(Int, dt / dt_stable)) @@ -1302,8 +1315,9 @@ active-reducer fraction without a measured activity fraction; `D-XRED` in `data/calibration/reference_d_requirements.csv` records that as blocked by this units error rather than by missing data. """ -function _assert_basis_gate(X_total::Float64, ack::Bool = false) - X_total == 1.0 && return nothing +function _assert_basis_gate(X_total::Float64, ack::Bool = false, + from_occupancy::Bool = false) + (X_total == 1.0 && !from_occupancy) && return nothing # THE ONE EXEMPTION. validate_serial.jl steps this path only to reproduce a # bit-for-bit CPM trajectory, and records no radiodialysis quantity: its CSV # carries CPM columns plus rd.m, whose ODE (dm/dt = -k_dam*Ddot_R*m) has no @@ -1567,6 +1581,10 @@ function run_simulation_coupled(params::CPMParams, rp::RadiolysisParams, c_ext = rp.c_ext, dt_rd = rp.dt_rd ) + # Declared where the occupancy basis is INSTALLED, and never + # cleared: c and s are path-dependent from here on, so provenance + # cannot be re-derived from the current X_total (rule 4). + rd.basis_from_occupancy = true end step_radiolysis!(rd, rp.dt_rd) @@ -1927,6 +1945,10 @@ function advance_window!(sim::CoupledSimulation, n_mcs::Int) basis_gate_ack = rp.basis_gate_ack, P0 = rp.P0, alpha_P = rp.alpha_P, k_dam = rp.k_dam, Ddot_R = rp.Ddot_R, c_ext = rp.c_ext, dt_rd = rp.dt_rd) + # Declared where the occupancy basis is INSTALLED, and never + # cleared: c and s are path-dependent from here on, so provenance + # cannot be re-derived from the current X_total (rule 4). + rd.basis_from_occupancy = true end step_radiolysis!(rd, rd.params.dt_rd) radial_to_3d!(sim.contaminant, rd.c, rd.r_grid, state.interior, N) diff --git a/biofilms_potts_jacc.jl b/biofilms_potts_jacc.jl index 04d9c3f..127efc2 100644 --- a/biofilms_potts_jacc.jl +++ b/biofilms_potts_jacc.jl @@ -282,8 +282,20 @@ mutable struct RadiolysisState m::Float64 t::Float64 params::RadiolysisParams + # Sticky basis provenance. NOT inferred from X_total: a coupled state can + # legitimately produce mean(X_tot) == 1.0 (every sampled interior site + # occupied), which would collide with the standalone default and slip the + # gate; and c/s are PATH-dependent, so once a step has run on an occupancy + # basis they stay gated even if X_total later returns to 1.0. Rule 4: the + # producer that installs the basis declares it, and the flag never clears. + basis_from_occupancy::Bool end +# Six-argument form: provenance defaults to false, so a standalone state is +# unmarked and only the coupled installers set it. +RadiolysisState(r_grid, c, s, m, t, params) = + RadiolysisState(r_grid, c, s, m, t, params, false) + function init_radiolysis(rp::RadiolysisParams; R::Float64 = 1.0) r_grid = collect(range(0.0, R, length = rp.Nr)) RadiolysisState(r_grid, zeros(rp.Nr), zeros(rp.Nr), 1.0, 0.0, rp) @@ -306,7 +318,8 @@ small. Correct the units to R = 1.0 cm and dt_stable becomes 0.132, n_sub = 4 at which point a port without this guard diverges. """ function step_radiolysis!(rd::RadiolysisState, dt::Float64) - _assert_basis_gate(rd.params.X_total, rd.params.basis_gate_ack) + _assert_basis_gate(rd.params.X_total, rd.params.basis_gate_ack, + rd.basis_from_occupancy) dr = rd.r_grid[2] - rd.r_grid[1] dt_stable = 0.4 * dr^2 / (2.0 * rd.params.D_eff) n_sub = max(1, ceil(Int, dt / dt_stable)) @@ -349,8 +362,9 @@ active-reducer fraction without a measured activity fraction; `D-XRED` in `data/calibration/reference_d_requirements.csv` records that as blocked by this units error rather than by missing data. """ -function _assert_basis_gate(X_total::Float64, ack::Bool = false) - X_total == 1.0 && return nothing +function _assert_basis_gate(X_total::Float64, ack::Bool = false, + from_occupancy::Bool = false) + (X_total == 1.0 && !from_occupancy) && return nothing # THE ONE EXEMPTION. validate_serial.jl steps this path only to reproduce a # bit-for-bit CPM trajectory, and records no radiodialysis quantity: its CSV # carries CPM columns plus rd.m, whose ODE (dm/dt = -k_dam*Ddot_R*m) has no @@ -504,6 +518,7 @@ function run_coupled(; N::Int = 40, n_cells_per_species::Int = 6, basis_gate_ack = rp0.basis_gate_ack, P0 = rp0.P0, alpha_P = rp0.alpha_P, k_dam = rp0.k_dam, Ddot_R = rp0.Ddot_R, c_ext = rp0.c_ext, dt_rd = rp0.dt_rd) + rd.basis_from_occupancy = true # see biofilms_potts.jl end step_radiolysis!(rd, rd.params.dt_rd) diff --git a/export_checkpoint.jl b/export_checkpoint.jl index 2c71a5b..d429e32 100644 --- a/export_checkpoint.jl +++ b/export_checkpoint.jl @@ -223,11 +223,17 @@ function export_restart_checkpoint(SR, sim, path) # fraction nor a reducer fraction. A consumer cannot tell that by # looking at the arrays, so the file says it. Asserted by # tests/checkpoint_io_tests.jl. - f["rd/basis_gate_blocked"] = sim.rd.params.X_total != 1.0 - f["rd/basis_gate_note"] = sim.rd.params.X_total != 1.0 ? + # Read from sticky provenance, not from the current X_total: a coupled + # state can hit mean(X_tot) == 1.0 by coincidence, and c/s stay gated + # once a step has run on an occupancy basis regardless of the value now. + gated = sim.rd.basis_from_occupancy || sim.rd.params.X_total != 1.0 + f["rd/basis_gate_blocked"] = gated + f["rd/basis_from_occupancy"] = sim.rd.basis_from_occupancy + f["rd/basis_gate_note"] = gated ? "RADIODIALYSIS: BLOCKED -- c and s were computed on a gated " * - "biomass basis (X_total = $(sim.rd.params.X_total), a site-occupancy " * - "mean). No magnitude read from them is a claim about a biofilm." : + "biomass basis (X_total = $(sim.rd.params.X_total), " * + "basis_from_occupancy = $(sim.rd.basis_from_occupancy)). " * + "No magnitude read from them is a claim about a biofilm." : "X_total == 1.0, the standalone default; the basis gate does not apply." buf = IOBuffer() serialize(buf, sim.rng) @@ -324,9 +330,13 @@ function restore_restart_checkpoint(SR, path; allow_version_mismatch::Bool = fal read(a["membrane_dose_rate_Gy_s"]), read(f["fields/melanin_drive"])) + # Provenance is restored, not recomputed: a restart of a gated run is + # still gated, and the flag cannot be re-derived from X_total. + rd_prov = haskey(f, "rd/basis_from_occupancy") ? + read(f["rd/basis_from_occupancy"]) : false rd = SR.RadiolysisState(read(f["rd/r_grid"]), read(f["rd/c"]), read(f["rd/s"]), read(a["rd_m"]), - read(a["rd_t"]), rp) + read(a["rd_t"]), rp, rd_prov) rng = deserialize(IOBuffer(read(f["rng/serialized"]))) return SR.CoupledSimulation(state, rd, read(f["contaminant"]), rng, Int(read(a["sim_mcs"])), diff --git a/tests/radiodialysis_basis_gate.jl b/tests/radiodialysis_basis_gate.jl index 13e6450..cf5382e 100644 --- a/tests/radiodialysis_basis_gate.jl +++ b/tests/radiodialysis_basis_gate.jl @@ -80,10 +80,111 @@ @test_throws ErrorException SR.step_radiolysis!(rdx, 0.5) end + @testset "provenance is declared and sticky, not inferred from X_total" begin + # Codex on 2c93d1e. Inferring "gated" from X_total != 1.0 fails twice: + # a coupled state can hit mean(X_tot) == 1.0 by coincidence, and c/s are + # path-dependent so they stay gated after the value moves back. + rp1 = SR.RadiolysisParams(Nr = 20, Ddot_R = 1.0, c_ext = 1.0) + + # (a) THE COLLISION: X_total == 1.0, but the basis came from occupancy. + rd = SR.init_radiolysis(rp1; R = 10.0) + rd.basis_from_occupancy = true + @test rd.params.X_total == 1.0 # indistinguishable by value + @test_throws ErrorException SR.step_radiolysis!(rd, 0.5) + + # (b) STICKINESS: gated, stepped, then the value returns to the default. + rd2 = SR.init_radiolysis( + SR.RadiolysisParams(Nr = 20, X_total = 0.065, Ddot_R = 1.0, + c_ext = 1.0, basis_gate_ack = true); R = 10.0) + rd2.basis_from_occupancy = true + SR.step_radiolysis!(rd2, 0.5) # acked, so it runs + rd2.params = SR.RadiolysisParams(Nr = 20, X_total = 1.0, Ddot_R = 1.0, + c_ext = 1.0) # ack dropped, value reset + @test rd2.params.X_total == 1.0 + @test rd2.basis_from_occupancy # provenance survives + @test_throws ErrorException SR.step_radiolysis!(rd2, 0.5) + + # (c) and a genuinely standalone state is still unmarked and allowed. + rd3 = SR.init_radiolysis(rp1; R = 10.0) + @test rd3.basis_from_occupancy === false + SR.step_radiolysis!(rd3, 0.5) + @test all(isfinite, rd3.c) + end + + @testset "the coupled path actually sets provenance" begin + # (a) and (b) above set the flag by hand; this proves the real coupled + # loop sets it, so the guard is not protecting a field nobody writes. + params = SR.CPMParams(N = 20, n_cells_per_species = 2, + snapshot_interval = 100) + rp = SR.RadiolysisParams(Nr = 20, Ddot_R = 1.0, c_ext = 1.0, + basis_gate_ack = true) + sim = SR.init_coupled_simulation(params, rp; seed = 3) + @test sim.rd.basis_from_occupancy === false # before any step + SR.advance_window!(sim, 2) + @test sim.rd.basis_from_occupancy === true # after the installer runs + end + @testset "the ack census: exactly these sites, no more" begin - # An exemption that can be added silently is not an exemption, it is a - # default. Every site opening the gate is enumerated here, so widening - # it is a deliberate edit to this list rather than a line someone adds. + # Codex on 2c93d1e: the first version grepped for the literal + # "basis_gate_ack = true", so `basis_gate_ack=true`, extra spacing, a + # line break, or a field assignment all opened the gate while the census + # stayed green -- a check asserting a bound it could not enforce. + # + # This parses instead. _opens_gate walks the Julia AST for assignments + # to basis_gate_ack in any syntactic form. + # + # SCOPE, stated rather than implied: this is a STATIC check. A value + # computed at runtime and splatted in (`RadiolysisParams(; kw...)`) is + # not statically detectable and this cannot catch it. What it does + # establish is that no site opens the gate by writing the field + # directly, in any spelling. The known-bad fixtures below prove it + # catches each form rather than passing because it matches nothing. + function _opens_gate(ex) + hit = false + walk(e) = begin + if e isa Expr + if e.head === :kw || e.head === :(=) + lhs, rhs = e.args[1], e.args[2] + names = lhs === :basis_gate_ack || + (lhs isa Expr && lhs.head === :. && + lhs.args[2] == QuoteNode(:basis_gate_ack)) + # PROPAGATION is not OPENING. The reconstruction sites + # write `basis_gate_ack = rp.basis_gate_ack`, forwarding + # a flag someone else already justified; and `= false` + # is the closed default. Anything else -- a literal + # true, or a computed value whose provenance this cannot + # see -- counts as opening the gate and must be listed. + forwards = rhs isa Expr && rhs.head === :. && + rhs.args[2] == QuoteNode(:basis_gate_ack) + names && !forwards && rhs !== false && (hit = true) + end + foreach(walk, e.args) + end + end + walk(ex) + return hit + end + opens(src) = _opens_gate(Meta.parseall(src)) + + # The control: each of these bypassed the old literal grep. + @test opens("RadiolysisParams(basis_gate_ack=true)") + @test opens("RadiolysisParams(basis_gate_ack = true)") + @test opens("RadiolysisParams(Nr = 20,\n basis_gate_ack =\n true)") + @test opens("rd.params.basis_gate_ack = true") + @test opens("x = (basis_gate_ack = true,)") + # and it must not fire on the guard's own comparison or on prose, + # or the census would name every file that mentions the symbol. + @test !opens("ack === :determinism_only") + @test !opens("f(basis_gate_ack)") + @test !opens("# basis_gate_ack = true, in a comment") + @test !opens("s = \"basis_gate_ack = true\"") + # propagation forwards a flag already justified elsewhere; the closed + # default opens nothing. Neither should name a site. + @test !opens("RadiolysisParams(basis_gate_ack = rp.basis_gate_ack)") + @test !opens("basis_gate_ack::Bool = false") + # but a value this cannot trace IS an opening, and must be listed + @test opens("RadiolysisParams(basis_gate_ack = should_ack())") + expected = Set([ "validate_serial.jl", # CPM trajectory determinism "tests/genealogy_tests.jl", # legacy vs windowed API equivalence @@ -98,12 +199,13 @@ for fname in files endswith(fname, ".jl") || continue path = joinpath(root, fname) - # The ASSIGNMENT form only. The guard itself compares with - # `ack === :determinism_only`, and the docstrings name the - # symbol in prose; neither opens the gate for anyone. - if occursin("basis_gate_ack = true", read(path, String)) - push!(found, relpath(path, repo)) + src = read(path, String) + parsed = try + Meta.parseall(src) + catch + continue end + _opens_gate(parsed) && push!(found, relpath(path, repo)) end end @test found == expected From f49fe94051b860998f6e1f3d84da292480912629 Mon Sep 17 00:00:00 2001 From: aurascoper Date: Fri, 28 Aug 2026 22:43:03 -0500 Subject: [PATCH 10/13] Census: keyword shorthand is a bypass, not the disclaimed residual Codex on 207f292. `RadiolysisParams(; basis_gate_ack)` passes a same-named binding as a bare Symbol under an AST :parameters node, so a walker examining only :kw and :(=) never counted the site. The distinction matters and Codex drew it: the field name is STATICALLY PRESENT here, so this is not the runtime-splat residual the last commit disclaimed. I had scoped the limitation to values I could not see, and this one was in plain sight. Four fixtures, including Codex's own `function build(basis_gate_ack); RadiolysisParams(; basis_gate_ack); end`. Confirmed to bite by adding a real site in that form, which the census names. The :parameters branch errs toward NAMING a site: a keyword parameter merely declared with this name fires too. A census that fails loudly on something it cannot classify is right in the direction rule 3 asks for. Julia 2/34/55/68/48/9/14/22. Co-Authored-By: Claude Opus 5 (1M context) --- tests/radiodialysis_basis_gate.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/tests/radiodialysis_basis_gate.jl b/tests/radiodialysis_basis_gate.jl index cf5382e..cf6f599 100644 --- a/tests/radiodialysis_basis_gate.jl +++ b/tests/radiodialysis_basis_gate.jl @@ -157,6 +157,16 @@ forwards = rhs isa Expr && rhs.head === :. && rhs.args[2] == QuoteNode(:basis_gate_ack) names && !forwards && rhs !== false && (hit = true) + elseif e.head === :parameters + # KEYWORD SHORTHAND: `RadiolysisParams(; basis_gate_ack)` + # passes a same-named binding. The field name is + # statically present, so this is NOT the runtime-splat + # residual disclaimed below -- it is a real bypass + # (Codex on 207f292). A bare Symbol under :parameters. + # Errs toward NAMING a site: a keyword parameter + # declared with this name fires too, failing the census + # loudly rather than missing an exemption. + any(a -> a === :basis_gate_ack, e.args) && (hit = true) end foreach(walk, e.args) end @@ -184,6 +194,12 @@ @test !opens("basis_gate_ack::Bool = false") # but a value this cannot trace IS an opening, and must be listed @test opens("RadiolysisParams(basis_gate_ack = should_ack())") + # keyword shorthand: same-named binding, field name still statically + # present, so the census must see it + @test opens("RadiolysisParams(; basis_gate_ack)") + @test opens("function build(basis_gate_ack)\n RadiolysisParams(; basis_gate_ack)\nend") + @test opens("(; basis_gate_ack)") + @test !opens("RadiolysisParams(; Nr, D_eff)") expected = Set([ "validate_serial.jl", # CPM trajectory determinism From ba5c424026a93a1bdc4060decd12f0839df622c7 Mon Sep 17 00:00:00 2001 From: aurascoper Date: Fri, 28 Aug 2026 23:22:18 -0500 Subject: [PATCH 11/13] Refuse unknown checkpoint provenance, and hybrid geometry lists Two of Codex's three on f49fe94. The third (the census collapsing sites to filenames) is deferred to an issue. MISSING PROVENANCE IS UNKNOWN, NOT FALSE. restore_restart_checkpoint defaulted absent rd/basis_from_occupancy to false, which is rule 3 inside a gate: a checkpoint written before the field existed carries no provenance, and if its occupancy-derived X_total happened to land on 1.0 the restored state read as standalone and resumed integrating gated c/s. It now refuses and makes the caller establish the basis, the same shape as the existing allow_version_mismatch. Six assertions: the refusal, both explicit declarations honoured, and a current checkpoint still restoring untouched, so the refusal is not simply blocking everything. HYBRID PARAMETER LISTS. bc_coef dispatched on geom while the interior rows consumed whatever face weights parms carried, and nothing tied the two together. Codex's example: take a default_parms() result, set geom = "slab", supply k_L, and get a slab boundary over cylindrical diffusion -- a plausible operator rather than a refusal. This is the same defect as the original bc_coef finding one layer down; validating half a contract proves nothing about the other half. The weights are now validated against face_weights(r_grid, geom) rather than recomputed, because recomputing would discard what the caller wrote and hide the mistake instead of naming it. geom is checked first so an unrecognised value still gets the named refusal rather than face_weights()'s generic one. Check 16 covers both directions and a single corrupted weight, not only a wholesale swap, and requires both presets to still run. Confirmed to bite by deleting the guard: all three bad lists go FALSE while presets_still_run stays TRUE. R 16/16, Julia 2/34/55/68/54/9/14/22, numpy 3/3. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/model-contracts.yml | 2 +- analysis/verify_biofilm_depth_profile.R | 28 ++++++++++++++++++ biofilms_radiodialysis.R | 27 ++++++++++++++++++ export_checkpoint.jl | 30 +++++++++++++++++-- tests/checkpoint_io_tests.jl | 38 +++++++++++++++++++++++++ 5 files changed, 121 insertions(+), 4 deletions(-) diff --git a/.github/workflows/model-contracts.yml b/.github/workflows/model-contracts.yml index 8ac85f6..ca0c2df 100644 --- a/.github/workflows/model-contracts.yml +++ b/.github/workflows/model-contracts.yml @@ -154,7 +154,7 @@ jobs: - name: Assert the verifier actually ran everything shell: Rscript {0} run: | - EXPECTED_CHECKS <- 15L + EXPECTED_CHECKS <- 16L if (!file.exists("r-verify.json")) stop("no receipt: the verifier did not run to completion") r <- jsonlite::fromJSON("r-verify.json", simplifyVector = TRUE) diff --git a/analysis/verify_biofilm_depth_profile.R b/analysis/verify_biofilm_depth_profile.R index b87b810..104e133 100644 --- a/analysis/verify_biofilm_depth_profile.R +++ b/analysis/verify_biofilm_depth_profile.R @@ -167,6 +167,34 @@ local({ refused, accepts[1], accepts[2])) }) +# --- 1b2. cached weights must match the declared geometry ------------------ +# Codex on f49fe94. bc_coef dispatched on geom while the interior rows used +# whatever weights parms carried, so geom = "slab" over cylindrical weights +# assembled a hybrid operator instead of refusing. Both directions, and the +# well-formed presets must still run or "refuses everything" would pass. +local({ + hybrid <- function(base, geom) { p <- base; p$geom <- geom; p } + y_of <- function(p) c(rep(0, length(p$r_grid)), rep(0, length(p$r_grid)), 1) + caught <- function(p) tryCatch({ radiodialysis_rhs(0, y_of(p), p); FALSE }, + error = function(e) + grepl("cached face weights do not match geom", + conditionMessage(e), fixed = TRUE)) + # cylindrical parms relabelled slab (Codex's example), and the reverse + cyl_as_slab <- hybrid(default_parms(), "slab"); cyl_as_slab$k_L <- 1e-3 + slab_as_cyl <- hybrid(slab_parms(X_total = 1.0), "cylindrical") + slab_as_cyl$P0 <- 0.01; slab_as_cyl$alpha_P <- 0.02 + # and a single corrupted weight, not just a wholesale swap + one_bad <- slab_parms(X_total = 1.0); one_bad$w_plus[5] <- 1.001 + bad_ok <- caught(cyl_as_slab) && caught(slab_as_cyl) && caught(one_bad) + good_ok <- all(vapply(list(default_parms(), slab_parms(X_total = 1.0)), + function(p) tryCatch({ radiodialysis_rhs(0, y_of(p), p); TRUE }, + error = function(e) FALSE), logical(1))) + report(bad_ok && good_ok, + "cached face weights are validated against the declared geom", + sprintf("cyl_as_slab=%s slab_as_cyl=%s one_weight=%s presets_still_run=%s", + caught(cyl_as_slab), caught(slab_as_cyl), caught(one_bad), good_ok)) +}) + # --- 1c. the reducing basis is a FRACTION of X_total ----------------------- # Codex P1 on #19. X_red was fixed at 0.3 while X_total was caller-declared, # so U did not scale with the declared basis and a corrected low X_total put diff --git a/biofilms_radiodialysis.R b/biofilms_radiodialysis.R index e711568..cf8b2a9 100644 --- a/biofilms_radiodialysis.R +++ b/biofilms_radiodialysis.R @@ -58,6 +58,33 @@ radiodialysis_rhs <- function(t, y, parms) { Nr <- length(r_grid) dr <- r_grid[2] - r_grid[1] + # The weights are CACHED in parms while geom is declared separately, and + # nothing tied the two together: bc_coef dispatched on geom while the + # interior rows consumed whatever weights the list happened to carry. A + # caller editing a default_parms() result to geom = "slab" and supplying + # k_L got a plausible HYBRID -- slab boundary, cylindrical diffusion -- + # instead of a refusal (Codex on f49fe94). Same defect as the bc_coef + # dispatch, one layer down: validating half a contract proves nothing. + # + # Validated, not silently recomputed. Recomputing would quietly discard + # what the caller wrote and hide the mistake; refusing names it (rule 3). + # Node 1 is NA in both by construction, so it is excluded. + # + # geom is checked FIRST so an unrecognised value gets this function's named + # refusal rather than face_weights()'s generic stopifnot one line later. + if (!(geom %in% c("slab", "cylindrical"))) + stop("radiodialysis_rhs(): unsupported geom ", sQuote(geom), + ". Supported: 'slab', 'cylindrical'. An unrecognised geometry ", + "must refuse, not fall through to one of them.", call. = FALSE) + .w <- face_weights(r_grid, geom) + if (!isTRUE(all.equal(w_plus[-1], .w$w_plus[-1])) || + !isTRUE(all.equal(w_minus[-1], .w$w_minus[-1]))) + stop("radiodialysis_rhs(): cached face weights do not match geom = ", + sQuote(geom), ". The parameter list declares one geometry and ", + "carries another's weights, which would assemble a hybrid ", + "operator. Rebuild it with default_parms() or slab_parms() rather ", + "than editing geom in place.", call. = FALSE) + c_vec <- y[seq_len(Nr)] s_vec <- y[Nr + seq_len(Nr)] m_val <- y[2 * Nr + 1] diff --git a/export_checkpoint.jl b/export_checkpoint.jl index d429e32..33bf576 100644 --- a/export_checkpoint.jl +++ b/export_checkpoint.jl @@ -249,7 +249,8 @@ Rebuild a `CoupledSimulation` whose continuation is bit-identical to the unbroken run. Refuses a Julia-version mismatch unless overridden (the RNG byte stream is version-pinned). """ -function restore_restart_checkpoint(SR, path; allow_version_mismatch::Bool = false) +function restore_restart_checkpoint(SR, path; allow_version_mismatch::Bool = false, + declare_basis_from_occupancy::Union{Bool,Nothing} = nothing) h5open(path, "r") do f a = attributes(f) ver = read(a["julia_version"]) @@ -332,8 +333,31 @@ function restore_restart_checkpoint(SR, path; allow_version_mismatch::Bool = fal # Provenance is restored, not recomputed: a restart of a gated run is # still gated, and the flag cannot be re-derived from X_total. - rd_prov = haskey(f, "rd/basis_from_occupancy") ? - read(f["rd/basis_from_occupancy"]) : false + # + # MISSING IS UNKNOWN, NOT FALSE (Codex on f49fe94). A checkpoint written + # before this field existed carries no provenance, and defaulting it to + # false is rule 3 inside a gate: if that file's occupancy-derived + # X_total happened to land on 1.0, the restored state would read as + # standalone and resume integrating gated c/s. So refuse, and make the + # caller establish the basis explicitly -- the same shape as + # allow_version_mismatch above. + rd_prov = if haskey(f, "rd/basis_from_occupancy") + read(f["rd/basis_from_occupancy"]) + elseif declare_basis_from_occupancy !== nothing + declare_basis_from_occupancy + else + error(""" + restart checkpoint predates rd/basis_from_occupancy, so whether + its c and s were computed on a gated biomass basis is UNKNOWN, + not false. X_total alone cannot settle it: an occupancy basis + can coincide with the standalone default of 1.0. + + Re-run, or resume with declare_basis_from_occupancy=true/false + once you have established which it was. Passing `true` is + always the safe direction; it gates a state that may not have + needed it, rather than releasing one that did. + """) + end rd = SR.RadiolysisState(read(f["rd/r_grid"]), read(f["rd/c"]), read(f["rd/s"]), read(a["rd_m"]), read(a["rd_t"]), rp, rd_prov) diff --git a/tests/checkpoint_io_tests.jl b/tests/checkpoint_io_tests.jl index a7e84bf..f55fbaf 100644 --- a/tests/checkpoint_io_tests.jl +++ b/tests/checkpoint_io_tests.jl @@ -14,6 +14,44 @@ p = SR.CPMParams(N = 20, n_cells_per_species = 2, snapshot_interval = 100) rp = SR.RadiolysisParams(Nr = 20, Ddot_R = 1.0, c_ext = 1.0, basis_gate_ack = true) +# ---------- a checkpoint with no provenance is UNKNOWN, not ungated ---------- + +let + # Codex on f49fe94: defaulting missing provenance to false is rule 3 inside + # a gate. Simulate a pre-provenance file by deleting the dataset. + rpx = SR.RadiolysisParams(Nr = 20, Ddot_R = 1.0, c_ext = 1.0, + basis_gate_ack = true) + sim = SR.init_coupled_simulation(p, rpx; seed = 11) + SR.advance_window!(sim, 2) + legacy = joinpath(tmp, "legacy_no_provenance.h5") + export_restart_checkpoint(SR, sim, legacy) + h5open(legacy, "r+") do f + delete_object(f, "rd/basis_from_occupancy") + end + @test !h5open(g -> haskey(g, "rd/basis_from_occupancy"), legacy, "r") + + # refuses rather than assuming + err = try + restore_restart_checkpoint(SR, legacy); nothing + catch e + sprint(showerror, e) + end + @test err !== nothing + @test occursin("UNKNOWN", err) + + # and both explicit declarations are honoured + simT = restore_restart_checkpoint(SR, legacy; declare_basis_from_occupancy = true) + @test simT.rd.basis_from_occupancy === true + simF = restore_restart_checkpoint(SR, legacy; declare_basis_from_occupancy = false) + @test simF.rd.basis_from_occupancy === false + + # a CURRENT checkpoint still restores with no declaration needed, or the + # refusal above would just be blocking everything. + current = joinpath(tmp, "current_provenance.h5") + export_restart_checkpoint(SR, sim, current) + @test restore_restart_checkpoint(SR, current).rd.basis_from_occupancy === true +end + # ---------- the exported file declares a gated basis (rule 4) ---------- let From 3daef5f848203d409184db126b5bee116328edaa Mon Sep 17 00:00:00 2001 From: aurascoper Date: Sat, 29 Aug 2026 00:57:55 -0500 Subject: [PATCH 12/13] Ledger what the merge gate cannot tell apart, and stop sharing checkouts GATE-06. scripts/preflight_merge.sh decides on select(.isResolved | not) and nothing else, so resolved-because-fixed and resolved-because-deferred produce one signal. Demonstrated on this PR: thread PRRT_kwDONEeyC86dXE6w ("Count acknowledgement sites rather than files") is resolved as deferred to #21, not fixed, and the gate now reads that as addressed. Filed as #22, with a register that only downgrades an OPEN thread so `resolved` reverts to meaning one thing and an unverifiable deferral refuses by construction rather than by an else branch. Verdict restate, not requalify: the claim becomes false as written rather than true under a narrower scope, which is why GATE-01 used restate when the set of blocking cases changed. The row says two unflattering things on purpose -- that #22 would be the first gate change to make it refuse LESS, and that unlike GATE-01..05 this was not caught by external review, so it should not be dressed as a catch. CLAUDE.md gains a fourth section. Two commits landed on the wrong branch in one session on 2026-08-28: 0af4505 appeared on feat/slab-depth-geometry mid-session, and fix/ledger-guard-document-scope was checked out under the session so census commit 5b100d7 landed there instead, recovered as f49fe94. Both times `git push origin feat/slab-depth-geometry` printed "Everything up-to-date" while HEAD was a commit ahead, because the branch named in the push was not the branch checked out -- truthful and useless. A mitigation already existed, in a plan scoped to one unrelated piece of work, and a worktree was even open at the time (Biofilms-v11). The tooling was in use and the practice was not, which is the same shape as the pipefail rule living in one script while every instance of the defect was typed by hand. It is standing practice now, and carried to other repos by a user memory. test_claims_ledger.py 15 passed. Co-Authored-By: Claude Opus 5 (1M context) --- CLAUDE.md | 33 +++++++++++++++++++++++++++++++++ data/claims_ledger.csv | 1 + 2 files changed, 34 insertions(+) diff --git a/CLAUDE.md b/CLAUDE.md index bf88150..4587bd9 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -40,3 +40,36 @@ first. Check what the control actually contains before trusting what it reports. State what was run and what it returned. "Conservation to nine significant figures" was true of one check and false about what that check verified; the masses it compared differed by 3%. If a suite skipped, say it skipped. + +## Work in a worktree, not the shared checkout + +Hunter works in these trees at the same time you do. A shared checkout means a +shared branch pointer, and a branch can change under a running session between +one command and the next — so a commit lands wherever `HEAD` happens to point, +not where the work belongs. + +This happened **twice in one session** on 2026-08-28, during PR #19: + +- `0af4505`, Hunter's claims-ledger work, appeared on `feat/slab-depth-geometry` + mid-session. Recovered onto `fix/claims-ledger-delete-verdicts`. +- `fix/ledger-guard-document-scope` was checked out under the running session, + and census commit `5b100d7` landed on it instead of the feature branch. + Recovered as `f49fe94`. + +Both times the confirming command lied. `git push origin feat/slab-depth-geometry` +printed **`Everything up-to-date`** while local `HEAD` was a commit ahead — because +the branch *named* in the push was not the branch *checked out*, and the named one +genuinely had nothing new. **After committing, treat `Everything up-to-date` as +evidence the commit went somewhere else**, not as confirmation of anything. + +Start in a worktree: `EnterWorktree`, or `git worktree add`. Note that +`worktree.baseRef` defaults to `fresh`, which branches from `origin/master`; +stacked work needs `head`. That a worktree already exists is not the rule — one +did (`Biofilms-v11`, on `fix/figure-artifacts-and-v11`) while both collisions +happened. The tooling was in use and the practice was not, because the mitigation +had been written into one plan rather than adopted as standing practice. A +mitigation scoped to a single piece of work is not in force for the next session. + +Read the branch in the same snapshot as the commit — `git status -sb` alongside +`git log -1`, one command — for the same reason every other check here is taken +atomically. A push's exit status is not evidence of where a commit landed. diff --git a/data/claims_ledger.csv b/data/claims_ledger.csv index c0d65df..6368f3c 100644 --- a/data/claims_ledger.csv +++ b/data/claims_ledger.csv @@ -442,6 +442,7 @@ REFINE-10,coupling/scripts/subvoxel_refinement.py,ratio refusal ordering,"Every GATE-04,scripts/test_preflight_merge.sh,cursor contract,The pagination control exercises real pagination.,capability,,,false,true,scripts/test_preflight_merge.sh,code,requalify,The fake refuses a query missing $endCursor or pageInfo before emitting page two,"The fake emitted page two on seeing --paginate alone, so deleting $endCursor or the pageInfo selection left the control green while real pagination silently stopped - the same fail-open shape the pagination fix existed to close, reproduced inside its own test. Caught by external review (Codex, PR #12), the third finding against the merge gate. Removing either element now fails the control." REFD-08,calibration/scripts/reference_d_status.py,phrase override ordering,Each problem is classified by the field it is about.,capability,,,false,true,calibration/tests/test_reference_d_requirements.py,code,restate,Relation phrases are matched only AFTER the field is identified,"THE SAME HOLE ONE SIZE SMALLER. The scope-mismatch phrase was tested BEFORE reading the message head, so approval_source_id = 'does not match the conditions' produced a provenance-only refusal that the checklist reported as criterion 6, leaving criterion 8 met. An echoed identifier could still redirect the institutional checklist. Caught by external review (Codex, PR #12) on the fix for REFD-07. A phrase may now narrow a message already known to be about a field; it may never decide which field a message is about." GATE-05,scripts/test_preflight_merge.sh,cursor binding,The pagination fake enforces the GraphQL cursor contract.,capability,,,false,true,scripts/test_preflight_merge.sh,code,restate,"The fake requires after:$endCursor on the connection, not just the declaration","Checking that $endCursor appeared ANYWHERE was satisfied by the variable declaration alone: remove the after: binding from the connection and cursor updates cannot advance reviewThreads, so the real gate sticks on page one while the control stays green. Codex verified this by removing only the binding and observing all nine controls pass. Fourth finding against the merge gate." +GATE-06,scripts/preflight_merge.sh,resolution disposition,A resolved review thread means the finding was addressed.,capability,,,false,true,scripts/preflight_merge.sh:183,code,restate,"A committed deferrals register: a thread stays OPEN and a well-formed row naming an OPEN issue prints DEFERRED, while anything unverifiable refuses.","The gate decides on select(.isResolved | not) and nothing else, so resolved-because-fixed and resolved-because-deferred produce one signal. Resolving a deferred finding makes it print 'Clear to merge' while a real finding is outstanding. Demonstrated on this PR: thread PRRT_kwDONEeyC86dXE6w (comment 3885494705, 'Count acknowledgement sites rather than files') was resolved as deferred to issue #21, not fixed. Filed as #22, with a register that only downgrades an OPEN thread so that 'resolved' reverts to meaning exactly one thing and an unverifiable deferral falls back to refusing by construction rather than by an else branch. TWO THINGS THIS ROW MUST NOT DRESS UP. First, this would be the first change to the gate that makes it refuse LESS: every prior GATE row closed a fail-open, this opens a narrow typed one, and what holds it is three controls -- a row naming a different thread must refuse, and the open-issue-clears path must be exercised through the real gh, since an implementation that hard-refuses every deferral passes every other control. Second, provenance differs from GATE-01..GATE-05, which were all caught by external review (Codex). This one was found by reasoning about the gate while using it, not by a review catching a defect in flight." REFD-09,calibration/scripts/reference_d_status.py,condition id quoting,Each problem is classified by the field it is about.,capability,,,false,true,calibration/tests/test_reference_d_requirements.py,code,restate,"The prefix strip handles both repr quote styles, so the id cannot steer the checklist","THE IDENTIFIER WAS DATA THAT COULD STEER THE VERDICT. approval.problems formats the growth condition id with !r, and Python switches to double quotes when the value contains an apostrophe - so an ordinary id like ""O'Brien-1"" produced a prefix the classifier's regex could not strip. The head became the word 'condition', nothing classified, and an unset source id was reported as a scope failure plus an UNMAPPED refusal. Caught by external review (Codex, PR #12), the third successive narrowing of the same defect: values reaching a decision they should have no part in." REFD-10,calibration/biofilm_calibration/approval.py,Refusal.subject,Each problem is classified by the field it is about.,capability,4,narrowing fixes,false,true,calibration/scripts/reference_d_status.py,code,delete,approval.classified carries the subject; nothing parses the prose,"FOUR FIXES NARROWED ONE DEFECT WITHOUT CLOSING IT. The checklist recovered each refusal's field by stripping a 'growth condition ' prefix and taking the first token. It failed on an echoed field value matching a criterion pattern; then on a longer echoed value impersonating a relation phrase; then on an apostrophe flipping repr to double quotes; then on an id containing BOTH quote characters, which repr must escape and no delimiter pair can bracket. Each fix was correct and each left the shape intact one size smaller, because the defect was never the regex - it was asking a sentence built from unrestricted CSV data to say which field it was about. Caught across four rounds of external review (Codex, PR #12), which recommended the structural fix. approval.classified now returns Refusal(subject, text) and problems() is derived from it, so the two cannot disagree and a condition id is data again. Verified: fourteen single-field breakages give identical correct verdicts under gid 'GC1' and under gid O'Brien ""lab""." PILOT-11,coupling/scripts/openmc_nested_pilot.py,resolve_output_dir,The refusal is the only route to the canonical directory.,capability,,,false,true,coupling/tests/test_pilot_seeding.py,code,restate,"The guard is on the DESTINATION, not on the --publish flag","IT WAS NOT THE ONLY ROUTE. --publish is one way into data/calibration/; --outdir data/calibration is another, and it reached the canonical tables without passing any check at all - so a budget-exhausted run could still overwrite published evidence with partial rows, which is exactly the loss the refusal exists to prevent. Caught by external review (Codex, PR #12) on the fix from PILOT-09. A flag is a proxy for the thing that matters; the directory IS the thing that matters. The check now resolves the target and refuses it, along with any subdirectory of it, by whichever route." From 63be89bb2911511db6629cd9e92289af3cca6882 Mon Sep 17 00:00:00 2001 From: aurascoper Date: Sat, 29 Aug 2026 08:33:19 -0500 Subject: [PATCH 13/13] The worktree rule prescribed the defect it warns about Codex on 3daef5f, and it is pointing at the rule itself. The closing directive said to read the branch "in the same snapshot as the commit -- `git status -sb` alongside `git log -1`, one command". Those are two git invocations. Joining them with `&&` sequences them; it does not make them atomic. Another process moving HEAD between the two -- the exact scenario the section describes -- lets `git status -sb` report the old branch while `git log -1` reports the new commit, so the check can falsely pair a commit with the intended branch. One invocation that actually is one: git status --porcelain=v2 --branch reporting branch.oid and branch.head together from a single read. The paragraph now says that an earlier draft of itself prescribed the two-command check, because a rule that quietly corrects its own history teaches less than one that says where it was wrong. The same correction is applied to the agent-sessions-use-worktrees memory, which carried the same line. Worth recording plainly: this section was written in the same session that used `git status -sb; git log --oneline -1` as its verification pattern, and neither the writing nor the using noticed. External review did. Co-Authored-By: Claude Opus 5 (1M context) --- CLAUDE.md | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 4587bd9..3061c22 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -70,6 +70,13 @@ happened. The tooling was in use and the practice was not, because the mitigatio had been written into one plan rather than adopted as standing practice. A mitigation scoped to a single piece of work is not in force for the next session. -Read the branch in the same snapshot as the commit — `git status -sb` alongside -`git log -1`, one command — for the same reason every other check here is taken -atomically. A push's exit status is not evidence of where a commit landed. +Read the branch and the commit in **one invocation**: + + git status --porcelain=v2 --branch + +which reports `branch.oid` and `branch.head` together, from a single read. Two +git commands on one shell line are not one snapshot — `&&` sequences them, it +does not make them atomic, so `git status -sb` can report the old branch while +`git log -1` reports the new commit. That is this section's own failure, and an +earlier draft of this paragraph prescribed exactly that two-command check. +A push's exit status is not evidence of where a commit landed either.