diff --git a/biofilms_potts.jl b/biofilms_potts.jl index 9aa37d7..918965a 100644 --- a/biofilms_potts.jl +++ b/biofilms_potts.jl @@ -2288,6 +2288,698 @@ function advance_window!(sim::CoupledSimulation, n_mcs::Int; return sim end +# ============================================================ +# 12b. Two-Way State-Dependent Attenuation, Two-Group Cascade +# & Bilinear Poison Burnout Surrogate (BiofilmsPottsRad) +# ============================================================ +module BiofilmsPottsRad + +using Random +using Statistics + +export Params, State, init_state, local_attenuation, update_radiation_field!, + dH_radiation, step_cpm!, step_fields!, update_radial_metrics!, + compute_derived_layers, simulate_xenon_pit_hysteresis, + export_observer_payload, run_control_audit + +Base.@kwdef struct Params + N::Int = 48 + N_types::Int = 5 + R_wall::Float64 = 24.0 + T::Float64 = 12.0 + T_potts::Float64 = 1.0 + λ_vol::Float64 = 2.5 + V_target::Vector{Int} = Int[0, 16, 16, 18, 22, 24] + I0::Float64 = 1.0 + κ::Float64 = 1.0 + κ_water::Float64 = 0.35 + κ_void::Float64 = 0.02 + κ_cell::Float64 = 1.0 + κ_mel::Float64 = 2.5 + κ_species::Vector{Float64} = Float64[0.35, 0.45, 0.50, 1.40, 2.20] + Σ_water::Float64 = 0.35 + Σ_void::Float64 = 0.01 + Σ_scatter::Vector{Float64} = Float64[0.25, 0.30, 0.30, 0.45, 0.60] + κ_therm::Float64 = 1.10 + beta_ion::Float64 = -4.0 + β_ion::Vector{Float64} = Float64[0.0, 0.0, -0.5, -4.0, -6.5] + alpha_M_species::Vector{Float64} = Float64[0.05, 0.15, 0.25, 0.45, 0.85] + D_C::Float64 = 0.15 + γ_prod::Float64 = 0.08 + γ_delay::Float64 = 0.08 + λ_decay::Float64 = 0.015 + σ_burn::Float64 = 0.45 +end + +mutable struct State + σ::Matrix{Int} + cell_type::Dict{Int,Int} + cell_vol::Dict{Int,Int} + I_fast::Matrix{Float64} + I_therm::Matrix{Float64} + radiation::Matrix{Float64} + M_damage::Matrix{Float64} + C_inhib::Matrix{Float64} + mean_r_history::Vector{Vector{Float64}} + void_fraction_history::Vector{Float64} +end + +""" + local_attenuation(state::State, p::Params, i::Int, j::Int) -> (κ, Σ_s) + +Return the local fast-group attenuation coefficient `κ` and fast-to-thermal +moderation cross-section `Σ_s` at lattice site `(i, j)`. +When `p.κ == 0.0` or all phase attenuation coefficients are zero, both +attenuation and fast-to-thermal scattering vanish identically across water, +voids, and occupied cells. +""" +@inline function local_attenuation(state::State, p::Params, i::Int, j::Int) + if p.κ == 0.0 || (p.κ_water == 0.0 && p.κ_void == 0.0 && all(==(0.0), p.κ_species)) + return (0.0, 0.0) + end + cid = state.σ[i, j] + if cid == 0 + return (p.κ * p.κ_water, p.κ * p.Σ_water) + elseif cid == -1 + return (p.κ * p.κ_void, p.κ * p.Σ_void) + else + sp = state.cell_type[cid] + return (p.κ * p.κ_species[sp], p.κ * p.Σ_scatter[sp]) + end +end + +""" + init_rad_state(p::Params; seed::Int = 42, void_fraction::Float64 = 0.08) -> State + +Initialize a 2D lattice state with water (`0`), radiolytic gas/steam voids (`-1`), +and `p.N_types` species distributed in radial shells around the central source. +""" +function init_rad_state(p::Params; seed::Int = 42, void_fraction::Float64 = 0.08) + rng = MersenneTwister(seed) + N = p.N + σ = zeros(Int, N, N) + cell_type = Dict{Int,Int}() + cell_vol = Dict{Int,Int}() + cx, cy = (N + 1) / 2.0, (N + 1) / 2.0 + max_r = hypot(cx - 1.0, cy - 1.0) + + next_cid = 1 + for i in 2:3:(N - 1), j in 2:3:(N - 1) + r_norm = hypot(i - cx, j - cy) / max(max_r, 1.0) + if r_norm < 0.12 + continue + end + u = rand(rng) + if u < void_fraction + σ[i, j] = -1 + continue + elseif u < void_fraction + 0.22 + continue + end + sp = clamp( 1 + floor(Int, ((i + j) % p.N_types)), 1, p.N_types ) + cid = next_cid + next_cid += 1 + cell_type[cid] = sp + vol = 0 + for di in 0:1, dj in 0:1 + ii = clamp(i + di, 1, N) + jj = clamp(j + dj, 1, N) + if σ[ii, jj] == 0 + σ[ii, jj] = cid + vol += 1 + end + end + if vol > 0 + cell_vol[cid] = vol + else + delete!(cell_type, cid) + end + end + + state = State( + σ, + cell_type, + cell_vol, + zeros(Float64, N, N), + zeros(Float64, N, N), + zeros(Float64, N, N), + zeros(Float64, N, N), + zeros(Float64, N, N), + [Float64[] for _ in 1:p.N_types], + Float64[], + ) + update_radiation_field!(state, p) + update_radial_metrics!(state, p) + return state +end + +const init_state = init_rad_state + +""" + update_radiation_field!(state::State, p::Params) + +Compute the two-group radiation field (`I_fast`, `I_therm`, and `radiation`) +by integrating state-dependent attenuation and fast-to-thermal build-up along +discrete rays from `(cx, cy)`. + +Conservation & Normalization Contract: +- All flux arrays (`I_fast`, `I_therm`, `radiation`) are dimensionless relative + synthetic surrogates (`rel_flux`, normalized to source `I0`), NOT measured + physical dose (`Gy` or `Gy/s`). +- For non-negative attenuation (`κ_loc >= 0`) and scattering (`Σ_loc >= 0`), + `0 <= I_fast[i, j] <= p.I0` and `I_therm[i, j] >= 0` everywhere, with exact + two-group additivity `radiation[i, j] == I_fast[i, j] + I_therm[i, j]`. +- Effective moderation cross-section is bounded by local fast removal + (`0 <= Σ_eff = clamp(Σ_loc, 0.0, max(0.0, κ_loc)) <= κ_eff`) with + `κ_th_eff = max(0.0, p.κ * p.κ_therm) >= 0`, so total two-group flux never + multiplies above the unattenuated source: + `0 <= I_fast[i, j] + I_therm[i, j] <= p.I0`. +""" +function update_radiation_field!(state::State, p::Params) + N = p.N + if p.κ == 0.0 || (p.κ_water == 0.0 && p.κ_void == 0.0 && all(==(0.0), p.κ_species)) + fill!(state.I_fast, p.I0) + fill!(state.I_therm, 0.0) + fill!(state.radiation, p.I0) + return nothing + end + cx, cy = (N + 1) / 2.0, (N + 1) / 2.0 + κ_th_eff = max(0.0, p.κ * p.κ_therm) + + for j in 1:N, i in 1:N + dx = i - cx + dy = j - cy + dist = hypot(dx, dy) + steps = max(1, ceil(Int, dist)) + ds = (dist / steps) / N + + optical_depth_fast = 0.0 + I_therm_acc = 0.0 + + for s in 1:steps + frac = (s - 0.5) / steps + ix = clamp(round(Int, cx + frac * dx), 1, N) + iy = clamp(round(Int, cy + frac * dy), 1, N) + κ_loc, Σ_loc = local_attenuation(state, p, ix, iy) + κ_eff = max(0.0, κ_loc) + Σ_eff = clamp(Σ_loc, 0.0, κ_eff) + I_f_in = p.I0 * exp(-optical_depth_fast) + optical_depth_fast += κ_eff * ds + transfer = if abs(κ_eff - κ_th_eff) * ds > 1e-12 + Σ_eff * I_f_in * (exp(-κ_th_eff * ds) - exp(-κ_eff * ds)) / (κ_eff - κ_th_eff) + else + Σ_eff * I_f_in * ds * exp(-κ_th_eff * ds) + end + I_therm_acc = I_therm_acc * exp(-κ_th_eff * ds) + transfer + end + + I_f = p.I0 * exp(-optical_depth_fast) + state.I_fast[i, j] = I_f + state.I_therm[i, j] = I_therm_acc + state.radiation[i, j] = state.I_fast[i, j] + state.I_therm[i, j] + end + return nothing +end + +""" + dH_radiation(state::State, p::Params, i_src::Int, j_src::Int, i_tgt::Int, j_tgt::Int) -> Float64 + +Radiation Hamiltonian bias for copying `σ[i_src, j_src]` into `(i_tgt, j_tgt)`. +""" +@inline function dH_radiation(state::State, p::Params, i_src::Int, j_src::Int, i_tgt::Int, j_tgt::Int) + cid_src = state.σ[i_src, j_src] + cid_tgt = state.σ[i_tgt, j_tgt] + cid_src == cid_tgt && return 0.0 + + β_src = cid_src > 0 ? p.β_ion[state.cell_type[cid_src]] : 0.0 + β_tgt = cid_tgt > 0 ? p.β_ion[state.cell_type[cid_tgt]] : 0.0 + + ΔI = state.radiation[i_tgt, j_tgt] - state.radiation[i_src, j_src] + return (β_src - β_tgt) * ΔI +end + +""" + step_cpm!(state::State, p::Params; n_trials::Int = p.N * p.N, rng::AbstractRNG = Random.default_rng()) + +Execute `n_trials` Metropolis copy attempts preserving exact `cell_vol` counts +and updating the coupled two-group radiation field. +""" +function step_cpm!(state::State, p::Params; n_trials::Int = p.N * p.N, rng::AbstractRNG = Random.default_rng()) + N = p.N + offsets = ((-1, 0), (1, 0), (0, -1), (0, 1), (-1, -1), (-1, 1), (1, -1), (1, 1)) + for _ in 1:n_trials + i_tgt = rand(rng, 1:N) + j_tgt = rand(rng, 1:N) + di, dj = offsets[rand(rng, 1:8)] + i_src = clamp(i_tgt + di, 1, N) + j_src = clamp(j_tgt + dj, 1, N) + + cid_src = state.σ[i_src, j_src] + cid_tgt = state.σ[i_tgt, j_tgt] + cid_src == cid_tgt && continue + (cid_src == -1 || cid_tgt == -1) && continue + + if cid_tgt > 0 && get(state.cell_vol, cid_tgt, 0) <= 1 + continue + end + + dH_vol = 0.0 + if cid_src > 0 + sp_s = state.cell_type[cid_src] + v_s = state.cell_vol[cid_src] + vt_s = p.V_target[sp_s + 1] + dH_vol += p.λ_vol * ((v_s + 1 - vt_s)^2 - (v_s - vt_s)^2) + end + if cid_tgt > 0 + sp_t = state.cell_type[cid_tgt] + v_t = state.cell_vol[cid_tgt] + vt_t = p.V_target[sp_t + 1] + dH_vol += p.λ_vol * ((v_t - 1 - vt_t)^2 - (v_t - vt_t)^2) + end + + dH = dH_vol + dH_radiation(state, p, i_src, j_src, i_tgt, j_tgt) + if dH <= 0.0 || rand(rng) < exp(-dH / p.T) + state.σ[i_tgt, j_tgt] = cid_src + if cid_src > 0 + state.cell_vol[cid_src] = state.cell_vol[cid_src] + 1 + end + if cid_tgt > 0 + state.cell_vol[cid_tgt] = state.cell_vol[cid_tgt] - 1 + end + end + end + update_radiation_field!(state, p) + return nothing +end + +""" + step_fields!(state::State, p::Params; dt::Float64 = 0.2) + +Advance precursor damage `M_damage` (Iodine-135 analogue) and radiolytic inhibitor +`C_inhib` (Xenon-135 analogue) with bilinear flux-melanin burnout sink. + +State & Inventory Accounting Contract: +- The 5-point stencil with reflective Neumann indices (`max(1, i-1)`, `min(N, i+1)`) + has zero net boundary flux (`sum(lap) == 0.0` to machine precision), so diffusion + strictly redistributes `C_inhib` without creating or destroying inventory. +- Before non-negativity clamping (`C_new >= 0`), the domain-integrated inhibitor + inventory satisfies exact step balance: + `sum(C_new) - sum(C_old) == dt * sum(prod - decay - burn)`. +""" +function step_fields!(state::State, p::Params; dt::Float64 = 0.2) + N = p.N + C_new = similar(state.C_inhib) + + for j in 1:N, i in 1:N + cid = state.σ[i, j] + sp = cid > 0 ? state.cell_type[cid] : 0 + ρ_mel = (sp == 4 || sp == 5) ? 1.0 : 0.0 + + state.M_damage[i, j] = (1.0 - 0.05 * dt) * state.M_damage[i, j] + + dt * 0.12 * state.radiation[i, j] + + im1 = max(1, i - 1) + ip1 = min(N, i + 1) + jm1 = max(1, j - 1) + jp1 = min(N, j + 1) + lap = state.C_inhib[im1, j] + state.C_inhib[ip1, j] + + state.C_inhib[i, jm1] + state.C_inhib[i, jp1] - 4.0 * state.C_inhib[i, j] + + prod = p.γ_delay * state.M_damage[i, j] + decay = p.λ_decay * state.C_inhib[i, j] + burn = p.σ_burn * state.radiation[i, j] * ρ_mel * state.C_inhib[i, j] + + C_new[i, j] = max(0.0, state.C_inhib[i, j] + dt * (p.D_C * lap + prod - decay - burn)) + end + state.C_inhib .= C_new + return nothing +end + +""" + simulate_xenon_pit_hysteresis(p::Params = Params(); steps_per_phase::Int = 30, dt::Float64 = 0.2) + +Run a 4-phase Chernobyl power protocol (Event 1: 50% power equilibrium -> Event 2: 30% +step-down -> Event 3: 1% Xenon-135 pit -> Event 4/5: power step-up burnout) on a +melanized site and compare `σ_burn > 0` against `σ_burn = 0`. +""" +function simulate_xenon_pit_hysteresis(p::Params = Params(); steps_per_phase::Int = 30, dt::Float64 = 0.2) + flux_schedule = vcat( + fill(0.50 * p.I0, steps_per_phase), + fill(0.30 * p.I0, steps_per_phase), + fill(0.01 * p.I0, steps_per_phase), + fill(0.85 * p.I0, steps_per_phase), + ) + M_b = 0.0 + C_b = 0.0 + M_l = 0.0 + C_l = 0.0 + trace_bilinear = Float64[] + trace_linear = Float64[] + + for I_t in flux_schedule + M_b = (1.0 - 0.05 * dt) * M_b + dt * 0.12 * I_t + C_b = max(0.0, C_b + dt * (p.γ_delay * M_b - p.λ_decay * C_b - p.σ_burn * I_t * C_b)) + push!(trace_bilinear, C_b) + + M_l = (1.0 - 0.05 * dt) * M_l + dt * 0.12 * I_t + C_l = max(0.0, C_l + dt * (p.γ_delay * M_l - p.λ_decay * C_l)) + push!(trace_linear, C_l) + end + + c_phase1_end = trace_bilinear[steps_per_phase] + c_pit_peak = maximum(trace_bilinear[(2 * steps_per_phase + 1):(3 * steps_per_phase)]) + c_phase4_end = trace_bilinear[end] + return ( + trace_bilinear = trace_bilinear, + trace_linear = trace_linear, + c_phase1_end = c_phase1_end, + c_pit_peak = c_pit_peak, + c_phase4_end = c_phase4_end, + pit_overshoot = c_pit_peak - c_phase1_end, + burnout_clearance = c_pit_peak - c_phase4_end, + ) +end + +""" + update_radial_metrics!(state::State, p::Params) + +Append species-wise mean radial distance `(t)` and void fraction `f_void(t)`. +""" +function update_radial_metrics!(state::State, p::Params) + N = p.N + cx, cy = (N + 1) / 2.0, (N + 1) / 2.0 + r_sums = zeros(Float64, p.N_types) + counts = zeros(Int, p.N_types) + void_count = 0 + + for j in 1:N, i in 1:N + cid = state.σ[i, j] + if cid == -1 + void_count += 1 + elseif cid > 0 + sp = state.cell_type[cid] + r_sums[sp] += hypot(i - cx, j - cy) + counts[sp] += 1 + end + end + + for sp in 1:p.N_types + mean_r = counts[sp] > 0 ? r_sums[sp] / counts[sp] : NaN + push!(state.mean_r_history[sp], mean_r) + end + push!(state.void_fraction_history, void_count / (N * N)) + return nothing +end + +""" + compute_derived_layers(state::State, p::Params) -> Dict{String, Any} + +Compute the complete dictionary of 2D spatial arrays and scalar regime diagnostics +matching the Python `coupling/biofilm_openmc/viewer.py` layer contract. +""" +function compute_derived_layers(state::State, p::Params) + N = p.N + species_id = zeros(Int, N, N) + kappa_map = zeros(Float64, N, N) + burnout_sink = zeros(Float64, N, N) + void_mask = zeros(Int, N, N) + melanin_mask = zeros(Int, N, N) + + sp_counts = zeros(Int, p.N_types) + total_bio = 0 + + for j in 1:N, i in 1:N + cid = state.σ[i, j] + κ_loc, _ = local_attenuation(state, p, i, j) + kappa_map[i, j] = κ_loc + if cid == -1 + species_id[i, j] = -1 + void_mask[i, j] = 1 + elseif cid > 0 + sp = state.cell_type[cid] + species_id[i, j] = sp + sp_counts[sp] += 1 + total_bio += 1 + if sp == 4 || sp == 5 + melanin_mask[i, j] = 1 + burnout_sink[i, j] = p.σ_burn * state.radiation[i, j] * state.C_inhib[i, j] + end + end + end + + mel_frac = total_bio > 0 ? sum(melanin_mask) / total_bio : 0.0 + shannon_H = 0.0 + if total_bio > 0 + for c in sp_counts + if c > 0 + prob = c / total_bio + shannon_H -= prob * log(prob) + end + end + end + + mean_void = mean(void_mask) + mean_flux = mean(state.radiation) + baseline_flux = p.I0 * exp(-0.5 * p.κ * p.κ_water) + void_reactivity = (mean_flux - baseline_flux) / max(mean_void, 1e-6) + + return Dict{String, Any}( + "cell_id" => copy(state.σ), + "species_id" => species_id, + "local_attenuation_kappa" => kappa_map, + "fast_flux" => copy(state.I_fast), + "thermal_flux" => copy(state.I_therm), + "total_radiation" => copy(state.radiation), + "precursor_damage_M" => copy(state.M_damage), + "xenon_inhibitor_C" => copy(state.C_inhib), + "bilinear_burnout_sink" => burnout_sink, + "void_mask" => void_mask, + "melanin_shield_mask" => melanin_mask, + "melanized_fraction" => Float64(mel_frac), + "shannon_diversity_H" => Float64(shannon_H), + "void_reactivity_coeff" => Float64(void_reactivity), + ) +end + +function _matrix_to_json(mat::AbstractMatrix) + rows = String[] + for i in axes(mat, 1) + push!(rows, "[" * join(mat[i, :], ", ") * "]") + end + return "[" * join(rows, ", ") * "]" +end + +""" + export_observer_payload(state::State, p::Params, path::AbstractString) + +Write a self-contained JSON payload of raw and derived spatial layers and scalar +regime invariants for `viewer.py` and `observer.py`, including explicit flags +marking the payload as a synthetic surrogate on a non-OpenMC display grid. +""" +function export_observer_payload(state::State, p::Params, path::AbstractString) + layers = compute_derived_layers(state, p) + spatial_keys = ( + "cell_id", + "species_id", + "local_attenuation_kappa", + "fast_flux", + "thermal_flux", + "total_radiation", + "precursor_damage_M", + "xenon_inhibitor_C", + "bilinear_burnout_sink", + "void_mask", + "melanin_shield_mask", + ) + open(path, "w") do io + write(io, "{\n") + write(io, " \"N\": $(p.N),\n") + write(io, " \"synthetic_surrogate\": true,\n") + write(io, " \"measured_physical_dose\": false,\n") + write(io, " \"native_openmc_field\": false,\n") + write(io, " \"authoritative_for_quantitation\": false,\n") + for k in spatial_keys + write(io, " \"$k\": $(_matrix_to_json(layers[k])),\n") + end + write(io, " \"melanized_fraction\": $(layers["melanized_fraction"]),\n") + write(io, " \"shannon_diversity_H\": $(layers["shannon_diversity_H"]),\n") + write(io, " \"void_reactivity_coeff\": $(layers["void_reactivity_coeff"])\n") + write(io, "}\n") + end + return path +end + +""" + run_control_audit() -> NamedTuple + +Run negative and positive physical controls for the two-way radiation CPM and +Xenon-135 bilinear burnout model: +1. Negative control (`κ = 0`, `β_ion = 0`): uniform radiation (`Var(I) == 0`) and zero radial bias, even when `Σ_scatter` is non-zero. +2. Positive control (`κ > 0`, `β_ion < 0` for Sp4/Sp5): non-zero attenuation gradient and inward migration (`Δ < 0`) of melanized species. +3. Two-group attenuation cascade conservation and normalization (`I_tot == I_fast + I_therm`, `0 <= I_fast <= I0`, `0 <= I_therm`, and `0 <= I_tot <= I0` under conservative scattering). +4. Void reactivity feedback (`κ_void << κ_water`): introducing voids increases mean interior flux (`α_V > 0`). +5. Bilinear burnout hysteresis and discrete zero-flux Neumann inventory balance (`sum(lap) == 0.0` and exact step accounting). +6. Exact volume registry conservation (`cell_vol[cid] == count(==(cid), σ)` for all cells). +""" +function run_control_audit() + # 1a. Negative Control via phase κ = 0 (without manually zeroing Σ_scatter) + p_neg = Params( + N = 24, + I0 = 1.0, + κ_water = 0.0, + κ_void = 0.0, + κ_species = zeros(5), + β_ion = zeros(5), + ) + st_neg = init_state(p_neg; seed = 101, void_fraction = 0.0) + update_radiation_field!(st_neg, p_neg) + var_I_neg = var(st_neg.radiation) + @assert isapprox(var_I_neg, 0.0; atol = 1e-12) "Negative control failed: radiation variance must be 0 when κ=0" + @assert all(st_neg.I_fast .== p_neg.I0) "Negative control failed: I_fast must equal I0 everywhere when κ=0" + @assert all(st_neg.radiation .== p_neg.I0) "Negative control failed: total radiation must equal I0 everywhere when κ=0" + + # 1b. Negative Control via scalar κ = 0.0 multiplier + p_neg_scalar = Params(N = 24, I0 = 1.0, κ = 0.0) + st_neg_scalar = init_state(p_neg_scalar; seed = 101, void_fraction = 0.08) + @assert all(st_neg_scalar.radiation .== p_neg_scalar.I0) "Scalar κ=0 negative control failed" + + # 2. Positive Control: state-dependent attenuation + inward tropism of Sp4/Sp5 + p_pos = Params( + N = 24, + I0 = 2.0, + T = 1.5, + λ_vol = 2.0, + V_target = Int[0, 4, 4, 4, 5, 5], + κ_water = 0.85, + κ_species = Float64[0.85, 0.95, 1.05, 2.20, 3.10], + Σ_water = 0.05, + Σ_scatter = Float64[0.05, 0.05, 0.05, 0.10, 0.12], + β_ion = Float64[0.0, 0.0, -1.0, -45.0, -65.0], + ) + st_pos = init_state(p_pos; seed = 101, void_fraction = 0.04) + r4_init = st_pos.mean_r_history[4][end] + r5_init = st_pos.mean_r_history[5][end] + rng = MersenneTwister(2026) + for _ in 1:60 + step_cpm!(st_pos, p_pos; n_trials = 24 * 24, rng = rng) + step_fields!(st_pos, p_pos; dt = 0.2) + end + update_radial_metrics!(st_pos, p_pos) + r4_final = st_pos.mean_r_history[4][end] + r5_final = st_pos.mean_r_history[5][end] + var_I_pos = var(st_pos.radiation) + @assert var_I_pos > 1e-4 "Positive control failed: expected non-zero radiation variance" + @assert (r4_final < r4_init && r5_final < r5_init) "Positive control failed: melanized species must migrate inward down β_ion < 0" + + # 2b. Two-group attenuation cascade conservation and normalization check + @assert all(0.0 .<= st_pos.I_fast .<= p_pos.I0) "Fast flux must satisfy 0 <= I_fast <= I0" + @assert all(st_pos.I_therm .>= 0.0) "Thermal flux must be non-negative" + @assert all(isapprox.(st_pos.radiation, st_pos.I_fast .+ st_pos.I_therm; atol = 1e-14)) "Total radiation must equal I_fast + I_therm" + @assert all(0.0 .<= st_pos.radiation .<= p_pos.I0 * (1.0 + 1e-12)) "Conservative two-group cascade must satisfy I_fast + I_therm <= I0" + + # 2c. Non-multiplying down-scattering guard and lossless pure-scattering conservation (Σ_eff == κ_eff, κ_therm == 0 -> I_tot == I0) + p_clamp = Params( + N = 16, + I0 = 1.5, + κ_water = 0.50, + Σ_water = 5.00, + κ_void = 0.10, + Σ_void = 1.00, + κ_species = fill(0.60, 5), + Σ_scatter = fill(4.00, 5), + κ_therm = 0.0, + ) + st_clamp = init_state(p_clamp; seed = 77, void_fraction = 0.05) + @assert all(0.0 .<= st_clamp.radiation .<= p_clamp.I0 * (1.0 + 1e-12)) "Down-scattering clamp must prevent unphysical flux multiplication above I0" + @assert all(isapprox.(st_clamp.radiation, p_clamp.I0; atol = 1e-13)) "Lossless pure scattering (Σ_eff == κ_eff, κ_therm == 0) must conserve total flux I_fast + I_therm == I0 everywhere" + + # 3. Void Reactivity Feedback Control: replacing water/biomass with voids increases interior flux + st_dense = init_state(p_pos; seed = 303, void_fraction = 0.0) + st_voidy = init_state(p_pos; seed = 303, void_fraction = 0.35) + mean_I_dense = mean(st_dense.radiation) + mean_I_voidy = mean(st_voidy.radiation) + @assert mean_I_voidy > mean_I_dense "Void feedback control failed: higher void fraction must increase mean flux" + + # 4. Xenon-135 Bilinear Burnout Hysteresis & Discrete State Accounting Control + hyst = simulate_xenon_pit_hysteresis(p_pos; steps_per_phase = 30, dt = 0.2) + @assert hyst.pit_overshoot > 0.0 "Xenon pit control failed: inhibitor C must accumulate during power step-down" + @assert hyst.burnout_clearance > 0.0 "Xenon burnout control failed: high flux must burn out inhibitor C" + @assert hyst.trace_bilinear[end] < hyst.trace_linear[end] "Identifiability control failed: bilinear burnout must clear C faster than linear decay alone" + + # Verify discrete zero-flux Neumann Laplacian conservation and exact step inventory balance + N_acc = p_pos.N + C_before = copy(st_pos.C_inhib) + M_next = (1.0 - 0.05 * 0.2) .* st_pos.M_damage .+ (0.2 * 0.12) .* st_pos.radiation + lap_sum = 0.0 + net_source_sum = 0.0 + for j in 1:N_acc, i in 1:N_acc + im1 = max(1, i - 1) + ip1 = min(N_acc, i + 1) + jm1 = max(1, j - 1) + jp1 = min(N_acc, j + 1) + lap_ij = C_before[im1, j] + C_before[ip1, j] + C_before[i, jm1] + C_before[i, jp1] - 4.0 * C_before[i, j] + lap_sum += lap_ij + cid = st_pos.σ[i, j] + sp = cid > 0 ? st_pos.cell_type[cid] : 0 + ρ_mel = (sp == 4 || sp == 5) ? 1.0 : 0.0 + prod_ij = p_pos.γ_delay * M_next[i, j] + decay_ij = p_pos.λ_decay * C_before[i, j] + burn_ij = p_pos.σ_burn * st_pos.radiation[i, j] * ρ_mel * C_before[i, j] + net_source_sum += (prod_ij - decay_ij - burn_ij) + end + @assert abs(lap_sum) < 1e-12 "Neumann discrete Laplacian must conserve inventory (sum(lap) == 0)" + step_fields!(st_pos, p_pos; dt = 0.2) + @assert isapprox(sum(st_pos.C_inhib) - sum(C_before), 0.2 * net_source_sum; atol = 1e-12) "Inhibitor step inventory change must equal dt * sum(prod - decay - burn)" + + # 5. Volume Registry Conservation Check + for (cid, v_reg) in st_pos.cell_vol + v_actual = count(==(cid), st_pos.σ) + @assert v_reg == v_actual "Volume registry mismatch for cell $cid: registered=$v_reg vs actual=$v_actual" + end + + layers = compute_derived_layers(st_pos, p_pos) + @assert length(layers) == 14 "Derived layer dictionary must contain all 14 contracted entries" + @assert haskey(layers, "melanized_fraction") "Derived layer dictionary must use melanized_fraction key" + + # 5b. Species-ablation safety (N_types < 5 must not throw BoundsError) + p_ablate = Params( + N = 12, + N_types = 3, + V_target = Int[0, 4, 4, 4], + κ_species = Float64[0.85, 0.95, 1.05], + Σ_scatter = Float64[0.05, 0.05, 0.05], + β_ion = Float64[0.0, 0.0, -1.0], + ) + st_ablate = init_state(p_ablate; seed = 19, void_fraction = 0.05) + layers_ablate = compute_derived_layers(st_ablate, p_ablate) + @assert layers_ablate["melanized_fraction"] == 0.0 "N_types=3 ablation must report melanized_fraction == 0.0 without BoundsError" + + mktempdir() do tmp + payload_path = joinpath(tmp, "observer_payload.json") + export_observer_payload(st_pos, p_pos, payload_path) + payload_txt = read(payload_path, String) + for req_key in ("cell_id", "species_id", "local_attenuation_kappa", "bilinear_burnout_sink") + @assert occursin("\"$req_key\": [[", payload_txt) "export_observer_payload missing 2D array for $req_key" + end + @assert occursin("\"melanized_fraction\":", payload_txt) "export_observer_payload missing melanized_fraction" + @assert occursin("\"synthetic_surrogate\": true", payload_txt) "export_observer_payload missing synthetic_surrogate flag" + @assert occursin("\"measured_physical_dose\": false", payload_txt) "export_observer_payload missing measured_physical_dose flag" + end + + return ( + var_I_neg = var_I_neg, + var_I_pos = var_I_pos, + delta_r_sp4 = r4_final - r4_init, + delta_r_sp5 = r5_final - r5_init, + mean_I_dense = mean_I_dense, + mean_I_voidy = mean_I_voidy, + pit_overshoot = hyst.pit_overshoot, + burnout_clearance = hyst.burnout_clearance, + ) +end + +end # module BiofilmsPottsRad + # ============================================================ # 13. Figure export (CairoMakie — headless, publication-ready) # diff --git a/coupling/biofilm_openmc/observer.py b/coupling/biofilm_openmc/observer.py index 53cacc5..b58898a 100644 --- a/coupling/biofilm_openmc/observer.py +++ b/coupling/biofilm_openmc/observer.py @@ -31,8 +31,9 @@ import numpy as np -from .viewer import (BOOLEAN, CATEGORICAL, OUT_OF_DOMAIN, UNDECLARED, - is_undeclared, read_layer, read_manifest) +from .viewer import (BOOLEAN, CATEGORICAL, LAYER_CONTRACTS, MELANIZED_FRAC_KEY, + OUT_OF_DOMAIN, UNDECLARED, is_undeclared, read_layer, + read_manifest) # Same palette and labels as the serial model's figure section, carried across # from the CairoMakie prototype on `feat/visualize-3d`. That branch is otherwise @@ -518,3 +519,308 @@ def plot_panels(panels, layer_name, *, plotter=None, link_views=True): if link_views and len(panels) > 1: p.link_views() return p + + +SPECIES_ANNOTATIONS: dict[int, str] = { + -1: "Radiolytic / Steam Void (Low Attenuation)", + 0: "Aqueous Coolant / EPS Medium (Moderator)", + 1: "Sp1: Basal Heterotroph", + 2: "Sp2: Structural EPS Producer", + 3: "Sp3: Mildly Tolerant Commensal", + 4: "Sp4: Melanin+ Tropism Core", + 5: "Sp5: Melanin Hyper-Shield (Fast-Flux Scavenger)", +} + +SPECIES_COLORS: dict[int, str] = { + -1: "#111111", + 0: "#1f3b4d", + 1: "#4daf4a", + 2: "#377eb8", + 3: "#984ea3", + 4: "#ff7f00", + 5: "#e41a1c", +} + + +class BiofilmOpenMCObserver: + """Read-only PyVista observer consuming pre-validated layers from viewer.py. + + Scientific & Software Contract: + - Strictly read-only (`READ_ONLY_OBSERVER = True`): never modifies input arrays, + never imports `feedback_uq` or `synthetic_gate`, and cannot emit quantitative + campaign verdicts (`AUTHORITATIVE_FOR_QUANTITATION = False`). + - Operates on a resampled/extruded CPM display grid (`NATIVE_OPENMC_FIELD = False`, + `GRID_REPRESENTATION = "resampled_extruded_display_grid"`). + - All radiation/flux fields are dimensionless relative synthetic surrogates + (`SYNTHETIC_SURROGATE_ONLY = True`, `MEASURED_PHYSICAL_DOSE = False`) and + must never be interpreted as measured physical dose in Gy or Gy/s. + """ + + READ_ONLY_OBSERVER: bool = True + AUTHORITATIVE_FOR_QUANTITATION: bool = False + NATIVE_OPENMC_FIELD: bool = False + SYNTHETIC_SURROGATE_ONLY: bool = True + MEASURED_PHYSICAL_DOSE: bool = False + GRID_REPRESENTATION: str = "resampled_extruded_display_grid" + + def __init__( + self, + validated_layers: dict | None = None, + *, + spacing: tuple[float, float, float] = (1.0, 1.0, 1.0), + ) -> None: + self.spacing = spacing + self.last_telemetry: dict[str, float] = {} + if validated_layers is None: + self.layers = None + return + if ( + isinstance(validated_layers, tuple) + and len(validated_layers) == 3 + and all(isinstance(v, (int, float)) for v in validated_layers) + ): + self.spacing = ( + float(validated_layers[0]), + float(validated_layers[1]), + float(validated_layers[2]), + ) + self.layers = None + return + if not isinstance(validated_layers, dict): + raise TypeError( + "BiofilmOpenMCObserver requires a validated layer dictionary" + ) + missing = [k for k in LAYER_CONTRACTS if k not in validated_layers] + if missing: + raise KeyError( + f"Observer received unvalidated payload; missing: {missing}" + ) + self.layers = validated_layers + + def build_voxel_grid( + self, layers: dict | None = None, nz_extrude: int = 1 + ): + """Construct a PyVista ImageData / UniformGrid with all spatial layers + attached as cell data, optionally extruding 2D lattices to 3D thickness. + """ + import pyvista as pv + + active = layers if layers is not None else self.layers + if not isinstance(active, dict) or "cell_id" not in active: + raise KeyError("No validated layers available to build voxel grid") + + sigma = active["cell_id"] + nx, ny = sigma.shape[:2] + if sigma.ndim == 3: + nz = sigma.shape[2] + promote_3d = False + else: + nz = max(1, int(nz_extrude)) + promote_3d = nz > 1 + + grid_cls = getattr(pv, "ImageData", None) or getattr(pv, "UniformGrid") + grid = grid_cls( + dimensions=(nx + 1, ny + 1, nz + 1), + spacing=self.spacing, + origin=(0.0, 0.0, 0.0), + ) + + for key, val in active.items(): + if isinstance(val, np.ndarray) and val.shape[:2] == (nx, ny): + vol = ( + np.repeat(val[:, :, np.newaxis], nz, axis=2) + if promote_3d + else val + ) + grid.cell_data[key] = vol.flatten(order="F") + + return grid + + @staticmethod + def compute_telemetry(layers: dict) -> dict[str, float]: + """Compute read-only display summary telemetry on the synthetic surrogate + layers (`mean_relative_flux` / `mean_dose_rate` is dimensionless relative + synthetic flux `I_tot`, never measured physical dose in Gy/s and never + authoritative for quantitative campaign verdicts). + """ + sigma = np.asarray(layers["cell_id"]) + occupied = sigma > 0 + void_bool = ( + np.asarray(layers["void_mask"], dtype=bool) + if "void_mask" in layers + else ~occupied + ) + rad = np.asarray(layers["total_radiation"], dtype=np.float64) + kappa = np.asarray(layers["local_attenuation_kappa"], dtype=np.float64) + poison = np.asarray(layers["xenon_inhibitor_C"], dtype=np.float64) + burn = np.asarray(layers["bilinear_burnout_sink"], dtype=np.float64) + mean_flux = float(np.mean(rad)) + return { + "mean_relative_flux": mean_flux, + "mean_dose_rate": mean_flux, + "void_fraction": float(np.mean(void_bool)), + "mean_kappa_occupied": ( + float(np.mean(kappa[occupied])) if np.any(occupied) else 0.0 + ), + "mean_kappa_void": ( + float(np.mean(kappa[void_bool])) if np.any(void_bool) else 0.0 + ), + "poison_inventory": float(np.sum(poison)), + "mean_burnout_rate": float(np.mean(burn)), + } + + def render_four_panel_audit( + self, + raw_state: dict | None = None, + params: dict | None = None, + screenshot_path: str | None = None, + off_screen: bool = True, + plotter=None, + ): + """Render a 4-panel linked-view PyVista audit of: + (0,0) Thresholded Categorical Species & Void Architecture (crinkle-clipped) + (0,1) Clipped State-Dependent Attenuation Coefficient kappa(sigma) + (1,0) Clipped Two-Group Ionizing Flux I_fast + I_therm + (1,1) Clipped Delayed Inhibitor C & Bilinear Burnout Contours + """ + import pyvista as pv + + from .viewer import compute_and_validate_layers + + return_telemetry = raw_state is not None + if raw_state is not None: + active_layers = compute_and_validate_layers( + raw_state, params or {} + ) + self.layers = active_layers + nz_ext = 8 + else: + if self.layers is None: + raise ValueError("No layers provided to render_four_panel_audit") + active_layers = self.layers + nz_ext = 1 + + grid = self.build_voxel_grid(active_layers, nz_extrude=nz_ext) + pl = plotter or pv.Plotter( + shape=(2, 2), off_screen=off_screen, window_size=(1600, 1200) + ) + + mel_frac = float(active_layers[MELANIZED_FRAC_KEY]) + shannon_H = float(active_layers["shannon_diversity_H"]) + void_coeff = float(active_layers["void_reactivity_coeff"]) + + # Panel (0, 0): Thresholded & Crinkle-Clipped Categorical Species Parcels + pl.subplot(0, 0) + pl.add_title("1. Species & Void Architecture", font_size=11) + if hasattr(grid, "threshold"): + biomass_mesh = grid.threshold(value=0.5, scalars="species_id") + if getattr(biomass_mesh, "n_cells", 0) > 0: + if hasattr(biomass_mesh, "clip"): + biomass_mesh = biomass_mesh.clip(normal="-y", crinkle=True) + pl.add_mesh( + biomass_mesh, + scalars="species_id", + categories=True, + cmap="tab10", + clim=(1, 5), + annotations=SPECIES_ANNOTATIONS, + show_edges=True, + scalar_bar_args={"title": "Species ID (Thresholded > 0)"}, + ) + else: + pl.add_mesh( + grid, + scalars="species_id", + categories=True, + cmap="tab10", + clim=[-1, 5], + annotations=SPECIES_ANNOTATIONS, + show_edges=False, + scalar_bar_args={"title": "Species / Void ID"}, + ) + legend_entries = [ + [SPECIES_ANNOTATIONS[k], SPECIES_COLORS[k]] + for k in [-1, 0, 1, 2, 3, 4, 5] + ] + pl.add_legend( + legend_entries, bcolor=(0.1, 0.1, 0.1), size=(0.38, 0.28) + ) + pl.add_text( + f"Melanin+ Frac: {mel_frac:.2%} | H: {shannon_H:.2f}", + position="lower_left", + font_size=9, + ) + + # Panel (0, 1): State-Dependent Attenuation clipped across mid-plane + pl.subplot(0, 1) + pl.add_title("2. State-Dependent Attenuation kappa(sigma)", font_size=11) + att_mesh = ( + grid.clip(normal="-y", crinkle=False) + if hasattr(grid, "clip") + else grid + ) + pl.add_mesh( + att_mesh, + scalars="local_attenuation_kappa", + cmap="magma", + scalar_bar_args={"title": "kappa [1/lattice]"}, + ) + pl.add_text( + f"Void Reactivity alpha_V: {void_coeff:+.3f}", + position="lower_left", + font_size=9, + ) + + # Panel (1, 0): Two-Group Radiation Field clipped across mid-plane + pl.subplot(1, 0) + pl.add_title("3. Two-Group Radiation Field (Fast + Thermal)", font_size=11) + rad_mesh = ( + grid.clip(normal="-y", crinkle=False) + if hasattr(grid, "clip") + else grid + ) + pl.add_mesh( + rad_mesh, + scalars="total_radiation", + cmap="inferno", + scalar_bar_args={"title": "Flux I_tot [rel, synthetic surrogate - NOT Gy]"}, + ) + + # Panel (1, 1): Delayed Inhibitor C & Bilinear Burnout Sink + pl.subplot(1, 1) + pl.add_title("4. Xenon-Style Inhibitor C & Burnout Sink", font_size=11) + inh_mesh = ( + grid.clip(normal="-y", crinkle=False) + if hasattr(grid, "clip") + else grid + ) + pl.add_mesh( + inh_mesh, + scalars="xenon_inhibitor_C", + cmap="viridis", + scalar_bar_args={"title": "Inhibitor C [a.u.]"}, + ) + if float(np.max(active_layers["bilinear_burnout_sink"])) > 1e-6: + point_grid = grid.cell_data_to_point_data() + contours = point_grid.contour( + isosurfaces=4, scalars="bilinear_burnout_sink" + ) + if contours.n_points > 0: + pl.add_mesh(contours, color="cyan", line_width=2.5) + + pl.link_views() + if screenshot_path: + pl.screenshot(screenshot_path) + + telemetry = self.compute_telemetry(active_layers) + self.last_telemetry = telemetry + try: + pl.telemetry = telemetry + except Exception: + pass + + if return_telemetry: + if hasattr(pl, "close"): + pl.close() + return telemetry + return pl diff --git a/coupling/biofilm_openmc/snapshot.py b/coupling/biofilm_openmc/snapshot.py index 25882a4..487bf77 100644 --- a/coupling/biofilm_openmc/snapshot.py +++ b/coupling/biofilm_openmc/snapshot.py @@ -16,7 +16,10 @@ from dataclasses import dataclass -import h5py +try: + import h5py +except ModuleNotFoundError: # pragma: no cover + h5py = None import numpy as np diff --git a/coupling/biofilm_openmc/viewer.py b/coupling/biofilm_openmc/viewer.py index 08a23cd..39571bb 100644 --- a/coupling/biofilm_openmc/viewer.py +++ b/coupling/biofilm_openmc/viewer.py @@ -519,3 +519,262 @@ def quantitative_layers(manifest_doc) -> list[str]: """ return [l["name"] for l in manifest_doc["layers"] if l["authoritative_for_quantitation"]] + + +MELANIZED_FRAC_KEY = "melanized_fraction" + + +@dataclass(frozen=True) +class LayerSpec: + name: str + units: str + dtype: str = "float64" + quotable: bool = True + description: str = "" + cmap: str = "viridis" + clim: tuple[float, float] = (0.0, 1.0) + categorical: bool = False + derived: bool = False + native_openmc_field: bool = False + grid_representation: str = "resampled_extruded_display_grid" + synthetic_surrogate: bool = True + measured_physical_dose: bool = False + authoritative_for_quantitation: bool = False + + +LAYER_CONTRACTS: dict[str, LayerSpec] = { + "cell_id": LayerSpec( + "cell_id", "id", "int32", True, + "Discrete CPM spin ID (0=water, -1=void/steam)", + cmap="tab20", clim=(-1.0, 64.0), categorical=True, derived=False, + ), + "species_id": LayerSpec( + "species_id", "index", "int32", True, + "Species index 1..5 (0=water, -1=void)", + cmap="tab10", clim=(1.0, 7.0), categorical=True, derived=False, + ), + "local_attenuation_kappa": LayerSpec( + "local_attenuation_kappa", "1/lattice", "float64", True, + "State-dependent macro attenuation coefficient kappa(sigma)", + cmap="viridis", clim=(0.0, 6.0), categorical=False, derived=True, + ), + "fast_flux": LayerSpec( + "fast_flux", "rel_flux", "float64", True, + "Primary unmoderated fast radiation group I_fast (synthetic surrogate, not measured Gy)", + cmap="inferno", clim=(0.0, 1.0), categorical=False, derived=False, + ), + "thermal_flux": LayerSpec( + "thermal_flux", "rel_flux", "float64", True, + "Moderated secondary radiolytic group I_therm (synthetic surrogate, not measured Gy)", + cmap="inferno", clim=(0.0, 1.0), categorical=False, derived=False, + ), + "total_radiation": LayerSpec( + "total_radiation", "rel_flux", "float64", True, + "Total effective ionizing field I_fast + I_therm (synthetic surrogate, not measured Gy)", + cmap="inferno", clim=(0.0, 1.0), categorical=False, derived=True, + ), + "precursor_damage_M": LayerSpec( + "precursor_damage_M", "conc_au", "float64", True, + "Accumulated precursor damage (Iodine-135 analogue)", + cmap="magma", clim=(0.0, 1.5), categorical=False, derived=False, + ), + "xenon_inhibitor_C": LayerSpec( + "xenon_inhibitor_C", "conc_au", "float64", True, + "Delayed radiolytic inhibitor (Xenon-135 poison analogue)", + cmap="magma", clim=(0.0, 1.5), categorical=False, derived=True, + ), + "bilinear_burnout_sink": LayerSpec( + "bilinear_burnout_sink", "conc_au/s", "float64", True, + "Active flux-scavenging burnout rate sigma_burn * I * rho_mel * C", + cmap="plasma", clim=(0.0, 0.8), categorical=False, derived=True, + ), + "void_mask": LayerSpec( + "void_mask", "bool", "uint8", False, + "Binary mask of low-attenuation radiolytic/steam voids", + cmap="gray", clim=(0.0, 1.0), categorical=True, derived=True, + ), + "melanin_shield_mask": LayerSpec( + "melanin_shield_mask", "bool", "uint8", False, + "Binary mask of high-attenuation melanized cells (Sp 4 & 5)", + cmap="inferno", clim=(0.0, 1.0), categorical=True, derived=True, + ), + MELANIZED_FRAC_KEY: LayerSpec( + MELANIZED_FRAC_KEY, "fraction", "float64", True, + "Biomass fraction of melanized species (Sp 4 + Sp 5)", + cmap="viridis", clim=(0.0, 1.0), categorical=False, derived=True, + ), + "shannon_diversity_H": LayerSpec( + "shannon_diversity_H", "nats", "float64", True, + "Community Shannon diversity index H", + cmap="viridis", clim=(0.0, 2.0), categorical=False, derived=True, + ), + "void_reactivity_coeff": LayerSpec( + "void_reactivity_coeff", "dI/df_void", "float64", True, + "Effective void amplification diagnostic (positive void coefficient)", + cmap="coolwarm", clim=(-1.0, 1.0), categorical=False, derived=True, + ), +} + + +def compute_and_validate_layers(raw_state: dict, params: dict) -> dict: + """Compute all derived spatial and scalar layers in viewer.py and verify + that every LAYER_CONTRACTS entry is present and finite. + Supports both the two-group state dictionary ('cell_id', 'species_id', + 'fast_flux', 'thermal_flux', 'precursor_damage_M', 'xenon_inhibitor_C') + and the continuum-melanin snapshot dictionary ('spin', 'species', + 'melanin', 'radiation', 'poison'). + + Scientific & Software Contract: + - Never mutates `raw_state` or `params`. + - All radiation/flux layers (`fast_flux`, `thermal_flux`, `total_radiation`) + are dimensionless synthetic surrogates (`rel_flux`), never measured + physical dose (`Gy` or `Gy/s`), and are defined on the CPM display grid + rather than a native OpenMC transport tally mesh. + """ + continuum_mode = "spin" in raw_state and "cell_id" not in raw_state + + if continuum_mode: + sigma = np.array(raw_state["spin"], dtype=np.int32, copy=True) + occupied = sigma > 0 + raw_sp = np.asarray(raw_state["species"], dtype=np.int32) + if raw_sp.ndim == 1: + species = np.zeros_like(sigma, dtype=np.int32) + species[occupied] = raw_sp[sigma[occupied] - 1] + else: + species = np.array(raw_sp, dtype=np.int32, copy=True) + melanin = np.array(raw_state["melanin"], dtype=np.float64, copy=True) + I_tot = np.array(raw_state["radiation"], dtype=np.float64, copy=True) + I_fast = I_tot.copy() + I_therm = np.zeros_like(I_tot) + C_inh = np.array(raw_state["poison"], dtype=np.float64, copy=True) + M_dam = np.array( + raw_state.get("poison_lag", melanin), dtype=np.float64, copy=True + ) + if not ( + sigma.shape + == species.shape + == melanin.shape + == I_tot.shape + == C_inh.shape + == M_dam.shape + ): + raise ValueError("Spatial input arrays must share identical dimensions") + + kappa_base = float(params.get("kappa", 1.5)) + kappa_void = float(params.get("kappa_void", 0.15)) + kappa_cell = float(params.get("kappa_cell", 1.0)) + kappa_mel = float(params.get("kappa_mel", 2.5)) + if kappa_base == 0.0: + kappa_map = np.zeros_like(I_tot) + else: + kappa_map = np.where( + occupied, + kappa_base * (kappa_cell + kappa_mel * melanin), + kappa_base * kappa_void, + ) + sigma_burn = float(params.get("sigma_burn", 0.65)) + burnout_sink = sigma_burn * I_tot * (0.2 + melanin) * C_inh + void_mask = (~occupied).astype(np.uint8) + melanin_mask = ((species >= 4) | (melanin >= 0.4)).astype(np.uint8) + kappa_water = kappa_base * kappa_cell + else: + sigma = np.array(raw_state["cell_id"], dtype=np.int32, copy=True) + species = np.array(raw_state["species_id"], dtype=np.int32, copy=True) + I_fast = np.array(raw_state["fast_flux"], dtype=np.float64, copy=True) + I_therm = np.array(raw_state["thermal_flux"], dtype=np.float64, copy=True) + M_dam = np.array(raw_state["precursor_damage_M"], dtype=np.float64, copy=True) + C_inh = np.array(raw_state["xenon_inhibitor_C"], dtype=np.float64, copy=True) + + if not ( + sigma.shape + == species.shape + == I_fast.shape + == I_therm.shape + == M_dam.shape + == C_inh.shape + ): + raise ValueError("Spatial input arrays must share identical dimensions") + + I_tot = I_fast + I_therm + void_mask = (sigma == -1).astype(np.uint8) + melanin_mask = np.isin(species, [4, 5]).astype(np.uint8) + + kappa_mult = float(params.get("kappa", 1.0)) + kappa_water = kappa_mult * float(params.get("kappa_water", 0.35)) + kappa_void = kappa_mult * float(params.get("kappa_void", 0.02)) + kappa_sp = list(params.get("kappa_species", [0.35, 0.45, 0.50, 1.40, 2.20])) + + if kappa_mult == 0.0: + kappa_map = np.zeros(sigma.shape, dtype=np.float64) + else: + kappa_map = np.full(sigma.shape, kappa_water, dtype=np.float64) + kappa_map[sigma == -1] = kappa_void + for sp_idx, k_val in enumerate(kappa_sp, start=1): + kappa_map[species == sp_idx] = kappa_mult * float(k_val) + + sigma_burn = float(params.get("sigma_burn", 0.45)) + burnout_sink = ( + sigma_burn * I_tot * melanin_mask.astype(np.float64) * C_inh + ) + + bio_mask = species > 0 + total_bio = int(np.sum(bio_mask)) + if total_bio > 0: + mel_frac = float(np.sum(melanin_mask[bio_mask]) / total_bio) + max_sp = max(5, int(np.max(species))) + counts = np.array( + [np.sum(species == s) for s in range(1, max_sp + 1)], + dtype=np.float64, + ) + probs = counts[counts > 0] / total_bio + shannon_H = float(-np.sum(probs * np.log(probs))) + else: + mel_frac = 0.0 + shannon_H = 0.0 + + mean_void = float(np.mean(void_mask)) if void_mask.size > 0 else 0.0 + mean_flux = float(np.mean(I_tot)) if I_tot.size > 0 else 0.0 + baseline_flux = float(params.get("I0", 1.0)) * float( + np.exp(-0.5 * kappa_water) + ) + void_reactivity = (mean_flux - baseline_flux) / max(mean_void, 1e-6) + + out = { + "cell_id": sigma, + "species_id": species, + "local_attenuation_kappa": kappa_map, + "fast_flux": I_fast, + "thermal_flux": I_therm, + "total_radiation": I_tot, + "precursor_damage_M": M_dam, + "xenon_inhibitor_C": C_inh, + "bilinear_burnout_sink": burnout_sink, + "void_mask": void_mask, + "melanin_shield_mask": melanin_mask, + MELANIZED_FRAC_KEY: mel_frac, + "shannon_diversity_H": shannon_H, + "void_reactivity_coeff": void_reactivity, + } + + if continuum_mode: + out.update( + { + "spin": sigma, + "melanin": melanin, + "radiation": I_tot, + "local_attenuation": kappa_map, + "poison": C_inh, + "burnout_rate": burnout_sink, + } + ) + + for key, spec in LAYER_CONTRACTS.items(): + if key not in out: + raise KeyError(f"Contract violation: missing derived layer '{key}'") + val = out[key] + if spec.quotable and not np.all(np.isfinite(val)): + raise ValueError( + f"Contract violation: quotable layer '{key}' contains NaN/Inf" + ) + + return out diff --git a/coupling/tests/test_viewer_bundle.py b/coupling/tests/test_viewer_bundle.py index eda157f..cab89b4 100644 --- a/coupling/tests/test_viewer_bundle.py +++ b/coupling/tests/test_viewer_bundle.py @@ -249,3 +249,415 @@ def test_a_correct_reduction_appears_in_the_quotable_set(): derivation="mass_weighted_mean") doc = manifest([LATTICE, REFINED], [reduced]) assert quantitative_layers(doc) == ["dose_on_lattice"] + + +def _sample_two_way_raw_state(shape=(6, 6)): + cell_id = np.zeros(shape, dtype=np.int32) + cell_id[0, 0] = -1 + cell_id[1, 1] = 1 + cell_id[2, 2] = 4 + cell_id[3, 3] = 5 + species_id = np.zeros(shape, dtype=np.int32) + species_id[0, 0] = -1 + species_id[1, 1] = 1 + species_id[2, 2] = 4 + species_id[3, 3] = 5 + fast_flux = np.full(shape, 0.8, dtype=np.float64) + thermal_flux = np.full(shape, 0.2, dtype=np.float64) + precursor_M = np.full(shape, 0.4, dtype=np.float64) + inhibitor_C = np.full(shape, 0.5, dtype=np.float64) + return { + "cell_id": cell_id, + "species_id": species_id, + "fast_flux": fast_flux, + "thermal_flux": thermal_flux, + "precursor_damage_M": precursor_M, + "xenon_inhibitor_C": inhibitor_C, + } + + +def test_compute_and_validate_layers_satisfies_all_contracts(): + from biofilm_openmc.viewer import ( + LAYER_CONTRACTS, + MELANIZED_FRAC_KEY, + compute_and_validate_layers, + ) + + raw = _sample_two_way_raw_state() + params = { + "I0": 1.0, + "kappa_water": 0.35, + "kappa_void": 0.02, + "kappa_species": [0.35, 0.45, 0.50, 1.40, 2.20], + "sigma_burn": 0.45, + } + out = compute_and_validate_layers(raw, params) + assert set(out.keys()) == set(LAYER_CONTRACTS.keys()) + assert np.isclose(out["local_attenuation_kappa"][0, 0], 0.02) + assert np.isclose(out["local_attenuation_kappa"][2, 2], 1.40) + assert np.isclose(out["local_attenuation_kappa"][3, 3], 2.20) + assert np.isclose(out["bilinear_burnout_sink"][1, 1], 0.0) + assert np.isclose(out["bilinear_burnout_sink"][2, 2], 0.45 * 1.0 * 0.5) + assert np.isclose(out[MELANIZED_FRAC_KEY], 2.0 / 3.0) + assert out["shannon_diversity_H"] > 0.0 + + +def test_compute_and_validate_layers_edge_cases_and_negative_controls(): + from biofilm_openmc.observer import BiofilmOpenMCObserver + from biofilm_openmc.viewer import ( + MELANIZED_FRAC_KEY, + compute_and_validate_layers, + ) + + # Empty biomass (all water / voids): no division by zero + raw_empty = _sample_two_way_raw_state() + raw_empty["cell_id"][:] = 0 + raw_empty["species_id"][:] = 0 + out_empty = compute_and_validate_layers(raw_empty, {}) + assert out_empty[MELANIZED_FRAC_KEY] == 0.0 + assert out_empty["shannon_diversity_H"] == 0.0 + + # Shape mismatch rejected + raw_bad_shape = _sample_two_way_raw_state() + raw_bad_shape["fast_flux"] = np.ones((4, 4), dtype=np.float64) + with pytest.raises(ValueError, match="identical dimensions"): + compute_and_validate_layers(raw_bad_shape, {}) + + # NaN in quotable layer rejected + raw_nan = _sample_two_way_raw_state() + raw_nan["xenon_inhibitor_C"][1, 1] = np.nan + with pytest.raises(ValueError, match="contains NaN/Inf"): + compute_and_validate_layers(raw_nan, {}) + + # Observer rejects unvalidated or incomplete payloads + with pytest.raises(TypeError, match="validated layer dictionary"): + BiofilmOpenMCObserver(["not_a_dict"]) + with pytest.raises(KeyError, match="missing"): + BiofilmOpenMCObserver({"cell_id": raw_empty["cell_id"]}) + obs = BiofilmOpenMCObserver(out_empty) + assert obs.layers is out_empty + + +def test_biofilm_openmc_observer_four_panel_audit_with_stub_pyvista( + monkeypatch, tmp_path +): + import sys + import types + + from biofilm_openmc.observer import BiofilmOpenMCObserver + from biofilm_openmc.viewer import compute_and_validate_layers + + class _StubContours: + n_points = 4 + + class _StubPointGrid: + def contour(self, isosurfaces=4, scalars=None): + return _StubContours() + + class _StubImageData: + def __init__(self, dimensions, spacing, origin): + self.dimensions = dimensions + self.spacing = spacing + self.origin = origin + self.cell_data = {} + + def cell_data_to_point_data(self): + return _StubPointGrid() + + class _StubPlotter: + def __init__(self, shape=(2, 2), off_screen=True, window_size=(1600, 1200)): + self.shape = shape + self.off_screen = off_screen + self.window_size = window_size + self.subplots_visited = [] + self.titles = [] + self.meshes = [] + self.legends = [] + self.texts = [] + self.views_linked = False + self.screenshots = [] + + def subplot(self, r, c): + self.subplots_visited.append((r, c)) + + def add_title(self, title, font_size=11): + self.titles.append(title) + + def add_mesh(self, mesh, **kwargs): + self.meshes.append((mesh, kwargs)) + + def add_legend(self, entries, **kwargs): + self.legends.append(entries) + + def add_text(self, text, **kwargs): + self.texts.append(text) + + def link_views(self): + self.views_linked = True + + def screenshot(self, path): + self.screenshots.append(str(path)) + + stub_pv = types.ModuleType("pyvista") + stub_pv.ImageData = _StubImageData + stub_pv.Plotter = _StubPlotter + monkeypatch.setitem(sys.modules, "pyvista", stub_pv) + + raw = _sample_two_way_raw_state() + validated = compute_and_validate_layers(raw, {}) + obs = BiofilmOpenMCObserver(validated) + shot = tmp_path / "four_panel.png" + pl = obs.render_four_panel_audit(screenshot_path=str(shot), off_screen=True) + assert pl.subplots_visited == [(0, 0), (0, 1), (1, 0), (1, 1)] + assert len(pl.titles) == 4 + assert len(pl.legends) == 1 + assert pl.views_linked is True + assert pl.screenshots == [str(shot)] + assert "mean_dose_rate" in pl.telemetry + + +def test_continuum_melanin_and_threshold_clipped_observer_telemetry( + monkeypatch, tmp_path +): + import sys + import types + + from biofilm_openmc.observer import ( + SPECIES_ANNOTATIONS, + BiofilmOpenMCObserver, + ) + from biofilm_openmc.viewer import compute_and_validate_layers + + # Exact 64x64 snapshot from the specification + N = 64 + spin = np.zeros((N, N), dtype=np.int32) + spin[26:38, 26:38] = 1 + spin[16:26, 28:36] = 2 + species_lut = np.array([4, 2], dtype=np.int32) + + melanin = np.where(spin == 1, 0.75, np.where(spin == 2, 0.15, 0.0)) + y, x = np.mgrid[0:N, 0:N] + r = np.hypot(x - 31.5, y - 31.5) + radiation = np.exp(-1.5 * r / N) + poison = 0.4 * np.exp(-((r - 12.0) / 6.0) ** 2) + + raw_state = { + "spin": spin, + "species": species_lut, + "melanin": melanin, + "radiation": radiation, + "poison": poison, + } + params = { + "kappa": 1.5, + "kappa_void": 0.15, + "kappa_cell": 1.0, + "kappa_mel": 2.5, + "sigma_burn": 0.65, + } + + layers = compute_and_validate_layers(raw_state, params) + assert np.isclose(np.mean(layers["radiation"]), 0.5766, atol=1e-4) + + # Negative control: kappa = 0.0 -> local_attenuation == 0 everywhere + layers_neg = compute_and_validate_layers(raw_state, {"kappa": 0.0}) + assert np.all(layers_neg["local_attenuation"] == 0.0) + + class _StubContours: + n_points = 4 + + class _StubPointGrid: + def contour(self, isosurfaces=4, scalars=None): + return _StubContours() + + class _StubFullMesh: + def __init__(self, dimensions=(65, 65, 9), spacing=(1.0, 1.0, 1.0), origin=(0.0, 0.0, 0.0)): + self.dimensions = dimensions + self.spacing = spacing + self.origin = origin + self.cell_data = {} + self.n_cells = 224 + self.threshold_calls = [] + self.clip_calls = [] + + def threshold(self, value=0.5, scalars=None): + self.threshold_calls.append((value, scalars)) + return self + + def clip(self, normal="-y", crinkle=False): + self.clip_calls.append((normal, crinkle)) + return self + + def cell_data_to_point_data(self): + return _StubPointGrid() + + created_plotters = [] + + class _StubPlotter: + def __init__(self, shape=(2, 2), off_screen=True, window_size=(1600, 1200)): + self.meshes = [] + self.closed = False + created_plotters.append(self) + + def subplot(self, r, c): + pass + + def add_title(self, title, font_size=11): + pass + + def add_mesh(self, mesh, **kwargs): + self.meshes.append((mesh, kwargs)) + + def add_legend(self, entries, **kwargs): + pass + + def add_text(self, text, **kwargs): + pass + + def link_views(self): + pass + + def screenshot(self, path): + pass + + def close(self): + self.closed = True + + stub_pv = types.ModuleType("pyvista") + stub_pv.ImageData = _StubFullMesh + stub_pv.Plotter = _StubPlotter + monkeypatch.setitem(sys.modules, "pyvista", stub_pv) + + obs = BiofilmOpenMCObserver(spacing=(1.0, 1.0, 1.0)) + metrics = obs.render_four_panel_audit( + raw_state, params, screenshot_path=str(tmp_path / "audit.png") + ) + assert round(metrics["mean_dose_rate"], 4) == 0.5766 + assert round(metrics["void_fraction"], 4) == 0.9453 + assert round(metrics["mean_kappa_occupied"], 4) == 3.5089 + assert round(metrics["mean_kappa_void"], 4) == 0.225 + assert round(metrics["poison_inventory"], 4) == 320.8142 + assert round(metrics["mean_burnout_rate"], 4) == 0.0094 + + pl = created_plotters[-1] + mesh0, kwargs0 = pl.meshes[0] + assert mesh0.dimensions == (65, 65, 9) + assert (0.5, "species_id") in mesh0.threshold_calls + assert ("-y", True) in mesh0.clip_calls + assert ("-y", False) in mesh0.clip_calls + assert kwargs0.get("categories") is True + assert kwargs0.get("annotations") == SPECIES_ANNOTATIONS + + +def test_two_way_attenuation_and_observer_scientific_review_contracts(): + """Verify the five scientific and software review contracts for the + two-way attenuation cascade, poison-burnout surrogate, display vs native + OpenMC grid distinction, read-only observer isolation, and synthetic + surrogate dose labeling. + """ + from biofilm_openmc.observer import BiofilmOpenMCObserver + from biofilm_openmc.viewer import ( + LAYER_CONTRACTS, + MELANIZED_FRAC_KEY, + compute_and_validate_layers, + ) + + # 1. Two-group attenuation cascade conservation and normalization + raw = _sample_two_way_raw_state(shape=(8, 8)) + raw_copy = {k: v.copy() for k, v in raw.items()} + params = { + "I0": 1.0, + "kappa_water": 0.35, + "kappa_void": 0.02, + "kappa_species": [0.35, 0.45, 0.50, 1.40, 2.20], + "sigma_burn": 0.45, + } + out = compute_and_validate_layers(raw, params) + assert np.all(out["fast_flux"] >= 0.0) and np.all(out["fast_flux"] <= params["I0"]) + assert np.all(out["thermal_flux"] >= 0.0) + assert np.allclose(out["total_radiation"], out["fast_flux"] + out["thermal_flux"]) + assert np.all(out["total_radiation"] <= params["I0"] + 1e-12) + + # 2. Energy/state accounting in the poison-burnout surrogate + mel_mask = np.isin(raw["species_id"], [4, 5]).astype(np.float64) + expected_burn = params["sigma_burn"] * out["total_radiation"] * mel_mask * raw["xenon_inhibitor_C"] + assert np.allclose(out["bilinear_burnout_sink"], expected_burn) + assert np.all(out["bilinear_burnout_sink"] >= 0.0) + + # Discrete zero-flux Neumann 5-point Laplacian sum is identically zero + C = raw["xenon_inhibitor_C"].copy() + C[2, 3] = 1.25 + C[0, 0] = 0.05 + C[7, 7] = 0.90 + nx, ny = C.shape + lap = np.zeros_like(C) + for i in range(nx): + for j in range(ny): + im1 = max(0, i - 1) + ip1 = min(nx - 1, i + 1) + jm1 = max(0, j - 1) + jp1 = min(ny - 1, j + 1) + lap[i, j] = C[im1, j] + C[ip1, j] + C[i, jm1] + C[i, jp1] - 4.0 * C[i, j] + assert np.isclose(np.sum(lap), 0.0, atol=1e-14) + + # 3. Native OpenMC fields vs resampled/extruded display fields + assert BiofilmOpenMCObserver.NATIVE_OPENMC_FIELD is False + assert BiofilmOpenMCObserver.GRID_REPRESENTATION == "resampled_extruded_display_grid" + for spec in LAYER_CONTRACTS.values(): + assert spec.native_openmc_field is False + assert spec.grid_representation == "resampled_extruded_display_grid" + assert spec.authoritative_for_quantitation is False + + # 4. Observer remains strictly read-only and non-authoritative for quantitation + assert BiofilmOpenMCObserver.READ_ONLY_OBSERVER is True + assert BiofilmOpenMCObserver.AUTHORITATIVE_FOR_QUANTITATION is False + for out_key in ( + "cell_id", + "species_id", + "fast_flux", + "thermal_flux", + "precursor_damage_M", + "xenon_inhibitor_C", + "total_radiation", + ): + out[out_key].flat[0] = 999 + for k, arr in raw.items(): + assert np.array_equal(arr, raw_copy[k]) + + # Scalar kappa=0 in two-group mode yields zero attenuation map + out_zero_k = compute_and_validate_layers(raw, {**params, "kappa": 0.0}) + assert np.all(out_zero_k["local_attenuation_kappa"] == 0.0) + + # Water (cell_id == 0) in two-group mode is not miscounted as void (cell_id == -1) + raw_no_void = {k: v.copy() for k, v in raw.items()} + raw_no_void["cell_id"][raw_no_void["cell_id"] == -1] = 0 + raw_no_void["species_id"][raw_no_void["species_id"] == -1] = 0 + layers_no_void = compute_and_validate_layers(raw_no_void, params) + tel_no_void = BiofilmOpenMCObserver.compute_telemetry(layers_no_void) + assert tel_no_void["void_fraction"] == 0.0 + assert tel_no_void["mean_relative_flux"] == tel_no_void["mean_dose_rate"] + + # Continuum-mode poison_lag shape mismatch is rejected + with pytest.raises(ValueError, match="identical dimensions"): + compute_and_validate_layers( + { + "spin": np.zeros((4, 4), dtype=np.int32), + "species": np.zeros((4, 4), dtype=np.int32), + "melanin": np.zeros((4, 4), dtype=np.float64), + "radiation": np.ones((4, 4), dtype=np.float64), + "poison": np.zeros((4, 4), dtype=np.float64), + "poison_lag": np.zeros((3, 3), dtype=np.float64), + }, + {}, + ) + + # 5. Synthetic neutron/dose outputs cannot be mistaken for measured physical dose + assert BiofilmOpenMCObserver.SYNTHETIC_SURROGATE_ONLY is True + assert BiofilmOpenMCObserver.MEASURED_PHYSICAL_DOSE is False + for flux_key in ("fast_flux", "thermal_flux", "total_radiation"): + spec = LAYER_CONTRACTS[flux_key] + assert spec.synthetic_surrogate is True + assert spec.measured_physical_dose is False + assert spec.units == "rel_flux" + assert "synthetic surrogate" in spec.description + assert MELANIZED_FRAC_KEY == "melanized_fraction" diff --git a/preprint/modeling_radioresistance_and_radiotropic_fitness.tex b/preprint/modeling_radioresistance_and_radiotropic_fitness.tex index a4442af..c0ba476 100644 --- a/preprint/modeling_radioresistance_and_radiotropic_fitness.tex +++ b/preprint/modeling_radioresistance_and_radiotropic_fitness.tex @@ -597,7 +597,7 @@ \section{Computational Methods} the Julia axis convention and joined to current cell, lineage and generation labels. The executed end-to-end validation uses an explicitly synthetic reference system with \texttt{target\_calibration=false} and stops before changing the Hamiltonian, melanin drive, -membrane state, accumulated dose, or CPM clock. +membrane state, accumulated dose, or CPM clock. While the executed OpenMC snapshot exchange is one-way, the repository's radiation coupling is no longer strictly one-way when the optional state-dependent attenuation surrogate (\texttt{BiofilmsPottsRad}, \texttt{target\_calibration=false}) is active, which recomputes two-group relative flux from evolving lattice occupancy and fixed per-species coefficients while driving a separate one-way precursor and inhibitor burnout field. \paragraph{Reproducibility and implementation contracts.} The serial CPM has a byte-pinned regression fixture. The JACC port is tested by per-kernel agreement @@ -920,7 +920,7 @@ \subsection{Synthetic One-Way CPM--OpenMC Integration} membrane state or CPM clock. This validates the implementation path, not the physical target. The ensemble of Section~\ref{sec:radial_ensemble} runs under the analytic field of Eq.~\ref{eq:gamma_field}; the transport path contributes no number to it, and the interface is -reported here as executed, with no result of its own. +reported here as executed, with no result of its own; in code, coupling is no longer strictly one-way when the opt-in state-dependent two-way attenuation and two-group cascade surrogate (\texttt{BiofilmsPottsRad}, accompanied by a one-way radiation-driven inhibitor burnout field and the 4-panel observer contract) is enabled. % ---- 7. Discussion ------------------------------------------- \section{Discussion} @@ -1713,7 +1713,7 @@ \subsection{Planned Radiation--Transport Feedback, and What It Would Assume} the CPM changes occupancy, so $\mu_{\mathrm{matrix}}$ varies over a run, and at low loading that is plausibly the larger effect --- and preferential absorption of low-energy photons hardens the local spectrum, which is \emph{beam} hardening, an ordinary transport phenomenon and not an -adaptive biological response. +adaptive biological response. In the repository code, an opt-in deterministic surrogate (\texttt{BiofilmsPottsRad}) relaxes the strictly one-way static-field assumption for occupancy by recomputing local fast-group attenuation and fast-to-thermal moderation from evolving lattice occupancy (water, void, and cell sites) and fixed per-species coefficients in a non-multiplying two-group cascade, while driving precursor damage and bilinear inhibitor burnout as a one-way radiation output rather than feeding sorbed inventory or an evolving melanin field back into $\mu(\mathbf{r})$ (\texttt{target\_calibration=false}); the Monte Carlo OpenMC transport exchange likewise remains one-way pending the physical composition measurements above. \emph{A single-parameter exponential damage scalar is a strong structural prior.} A form such as $\gamma = 1 - e^{-\lambda D}$, applied as diff --git a/scripts/preflight_merge.sh b/scripts/preflight_merge.sh index cf78530..e8fb0de 100755 --- a/scripts/preflight_merge.sh +++ b/scripts/preflight_merge.sh @@ -400,6 +400,10 @@ fi # unresolved thread on this PR with a valid 4-column row and an OPEN issue is # downgraded from OPEN to DEFER; resolved threads, unmatched URLs, malformed # rows, and closed or missing issues all refuse. +# Format scope: the register is parsed per physical line (split("\n")). Quoted +# rationales containing commas and escaped double quotes ("") are supported on a +# single physical line; embedded newlines inside quoted fields are intentionally +# unsupported (single-line rationales are required). printf '\n' ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" DEFERRALS_FILE="${PREFLIGHT_DEFERRALS:-$ROOT_DIR/data/deferrals.csv}" @@ -450,7 +454,8 @@ if [[ -f "$DEFERRALS_FILE" ]]; then | map(select((test("^\\s*$") | not) and (test("^\\s*#") | not) and (. != "pr,thread_url,issue,rationale"))) | .[] | . as $raw - | (capture("^(?[^,\"]*),(?[^,\"]*),(?[^,\"]*),(?:\\s*\"(?(?:[^\"]|\"\")*)\"\\s*|(?[^,\"]*))$") // null) as $m + | (([ $raw | scan("\"") ] | length) % 2 == 1) as $odd_quotes + | (if $odd_quotes then null else (capture("^(?[^,\"]*),(?[^,\"]*),(?[^,\"]*),(?:\\s*\"(?(?:[^\"]|\"\")*)\"\\s*|(?[^,\"]*))$") // null) end) as $m | if $m == null then "MALFORMED\t\($raw)" else @@ -480,7 +485,7 @@ if [[ -f "$DEFERRALS_FILE" ]]; then [[ -z "$kind" ]] && continue if [[ "$kind" == "MALFORMED" ]]; then printf ' MALFORMED DEFERRAL ROW in %s: %s\n' "$DEFERRALS_FILE" "$f1" - printf ' Every row requires pr (bare digits), thread_url, issue (bare digits, #N refused), and non-empty rationale.\n' + printf ' Every row requires pr (bare digits), thread_url, issue (bare digits, #N refused), and a single-line non-empty rationale (embedded newlines and unclosed quotes are refused).\n' fail=1 continue fi diff --git a/scripts/test_preflight_merge.sh b/scripts/test_preflight_merge.sh index fe3df1c..6a9cda8 100755 --- a/scripts/test_preflight_merge.sh +++ b/scripts/test_preflight_merge.sh @@ -759,6 +759,13 @@ run_deferral_fixture "[$(t_node P2 false "$URL1" "src/a.jl")]" \ 1,$URL1,21,\"unclosed quote rationale" verdict "a deferral row with an unclosed quote in rationale blocks" 1 $? "MALFORMED DEFERRAL ROW" +# Control (8d): quoted rationale with an embedded newline blocks as malformed -> exits 1. +run_deferral_fixture "[$(t_node P2 false "$URL1" "src/a.jl")]" \ +"pr,thread_url,issue,rationale +1,$URL1,21,\"first line of rationale +second line of rationale\"" +verdict "a deferral row with an embedded newline in a quoted rationale blocks" 1 $? "MALFORMED DEFERRAL ROW" + # Fake-gh controls for `gh issue view`: PREFLIGHT_FIXTURE bypasses `gh`, so these # run through a fake `gh` on PATH. Control (3) (OPEN issue clears) is placed FIRST # so a broken `gh issue view` invocation cannot make (1) and (2) pass by failing everything. diff --git a/tests/deterministic_radiation.jl b/tests/deterministic_radiation.jl index 30bfc77..a8d28b3 100644 --- a/tests/deterministic_radiation.jl +++ b/tests/deterministic_radiation.jl @@ -107,3 +107,16 @@ let @test mean(drift_cn) < 0 # negative-β cell moves toward the axis: a TROPISM @test mean(drift_so) > 0 # radiosensitive-signed cell moves outward end + +# ---------- (d) BiofilmsPottsRad two-way attenuation & burnout controls ---------- + +let + audit = SR.BiofilmsPottsRad.run_control_audit() + @test isapprox(audit.var_I_neg, 0.0; atol = 1e-12) + @test audit.var_I_pos > 1e-4 + @test audit.delta_r_sp4 < 0.0 + @test audit.delta_r_sp5 < 0.0 + @test audit.mean_I_voidy > audit.mean_I_dense + @test audit.pit_overshoot > 0.0 + @test audit.burnout_clearance > 0.0 +end