Unify mortar/concrete optimization, consolidate constraints, 100% test coverage

- Rename BOxCrete_data.csv → boxcrete_data.csv (lowercase, remove clay columns)
- Remove compressive_strength.csv (subsumed by boxcrete data)
- Rename strength_model_tutorial.ipynb → prediction_and_optimization_tutorial.ipynb
- Update all constants to boxcrete column names (with units)
- Consolidate get_mortar_bounds + get_bounds → single get_bounds with
  MORTAR_BOUNDS_DICT / CONCRETE_BOUNDS_DICT presets
- Consolidate get_mortar_constraints + get_concrete_constraints → single
  get_constraints with equality_sums parameter and MORTAR_CONSTRAINTS /
  CONCRETE_CONSTRAINTS presets
- Inline thin constraint wrapper functions into get_constraints
- Generalize get_model_list(fixed_features=None) for mortar/concrete
  optimization — backward compatible, wraps GWP model when non-Time
  features are fixed
- Add reduce_to_optimization_space for dimension reduction of bounds
  and constraints when using FixedFeatureModel
- Add boxcrete/plotting.py with plot_strength_curve utility
- Parameterize notebooks with optimization_mode env var (mortar/concrete)
- Fix demo notebook hardcoded tensors → data-driven column lookup
- CI executes notebooks in both mortar and concrete modes
- Add docstrings to all public methods and classes
- 79 tests, 100% coverage

Author: SebastianAment

.github/workflows/notebooks.yml

@@ -58,18 +58,23 @@ jobs:
         run: |
           python -m ipykernel install --user --name python3
 
-      - name: Execute notebooks
+      - name: Execute notebooks (mortar and concrete modes)
         run: |
-          python - << 'PYEOF'
-          import subprocess, sys
+          for mode in mortar concrete; do
+            echo "========================================"
+            echo "Running notebooks with BOXCRETE_OPTIMIZATION_MODE=$mode"
+            echo "========================================"
+            BOXCRETE_OPTIMIZATION_MODE=$mode python - << 'PYEOF'
+          import subprocess, sys, os
           from pathlib import Path
 
+          mode = os.environ.get("BOXCRETE_OPTIMIZATION_MODE", "concrete")
           notebooks = sorted(Path("notebooks").glob("*.ipynb"))
           failed = []
 
           for nb in notebooks:
               print(f"{'=' * 40}")
-              print(f"Executing: {nb}")
+              print(f"Executing ({mode}): {nb}")
               print(f"{'=' * 40}")
 
               result = subprocess.run(
@@ -87,19 +92,20 @@ jobs:
 
               if result.returncode != 0:
                   failed.append(str(nb))
-                  print(f"❌ Failed: {nb}")
+                  print(f"❌ Failed ({mode}): {nb}")
               else:
-                  print(f"✅ Successfully executed: {nb}")
+                  print(f"✅ Successfully executed ({mode}): {nb}")
               print()
 
           if failed:
-              print(f"\n{len(failed)} notebook(s) failed:")
+              print(f"\n{len(failed)} notebook(s) failed in {mode} mode:")
               for f in failed:
                   print(f"  - {f}")
               sys.exit(1)
           else:
-              print(f"\nAll {len(notebooks)} notebook(s) executed successfully.")
+              print(f"\nAll {len(notebooks)} notebook(s) executed successfully in {mode} mode.")
           PYEOF
+          done
 
       - name: Upload executed notebooks as artifacts
         uses: actions/upload-artifact@v4

.gitignore

@@ -56,5 +56,17 @@ Thumbs.db
 *.tar.gz
 *.zip
 
+# Legacy root-level modules (superseded by boxcrete/ package)
+/models.py
+/utils.py
+/input_transform.py
+/doodles.py
+/posterior_mean_pareto.py
+
+# Notebook-generated output files
+ConcreteFormulae.csv
+ConcreteFormula_*.png
+GWPvsStrength.png
+
 # Planning files
 *.plan.md

README.md

@@ -1,7 +1,7 @@
 # BOxCrete: A Bayesian Optimization open-source AI Model for Concrete Mix Design & Optimization
 
-Concrete, the second most widely used material in the world, accounts for **6–8% of global anthropogenic CO₂ emissions**, largely due to Portland cement production (~0.8 tons CO₂ per ton of cement). Partial replacement with Supplementary Cementitious Materials (SCMs) such as fly ash, slag, and natural pozzolan reduces embodied carbon and often improves durability, but high SCM usage makes compressive strength a highly nonlinear function of multiple interacting mix parameters, rendering traditional design empirical and trial-and-error driven. To systematically navigate this complex composition space, data-driven frameworks are needed. 
-Here, we introduce BOxCrete, an open-source Bayesian optimization framework for probabilistic strength curve prediction and sustainable mix design. 
+Concrete, the second most widely used material in the world, accounts for **6–8% of global anthropogenic CO₂ emissions**, largely due to Portland cement production (~0.8 tons CO₂ per ton of cement). Partial replacement with Supplementary Cementitious Materials (SCMs) such as fly ash, slag, and natural pozzolan reduces embodied carbon and often improves durability, but high SCM usage makes compressive strength a highly nonlinear function of multiple interacting mix parameters, rendering traditional design empirical and trial-and-error driven. To systematically navigate this complex composition space, data-driven frameworks are needed.
+Here, we introduce BOxCrete, an open-source Bayesian optimization framework for probabilistic strength curve prediction and sustainable mix design.
 We invite researchers and practitioners of both machine learning and civil engineering
 to collaborate on discovering more sustainable concrete formulations that are applicable
 to a wide array of construction projects, at scale.
@@ -13,9 +13,11 @@ This repository contains probabilistic models and data for the
 1) Compressive strength of concrete and mortar mixes
 2) The associated global warming potential (GWP)
 
-as a function of their composition, consisting of
-cement, slag, water, to name a few basic ingredients.
-See `boxcrete/models.py` for implementation details.
+as a function of their composition, consisting of cement, fly ash, slag, fine and coarse aggregate, admixtures, and water, to name a few basic ingredients. See `boxcrete/models.py` for implementation details.
+
+### Included Datasets
+
+- **BOxCrete data** (`data/boxcrete_data.csv`): Combined mortar and concrete mix compositions with strength measurements at multiple curing ages, GWP values, and multiple material sources. This is the single unified dataset used for all model training.
 
 ## Installation
 
@@ -45,7 +47,7 @@ pip install -e ".[notebooks]"
 
 ```python
 from boxcrete.models import SustainableConcreteModel
-from boxcrete.utils import load_concrete_strength, get_mortar_bounds
+from boxcrete.utils import load_concrete_strength, get_bounds
 ```
 
 The models can be used for a variety of tasks, including but not limited to
@@ -58,7 +60,7 @@ The models can be used for a variety of tasks, including but not limited to
 
 ## Compressive Strength Model
 
-The `SustainableConcreteModel` in ['BOxCrete_models.py'](BOxCrete_models.py) includes a strength_model that predicts the evolution of compressive strength as a function of mixture composition. A tutorial is provided in [notebooks/BOxCrete Concrete Strength Prediction for GitHub.ipynb](<notebooks/BOxCrete Concrete Strength Prediction for GitHub.ipynb>), which demonstrates how the model can be used to predict the full strength development curve for any user-specified mix. The model is based on Gaussian Process (GP) regression and incorporates custom modeling steps to ensure physically consistent strength evolution and calibrated uncertainty. Example strength curve predictions generated using the notebook are shown in the figures below.
+The `SustainableConcreteModel` in [`boxcrete/models.py`](boxcrete/models.py) includes a strength_model that predicts the evolution of compressive strength as a function of mixture composition. A demo is provided in [`notebooks/strength_curve_prediction_demo.ipynb`](notebooks/strength_curve_prediction_demo.ipynb), which demonstrates how the model can be used to predict the full strength development curve for any user-specified mix. A comprehensive tutorial covering prediction, calibration, Pareto frontiers, and gradient-based experimental design is available in [`notebooks/prediction_and_optimization_tutorial.ipynb`](notebooks/prediction_and_optimization_tutorial.ipynb). The model is based on Gaussian Process (GP) regression and incorporates custom modeling steps to ensure physically consistent strength evolution and calibrated uncertainty. Example strength curve predictions generated using the notebook are shown in the figures below.
 
 <p align="center">
   <img src="fig/Picture1.png">
@@ -80,7 +82,7 @@ Further, when trained on the mortar and concrete mix strength data contained in
 
 ## Experimental Design
 
-The probabilistic compressive strength model can also be used to design new concrete mixtures that achieve optimal trade-offs between mechanical performance and environmental impact. In particular, the framework enables multi-objective optimization of early-age (1-day) and later-age (28-day) compressive strength alongside Global Warming Potential (GWP). By systematically exploring the composition space, BOxCrete can generate candidate mixes that balance structural performance requirements with carbon reduction targets. 
+The probabilistic compressive strength model can also be used to design new concrete mixtures that achieve optimal trade-offs between mechanical performance and environmental impact. In particular, the framework enables multi-objective optimization of early-age (1-day) and later-age (28-day) compressive strength alongside Global Warming Potential (GWP). By systematically exploring the composition space, BOxCrete can generate candidate mixes that balance structural performance requirements with carbon reduction targets.
 
 As illustrated in the figure below, the model identifies a Pareto front capturing the trade-off between 1-day strength, 28-day strength, and GWP across candidate mixtures.
 

boxcrete/__init__.py

@@ -11,33 +11,47 @@
     fit_strength_gp,
     get_strength_gp_input_transform,
 )
+from boxcrete.plotting import plot_strength_curve
 from boxcrete.utils import (
-    DEFAULT_DATA_PATH,
+    CONCRETE_BOUNDS_DICT,
+    CONCRETE_CONSTRAINTS,
+    DATA_PATH,
+    DEFAULT_BOUNDS_DICT,
     DEFAULT_X_COLUMNS,
     DEFAULT_Y_COLUMNS,
     DEFAULT_YSTD_COLUMNS,
+    MORTAR_BOUNDS_DICT,
+    MORTAR_CONSTRAINTS,
     SustainableConcreteDataset,
+    get_bounds,
+    get_constraints,
     get_day_zero_data,
-    get_mortar_bounds,
-    get_mortar_constraints,
     get_reference_point,
     load_concrete_strength,
+    reduce_to_optimization_space,
 )
 
 __all__ = [
-    "FixedFeatureModel",
-    "SustainableConcreteModel",
-    "fit_gwp_gp",
-    "fit_strength_gp",
-    "get_strength_gp_input_transform",
-    "DEFAULT_DATA_PATH",
+    "CONCRETE_BOUNDS_DICT",
+    "CONCRETE_CONSTRAINTS",
+    "DATA_PATH",
+    "DEFAULT_BOUNDS_DICT",
     "DEFAULT_X_COLUMNS",
     "DEFAULT_Y_COLUMNS",
     "DEFAULT_YSTD_COLUMNS",
+    "FixedFeatureModel",
+    "MORTAR_BOUNDS_DICT",
+    "MORTAR_CONSTRAINTS",
+    "SustainableConcreteModel",
     "SustainableConcreteDataset",
+    "fit_gwp_gp",
+    "fit_strength_gp",
+    "get_bounds",
+    "get_constraints",
     "get_day_zero_data",
-    "get_mortar_bounds",
-    "get_mortar_constraints",
     "get_reference_point",
+    "get_strength_gp_input_transform",
     "load_concrete_strength",
+    "plot_strength_curve",
+    "reduce_to_optimization_space",
 ]

boxcrete/models.py

@@ -33,6 +33,15 @@
 
 
 class SustainableConcreteModel:
+    """Multi-output model that jointly predicts GWP and compressive strength.
+
+    The model consists of a GWP model (independent of curing time) and a strength
+    model (dependent on composition *and* time).  At optimisation time the strength
+    model is sliced at each of the ``strength_days`` via ``FixedFeatureModel`` to
+    produce a ``ModelList`` that maps composition only to ``[GWP, 1-day strength,
+    28-day strength, ...]``.
+    """
+
     def __init__(
         self,
         strength_days: list[int],
@@ -100,37 +109,82 @@ def _set_d(self, d: int) -> None:
         if self.d is None:
             self.d = d
 
-    def get_model_list(self) -> ModelList:
-        """Returns a ModelList modeling the GWP and compressive strength objectives as a function
-        of composition only.
-        Converts the strength and gwp models into a model list of independent models for gwp,
-        and x-day strengths, by fixing the time input of the strength model at 1 and 28 days.
+    def get_model_list(
+        self, fixed_features: dict[int, float] | None = None
+    ) -> ModelList:
+        """Returns a ``ModelList`` modelling GWP and compressive strength as a
+        function of composition only.
+
+        Converts the strength and GWP models into a model list of independent
+        models for GWP and *x*-day strengths by fixing the time input of the
+        strength model at each ``strength_day``.
+
+        Args:
+            fixed_features: Optional mapping from input column **index** to a
+                fixed value.  When provided these features are fixed *in
+                addition to* the Time dimension for the strength models, and
+                the non-Time entries are also applied to the GWP model via
+                ``FixedFeatureModel``.  Useful for fixing e.g.
+                ``Coarse Aggregates = 0`` in mortar mode.
+
+        Returns:
+            A ``ModelList`` with ``1 + len(strength_days)`` sub-models.
+
+        Raises:
+            ValueError: If the model has not been fitted yet.
         """
         if self.d is None:
             raise ValueError("Model not fit yet.")
-        models = [
-            self.gwp_model,
-            *(
+
+        time_idx = self.d - 1  # last column is Time
+
+        if fixed_features is None:
+            gwp_model = self.gwp_model
+        else:
+            non_time = {k: v for k, v in fixed_features.items() if k != time_idx}
+            if non_time:
+                gwp_indices = sorted(non_time.keys())
+                gwp_values = [non_time[i] for i in gwp_indices]
+                gwp_model = FixedFeatureModel(
+                    base_model=self.gwp_model,
+                    dim=self.d - 1,  # GWP input has no Time
+                    indices=gwp_indices,
+                    values=gwp_values,
+                )
+            else:
+                gwp_model = self.gwp_model
+
+        strength_models = []
+        for day in self.strength_days:
+            indices = [time_idx]
+            values: list[float] = [float(day)]
+            if fixed_features is not None:
+                for idx, val in sorted(fixed_features.items()):
+                    if idx != time_idx:
+                        indices.append(idx)
+                        values.append(val)
+            strength_models.append(
                 FixedFeatureModel(
                     base_model=self.strength_model,
                     dim=self.d,
-                    indices=[self.d - 1],
-                    values=[day],
+                    indices=indices,
+                    values=values,
                 )
-                for day in self.strength_days
-            ),
-        ]
-        model = ModelList(*models)
-        return model  # for use with multi-objective optimization
+            )
+
+        return ModelList(gwp_model, *strength_models)
 
     # TODO: add plot_strength_curve utility for visualizing predicted strength
     # curves as a function of time for a given composition.
 
 
 class FixedFeatureModel(Model):
-    # advantage: only need to implement posterior for it to work with qNEHI
-    # disadvantage: makes the strength outputs independent (IDEA: could add joint model)
-    # TODO: check that these are appended before the InputTransforms are applied, not after.
+    """Wraps a GP model to fix a subset of inputs to constant values.
+
+    At evaluation time the fixed features are spliced back into the input
+    tensor before delegating to the ``base_model``.
+    """
+
     def __init__(
         self,
         base_model: Model,
@@ -202,9 +256,19 @@ def posterior(self, X: Tensor, *args, **kwargs) -> Posterior:
 
     @property
     def num_outputs(self) -> int:
-        return self.base_model.num_outputs  # need to adjust if we batch fixed features
+        """The number of outputs of the base model."""
+        return self.base_model.num_outputs
 
     def subset_output(self, idcs: list[int]) -> FixedFeatureModel:
+        """Returns a new ``FixedFeatureModel`` whose base model is subset to
+        the given output indices.
+
+        Args:
+            idcs: Output indices to keep.
+
+        Returns:
+            A ``FixedFeatureModel`` wrapping the subset base model.
+        """
         return FixedFeatureModel(
             base_model=self.base_model.subset_output(idcs),
             dim=self._dim,