Add interactive GP concrete strength explorer (BOxCrete website) + PartialFixedNoiseLikelihood

.github/workflows/js-sync.yml

@@ -0,0 +1,40 @@
+# Verifies that the JavaScript GP implementation matches Python reference predictions.
+# Catches drift if model parameters or JS code are changed without re-syncing.
+
+name: JS Model Sync
+
+# Restrict GITHUB_TOKEN to the minimum needed: read-only access to repo
+# contents (required by actions/checkout). The workflow never writes to the
+# repo, comments on PRs, deploys, or interacts with any other GitHub APIs.
+permissions:
+  contents: read
+
+on:
+  push:
+    branches: [main, master]
+    paths:
+      - 'docs/gp.mjs'
+      - 'docs/model/**'
+      - 'test/test_js_gp.mjs'
+  pull_request:
+    branches: [main, master]
+    paths:
+      - 'docs/gp.mjs'
+      - 'docs/model/**'
+      - 'test/test_js_gp.mjs'
+
+env:
+  FORCE_JAVASCRIPT_ACTIONS_TO_NODE24: true
+
+jobs:
+  js-model-sync:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: actions/setup-node@v4
+        with:
+          node-version: '20'
+
+      - name: Run JS GP sync test
+        run: node test/test_js_gp.mjs

boxcrete/__init__.py

@@ -10,6 +10,7 @@
     fit_slump_gp,
     fit_strength_gp,
     get_strength_gp_input_transform,
+    PartialFixedNoiseLikelihood,
     SustainableConcreteModel,
 )
 from boxcrete.plotting import (
@@ -72,6 +73,7 @@
     "MORTAR_BOUNDS_DICT",
     "MORTAR_CONSTRAINTS",
     "MORTAR_REFERENCE_POINT",
+    "PartialFixedNoiseLikelihood",
     "SLUMP_DISPLAY_SCALE",
     "SLUMP_Y_COLUMNS",
     "STRENGTH_DISPLAY_SCALE",

boxcrete/models.py

@@ -34,6 +34,7 @@
 from gpytorch.kernels import MaternKernel, RBFKernel, ScaleKernel
 from gpytorch.likelihoods import GaussianLikelihood
 from gpytorch.mlls import ExactMarginalLogLikelihood
+from linear_operator.operators import DiagLinearOperator
 from torch import Tensor
 
 # Indices into DEFAULT_X_COLUMNS (without Time) for derived feature computation
@@ -415,6 +416,68 @@ def get_model_dict(
         return dict(zip(self.model_names, model_list.models))
 
 
+class PartialFixedNoiseLikelihood(GaussianLikelihood):
+    """Gaussian likelihood that learns noise for real observations while applying
+    fixed near-zero noise to pseudo-observations.
+
+    This enables conditioning the GP to pass through pseudo-observations (e.g.,
+    zero strength at time zero) with high certainty, while still learning the
+    observation noise for real data points via marginal likelihood optimization.
+
+    Args:
+        n_real: Number of real observations (must come first in training data).
+        n_pseudo: Number of pseudo-observations (must come last in training data).
+        pseudo_noise: Fixed noise variance for pseudo-observations.
+        **kwargs: Additional keyword arguments passed to GaussianLikelihood
+            (e.g., noise_constraint).
+    """
+
+    def __init__(
+        self,
+        n_real: int,
+        n_pseudo: int,
+        pseudo_noise: float = 1e-6,
+        **kwargs,
+    ):
+        super().__init__(**kwargs)
+        self._n_real = n_real
+        self._n_pseudo = n_pseudo
+        self._pseudo_noise = pseudo_noise
+
+    @property
+    def n_real(self) -> int:
+        return self._n_real
+
+    @property
+    def n_pseudo(self) -> int:
+        return self._n_pseudo
+
+    @property
+    def pseudo_noise(self) -> float:
+        return self._pseudo_noise
+
+    def _shaped_noise_covar(self, base_shape, *params, **kwargs):
+        n = base_shape[-1]
+        noise = self.noise_covar.noise.squeeze()  # learned scalar noise
+
+        if n == self._n_real + self._n_pseudo:
+            # Training: learned noise for real obs, fixed for pseudo-obs
+            diag = torch.cat(
+                [
+                    noise.expand(self._n_real),
+                    torch.full(
+                        (self._n_pseudo,),
+                        self._pseudo_noise,
+                        device=noise.device,
+                        dtype=noise.dtype,
+                    ),
+                ]
+            )
+            return DiagLinearOperator(diag)
+        # Prediction at test points: use learned noise
+        return super()._shaped_noise_covar(base_shape, *params, **kwargs)
+
+
 def fit_strength_gp(
     X: Tensor,
     Y: Tensor,
@@ -443,6 +506,8 @@ def fit_strength_gp(
 
     # add data to condition GP to be zero at day zero
     X_0, Y_0, Yvar_0 = get_day_zero_data(X=X, bounds=X_bounds, n=128)
+    n_real = X.shape[0]
+    n_pseudo = X_0.shape[0]
     X = torch.cat((X, X_0), dim=0)
     Y = torch.cat((Y, Y_0), dim=0)
     Yvar = torch.cat((Yvar, Yvar_0), dim=0)
@@ -484,8 +549,11 @@ def fit_strength_gp(
     if use_fixed_noise:
         model_kwargs["train_Yvar"] = Yvar
     else:
-        model_kwargs["likelihood"] = GaussianLikelihood(
-            noise_constraint=LogTransformedInterval(1e-6, 1.0, initial_value=1e-1)
+        model_kwargs["likelihood"] = PartialFixedNoiseLikelihood(
+            n_real=n_real,
+            n_pseudo=n_pseudo,
+            pseudo_noise=1e-6,
+            noise_constraint=LogTransformedInterval(1e-6, 1.0, initial_value=1e-1),
         )
     model = SingleTaskGP(**model_kwargs)
     mll = ExactMarginalLogLikelihood(model.likelihood, model)

boxcrete/utils.py

@@ -929,31 +929,26 @@ def get_subset_sum_tensors(
 # by `fit_gwp_model`.
 # Derived via per-class least-squares regression on training data (which
 # stores -GWP). Magnitudes here are the absolute emission factors.
-DEFAULT_GWP_COEFFICIENTS: dict[int, dict[str, tuple[float, float]]] = {
+DEFAULT_GWP_COEFFICIENTS = {
     0: {  # Material Source 0
-        "Cement (kg/m3)": (0.762613, 0.000384),
-        "Fly Ash (kg/m3)": (0.029577, 0.000457),
-        "Slag (kg/m3)": (0.085921, 0.000326),
-        "Water (kg/m3)": (0.001829, 0.001177),
-        # HRWR GWP: consistent with published EPDs for polycarboxylate-based
-        # superplasticizers (1.5-5.0 kg CO₂/kg range; e.g. BASF MasterGlenium,
-        # Sika ViscoCrete). Tight std confirms deterministic upstream formula.
-        "HRWR (kg/m3)": (3.184316, 0.016840),
-        "Fine Aggregate (kg/m3)": (0.002762, 0.000113),
-        "Coarse Aggregates (kg/m3)": (0.003895, 0.000121),
-        "Temp (C)": (0.002967, 0.005397),
+        "Cement (kg/m3)": (0.762610, 0.000365),
+        "Fly Ash (kg/m3)": (0.029601, 0.000432),
+        "Slag (kg/m3)": (0.085926, 0.000310),
+        "Water (kg/m3)": (-0.001765, 0.001114),
+        "HRWR (kg/m3)": (3.184692, 0.016001),
+        "Fine Aggregate (kg/m3)": (0.002788, 0.000097),
+        "Coarse Aggregates (kg/m3)": (0.003910, 0.000112),
     },
     1: {  # Material Source 1
-        "Cement (kg/m3)": (0.774398, 0.007709),
-        "Fly Ash (kg/m3)": (0.036826, 0.005956),
-        "Slag (kg/m3)": (0.094776, 0.007031),
-        "Water (kg/m3)": (0.001752, 0.020415),
-        # HRWR GWP: consistent with Source 0 (3.18 vs 3.10 kg CO₂/kg).
+        "Cement (kg/m3)": (0.773814, 0.007240),
+        "Fly Ash (kg/m3)": (0.035681, 0.005542),
+        "Slag (kg/m3)": (0.092849, 0.006469),
+        "Water (kg/m3)": (0.003864, 0.019144),
+        # HRWR GWP: consistent with Source 0 (3.18 vs 3.15 kg CO₂/kg).
         # See comment above for EPD references.
-        "HRWR (kg/m3)": (3.102591, 0.531138),
-        "Fine Aggregate (kg/m3)": (0.002823, 0.003714),
-        "Coarse Aggregates (kg/m3)": (0.001063, 0.003163),
-        "Temp (C)": (0.072522, 0.054983),
+        "HRWR (kg/m3)": (3.151231, 0.498391),
+        "Fine Aggregate (kg/m3)": (0.002513, 0.003486),
+        "Coarse Aggregates (kg/m3)": (-0.000039, 0.002870),
     },
 }
 
@@ -1016,31 +1011,40 @@ def make_linear_coefficients(
     return means, variances
 
 
-def get_day_zero_data(X: Tensor, bounds: Tensor | None, n: int = 128):
-    """Computes a tensor of n sobol points that satisfy the bounds, appended with a
-    zeros tensor. Useful to condition the strength GP to be zero at day zero.
+def get_day_zero_data(X: Tensor, bounds: Tensor | None = None, n: int = 128):
+    """Generates pseudo-observations at time=0 for conditioning the GP to predict
+    zero strength at day zero.
+
+    Uses the unique compositions from the training data (without time) to ensure
+    the constraint is enforced at all observed mix designs. If the number of unique
+    compositions exceeds n, a random subset is selected.
 
     Args:
-        X: The input tensor.
-        bounds: The bounds of the input tensor. If None, will be inferred from X.
-        n: The number of sobol points to generate.
+        X: The input tensor (n_train x d), where the last column is time.
+        bounds: Unused, kept for API compatibility. Will be removed in a future version.
+        n: Maximum number of pseudo-observations. If there are fewer unique
+            compositions than n, all unique compositions are used.
 
     Returns:
-        A tensor of n sobol points that satisfy the bounds, appended with a zeros
-        tensor, corresponding to the strength at day zero.
+        A tuple (X_0, Y_0, Yvar_0) of pseudo-observations at time=0.
     """
-    if bounds is None:
-        bounds = torch.stack((X.amin(dim=0), X.amax(dim=0)))
+    # Use unique observed compositions (without time)
+    unique_comps = torch.unique(X[:, :-1], dim=0)
+    n_unique = unique_comps.shape[0]
 
-    d = bounds.shape[-1]
-    sobol_engine = torch.quasirandom.SobolEngine(dimension=(d - 1))  # excluding time
-    X_0 = sobol_engine.draw(n)
-    X_0 = torch.cat((X_0, torch.zeros(n, 1)), dim=-1)  # append time (zero)
-    a, b = bounds[0], bounds[1]
-    X_0 = (b - a) * X_0 + a  # scaling according to bounds
-    X_0[:, -1] = 0.0  # explicitly set time to zero (day zero conditioning)
-    Y_0 = torch.zeros(n, 1)  #  zero strength
-    Yvar_0 = torch.full((n, 1), 1e-4)  #  with large certainty
+    if n_unique <= n:
+        # Use all unique compositions
+        X_comps = unique_comps
+    else:
+        # Random subset of unique compositions
+        perm = torch.randperm(n_unique)[:n]
+        X_comps = unique_comps[perm]
+
+    n_out = X_comps.shape[0]
+    # Append time=0
+    X_0 = torch.cat((X_comps, torch.zeros(n_out, 1, dtype=X.dtype)), dim=-1)
+    Y_0 = torch.zeros(n_out, 1, dtype=X.dtype)
+    Yvar_0 = torch.full((n_out, 1), 1e-4, dtype=X.dtype)
     return X_0, Y_0, Yvar_0