Three-class Material Source strength GP

The BOxCrete dataset spans three physically distinct material sets, not
two: Set 1 mortars (Amrize 1L / Class C fly ash, no coarse aggregate),
Set 2 concretes (Heidelberg 1L / Class C), and Set 3 concretes (Amrize
1L / Class F). The deployed model previously treated Material Source as
a continuous ARD coordinate, which forces a single lengthscale to span
both the Set 1 <-> Set 2 and the Set 1+2 <-> Set 3 gaps -- ill-defined
when the axis is really an unordered categorical.

Model
- Material Source becomes a 3-level categorical. The source-aware branch
  is now ScaleKernel(CategoricalKernel(source) * Matern(rest)), a Bonilla
  et al. (2008) intrinsic-coregionalization factorisation: one shared
  Matern shape plus a single learned cross-class correlation
  rho = exp(-1/ell). Full kernel is blind_matern + source_branch +
  additive_rbf_time, wrapped in TimeGatedKernel (gate tau = 0.10) so the
  prior variance at t=0 is exactly zero and f(x, 0) = 0 holds
  structurally, without day-zero anchor pseudo-observations.
- Training minimises a combined block-LOO + marginal-likelihood
  objective (boxcrete/block_loo.py). Blocks are unique compositions, so
  the loss measures held-out-composition curve prediction -- what the
  explorer and the BO loop actually consume. The MLL term contributes
  the -0.5 log|K| complexity penalty that block-LOO structurally lacks,
  keeping curves smooth between training ages.
- ComposedLengthscalePrior lets the within-group shrinkage prior
  (Cement/FA/Slag tied, Fine/Coarse Aggregate tied) compose with the
  per-dim lognormal baseline.

Data and naming
- boxcrete/mix_naming.py decodes the 3 classes to canonical mix names
  from composition alone, backed by _mix_naming_table.csv: M1..M69
  (69 Set 1 mortars), C1..C27 (27 Set 2 concretes), C28..C80 (53 Set 3
  concretes), 149 mixes total.
- data/boxcrete_data.csv carries the Set 3 rows and the 3-class Material
  Source encoding. Clay-using mortars M75/76/77 are dropped.
- docs/materials_background.md documents the chemistry behind each class.

Explorer
- docs/ui.mjs exposes Material Source as a 3-way selector; gp.mjs and
  gp_v2_fast.mjs implement the categorical branch for JS-Python parity.
- Regenerated docs/model/* artifacts (strength_model.pt, test_vectors,
  compositions, mix_analyses) via experiments/regenerate_all_artifacts.sh.
- docs/model/mix_analyses.json is re-authored for the 3-class catalog.
  The renumbering left the previous narratives keyed to a catalog order
  that no longer exists (only 31 of 144 entries still index-aligned, and
  70 changed Material Source label), so they could not be carried over.
  Each of the 149 entries now names its canonical mix, states its
  verified mix-design figures, and interprets its behaviour against
  sibling mixes in the same designed series. The entries deliberately do
  not restate measured strength points, embodied carbon, or Pareto
  status: the explorer already overlays observations on the strength
  curve, gives GWP/cost/W-B their own readouts, and shows Pareto
  membership as a pill, so repeating them only padded every entry.
  regenerate_all_artifacts.sh no longer invokes
  docs/generate_mix_analyses.py, whose templated fallback output would
  overwrite the authored prose -- this restores the contract that
  script's own docstring already specified.

Tests
- test_composed_prior.py, test_mix_naming.py cover the new modules.
- test_kernel_layout.py and test_lengthscale_identifiability.py updated
  for the categorical branch.
- Full suite: 270 passing, 100% line coverage on boxcrete/.

Author: SebastianAment

boxcrete/_mix_naming_table.csv

@@ -0,0 +1,150 @@
+legacy_int_name,canonical_name,material_class,source_evidence
+M1,M1,0,boxcrete_all direct (v6)
+M2,M2,0,boxcrete_all direct (v6)
+M3,M3,0,boxcrete_all direct (v6)
+M4,M4,0,boxcrete_all direct (v6)
+M5,M5,0,boxcrete_all direct (v6)
+M6,M6,0,boxcrete_all direct (v6)
+M7,M7,0,boxcrete_all direct (v6)
+M8,M8,0,boxcrete_all direct (v6)
+M9,M9,0,boxcrete_all direct (v6)
+M10,M10,0,boxcrete_all direct (v6)
+M11,M11,0,boxcrete_all direct (v6)
+M12,M12,0,boxcrete_all direct (v6)
+M13,M13,0,boxcrete_all direct (v6)
+M14,M14,0,boxcrete_all direct (v6)
+M15,M15,0,boxcrete_all direct (v6)
+M16,M16,0,boxcrete_all direct (v6)
+M17,M17,0,boxcrete_all direct (v6)
+M18,M18,0,boxcrete_all direct (v6)
+M19,M19,0,boxcrete_all direct (v6)
+M20,M20,0,boxcrete_all direct (v6)
+M21,M21,0,boxcrete_all direct (v6)
+M22,M22,0,boxcrete_all direct (v6)
+M23,M23,0,boxcrete_all direct (v6)
+M24,M24,0,boxcrete_all direct (v6)
+M25,M25,0,boxcrete_all direct (v6)
+M26,M26,0,boxcrete_all direct (v6)
+M27,M27,0,boxcrete_all direct (v6)
+M28,M28,0,boxcrete_all direct (v6)
+M29,M29,0,boxcrete_all direct (v6)
+M30,M30,0,boxcrete_all direct (v6)
+M31,M31,0,boxcrete_all direct (v6)
+M32,M32,0,boxcrete_all direct (v6)
+M33,M33,0,boxcrete_all direct (v6)
+M34,M34,0,boxcrete_all direct (v6)
+M35,M35,0,boxcrete_all direct (v6)
+M36,M36,0,boxcrete_all direct (v6)
+M37,M37,0,boxcrete_all direct (v6)
+M38,M38,0,boxcrete_all direct (v6)
+M39,M39,0,boxcrete_all direct (v6)
+M40,M40,0,boxcrete_all direct (v6)
+M41,M41,0,boxcrete_all direct (v6)
+M42,M42,0,boxcrete_all direct (v6)
+M43,M43,0,boxcrete_all direct (v6)
+M44,M44,0,boxcrete_all direct (v6)
+M45,M45,0,boxcrete_all direct (v6)
+M46,M46,0,boxcrete_all direct (v6)
+M47,M47,0,boxcrete_all direct (v6)
+M48,M48,0,boxcrete_all direct (v6)
+M49,M49,0,boxcrete_all direct (v6)
+M50,M50,0,boxcrete_all direct (v6)
+M51,M51,0,boxcrete_all direct (v6)
+M52,M52,0,boxcrete_all direct (v6)
+M53,M53,0,boxcrete_all direct (v6)
+M54,M54,0,boxcrete_all direct (v6)
+M55,M55,0,boxcrete_all direct (v6)
+M56,M56,0,boxcrete_all direct (v6)
+M57,M57,0,boxcrete_all direct (v6)
+M58,M58,0,boxcrete_all direct (v6)
+M59,M59,0,boxcrete_all direct (v6)
+M60,M60,0,boxcrete_all direct (v6)
+M61,M61,0,boxcrete_all direct (v6)
+M62,M62,0,boxcrete_all direct (v6)
+M63,M63,0,boxcrete_all direct (v6)
+M64,M64,0,boxcrete_all direct (v6)
+M65,M65,0,boxcrete_all direct (v6)
+M66,M66,0,boxcrete_all direct (v6)
+M67,M67,0,boxcrete_all direct (v6)
+M68,M68,0,boxcrete_all direct (v6)
+M69,M69,0,boxcrete_all direct (v6)
+C1,C1,1,boxcrete_all direct (v6)
+C2,C2,1,boxcrete_all direct (v6)
+C3,C3,1,boxcrete_all direct (v6)
+C4,C4,1,boxcrete_all direct (v6)
+C5,C5,1,boxcrete_all direct (v6)
+C6,C6,1,boxcrete_all direct (v6)
+C7,C7,1,boxcrete_all direct (v6)
+C8,C8,1,boxcrete_all direct (v6)
+C9,C9,1,boxcrete_all direct (v6)
+C10,C10,1,boxcrete_all direct (v6)
+C11,C11,1,boxcrete_all direct (v6)
+C12,C12,1,boxcrete_all direct (v6)
+C13,C13,1,boxcrete_all direct (v6)
+C14,C14,1,boxcrete_all direct (v6)
+C15,C15,1,boxcrete_all direct (v6)
+C16,C16,1,boxcrete_all direct (v6)
+C17,C17,1,boxcrete_all direct (v6)
+C18,C18,1,boxcrete_all direct (v6)
+C19,C19,1,boxcrete_all direct (v6)
+C20,C20,1,boxcrete_all direct (v6)
+C21,C21,1,boxcrete_all direct (v6)
+C22,C22,1,boxcrete_all direct (v6)
+C23,C23,1,boxcrete_all direct (v6)
+C24,C24,1,boxcrete_all direct (v6)
+C25,C25,1,boxcrete_all direct (v6)
+C26,C26,1,boxcrete_all direct (v6)
+C27,C27,1,boxcrete_all direct (v6)
+C28,C28,2,boxcrete_all direct (v6)
+C29,C29,2,boxcrete_all direct (v6)
+C30,C30,2,boxcrete_all direct (v6)
+C31,C31,2,boxcrete_all direct (v6)
+C32,C32,2,boxcrete_all direct (v6)
+C33,C33,2,boxcrete_all direct (v6)
+C34,C34,2,boxcrete_all direct (v6)
+C35,C35,2,boxcrete_all direct (v6)
+C36,C36,2,boxcrete_all direct (v6)
+C37,C37,2,boxcrete_all direct (v6)
+C38,C38,2,boxcrete_all direct (v6)
+C39,C39,2,boxcrete_all direct (v6)
+C40,C40,2,boxcrete_all direct (v6)
+C41,C41,2,boxcrete_all direct (v6)
+C42,C42,2,boxcrete_all direct (v6)
+C43,C43,2,boxcrete_all direct (v6)
+C44,C44,2,boxcrete_all direct (v6)
+C45,C45,2,boxcrete_all direct (v6)
+C46,C46,2,boxcrete_all direct (v6)
+C47,C47,2,boxcrete_all direct (v6)
+C48,C48,2,boxcrete_all direct (v6)
+C49,C49,2,boxcrete_all direct (v6)
+C50,C50,2,boxcrete_all direct (v6)
+C51,C51,2,boxcrete_all direct (v6)
+C52,C52,2,boxcrete_all direct (v6)
+C53,C53,2,boxcrete_all direct (v6)
+C54,C54,2,boxcrete_all direct (v6)
+C55,C55,2,boxcrete_all direct (v6)
+C56,C56,2,boxcrete_all direct (v6)
+C57,C57,2,boxcrete_all direct (v6)
+C58,C58,2,boxcrete_all direct (v6)
+C59,C59,2,boxcrete_all direct (v6)
+C60,C60,2,boxcrete_all direct (v6)
+C61,C61,2,boxcrete_all direct (v6)
+C62,C62,2,boxcrete_all direct (v6)
+C63,C63,2,boxcrete_all direct (v6)
+C64,C64,2,boxcrete_all direct (v6)
+C65,C65,2,boxcrete_all direct (v6)
+C66,C66,2,boxcrete_all direct (v6)
+C67,C67,2,boxcrete_all direct (v6)
+C68,C68,2,boxcrete_all direct (v6)
+C69,C69,2,boxcrete_all direct (v6)
+C70,C70,2,boxcrete_all direct (v6)
+C71,C71,2,boxcrete_all direct (v6)
+C72,C72,2,boxcrete_all direct (v6)
+C73,C73,2,boxcrete_all direct (v6)
+C74,C74,2,boxcrete_all direct (v6)
+C75,C75,2,boxcrete_all direct (v6)
+C76,C76,2,boxcrete_all direct (v6)
+C77,C77,2,boxcrete_all direct (v6)
+C78,C78,2,boxcrete_all direct (v6)
+C79,C79,2,boxcrete_all direct (v6)
+C80,C80,2,boxcrete_all direct (v6)

boxcrete/block_loo.py

@@ -0,0 +1,165 @@
+# Copyright (c) Meta Platforms, Inc. and affiliates.
+#
+# This source code is licensed under the MIT license found in the
+# LICENSE file in the root directory of this source tree.
+
+"""Block-leave-one-out (block-LOO) training objective for the strength GP.
+
+The deployed strength model is trained by minimising the **block-LOO
+predictive negative-log-likelihood** rather than the marginal likelihood
+(MLL). Each block is one unique composition (rows sharing the leading
+composition dims), so block-LOO measures how well the model predicts a
+*held-out composition's* strength curve — the quantity the explorer and
+Bayesian-optimisation loop actually rely on. Training directly on this
+objective (from initialisation, no MLL warm-up) gave the best and most
+stable held-out RMSE in the design-exploration benchmark.
+
+Closed form (Sundararajan-Keerthi block-inverse identity): with
+``K = K_kernel + sigma^2 I``, ``alpha = K^{-1}(y - mu_prior)`` and, for a
+block ``B``,
+
+    residual_B = (K^{-1}_BB)^{-1} @ alpha_B
+    var_B      = diag((K^{-1}_BB)^{-1})     # includes aleatoric noise
+
+so the per-row loss is the Gaussian NLL ``0.5*(log 2pi + log var + r^2/var)``.
+Day-zero anchor pseudo-rows (indices ``[n_real:]``) join their
+composition's block — holding a composition out also holds out its anchor
+(no leak) — but are not scored.
+
+The deployed model is trained on a **single combined objective** that
+blends this block-LOO loss with a fraction of the marginal likelihood
+(see :func:`train_block_loo`): the MLL term supplies the ``-0.5 log|K|``
+complexity penalty that block-LOO structurally lacks, keeping the strength
+curves smooth/monotone between the training ages at negligible block-LOO
+cost.
+"""
+
+from __future__ import annotations
+
+import math
+
+import torch
+from botorch.models import SingleTaskGP
+from gpytorch.utils.cholesky import psd_safe_cholesky
+
+__all__ = ["block_loo_loss", "train_block_loo"]
+
+
+def _model_prior_log_prob(model: SingleTaskGP) -> torch.Tensor:
+    """Sum of ``log p(theta)`` over every prior registered on the model.
+
+    Including this term makes block-LOO training maximise the posterior
+    ``log p(theta | y)`` rather than the bare predictive likelihood, which
+    keeps the lengthscale-estimation problem well-posed (the within-group
+    shrinkage prior stays active). ``named_priors()`` yields
+    ``(name, parent_module, prior, closure, setting_closure)``; the closure
+    takes the parent module (the kernel the prior was registered on).
+    """
+    total = torch.zeros((), dtype=torch.double)
+    for _name, parent_module, prior, closure, _setting in model.named_priors():
+        contrib = prior.log_prob(closure(parent_module)).sum()
+        total = total.to(contrib) + contrib
+    return total
+
+
+def block_loo_loss(
+    model: SingleTaskGP,
+    n_real: int,
+    *,
+    n_composition_dims: int = 9,
+) -> torch.Tensor:
+    """Differentiable mean block-LOO negative-log-likelihood per real row.
+
+    Args:
+        model: a ``SingleTaskGP`` whose hyperparameters require grad.
+        n_real: number of real rows; anchor pseudo-rows live at
+            ``[n_real:]`` and are held out with their block but not scored.
+        n_composition_dims: number of leading input dims that define a
+            block (default 9 = composition + temperature).
+
+    Returns a scalar tensor (minimise it to improve held-out calibration).
+    """
+    train_X = model.train_inputs[0]
+    train_Y = model.train_targets
+    prior = model(train_X)
+    noisy = model.likelihood(prior, train_X)
+    K = noisy.lazy_covariance_matrix.to_dense()
+
+    n = K.shape[-1]
+    L = psd_safe_cholesky(K)
+    Y_t = train_Y.unsqueeze(-1) if train_Y.dim() == 1 else train_Y
+    residuals = Y_t - prior.mean.unsqueeze(-1)
+    alpha = torch.cholesky_solve(residuals, L).squeeze(-1)
+    K_inv = torch.cholesky_solve(torch.eye(n, dtype=K.dtype, device=K.device), L)
+
+    fingerprints = train_X[..., :n_composition_dims]
+    unique_fp, inverse = torch.unique(fingerprints, dim=0, return_inverse=True)
+
+    log_two_pi = math.log(2.0 * math.pi)
+    total_nll = torch.zeros((), dtype=K.dtype, device=K.device)
+    n_scored = 0
+    for g in range(unique_fp.shape[0]):
+        idx = (inverse == g).nonzero(as_tuple=True)[0]
+        K_inv_BB_inv = torch.linalg.inv(K_inv[idx][:, idx])
+        block_residual = K_inv_BB_inv @ alpha[idx]
+        block_var = torch.diagonal(K_inv_BB_inv).clamp_min(1e-12)
+        block_nll = 0.5 * (
+            log_two_pi + torch.log(block_var) + (block_residual**2) / block_var
+        )
+        real_mask = idx < n_real
+        total_nll = total_nll + block_nll[real_mask].sum()
+        n_scored += int(real_mask.sum())
+
+    nll_per_row = total_nll / n_scored
+    return nll_per_row - _model_prior_log_prob(model) / n_scored
+
+
+def train_block_loo(
+    model: SingleTaskGP,
+    n_real: int,
+    *,
+    mll_weight: float = 0.5,
+    max_iter: int = 150,
+    lr: float = 0.1,
+) -> float:
+    """Train the model's hyperparameters from their current initialisation
+    on a single combined objective, minimised with LBFGS::
+
+        L(theta) = (1 - mll_weight) * block_loo_NLL(theta)
+                   +      mll_weight * MLL_NLL(theta)
+
+    The block-LOO term optimises held-out predictive calibration (the
+    deployment metric), but is scored only at the training ages, so on its
+    own it leaves the predictive mean free to *oscillate between* those ages
+    (non-monotone strength curves). The marginal-likelihood (MLL) term adds
+    back the ``-0.5 log|K|`` Occam/complexity penalty that block-LOO
+    structurally lacks; a modest ``mll_weight`` restores smooth, monotone
+    curves at negligible block-LOO cost. Trained from a freshly-constructed
+    model (no MLL warm-up) this is a genuine single-objective fit, not a
+    two-stage MLL-then-refine. Returns the final (pure) block-LOO loss.
+    """
+    from gpytorch.mlls import ExactMarginalLogLikelihood
+
+    model.train()
+    mll = ExactMarginalLogLikelihood(model.likelihood, model)
+    params = [p for p in model.parameters() if p.requires_grad]
+    optimizer = torch.optim.LBFGS(
+        params,
+        lr=lr,
+        max_iter=max_iter,
+        tolerance_grad=1e-6,
+        tolerance_change=1e-9,
+        line_search_fn="strong_wolfe",
+    )
+
+    def closure():
+        optimizer.zero_grad()
+        output = model(*model.train_inputs)
+        mll_nll = -mll(output, model.train_targets)
+        loss = (1.0 - mll_weight) * block_loo_loss(model, n_real) + mll_weight * mll_nll
+        loss.backward()
+        return loss
+
+    optimizer.step(closure)
+    with torch.no_grad():
+        return block_loo_loss(model, n_real).item()

boxcrete/features.py

@@ -8,7 +8,7 @@
 Public extension surface for users / variant authors:
 
   * :data:`F5_ALLLOG_FEATURES` — the deployed V2 strength GP's 7-feature set.
-  * :data:`GATE_TAU` — the production time-gate timescale (0.05).
+  * :data:`GATE_TAU` — the production time-gate timescale (0.10).
   * :data:`IDX` — column-name → integer-index map for the 10 raw input
     dims (Cement, Fly Ash, Slag, Water, HRWR, Fine, Coarse, Source,
     Temperature, Time).
@@ -38,13 +38,25 @@
 from boxcrete.utils import DEFAULT_X_COLUMNS
 
 # Time gate constant: h(t) = 1 - exp(-t / GATE_TAU).
-# tau=0.05 (post-input-transform time units) was found to be optimal in
-# the τ-sweep; see §4.5 of the benchmark.
-GATE_TAU = 0.05
-
-
-# Default feature set for the V2 strength GP. See STRENGTH_GP_BENCHMARK.md
-# §4.3 for the per-feature ablation showing all 5 log-transforms help.
+# tau=0.10 (post-input-transform time units). Chosen from a τ-sweep over
+# {0.05, 0.10, 0.15, 0.20} as the best point that keeps EXACT JS↔Python
+# parity. Raising τ reduces the early-time (t<1 day) negative overshoot in
+# the data-free window before the first (t=1 day) measurement — the min
+# predicted strength there improves from −1957 psi (τ=0.05) to −1302 psi
+# (τ=0.10), −975 (0.15), −809 (0.2) — but it also grows the fitted
+# outputscale and worsens the training-kernel conditioning, which makes the
+# in-browser JS posterior drift from the Python posterior. That drift stays
+# within the tight JS-port test (test/test_js_gp.mjs, rtol=1e-4/atol=1e-2)
+# only for τ ≤ 0.10; τ ≥ 0.15 breaks it (mean/variance drift up to ~0.5 psi
+# at 0.15, ~3.5 psi at 0.2). Held-out block-LOO loss is flat across the
+# sweep (−1.6509 at 0.10 vs −1.6529 best at 0.15, within noise), so τ=0.10
+# captures ~1/3 of the achievable overshoot reduction at zero parity cost.
+# Revisit (with a parity fix) if more early-window damping is needed.
+GATE_TAU = 0.1
+
+
+# Default feature set for the V2 strength GP. The per-feature ablation
+# found that all 5 log-transforms help.
 F5_ALLLOG_FEATURES = (
     "wb_ratio",
     "scm_frac",