Fix spectralog benchmark cases

benchmark_cases/README.md

@@ -23,7 +23,10 @@ Current tracked campaigns:
   GAMS profile: DICOPT for transformed/local MINLP roles, IPOPTH for NLP roles,
   and Gurobi for MIP roles. Rows may override the direct GAMS solver or GDPopt
   role solvers when a model/strategy is known to produce a linear MIP under this
-  profile.
+  profile. The `spectralog`/`gdpopt.enumerate` row is intentionally omitted
+  because the model has 30 independent disjunctions, which implies `2^30`
+  discrete realizations; Pyomo currently materializes that enumeration list
+  before enforcing `time_limit` (Pyomo/pyomo#3953).
 
 Run it with:
 

benchmark_cases/pr58_local.csv

@@ -99,7 +99,6 @@ small_batch,gdpopt.lbb,gams-local,dicopt,3600,local,ipopth,gurobi,dicopt,dicopt
 small_batch,gdpopt.ric,gams-local,dicopt,3600,local,ipopth,gurobi,dicopt,dicopt
 spectralog,gdp.bigm,gams-local,dicopt,3600,local,ipopth,gurobi,dicopt,dicopt
 spectralog,gdp.hull,gams-local,dicopt,3600,local,ipopth,gurobi,dicopt,dicopt
-spectralog,gdpopt.enumerate,gams-local,dicopt,3600,local,ipopth,gurobi,dicopt,dicopt
 spectralog,gdpopt.loa,gams-local,dicopt,3600,local,ipopth,gurobi,dicopt,dicopt
 spectralog,gdpopt.gloa,gams-local,dicopt,3600,local,ipopth,gurobi,dicopt,dicopt
 spectralog,gdpopt.lbb,gams-local,dicopt,3600,local,ipopth,gurobi,dicopt,dicopt

gdplib/spectralog/spectralog.py

@@ -11,6 +11,7 @@
 """
 
 from pyomo.environ import *
+from pyomo.environ import value
 from pyomo.gdp import *
 from pyomo.core.expr.logical_expr import *
 from pyomo.core.plugins.transform.logical_to_linear import (
@@ -19,6 +20,26 @@
 from six import StringIO
 import pandas as pd
 
+P_UPPER_BOUND = 1000
+
+
+def _val_bounds(m, j):
+    """Bounds for each squared residual term under the source-scale P bounds.
+
+    Assumes the committed data's identity ``R`` matrix (which makes ``val[j]``
+    a sum of squared residuals) and nonnegative ``A``/``C`` data; revisit if
+    ``m.R`` ever gains off-diagonal entries.
+    """
+    # Independent of compound k: the largest modeled signal sum(P * A) over
+    # the P box [0, P_UPPER_BOUND].
+    max_modeled_signal = P_UPPER_BOUND * sum(value(m.A[i, j]) for i in m.wave_number)
+    residual_bound = 0
+    for k in m.compounds:
+        concentration = value(m.C[k, j]) / 100
+        max_residual = max(abs(concentration), abs(concentration - max_modeled_signal))
+        residual_bound += max_residual**2
+    return 0, residual_bound
+
 
 def build_model():
     """
@@ -106,7 +127,9 @@ def build_model():
     )
 
     m.val = Var(
-        m.spectra_data, doc="Calculated objective values for each spectra data point"
+        m.spectra_data,
+        bounds=_val_bounds,
+        doc="Calculated objective values for each spectra data point",
     )
     m.ent = Var(
         m.compounds,
@@ -122,7 +145,7 @@ def build_model():
     m.P = Var(
         m.compounds,
         m.wave_number,
-        bounds=(0, 1000),
+        bounds=(0, P_UPPER_BOUND),
         doc="Continuous variables estimating the concentration level of each compound at each wave number.",
     )
 
@@ -146,7 +169,7 @@ def d(m, k, i):
         """
         return [
             [
-                m.P[k, i] <= 1000,
+                m.P[k, i] <= P_UPPER_BOUND,
                 m.P[k, i] >= 0,
                 m.ent[k, i] == 1,
             ],  # Conditions for the compound being active

tests/test_benchmark.py

@@ -23,6 +23,8 @@
     summarize_warning_rows,
 )
 
+PR58_LOCAL_OMITTED_CASES = {("spectralog", "gdpopt.enumerate")}
+
 
 def test_gams_solve_options_use_requested_timelimit():
     assert _gams_solve_options(123) == [
@@ -195,7 +197,14 @@ def test_committed_pr58_local_cases_cover_default_matrix():
 
     cases = read_benchmark_cases_file(cases_file, defaults)
 
-    assert len(cases) == len(PR58_BENCHMARK_INSTANCES) * len(DEFAULT_STRATEGIES)
+    expected_cases = {
+        (instance, strategy)
+        for instance in PR58_BENCHMARK_INSTANCES
+        for strategy in DEFAULT_STRATEGIES
+    } - PR58_LOCAL_OMITTED_CASES
+
+    assert len(cases) == len(expected_cases)
+    assert {(case.instance, case.strategy) for case in cases} == expected_cases
     assert {case.instance for case in cases} == set(PR58_BENCHMARK_INSTANCES)
     assert {case.strategy for case in cases} == set(DEFAULT_STRATEGIES)
     assert {case.solver_profile for case in cases} == {"gams-local"}

tests/test_spectralog.py

@@ -0,0 +1,38 @@
+import pytest
+import pyomo.environ as pyo
+from pyomo.gdp import Disjunct, Disjunction
+
+from gdplib.spectralog import build_model
+
+EXPECTED_VAL_UPPER_BOUNDS = {
+    1: 157341.7613,
+    2: 807704.4181,
+    3: 511037.7934,
+    4: 476646.6045,
+    5: 802891.9736,
+    6: 569936.876,
+    7: 429417.6493,
+    8: 186408.5249,
+}
+
+
+def test_spectralog_objective_value_variables_have_source_bounds():
+    model = build_model()
+
+    for j, upper_bound in EXPECTED_VAL_UPPER_BOUNDS.items():
+        assert model.val[j].lb == 0
+        assert model.val[j].ub == pytest.approx(upper_bound)
+
+
+@pytest.mark.parametrize("transformation", ["gdp.bigm", "gdp.hull"])
+def test_spectralog_reformulates_with_supported_gdp_transformations(transformation):
+    model = build_model()
+
+    pyo.TransformationFactory(transformation).apply_to(model)
+
+    # Forward guard: the model currently declares no LogicalConstraint (its
+    # logic uses BooleanVar + associate_binary_var), so this only fires if
+    # future edits add one that a transformation fails to deactivate.
+    assert not any(model.component_data_objects(pyo.LogicalConstraint, active=True))
+    assert not any(model.component_data_objects(Disjunction, active=True))
+    assert not any(model.component_data_objects(Disjunct, active=True))