diff --git a/ortools/sat/BUILD.bazel b/ortools/sat/BUILD.bazel index c1a151b222..7ba3031313 100644 --- a/ortools/sat/BUILD.bazel +++ b/ortools/sat/BUILD.bazel @@ -205,7 +205,8 @@ cc_library( ":util", "//ortools/algorithms:binary_search", "//ortools/util:sorted_interval_list", - "@com_google_absl//absl/log", + "@com_google_absl//absl/functional:bind_front", + "@com_google_absl//absl/log", "@com_google_absl//absl/random", "@com_google_absl//absl/strings", "@com_google_absl//absl/time", diff --git a/ortools/sat/constraint_violation.cc b/ortools/sat/constraint_violation.cc index e8b8d68786..00f0fc5be6 100644 --- a/ortools/sat/constraint_violation.cc +++ b/ortools/sat/constraint_violation.cc @@ -1284,11 +1284,13 @@ void AddCircuitFlowConstraints(LinearIncrementalEvaluator& linear_evaluator, LsEvaluator::LsEvaluator(const CpModelProto& cp_model) : cp_model_(cp_model) { var_to_constraints_.resize(cp_model_.variables_size()); jump_value_optimal_.resize(cp_model_.variables_size(), true); + num_violated_constraint_per_var_.assign(cp_model_.variables_size(), 0); std::vector ignored_constraints(cp_model_.constraints_size(), false); std::vector additional_constraints; CompileConstraintsAndObjective(ignored_constraints, additional_constraints); BuildVarConstraintGraph(); + pos_in_violated_constraints_.assign(NumEvaluatorConstraints(), -1); } LsEvaluator::LsEvaluator( @@ -1297,8 +1299,10 @@ LsEvaluator::LsEvaluator( : cp_model_(cp_model) { var_to_constraints_.resize(cp_model_.variables_size()); jump_value_optimal_.resize(cp_model_.variables_size(), true); + num_violated_constraint_per_var_.assign(cp_model_.variables_size(), 0); CompileConstraintsAndObjective(ignored_constraints, additional_constraints); BuildVarConstraintGraph(); + pos_in_violated_constraints_.assign(NumEvaluatorConstraints(), -1); } void LsEvaluator::BuildVarConstraintGraph() { @@ -1497,7 +1501,11 @@ void LsEvaluator::CompileConstraintsAndObjective( bool LsEvaluator::ReduceObjectiveBounds(int64_t lb, int64_t ub) { if (!cp_model_.has_objective()) return false; - return linear_evaluator_.ReduceBounds(/*c=*/0, lb, ub); + if (linear_evaluator_.ReduceBounds(/*c=*/0, lb, ub)) { + UpdateViolatedList(0); + return true; + } + return false; } void LsEvaluator::OverwriteCurrentSolution(absl::Span solution) { @@ -1667,27 +1675,6 @@ double LsEvaluator::WeightedViolationDelta(absl::Span weights, return linear_evaluator_.WeightedViolationDelta(weights, var, delta) + WeightedNonLinearViolationDelta(weights, var, delta); } -// TODO(user): Speed-up: -// - Use a row oriented representation of the model (could reuse the Apply -// methods on proto). -// - Maintain the list of violated constraints ? -std::vector LsEvaluator::VariablesInViolatedConstraints() const { - std::vector variables; - for (int var = 0; var < cp_model_.variables_size(); ++var) { - if (linear_evaluator_.AppearsInViolatedConstraints(var)) { - variables.push_back(var); - } else { - for (const int ct_index : var_to_constraints_[var]) { - if (constraints_[ct_index]->violation() > 0) { - variables.push_back(var); - break; - } - } - } - } - - return variables; -} bool LsEvaluator::VariableOnlyInLinearConstraintWithConvexViolationChange( int var) const { @@ -1695,6 +1682,7 @@ bool LsEvaluator::VariableOnlyInLinearConstraintWithConvexViolationChange( } void LsEvaluator::RecomputeViolatedList(bool linear_only) { + num_violated_constraint_per_var_.assign(cp_model_.variables_size(), 0); violated_constraints_.clear(); pos_in_violated_constraints_.assign(NumEvaluatorConstraints(), -1); @@ -1713,6 +1701,9 @@ void LsEvaluator::UpdateViolatedList(const int c) { if (pos != -1) return; pos_in_violated_constraints_[c] = violated_constraints_.size(); violated_constraints_.push_back(c); + for (const int v : ConstraintToVars(c)) { + num_violated_constraint_per_var_[v] += 1; + } return; } @@ -1723,6 +1714,9 @@ void LsEvaluator::UpdateViolatedList(const int c) { violated_constraints_[pos] = last; pos_in_violated_constraints_[c] = -1; violated_constraints_.pop_back(); + for (const int v : ConstraintToVars(c)) { + num_violated_constraint_per_var_[v] -= 1; + } } } // namespace sat diff --git a/ortools/sat/constraint_violation.h b/ortools/sat/constraint_violation.h index f850050a4f..80dbb16c35 100644 --- a/ortools/sat/constraint_violation.h +++ b/ortools/sat/constraint_violation.h @@ -359,19 +359,19 @@ class LsEvaluator { return &linear_evaluator_; } - // Returns the set of variables appearing in a violated constraint. - std::vector VariablesInViolatedConstraints() const; - // List of the currently violated constraints. // - It is initialized by RecomputeViolatedList() // - And incrementally maintained by UpdateVariableValue() // - // The order depends on the algorithm used and shouldn't be + // The order depends on the algorithm used and shouldn't be relied on. void RecomputeViolatedList(bool linear_only); const std::vector& ViolatedConstraints() const { return violated_constraints_; } - + // Returns the number of constraints in ViolatedConstraints containing `var`. + int NumViolatedConstraintsForVar(int var) const { + return num_violated_constraint_per_var_[var]; + } // Indicates if the computed jump value is always the best choice. bool VariableOnlyInLinearConstraintWithConvexViolationChange(int var) const; @@ -435,6 +435,7 @@ class LsEvaluator { // - violated_constraints_[pos_in_violated_constraints_[c]] = c std::vector pos_in_violated_constraints_; std::vector violated_constraints_; + std::vector num_violated_constraint_per_var_; // Constraint index and violation delta for the last update. std::vector> last_update_violation_changes_; diff --git a/ortools/sat/docs/README.md b/ortools/sat/docs/README.md index aed756a545..4f15cbbd24 100644 --- a/ortools/sat/docs/README.md +++ b/ortools/sat/docs/README.md @@ -1,13 +1,9 @@ -| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) | -| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- | - +[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) +----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | --------------------------------------------------------------------------------------- # Using the CP-SAT solver https://developers.google.com/optimization/cp/cp_solver - - - ## Documentation structure This document presents modeling recipes for the CP-SAT solver. diff --git a/ortools/sat/docs/boolean_logic.md b/ortools/sat/docs/boolean_logic.md index ab2b938f65..e581c04be6 100644 --- a/ortools/sat/docs/boolean_logic.md +++ b/ortools/sat/docs/boolean_logic.md @@ -1,13 +1,9 @@ -| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) | -| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- | - +[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) +----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | --------------------------------------------------------------------------------------- # Boolean logic recipes for the CP-SAT solver. https://developers.google.com/optimization/ - - - ## Introduction The CP-SAT solver can express Boolean variables and constraints. A **Boolean diff --git a/ortools/sat/docs/channeling.md b/ortools/sat/docs/channeling.md index a00049c049..7dcf50f3f2 100644 --- a/ortools/sat/docs/channeling.md +++ b/ortools/sat/docs/channeling.md @@ -1,12 +1,9 @@ -| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) | -| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- | - +[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) +----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | --------------------------------------------------------------------------------------- # Channeling constraints https://developers.google.com/optimization/ - - A *channeling constraint* links variables inside a model. They're used when you want to express a complicated relationship between variables, such as "if this variable satisfies a condition, force another variable to a particular value". diff --git a/ortools/sat/docs/integer_arithmetic.md b/ortools/sat/docs/integer_arithmetic.md index 63c9ecdc47..11967ff53c 100644 --- a/ortools/sat/docs/integer_arithmetic.md +++ b/ortools/sat/docs/integer_arithmetic.md @@ -1,12 +1,9 @@ -| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) | -| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- | - +[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) +----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | --------------------------------------------------------------------------------------- # Integer arithmetic recipes for the CP-SAT solver. https://developers.google.com/optimization/ - - ## Introduction The CP-SAT solver can express integer variables and constraints. diff --git a/ortools/sat/docs/model.md b/ortools/sat/docs/model.md index 62fd2ff9ec..9cfde48e56 100644 --- a/ortools/sat/docs/model.md +++ b/ortools/sat/docs/model.md @@ -1,12 +1,9 @@ -| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) | -| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- | - +[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) +----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | --------------------------------------------------------------------------------------- # Model manipulation https://developers.google.com/optimization/ - - ## Introduction In all languages, the CpModel class is a thin wrapper around a diff --git a/ortools/sat/docs/scheduling.md b/ortools/sat/docs/scheduling.md index 2d412b0912..a7ee164bbf 100644 --- a/ortools/sat/docs/scheduling.md +++ b/ortools/sat/docs/scheduling.md @@ -1,12 +1,9 @@ -| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) | -| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- | - +[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) +----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | --------------------------------------------------------------------------------------- # Scheduling recipes for the CP-SAT solver. https://developers.google.com/optimization/ - - ## Introduction Scheduling in Operations Research involves problems of tasks, resources and @@ -754,7 +751,7 @@ def main(): index=tasks_df.index, starts=starts, sizes=tasks_df.duration, - performed_literals=performed, + are_present=performed, ) # Set up the profile. We use fixed (intervals, demands) to fill in the space @@ -795,7 +792,7 @@ def main(): else: print(f"task {task} is not performed") elif status == cp_model.INFEASIBLE: - print("The problem is infeasible") + print("No solution found") else: print("Something is wrong, check the status and the log of the solve") @@ -808,7 +805,7 @@ if __name__ == "__main__": ## Ranking tasks in a disjunctive resource -To rank intervals in a NoOverlap constraint, we will count the number of +To rank intervals in a no_overlap constraint, we will count the number of performed intervals that precede each interval. This is slightly complicated if some interval variables are optional. To @@ -1455,8 +1452,8 @@ public class RankingSampleSat ## Ranking tasks in a disjunctive resource with a circuit constraint. -To rank intervals in a NoOverlap constraint, we will use a circuit constraint to -perform the transitive reduction from precedences to successors. +To rank intervals in a no_overlap constraint, we will use a circuit constraint +to perform the transitive reduction from precedences to successors. This is slightly complicated if some interval variables are optional, and we need to take into account the case where no task is performed. @@ -1728,10 +1725,10 @@ SchedulingWithCalendarSampleSat() ## Detecting if two intervals overlap. -We want a Boolean variable to be true iff two intervals overlap. To enforce -this, we will create 3 Boolean variables, link two of them to the relative -positions of the two intervals, and define the third one using the other two -Boolean variables. +We want a Boolean variable to be true if and only if two intervals overlap. To +enforce this, we will create 3 Boolean variables, link two of them to the +relative positions of the two intervals, and define the third one using the +other two Boolean variables. There are two ways of linking the three Boolean variables. The first version uses one clause and two implications. Propagation will be faster using this @@ -1855,5 +1852,3 @@ OverlappingIntervals() ## Precedences between intervals ## Convex hull of a set of intervals - -## Reservoir constraint diff --git a/ortools/sat/docs/solver.md b/ortools/sat/docs/solver.md index aaf6f991c7..019d362141 100644 --- a/ortools/sat/docs/solver.md +++ b/ortools/sat/docs/solver.md @@ -1,12 +1,9 @@ -| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) | -| ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- | - +[home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Troubleshooting](troubleshooting.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) +----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | ------------------------------------- | --------------------------------------------------------------------------------------- # Solving a CP-SAT model https://developers.google.com/optimization/ - - ## Changing the parameters of the solver The SatParameters protobuf encapsulates the set of parameters of a CP-SAT diff --git a/ortools/sat/docs/troubleshooting.md b/ortools/sat/docs/troubleshooting.md new file mode 100644 index 0000000000..04dc6e36eb --- /dev/null +++ b/ortools/sat/docs/troubleshooting.md @@ -0,0 +1,268 @@ +# Troubleshooting + +## Enable logging + +CP-SAT supports extensive logging of its internals. + +To enable it, just set the parameter `log_search_progress` to true. + +A good explanation of the log can be found in the +[cpsat-primer](https://github.com/d-krupke/cpsat-primer/blob/main/understanding_the_log.md). + +## Improving performance with multiple workers + +CP-SAT is built with parallelism in mind. While single worker is a good +solution, the solver is tuned to reached its best performance with multiple +workers. + +We distinguish multiple thresholds: + +- **[8 workers]** This is the minimal setting for parallel search. It blends + workers with different linear relaxations (none, default, maximal), core + based search if applicable, a quick_restart subsolver, a feasibility_jump + first solution subsolver, and dedicated Large Neighborhood Search subsolvers + to find improving solutions. +- **[16 workers]** Bumping to 16 workers adds a continuous probing subsolver, + more first solution subsolvers (random search and feasibility_jump), two + dual subsolvers dedicated to improving the lower bound of the objective + (when minimizing), and more LNS subsolvers. +- **[24 workers]** Adds more dual subsolvers, no_lp and max_lp variations of + previous subsolvers, more feasibility_jump first solution subsolvers, and + more LNS subsolvers. +- **[32 workers and more]** Adds more LNS workers, and more first solution + subsolvers. + +## Solve status + +Solving a problem yields the following possible status (CpSolverStatus): + +- **[UNKNOWN]** The status of the model is still unknown. A search limit has + been reached before any of the statuses below could be determined. +- **[MODEL_INVALID]** The given CpModelProto didn't pass the validation step. + You can get a detailed error by calling `ValidateCpModel(model_proto)`in + C++, or `model.Validate()` in other languages. +- **[FEASIBLE]** A feasible solution has been found. But the search was + stopped before we could prove optimality or before we enumerated all + solutions of a feasibility problem (if asked). +- **[INFEASIBLE]** The problem has been proven infeasible. +- **[OPTIMAL]** An optimal feasible solution has been found. More generally, + this status represent a success. So we also return OPTIMAL if we find a + solution for a pure feasibility problem or if a gap limit has been specified + and we return a solution within this limit. In the case where we need to + return all the feasible solution, this status will only be returned if + solve() enumerated all of them; If we stopped before, we will return + FEASIBLE. + +## Debugging infeasible models + +Oftentimes, solving a model yields an infeasibility status. While this can be a +valid answer, it can happen because of either a data issue, or a model issue. + +Here a simple common sense processes to help diagnose the source of the +infeasibility. + +- reduce the model size if this is a parameter +- remove all constraints while keeping the problem infeasible +- play with the domain of the variables to enlarge them +- inject a known feasible solution and try to find where it breaks +- add assumptions to check the soundness of your data (capacity is >= 0, no + crazy ranges of values, ...) + +## Using assumptions to explain infeasibility + +CP-SAT implements a system of assumptions to find the root cause of an +infeasible problem. To use it, one must instrument constraints with literals and +add these literals to the set of assumptions. + +Note that this set is minimized but not guaranteed to be minimal. To find a +minimal unsatisfiable set (MUS), you must minimize the (weighted) sum of these +assumption literals instead of using the assumptions mechanism. + +Finally, be aware that solving with assumptions is not compatible with +parallelism. Therefore, the number of workers must be set to 1. + +### Python code sample + +```python +#!/usr/bin/env python3 +"""Code sample that solves a model and gets the infeasibility assumptions.""" +from ortools.sat.python import cp_model + + +def main(): + """Showcases assumptions.""" + # Creates the model. + model = cp_model.CpModel() + + # Creates the variables. + x = model.NewIntVar(0, 10, "x") + y = model.NewIntVar(0, 10, "y") + z = model.NewIntVar(0, 10, "z") + a = model.NewBoolVar("a") + b = model.NewBoolVar("b") + c = model.NewBoolVar("c") + + # Creates the constraints. + model.Add(x > y).OnlyEnforceIf(a) + model.Add(y > z).OnlyEnforceIf(b) + model.Add(z > x).OnlyEnforceIf(c) + + # Add assumptions + model.AddAssumptions([a, b, c]) + + # Creates a solver and solves. + solver = cp_model.CpSolver() + status = solver.Solve(model) + + # Print solution. + print(f"Status = {solver.StatusName(status)}") + if status == cp_model.INFEASIBLE: + print( + "SufficientAssumptionsForInfeasibility = " + f"{solver.SufficientAssumptionsForInfeasibility()}" + ) + + +if __name__ == "__main__": + main() +``` + +### C++ code samples + +```cpp +#include + +#include "absl/types/span.h" +#include "ortools/base/logging.h" +#include "ortools/sat/cp_model.h" +#include "ortools/sat/cp_model.pb.h" +#include "ortools/sat/cp_model_solver.h" +#include "ortools/util/sorted_interval_list.h" + +namespace operations_research { +namespace sat { + +void AssumptionsSampleSat() { + CpModelBuilder cp_model; + + const Domain domain(0, 10); + const IntVar x = cp_model.NewIntVar(domain).WithName("x"); + const IntVar y = cp_model.NewIntVar(domain).WithName("y"); + const IntVar z = cp_model.NewIntVar(domain).WithName("z"); + const BoolVar a = cp_model.NewBoolVar().WithName("a"); + const BoolVar b = cp_model.NewBoolVar().WithName("b"); + const BoolVar c = cp_model.NewBoolVar().WithName("c"); + + cp_model.AddGreaterThan(x, y).OnlyEnforceIf(a); + cp_model.AddGreaterThan(y, z).OnlyEnforceIf(b); + cp_model.AddGreaterThan(z, x).OnlyEnforceIf(c); + + // Add assumptions + cp_model.AddAssumptions({a, b, c}); + + // Solving part. + const CpSolverResponse response = Solve(cp_model.Build()); + + // Print solution. + LOG(INFO) << CpSolverResponseStats(response); + if (response.status() == CpSolverStatus::INFEASIBLE) { + for (const int index : + response.sufficient_assumptions_for_infeasibility()) { + LOG(INFO) << index; + } + } +} +} // namespace sat +} // namespace operations_research + +int main(int argc, char** argv) { + operations_research::sat::AssumptionsSampleSat(); + return EXIT_SUCCESS; +} +``` + +### Java code samples + +```java +package com.google.ortools.sat.samples; +import com.google.ortools.Loader; +import com.google.ortools.sat.CpModel; +import com.google.ortools.sat.CpSolver; +import com.google.ortools.sat.CpSolverStatus; +import com.google.ortools.sat.IntVar; +import com.google.ortools.sat.Literal; + +/** Minimal CP-SAT example to showcase assumptions. */ +public class AssumptionsSampleSat { + public static void main(String[] args) { + Loader.loadNativeLibraries(); + // Create the model. + CpModel model = new CpModel(); + + // Create the variables. + IntVar x = model.newIntVar(0, 10, "x"); + IntVar y = model.newIntVar(0, 10, "y"); + IntVar z = model.newIntVar(0, 10, "z"); + Literal a = model.newBoolVar("a"); + Literal b = model.newBoolVar("b"); + Literal c = model.newBoolVar("c"); + + // Creates the constraints. + model.addGreaterThan(x, y).onlyEnforceIf(a); + model.addGreaterThan(y, z).onlyEnforceIf(b); + model.addGreaterThan(z, x).onlyEnforceIf(c); + + // Add assumptions + model.addAssumptions(new Literal[] {a, b, c}); + + // Create a solver and solve the model. + CpSolver solver = new CpSolver(); + CpSolverStatus status = solver.solve(model); + + // Print solution. + // Check that the problem is infeasible. + if (status == CpSolverStatus.INFEASIBLE) { + System.out.println(solver.sufficientAssumptionsForInfeasibility()); + } + } + + private AssumptionsSampleSat() {} +} +``` + +### C\# code samples + +```cs +using System; +using Google.OrTools.Sat; + +public class AssumptionsSampleSat +{ + static void Main() + { + // Creates the model. + CpModel model = new CpModel(); + + // Creates the variables. + IntVar x = model.NewIntVar(0, 10, "x"); + IntVar y = model.NewIntVar(0, 10, "y"); + IntVar z = model.NewIntVar(0, 10, "z"); + ILiteral a = model.NewBoolVar("a"); + ILiteral b = model.NewBoolVar("b"); + ILiteral c = model.NewBoolVar("c"); + + // Creates the constraints. + model.Add(x > y).OnlyEnforceIf(a); + model.Add(y > z).OnlyEnforceIf(b); + model.Add(z > x).OnlyEnforceIf(c); + + // Add assumptions + model.AddAssumptions(new ILiteral[] { a, b, c }); + + // Creates a solver and solves the model. + CpSolver solver = new CpSolver(); + CpSolverStatus status = solver.Solve(model); + Console.WriteLine(solver.SufficientAssumptionsForInfeasibility()); + } +} +``` diff --git a/ortools/sat/feasibility_jump.cc b/ortools/sat/feasibility_jump.cc index 08bad19de5..b32ec5074e 100644 --- a/ortools/sat/feasibility_jump.cc +++ b/ortools/sat/feasibility_jump.cc @@ -13,6 +13,8 @@ #include "ortools/sat/feasibility_jump.h" +#include + #include #include #include @@ -25,6 +27,7 @@ #include #include +#include "absl/functional/any_invocable.h" #include "absl/functional/function_ref.h" #include "absl/log/check.h" #include "absl/random/bit_gen_ref.h" @@ -52,21 +55,54 @@ namespace operations_research::sat { namespace { // How much do we discount moves we might fix later. constexpr double kCompoundDiscount = 1. / 1024; - -std::pair FindBestValue(Domain domain, int64_t current_value, - absl::FunctionRef f) { - std::pair result = std::make_pair(current_value, 0.0); - domain = domain.IntersectionWith( - Domain(current_value, current_value).Complement()); - for (int i = 0; i < domain.NumIntervals(); ++i) { - const auto& [min, max] = domain[i]; - const auto& [val, score] = RangeConvexMinimum(result, min, max + 1, f); - if (score < result.second) result = std::make_pair(val, score); - } - return result; -} } // namespace +JumpTable::JumpTable( + absl::AnyInvocable(int)> compute_jump) + : compute_jump_(std::move(compute_jump)) {} + +void JumpTable::RecomputeAll(int num_variables) { + deltas_.resize(num_variables); + scores_.resize(num_variables); + needs_recomputation_.assign(num_variables, true); +} + +void JumpTable::SetJump(int var, int64_t delta, double score) { + deltas_[var] = delta; + scores_[var] = score; + needs_recomputation_[var] = false; +} + +void JumpTable::Recompute(int var) { needs_recomputation_[var] = true; } + +bool JumpTable::PossiblyGood(int var) const { + return needs_recomputation_[var] || scores_[var] < 0; +} + +bool JumpTable::JumpIsUpToDate(int var) { + const auto& [delta, score] = compute_jump_(var); + if (delta != deltas_[var]) { + LOG(ERROR) << "Incorrect delta for var " << var << ": " << deltas_[var] + << " (should be " << delta << ")"; + } + bool score_ok = true; + if (abs(score - scores_[var]) / std::max(abs(score), 1.0) > 1e-2) { + score_ok = false; + LOG(ERROR) << "Incorrect score for var " << var << ": " << scores_[var] + << " (should be " << score << ") " + << " delta = " << delta; + } + return delta == deltas_[var] && score_ok; +} + +std::pair JumpTable::GetJump(int var) { + if (needs_recomputation_[var]) { + needs_recomputation_[var] = false; + std::tie(deltas_[var], scores_[var]) = compute_jump_(var); + } + return std::make_pair(deltas_[var], scores_[var]); +} + FeasibilityJumpSolver::~FeasibilityJumpSolver() { if (!VLOG_IS_ON(1)) return; std::vector> stats; @@ -101,8 +137,16 @@ void FeasibilityJumpSolver::Initialize() { ReadDomainFromProto(linear_model_->model_proto().variables(v)); var_has_two_values_[v] = var_domains_[v].HasTwoValues(); } + vars_to_scan_.reserve(num_variables); + in_vars_to_scan_.assign(num_variables, false); move_ = std::make_unique(evaluator_.get(), num_variables); + var_occurs_in_non_linear_constraint_.resize(num_variables); + for (int c = 0; c < evaluator_->NumNonLinearConstraints(); ++c) { + for (int v : evaluator_->GeneralConstraintToVars(c)) { + var_occurs_in_non_linear_constraint_[v] = true; + } + } } namespace { @@ -272,7 +316,7 @@ std::string FeasibilityJumpSolver::OneLineStats() const { const std::string non_solution_str = num_infeasible_cts == 0 ? "" - : absl::StrCat(" #good_lin_moves:", FormatCounter(good_jumps_.size()), + : absl::StrCat(" #good_moves:", FormatCounter(vars_to_scan_.size()), " #inf_cts:", FormatCounter(evaluator_->NumInfeasibleConstraints())); @@ -409,6 +453,7 @@ std::function FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) { if (params_.violation_ls_use_compound_moves()) { use_compound_moves_ = absl::Bernoulli(random_, 0.25); if (use_compound_moves_) { + compound_move_max_discrepancy_ = 0; compound_weight_changed_.clear(); in_compound_weight_changed_.assign(weights_.size(), false); compound_weights_.assign(weights_.size(), kCompoundDiscount); @@ -450,124 +495,140 @@ std::function FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) { }; } -bool FeasibilityJumpSolver::IsGood(int var) const { - if (jump_scores_[var] < 0.0) return true; - if (jump_scores_[var] > 0.0) return false; - return evaluator_->ObjectiveDelta(var, jump_deltas_[var]) < 0; +double FeasibilityJumpSolver::ComputeScore( + absl::Span scan_weights, int var, int64_t delta, + bool linear_only) const { + constexpr double kEpsilon = 1.0 / std::numeric_limits::max(); + double score = evaluator_->LinearEvaluator().WeightedViolationDelta( + scan_weights, var, delta); + if (!linear_only) { + score += + evaluator_->WeightedNonLinearViolationDelta(scan_weights, var, delta); + } + score += kEpsilon * evaluator_->ObjectiveDelta(var, delta); + return score; } -void FeasibilityJumpSolver::RecomputeJump(int var) { +std::pair FeasibilityJumpSolver::ComputeLinearJump(int var) { const std::vector& solution = evaluator_->current_solution(); - ++num_linear_evals_; - jump_need_recomputation_[var] = false; if (var_domains_[var].IsFixed()) { - jump_deltas_[var] = 0; - jump_scores_[var] = 0.0; - return; + return std::make_pair(0l, 0.0); } + + ++num_linear_evals_; const LinearIncrementalEvaluator& linear_evaluator = evaluator_->LinearEvaluator(); if (var_has_two_values_[var]) { const int64_t min_value = var_domains_[var].Min(); const int64_t max_value = var_domains_[var].Max(); - jump_deltas_[var] = solution[var] == min_value ? max_value - min_value - : min_value - max_value; - jump_scores_[var] = linear_evaluator.WeightedViolationDelta( - weights_, var, jump_deltas_[var]); - } else { - // In practice, after a few iterations, the chance of finding an improving - // move is slim, and we can test that fairly easily with at most two - // queries! - // - // Tricky/Annoying: if the value is not in the domain, we returns it. - const int64_t p1 = var_domains_[var].ValueAtOrBefore(solution[var] - 1); - const int64_t p2 = var_domains_[var].ValueAtOrAfter(solution[var] + 1); + const int64_t delta = solution[var] == min_value ? max_value - min_value + : min_value - max_value; + return std::make_pair(delta, + ComputeScore(var, delta, /*linear_only=*/true)); + } - std::pair best_jump; - const double v1 = var_domains_[var].Contains(p1) - ? linear_evaluator.WeightedViolationDelta( - weights_, var, p1 - solution[var]) - : std::numeric_limits::infinity(); - if (v1 < 0.0) { - // Point p1 is improving. Look for best before it. - // Note that we can exclude all point after solution[var] since it is - // worse and we assume convexity. + // In practice, after a few iterations, the chance of finding an improving + // move is slim, and we can test that fairly easily with at most two + // queries! + // + // Tricky/Annoying: if the value is not in the domain, we returns it. + const int64_t p1 = var_domains_[var].ValueAtOrBefore(solution[var] - 1); + const int64_t p2 = var_domains_[var].ValueAtOrAfter(solution[var] + 1); + + std::pair best_jump; + const double v1 = + var_domains_[var].Contains(p1) + ? ComputeScore(var, p1 - solution[var], /*linear_only=*/true) + : std::numeric_limits::infinity(); + if (v1 < 0.0) { + // Point p1 is improving. Look for best before it. + // Note that we can exclude all point after solution[var] since it is + // worse and we assume convexity. + const Domain dom = var_domains_[var].IntersectionWith( + Domain(std::numeric_limits::min(), p1 - 1)); + if (dom.IsEmpty()) { + best_jump = {p1, v1}; + } else { + tmp_breakpoints_ = + linear_evaluator.SlopeBreakpoints(var, solution[var], dom); + best_jump = ConvexMinimum( + /*is_to_the_right=*/true, {p1, v1}, tmp_breakpoints_, + [this, var, &solution](int64_t jump_value) { + return ComputeScore(var, jump_value - solution[var], + /*linear_only=*/true); + }); + } + } else { + const double v2 = + var_domains_[var].Contains(p2) + ? ComputeScore(var, p2 - solution[var], /*linear_only=*/true) + : std::numeric_limits::infinity(); + if (v2 < 0.0) { + // Point p2 is improving. Look for best after it. + // Similarly, we exclude the other points by convexity. const Domain dom = var_domains_[var].IntersectionWith( - Domain(std::numeric_limits::min(), p1 - 1)); + Domain(p2 + 1, std::numeric_limits::max())); if (dom.IsEmpty()) { - best_jump = {p1, v1}; + best_jump = {p2, v2}; } else { tmp_breakpoints_ = linear_evaluator.SlopeBreakpoints(var, solution[var], dom); best_jump = ConvexMinimum( - /*is_to_the_right=*/true, {p1, v1}, tmp_breakpoints_, - [var, this, &linear_evaluator, &solution](int64_t jump_value) { - return linear_evaluator.WeightedViolationDelta( - weights_, var, jump_value - solution[var]); + /*is_to_the_right=*/false, {p2, v2}, tmp_breakpoints_, + [this, var, &solution](int64_t jump_value) { + return ComputeScore(var, jump_value - solution[var], + /*linear_only=*/true); }); } } else { - const double v2 = var_domains_[var].Contains(p2) - ? linear_evaluator.WeightedViolationDelta( - weights_, var, p2 - solution[var]) - : std::numeric_limits::infinity(); - if (v2 < 0.0) { - // Point p2 is improving. Look for best after it. - // Similarly, we exclude the other points by convexity. - const Domain dom = var_domains_[var].IntersectionWith( - Domain(p2 + 1, std::numeric_limits::max())); - if (dom.IsEmpty()) { - best_jump = {p2, v2}; - } else { - tmp_breakpoints_ = - linear_evaluator.SlopeBreakpoints(var, solution[var], dom); - best_jump = ConvexMinimum( - /*is_to_the_right=*/false, {p2, v2}, tmp_breakpoints_, - [var, this, &linear_evaluator, &solution](int64_t jump_value) { - return linear_evaluator.WeightedViolationDelta( - weights_, var, jump_value - solution[var]); - }); - } + // We have no improving point, result is either p1 or p2. This is the + // most common scenario, and require no breakpoint computation! + // Choose the direction which increases violation the least, + // disambiguating by best objective. + if (v1 < v2) { + best_jump = {p1, v1}; } else { - // We have no improving point, result is either p1 or p2. This is the - // most common scenario, and require no breakpoint computation! - // Choose the direction which increases violation the least, - // disambiguating by best objective. - if (v1 < v2 || (v1 == v2 && evaluator_->ObjectiveDelta( - var, p1 - solution[var]) < 0)) { - best_jump = {p1, v1}; - } else { - best_jump = {p2, v2}; - } + best_jump = {p2, v2}; } } - - DCHECK_NE(best_jump.first, solution[var]); - jump_deltas_[var] = best_jump.first - solution[var]; - jump_scores_[var] = best_jump.second; } + DCHECK_NE(best_jump.first, solution[var]); + return std::make_pair(best_jump.first - solution[var], best_jump.second); } -void FeasibilityJumpSolver::RecomputeAllJumps() { - const int num_variables = static_cast(var_domains_.size()); - jump_deltas_.resize(num_variables); - jump_scores_.resize(num_variables); - jump_need_recomputation_.assign(num_variables, true); - - in_good_jumps_.assign(num_variables, false); - good_jumps_.clear(); - - for (int var = 0; var < num_variables; ++var) { - RecomputeJump(var); - if (IsGood(var)) { - in_good_jumps_[var] = true; - good_jumps_.push_back(var); - } +std::pair FeasibilityJumpSolver::ComputeGeneralJump(int var) { + if (!var_occurs_in_non_linear_constraint_[var]) { + return ComputeLinearJump(var); } + Domain domain = var_domains_[var]; + if (domain.IsFixed()) return std::make_pair(0, 0.0); + + ++num_general_evals_; + const int64_t current_value = evaluator_->current_solution()[var]; + domain = domain.IntersectionWith( + Domain(current_value, current_value).Complement()); + std::pair result = RangeConvexMinimum( + domain[0].start - current_value, domain[0].end - current_value + 1, + [&](int64_t delta) -> double { + return ComputeScore(var, delta, /*linear_only=*/false); + }); + for (int i = 1; i < domain.NumIntervals(); ++i) { + const int64_t min_delta = domain[i].start - current_value; + const int64_t max_delta = domain[i].end - current_value; + const auto& [delta, score] = RangeConvexMinimum( + result, min_delta, max_delta + 1, [&](int64_t delta) -> double { + return ComputeScore(var, delta, /*linear_only=*/false); + }); + if (score < result.second) result = std::make_pair(delta, score); + } + DCHECK(domain.Contains(current_value + result.first)) + << current_value << "+" << result.first << " not in domain " + << domain.ToString(); + return result; } -void FeasibilityJumpSolver::UpdateViolatedConstraintWeights() { +void FeasibilityJumpSolver::UpdateViolatedConstraintWeights(JumpTable& jumps) { ++num_weight_updates_; // Because we update the weight incrementally, it is better to not have a @@ -576,6 +637,7 @@ void FeasibilityJumpSolver::UpdateViolatedConstraintWeights() { // DCHECK(JumpIsUpToDate(var)) will fail more often. const double kMaxWeight = 1e10; const double kBumpFactor = 1.0 / params_.feasibility_jump_decay(); + const int num_variables = var_domains_.size(); if (use_decay_) { bump_value_ *= kBumpFactor; } @@ -597,7 +659,7 @@ void FeasibilityJumpSolver::UpdateViolatedConstraintWeights() { weights_[c] *= factor; if (use_compound_moves_) compound_weights_[c] *= factor; } - RecomputeAllJumps(); + jumps.RecomputeAll(num_variables); return; } @@ -610,10 +672,18 @@ void FeasibilityJumpSolver::UpdateViolatedConstraintWeights() { LinearIncrementalEvaluator* linear_evaluator = evaluator_->MutableLinearEvaluator(); linear_evaluator->ClearAffectedVariables(); - for_weight_update_.resize(jump_scores_.size()); + for_weight_update_.resize(num_variables); for (const int c : evaluator_->ViolatedConstraints()) { - linear_evaluator->UpdateScoreOnWeightUpdate( - c, jump_deltas_, absl::MakeSpan(for_weight_update_)); + // TODO(user): We can probably handle compound moves too. + if (c < evaluator_->NumLinearConstraints() && !use_compound_moves_) { + linear_evaluator->UpdateScoreOnWeightUpdate( + c, jumps.Deltas(), absl::MakeSpan(for_weight_update_)); + } else { + for (const int v : evaluator_->ConstraintToVars(c)) { + jumps.Recompute(v); + AddVarToScan(jumps, v); + } + } } // Recompute the affected jumps. @@ -621,65 +691,32 @@ void FeasibilityJumpSolver::UpdateViolatedConstraintWeights() { for (const int var : linear_evaluator->VariablesAffectedByLastUpdate()) { // Apply the delta. // - // TODO(user): We could compute the minimal bump that lead to a good move. - // That might change depending on the jump value though, so we can only - // do that easily for Boolean I think. - jump_scores_[var] += bump_value_ * for_weight_update_[var]; - - // We don't need to recompute score of binary variable, it should - // already be correct. - if (!jump_need_recomputation_[var] && var_has_two_values_[var]) { - DCHECK(JumpIsUpToDate(var)); - if (IsGood(var) && !in_good_jumps_[var]) { - in_good_jumps_[var] = true; - good_jumps_.push_back(var); - } - continue; - } - - // This jump might be good, so we need to add it to the queue so it can be - // evaluated when choosing the next jump. - jump_need_recomputation_[var] = true; - if (!in_good_jumps_[var]) { - in_good_jumps_[var] = true; - good_jumps_.push_back(var); + // TODO(user): We could compute the minimal bump that would lead to a + // good move. That might change depending on the jump value though, so + // we can only do that easily for Booleans. + // TODO(user): We can probably handle compound moves too. + if (!var_has_two_values_[var] || use_compound_moves_) { + jumps.Recompute(var); + } else { + // We may know the correct score for binary vars. + jumps.MutableScores()[var] += bump_value_ * for_weight_update_[var]; } + AddVarToScan(jumps, var); } } -// Important: This is for debugging, but unfortunately it currently change the -// deterministic time and so the overall algo behavior. -bool FeasibilityJumpSolver::JumpIsUpToDate(int var) { - // With decay, we have a wide magnitude of weight, so the incremental update - // of the floating point score can be widely off. Here we just check that - // without decay (integer weight) the update should be reasonably precise. - if (use_decay_) return true; - - const int64_t old_jump_delta = jump_deltas_[var]; - const double old_score = jump_scores_[var]; - RecomputeJump(var); - CHECK_EQ(jump_deltas_[var], old_jump_delta); // No change - const double relative = - std::max({std::abs(jump_scores_[var]), std::abs(old_score), 1.0}); - const bool result = std::abs(jump_scores_[var] - old_score) / relative < 1e-2; - - // To keep the same behavior as in -c opt, we restore the old value. - jump_scores_[var] = old_score; - return result; -} - bool FeasibilityJumpSolver::DoSomeLinearIterations() { - RecomputeAllJumps(); - evaluator_->RecomputeViolatedList(/*linear_only=*/true); - if (VLOG_IS_ON(1)) { shared_response_->LogMessageWithThrottling(name(), OneLineStats()); } - // TODO(user): It should be possible to support compound moves with the - // specialized linear code, but lets keep it simpler for now. + // TODO(user): It should be possible to support compound moves with + // the specialized linear code, but lets keep it simpler for now. if (use_compound_moves_) return true; + evaluator_->RecomputeViolatedList(/*linear_only=*/true); + RecomputeVarsToScan(linear_jumps_); + // Do a batch of a given number of loop here. // Outer loop: when no more greedy moves, update the weight. const int kBatchSize = 10000; @@ -687,90 +724,41 @@ bool FeasibilityJumpSolver::DoSomeLinearIterations() { for (int loop = 0; loop < kBatchSize; ++loop) { // Inner loop: greedy descent. for (; loop < kBatchSize; ++loop) { - // Test the shared limit not too often. - // - // TODO(user): depending on the size of the problem that might be too - // little, use deterministic time instead. - if (loop % 100 == 0) { - if (shared_time_limit_->LimitReached()) { - return false; - } - } - // Take the best jump score amongst some random candidates. // It is okay if we pick twice the same, we don't really care. - int best_var = -1; - int best_index = -1; - int64_t best_delta = 0; - double best_score = 0.0; - int64_t best_obj_delta = 0; - int num_improving_jump_tested = 0; - while (!good_jumps_.empty() && num_improving_jump_tested < 5) { - const int index = absl::Uniform( - random_, 0, static_cast(good_jumps_.size())); - const int var = good_jumps_[index]; - - // We lazily update the jump value. - if (jump_need_recomputation_[var]) { - RecomputeJump(var); - } else { - DCHECK(JumpIsUpToDate(var)); - } - - if (!IsGood(var)) { - // Lazily remove. - in_good_jumps_[var] = false; - good_jumps_[index] = good_jumps_.back(); - good_jumps_.pop_back(); - if (best_index == good_jumps_.size()) best_index = index; - continue; - } - - ++num_improving_jump_tested; - const int64_t obj_delta = - evaluator_->ObjectiveDelta(var, jump_deltas_[var]); - if (std::make_pair(jump_scores_[var], obj_delta) < - std::make_pair(best_score, best_obj_delta)) { - best_var = var; - best_index = index; - best_delta = jump_deltas_[var]; - best_score = jump_scores_[var]; - best_obj_delta = obj_delta; - } + int best_var; + int64_t best_value; + double best_score; + if (!ScanRelevantVariables(/*num_to_scan=*/5, linear_jumps_, &best_var, + &best_value, &best_score)) { + break; } - - if (good_jumps_.empty()) break; - DCHECK_NE(best_var, -1); - DCHECK_NE(best_index, -1); + const int64_t current_value = solution[best_var]; // Perform the move. ++num_linear_moves_; - const int64_t best_value = solution[best_var] + best_delta; evaluator_->UpdateLinearScores(best_var, best_value, weights_, - jump_deltas_, - absl::MakeSpan(jump_scores_)); + linear_jumps_.Deltas(), + linear_jumps_.MutableScores()); evaluator_->UpdateVariableValue(best_var, best_value); - // We already know the score of undoing the move we just did, and we know - // this move is bad, so we can remove it from good_jumps_ right away. - jump_deltas_[best_var] = -jump_deltas_[best_var]; - jump_scores_[best_var] = -best_score; if (var_has_two_values_[best_var]) { - CHECK_EQ(good_jumps_[best_index], best_var); - in_good_jumps_[best_var] = false; - good_jumps_[best_index] = good_jumps_.back(); - good_jumps_.pop_back(); + // We already know the score of undoing the move we just did, and that + // this is optimal. + linear_jumps_.SetJump(best_var, current_value - best_value, + -best_score); } else { - jump_need_recomputation_[best_var] = true; + linear_jumps_.Recompute(best_var); } - MarkJumpsThatNeedsToBeRecomputed(best_var); + MarkJumpsThatNeedToBeRecomputed(best_var, linear_jumps_); } + if (time_limit_crossed_) return false; // We will update the weight unless the queue is non-empty. - if (good_jumps_.empty()) { + if (vars_to_scan_.empty()) { // Note that we only count linear constraint as violated here. if (evaluator_->ViolatedConstraints().empty()) return true; - UpdateViolatedConstraintWeights(); + UpdateViolatedConstraintWeights(linear_jumps_); } } return false; @@ -789,29 +777,22 @@ bool FeasibilityJumpSolver::DoSomeLinearIterations() { // TODO(user): For non-Boolean, we could easily detect if a non-improving // score cannot become improving. We don't need to add such variable to // the queue. -void FeasibilityJumpSolver::MarkJumpsThatNeedsToBeRecomputed(int changed_var) { +void FeasibilityJumpSolver::MarkJumpsThatNeedToBeRecomputed(int changed_var, + JumpTable& jumps) { for (const int var : evaluator_->VariablesAffectedByLastLinearUpdate()) { - if (var == changed_var) continue; - if (jump_need_recomputation_[var]) { - DCHECK(in_good_jumps_[var]); - continue; + // TODO(user): We can probably handle compound moves too. + if (var != changed_var && + (!var_has_two_values_[var] || use_compound_moves_)) { + jumps.Recompute(var); } - - // We don't need to recompute score of binary variable, it should - // already be correct. - if (var_has_two_values_[var]) { - DCHECK(JumpIsUpToDate(var)) << var << " " << jump_scores_[var]; - if (IsGood(var) && !in_good_jumps_[var]) { - in_good_jumps_[var] = true; - good_jumps_.push_back(var); - } - continue; - } - - jump_need_recomputation_[var] = true; - if (!in_good_jumps_[var]) { - in_good_jumps_[var] = true; - good_jumps_.push_back(var); + AddVarToScan(jumps, var); + } + for (const auto& [c, violation_delta] : + evaluator_->last_update_violation_changes()) { + if (c < evaluator_->NumLinearConstraints()) continue; + for (const int var : evaluator_->ConstraintToVars(c)) { + if (var != changed_var) jumps.Recompute(var); + AddVarToScan(jumps, var); } } } @@ -821,105 +802,93 @@ bool FeasibilityJumpSolver::DoSomeGeneralIterations() { return true; } const std::vector& solution = evaluator_->current_solution(); - // Non-linear constraints are not evaluated in the linear phase. evaluator_->UpdateAllNonLinearViolations(); evaluator_->RecomputeViolatedList(/*linear_only=*/false); - RecomputeVarsToScan(); - if (use_compound_moves_) { - compound_move_beam_search_width_ = absl::Bernoulli(random_, 0.75) ? 1 : 2; - } + RecomputeVarsToScan(general_jumps_); auto effort = [&]() { - return num_general_evals_ + num_weight_updates_ + num_general_moves_; + return num_general_evals_ + num_linear_evals_ + num_weight_updates_ + + num_general_moves_; }; - const int64_t effort_limit = effort() + 100000; + const int64_t effort_limit = effort() + 20000; while (effort() < effort_limit) { int var; int64_t value; double score; - bool time_limit_crossed = false; - DCHECK(!move_->IsImproving()); - const bool found_move = - ScanRelevantVariables(&var, &value, &score, &time_limit_crossed); + const bool found_move = ScanRelevantVariables( + /*num_to_scan=*/3, general_jumps_, &var, &value, &score); const bool backtrack = !found_move && move_->Backtrack(&var, &value, &score); if (found_move || backtrack) { - const int64_t prev_value = solution[var]; // Perform the move. - num_general_moves_++; - + ++num_general_moves_; + CHECK_NE(var, -1) << var << " " << found_move << " " << backtrack; + const int64_t prev_value = solution[var]; + DCHECK_NE(prev_value, value); // Update the linear part. - evaluator_->UpdateLinearScores(var, value, weights_, jump_deltas_, - absl::MakeSpan(jump_scores_)); - jump_scores_[var] = -score; - jump_deltas_[var] = -jump_deltas_[var]; - - // This score might include non-linear weights, so may be wrong. - // We don't actually use it, but if we make things more incremental - // across batches we may want this in future and it stops - // UpdateViolatedConstraintWeights() from DCHECKing. - jump_need_recomputation_[var] = MustRecomputeJumpOnGeneralUpdate(var); - + evaluator_->UpdateLinearScores(var, value, weights_, + general_jumps_.Deltas(), + general_jumps_.MutableScores()); // Update the non-linear part. Note it also commits the move. evaluator_->UpdateNonLinearViolations(var, value); evaluator_->UpdateVariableValue(var, value); - for (const auto& [c, violation_delta] : - evaluator_->last_update_violation_changes()) { - if (violation_delta == 0) continue; - for (const int var : evaluator_->ConstraintToVars(c)) { - if (violation_delta == evaluator_->Violation(c)) { - num_violated_constraints_per_var_[var] += 1; - } else if (violation_delta < 0 && !evaluator_->IsViolated(c)) { - num_violated_constraints_per_var_[var] -= 1; - } - if (use_compound_moves_ && !in_compound_weight_changed_[c]) { + // TODO(user): Handle compound moves? + if (use_compound_moves_ && !backtrack) { + // `!backtrack` is just an optimisation - we can never break any new + // constraints on backtrack, so we can never change any + // compound_weight_. + for (const auto& [c, violation_delta] : + evaluator_->last_update_violation_changes()) { + if (violation_delta == 0) continue; + if (violation_delta > 0 && compound_weights_[c] != weights_[c]) { compound_weights_[c] = weights_[c]; - compound_weight_changed_.push_back(c); + if (!in_compound_weight_changed_[c]) { + in_compound_weight_changed_[c] = true; + compound_weight_changed_.push_back(c); + } + for (const int v : evaluator_->ConstraintToVars(c)) { + general_jumps_.Recompute(v); + // Vars will be added in MarkJumpsThatNeedToBeRecomputed. + } + } else if (violation_delta < 0 && !in_compound_weight_changed_[c] && + !evaluator_->IsViolated(c)) { + DCHECK_EQ(compound_weights_[c], weights_[c]); in_compound_weight_changed_[c] = true; + compound_weight_changed_.push_back(c); } } } - - // We call AddVarToScan() after num_violated_constraints_per_var_ has - // been computed. - for (const int v : evaluator_->VariablesAffectedByLastLinearUpdate()) { - jump_need_recomputation_[v] = MustRecomputeJumpOnGeneralUpdate(v); - AddVarToScan(v); + if (!use_decay_) { + DCHECK_EQ(-score, ComputeScore(var, prev_value - value, false)); } - for (const int general_c : evaluator_->VarToGeneralConstraints(var)) { - for (const int v : evaluator_->GeneralConstraintToVars(general_c)) { - AddVarToScan(v); - } + if (var_has_two_values_[var]) { + // We already know the score of the only possible move (undoing what we + // just did). + general_jumps_.SetJump(var, prev_value - value, -score); + } else { + general_jumps_.Recompute(var); } - DCHECK(SlowCheckNumViolatedConstraints()); - + MarkJumpsThatNeedToBeRecomputed(var, general_jumps_); if (use_compound_moves_ && !backtrack) { // Make sure we can undo the move. move_->Push(var, prev_value, score); - if (move_->IsImproving()) { - if (move_->Size() > 1) { - num_compound_moves_ += move_->Size(); - } + if (move_->Score() < 0) { + num_compound_moves_ += move_->Size(); move_->Clear(); - ResetChangedCompoundWeights(); + compound_move_max_discrepancy_ = 0; } } continue; - } else if (time_limit_crossed) { + } else if (time_limit_crossed_) { return false; } DCHECK_EQ(move_->Size(), 0); if (evaluator_->ViolatedConstraints().empty()) return true; - UpdateViolatedConstraintWeights(); - - DCHECK(SlowCheckNumViolatedConstraints()); - // Constraints with increased weight may lead to new negative score moves. - for (const int c : evaluator_->ViolatedConstraints()) { - for (const int v : evaluator_->ConstraintToVars(c)) { - AddVarToScan(v); - } + if (use_compound_moves_) ResetChangedCompoundWeights(); + if (!use_compound_moves_ || ++compound_move_max_discrepancy_ > 2) { + compound_move_max_discrepancy_ = 0; + UpdateViolatedConstraintWeights(general_jumps_); } - ResetChangedCompoundWeights(); } return false; } @@ -929,136 +898,135 @@ void FeasibilityJumpSolver::ResetChangedCompoundWeights() { DCHECK_EQ(move_->Size(), 0); for (const int c : compound_weight_changed_) { in_compound_weight_changed_[c] = false; - compound_weights_[c] = weights_[c]; - if (!evaluator_->IsViolated(c)) { - compound_weights_[c] *= kCompoundDiscount; - for (const int var : evaluator_->ConstraintToVars(c)) { - AddVarToScan(var); - } + double expected_weight = + (evaluator_->IsViolated(c) ? 1.0 : kCompoundDiscount) * weights_[c]; + if (compound_weights_[c] == expected_weight) continue; + compound_weights_[c] = expected_weight; + for (const int var : evaluator_->ConstraintToVars(c)) { + general_jumps_.Recompute(var); + AddVarToScan(general_jumps_, var); } } compound_weight_changed_.clear(); } -bool FeasibilityJumpSolver::ShouldAcceptUnitMove(int var, int64_t delta, - double score) { - if (score < 0) return true; - return score == 0 && evaluator_->ObjectiveDelta(var, delta) < 0; -} - -bool FeasibilityJumpSolver::ShouldExtendCompoundMove(int var, int64_t delta, - double score, +bool FeasibilityJumpSolver::ShouldExtendCompoundMove(double score, double novelty) { if (move_->Score() + score - std::max(novelty, 0.0) < 0) { return true; } - if (move_->Score() + score == 0 && - move_->ObjectiveDelta() + evaluator_->ObjectiveDelta(var, delta) < 0) { - return true; - } return score < move_->BestChildScore(); } -bool FeasibilityJumpSolver::ScanRelevantVariables(int* improving_var, - int64_t* improving_value, - double* improving_score, - bool* time_limit_crossed) { - if (move_->NumChildrenExplored() >= compound_move_beam_search_width_) { +bool FeasibilityJumpSolver::ScanRelevantVariables(int num_to_scan, + JumpTable& jumps, + int* best_var, + int64_t* best_value, + double* best_score) { + if (move_->Discrepancy() > compound_move_max_discrepancy_) { return false; } - DCHECK_GE(move_->Score(), 0); - absl::Span solution = evaluator_->current_solution(); - absl::Span scan_weights = - use_compound_moves_ ? compound_weights_ : weights_; - while (!vars_to_scan_.empty()) { - const int idx = absl::Uniform(random_, 0, vars_to_scan_.size()); - const int var = vars_to_scan_[idx]; - vars_to_scan_[idx] = vars_to_scan_.back(); - vars_to_scan_.pop_back(); - in_vars_to_scan_[var] = false; - // Skip evaluating `var` if it cannot have an improving move. - if (!ShouldScan(var)) continue; + double best_scan_score = 0.0; + int num_good = 0; + int best_index = -1; + *best_var = -1; + *best_score = 0.0; - int64_t new_value; - double scan_score; - const int64_t current_value = solution[var]; - if (!use_compound_moves_ && - evaluator_->VariableOnlyInLinearConstraintWithConvexViolationChange( - var)) { - // We lazily update the jump value. - if (jump_need_recomputation_[var]) { - // We don't use good_jumps_ and in_good_jumps_ here, so we don't update - // them. - RecomputeJump(var); - } else { - DCHECK(JumpIsUpToDate(var)); - } - new_value = current_value + jump_deltas_[var]; - scan_score = jump_scores_[var]; - } else { - std::tie(new_value, scan_score) = FindBestValue( - var_domains_[var], current_value, [&](int64_t value) -> double { - // Check the time limit periodically. - if (++num_general_evals_ % 1000 == 0 && - shared_time_limit_ != nullptr && - shared_time_limit_->LimitReached()) { - *time_limit_crossed = true; - } - if (*time_limit_crossed) return 0.0; - return evaluator_->WeightedViolationDelta(scan_weights, var, - value - current_value); - }); + auto remove_var_to_scan_at_index = [&](int index) { + in_vars_to_scan_[vars_to_scan_[index]] = false; + vars_to_scan_[index] = vars_to_scan_.back(); + vars_to_scan_.pop_back(); + if (best_index == vars_to_scan_.size()) { + best_index = index; } - if (*time_limit_crossed) return false; - const int64_t delta = new_value - current_value; - if (delta == 0) continue; - const double score = use_compound_moves_ - ? evaluator_->WeightedViolationDelta( - weights_, var, new_value - current_value) - : scan_score; + }; + while (!vars_to_scan_.empty() && num_good < num_to_scan) { + const int index = absl::Uniform(random_, 0, vars_to_scan_.size()); + const int var = vars_to_scan_[index]; + DCHECK_GE(var, 0); + DCHECK(in_vars_to_scan_[var]); + + if (!ShouldScan(jumps, var)) { + remove_var_to_scan_at_index(index); + continue; + } + const auto [delta, scan_score] = jumps.GetJump(var); + if ((num_general_evals_ + num_linear_evals_) % 100 == 0 && + shared_time_limit_ != nullptr && shared_time_limit_->LimitReached()) { + time_limit_crossed_ = true; + return false; + } + if (time_limit_crossed_) return false; + const int64_t current_value = evaluator_->current_solution()[var]; + DCHECK(var_domains_[var].Contains(current_value + delta)); + DCHECK(!var_domains_[var].IsFixed()); + // Note that this will likely fail if you use decaying weights as they + // will have large magnitudes and the incremental update will be + // imprecise. + DCHECK(use_decay_ || jumps.JumpIsUpToDate(var)) + << var << " " << var_domains_[var].ToString(); + if (scan_score >= 0) { + remove_var_to_scan_at_index(index); + continue; + } + double score = scan_score; if (use_compound_moves_) { - if (!ShouldExtendCompoundMove(var, delta, score, score - scan_score)) { + // We only use compound moves in general iterations. + score = ComputeScore(weights_, var, delta, /*linear_only=*/false); + if (!ShouldExtendCompoundMove(score, score - scan_score)) { + remove_var_to_scan_at_index(index); continue; } - } else { - if (!ShouldAcceptUnitMove(var, delta, score)) continue; } - *improving_var = var; - *improving_value = new_value; - *improving_score = score; + + ++num_good; + if (scan_score < best_scan_score) { + CHECK_NE(delta, 0) << score; + *best_var = var; + *best_value = current_value + delta; + *best_score = score; + best_scan_score = scan_score; + best_index = index; + } + } + if (best_index != -1) { + DCHECK_GT(num_good, 0); + DCHECK_GE(*best_var, 0); + DCHECK_EQ(*best_var, vars_to_scan_[best_index]); + remove_var_to_scan_at_index(best_index); return true; } return false; } -bool FeasibilityJumpSolver::MustRecomputeJumpOnGeneralUpdate(int var) const { - return !evaluator_->VariableOnlyInLinearConstraintWithConvexViolationChange( - var) || - !var_has_two_values_[var]; -} - -void FeasibilityJumpSolver::AddVarToScan(int var) { - if (in_vars_to_scan_[var] || !ShouldScan(var)) return; +void FeasibilityJumpSolver::AddVarToScan(const JumpTable& jumps, int var) { + DCHECK_GE(var, 0); + if (in_vars_to_scan_[var] || !ShouldScan(jumps, var)) return; vars_to_scan_.push_back(var); in_vars_to_scan_[var] = true; } -bool FeasibilityJumpSolver::ShouldScan(int var) const { - if (move_->OnStack(var) || var_domains_[var].IsFixed()) return false; - if (num_violated_constraints_per_var_[var] > 0) return true; +bool FeasibilityJumpSolver::ShouldScan(const JumpTable& jumps, int var) const { + DCHECK_GE(var, 0); + if (var_domains_[var].IsFixed()) return false; + if (!jumps.PossiblyGood(var)) return false; + if (move_->OnStack(var)) return false; + if (evaluator_->NumViolatedConstraintsForVar(var) > 0) return true; const int64_t value = evaluator_->current_solution()[var]; // Return true iff var is has a better objective value in its domain. - return evaluator_->ObjectiveDelta(var, var_domains_[var].Max() - value) < 0 || - evaluator_->ObjectiveDelta(var, var_domains_[var].Min() - value) < 0; + return evaluator_->ObjectiveDelta(var, var_domains_[var].Min() - value) < 0 || + evaluator_->ObjectiveDelta(var, var_domains_[var].Max() - value) < 0; } -void FeasibilityJumpSolver::RecomputeVarsToScan() { - num_violated_constraints_per_var_.assign(var_domains_.size(), 0); - in_vars_to_scan_.assign(evaluator_->current_solution().size(), false); +void FeasibilityJumpSolver::RecomputeVarsToScan(JumpTable& jumps) { + const int num_variables = var_domains_.size(); + jumps.RecomputeAll(num_variables); + in_vars_to_scan_.assign(num_variables, false); + vars_to_scan_.clear(); + DCHECK(SlowCheckNumViolatedConstraints()); for (const int c : evaluator_->ViolatedConstraints()) { for (const int v : evaluator_->ConstraintToVars(c)) { - num_violated_constraints_per_var_[v] += 1; - AddVarToScan(v); + AddVarToScan(jumps, v); } } } @@ -1072,15 +1040,11 @@ bool FeasibilityJumpSolver::SlowCheckNumViolatedConstraints() const { } } for (int v = 0; v < result.size(); ++v) { - CHECK_EQ(result[v], num_violated_constraints_per_var_[v]); + CHECK_EQ(result[v], evaluator_->NumViolatedConstraintsForVar(v)); } return true; } -bool CompoundMoveBuilder::IsImproving() const { - return Score() < 0 || (Score() == 0 && ObjectiveDelta() < 0); -} - void CompoundMoveBuilder::Clear() { for (const UnitMove& move : stack_) { var_on_stack_[move.var] = false; @@ -1101,25 +1065,26 @@ bool CompoundMoveBuilder::Backtrack(int* var, int64_t* value, double* score) { var_on_stack_[*var] = false; stack_.pop_back(); DCHECK_NE(*value, evaluator_->current_solution()[*var]); + if (!stack_.empty()) { + ++stack_.back().discrepancy; + } return true; } void CompoundMoveBuilder::Push(int var, int64_t prev_value, double score) { DCHECK_NE(prev_value, evaluator_->current_solution()[var]); - const double obj_delta = evaluator_->ObjectiveDelta( - var, evaluator_->current_solution()[var] - prev_value); DCHECK(!var_on_stack_[var]); if (!stack_.empty()) { stack_.back().best_child_score = std::min(stack_.back().best_child_score, score); - ++stack_.back().num_children; } var_on_stack_[var] = true; - stack_.push_back( - {.var = var, - .prev_value = prev_value, - .score = -score, - .cumulative_score = Score() + score, - .cumulative_objective_delta = ObjectiveDelta() + obj_delta}); + stack_.push_back({ + .var = var, + .prev_value = prev_value, + .score = -score, + .cumulative_score = Score() + score, + .discrepancy = Discrepancy(), + }); } } // namespace operations_research::sat diff --git a/ortools/sat/feasibility_jump.h b/ortools/sat/feasibility_jump.h index 87e3111310..d618c37b98 100644 --- a/ortools/sat/feasibility_jump.h +++ b/ortools/sat/feasibility_jump.h @@ -20,9 +20,11 @@ #include #include #include +#include #include -#include "absl/time/time.h" +#include "absl/functional/any_invocable.h" +#include "absl/functional/bind_front.h" #include "absl/types/span.h" #include "ortools/sat/constraint_violation.h" #include "ortools/sat/linear_model.h" @@ -36,6 +38,55 @@ namespace operations_research::sat { class CompoundMoveBuilder; +// This class lazily caches the results of `compute_jump(var)` which returns a +// pair. +// Variables' scores can be manually modified using MutableScores (if the +// optimal jump is known not to change), or marked for recomputation on the next +// call to GetJump(var) by calling Recompute. +class JumpTable { + public: + explicit JumpTable( + absl::AnyInvocable(int)> compute_jump); + void RecomputeAll(int num_variables); + + // Gets the current jump delta and score, recomputing if necessary. + std::pair GetJump(int var); + // If the new optimum value and score is known, users can update it directly. + // e.g. after weight rescaling, or after changing a binary variable. + void SetJump(int var, int64_t delta, double score); + // Recompute the jump for `var` when `GetJump(var)` is next called. + void Recompute(int var); + // Returns true if the jump score for `var` might be negative. + bool PossiblyGood(int var) const; + + // Advanced usage, allows users to read possibly stale deltas for incremental + // score updates. + absl::Span Deltas() const { + return absl::MakeConstSpan(deltas_); + } + absl::Span Scores() const { + return absl::MakeConstSpan(scores_); + } + + absl::Span MutableScores() { return absl::MakeSpan(scores_); } + + // Note if you have very high weights (e.g. when using decay), the tolerances + // in this function are likely too tight. + bool JumpIsUpToDate(int var); // For debugging and testing. + + private: + absl::AnyInvocable(int)> compute_jump_; + + // For each variable, we store: + // - A jump delta which represents a change in variable value: + // new_value = old + delta, which is non-zero except for fixed variables. + // - The associated weighted feasibility violation change if we take this + // jump. + std::vector deltas_; + std::vector scores_; + std::vector needs_recomputation_; +}; + // Implements and heuristic similar to the one described in the paper: // "Feasibility Jump: an LP-free Lagrangian MIP heuristic", Bjørnar // Luteberget, Giorgio Sartor, 2023, Mathematical Programming Computation. @@ -62,7 +113,11 @@ class FeasibilityJumpSolver : public SubSolver { shared_response_(shared_response), shared_bounds_(shared_bounds), shared_stats_(shared_stats), - random_(params_) {} + random_(params_), + linear_jumps_( + absl::bind_front(&FeasibilityJumpSolver::ComputeLinearJump, this)), + general_jumps_(absl::bind_front( + &FeasibilityJumpSolver::ComputeGeneralJump, this)) {} // If VLOG_IS_ON(1), it will export a bunch of statistics. ~FeasibilityJumpSolver() override; @@ -103,51 +158,60 @@ class FeasibilityJumpSolver : public SubSolver { void PerturbateCurrentSolution(); std::string OneLineStats() const; - // Linear only. - bool JumpIsUpToDate(int var); // Debug. - void RecomputeJump(int var); - void RecomputeAllJumps(); - void MarkJumpsThatNeedsToBeRecomputed(int changed_var); + // Returns the weighted violation delta plus epsilon * the objective delta. + double ComputeScore(absl::Span scan_weights, int var, + int64_t delta, bool linear_only) const; + + // As above, but uses appropriate scan weights based on the current value of + // `use_compound_moves_`. + double ComputeScore(int var, int64_t delta, bool linear_only) const { + return ComputeScore(use_compound_moves_ ? compound_weights_ : weights_, var, + delta, linear_only); + } + // Computes the optimal value for variable v, considering only the violation + // of linear constraints. + std::pair ComputeLinearJump(int var); + + // Computes the optimal value for variable v, considering all constraints + // (assuming violation functions are convex). + std::pair ComputeGeneralJump(int var); + + // Marks all variables whose jump value may have changed due to the last + // update, except for `changed var`. + void MarkJumpsThatNeedToBeRecomputed(int changed_var, JumpTable& jumps); // Moves. bool DoSomeLinearIterations(); bool DoSomeGeneralIterations(); // Returns true if an improving move was found. - bool ScanRelevantVariables(int* improving_var, int64_t* improving_value, - double* improving_score, bool* time_limit_crossed); + bool ScanRelevantVariables(int num_to_scan, JumpTable& jumps, int* var, + int64_t* value, double* score); // Increases the weight of the currently infeasible constraints. - void UpdateViolatedConstraintWeights(); + // Ensures jumps remains consistent. + void UpdateViolatedConstraintWeights(JumpTable& jumps); - // Returns true if changing this variable to its jump value reduces weighted - // violation, or has no impact on weighted violation but reduces the - // objective. - bool IsGood(int var) const; + void UpdateNumViolatedConstraintsPerVar(); - void RecomputeVarsToScan(); + void RecomputeVarsToScan(JumpTable&); // Returns true if it is possible that `var` may have value that reduces // weighted violation or improve the objective. // Note that this is independent of the actual weights used. - bool ShouldScan(int var) const; - void AddVarToScan(int var); - bool MustRecomputeJumpOnGeneralUpdate(int var) const; + bool ShouldScan(const JumpTable& jumps, int var) const; + void AddVarToScan(const JumpTable&, int var); - // Resets the weights used to find compound moves. Only modifies weights that - // have changed, but the invariant below holds for all constraints: + // Resets the weights used to find compound moves. + // Ensures the following invariant holds afterwards: // compound_weights[c] = weights_[c] if c is violated, and epsilon * // weights_[c] otherwise. void ResetChangedCompoundWeights(); - // Returns true if we should accept a single variable move. - bool ShouldAcceptUnitMove(int var, int64_t delta, double score); - // Returns true if we should push this change to move_. // `novelty` is the total discount applied to the score due to using // `cumulative_weights_`, should always be positive (modulo floating-point // errors). - bool ShouldExtendCompoundMove(int var, int64_t delta, double score, - double novelty); + bool ShouldExtendCompoundMove(double score, double novelty); // Validates each element in num_violated_constraints_per_var_ against // evaluator_->ViolatedConstraints. @@ -168,19 +232,15 @@ class FeasibilityJumpSolver : public SubSolver { bool is_initialized_ = false; std::atomic model_is_supported_ = true; std::atomic task_generated_ = false; + bool time_limit_crossed_ = false; std::unique_ptr evaluator_; std::vector var_domains_; std::vector var_has_two_values_; + std::vector var_occurs_in_non_linear_constraint_; - // For each variable, we store: - // - A jump delta which represent a change in variable value: - // new_value = old + delta, which is non-zero except for fixed variable. - // - The associated weighted feasibility violation change if we take this - // jump. - std::vector jump_need_recomputation_; - std::vector jump_deltas_; - std::vector jump_scores_; + JumpTable linear_jumps_; + JumpTable general_jumps_; std::vector for_weight_update_; // The score of a solution is just the sum of infeasibility of each @@ -199,11 +259,9 @@ class FeasibilityJumpSolver : public SubSolver { // weight like for SAT activities. double bump_value_ = 1.0; - // Sparse list of jumps with a potential delta < 0.0. - // If jump_need_recomputation_[var] is true, we lazily recompute the exact - // delta as we randomly pick variables from here. - std::vector in_good_jumps_; - std::vector good_jumps_; + // A list of variables that might be relevant to check for improving jumps. + std::vector in_vars_to_scan_; + std::vector vars_to_scan_; // We restart each time our local deterministic time crosses this. double dtime_restart_threshold_ = 0.0; @@ -219,9 +277,14 @@ class FeasibilityJumpSolver : public SubSolver { // Each time we reset the weights, randomly decide if we will use compound // moves or not. bool use_compound_moves_ = false; - // We randomize the beam width each batch, choosing either 1 or 2, biased - // towards 1. - int compound_move_beam_search_width_ = 1; + + // Limit the discrepancy in compound move search (i.e. limit the number of + // backtracks to any ancestor of the current leaf). This is set to 0 whenever + // a new incumbent is found or weights are updated, and increased at fixed + // point. + // Weights are only increased if no moves are found with discrepancy 2. + // Empirically we have seen very few moves applied with discrepancy > 2. + int compound_move_max_discrepancy_ = 0; // Statistics int64_t num_batches_ = 0; @@ -235,14 +298,9 @@ class FeasibilityJumpSolver : public SubSolver { int64_t num_solutions_imported_ = 0; int64_t num_weight_updates_ = 0; - // A list of variables that might be relevant to check in general iterations. - std::vector in_vars_to_scan_; - std::vector vars_to_scan_; std::unique_ptr move_; // Counts the number of violated constraints each var is in. - // Only maintained in general iterations as jump_score_ is used for filtering - // moves in linear iterations. std::vector num_violated_constraints_per_var_; // Info on the last solution loaded. @@ -281,29 +339,15 @@ class CompoundMoveBuilder { return stack_.empty() ? 0.0 : stack_.back().cumulative_score; } - // Returns the sum of objective deltas from the root. - double ObjectiveDelta() const { - return stack_.empty() ? 0.0 : stack_.back().cumulative_objective_delta; - } - - // Returns the minimum score of any child of the current node in the stack, - // with a maximum of 0. - // NB: We assume all improving moves will immediately be committed, so this - // always returns 0 when the stack is empty. double BestChildScore() const { return stack_.empty() ? 0.0 : stack_.back().best_child_score; } - // Returns the number of times a child has been pushed to the current node. - // NB: Always returns 0 when the stack is empty. - int NumChildrenExplored() const { - return stack_.empty() ? 0 : stack_.back().num_children; + // Returns the number of backtracks to any ancestor of the current leaf. + int Discrepancy() const { + return stack_.empty() ? 0 : stack_.back().discrepancy; } - // Returns true if this move reduces weighted violation, or improves the - // objective without increasing violation. - bool IsImproving() const; - // Returns the number of backtracking moves that have been applied. int NumBacktracks() const { return num_backtracks_; } @@ -316,13 +360,12 @@ class CompoundMoveBuilder { // Instead of tracking this on the compound move, we store the partial sums // to avoid numerical issues causing negative scores after backtracking. double cumulative_score; - double cumulative_objective_delta; // Used to avoid infinite loops, this tracks the best score of any immediate // children (and not deeper descendants) to avoid re-exploring the same // prefix. double best_child_score = 0.0; - int num_children = 0; + int discrepancy = 0; }; LsEvaluator* evaluator_; std::vector var_on_stack_; diff --git a/ortools/sat/python/cp_model.py b/ortools/sat/python/cp_model.py index 3f127dfebc..a870ed8933 100644 --- a/ortools/sat/python/cp_model.py +++ b/ortools/sat/python/cp_model.py @@ -2146,10 +2146,11 @@ class CpModel: """Creates an optional interval var from start, size, end, and is_present. An optional interval variable is a constraint, that is itself used in other - constraints like NoOverlap. This constraint is protected by an is_present + constraints like NoOverlap. This constraint is protected by a presence literal that indicates if it is active or not. - Internally, it ensures that `is_present` implies `start + size == end`. + Internally, it ensures that `is_present` implies `start + size == + end`. Args: start: The start of the interval. It can be an integer value, or an @@ -2193,7 +2194,7 @@ class CpModel: starts: Union[LinearExprT, pd.Series], sizes: Union[LinearExprT, pd.Series], ends: Union[LinearExprT, pd.Series], - performed_literals: Union[LiteralT, pd.Series], + are_present: Union[LiteralT, pd.Series], ) -> pd.Series: """Creates a series of interval variables with the given name. @@ -2209,6 +2210,9 @@ class CpModel: ends (Union[LinearExprT, pd.Series]): The ends of each interval in the set. If a `pd.Series` is passed in, it will be based on the corresponding values of the pd.Series. + are_present (Union[LiteralT, pd.Series]): The performed literal of each + interval in the set. If a `pd.Series` is passed in, it will be based on + the corresponding values of the pd.Series. Returns: pd.Series: The interval variable set indexed by its corresponding @@ -2227,9 +2231,7 @@ class CpModel: starts = _ConvertToLinearExprSeriesAndValidateIndex(starts, index) sizes = _ConvertToLinearExprSeriesAndValidateIndex(sizes, index) ends = _ConvertToLinearExprSeriesAndValidateIndex(ends, index) - performed_literals = _ConvertToLiteralSeriesAndValidateIndex( - performed_literals, index - ) + are_present = _ConvertToLiteralSeriesAndValidateIndex(are_present, index) interval_array = [] for i in index: @@ -2238,14 +2240,18 @@ class CpModel: start=starts[i], size=sizes[i], end=ends[i], - is_present=performed_literals[i], + is_present=are_present[i], name=f"{name}[{i}]", ) ) return pd.Series(index=index, data=interval_array) def NewOptionalFixedSizeIntervalVar( - self, start: LinearExprT, size: IntegralT, is_present: LiteralT, name: str + self, + start: LinearExprT, + size: IntegralT, + is_present: LiteralT, + name: str, ) -> IntervalVar: """Creates an interval variable from start, and a fixed size. @@ -2273,7 +2279,12 @@ class CpModel: ) is_present_index = self.GetOrMakeBooleanIndex(is_present) return IntervalVar( - self.__model, start_expr, size_expr, end_expr, is_present_index, name + self.__model, + start_expr, + size_expr, + end_expr, + is_present_index, + name, ) def NewOptionalFixedSizeIntervalVarSeries( @@ -2282,7 +2293,7 @@ class CpModel: index: pd.Index, starts: Union[LinearExprT, pd.Series], sizes: Union[IntegralT, pd.Series], - performed_literals: Union[LiteralT, pd.Series], + are_present: Union[LiteralT, pd.Series], ) -> pd.Series: """Creates a series of interval variables with the given name. @@ -2295,9 +2306,9 @@ class CpModel: sizes (Union[IntegralT, pd.Series]): The fixed size of each interval in the set. If a `pd.Series` is passed in, it will be based on the corresponding values of the pd.Series. - performed_literals (Union[LiteralT, pd.Series]): The performed literal of - each interval in the set. If a `pd.Series` is passed in, it will be - based on the corresponding values of the pd.Series. + are_present (Union[LiteralT, pd.Series]): The performed literal of each + interval in the set. If a `pd.Series` is passed in, it will be based on + the corresponding values of the pd.Series. Returns: pd.Series: The interval variable set indexed by its corresponding @@ -2315,16 +2326,14 @@ class CpModel: starts = _ConvertToLinearExprSeriesAndValidateIndex(starts, index) sizes = _ConvertToIntegralSeriesAndValidateIndex(sizes, index) - performed_literals = _ConvertToLiteralSeriesAndValidateIndex( - performed_literals, index - ) + are_present = _ConvertToLiteralSeriesAndValidateIndex(are_present, index) interval_array = [] for i in index: interval_array.append( self.NewOptionalFixedSizeIntervalVar( start=starts[i], size=sizes[i], - is_present=performed_literals[i], + is_present=are_present[i], name=f"{name}[{i}]", ) ) diff --git a/ortools/sat/python/cp_model_test.py b/ortools/sat/python/cp_model_test.py index 934fd370fe..1afcbec251 100644 --- a/ortools/sat/python/cp_model_test.py +++ b/ortools/sat/python/cp_model_test.py @@ -1444,6 +1444,64 @@ class CpModelTest(absltest.TestCase): solution = solver.BooleanValues(x) self.assertTrue((solution.values == [False, True, False]).all()) + def testFixedSizeIntervalVarSeries(self): + print("testFixedSizeIntervalVarSeries") + df = pd.DataFrame([2, 4, 6], columns=["size"]) + model = cp_model.CpModel() + starts = model.NewIntVarSeries( + name="starts", index=df.index, lower_bounds=0, upper_bounds=5 + ) + presences = model.NewBoolVarSeries(name="rresences", index=df.index) + fixed_size_intervals = model.NewFixedSizeIntervalVarSeries( + name="fixed_size_intervals", + index=df.index, + starts=starts, + sizes=df.size, + ) + opt_fixed_size_intervals = model.NewOptionalFixedSizeIntervalVarSeries( + name="fixed_size_intervals", + index=df.index, + starts=starts, + sizes=df.size, + are_present=presences, + ) + model.AddNoOverlap( + fixed_size_intervals.to_list() + opt_fixed_size_intervals.to_list() + ) + self.assertLen(model.Proto().constraints, 7) + + def testIntervalVarSeries(self): + print("testIntervalVarSeries") + df = pd.DataFrame([2, 4, 6], columns=["size"]) + model = cp_model.CpModel() + starts = model.NewIntVarSeries( + name="starts", index=df.index, lower_bounds=0, upper_bounds=5 + ) + sizes = model.NewIntVarSeries( + name="sizes", index=df.index, lower_bounds=2, upper_bounds=4 + ) + ends = model.NewIntVarSeries( + name="ends", index=df.index, lower_bounds=0, upper_bounds=10 + ) + presences = model.NewBoolVarSeries(name="rresences", index=df.index) + intervals = model.NewIntervalVarSeries( + name="fixed_size_intervals", + index=df.index, + starts=starts, + sizes=sizes, + ends=ends, + ) + opt_intervals = model.NewOptionalIntervalVarSeries( + name="fixed_size_intervals", + index=df.index, + starts=starts, + sizes=sizes, + ends=ends, + are_present=presences, + ) + model.AddNoOverlap(intervals.to_list() + opt_intervals.to_list()) + self.assertLen(model.Proto().constraints, 13) + if __name__ == "__main__": absltest.main() diff --git a/ortools/sat/samples/assignment_sat.py b/ortools/sat/samples/assignment_sat.py index 4f84089722..f8cbfe7f3b 100644 --- a/ortools/sat/samples/assignment_sat.py +++ b/ortools/sat/samples/assignment_sat.py @@ -92,8 +92,10 @@ def main(): selected = data.loc[solver.BooleanValues(x).loc[lambda x: x].index] for unused_index, row in selected.iterrows(): print(f"{row.task} assigned to {row.worker} with a cost of {row.cost}") + elif status == cp_model.INFEASIBLE: + print("No solution found") else: - print("No solution found.") + print("Something is wrong, check the status and the log of the solve") # [END print_solution] diff --git a/ortools/sat/samples/bin_packing_sat.py b/ortools/sat/samples/bin_packing_sat.py index 2e6d761dd7..0efeb0315a 100644 --- a/ortools/sat/samples/bin_packing_sat.py +++ b/ortools/sat/samples/bin_packing_sat.py @@ -105,13 +105,13 @@ def main(): # [END objective] # [START solve] - # Create the solver with the CP-SAT backend, and solve the model. + # Create the solver and solve the model. solver = cp_model.CpSolver() status = solver.Solve(model) # [END solve] # [START print_solution] - if status == cp_model.OPTIMAL: + if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE: print(f"Number of bins used = {solver.ObjectiveValue()}") x_values = solver.BooleanValues(x) @@ -129,8 +129,10 @@ def main(): print(f"Total packed weight: {items.weight.sum()}") print() print(f"Time = {solver.WallTime()} seconds") + elif status == cp_model.INFEASIBLE: + print("No solution found") else: - print("The problem does not have an optimal solution.") + print("Something is wrong, check the status and the log of the solve") # [END print_solution] diff --git a/ortools/sat/samples/cumulative_variable_profile_sample_sat.py b/ortools/sat/samples/cumulative_variable_profile_sample_sat.py index eeeeb19aaf..49484443a0 100644 --- a/ortools/sat/samples/cumulative_variable_profile_sample_sat.py +++ b/ortools/sat/samples/cumulative_variable_profile_sample_sat.py @@ -117,7 +117,7 @@ def main(): index=tasks_df.index, starts=starts, sizes=tasks_df.duration, - performed_literals=performed, + are_present=performed, ) # [END variables] @@ -166,7 +166,7 @@ def main(): else: print(f"task {task} is not performed") elif status == cp_model.INFEASIBLE: - print("The problem is infeasible") + print("No solution found") else: print("Something is wrong, check the status and the log of the solve") # [END print_solution] diff --git a/ortools/sat/synchronization.h b/ortools/sat/synchronization.h index 40c32c85da..fedac37445 100644 --- a/ortools/sat/synchronization.h +++ b/ortools/sat/synchronization.h @@ -400,7 +400,7 @@ class SharedResponseManager { CpSolverStatus synchronized_best_status_ ABSL_GUARDED_BY(mutex_) = CpSolverStatus::UNKNOWN; std::vector unsat_cores_ ABSL_GUARDED_BY(mutex_); - SharedSolutionRepository solutions_ ABSL_GUARDED_BY(mutex_); + SharedSolutionRepository solutions_; // Thread-safe. int num_solutions_ ABSL_GUARDED_BY(mutex_) = 0; int64_t inner_objective_lower_bound_ ABSL_GUARDED_BY(mutex_) =