diff --git a/pyproject.toml b/pyproject.toml index a1e4bfed..d49ccccc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -14,6 +14,12 @@ exclude = ["build*"] [tool.uv.sources] simopt-extensions = { git = "https://github.com/cenwangumass/simopt-extensions", rev = "b40a9c8" } +simopt-dsl = { workspace = true } + +[tool.uv.workspace] +members = [ + "simopt-dsl", +] [project] name = "simoptlib" @@ -50,6 +56,7 @@ dependencies = [ "scipy>=1.16.3", "seaborn>=0.13.2", "simpy>=4.1.2", + "simopt-dsl>=0.1.0", "structlog>=25.5.0", ] diff --git a/simopt-dsl/.gitignore b/simopt-dsl/.gitignore new file mode 100644 index 00000000..5fed5c72 --- /dev/null +++ b/simopt-dsl/.gitignore @@ -0,0 +1,174 @@ +### Python ### +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# poetry +# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control. +# This is especially recommended for binary packages to ensure reproducibility, and is more +# commonly ignored for libraries. +# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control +#poetry.lock + +# pdm +# Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control. +#pdm.lock +# pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it +# in version control. +# https://pdm.fming.dev/#use-with-ide +.pdm.toml + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + +# PyCharm +# JetBrains specific template is maintained in a separate JetBrains.gitignore that can +# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore +# and can be added to the global gitignore or merged into this file. For a more nuclear +# option (not recommended) you can uncomment the following to ignore the entire idea folder. +#.idea/ + +### Python Patch ### +# Poetry local configuration file - https://python-poetry.org/docs/configuration/#local-configuration +poetry.toml + +# ruff +.ruff_cache/ + +# LSP config files +pyrightconfig.json + +# uv +uv.lock diff --git a/simopt-dsl/.python-version b/simopt-dsl/.python-version new file mode 100644 index 00000000..24ee5b1b --- /dev/null +++ b/simopt-dsl/.python-version @@ -0,0 +1 @@ +3.13 diff --git a/simopt-dsl/README.md b/simopt-dsl/README.md new file mode 100644 index 00000000..e69de29b diff --git a/simopt-dsl/pyproject.toml b/simopt-dsl/pyproject.toml new file mode 100644 index 00000000..9d635e03 --- /dev/null +++ b/simopt-dsl/pyproject.toml @@ -0,0 +1,56 @@ +[project] +name = "simopt-dsl" +version = "0.1.0" +description = "A declarative modeling language for simulation optimization." +readme = "README.md" +authors = [ + { name = "Cen Wang", email = "cenwang@umass.edu" } +] +requires-python = ">=3.11, <3.14" +dependencies = [ + "numpy>=2.0.0", + "sympy>=1.14.0", +] + +[build-system] +requires = ["uv_build>=0.12.4,<0.13.0"] +build-backend = "uv_build" + +[dependency-groups] +dev = [ + "pytest>=9.1.1", +] + +[tool.ruff] +line-length = 100 + +[tool.ruff.format] +skip-magic-trailing-comma = true + +[tool.ruff.lint] +isort.split-on-trailing-comma = false + +select = [ + "E", # pycodestyle (error) + "F", # pyflakes + "B", # bugbear + "B9", + "C4", # flake8-comprehensions + "SIM", # flake8-simplify + "I", # isort + "UP", # pyupgrade + "PIE", # flake8-pie + "PGH", # pygrep-hooks + "PYI", # flake8-pyi + "RUF", + "S602", # flake8-bandit: subprocess-popen-with-shell-equals-true +] + +ignore = [ + # only relevant if you run a script with `python -0`, + # which seems unlikely for any of the scripts in this repo + "B011", + # Leave it to the formatter to split long lines and + # the judgement of all of us. + "E501" +] diff --git a/simopt-dsl/src/simopt_dsl/__init__.py b/simopt-dsl/src/simopt_dsl/__init__.py new file mode 100644 index 00000000..75bd59e9 --- /dev/null +++ b/simopt-dsl/src/simopt_dsl/__init__.py @@ -0,0 +1,8 @@ +"""Declarative modeling primitives for simulation optimization.""" + +from simopt_dsl.expressions import mean, sum +from simopt_dsl.model import Model +from simopt_dsl.simulation import Simulation +from simopt_dsl.variables import Variable, VectorVariable + +__all__ = ["Model", "Simulation", "Variable", "VectorVariable", "mean", "sum"] diff --git a/simopt-dsl/src/simopt_dsl/expressions.py b/simopt-dsl/src/simopt_dsl/expressions.py new file mode 100644 index 00000000..6a356488 --- /dev/null +++ b/simopt-dsl/src/simopt_dsl/expressions.py @@ -0,0 +1,201 @@ +"""Expressions used to declare objectives and constraints.""" + +from __future__ import annotations + +from collections.abc import Iterable +from dataclasses import dataclass +from typing import Any + + +class Expression: + """Base class for scalar, replication-level expressions.""" + + def evaluate(self, context: EvaluationContext) -> float: + raise NotImplementedError + + def __add__(self, other: object) -> Expression | AggregateExpression: + return _binary_expression("+", self, other) + + def __radd__(self, other: object) -> Expression | AggregateExpression: + return _binary_expression("+", other, self) + + def __sub__(self, other: object) -> Expression | AggregateExpression: + return _binary_expression("-", self, other) + + def __rsub__(self, other: object) -> Expression | AggregateExpression: + return _binary_expression("-", other, self) + + def __mul__(self, other: object) -> Expression | AggregateExpression: + return _binary_expression("*", self, other) + + def __rmul__(self, other: object) -> Expression | AggregateExpression: + return _binary_expression("*", other, self) + + def __truediv__(self, other: object) -> Expression | AggregateExpression: + return _binary_expression("/", self, other) + + def __rtruediv__(self, other: object) -> Expression | AggregateExpression: + return _binary_expression("/", other, self) + + def __neg__(self) -> Expression: + return BinaryExpression("*", Constant(-1.0), self) + + def __le__(self, other: object) -> Constraint: + return Constraint(self, "<=", as_expression(other)) + + def __ge__(self, other: object) -> Constraint: + return Constraint(self, ">=", as_expression(other)) + + +@dataclass(frozen=True) +class Constant(Expression): + value: float + + def evaluate(self, context: EvaluationContext) -> float: + return self.value + + +@dataclass(frozen=True) +class BinaryExpression(Expression): + operator: str + left: Expression + right: Expression + + def evaluate(self, context: EvaluationContext) -> float: + return apply_binary_operator( + self.operator, self.left.evaluate(context), self.right.evaluate(context) + ) + + +@dataclass(frozen=True) +class Constraint: + left: Expression + sense: str + right: Expression + + def satisfied(self, context: EvaluationContext) -> bool: + left = self.left.evaluate(context) + right = self.right.evaluate(context) + if self.sense == "<=": + return left <= right + if self.sense == ">=": + return left >= right + raise ValueError(f"unknown constraint sense {self.sense!r}") + + def residual(self) -> Expression: + """Return an expression whose feasible values are nonpositive.""" + if self.sense == "<=": + return BinaryExpression("-", self.left, self.right) + if self.sense == ">=": + return BinaryExpression("-", self.right, self.left) + raise ValueError(f"unknown constraint sense {self.sense!r}") + + +class AggregateExpression: + """Base class for expressions estimated over simulation replications.""" + + def __add__(self, other: object) -> AggregateExpression: + return BinaryAggregateExpression("+", self, as_aggregate_expression(other)) + + def __radd__(self, other: object) -> AggregateExpression: + return BinaryAggregateExpression("+", as_aggregate_expression(other), self) + + def __sub__(self, other: object) -> AggregateExpression: + return BinaryAggregateExpression("-", self, as_aggregate_expression(other)) + + def __rsub__(self, other: object) -> AggregateExpression: + return BinaryAggregateExpression("-", as_aggregate_expression(other), self) + + def __mul__(self, other: object) -> AggregateExpression: + return BinaryAggregateExpression("*", self, as_aggregate_expression(other)) + + def __rmul__(self, other: object) -> AggregateExpression: + return BinaryAggregateExpression("*", as_aggregate_expression(other), self) + + def __truediv__(self, other: object) -> AggregateExpression: + return BinaryAggregateExpression("/", self, as_aggregate_expression(other)) + + def __rtruediv__(self, other: object) -> AggregateExpression: + return BinaryAggregateExpression("/", as_aggregate_expression(other), self) + + def __neg__(self) -> AggregateExpression: + return BinaryAggregateExpression("*", Constant(-1.0), self) + + +@dataclass(frozen=True) +class Mean(AggregateExpression): + expression: Expression + + +@dataclass(frozen=True) +class BinaryAggregateExpression(AggregateExpression): + operator: str + left: Expression | AggregateExpression + right: Expression | AggregateExpression + + +@dataclass +class EvaluationContext: + variables: dict[Any, float] + metrics: dict[tuple[str, str], Any] + metric_derivatives: dict[tuple[str, str, tuple[int, ...], str], float] + + def __init__(self, variables: dict[Any, float]) -> None: + self.variables = variables + self.metrics = {} + self.metric_derivatives = {} + + +def mean(expression: object) -> AggregateExpression: + """Return the replication mean of a scalar expression.""" + return Mean(as_expression(expression)) + + +def sum(expressions: Iterable[object]) -> Expression: + """Return the sum of scalar expressions, or zero when empty.""" + iterator = iter(expressions) + total = as_expression(next(iterator, 0.0)) + for expression in iterator: + total = BinaryExpression("+", total, as_expression(expression)) + return total + + +def as_expression(value: object) -> Expression: + if isinstance(value, Expression): + return value + if isinstance(value, (int, float)): + return Constant(float(value)) + raise TypeError(f"expected an expression or number, got {type(value).__name__}") + + +def as_aggregate_expression(value: object) -> Expression | AggregateExpression: + if isinstance(value, (Expression, AggregateExpression)): + return value + if isinstance(value, (int, float)): + return Constant(float(value)) + raise TypeError(f"expected an expression, statistic, or number, got {type(value).__name__}") + + +def apply_binary_operator(operator: str, left: float, right: float) -> float: + if operator == "+": + return left + right + if operator == "-": + return left - right + if operator == "*": + return left * right + if operator == "/": + return left / right + raise ValueError(f"unknown operator {operator!r}") + + +def _binary_expression( + operator: str, left: object, right: object +) -> Expression | AggregateExpression: + if isinstance(left, AggregateExpression) or isinstance(right, AggregateExpression): + return BinaryAggregateExpression( + operator, as_aggregate_expression(left), as_aggregate_expression(right) + ) + return BinaryExpression(operator, as_expression(left), as_expression(right)) + + +__all__ = ["mean", "sum"] diff --git a/simopt-dsl/src/simopt_dsl/model.py b/simopt-dsl/src/simopt_dsl/model.py new file mode 100644 index 00000000..e6f94731 --- /dev/null +++ b/simopt-dsl/src/simopt_dsl/model.py @@ -0,0 +1,552 @@ +"""Simulation-optimization model declarations and replication evaluation.""" + +from __future__ import annotations + +from collections.abc import Callable, Iterable, Mapping, Sequence +from dataclasses import dataclass +from math import isinf +from numbers import Real +from typing import Any, cast + +import sympy as sp + +from simopt_dsl.expressions import ( + AggregateExpression, + BinaryAggregateExpression, + BinaryExpression, + Constant, + Constraint, + EvaluationContext, + Expression, + Mean, + apply_binary_operator, + mean, +) +from simopt_dsl.simulation import Metric, Simulation, SimulationResult, component_items +from simopt_dsl.variables import DecisionVariable, Variable, VectorVariable, components + +Number = int | float + + +@dataclass(frozen=True) +class StochasticConstraintEvaluation: + value: float + gradient: tuple[float, ...] | None + + +@dataclass(frozen=True) +class ReplicationEvaluation: + objective: float + objective_gradient: tuple[float, ...] | None + stochastic_constraints: tuple[StochasticConstraintEvaluation, ...] = () + + +class Model: + """A declarative simulation-optimization model.""" + + def __init__(self, name: str = "") -> None: + if not isinstance(name, str): + raise TypeError("model name must be a string") + self.name = name + self.variables: list[Variable] = [] + self.simulations: list[Simulation] = [] + self.deterministic_constraints: list[tuple[str, Constraint]] = [] + self.stochastic_constraints: list[tuple[str, Constraint]] = [] + self.objective: AggregateExpression | None = None + self.objective_sense = "minimize" + self.n_rngs = 1 + self._declared_variables: list[DecisionVariable] = [] + self._evaluation_plan: _EvaluationPlan | None = None + + def __getstate__(self) -> dict[str, Any]: + """Return serializable state without process-local compiled callables.""" + state = self.__dict__.copy() + state["_evaluation_plan"] = None + return state + + def add_continuous_variable( + self, lb: Number, ub: Number, name: str | None = None, initial: Number | None = None + ) -> Variable: + """Add a scalar continuous decision variable.""" + return self._add_scalar_variable(lb, ub, name, initial, integer=False) + + def add_continuous_vector( + self, + lb: Number | Iterable[Number], + ub: Number | Iterable[Number], + name: str | None = None, + shape: int | tuple[int] | None = None, + initial: Number | Iterable[Number] | None = None, + ) -> VectorVariable: + """Add a one-dimensional continuous decision variable.""" + return self._add_vector_variable(lb, ub, name, shape, initial, integer=False) + + def add_integer_variable( + self, lb: Number, ub: Number, name: str | None = None, initial: Number | None = None + ) -> Variable: + """Add a scalar integer decision variable.""" + return self._add_scalar_variable(lb, ub, name, initial, integer=True) + + def add_integer_vector( + self, + lb: Number | Iterable[Number], + ub: Number | Iterable[Number], + name: str | None = None, + shape: int | tuple[int] | None = None, + initial: Number | Iterable[Number] | None = None, + ) -> VectorVariable: + """Add a one-dimensional integer decision variable.""" + return self._add_vector_variable(lb, ub, name, shape, initial, integer=True) + + def add_linear_constraint(self, constraint: Constraint, name: str = "") -> None: + """Add a deterministic constraint.""" + self._validate_constraint(constraint, name) + self.deterministic_constraints.append((name, constraint)) + self._invalidate_evaluation_plan() + + def add_stochastic_constraint(self, constraint: Constraint, name: str = "") -> None: + """Add a replication-level stochastic constraint.""" + self._validate_constraint(constraint, name) + self.stochastic_constraints.append((name, constraint)) + self._invalidate_evaluation_plan() + + def add_simulation( + self, + name: str | None = None, + run: Callable[..., SimulationResult] | None = None, + decisions: Mapping[str, DecisionVariable] | None = None, + n_rngs: int = 1, + ) -> Simulation: + """Add a simulation callback and bind its decision variables.""" + if run is None or not callable(run): + raise TypeError("simulation run callable is required") + if decisions is None: + raise TypeError("simulation decisions are required") + if not isinstance(n_rngs, int) or isinstance(n_rngs, bool) or n_rngs <= 0: + raise ValueError("simulation must use at least one RNG") + simulation_name = self._resolve_name(name, "simulation", self.simulations) + + decision_variables: dict[str, DecisionVariable] = {} + seen_components: set[int] = set() + for decision_name, variable in decisions.items(): + if not isinstance(decision_name, str): + raise TypeError("simulation decision names must be strings") + if not decision_name: + raise ValueError("simulation decision names cannot be empty") + if not isinstance(variable, (Variable, VectorVariable)): + raise TypeError("simulation decisions must be decision variables") + if not any(variable is declared for declared in self._declared_variables): + raise ValueError(f"simulation decision {decision_name!r} is not in this model") + variable_components = components(variable) + if any(id(component) in seen_components for component in variable_components): + raise ValueError("simulation decision variable components must be unique") + decision_variables[decision_name] = variable + seen_components.update(id(component) for component in variable_components) + + simulation = Simulation(simulation_name, run, decision_variables, n_rngs) + self.simulations.append(simulation) + self.n_rngs = max(self.n_rngs, n_rngs) + self._invalidate_evaluation_plan() + return simulation + + def maximize(self, objective: AggregateExpression | Expression) -> None: + """Set the model's maximization objective.""" + self._set_objective(objective, "maximize") + + def minimize(self, objective: AggregateExpression | Expression) -> None: + """Set the model's minimization objective.""" + self._set_objective(objective, "minimize") + + def initial_vector(self) -> tuple[float, ...]: + """Return initial values in flattened solver order.""" + return tuple(variable.initial for variable in self.variables) + + def lower_bounds(self) -> tuple[float, ...]: + """Return lower bounds in flattened solver order.""" + return tuple(variable.lb for variable in self.variables) + + def upper_bounds(self) -> tuple[float, ...]: + """Return upper bounds in flattened solver order.""" + return tuple(variable.ub for variable in self.variables) + + def unpack_vector(self, values: Iterable[float]) -> dict[str, float | tuple[float, ...]]: + """Regroup a flat solver vector by declared scalar and vector variables.""" + flat_values = tuple(float(value) for value in values) + if len(flat_values) != len(self.variables): + raise ValueError("solution dimension does not match model variables") + + unpacked: dict[str, float | tuple[float, ...]] = {} + offset = 0 + for variable in self._declared_variables: + next_offset = offset + len(components(variable)) + component_values = flat_values[offset:next_offset] + unpacked[variable.name] = ( + component_values[0] if isinstance(variable, Variable) else component_values + ) + offset = next_offset + return unpacked + + def run_replication( + self, values: Iterable[float], rngs: Sequence[Any] + ) -> ReplicationEvaluation: + """Evaluate one replication at a solution using caller-owned RNGs.""" + plan = self._get_evaluation_plan() + flat_values = tuple(float(value) for value in values) + if len(flat_values) != len(self.variables): + raise ValueError("solution dimension does not match model variables") + for variable, value in zip(self.variables, flat_values, strict=True): + if not variable.lb <= value <= variable.ub: + raise ValueError("solution variable value must be within its bounds") + if variable.integer: + _validate_integer(value, "integer variable value") + if len(rngs) < self.n_rngs: + raise ValueError("not enough RNGs for model simulations") + + context = EvaluationContext(dict(zip(self.variables, flat_values, strict=True))) + for _, constraint in self.deterministic_constraints: + if not constraint.satisfied(context): + raise ValueError("solution violates a deterministic constraint") + + for simulation in self.simulations: + evaluation = simulation.evaluate(context, rngs[: simulation.n_rngs]) + for metric_name, metric_value in evaluation.metrics.items(): + context.metrics[(simulation.name, metric_name)] = metric_value + for key, derivative in evaluation.derivatives.items(): + metric_name, metric_indices, decision_name = key + context.metric_derivatives[ + (simulation.name, metric_name, metric_indices, decision_name) + ] = derivative + + objective_value = _evaluate_aggregate(plan.objective, context) + objective_gradient = plan.estimate_gradient(plan.objective, context) + stochastic_constraints = tuple( + StochasticConstraintEvaluation( + residual.evaluate(context), plan.estimate_gradient(residual, context) + ) + for residual in plan.stochastic_constraint_residuals + ) + return ReplicationEvaluation(objective_value, objective_gradient, stochastic_constraints) + + def _add_scalar_variable( + self, lb: Number, ub: Number, name: str | None, initial: Number | None, *, integer: bool + ) -> Variable: + lower = _coerce_number(lb, "variable lower bound") + upper = _coerce_number(ub, "variable upper bound") + if integer: + _validate_integer_bound(lower, "integer variable lower bound") + _validate_integer_bound(upper, "integer variable upper bound") + if lower >= upper: + raise ValueError("variable lower bound must be less than upper bound") + variable_name = self._resolve_name(name, "variable", self._declared_variables) + initial_value = lower if initial is None else _coerce_number(initial, "initial value") + if integer: + _validate_integer(initial_value, "integer variable initial value") + _validate_initial(initial_value, lower, upper) + variable = Variable(variable_name, lower, upper, initial_value, integer) + self.variables.append(variable) + self._declared_variables.append(variable) + self._invalidate_evaluation_plan() + return variable + + def _add_vector_variable( + self, + lb: Number | Iterable[Number], + ub: Number | Iterable[Number], + name: str | None, + shape: int | tuple[int] | None, + initial: Number | Iterable[Number] | None, + *, + integer: bool, + ) -> VectorVariable: + if shape is None: + raise TypeError("vector variable shape is required") + size = _vector_size(shape) + lower = _vector_values(lb, size, "lower bound") + upper = _vector_values(ub, size, "upper bound") + initial_values = ( + lower if initial is None else _vector_values(initial, size, "initial value") + ) + if integer: + for lower_value in lower: + _validate_integer_bound(lower_value, "integer variable lower bound") + for upper_value in upper: + _validate_integer_bound(upper_value, "integer variable upper bound") + for initial_value in initial_values: + _validate_integer(initial_value, "integer variable initial value") + if any( + lower_value >= upper_value + for lower_value, upper_value in zip(lower, upper, strict=True) + ): + raise ValueError("each vector variable lower bound must be less than its upper bound") + for initial_value, lower_value, upper_value in zip( + initial_values, lower, upper, strict=True + ): + _validate_initial(initial_value, lower_value, upper_value) + + variable_name = self._resolve_name(name, "variable", self._declared_variables) + scalar_components = tuple( + Variable( + f"{variable_name}[{index}]", + lower[index], + upper[index], + initial_values[index], + integer, + ) + for index in range(size) + ) + variable = VectorVariable(variable_name, scalar_components) + self.variables.extend(scalar_components) + self._declared_variables.append(variable) + self._invalidate_evaluation_plan() + return variable + + def _set_objective(self, objective: AggregateExpression | Expression, sense: str) -> None: + if not isinstance(objective, (AggregateExpression, Expression)): + raise TypeError("objective must be an expression") + self.objective = ( + objective if isinstance(objective, AggregateExpression) else mean(objective) + ) + self.objective_sense = sense + self._invalidate_evaluation_plan() + + def _resolve_name(self, name: str | None, prefix: str, existing: Iterable[Any]) -> str: + existing_names = {item.name for item in existing} + if name is None or name == "": + index = 1 + while f"{prefix}_{index}" in existing_names: + index += 1 + return f"{prefix}_{index}" + if not isinstance(name, str): + raise TypeError(f"{prefix} name must be a string") + if name in existing_names: + raise ValueError(f"duplicate {prefix} name {name!r}") + return name + + @staticmethod + def _validate_constraint(constraint: Constraint, name: str) -> None: + if not isinstance(constraint, Constraint): + raise TypeError("constraint must be created with <= or >=") + if not isinstance(name, str): + raise TypeError("constraint name must be a string") + + def _invalidate_evaluation_plan(self) -> None: + self._evaluation_plan = None + + def _get_evaluation_plan(self) -> _EvaluationPlan: + if self._evaluation_plan is None: + self._evaluation_plan = _EvaluationPlan(self) + return self._evaluation_plan + + +class _UnsupportedGradient(Exception): + pass + + +class _SymbolicGradient: + def __init__( + self, variables: Sequence[Variable], expression: Expression | AggregateExpression + ) -> None: + self.variables = variables + self.variable_symbols = { + variable: sp.Symbol(f"v_{index}", real=True) for index, variable in enumerate(variables) + } + self.metric_applications: dict[tuple[str, str, tuple[int, ...]], sp.Expr] = {} + self.metric_derivative_keys: dict[sp.Expr, tuple[str, str, tuple[int, ...], str]] = {} + + sample_expression = self._expression(expression) + gradient_expressions = tuple( + sp.simplify(sp.diff(sample_expression, self.variable_symbols[variable])) + for variable in variables + ) + arguments = ( + *self.variable_symbols.values(), + *self.metric_applications.values(), + *self.metric_derivative_keys, + ) + self._compiled_gradient: Callable[..., Any] = sp.lambdify( + arguments, gradient_expressions, modules="math", dummify=True + ) + + def estimate(self, context: EvaluationContext) -> tuple[float, ...] | None: + arguments = self._argument_values(context) + if arguments is None: + return None + try: + return tuple(float(value) for value in self._compiled_gradient(*arguments)) + except (TypeError, ValueError, ZeroDivisionError, OverflowError): + return None + + def _argument_values(self, context: EvaluationContext) -> tuple[float, ...] | None: + values = [context.variables[variable] for variable in self.variables] + for simulation_name, metric_name, metric_indices in self.metric_applications: + value = context.metrics[(simulation_name, metric_name)] + for index in metric_indices: + value = value[index] + values.append(float(value)) + for derivative_key in self.metric_derivative_keys.values(): + if derivative_key not in context.metric_derivatives: + return None + values.append(context.metric_derivatives[derivative_key]) + return tuple(values) + + def _expression(self, expression: Expression | AggregateExpression) -> sp.Expr: + if isinstance(expression, Constant): + return sp.Float(expression.value) + if isinstance(expression, Variable): + return self.variable_symbols[expression] + if isinstance(expression, Metric): + return self._metric_application(expression) + if isinstance(expression, BinaryExpression): + return self._binary_expression(expression.operator, expression.left, expression.right) + if isinstance(expression, Mean): + return self._expression(expression.expression) + if isinstance(expression, BinaryAggregateExpression): + return self._binary_expression(expression.operator, expression.left, expression.right) + raise _UnsupportedGradient( + f"unsupported objective expression {type(expression).__name__!r}" + ) + + def _binary_expression( + self, + operator: str, + left: Expression | AggregateExpression, + right: Expression | AggregateExpression, + ) -> sp.Expr: + left_expression = self._expression(left) + right_expression = self._expression(right) + if operator == "+": + return left_expression + right_expression + if operator == "-": + return left_expression - right_expression + if operator == "*": + return left_expression * right_expression + if operator == "/": + return left_expression / right_expression + raise ValueError(f"unknown operator {operator!r}") + + def _metric_application(self, metric: Metric) -> sp.Expr: + key = (metric.simulation.name, metric.name, metric.indices) + if key not in self.metric_applications: + function: Any = sp.Function(f"metric_{len(self.metric_applications)}") + items = component_items(metric.simulation.decisions) + application = function(*(self.variable_symbols[component] for _, component in items)) + self.metric_applications[key] = application + for component_name, component in items: + derivative = sp.Derivative(application, self.variable_symbols[component]) + self.metric_derivative_keys[derivative] = ( + metric.simulation.name, + metric.name, + metric.indices, + component_name, + ) + return self.metric_applications[key] + + +class _EvaluationPlan: + """Compiled, reusable objective and constraint evaluation structure.""" + + def __init__(self, model: Model) -> None: + if model.objective is None: + raise ValueError("model has no objective") + self.variables = tuple(model.variables) + self.objective = model.objective + self.stochastic_constraint_residuals = tuple( + constraint.residual() for _, constraint in model.stochastic_constraints + ) + self._gradients: dict[ + int, tuple[Expression | AggregateExpression, _SymbolicGradient | None] + ] = {} + self._gradient(self.objective) + for residual in self.stochastic_constraint_residuals: + self._gradient(residual) + + def estimate_gradient( + self, expression: Expression | AggregateExpression, context: EvaluationContext + ) -> tuple[float, ...] | None: + gradient = self._gradient(expression) + return None if gradient is None else gradient.estimate(context) + + def _gradient(self, expression: Expression | AggregateExpression) -> _SymbolicGradient | None: + key = id(expression) + cached = self._gradients.get(key) + if cached is not None and cached[0] is expression: + return cached[1] + try: + gradient = _SymbolicGradient(self.variables, expression) + except _UnsupportedGradient: + gradient = None + self._gradients[key] = (expression, gradient) + return gradient + + +def _evaluate_aggregate( + expression: Expression | AggregateExpression, context: EvaluationContext +) -> float: + if isinstance(expression, Mean): + return expression.expression.evaluate(context) + if isinstance(expression, BinaryAggregateExpression): + return apply_binary_operator( + expression.operator, + _evaluate_aggregate(expression.left, context), + _evaluate_aggregate(expression.right, context), + ) + if isinstance(expression, Expression): + return expression.evaluate(context) + raise TypeError(f"unsupported aggregate objective {type(expression).__name__!r}") + + +def _coerce_number(value: object, parameter_name: str) -> float: + if not isinstance(value, Real) or isinstance(value, bool): + raise TypeError(f"{parameter_name} must be numeric") + return float(value) + + +def _validate_initial(initial: float, lower: float, upper: float) -> None: + if not lower <= initial <= upper: + raise ValueError("variable initial value must be within its bounds") + + +def _validate_integer(value: float, parameter_name: str) -> None: + if value % 1 != 0: + raise ValueError(f"{parameter_name} must be an integer") + + +def _validate_integer_bound(value: float, parameter_name: str) -> None: + if not isinf(value): + _validate_integer(value, parameter_name) + + +def _vector_size(shape: int | tuple[int]) -> int: + if isinstance(shape, bool): + raise TypeError("vector variable size must be a positive integer") + if isinstance(shape, int): + size = shape + else: + dimensions = tuple(shape) + if len(dimensions) != 1: + raise ValueError("vector variable shape must have exactly one dimension") + size = dimensions[0] + if not isinstance(size, int) or isinstance(size, bool) or size <= 0: + raise ValueError("vector variable size must be a positive integer") + return size + + +def _vector_values( + value: Number | Iterable[Number], size: int, parameter_name: str +) -> tuple[float, ...]: + if isinstance(value, Real) and not isinstance(value, bool): + return (float(value),) * size + if isinstance(value, (str, bytes)): + raise TypeError(f"vector variable {parameter_name} must be numeric") + try: + raw_values = tuple(cast(Iterable[Number], value)) + except TypeError as exc: + raise TypeError(f"vector variable {parameter_name} must be numeric") from exc + if len(raw_values) != size: + raise ValueError( + f"vector variable {parameter_name} must have {size} values, got {len(raw_values)}" + ) + return tuple(_coerce_number(item, f"vector variable {parameter_name}") for item in raw_values) + + +__all__ = ["Model"] diff --git a/simopt-dsl/src/simopt_dsl/py.typed b/simopt-dsl/src/simopt_dsl/py.typed new file mode 100644 index 00000000..e69de29b diff --git a/simopt-dsl/src/simopt_dsl/simulation.py b/simopt-dsl/src/simopt_dsl/simulation.py new file mode 100644 index 00000000..57402691 --- /dev/null +++ b/simopt-dsl/src/simopt_dsl/simulation.py @@ -0,0 +1,208 @@ +"""Simulation callbacks and result normalization.""" + +from __future__ import annotations + +from collections.abc import Callable, Mapping, Sequence +from dataclasses import dataclass +from numbers import Real +from typing import Any + +import numpy as np + +from simopt_dsl.expressions import EvaluationContext, Expression +from simopt_dsl.variables import DecisionVariable, Variable, components + +SimulationResult = tuple[Mapping[str, Any], Mapping[str, Mapping[str, Any]]] + + +@dataclass(frozen=True) +class SimulationEvaluation: + metrics: dict[str, Any] + derivatives: dict[tuple[str, tuple[int, ...], str], float] + + +@dataclass(frozen=True) +class Metric(Expression): + simulation: Simulation + name: str + indices: tuple[int, ...] = () + + def evaluate(self, context: EvaluationContext) -> float: + value = context.metrics[(self.simulation.name, self.name)] + for index in self.indices: + value = value[index] + try: + return float(value) + except (TypeError, ValueError) as exc: + suffix = "".join(f"[{index}]" for index in self.indices) + raise TypeError( + f"simulation metric {self.name!r}{suffix} is not scalar; " + "select a scalar component with metric[index]" + ) from exc + + def __getitem__(self, index: int) -> Metric: + if not isinstance(index, int): + raise TypeError("simulation metric indices must be integers") + return Metric(self.simulation, self.name, (*self.indices, index)) + + +@dataclass +class Simulation: + """A simulation callback and its decision-variable bindings.""" + + name: str + run: Callable[..., SimulationResult] + decisions: dict[str, DecisionVariable] + n_rngs: int = 1 + + def metric(self, name: str) -> Metric: + """Return a scalar expression referring to a simulation response.""" + if not isinstance(name, str): + raise TypeError("simulation metric name must be a string") + if not name: + raise ValueError("simulation metric name cannot be empty") + return Metric(self, name) + + def evaluate(self, context: EvaluationContext, rngs: Sequence[Any]) -> SimulationEvaluation: + """Evaluate the callback once and normalize its metrics and derivatives.""" + decisions = {name: variable.evaluate(context) for name, variable in self.decisions.items()} + result = self.run(decisions, rngs) + return _coerce_evaluation(result, self.decisions) + + +def component_items(decisions: Mapping[str, DecisionVariable]) -> tuple[tuple[str, Variable], ...]: + items: list[tuple[str, Variable]] = [] + seen_names: set[str] = set() + for decision_name, decision in decisions.items(): + decision_components = components(decision) + component_names = ( + (decision_name,) + if isinstance(decision, Variable) + else tuple(f"{decision_name}[{index}]" for index in range(len(decision_components))) + ) + for component_name, component in zip(component_names, decision_components, strict=True): + if component_name in seen_names: + raise ValueError( + f"simulation decision component name {component_name!r} is ambiguous" + ) + seen_names.add(component_name) + items.append((component_name, component)) + return tuple(items) + + +def _coerce_evaluation( + result: Any, decisions: Mapping[str, DecisionVariable] +) -> SimulationEvaluation: + if not isinstance(result, tuple) or len(result) != 2: + raise TypeError("simulation must return a (responses, gradients) pair") + raw_metrics, raw_derivatives = result + if not isinstance(raw_metrics, Mapping): + raise TypeError("simulation response payload must be a mapping") + if not isinstance(raw_derivatives, Mapping): + raise TypeError("simulation gradient payload must be a mapping") + + metrics: dict[str, Any] = {} + for name, value in raw_metrics.items(): + if not isinstance(name, str): + raise TypeError("simulation metric names must be strings") + metrics[name] = value + + return SimulationEvaluation(metrics, _coerce_derivatives(raw_derivatives, decisions, metrics)) + + +def _coerce_derivatives( + raw: Any, decisions: Mapping[str, DecisionVariable], metrics: Mapping[str, Any] +) -> dict[tuple[str, tuple[int, ...], str], float]: + if not isinstance(raw, Mapping): + raise TypeError("simulation derivatives must be a mapping") + + derivatives: dict[tuple[str, tuple[int, ...], str], float] = {} + for metric_name, derivative_values in raw.items(): + if not isinstance(metric_name, str): + raise TypeError("gradient response names must be strings") + metric_shape = _metric_shape(metrics, metric_name) + if not isinstance(derivative_values, Mapping): + raise TypeError(f"gradients for response {metric_name!r} must be a mapping") + for decision_name, value in derivative_values.items(): + if not isinstance(decision_name, str): + raise TypeError("gradient decision names must be strings") + if decision_name in decisions: + _store_decision_derivatives( + derivatives, + metric_name, + metric_shape, + decision_name, + decisions[decision_name], + value, + ) + return derivatives + + +def _store_decision_derivatives( + derivatives: dict[tuple[str, tuple[int, ...], str], float], + metric_name: str, + metric_shape: tuple[int, ...], + decision_name: str, + decision: DecisionVariable, + raw_values: Any, +) -> None: + items = component_items({decision_name: decision}) + ordered_values = _flatten_values(raw_values) + indices = _metric_indices(metric_shape) + expected = len(indices) * len(items) + if len(ordered_values) != expected: + raise ValueError( + f"derivative for decision {decision_name!r} must have " + f"{expected} component(s), got {len(ordered_values)}" + ) + offset = 0 + for metric_index in indices: + for component_name, _ in items: + derivatives[(metric_name, metric_index, component_name)] = ordered_values[offset] + offset += 1 + + +def _metric_shape(metrics: Mapping[str, Any], metric_name: str) -> tuple[int, ...]: + if metric_name not in metrics: + raise ValueError(f"derivative provided for unknown metric {metric_name!r}") + value = metrics[metric_name] + if _is_scalar(value): + return () + try: + array = np.asarray(value) + except (TypeError, ValueError) as exc: + raise TypeError(f"simulation metric {metric_name!r} must be scalar or array-like") from exc + if array.dtype == object: + raise TypeError(f"simulation metric {metric_name!r} must have a rectangular shape") + return tuple(int(size) for size in array.shape) + + +def _metric_indices(shape: tuple[int, ...]) -> tuple[tuple[int, ...], ...]: + if not shape: + return ((),) + return tuple(tuple(int(index) for index in indices) for indices in np.ndindex(shape)) + + +def _is_scalar(value: Any) -> bool: + if isinstance(value, Real): + return True + try: + float(value) + except (TypeError, ValueError): + return False + return True + + +def _flatten_values(raw: Any) -> tuple[float, ...]: + if _is_scalar(raw): + return (float(raw),) + if isinstance(raw, (str, bytes)): + raise TypeError("derivative values must contain ordered numeric values") + try: + values = tuple(raw) + except TypeError as exc: + raise TypeError("derivative values must contain ordered numeric values") from exc + return tuple(item for value in values for item in _flatten_values(value)) + + +__all__ = ["Simulation"] diff --git a/simopt-dsl/src/simopt_dsl/variables.py b/simopt-dsl/src/simopt_dsl/variables.py new file mode 100644 index 00000000..bb18cdfe --- /dev/null +++ b/simopt-dsl/src/simopt_dsl/variables.py @@ -0,0 +1,77 @@ +"""Scalar and vector decision variables.""" + +from __future__ import annotations + +from collections.abc import Iterator +from dataclasses import dataclass, field +from typing import overload +from uuid import UUID, uuid4 + +from simopt_dsl.expressions import EvaluationContext, Expression + + +@dataclass(frozen=True) +class Variable(Expression): + """A scalar decision variable.""" + + name: str + lb: float + ub: float + initial: float + integer: bool = False + _id: UUID = field(default_factory=uuid4, init=False, repr=False) + + def evaluate(self, context: EvaluationContext) -> float: + return context.variables[self] + + +@dataclass(frozen=True) +class VectorVariable: + """A one-dimensional decision variable with scalar expression components.""" + + name: str + components: tuple[Variable, ...] + + @property + def shape(self) -> tuple[int]: + return (len(self.components),) + + @property + def lb(self) -> tuple[float, ...]: + return tuple(component.lb for component in self.components) + + @property + def ub(self) -> tuple[float, ...]: + return tuple(component.ub for component in self.components) + + @property + def initial(self) -> tuple[float, ...]: + return tuple(component.initial for component in self.components) + + def evaluate(self, context: EvaluationContext) -> tuple[float, ...]: + return tuple(component.evaluate(context) for component in self.components) + + def __len__(self) -> int: + return len(self.components) + + def __iter__(self) -> Iterator[Variable]: + return iter(self.components) + + @overload + def __getitem__(self, index: int) -> Variable: ... + + @overload + def __getitem__(self, index: slice) -> tuple[Variable, ...]: ... + + def __getitem__(self, index: int | slice) -> Variable | tuple[Variable, ...]: + return self.components[index] + + +DecisionVariable = Variable | VectorVariable + + +def components(decision: DecisionVariable) -> tuple[Variable, ...]: + return (decision,) if isinstance(decision, Variable) else decision.components + + +__all__ = ["Variable", "VectorVariable"] diff --git a/simopt-dsl/tests/test_model.py b/simopt-dsl/tests/test_model.py new file mode 100644 index 00000000..6cf1e019 --- /dev/null +++ b/simopt-dsl/tests/test_model.py @@ -0,0 +1,17 @@ +"""Tests for model declarations and replication evaluation.""" + +from simopt_dsl import Model + + +def test_vector_component_does_not_alias_scalar_with_same_name() -> None: + model = Model() + vector = model.add_continuous_vector(0, 10, name="x", shape=1) + scalar = model.add_continuous_variable(0, 10, name="x[0]") + model.minimize(vector[0]) + + assert vector[0] != scalar + + evaluation = model.run_replication((2, 7), [object()]) + + assert evaluation.objective == 2 + assert evaluation.objective_gradient == (1, 0) diff --git a/simopt/dsl.py b/simopt/dsl.py new file mode 100644 index 00000000..14a1a78e --- /dev/null +++ b/simopt/dsl.py @@ -0,0 +1,5 @@ +"""SimOpt's public interface to the modeling language.""" + +from simopt_dsl import Model, Simulation, Variable, VectorVariable, mean, sum # noqa: A004 + +__all__ = ["Model", "Simulation", "Variable", "VectorVariable", "mean", "sum"] diff --git a/simopt/models/ambulance.py b/simopt/models/ambulance.py index b4c0d09b..aa91b319 100644 --- a/simopt/models/ambulance.py +++ b/simopt/models/ambulance.py @@ -10,12 +10,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import Beta, Exp @@ -288,45 +287,20 @@ class AmbulanceMinAvgResponse(Problem): } model_decision_factors: ClassVar[set[str]] = {"variable_locs"} - @property @override - def dim(self) -> int: - return int(2 * self.model.factors["variable_base_count"]) - - @property - @override - def lower_bounds(self) -> tuple: - return tuple(0.0 for _ in range(self.dim)) - - @property - @override - def upper_bounds(self) -> tuple: - return tuple(20.0 for _ in range(self.dim)) + def build(self) -> dsl.Model: + problem = dsl.Model() + initial = self.factors["initial_solution"] + dim = len(initial) + location = problem.add_continuous_vector(lb=0.0, ub=20.0, shape=(dim,), initial=initial) + simulation = self.add_simulation(problem, {"variable_locs": location}) + problem.minimize(dsl.mean(simulation.metric("avg_response_time"))) + return problem @override def vector_to_factor_dict(self, vector: tuple) -> dict: return {"variable_locs": list(vector)} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - # 1. Run the simulation - responses, gradients = self.model.replicate(model_factors, rngs) - - # 2. Construct the Objective - # Since this problem has no deterministic cost component, - # deterministic values are 0. - objectives = [ - Objective( - stochastic=responses["avg_response_time"], - stochastic_gradients=gradients["avg_response_time"]["variable_locs"], - deterministic=0.0, - deterministic_gradients=(0.0,) * self.dim, - ) - ] - - # 3. Return result - return RepResult(objectives=objectives) - @override def check_deterministic_constraints(self, _x: tuple) -> bool: return len(_x) == self.dim and all(0 <= xi <= 20 for xi in _x) diff --git a/simopt/models/amusementpark.py b/simopt/models/amusementpark.py index f5d54c30..75e071de 100644 --- a/simopt/models/amusementpark.py +++ b/simopt/models/amusementpark.py @@ -9,12 +9,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import Exp, Gamma, WeightedChoice @@ -440,28 +439,27 @@ class AmusementParkMinDepart(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"queue_capacities"} - @property - def dim(self) -> int: - return self.model.factors["number_attractions"] - - @property - def lower_bounds(self) -> tuple: - return (0,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (self.model.factors["park_capacity"],) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + queue_capacities = problem.add_integer_vector( + lb=0, + ub=self.model.factors["park_capacity"], + shape=(self.model.factors["number_attractions"],), + initial=tuple(self.factors["initial_solution"]), + ) + problem.add_linear_constraint( + dsl.sum(queue_capacities) <= self.model.factors["park_capacity"] + ) + simulation = self.add_simulation(problem, {"queue_capacities": queue_capacities}) + problem.minimize(dsl.mean(simulation.metric("total_departed"))) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict[str, tuple]: return { "queue_capacities": vector[:], } - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - return RepResult(objectives=[Objective(stochastic=responses["total_departed"])]) - def check_deterministic_constraints(self, x: tuple) -> bool: # Check box constraints. if not super().check_deterministic_constraints(x): diff --git a/simopt/models/chessmm.py b/simopt/models/chessmm.py index edb81d67..9c356173 100644 --- a/simopt/models/chessmm.py +++ b/simopt/models/chessmm.py @@ -12,13 +12,11 @@ from scipy import special from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, - StochasticConstraint, VariableType, ) from simopt.input_models import Exp, InputModel @@ -248,34 +246,22 @@ class ChessAvgDifference(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"allowable_diff"} - @property - def dim(self) -> int: - return 1 - - @property - def lower_bounds(self) -> tuple: - return (0,) - - @property - def upper_bounds(self) -> tuple: - return (2400,) + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + allowable_diff = problem.add_continuous_variable( + lb=0.0, ub=2400.0, initial=self.factors["initial_solution"][0] + ) + simulation = self.add_simulation(problem, {"allowable_diff": allowable_diff}) + problem.minimize(dsl.mean(simulation.metric("avg_diff"))) + problem.add_stochastic_constraint( + simulation.metric("avg_wait_time") <= self.factors["upper_time"] + ) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"allowable_diff": vector[0]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - return RepResult( - objectives=[Objective(stochastic=responses["avg_diff"])], - stochastic_constraints=[ - StochasticConstraint( - stochastic=responses["avg_wait_time"], - deterministic=-1 * self.factors["upper_time"], - ) - ], - ) - def check_deterministic_constraints(self, x: tuple) -> bool: return all(x_val > 0 for x_val in x) diff --git a/simopt/models/cntnv.py b/simopt/models/cntnv.py index 77af925b..84609dec 100644 --- a/simopt/models/cntnv.py +++ b/simopt/models/cntnv.py @@ -9,12 +9,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import InputModel @@ -254,33 +253,19 @@ class CntNVMaxProfit(Problem): } model_decision_factors: ClassVar[set[str]] = {"order_quantity"} - @property - def dim(self) -> int: - return 1 - - @property - def lower_bounds(self) -> tuple: - return (0,) - - @property - def upper_bounds(self) -> tuple: - return (np.inf,) + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + order_quantity = problem.add_continuous_variable( + lb=0.0, ub=np.inf, initial=self.factors["initial_solution"][0] + ) + simulation = self.add_simulation(problem, {"order_quantity": order_quantity}) + problem.maximize(dsl.mean(simulation.metric("profit"))) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"order_quantity": vector[0]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, gradients = self.model.replicate(model_factors, rngs) - return RepResult( - objectives=[ - Objective( - stochastic=responses["profit"], - stochastic_gradients=gradients["profit"]["order_quantity"], - ) - ], - ) - def check_deterministic_constraints(self, x: tuple) -> bool: return x[0] > 0 diff --git a/simopt/models/contam.py b/simopt/models/contam.py index 238afd70..c3b3bbc0 100644 --- a/simopt/models/contam.py +++ b/simopt/models/contam.py @@ -8,13 +8,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, - StochasticConstraint, VariableType, ) from simopt.input_models import Beta @@ -365,40 +363,46 @@ class ContaminationTotalCostDisc(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"prev_decision"} - @property - def dim(self) -> int: - return self.model.factors["stages"] - - @property - def lower_bounds(self) -> tuple: - return (0,) * self.model.factors["stages"] - - @property - def upper_bounds(self) -> tuple: - return (1,) * self.model.factors["stages"] + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + stages = self.model.factors["stages"] + prev_decision = problem.add_integer_vector( + lb=0, ub=1, shape=(stages,), initial=tuple(self.factors["initial_solution"]) + ) + + def run( + decisions: dict[str, float | tuple[float, ...]], + rngs: list[MRG32k3a], + ) -> tuple[dict, dict]: + prevention = decisions["prev_decision"] + if not isinstance(prevention, tuple): + raise TypeError("prev_decision must be a vector") + decision_factors = {"prev_decision": prevention} + responses, _ = self.model.replicate(self.model.factors | decision_factors, rngs) + under_control = np.asarray(responses["level"]) <= np.asarray( + self.factors["upper_thres"] + ) + return {"under_control": under_control.astype(float)}, {} + + simulation = problem.add_simulation( + run=run, decisions={"prev_decision": prev_decision}, n_rngs=self.model.n_rngs + ) + problem.minimize( + sum( + cost * prevention + for cost, prevention in zip(self.factors["prev_cost"], prev_decision, strict=True) + ) + ) + for stage, error_probability in enumerate(self.factors["error_prob"]): + problem.add_stochastic_constraint( + (1.0 - error_probability) - simulation.metric("under_control")[stage] <= 0.0 + ) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"prev_decision": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - x = tuple(model_factors["prev_decision"]) - responses, _ = self.model.replicate(model_factors, rngs) - objectives = [ - Objective( - stochastic=0.0, - deterministic=np.dot(self.factors["prev_cost"], x), - deterministic_gradients=self.factors["prev_cost"], - ) - ] - under_control = responses["level"] <= self.factors["upper_thres"] - error_prob = self.factors["error_prob"] - stochastic_constraints = [ - StochasticConstraint(stochastic=-1 * under_control[i], deterministic=1 - error_prob[i]) - for i in range(len(under_control)) - ] - return RepResult(objectives=objectives, stochastic_constraints=stochastic_constraints) - def check_deterministic_constraints(self, x: tuple) -> bool: return all(0 <= u <= 1 for u in x) @@ -424,54 +428,46 @@ class ContaminationTotalCostCont(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"prev_decision"} - @property - def dim(self) -> int: - return self.model.factors["stages"] - - @property - def lower_bounds(self) -> tuple: - return (0,) * self.model.factors["stages"] - - @property - def upper_bounds(self) -> tuple: - return (1,) * self.model.factors["stages"] - - # # TODO: figure out how Problem.check_simulatable_factors() works - # def check_simulatable_factors(self) -> bool: - # lower_len = len(self.lower_bounds) - # upper_len = len(self.upper_bounds) - # if lower_len != upper_len or lower_len != self.dim: - # error_msg = ( - # f"Lower bounds: {lower_len}, " - # f"Upper bounds: {upper_len}, " - # f"Dim: {self.dim}" - # ) - # raise ValueError(error_msg) - # return True + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + stages = self.model.factors["stages"] + prev_decision = problem.add_continuous_vector( + lb=0.0, ub=1.0, shape=(stages,), initial=tuple(self.factors["initial_solution"]) + ) + + def run( + decisions: dict[str, float | tuple[float, ...]], + rngs: list[MRG32k3a], + ) -> tuple[dict, dict]: + prevention = decisions["prev_decision"] + if not isinstance(prevention, tuple): + raise TypeError("prev_decision must be a vector") + decision_factors = {"prev_decision": prevention} + responses, _ = self.model.replicate(self.model.factors | decision_factors, rngs) + under_control = np.asarray(responses["level"]) <= np.asarray( + self.factors["upper_thres"] + ) + return {"under_control": under_control.astype(float)}, {} + + simulation = problem.add_simulation( + run=run, decisions={"prev_decision": prev_decision}, n_rngs=self.model.n_rngs + ) + problem.minimize( + sum( + cost * prevention + for cost, prevention in zip(self.factors["prev_cost"], prev_decision, strict=True) + ) + ) + for stage, error_probability in enumerate(self.factors["error_prob"]): + problem.add_stochastic_constraint( + (1.0 - error_probability) - simulation.metric("under_control")[stage] <= 0.0 + ) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"prev_decision": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - x = tuple(model_factors["prev_decision"]) - responses, _ = self.model.replicate(model_factors, rngs) - deterministic_cost = np.dot(self.factors["prev_cost"], x) - objectives = [ - Objective( - stochastic=0.0, - deterministic=deterministic_cost, - deterministic_gradients=self.factors["prev_cost"], - ) - ] - under_control = responses["level"] <= self.factors["upper_thres"] - error_prob = self.factors["error_prob"] - stochastic_constraints = [ - StochasticConstraint(stochastic=-1 * under_control[i], deterministic=1 - error_prob[i]) - for i in range(len(under_control)) - ] - return RepResult(objectives=objectives, stochastic_constraints=stochastic_constraints) - def check_deterministic_constraints(self, x: tuple) -> bool: return all(0 <= u <= 1 for u in x) diff --git a/simopt/models/dualsourcing.py b/simopt/models/dualsourcing.py index d876d8bd..2808e22a 100644 --- a/simopt/models/dualsourcing.py +++ b/simopt/models/dualsourcing.py @@ -9,12 +9,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import InputModel @@ -305,17 +304,24 @@ class DualSourcingMinCost(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"order_level_exp", "order_level_reg"} - @property - def dim(self) -> int: - return 2 - - @property - def lower_bounds(self) -> tuple: - return (0, 0) - - @property - def upper_bounds(self) -> tuple: - return (np.inf, np.inf) + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + initial_solution = self.factors["initial_solution"] + order_level_exp = problem.add_integer_variable(lb=0, ub=np.inf, initial=initial_solution[0]) + order_level_reg = problem.add_integer_variable(lb=0, ub=np.inf, initial=initial_solution[1]) + + simulation = self.add_simulation( + problem, {"order_level_exp": order_level_exp, "order_level_reg": order_level_reg} + ) + problem.minimize( + dsl.mean( + simulation.metric("average_ordering_cost") + + simulation.metric("average_penalty_cost") + + simulation.metric("average_holding_cost") + ) + ) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return { @@ -323,19 +329,6 @@ def vector_to_factor_dict(self, vector: tuple) -> dict: "order_level_reg": vector[1], } - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - return RepResult( - objectives=[ - Objective( - stochastic=responses["average_ordering_cost"] - + responses["average_penalty_cost"] - + responses["average_holding_cost"], - ) - ], - ) - def check_deterministic_constraints(self, x: tuple) -> bool: return x[0] >= 0 and x[1] >= 0 diff --git a/simopt/models/dynamnews.py b/simopt/models/dynamnews.py index e1190d19..1a306076 100644 --- a/simopt/models/dynamnews.py +++ b/simopt/models/dynamnews.py @@ -10,12 +10,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import InputModel @@ -300,27 +299,22 @@ class DynamNewsMaxProfit(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"init_level"} - @property - def dim(self) -> int: - return self.model.factors["num_prod"] - - @property - def lower_bounds(self) -> tuple: - return (0,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (np.inf,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + init_level = problem.add_continuous_vector( + lb=0.0, + ub=np.inf, + shape=(self.model.factors["num_prod"],), + initial=self.factors["initial_solution"], + ) + simulation = self.add_simulation(problem, {"init_level": init_level}) + problem.maximize(dsl.mean(simulation.metric("profit"))) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"init_level": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - objectives = [Objective(stochastic=responses["profit"])] - return RepResult(objectives=objectives) - def check_deterministic_constraints(self, x: tuple) -> bool: return all(x[j] > 0 for j in range(self.dim)) diff --git a/simopt/models/ermexample.py b/simopt/models/ermexample.py index 36fae180..540ac037 100644 --- a/simopt/models/ermexample.py +++ b/simopt/models/ermexample.py @@ -12,12 +12,11 @@ from pydantic import BaseModel, Field from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import InputModel @@ -95,7 +94,6 @@ def __init__(self, fixed_factors: dict | None = None) -> None: super().__init__(fixed_factors) self.resample_model = FileInputModel("workshop/erm_data.npy") - def replicate(self, factors: dict, rngs: list[MRG32k3a]) -> tuple[dict, dict]: """Evaluate the squared error loss of a single observation. @@ -142,9 +140,7 @@ def optimal_value(self) -> float | None: # noqa: D102 x = all_data[:, 0] y = all_data[:, 1] optbeta1, optbeta0 = np.polyfit(x, y, 1) - opttrainingmse = np.mean( - [(yy - optbeta0 - optbeta1 * xx) ** 2 for (xx, yy) in zip(x, y)] - ) + opttrainingmse = np.mean([(yy - optbeta0 - optbeta1 * xx) ** 2 for (xx, yy) in zip(x, y)]) return opttrainingmse @property @@ -156,32 +152,17 @@ def optimal_solution(self) -> tuple | None: # noqa: D102 optbeta1, optbeta0 = np.polyfit(x, y, 1) return (optbeta0, optbeta1) - @property - def dim(self) -> int: # noqa: D102 - return 2 - - @property - def lower_bounds(self) -> tuple: # noqa: D102 - return (-np.inf,) * self.dim - - @property - def upper_bounds(self) -> tuple: # noqa: D102 - return (np.inf,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + beta = problem.add_continuous_vector(lb=-np.inf, ub=np.inf, shape=(2,), initial=self.factors["initial_solution"]) + simulation = self.add_simulation(problem, {"beta": beta}) + problem.minimize(dsl.mean(simulation.metric("sq_error_loss"))) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: # noqa: D102 return {"beta": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: # noqa: D102 - responses, gradients = self.model.replicate(model_factors, rngs) - objectives = [ - Objective( - stochastic=responses["sq_error_loss"], - stochastic_gradients=gradients["sq_error_loss"]["beta"], - ) - ] - return RepResult(objectives=objectives) - def get_random_solution(self, rand_sol_rng: MRG32k3a) -> tuple: # noqa: D102 # beta = tuple([rand_sol_rng.uniform(-2, 2) for _ in range(self.dim)]) beta = tuple( diff --git a/simopt/models/example.py b/simopt/models/example.py index fdccb717..68b3cc8d 100644 --- a/simopt/models/example.py +++ b/simopt/models/example.py @@ -12,12 +12,11 @@ from pydantic import BaseModel, Field from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import Normal @@ -128,32 +127,20 @@ def optimal_solution(self) -> tuple | None: # TODO: figure out what f is return (0,) * self.dim - @property - def dim(self) -> int: - return len(self.factors["initial_solution"]) - - @property - def lower_bounds(self) -> tuple: - return (-np.inf,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (np.inf,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + initial_solution = self.factors["initial_solution"] + x = problem.add_continuous_vector( + lb=-np.inf, ub=np.inf, shape=(len(initial_solution),), initial=initial_solution + ) + simulation = self.add_simulation(problem, {"x": x}) + problem.minimize(dsl.mean(simulation.metric("est_f(x)"))) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"x": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, gradients = self.model.replicate(model_factors, rngs) - objectives = [ - Objective( - stochastic=responses["est_f(x)"], - stochastic_gradients=gradients["est_f(x)"]["x"], - ) - ] - return RepResult(objectives=objectives) - def get_random_solution(self, rand_sol_rng: MRG32k3a) -> tuple: # x = tuple([rand_sol_rng.uniform(-2, 2) for _ in range(self.dim)]) return tuple( @@ -258,26 +245,19 @@ def optimal_value(self) -> float | None: def optimal_solution(self) -> tuple | None: return (1, 2, 3, 4) - @property - def dim(self) -> int: - return 4 - - @property - def lower_bounds(self) -> tuple: - return (-4,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (4,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + initial_solution = self.factors["initial_solution"] + x = problem.add_integer_vector( + lb=-4, ub=4, shape=(len(initial_solution),), initial=initial_solution + ) + simulation = self.add_simulation(problem, {"x": x}) + problem.minimize(dsl.mean(simulation.metric("est_f(x)"))) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"x": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - objectives = [Objective(stochastic=responses["est_f(x)"])] - return RepResult(objectives=objectives) - def get_random_solution(self, rand_sol_rng: MRG32k3a) -> tuple: return tuple(rand_sol_rng.randint(-4, 4) for _ in range(self.dim)) diff --git a/simopt/models/facilitysizing.py b/simopt/models/facilitysizing.py index 90e77bac..cf166773 100644 --- a/simopt/models/facilitysizing.py +++ b/simopt/models/facilitysizing.py @@ -9,13 +9,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, - StochasticConstraint, VariableType, ) from simopt.input_models import InputModel @@ -299,40 +297,28 @@ class FacilitySizingTotalCost(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"capacity"} - @property - def dim(self) -> int: - return self.model.factors["n_fac"] - - @property - def lower_bounds(self) -> tuple: - return (0,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (np.inf,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + capacity = problem.add_continuous_vector( + lb=0.0, + ub=np.inf, + shape=(self.model.factors["n_fac"],), + initial=self.factors["initial_solution"], + ) + simulation = self.add_simulation(problem, {"capacity": capacity}) + installation_cost = sum( + cost * capacity[index] for index, cost in enumerate(self.factors["installation_costs"]) + ) + problem.minimize(installation_cost) + problem.add_stochastic_constraint( + simulation.metric("stockout_flag") <= self.factors["epsilon"] + ) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"capacity": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - x = tuple(model_factors["capacity"]) - objectives = [ - Objective( - stochastic=0.0, - deterministic=np.dot(self.factors["installation_costs"], x), - deterministic_gradients=self.factors["installation_costs"], - ) - ] - stochastic_constraints = [ - StochasticConstraint( - stochastic=responses["stockout_flag"], - deterministic=-self.factors["epsilon"], - ) - ] - return RepResult(objectives=objectives, stochastic_constraints=stochastic_constraints) - def get_random_solution(self, rand_sol_rng: MRG32k3a) -> tuple: cov_matrix = np.diag([x**2 for x in self.factors["initial_solution"]]) x = rand_sol_rng.mvnormalvariate( @@ -363,28 +349,26 @@ class FacilitySizingMaxService(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"capacity"} - @property - def dim(self) -> int: - return self.model.factors["n_fac"] - - @property - def lower_bounds(self) -> tuple: - return (0,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (np.inf,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + capacity = problem.add_continuous_vector( + lb=0.0, + ub=np.inf, + shape=(self.model.factors["n_fac"],), + initial=self.factors["initial_solution"], + ) + simulation = self.add_simulation(problem, {"capacity": capacity}) + problem.maximize(1 - dsl.mean(simulation.metric("stockout_flag"))) + installation_cost = sum( + cost * capacity[index] for index, cost in enumerate(self.factors["installation_costs"]) + ) + problem.add_linear_constraint(installation_cost <= self.factors["installation_budget"]) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"capacity": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - service_value = 1 - responses["stockout_flag"] - objectives = [Objective(stochastic=service_value)] - return RepResult(objectives=objectives) - def check_deterministic_constraints(self, x: tuple) -> bool: # Check budget constraint budget_feasible = ( diff --git a/simopt/models/fixedsan.py b/simopt/models/fixedsan.py index b133cbcb..44bc4c80 100644 --- a/simopt/models/fixedsan.py +++ b/simopt/models/fixedsan.py @@ -8,12 +8,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import Exp @@ -243,17 +242,22 @@ class FixedSANLongestPath(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"arc_means"} - @property - def dim(self) -> int: - return self.model.factors["num_arcs"] - - @property - def lower_bounds(self) -> tuple: - return (1e-2,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (np.inf,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + arc_means = problem.add_continuous_vector( + lb=1e-2, + ub=np.inf, + shape=(self.model.factors["num_arcs"],), + initial=tuple(self.factors["initial_solution"]), + ) + simulation = self.add_simulation(problem, {"arc_means": arc_means}) + deterministic_cost = sum( + cost / arc_mean + for cost, arc_mean in zip(self.factors["arc_costs"], arc_means, strict=True) + ) + problem.minimize(dsl.mean(simulation.metric("longest_path_length")) + deterministic_cost) + return problem def check_arc_costs(self) -> bool: """Check if all arc costs are positive and match the number of arcs.""" @@ -264,20 +268,6 @@ def check_arc_costs(self) -> bool: def vector_to_factor_dict(self, vector: tuple) -> dict: return {"arc_means": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, gradients = self.model.replicate(model_factors, rngs) - x = tuple(model_factors["arc_means"]) - objectives = [ - Objective( - stochastic=responses["longest_path_length"], - stochastic_gradients=gradients["longest_path_length"]["arc_means"], - deterministic=np.sum(np.array(self.factors["arc_costs"]) / np.array(x)), - deterministic_gradients=-np.array(self.factors["arc_costs"]) / (np.array(x) ** 2), - ) - ] - return RepResult(objectives=objectives) - def check_deterministic_constraints(self, x: tuple) -> bool: return all(x_i >= 0 for x_i in x) diff --git a/simopt/models/hotel.py b/simopt/models/hotel.py index e176ec13..36ea6d0e 100644 --- a/simopt/models/hotel.py +++ b/simopt/models/hotel.py @@ -10,12 +10,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import Exp @@ -354,36 +353,25 @@ class HotelRevenue(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"booking_limits"} - @property - def dim(self) -> int: - return self.model.factors["num_products"] - - @property - def lower_bounds(self) -> tuple: - return (0,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (self.model.factors["num_rooms"],) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + booking_limits = problem.add_integer_vector( + lb=0, + ub=self.model.factors["num_rooms"], + shape=(self.model.factors["num_products"],), + initial=self.factors["initial_solution"], + ) - # # TODO: figure out how Problem.check_simulatable_factors() works - # def check_simulatable_factors(self) -> bool: - # return not ( - # len(self.lower_bounds) != self.dim or len(self.upper_bounds) != self.dim - # ) + simulation = self.add_simulation(problem, {"booking_limits": booking_limits}) + problem.maximize(dsl.mean(simulation.metric("revenue"))) + return problem + @override def vector_to_factor_dict(self, vector: tuple) -> dict: return {"booking_limits": vector[:]} @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - objectives = [Objective(stochastic=responses["revenue"])] - return RepResult(objectives=objectives) - - def check_deterministic_constraints(self, _x: tuple) -> bool: - return True - def get_random_solution(self, rand_sol_rng: MRG32k3a) -> tuple: return tuple( [rand_sol_rng.randint(0, self.model.factors["num_rooms"]) for _ in range(self.dim)] diff --git a/simopt/models/ironore.py b/simopt/models/ironore.py index 1721314f..67969caa 100644 --- a/simopt/models/ironore.py +++ b/simopt/models/ironore.py @@ -14,12 +14,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import InputModel @@ -362,17 +361,25 @@ class IronOreMaxRev(Problem): "price_sell", } - @property - def dim(self) -> int: - return 4 - - @property - def lower_bounds(self) -> tuple: - return (0,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (np.inf,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + initial_solution = self.factors["initial_solution"] + price_prod = problem.add_continuous_variable(lb=0.0, ub=np.inf, initial=initial_solution[0]) + inven_stop = problem.add_integer_variable(lb=0, ub=np.inf, initial=initial_solution[1]) + price_stop = problem.add_continuous_variable(lb=0.0, ub=np.inf, initial=initial_solution[2]) + price_sell = problem.add_continuous_variable(lb=0.0, ub=np.inf, initial=initial_solution[3]) + simulation = self.add_simulation( + problem, + { + "price_prod": price_prod, + "inven_stop": inven_stop, + "price_stop": price_stop, + "price_sell": price_sell, + }, + ) + problem.maximize(dsl.mean(simulation.metric("total_profit"))) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return { @@ -382,12 +389,6 @@ def vector_to_factor_dict(self, vector: tuple) -> dict: "price_sell": vector[3], } - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - objectives = [Objective(stochastic=responses["total_profit"])] - return RepResult(objectives=objectives) - def get_random_solution(self, rand_sol_rng: MRG32k3a) -> tuple: # return ( # rand_sol_rng.randint(70, 90), @@ -425,17 +426,23 @@ class IronOreMaxRevCnt(Problem): "price_sell", } - @property - def dim(self) -> int: - return 3 - - @property - def lower_bounds(self) -> tuple: - return (0.0,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (np.inf,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + initial_solution = self.factors["initial_solution"] + price_prod = problem.add_continuous_variable(lb=0.0, ub=np.inf, initial=initial_solution[0]) + price_stop = problem.add_continuous_variable(lb=0.0, ub=np.inf, initial=initial_solution[1]) + price_sell = problem.add_continuous_variable(lb=0.0, ub=np.inf, initial=initial_solution[2]) + simulation = self.add_simulation( + problem, + { + "price_prod": price_prod, + "price_stop": price_stop, + "price_sell": price_sell, + }, + ) + problem.maximize(dsl.mean(simulation.metric("total_profit"))) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return { @@ -444,12 +451,6 @@ def vector_to_factor_dict(self, vector: tuple) -> dict: "price_sell": vector[2], } - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - objectives = [Objective(stochastic=responses["total_profit"])] - return RepResult(objectives=objectives) - def check_deterministic_constraints(self, x: tuple) -> bool: return x[0] >= 0 and x[1] >= 0 and x[2] >= 0 diff --git a/simopt/models/mm1queue.py b/simopt/models/mm1queue.py index 0361b117..94c1168f 100644 --- a/simopt/models/mm1queue.py +++ b/simopt/models/mm1queue.py @@ -11,12 +11,11 @@ from pydantic import BaseModel, Field from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import Exp @@ -289,35 +288,21 @@ class MM1MinMeanSojournTime(Problem): model_default_factors: ClassVar[dict] = {"warmup": 50, "people": 200} model_decision_factors: ClassVar[set[str]] = {"mu"} - @property - def dim(self) -> int: - return 1 - - @property - def lower_bounds(self) -> tuple: - return (0,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (np.inf,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + mu = problem.add_continuous_variable( + lb=0.0, ub=np.inf, initial=self.factors["initial_solution"][0] + ) + simulation = self.add_simulation(problem, {"mu": mu}) + problem.minimize( + dsl.mean(simulation.metric("avg_sojourn_time")) + self.factors["cost"] * mu * mu + ) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"mu": vector[0]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, gradients = self.model.replicate(model_factors, rngs) - x = (model_factors["mu"],) - objectives = [ - Objective( - stochastic=responses["avg_sojourn_time"], - stochastic_gradients=gradients["avg_sojourn_time"]["mu"], - deterministic=self.factors["cost"] * (x[0] ** 2), - deterministic_gradients=2 * self.factors["cost"] * x[0], - ) - ] - return RepResult(objectives=objectives) - def get_random_solution(self, rand_sol_rng: MRG32k3a) -> tuple: # Generate an Exponential(rate = 1/3) r.v. return (rand_sol_rng.expovariate(1 / 3),) diff --git a/simopt/models/network.py b/simopt/models/network.py index 3ccb4e3b..cf2c8b10 100644 --- a/simopt/models/network.py +++ b/simopt/models/network.py @@ -12,12 +12,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import Exp, InputModel, Triangular @@ -345,27 +344,23 @@ class NetworkMinTotalCost(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"process_prob"} - @property - def dim(self) -> int: - return self.model.factors["n_networks"] - - @property - def lower_bounds(self) -> tuple: - return (0,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (1,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + n_networks = self.model.factors["n_networks"] + process_prob = problem.add_continuous_vector( + lb=0.0, ub=1.0, shape=(n_networks,), initial=self.factors["initial_solution"] + ) + total_probability = dsl.sum(process_prob) + problem.add_linear_constraint(total_probability <= 1.0 + 1e-10) + problem.add_linear_constraint(total_probability >= 1.0 - 1e-10) + simulation = self.add_simulation(problem, {"process_prob": process_prob}) + problem.minimize(dsl.mean(simulation.metric("total_cost"))) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"process_prob": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - objectives = [Objective(stochastic=responses["total_cost"])] - return RepResult(objectives=objectives) - def check_deterministic_constraints(self, x: tuple) -> bool: # Check box constraints. box_feasible = super().check_deterministic_constraints(x) diff --git a/simopt/models/paramesti.py b/simopt/models/paramesti.py index 64a3c542..bd5b8b50 100644 --- a/simopt/models/paramesti.py +++ b/simopt/models/paramesti.py @@ -9,12 +9,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import Gamma @@ -164,27 +163,19 @@ def optimal_solution(self) -> tuple | None: return tuple(solution) return solution - @property - def dim(self) -> int: - return 2 - - @property - def lower_bounds(self) -> tuple: - return (0.1,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (10,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + x = problem.add_continuous_vector( + lb=0.1, ub=10.0, shape=(2,), initial=self.factors["initial_solution"] + ) + simulation = self.add_simulation(problem, {"x": x}) + problem.maximize(dsl.mean(simulation.metric("loglik"))) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"x": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - objectives = [Objective(stochastic=responses["loglik"])] - return RepResult(objectives=objectives) - def check_deterministic_constraints(self, _x: tuple) -> bool: return True diff --git a/simopt/models/rmitd.py b/simopt/models/rmitd.py index efd59379..feb28336 100644 --- a/simopt/models/rmitd.py +++ b/simopt/models/rmitd.py @@ -8,12 +8,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import InputModel @@ -266,30 +265,39 @@ class RMITDMaxRevenue(Problem): "reservation_qtys", } - @property - def dim(self) -> int: - return 3 - - @property - def lower_bounds(self) -> tuple: - return (0,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (np.inf,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + initial_solution = tuple(self.factors["initial_solution"]) + initial_inventory = problem.add_integer_variable( + lb=0, ub=np.inf, initial=initial_solution[0] + ) + reservation_qtys = problem.add_integer_vector( + lb=0, + ub=np.inf, + shape=(self.model.factors["time_horizon"] - 1,), + initial=initial_solution[1:], + ) + problem.add_linear_constraint(initial_inventory >= reservation_qtys[0]) + for idx in range(len(reservation_qtys) - 1): + problem.add_linear_constraint(reservation_qtys[idx] >= reservation_qtys[idx + 1]) + + simulation = self.add_simulation( + problem, + { + "initial_inventory": initial_inventory, + "reservation_qtys": reservation_qtys, + }, + ) + problem.maximize(dsl.mean(simulation.metric("revenue"))) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return { "initial_inventory": vector[0], - "reservation_qtys": list(vector[0:]), + "reservation_qtys": list(vector[1:]), } - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - objectives = [Objective(stochastic=responses["revenue"])] - return RepResult(objectives=objectives) - def check_deterministic_constraints(self, x: tuple) -> bool: return all(x[idx] >= x[idx + 1] for idx in range(self.dim - 1)) diff --git a/simopt/models/san.py b/simopt/models/san.py index 8040e087..35d9a63f 100644 --- a/simopt/models/san.py +++ b/simopt/models/san.py @@ -9,13 +9,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, - StochasticConstraint, VariableType, ) from simopt.input_models import Exp @@ -247,7 +245,7 @@ def replicate(self, factors: dict, rngs: list[MRG32k3a]) -> tuple[dict, dict]: # If an arc is not on the longest path, the component of the gradient is zero. arc_to_index = {arc: i for i, arc in enumerate(arcs)} - grads = [] + grads = np.zeros((num_nodes, len(arcs))) for node in topo_order: gradient = np.zeros(len(arcs)) current = node @@ -260,7 +258,7 @@ def replicate(self, factors: dict, rngs: list[MRG32k3a]) -> tuple[dict, dict]: current = backtrack backtrack = int(prev[backtrack - 1]) - grads.append(gradient) + grads[node - 1] = gradient # Compose responses and gradients. responses = { @@ -269,13 +267,9 @@ def replicate(self, factors: dict, rngs: list[MRG32k3a]) -> tuple[dict, dict]: "topo_order": topo_order, } gradients = { - response_key: { - factor_key: np.zeros(len(self.specifications)) for factor_key in self.specifications - } - for response_key in responses + "longest_path_length": {"arc_means": grads[-1]}, + "longest_path_to_all_nodes": {"arc_means": grads}, } - gradients["longest_path_length"]["arc_means"] = grads[-1] - gradients["longest_path_to_all_nodes"]["arc_means"] = np.array(grads) return responses, gradients @@ -297,35 +291,26 @@ class SANLongestPath(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"arc_means"} - @property - def dim(self) -> int: - return len(self.model.factors["arcs"]) - - @property - def lower_bounds(self) -> tuple: - return (1e-2,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (np.inf,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + arc_means = problem.add_continuous_vector( + lb=1e-2, + ub=np.inf, + shape=(len(self.model.factors["arcs"]),), + initial=tuple(self.factors["initial_solution"]), + ) + simulation = self.add_simulation(problem, {"arc_means": arc_means}) + deterministic_cost = sum( + cost / arc_mean + for cost, arc_mean in zip(self.factors["arc_costs"], arc_means, strict=True) + ) + problem.minimize(dsl.mean(simulation.metric("longest_path_length")) + deterministic_cost) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"arc_means": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, gradients = self.model.replicate(model_factors, rngs) - x = tuple(model_factors["arc_means"]) - objectives = [ - Objective( - stochastic=responses["longest_path_length"], - stochastic_gradients=gradients["longest_path_length"]["arc_means"], - deterministic=np.sum(np.array(self.factors["arc_costs"]) / np.array(x)), - deterministic_gradients=-np.array(self.factors["arc_costs"]) / (np.array(x) ** 2), - ) - ] - return RepResult(objectives=objectives) - def check_deterministic_constraints(self, x: tuple) -> bool: return all(x_i >= 0 for x_i in x) @@ -406,62 +391,32 @@ class SANLongestPathStochastic(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"arc_means"} - @property - def dim(self) -> int: - return len(self.model.factors["arcs"]) - - @property - def lower_bounds(self) -> tuple: - return (1e-2,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (100.0,) * self.dim - - def vector_to_factor_dict(self, vector: tuple) -> dict: - return {"arc_means": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, gradients = self.model.replicate(model_factors, rngs) - x = tuple(model_factors["arc_means"]) - - objectives = [ - Objective( - stochastic=responses["longest_path_length"], - stochastic_gradients=gradients["longest_path_length"]["arc_means"], - deterministic=np.sum(np.array(self.factors["arc_costs"]) / np.array(x)), - deterministic_gradients=-np.array(self.factors["arc_costs"]) / (np.array(x) ** 2), - ) - ] - - topo_order = responses["topo_order"] - longest_path_nodes = responses["longest_path_to_all_nodes"] - node_positions = {node: idx for idx, node in enumerate(topo_order)} - arc_gradients = gradients["longest_path_to_all_nodes"]["arc_means"] - constraint_limits = self.factors["length_to_node_constraint"] - constraint_nodes = self.factors["constraint_nodes"] - - stochastic_constraints = [] - for i, const_node in enumerate(constraint_nodes): - idx = node_positions[const_node] - stochastic_value = longest_path_nodes[idx] - stochastic_grad = arc_gradients[idx] - deterministic_value = -constraint_limits[i] - deterministic_grad = [0.0] * self.dim - stochastic_constraints.append( - StochasticConstraint( - stochastic_value, - stochastic_grad, - deterministic_value, - deterministic_grad, - ) + def build(self) -> dsl.Model: + problem = dsl.Model() + arc_means = problem.add_continuous_vector( + lb=1e-2, + ub=100.0, + shape=(len(self.model.factors["arcs"]),), + initial=tuple(self.factors["initial_solution"]), + ) + constraint_nodes = tuple(self.factors["constraint_nodes"]) + simulation = self.add_simulation(problem, {"arc_means": arc_means}) + deterministic_cost = sum( + cost / arc_mean + for cost, arc_mean in zip(self.factors["arc_costs"], arc_means, strict=True) + ) + problem.minimize(dsl.mean(simulation.metric("longest_path_length")) + deterministic_cost) + for node, limit in zip( + constraint_nodes, self.factors["length_to_node_constraint"], strict=True + ): + problem.add_stochastic_constraint( + simulation.metric("longest_path_to_all_nodes")[node - 1] <= limit ) + return problem - return RepResult( - objectives=objectives, - stochastic_constraints=stochastic_constraints, - ) + def vector_to_factor_dict(self, vector: tuple) -> dict: + return {"arc_means": vector[:]} def check_deterministic_constraints(self, x: tuple) -> bool: return all(x_i >= 0 for x_i in x) diff --git a/simopt/models/sscont.py b/simopt/models/sscont.py index 516a1bda..94bc6da8 100644 --- a/simopt/models/sscont.py +++ b/simopt/models/sscont.py @@ -9,6 +9,7 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, @@ -316,17 +317,25 @@ class SSContMinCost(Problem): model_default_factors: ClassVar[dict] = {"demand_mean": 100.0, "lead_mean": 6.0} model_decision_factors: ClassVar[set[str]] = {"s", "S"} - @property - def dim(self) -> int: - return 2 - - @property - def lower_bounds(self) -> tuple: - return (0,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + reorder_point = problem.add_continuous_variable( + lb=0.0, ub=np.inf, initial=self.factors["initial_solution"][0] + ) + order_gap = problem.add_continuous_variable( + lb=0.0, ub=np.inf, initial=self.factors["initial_solution"][1] + ) - @property - def upper_bounds(self) -> tuple: - return (np.inf,) * self.dim + simulation = self.add_simulation(problem, {"s": reorder_point, "order_gap": order_gap}) + problem.minimize( + dsl.mean( + simulation.metric("avg_backorder_costs") + + simulation.metric("avg_order_costs") + + simulation.metric("avg_holding_costs") + ) + ) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"s": vector[0], "S": vector[0] + vector[1]} diff --git a/simopt/models/tableallocation.py b/simopt/models/tableallocation.py index 8703cee6..23022c05 100644 --- a/simopt/models/tableallocation.py +++ b/simopt/models/tableallocation.py @@ -10,12 +10,11 @@ from pydantic import BaseModel, Field, model_validator from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.base import ( ConstraintType, Model, - Objective, Problem, - RepResult, VariableType, ) from simopt.input_models import Exp, Poisson, Uniform, WeightedChoice @@ -288,27 +287,29 @@ class TableAllocationMaxRev(Problem): model_default_factors: ClassVar[dict] = {} model_decision_factors: ClassVar[set[str]] = {"num_tables"} - @property - def dim(self) -> int: - return 4 - - @property - def lower_bounds(self) -> tuple: - return (0,) * self.dim - - @property - def upper_bounds(self) -> tuple: - return (np.inf,) * self.dim + @override + def build(self) -> dsl.Model: + problem = dsl.Model() + num_tables = problem.add_integer_vector( + lb=0, + ub=np.inf, + shape=(len(self.model.factors["table_cap"]),), + initial=tuple(self.factors["initial_solution"]), + ) + allocated_capacity = sum( + table_capacity * table_count + for table_capacity, table_count in zip( + self.model.factors["table_cap"], num_tables, strict=True + ) + ) + problem.add_linear_constraint(allocated_capacity <= self.model.factors["capacity"]) + simulation = self.add_simulation(problem, {"num_tables": num_tables}) + problem.maximize(dsl.mean(simulation.metric("total_revenue"))) + return problem def vector_to_factor_dict(self, vector: tuple) -> dict: return {"num_tables": vector[:]} - @override - def replicate(self, model_factors: dict, rngs: list[MRG32k3a]) -> RepResult: - responses, _ = self.model.replicate(model_factors, rngs) - objectives = [Objective(stochastic=responses["total_revenue"])] - return RepResult(objectives=objectives) - def check_deterministic_constraints(self, x: tuple) -> bool: return ( np.sum(np.multiply(self.model.factors["table_cap"], x)) diff --git a/simopt/problem.py b/simopt/problem.py index d0e77c41..a83af911 100644 --- a/simopt/problem.py +++ b/simopt/problem.py @@ -1,7 +1,7 @@ """Base classes for simulation optimization problems and models.""" from abc import ABC, abstractmethod -from collections.abc import Iterable +from collections.abc import Iterable, Mapping from typing import ClassVar import numpy as np @@ -9,6 +9,7 @@ from pydantic import BaseModel from mrg32k3a.mrg32k3a import MRG32k3a +from simopt import dsl from simopt.model import Model from simopt.problem_types import ConstraintType, VariableType from simopt.utils import get_specifications @@ -213,6 +214,7 @@ def __init__( # Set the model self.model = self.model_class(model_factors) + self._optimization_problem: dsl.Model | None = None self.rng_list: list[MRG32k3a] = [] def __eq__(self, other: object) -> bool: @@ -261,22 +263,38 @@ def optimal_solution(self) -> tuple | None: return None @property - @abstractmethod def dim(self) -> int: """Number of decision variables.""" - raise NotImplementedError + return len(self.optimization_problem.variables) @property - @abstractmethod - def lower_bounds(self) -> tuple[float, ...]: + def lower_bounds(self) -> tuple[int | float, ...]: """Lower bound for each decision variable.""" - raise NotImplementedError + problem = self.optimization_problem + return tuple( + int(bound) if variable.integer and np.isfinite(bound) else bound + for variable, bound in zip(problem.variables, problem.lower_bounds(), strict=True) + ) @property - @abstractmethod - def upper_bounds(self) -> tuple[float, ...]: + def upper_bounds(self) -> tuple[int | float, ...]: """Upper bound for each decision variable.""" - raise NotImplementedError + problem = self.optimization_problem + return tuple( + int(bound) if variable.integer and np.isfinite(bound) else bound + for variable, bound in zip(problem.variables, problem.upper_bounds(), strict=True) + ) + + @property + def optimization_problem(self) -> dsl.Model: + """Optimization problem expressed with the SimOpt DSL.""" + if self._optimization_problem is None: + self._optimization_problem = self.build() + return self._optimization_problem + + def build(self) -> dsl.Model: + """Build an optimization problem with the SimOpt DSL.""" + raise NotImplementedError(f"{type(self).__name__} does not build a SimOpt DSL problem") @classproperty def compatibility(cls) -> str: # noqa: N805 @@ -320,6 +338,44 @@ def vector_to_factor_dict(self, vector: tuple) -> dict: """ raise NotImplementedError + def add_simulation( + self, + target: dsl.Model, + decisions: Mapping[str, dsl.Variable | dsl.VectorVariable], + *, + name: str = "", + ) -> dsl.Simulation: + """Add this problem's simulation model to the DSL model. + + The problem supplies the replication adapter and RNG count so problem + implementations only need to declare the simulation's decisions. + + Args: + target: DSL model to which the simulation is added. + decisions: Decision variables passed to the simulation. + name: Name of the simulation. + + Returns: + The simulation added to ``target``. + """ + + def run( + decisions: dict[str, float | tuple[float, ...]], rngs: list[MRG32k3a] + ) -> tuple[dict, dict]: + decision_vector = [] + for decision in decisions.values(): + if isinstance(decision, tuple): + decision_vector.extend(decision) + else: + decision_vector.append(decision) + decision_vector = tuple(decision_vector) + decision_factors = self.vector_to_factor_dict(decision_vector) + return self.model.replicate(self.model.factors | decision_factors, rngs) + + return target.add_simulation( + name=name, run=run, decisions=decisions, n_rngs=self.model.n_rngs + ) + def check_deterministic_constraints(self, x: tuple, /) -> bool: """Check if a solution `x` satisfies the problem's deterministic constraints. @@ -350,15 +406,70 @@ def get_random_solution(self, rand_sol_rng: MRG32k3a) -> tuple: """ raise NotImplementedError - @abstractmethod def replicate(self, model_factors: dict, rngs: list[MRG32k3a], /) -> RepResult: """Replicate the problem for the supplied model factors. + The default implementation evaluates the DSL-defined objective. It + supports DSL decisions whose names and scalar/vector structure match + the corresponding model factors. Problems with a different coordinate + mapping can override this method. + Args: model_factors (dict): Complete model factors used for the replication. rngs (list[MRG32k3a]): RNGs used to drive the simulation. """ - raise NotImplementedError + problem = self.optimization_problem + values_by_variable_id: dict[int, float] = {} + + for simulation in problem.simulations: + for decision_name, decision in simulation.decisions.items(): + if decision_name not in model_factors: + raise ValueError(f"DSL decision {decision_name!r} has no matching model factor") + + if isinstance(decision, dsl.Variable): + components = (decision,) + values = (model_factors[decision_name],) + else: + components = decision.components + values = tuple(model_factors[decision_name]) + + if len(values) != len(components): + raise ValueError( + f"model factor {decision_name!r} must have " + f"{len(components)} component(s), got {len(values)}" + ) + + for component, value in zip(components, values, strict=True): + values_by_variable_id[id(component)] = float(value) + + missing_variables = [ + variable.name + for variable in problem.variables + if id(variable) not in values_by_variable_id + ] + if missing_variables: + raise ValueError(f"DSL variables are not bound to model factors: {missing_variables!r}") + + decision_vector = tuple( + values_by_variable_id[id(variable)] for variable in problem.variables + ) + result = problem.run_replication(decision_vector, rngs) + return RepResult( + objectives=[ + Objective( + stochastic=result.objective, + stochastic_gradients=result.objective_gradient, + ) + ], + stochastic_constraints=[ + StochasticConstraint( + stochastic=constraint.value, + stochastic_gradients=constraint.gradient, + ) + for constraint in result.stochastic_constraints + ] + or None, + ) def simulate(self, solution: "Solution", num_macroreps: int = 1) -> None: """Simulate `m` i.i.d. replications at solution `x`. diff --git a/test/expected_results/FIXEDSAN1_ADAM.pickle.zst b/test/expected_results/FIXEDSAN1_ADAM.pickle.zst index b47d4eb9..527bfef5 100644 Binary files a/test/expected_results/FIXEDSAN1_ADAM.pickle.zst and b/test/expected_results/FIXEDSAN1_ADAM.pickle.zst differ diff --git a/test/expected_results/FIXEDSAN1_ASTRODF_darwin_arm64.pickle.zst b/test/expected_results/FIXEDSAN1_ASTRODF_darwin_arm64.pickle.zst index f42ef307..4ae0dc9f 100644 Binary files a/test/expected_results/FIXEDSAN1_ASTRODF_darwin_arm64.pickle.zst and b/test/expected_results/FIXEDSAN1_ASTRODF_darwin_arm64.pickle.zst differ diff --git a/test/expected_results/RMITD1_RNDSRCH.pickle.zst b/test/expected_results/RMITD1_RNDSRCH.pickle.zst index 91851896..dd4dce86 100644 Binary files a/test/expected_results/RMITD1_RNDSRCH.pickle.zst and b/test/expected_results/RMITD1_RNDSRCH.pickle.zst differ diff --git a/test/expected_results/SAN1_ADAM.pickle.zst b/test/expected_results/SAN1_ADAM.pickle.zst index 057c2545..7cf3ac0c 100644 Binary files a/test/expected_results/SAN1_ADAM.pickle.zst and b/test/expected_results/SAN1_ADAM.pickle.zst differ