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.

Explorer performance
- predictStrengthCurveV2 gains a `meanOnly` option that skips the
  Cholesky/dtrsm variance solve and the [n x nTimes] K allocation. The
  mean is already accumulated during the kernel build, so output is
  bit-identical; only discarded work is removed. Routed the two callers
  that render no uncertainty band -- the dashed preview curve and the
  scatter position marker -- through it. Preview at 48 points: 5.94 ms
  -> 2.86 ms (2.1x).
- drawStrengthCurve's `isInteracting` now also covers preview activity.
  Previously, while the preview curve animated at 60 fps the main curve
  stayed on the 64-point grid because neither a slider nor a composition
  animation was in flight, making preview-settling the most expensive
  frame in the app. Behaviour at rest is unchanged (showPreview is false,
  so the expression evaluates exactly as before and the curve is still
  sampled at 64 points).
- Measured per frame, n_train=670, WASM active, desktop arm64:
  scatter hover / preview settling 14.08 ms -> 6.91 ms (51% less);
  animated transition 10.07 ms -> 6.91 ms (31% less); sustained slider
  drag is unaffected (4.13 -> 4.04 ms) because the preview does not run
  during a drag.

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",