Feature/example2

pyproject.toml

@@ -34,11 +34,12 @@ classifiers = [
 ]
 dependencies = [
     "boltons>=25.0.0",
+    "colorama>=0.4.6",
     "cvxpy>=1.7.3",
     "joblib>=1.5.2",
     "matplotlib>=3.10.7",
     "mrg32k3a[rust]>=2.0.0",
-    "numpy>=2.3.4",
+    "numpy>=2.3.4,<2.4",
     "pandas>=2.3.3",
     "pillow>=12.0.0",
     "pydantic>=2.12.3",

scripts/generate_experiment_results.py

@@ -5,13 +5,16 @@
 from pathlib import Path
 
 import zstandard as zstd
+from colorama import Fore, init
 
 # Append the parent directory (simopt package) to the system path
 sys.path.append(str(Path(__file__).resolve().parent.parent))
 
 import simopt.directory as directory
 from simopt.experiment_base import ProblemSolver, post_normalize
 
+init(autoreset=True)
+
 # Workaround for AutoAPI
 problem_directory = directory.problem_directory
 solver_directory = directory.solver_directory
@@ -24,6 +27,10 @@
 EXPECTED_RESULTS_DIR = HOME_DIR / "test" / "expected_results"
 
 
+def _color_text(text: str, color: str | int) -> str:
+    return color + text  # type: ignore
+
+
 # Based off the similar function in simopt/experiment_base.py
 def is_compatible(problem_name: str, solver_name: str) -> bool:
     """Check if a solver is compatible with a problem.
@@ -100,6 +107,7 @@ def create_test(problem_name: str, solver_name: str) -> None:
 
 def main() -> None:
     """Create test cases for all compatible problem-solver pairs."""
+    skip_problems = {"ERM-EXAMPLE-1"}
     # Create a list of compatible problem-solver pairs
     compatible_pairs = [
         (problem_name, solver_name)
@@ -123,7 +131,10 @@ def main() -> None:
         results_filename = f"{file_problem_name}_{file_solver_name}.pickle.zst"
         # If file exists, skip it
         if results_filename in existing_results:
-            print(f"Test for {pair} already exists")
+            print(_color_text(f"Test for {pair} already exists", Fore.GREEN))
+            continue
+        if problem_name in skip_problems:
+            print(_color_text(f"Skipping test for {pair}", Fore.YELLOW))
             continue
         # If file doesn't exist, create it
         print(f"Creating test for {pair}")

simopt/bootstrap.py

@@ -321,41 +321,45 @@ def functional_of_curves(
                 for curves in solver_1_curves
             ]
         ),
-        PlotType.DIFFERENCE_OF_CDF_SOLVABILITY: lambda: curve_utils.difference_of_curves(  # noqa: E501
-            curve_utils.mean_of_curves(
-                [
-                    curve_utils.cdf_of_curves_crossing_times(
-                        curves, threshold=solve_tol
-                    )
-                    for curves in solver_1_curves
-                ]
-            ),
-            curve_utils.mean_of_curves(
-                [
-                    curve_utils.cdf_of_curves_crossing_times(
-                        curves, threshold=solve_tol
-                    )
-                    for curves in solver_2_curves  # type: ignore
-                ]
-            ),
+        PlotType.DIFFERENCE_OF_CDF_SOLVABILITY: lambda: (
+            curve_utils.difference_of_curves(
+                curve_utils.mean_of_curves(
+                    [
+                        curve_utils.cdf_of_curves_crossing_times(
+                            curves, threshold=solve_tol
+                        )
+                        for curves in solver_1_curves
+                    ]
+                ),
+                curve_utils.mean_of_curves(
+                    [
+                        curve_utils.cdf_of_curves_crossing_times(
+                            curves, threshold=solve_tol
+                        )
+                        for curves in solver_2_curves  # type: ignore
+                    ]
+                ),
+            )
         ),
-        PlotType.DIFFERENCE_OF_QUANTILE_SOLVABILITY: lambda: curve_utils.difference_of_curves(  # noqa: E501
-            curve_utils.mean_of_curves(
-                [
-                    curve_utils.quantile_cross_jump(
-                        curves, threshold=solve_tol, beta=beta
-                    )
-                    for curves in solver_1_curves
-                ]
-            ),
-            curve_utils.mean_of_curves(
-                [
-                    curve_utils.quantile_cross_jump(
-                        curves, threshold=solve_tol, beta=beta
-                    )
-                    for curves in solver_2_curves  # type: ignore
-                ]
-            ),
+        PlotType.DIFFERENCE_OF_QUANTILE_SOLVABILITY: lambda: (
+            curve_utils.difference_of_curves(
+                curve_utils.mean_of_curves(
+                    [
+                        curve_utils.quantile_cross_jump(
+                            curves, threshold=solve_tol, beta=beta
+                        )
+                        for curves in solver_1_curves
+                    ]
+                ),
+                curve_utils.mean_of_curves(
+                    [
+                        curve_utils.quantile_cross_jump(
+                            curves, threshold=solve_tol, beta=beta
+                        )
+                        for curves in solver_2_curves  # type: ignore
+                    ]
+                ),
+            )
         ),
         PlotType.MEAN_FEASIBILITY_PROGRESS: lambda: curve_utils.mean_of_curves(
             single_curves

simopt/experiment/run_solver.py

@@ -2,6 +2,7 @@
 
 import logging
 import time
+from copy import deepcopy
 
 import pandas as pd
 from joblib import Parallel, delayed
@@ -111,13 +112,17 @@ def run_solver(
     logging.info(f"Running solver {solver.name} on problem {problem.name}.")
     logging.debug("Starting macroreplications")
 
+    # TODO: Long-term fix is to make Solver/Problem immutable (stateless)
+    # so macroreps can share instances safely without deepcopy.
     if n_jobs == 1:
         results: list[tuple] = [
-            _run_mrep(solver, problem, i) for i in range(n_macroreps)
+            _run_mrep(deepcopy(solver), deepcopy(problem), i)
+            for i in range(n_macroreps)
         ]
     else:
         results: list[tuple] = Parallel(n_jobs=n_jobs)(
-            delayed(_run_mrep)(solver, problem, i) for i in range(n_macroreps)
+            delayed(_run_mrep)(deepcopy(solver), deepcopy(problem), i)
+            for i in range(n_macroreps)
         )
 
     dfs = []

simopt/models/amusementpark.py

@@ -245,8 +245,9 @@ class AmusementParkMinDepartConfig(BaseModel):
     initial_solution: Annotated[
         tuple[int, ...],
         Field(
-            default_factory=lambda: (PARK_CAPACITY - NUM_ATTRACTIONS + 1,)
-            + (1,) * (NUM_ATTRACTIONS - 1),
+            default_factory=lambda: (
+                (PARK_CAPACITY - NUM_ATTRACTIONS + 1,) + (1,) * (NUM_ATTRACTIONS - 1)
+            ),
             description="Initial solution from which solvers start.",
         ),
     ]

simopt/models/example.py

@@ -167,3 +167,126 @@
                 factorized=False,
             )
         )
+
+
+class Example2ModelConfig(BaseModel):
+    """Configuration model for Example-2 simulation.
+
+    A model that is a deterministic quadratic function evaluated with noise.
+    """
+
+    x: Annotated[
+        tuple[int, ...],
+        Field(
+            default=(0, 0, 0, 0),
+            description="point to evaluate",
+        ),
+    ]
+
+
+class Example2ProblemConfig(BaseModel):
+    """Configuration model for Example-2 Problem.
+
+    Base class to implement simulation-optimization problems.
+    """
+
+    initial_solution: Annotated[
+        tuple[int, ...],
+        Field(
+            default=(0, 0, 0, 0),
+            description="initial solution",
+        ),
+    ]
+    budget: Annotated[
+        int,
+        Field(
+            default=1000,
+            description="max # of replications for a solver to take",
+            gt=0,
+            json_schema_extra={"isDatafarmable": False},
+        ),
+    ]
+
+
+class Example2Model(Model):
+    """A model that is a deterministic quadratic function evaluated with noise."""
+
+    class_name_abbr: ClassVar[str] = "EXAMPLE-2-MODEL"
+    class_name: ClassVar[str] = "Quadratic Function + Noise (Discrete)"
+    config_class: ClassVar[type[BaseModel]] = Example2ModelConfig
+    n_rngs: ClassVar[int] = 1
+    n_responses: ClassVar[int] = 1
+
+    def __init__(self, fixed_factors: dict | None = None) -> None:
+        """Initialize the model.
+
+        Args:
+            fixed_factors (dict | None): fixed factors of the model.
+                If None, use default values.
+        """
+        # Let the base class handle default arguments.
+        super().__init__(fixed_factors)
+        self.noise_model = Normal()
+
+    def before_replicate(self, rng_list: list[MRG32k3a]) -> None:  # noqa: D102
+        self.noise_model.set_rng(rng_list[0])
+
+    def replicate(self) -> tuple[dict, dict]:
+        """Evaluate a quadratic function f(x) with stochastic noise."""
+        x = np.array(self.factors["x"])
+        target = np.array([1, 2, 3, 4])
+        fn_eval_at_x = np.sum((x - target) ** 2) + self.noise_model.random()
+
+        responses = {"est_f(x)": fn_eval_at_x}
+        return responses, {}
+
+
+class Example2Problem(Problem):
+    """Discrete quadratic minimization example with noise."""
+
+    class_name_abbr: ClassVar[str] = "EXAMPLE-2"
+    class_name: ClassVar[str] = "Min Quadratic Function + Noise (Discrete)"
+    config_class: ClassVar[type[BaseModel]] = Example2ProblemConfig
+    model_class: ClassVar[type[Model]] = Example2Model
+    n_objectives: ClassVar[int] = 1
+    n_stochastic_constraints: ClassVar[int] = 0
+    minmax: ClassVar[tuple[int, ...]] = (-1,)
+    constraint_type: ClassVar[ConstraintType] = ConstraintType.UNCONSTRAINED
+    variable_type: ClassVar[VariableType] = VariableType.DISCRETE
+    gradient_available: ClassVar[bool] = False
+    model_default_factors: ClassVar[dict] = {}
+    model_decision_factors: ClassVar[set[str]] = {"x"}
+
+    @property
+    def optimal_value(self) -> float | None:  # noqa: D102
+        return 0.0
+
+    @property
+    def optimal_solution(self) -> tuple | None:  # noqa: D102
+        return (1, 2, 3, 4)
+
+    @property
+    def dim(self) -> int:  # noqa: D102
+        return 4
+
+    @property
+    def lower_bounds(self) -> tuple:  # noqa: D102
+        return (-4,) * self.dim
+
+    @property
+    def upper_bounds(self) -> tuple:  # noqa: D102
+        return (4,) * self.dim
+
+    def vector_to_factor_dict(self, vector: tuple) -> dict:  # noqa: D102
+        return {"x": vector[:]}
+
+    def factor_dict_to_vector(self, factor_dict: dict) -> tuple:  # noqa: D102
+        return tuple(factor_dict["x"])
+
+    def replicate(self, _x: tuple) -> RepResult:  # noqa: D102
+        responses, _ = self.model.replicate()
+        objectives = [Objective(stochastic=responses["est_f(x)"])]
+        return RepResult(objectives=objectives)
+
+    def get_random_solution(self, rand_sol_rng: MRG32k3a) -> tuple:  # noqa: D102
+        return tuple(rand_sol_rng.randint(-4, 4) for _ in range(self.dim))