From 5b8cdcd6bf4699062938ab3a0ff36be136ca3c71 Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Thu, 27 Jul 2023 08:50:52 -0700 Subject: [PATCH] [ModelBuilder] big update on the code; remove hash table, use c++ canonicalization of linear expressions; remove code duplication --- .../proto_solver/scip_proto_solver.cc | 6 - ortools/linear_solver/python/model_builder.py | 443 +++++++----------- .../python/model_builder_helper.cc | 34 +- .../python/model_builder_test.py | 198 ++++++-- 4 files changed, 359 insertions(+), 322 deletions(-) diff --git a/ortools/linear_solver/proto_solver/scip_proto_solver.cc b/ortools/linear_solver/proto_solver/scip_proto_solver.cc index 60b35311cd4..0304da7a6ba 100644 --- a/ortools/linear_solver/proto_solver/scip_proto_solver.cc +++ b/ortools/linear_solver/proto_solver/scip_proto_solver.cc @@ -73,13 +73,7 @@ absl::StatusOr ScipInfClamp(SCIP* scip, const double d) { return absl::InvalidArgumentError("Cannot clamp a NaN."); } const double kScipInf = SCIPinfinity(scip); - if (d > -kScipInf && d < kScipInf) { - return d; - } const double clamped = std::clamp(d, -kScipInf, kScipInf); - LOG_EVERY_N_SEC(WARNING, 5) << "A bound was clamped to SCIP's 'infinite' " - "value for safety; this may affect results. " - << "Was: " << d << ", is now: " << clamped; return clamped; } diff --git a/ortools/linear_solver/python/model_builder.py b/ortools/linear_solver/python/model_builder.py index f7b398b5e78..e601031d2bc 100644 --- a/ortools/linear_solver/python/model_builder.py +++ b/ortools/linear_solver/python/model_builder.py @@ -32,12 +32,11 @@ rather than for solving specific optimization problems. """ import abc -import collections import dataclasses import math import numbers import typing -from typing import Callable, Mapping, Optional, Sequence, Union, cast +from typing import Callable, List, Optional, Sequence, Tuple, Union, cast import numpy as np from numpy import typing as npt @@ -61,6 +60,9 @@ SolveStatus = mbh.SolveStatus +# pylint: disable=protected-access + + class LinearExpr(metaclass=abc.ABCMeta): """Holds an linear expression. @@ -126,7 +128,7 @@ def weighted_sum( # pytype: disable=annotation-type-mismatch # numpy-scalars coefficients: Sequence[NumberT], *, constant: NumberT = 0.0, - ) -> Union[NumberT, "_WeightedSum"]: + ) -> Union[NumberT, "_LinearExpression"]: """Creates `sum(expressions[i] * coefficients[i]) + constant`. It can perform simple simplifications and returns different object, @@ -138,7 +140,7 @@ def weighted_sum( # pytype: disable=annotation-type-mismatch # numpy-scalars constant: a numerical constant. Returns: - a _WeightedSum instance or a numerical constant. + a _LinearExpression instance or a numerical constant. """ if len(expressions) != len(coefficients): raise ValueError( @@ -148,46 +150,9 @@ def weighted_sum( # pytype: disable=annotation-type-mismatch # numpy-scalars checked_constant: np.double = mbn.assert_is_a_number(constant) if not expressions: return checked_constant - - # Collect sub-arrays to concatenate. - indices = [] - coeffs = [] - helper = None - for e, c in zip(expressions, coefficients): - if mbn.is_zero(c): - continue - - if mbn.is_a_number(e): - checked_constant += np.double(c * e) - elif isinstance(e, Variable): - if not helper: - helper = e.helper - indices.append(np.array([e.index], dtype=np.int32)) - coeffs.append(np.array([c], dtype=np.double)) - elif isinstance(e, _WeightedSum): - if not helper: - helper = e.helper - checked_constant += np.double(c * e.constant) - indices.append(e.variable_indices) - coeffs.append(e.coefficients * c) - elif isinstance(e, LinearExpr): - expr = _as_flat_linear_expression(e) - # pylint: disable=protected-access - checked_constant += np.double(c * expr._offset) - for variable, coeff in expr._terms.items(): - if not helper: - helper = variable.helper - indices.append(np.array([variable.index], dtype=np.int32)) - coeffs.append(np.array([coeff * c], dtype=np.double)) - - if helper: - return _WeightedSum( - helper=helper, - variable_indices=np.concatenate(indices, axis=0), - coefficients=np.concatenate(coeffs, axis=0), - constant=checked_constant, - ) - return checked_constant + return _sum_as_flat_linear_expression( + to_process=list(zip(expressions, coefficients)), offset=checked_constant + ) @classmethod def term( # pytype: disable=annotation-type-mismatch # numpy-scalars @@ -218,22 +183,10 @@ def term( # pytype: disable=annotation-type-mismatch # numpy-scalars return expression if mbn.is_a_number(expression): return np.double(expression) * checked_coefficient + checked_constant - if isinstance(expression, Variable): - return _WeightedSum( - helper=expression.helper, - variable_indices=np.array([expression.index], dtype=np.int32), - coefficients=np.array([checked_coefficient], dtype=np.double), - constant=checked_constant, - ) - if isinstance(expression, _WeightedSum): - return _WeightedSum( - helper=expression.helper, - variable_indices=np.copy(expression.variable_indices), - coefficients=expression.coefficients * checked_coefficient, - constant=expression.constant * checked_coefficient + checked_constant, - ) if isinstance(expression, LinearExpr): - return expression * checked_coefficient + checked_constant + return _as_flat_linear_expression( + expression * checked_coefficient + checked_constant + ) raise TypeError(f"Unknown expression {expression!r} of type {type(expression)}") def __hash__(self): @@ -280,89 +233,6 @@ def __le__(self, arg: LinearExprT) -> "BoundedLinearExpression": ) # pytype: disable=wrong-arg-types # numpy-scalars -class _WeightedSum(LinearExpr): - """Represents sum(ai * xi) + b.""" - - def __init__( - self, - helper: mbh.ModelBuilderHelper, - *, - variable_indices: npt.NDArray[np.int32], - coefficients: npt.NDArray[np.double], - constant: np.double = np.double(0.0), - ): - super().__init__() - self.__helper: mbh.ModelBuilderHelper = helper - self.__variable_indices: npt.NDArray[np.int32] = variable_indices - self.__coefficients: npt.NDArray[np.double] = mbn.assert_is_a_number_array( - coefficients - ) - self.__constant: np.double = constant - - def multiply_by(self, arg: NumberT) -> LinearExprT: - if mbn.is_zero(arg): - return 0.0 # pytype: disable=bad-return-type # numpy-scalars - if self.__variable_indices.size > 0: - return _WeightedSum( - helper=self.__helper, - variable_indices=np.copy(self.__variable_indices), - coefficients=self.__coefficients * arg, - constant=self.__constant * arg, - ) - else: - return self.constant * arg - - @property - def helper(self) -> mbh.ModelBuilderHelper: - return self.__helper - - @property - def variable_indices(self) -> npt.NDArray[np.int32]: - return self.__variable_indices - - @property - def coefficients(self) -> npt.NDArray[np.double]: - return self.__coefficients - - @property - def constant(self) -> np.double: - return self.__constant - - def pretty_string(self) -> str: - """Pretty print a linear expression into a string.""" - output: str = "" - for index, coeff in zip(self.variable_indices, self.coefficients): - var_name = Variable(self.helper, index, None, None, None).name - - if not output and mbn.is_one(coeff): - output = var_name - elif not output and mbn.is_minus_one(coeff): - output = f"-{var_name}" - elif not output: - output = f"{coeff} * {var_name}" - elif mbn.is_one(coeff): - output += f" + {var_name}" - elif mbn.is_minus_one(coeff): - output += f" - {var_name}" - elif coeff > 0.0: - output += f" + {coeff} * {var_name}" - elif coeff < 0.0: - output += " - {-coeff} * {var_name}" - if self.constant > 0: - output += f" + {self.constant}" - elif self.constant < 0: - output += f" - {-self.constant}" - if not output: - output = "0.0" - return output - - def __repr__(self): - return ( - f"WeightedSum(indices = {self.variable_indices}, coefficients =" - f" {self.coefficients}, constant = {self.constant})" - ) - - class Variable(LinearExpr): """A variable (continuous or integral). @@ -494,11 +364,6 @@ def __eq__(self, arg: Optional[LinearExprT]) -> ConstraintT: def __hash__(self): return hash((self.__helper, self.__index)) - def multiply_by(self, arg: NumberT) -> LinearExprT: - return LinearExpr.weighted_sum( - [self], [arg], constant=0.0 - ) # pytype: disable=wrong-arg-types # numpy-scalars - class _BoundedLinearExpr(metaclass=abc.ABCMeta): """Interface for types that can build bounded linear (boolean) expressions. @@ -525,10 +390,10 @@ def _add_linear_constraint( """ -def _add_linear_constraint( - constraint: Union[bool, _BoundedLinearExpr], +def _add_linear_constraint_to_helper( + bounded_expr: Union[bool, _BoundedLinearExpr], helper: mbh.ModelBuilderHelper, - name: str, + name: Optional[str], ): """Creates a new linear constraint in the helper. @@ -536,7 +401,7 @@ def _add_linear_constraint( BoundedLinearExpressions). Args: - constraint: The constraint to be created. + bounded_expr: The bounded expression used to create the constraint. helper: The helper to create the constraint. name: The name of the constraint to be created. @@ -546,22 +411,23 @@ def _add_linear_constraint( Raises: TypeError: If constraint is an invalid type. """ - if isinstance(constraint, bool): + if isinstance(bounded_expr, bool): c = LinearConstraint(helper) - helper.set_constraint_name(c.index, name) - if constraint: - # constraint that is always feasible: -inf <= nothing <= inf - helper.set_constraint_lower_bound(c.index, -math.inf) - helper.set_constraint_upper_bound(c.index, math.inf) + if name is not None: + helper.set_constraint_name(c.index, name) + if bounded_expr: + # constraint that is always feasible: 0.0 <= nothing <= 0.0 + helper.set_constraint_lower_bound(c.index, 0.0) + helper.set_constraint_upper_bound(c.index, 0.0) else: - # constraint that is always infeasible: -1 <= nothing <= -1 - helper.set_constraint_lower_bound(c.index, -1) + # constraint that is always infeasible: +oo <= nothing <= -oo + helper.set_constraint_lower_bound(c.index, 1) helper.set_constraint_upper_bound(c.index, -1) return c - if isinstance(constraint, _BoundedLinearExpr): + if isinstance(bounded_expr, _BoundedLinearExpr): # pylint: disable=protected-access - return constraint._add_linear_constraint(helper, name) - raise TypeError("invalid type={}".format(type(constraint))) + return bounded_expr._add_linear_constraint(helper, name) + raise TypeError("invalid type={}".format(type(bounded_expr))) @dataclasses.dataclass(repr=False, eq=False, frozen=True) @@ -646,17 +512,19 @@ def __bool__(self) -> bool: ) def _add_linear_constraint( - self, helper: mbh.ModelBuilderHelper, name: str + self, helper: mbh.ModelBuilderHelper, name: Optional[str] ) -> "LinearConstraint": c = LinearConstraint(helper) - expr = _as_flat_linear_expression(self.__expr) + flat_expr = _as_flat_linear_expression(self.__expr) # pylint: disable=protected-access - for variable, coeff in expr._terms.items(): - helper.add_term_to_constraint(c.index, variable.index, coeff) - helper.set_constraint_lower_bound(c.index, self.__lb - expr._offset) - helper.set_constraint_upper_bound(c.index, self.__ub - expr._offset) + helper.add_terms_to_constraint( + c.index, flat_expr._variable_indices, flat_expr._coefficients + ) + helper.set_constraint_lower_bound(c.index, self.__lb - flat_expr._offset) + helper.set_constraint_upper_bound(c.index, self.__ub - flat_expr._offset) # pylint: enable=protected-access - helper.set_constraint_name(c.index, name) + if name is not None: + helper.set_constraint_name(c.index, name) return c @@ -716,18 +584,38 @@ def name(self) -> str: def name(self, name: str) -> None: return self.__helper.set_constraint_name(self.__index, name) + def is_always_false(self) -> bool: + """Returns True if the constraint is always false. + + Usually, it means that it was created by model.add(False) + """ + return self.lower_bound > self.upper_bound + def __str__(self): return self.name def __repr__(self): - return self.__str__() + return ( + f"LinearConstraint({self.name}, lb={self.lower_bound}," + f" ub={self.upper_bound}," + f" var_indices={self.helper.constraint_var_indices(self.index)}," + f" coefficients={self.helper.constraint_coefficients(self.index)})" + ) def set_coefficient(self, var: Variable, coeff: NumberT) -> None: """Sets the coefficient of the variable in the constraint.""" + if self.is_always_false(): + raise ValueError( + f"Constraint {self.index} is always false and cannot be modified" + ) self.__helper.set_constraint_coefficient(self.__index, var.index, coeff) def add_term(self, var: Variable, coeff: NumberT) -> None: """Adds var * coeff to the constraint.""" + if self.is_always_false(): + raise ValueError( + f"Constraint {self.index} is always false and cannot be modified" + ) self.__helper.safe_add_term_to_constraint(self.__index, var.index, coeff) @@ -1165,28 +1053,13 @@ def add_linear_constraint( # pytype: disable=annotation-type-mismatch # numpy- self.__helper.set_constraint_lower_bound(ct.index, lb) self.__helper.set_constraint_upper_bound(ct.index, ub) self.__helper.add_term_to_constraint(ct.index, linear_expr.index, 1.0) - elif isinstance(linear_expr, _WeightedSum): - self.__helper.set_constraint_lower_bound( - ct.index, lb - linear_expr.constant - ) - self.__helper.set_constraint_upper_bound( - ct.index, ub - linear_expr.constant - ) - self.__helper.add_terms_to_constraint( - ct.index, linear_expr.variable_indices, linear_expr.coefficients - ) elif isinstance(linear_expr, LinearExpr): flat_expr = _as_flat_linear_expression(linear_expr) # pylint: disable=protected-access self.__helper.set_constraint_lower_bound(ct.index, lb - flat_expr._offset) self.__helper.set_constraint_upper_bound(ct.index, ub - flat_expr._offset) - variable_indices = [] - coefficients = [] - for variable, coeff in flat_expr._terms.items(): - variable_indices.append(variable.index) - coefficients.append(coeff) self.__helper.add_terms_to_constraint( - ct.index, variable_indices, coefficients + ct.index, flat_expr._variable_indices, flat_expr._coefficients ) else: raise TypeError( @@ -1206,38 +1079,31 @@ def add( Returns: An instance of the `Constraint` class. + + Note that a special treatment is done when the argument does not contain any + variable, and thus evaluates to True or False. + + model.add(True) will create a constraint 0 <= empty sum <= 0 + + model.add(False) will create a constraint inf <= empty sum <= -inf + + you can check the if a constraint is always false (lb=inf, ub=-inf) by + calling LinearConstraint.is_always_false() """ - if isinstance(ct, BoundedLinearExpression): - return self.add_linear_constraint( - ct.expression, ct.lower_bound, ct.upper_bound, name - ) - elif isinstance(ct, VarEqVar): - new_ct = LinearConstraint(self.__helper) - new_ct.lower_bound = 0.0 - new_ct.upper_bound = 0.0 - new_ct.add_term( - ct.left, 1.0 - ) # pytype: disable=wrong-arg-types # numpy-scalars - new_ct.add_term( - ct.right, -1.0 - ) # pytype: disable=wrong-arg-types # numpy-scalars - return new_ct + if isinstance(ct, _BoundedLinearExpr): + return ct._add_linear_constraint(self.__helper, name) + elif isinstance(ct, bool): + return _add_linear_constraint_to_helper(ct, self.__helper, name) elif isinstance(ct, pd.Series): return pd.Series( index=ct.index, data=[ - _add_linear_constraint(expr, self.__helper, f"{name}[{i}]") + _add_linear_constraint_to_helper( + expr, self.__helper, f"{name}[{i}]" + ) for (i, expr) in zip(ct.index, ct) ], ) - elif ct and isinstance(ct, bool): - return self.add_linear_constraint( - linear_expr=0.0 - ) # Evaluate to True. # pytype: disable=wrong-arg-types # numpy-scalars - elif not ct and isinstance(ct, bool): - return self.add_linear_constraint( - 1.0, 0.0, 0.0 - ) # Evaluate to False. # pytype: disable=wrong-arg-types # numpy-scalars else: raise TypeError("Not supported: ModelBuilder.Add(" + str(ct) + ")") @@ -1260,21 +1126,13 @@ def __optimize(self, linear_expr: LinearExprT, maximize: bool) -> None: self.helper.set_objective_offset(linear_expr) elif isinstance(linear_expr, Variable): self.helper.set_var_objective_coefficient(linear_expr.index, 1.0) - elif isinstance(linear_expr, _WeightedSum): - self.helper.set_objective_offset(linear_expr.constant) - self.helper.set_objective_coefficients( - linear_expr.variable_indices, linear_expr.coefficients - ) elif isinstance(linear_expr, LinearExpr): flat_expr = _as_flat_linear_expression(linear_expr) # pylint: disable=protected-access self.helper.set_objective_offset(flat_expr._offset) - variable_indices = [] - coefficients = [] - for variable, coeff in flat_expr._terms.items(): - variable_indices.append(variable.index) - coefficients.append(coeff) - self.helper.set_objective_coefficients(variable_indices, coefficients) + self.helper.set_objective_coefficients( + flat_expr._variable_indices, flat_expr._coefficients + ) else: raise TypeError( f"Not supported: ModelBuilder.minimize/maximize({linear_expr})" @@ -1293,6 +1151,7 @@ def objective_expression(self) -> "_LinearExpression": sum( variable * self.__helper.var_objective_coefficient(variable.index) for variable in self.get_variables() + if self.__helper.var_objective_coefficient(variable.index) != 0.0 ) + self.__helper.objective_offset() ) @@ -1392,20 +1251,12 @@ def value(self, expr: LinearExprT) -> np.double: return expr elif isinstance(expr, Variable): return self.__solve_helper.var_value(expr.index) - elif isinstance(expr, _WeightedSum): - return self.__solve_helper.expression_value( - expr.variable_indices, expr.coefficients, expr.constant - ) elif isinstance(expr, LinearExpr): flat_expr = _as_flat_linear_expression(expr) - variable_indices = [] - coefficients = [] - # pylint: disable=protected-access - for variable, coeff in flat_expr._terms.items(): - variable_indices.append(variable.index) - coefficients.append(coeff) return self.__solve_helper.expression_value( - variable_indices, coefficients, flat_expr._offset + flat_expr._variable_indices, + flat_expr._coefficients, + flat_expr._offset, ) else: raise TypeError(f"Unknown expression {expr!r} of type {type(expr)}") @@ -1533,76 +1384,122 @@ def user_time(self) -> np.double: class _LinearExpression(LinearExpr): """For variables x, an expression: offset + sum_{i in I} coeff_i * x_i.""" - __slots__ = ("_terms", "_offset") + __slots__ = ("_variable_indices", "_coefficients", "_offset", "_helper") - _terms: Mapping["Variable", float] + _variable_indices: npt.NDArray[np.int32] + _coefficients: npt.NDArray[np.double] _offset: float + _helper: Optional[mbh.ModelBuilderHelper] + + @property + def variable_indices(self) -> npt.NDArray[np.int32]: + return self._variable_indices + + @property + def coefficients(self) -> npt.NDArray[np.double]: + return self._coefficients + + @property + def constant(self) -> float: + return self._offset + + @property + def helper(self) -> Optional[mbh.ModelBuilderHelper]: + return self._helper def __repr__(self): return self.__str__() def __str__(self): + if self._helper is None: + return str(self._offset) + result = [] - if self._offset != 0.0: - result.append(str(self._offset)) - sorted_keys = sorted(self._terms.keys(), key=str) - num_displayed_terms = 0 - for variable in sorted_keys: - if num_displayed_terms == _MAX_LINEAR_EXPRESSION_REPR_TERMS: + for index, coeff in zip(self.variable_indices, self.coefficients): + if len(result) >= _MAX_LINEAR_EXPRESSION_REPR_TERMS: result.append(" + ...") break - coefficient = self._terms[variable] - if coefficient == 0.0: - continue - if result: - if coefficient > 0: - result.append(" + ") - else: - result.append(" - ") - elif coefficient < 0: - result.append("- ") - if abs(coefficient) != 1.0: - result.append(f"{abs(coefficient)} * ") - result.append(f"{variable}") - num_displayed_terms += 1 + var_name = Variable(self._helper, index, None, None, None).name + if not result and mbn.is_one(coeff): + result.append(var_name) + elif not result and mbn.is_minus_one(coeff): + result.append(f"-{var_name}") + elif not result: + result.append(f"{coeff} * {var_name}") + elif mbn.is_one(coeff): + result.append(f" + {var_name}") + elif mbn.is_minus_one(coeff): + result.append(f" - {var_name}") + elif coeff > 0.0: + result.append(f" + {coeff} * {var_name}") + elif coeff < 0.0: + result.append(f" - {-coeff} * {var_name}") + + if not result: + return f"{self.constant}" + if self.constant > 0: + result.append(f" + {self.constant}") + elif self.constant < 0: + result.append(f" - {-self.constant}") return "".join(result) -def _as_flat_linear_expression(base_expr: LinearExprT) -> _LinearExpression: - """Converts floats, ints and Linear objects to a LinearExpression.""" - # pylint: disable=protected-access - if isinstance(base_expr, _LinearExpression): - return base_expr - terms = collections.defaultdict(lambda: 0.0) - offset: float = 0.0 - to_process = [(base_expr, 1.0)] +def _sum_as_flat_linear_expression( + to_process: List[Tuple[LinearExprT, float]], offset: float = 0.0 +) -> _LinearExpression: + """Creates a _LinearExpression as the sum of terms.""" + indices = [] + coeffs = [] + helper = None while to_process: # Flatten AST of LinearTypes. expr, coeff = to_process.pop() if isinstance(expr, _Sum): to_process.append((expr._left, coeff)) to_process.append((expr._right, coeff)) elif isinstance(expr, Variable): - terms[expr] += coeff + indices.append([expr.index]) + coeffs.append([coeff]) + if helper is None: + helper = expr.helper elif mbn.is_a_number(expr): offset += coeff * cast(NumberT, expr) elif isinstance(expr, _Product): to_process.append((expr._expression, coeff * expr._coefficient)) elif isinstance(expr, _LinearExpression): offset += coeff * expr._offset - for variable, variable_coefficient in expr._terms.items(): - terms[variable] += coeff * variable_coefficient - elif isinstance(expr, _WeightedSum): - offset += coeff * expr.constant - for variable_index, variable_coefficient in zip( - expr.variable_indices, expr.coefficients - ): - variable = Variable(expr.helper, variable_index, None, None, None) - terms[variable] += coeff * variable_coefficient + if expr._helper is not None: + indices.append(expr.variable_indices) + coeffs.append(np.multiply(expr.coefficients, coeff)) + if helper is None: + helper = expr._helper else: raise TypeError( "Unrecognized linear expression: " + str(expr) + f" {type(expr)}" ) - return _LinearExpression(terms, offset) + + if helper is not None: + all_indices: npt.NDArray[np.int32] = np.concatenate(indices, axis=0) + all_coeffs: npt.NDArray[np.double] = np.concatenate(coeffs, axis=0) + sorted_indices, sorted_coefficients = helper.sort_and_regroup_terms( + all_indices, all_coeffs + ) + return _LinearExpression(sorted_indices, sorted_coefficients, offset, helper) + else: + assert not indices + assert not coeffs + return _LinearExpression( + _variable_indices=np.zeros(dtype=np.int32, shape=[0]), + _coefficients=np.zeros(dtype=np.double, shape=[0]), + _offset=offset, + _helper=None, + ) + + +def _as_flat_linear_expression(base_expr: LinearExprT) -> _LinearExpression: + """Converts floats, ints and Linear objects to a LinearExpression.""" + if isinstance(base_expr, _LinearExpression): + return base_expr + return _sum_as_flat_linear_expression(to_process=[(base_expr, 1.0)], offset=0.0) @dataclasses.dataclass(repr=False, eq=False, frozen=True) diff --git a/ortools/linear_solver/python/model_builder_helper.cc b/ortools/linear_solver/python/model_builder_helper.cc index b270a64a777..190707af272 100644 --- a/ortools/linear_solver/python/model_builder_helper.cc +++ b/ortools/linear_solver/python/model_builder_helper.cc @@ -121,7 +121,7 @@ void BuildModelFromSparseData( } std::vector> SortedGroupedTerms( - const std::vector& indices, const std::vector& coefficients) { + absl::Span indices, absl::Span coefficients) { CHECK_EQ(indices.size(), coefficients.size()); std::vector> terms; terms.reserve(indices.size()); @@ -139,14 +139,15 @@ std::vector> SortedGroupedTerms( }); int pos = 0; for (int i = 0; i < terms.size(); ++i) { - if (i == 0 || terms[i].first != terms[i - 1].first) { - if (i != pos) { - terms[pos] = terms[i]; - } - pos++; - } else { - terms[pos].second += terms[i].second; + const int var = terms[i].first; + double coeff = terms[i].second; + while (i + 1 < terms.size() && terms[i + 1].first == var) { + coeff += terms[i + 1].second; + ++i; } + if (coeff == 0.0) continue; + terms[pos] = {var, coeff}; + ++pos; } terms.resize(pos); return terms; @@ -341,7 +342,22 @@ PYBIND11_MODULE(model_builder_helper, m) { .def("set_maximize", &ModelBuilderHelper::SetMaximize, arg("maximize")) .def("set_objective_offset", &ModelBuilderHelper::SetObjectiveOffset, arg("offset")) - .def("objective_offset", &ModelBuilderHelper::ObjectiveOffset); + .def("objective_offset", &ModelBuilderHelper::ObjectiveOffset) + .def("sort_and_regroup_terms", + [](ModelBuilderHelper* helper, py::array_t indices, + py::array_t coefficients) { + const std::vector> terms = + SortedGroupedTerms(indices, coefficients); + std::vector sorted_indices; + std::vector sorted_coefficients; + sorted_indices.reserve(terms.size()); + sorted_coefficients.reserve(terms.size()); + for (const auto& [i, c] : terms) { + sorted_indices.push_back(i); + sorted_coefficients.push_back(c); + } + return std::make_pair(sorted_indices, sorted_coefficients); + }); py::enum_(m, "SolveStatus") .value("OPTIMAL", SolveStatus::OPTIMAL) diff --git a/ortools/linear_solver/python/model_builder_test.py b/ortools/linear_solver/python/model_builder_test.py index f8ee5a2858f..2c7ae548c3e 100644 --- a/ortools/linear_solver/python/model_builder_test.py +++ b/ortools/linear_solver/python/model_builder_test.py @@ -30,6 +30,16 @@ from ortools.linear_solver.python import model_builder as mb +def build_dict(expr: mb.LinearExprT) -> Dict[mb.Variable, float]: + res = {} + for i, c in zip(expr.variable_indices, expr.coefficients): + if not c: + continue + var = mb.Variable(expr.helper, lb=i, ub=None, is_integral=None, name=None) + res[var] = c + return res + + class ModelBuilderTest(absltest.TestCase): # Number of decimal places to use for numerical tolerance for # checking primal, dual, objective values and other values. @@ -204,7 +214,7 @@ def test_class_api(self): np.array([1, 1, 1], dtype=np.double), e1.coefficients ) self.assertEqual(e1.constant, 0.0) - self.assertEqual(e1.pretty_string(), "x + y + z") + self.assertEqual(e1.__str__(), "x + y + z") e2 = mb.LinearExpr.sum([e1, 4.0]) np_testing.assert_array_equal(expected_vars, e2.variable_indices) @@ -212,7 +222,7 @@ def test_class_api(self): np.array([1, 1, 1], dtype=np.double), e2.coefficients ) self.assertEqual(e2.constant, 4.0) - self.assertEqual(e2.pretty_string(), "x + y + z + 4.0") + self.assertEqual(e2.__str__(), "x + y + z + 4.0") e3 = mb.LinearExpr.term(e2, 2) np_testing.assert_array_equal(expected_vars, e3.variable_indices) @@ -220,7 +230,7 @@ def test_class_api(self): np.array([2, 2, 2], dtype=np.double), e3.coefficients ) self.assertEqual(e3.constant, 8.0) - self.assertEqual(e3.pretty_string(), "2.0 * x + 2.0 * y + 2.0 * z + 8.0") + self.assertEqual(e3.__str__(), "2.0 * x + 2.0 * y + 2.0 * z + 8.0") e4 = mb.LinearExpr.weighted_sum([x, t], [-1, 1], constant=2) np_testing.assert_array_equal( @@ -230,7 +240,7 @@ def test_class_api(self): np.array([-1, 1], dtype=np.double), e4.coefficients ) self.assertEqual(e4.constant, 2.0) - self.assertEqual(e4.pretty_string(), "-x + t + 2.0") + self.assertEqual(e4.__str__(), "-x + t + 2.0") e4b = mb.LinearExpr.weighted_sum([e4 * 3], [1]) np_testing.assert_array_equal( @@ -240,17 +250,17 @@ def test_class_api(self): np.array([-3, 3], dtype=np.double), e4b.coefficients ) self.assertEqual(e4b.constant, 6.0) - self.assertEqual(e4b.pretty_string(), "-3.0 * x + 3.0 * t + 6.0") + self.assertEqual(e4b.__str__(), "-3.0 * x + 3.0 * t + 6.0") e5 = mb.LinearExpr.sum([e1, -3, e4]) np_testing.assert_array_equal( - np.array([0, 1, 2, 0, 3], dtype=np.int32), e5.variable_indices + np.array([1, 2, 3], dtype=np.int32), e5.variable_indices ) np_testing.assert_array_equal( - np.array([1, 1, 1, -1, 1], dtype=np.double), e5.coefficients + np.array([1, 1, 1], dtype=np.double), e5.coefficients ) self.assertEqual(e5.constant, -1.0) - self.assertEqual(e5.pretty_string(), "x + y + z - x + t - 1.0") + self.assertEqual(e5.__str__(), "y + z + t - 1.0") e6 = mb.LinearExpr.term(x, 2.0, constant=1.0) np_testing.assert_array_equal( @@ -375,6 +385,40 @@ def test_vareqvar(self): self.assertEqual(ct.left.index, x.index) self.assertEqual(ct.right.index, y.index) + def test_create_false_ct(self): + # Create the model. + model = mb.ModelBuilder() + x = model.new_num_var(0.0, math.inf, "x") + + ct = model.add(False) + self.assertTrue(ct.is_always_false()) + self.assertRaises(ValueError, ct.add_term, x, 1) + + model.maximize(x) + + solver = mb.ModelSolver("glop") + status = solver.solve(model) + self.assertEqual(status, mb.SolveStatus.INFEASIBLE) + + def test_create_true_ct(self): + # Create the model. + model = mb.ModelBuilder() + x = model.new_num_var(0.0, 5.0, "x") + + ct = model.add(True) + self.assertEqual(ct.lower_bound, 0.0) + self.assertEqual(ct.upper_bound, 0.0) + ct.add_term(var=x, coeff=1) + + model.maximize(x) + + solver = mb.ModelSolver("glop") + status = solver.solve(model) + + self.assertEqual(status, mb.SolveStatus.OPTIMAL) + # Note that ct is binding. + self.assertEqual(0.0, solver.objective_value) + class InternalHelperTest(absltest.TestCase): def test_anonymous_variables(self): @@ -433,22 +477,22 @@ def setUp(self): dict( testcase_name="x[0] + 5", expr=lambda x, y: x[0] + 5, - expected_repr="5.0 + x[0]", + expected_repr="x[0] + 5.0", ), dict( testcase_name="x[0] - 5", expr=lambda x, y: x[0] - 5, - expected_repr="-5.0 + x[0]", + expected_repr="x[0] - 5.0", ), dict( testcase_name="5 - x[0]", expr=lambda x, y: 5 - x[0], - expected_repr="5.0 - x[0]", + expected_repr="-x[0] + 5.0", ), dict( testcase_name="5 + x[0]", expr=lambda x, y: 5 + x[0], - expected_repr="5.0 + x[0]", + expected_repr="x[0] + 5.0", ), dict( testcase_name="x[0] + y[0]", @@ -458,12 +502,12 @@ def setUp(self): dict( testcase_name="x[0] + y[0] + 5", expr=lambda x, y: x[0] + y[0] + 5, - expected_repr="5.0 + x[0] + y[0]", + expected_repr="x[0] + y[0] + 5.0", ), dict( testcase_name="5 + x[0] + y[0]", expr=lambda x, y: 5 + x[0] + y[0], - expected_repr="5.0 + x[0] + y[0]", + expected_repr="x[0] + y[0] + 5.0", ), dict( testcase_name="5 + x[0] - x[0]", @@ -473,7 +517,7 @@ def setUp(self): dict( testcase_name="5 + x[0] - y[0]", expr=lambda x, y: 5 + x[0] - y[0], - expected_repr="5.0 + x[0] - y[0]", + expected_repr="x[0] - y[0] + 5.0", ), dict( testcase_name="x.sum()", @@ -483,18 +527,23 @@ def setUp(self): dict( testcase_name="x.add(y, fill_value=0).sum() + 5", expr=lambda x, y: x.add(y, fill_value=0).sum() + 5, - expected_repr="5.0 + x[0] + x[1] + x[2] + y[0] + y[1] + ...", + expected_repr="x[0] + x[1] + x[2] + y[0] + y[1] + ... + 5.0", + ), + dict( + testcase_name="sum(x, y + 5)", + expr=lambda x, y: mb.LinearExpr.sum([x.sum(), y.sum() + 5]), + expected_repr="x[0] + x[1] + x[2] + y[0] + y[1] + ... + 5.0", ), # Product dict( testcase_name="- x.sum()", expr=lambda x, y: -x.sum(), - expected_repr="- x[0] - x[1] - x[2]", + expected_repr="-x[0] - x[1] - x[2]", ), dict( testcase_name="5 - x.sum()", expr=lambda x, y: 5 - x.sum(), - expected_repr="5.0 - x[0] - x[1] - x[2]", + expected_repr="-x[0] - x[1] - x[2] + 5.0", ), dict( testcase_name="x.sum() / 2.0", @@ -935,8 +984,8 @@ class ModelBuilderLinearConstraintsTest(parameterized.TestCase): name="true", bounded_exprs=lambda x, y: True, constraint_count=1, - lower_bounds=[-math.inf], - upper_bounds=[math.inf], + lower_bounds=[0.0], + upper_bounds=[0.0], expression_terms=lambda x, y: [{}], expression_offsets=[0], ), @@ -945,8 +994,8 @@ class ModelBuilderLinearConstraintsTest(parameterized.TestCase): name="true", bounded_exprs=lambda x, y: pd.Series(True), constraint_count=1, - lower_bounds=[-math.inf], - upper_bounds=[math.inf], + lower_bounds=[0.0], + upper_bounds=[0.0], expression_terms=lambda x, y: [{}], expression_offsets=[0], ), @@ -955,8 +1004,8 @@ class ModelBuilderLinearConstraintsTest(parameterized.TestCase): name="false", bounded_exprs=lambda x, y: False, constraint_count=1, - lower_bounds=[-1], - upper_bounds=[-1], + lower_bounds=[1.0], + upper_bounds=[-1.0], expression_terms=lambda x, y: [{}], expression_offsets=[0], ), @@ -965,8 +1014,8 @@ class ModelBuilderLinearConstraintsTest(parameterized.TestCase): name="false", bounded_exprs=lambda x, y: pd.Series(False), constraint_count=1, - lower_bounds=[-1], - upper_bounds=[-1], + lower_bounds=[1.0], + upper_bounds=[-1.0], expression_terms=lambda x, y: [{}], expression_offsets=[0], ), @@ -1047,7 +1096,7 @@ class ModelBuilderLinearConstraintsTest(parameterized.TestCase): constraint_count=1, lower_bounds=[0], upper_bounds=[0], - expression_terms=lambda x, y: [{x[0]: 0}], + expression_terms=lambda x, y: [{}], expression_offsets=[0], ), dict( @@ -1057,7 +1106,7 @@ class ModelBuilderLinearConstraintsTest(parameterized.TestCase): constraint_count=1, lower_bounds=[0], upper_bounds=[0], - expression_terms=lambda x, y: [{x[0]: 0}], + expression_terms=lambda x, y: [{}], expression_offsets=[0], ), dict( @@ -1384,7 +1433,7 @@ def test_get_linear_constraint_expressions( expr_terms = expression_terms(x, y) self.assertLen(linear_constraint_expressions, len(expr_terms)) for expr, expr_term in zip(linear_constraint_expressions, expr_terms): - self.assertDictEqual(expr._terms, expr_term) + self.assertDictEqual(build_dict(expr), expr_term) self.assertSequenceAlmostEqual( [expr._offset for expr in linear_constraint_expressions], expression_offsets, @@ -1423,10 +1472,14 @@ def assertLinearExpressionAlmostEqual( expr2: mb._LinearExpression, ) -> None: """Test that the two linear expressions are almost equal.""" - for variable, coeff in expr1._terms.items(): - self.assertAlmostEqual(expr2._terms.get(variable, 0), coeff) - for variable, coeff in expr2._terms.items(): - self.assertAlmostEqual(expr1._terms.get(variable, 0), coeff) + self.assertEqual(len(expr1.variable_indices), len(expr2.variable_indices)) + if expr1.variable_indices: + self.assertSequenceEqual(expr1.variable_indices, expr2.variable_indices) + self.assertSequenceAlmostEqual( + expr1.coefficients, expr2.coefficients, places=5 + ) + else: + self.assertEmpty(expr2.coefficients) self.assertAlmostEqual(expr1._offset, expr2._offset) @parameterized.product( @@ -1462,7 +1515,7 @@ def test_set_new_objective(self): # Set and check for old objective. model.maximize(old_objective_expression) got_objective_expression = model.objective_expression() - for var_coeff in got_objective_expression._terms.values(): + for var_coeff in got_objective_expression.coefficients: self.assertAlmostEqual(var_coeff, 0) self.assertAlmostEqual(got_objective_expression._offset, 1) @@ -1859,6 +1912,8 @@ def test_get_objective_value( else: self.assertAlmostEqual(model_solver.objective_value, 0) + +class ModelBuilderExamplesTest(absltest.TestCase): def test_simple_problem(self): # max 5x1 + 4x2 + 3x3 # s.t 2x1 + 3x2 + x3 <= 5 @@ -1902,6 +1957,81 @@ def test_simple_problem(self): self.assertAlmostEqual(0, solver.dual_value((x[1]))) self.assertAlmostEqual(1, solver.dual_value((x[2]))) + def test_graph_k_color(self): + # Assign a color to each vertex of graph, s.t that no two adjacent + # vertices share the same color. Assume N vertices, and max number of + # colors is K + # Consider graph with edges: + # Edge 1: (0,1) + # Edge 2: (0,2) + # Edge 3: (1,3) + # Trying to color graph with at most 3 colors (0, 1, 2) + # Two sets of variables: + # x - pandas series representing coloring status of nodes + # y - pandas series indicating whether color has been used + # Min: y0 + y1 + y2 + # s.t: Every vertex must be assigned exactly one color + # if two vertices are adjacent they cannot have the same color + + model = mb.ModelBuilder() + num_colors = 3 + num_nodes = 4 + x = model.new_var_series( + "x", + pd.MultiIndex.from_product( + (range(num_nodes), range(num_colors)), + names=["node", "color"], + ), + lower_bounds=0, + upper_bounds=1, + is_integral=True, + ) + y = model.new_var_series( + "y", + pd.Index(range(num_colors), name="color"), + lower_bounds=0, + upper_bounds=1, + is_integral=True, + ) + model.minimize(y.dot([1, 1, 1])) + # Every vertex must be assigned exactly one color + model.add(x.groupby("node").sum().apply(lambda expr: expr == 1)) + # If a vertex i is assigned to a color j then color j is used: + # namely, we re-arrange the terms to express: "x - y <= 0". + model.add(x.sub(y, fill_value=0).apply(lambda expr: expr <= 0)) + # if two vertices are adjacent they cannot have the same color j + for j in range(num_colors): + model.add(x[0, j] + x[1, j] <= 1) + model.add(x[0, j] + x[2, j] <= 1) + model.add(x[1, j] + x[3, j] <= 1) + solver = mb.ModelSolver("sat") + run = solver.solve(model) + self.assertEqual(run, mb.SolveStatus.OPTIMAL) + self.assertEqual(solver.objective_value, 2) + + def test_knapsack_problem(self): + # Basic Knapsack Problem: Given N items, + # with N different weights and values, find the maximum possible value of + # the Items while meeting a weight requirement + # We have 3 items: + # Item x1: Weight = 10, Value = 60 + # Item x2: Weight = 20, Value = 100 + # Item x3: Weight = 30, Value = 120 + # Max: 60x1 + 100x2 + 120x3 + # s.t: 10x1 + 20x2 + 30x3 <= 50 + model = mb.ModelBuilder() + x = model.new_bool_var_series("x", pd.Index(range(3))) + self.assertLen(model.get_variables(), 3) + model.maximize(x.dot([60, 100, 120])) + model.add(x.dot([10, 20, 30]) <= 50) + self.assertLen(model.get_linear_constraints(), 1) + solver = mb.ModelSolver("sat") + run = solver.solve(model) + self.assertEqual(run, mb.SolveStatus.OPTIMAL) + i = solver.values(model.get_variables()) + self.assertEqual(solver.objective_value, 220) + self.assertSequenceAlmostEqual(i, [0, 1, 1]) + if __name__ == "__main__": absltest.main()