From dc3d9ccf846cbe5fe0aa7f9becb22a13cd3b0541 Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Thu, 5 Dec 2019 16:36:11 +0100 Subject: [PATCH] internal improvemensts --- ortools/base/random.h | 2 +- ortools/glop/lp_solver.cc | 11 +- ortools/glop/preprocessor.h | 3 - ortools/graph/perfect_matching.h | 2 + ortools/linear_solver/cbc_interface.cc | 4 +- ortools/linear_solver/gurobi_proto_solver.cc | 2 +- ortools/linear_solver/linear_solver.h | 4 +- .../linear_solver/samples/SimpleMipProgram.cs | 6 +- .../linear_solver/samples/mip_var_array.cc | 16 +- ortools/linear_solver/scip_proto_solver.cc | 2 +- ortools/sat/cp_model_expand.cc | 4 +- ortools/sat/cp_model_lns.cc | 171 +++++++++- ortools/sat/cp_model_lns.h | 100 ++++-- ortools/sat/cp_model_solver.cc | 113 +++---- ortools/sat/cuts.cc | 318 +++++++++++++----- ortools/sat/cuts.h | 7 +- ortools/sat/integer.h | 4 +- ortools/sat/linear_constraint_manager.cc | 6 + ortools/sat/linear_programming_constraint.cc | 165 ++++++++- ortools/sat/linear_programming_constraint.h | 1 + ortools/sat/python/cp_model.py | 36 +- ortools/sat/samples/default.profraw | Bin 2056 -> 0 bytes ortools/sat/sat_parameters.proto | 2 +- ortools/sat/swig_helper.h | 1 + 24 files changed, 753 insertions(+), 227 deletions(-) delete mode 100644 ortools/sat/samples/default.profraw diff --git a/ortools/base/random.h b/ortools/base/random.h index 643fc2ac8d..888d3d65f3 100644 --- a/ortools/base/random.h +++ b/ortools/base/random.h @@ -22,7 +22,7 @@ namespace operations_research { -// ACM minimal standard random number generator based on std::mt19937. +// A wrapper around std::mt19937. Called ACMRandom for historical reasons. class ACMRandom { public: explicit ACMRandom(int32 seed) : generator_(seed) {} diff --git a/ortools/glop/lp_solver.cc b/ortools/glop/lp_solver.cc index a3090f445e..d99ef512fd 100644 --- a/ortools/glop/lp_solver.cc +++ b/ortools/glop/lp_solver.cc @@ -600,7 +600,16 @@ bool LPSolver::IsProblemSolutionConsistent( ++num_basic_variables; break; case VariableStatus::FIXED_VALUE: - if (lb != ub || value != lb) { + // TODO(user): Because of scaling, it is possible that a FIXED_VALUE + // status (only reserved for the exact lb == ub case) is now set for a + // variable where (ub == lb + epsilon). So we do not check here that the + // two bounds are exactly equal. The best is probably to remove the + // FIXED status from the API completely and report one of AT_LOWER_BOUND + // or AT_UPPER_BOUND instead. This also allows to indicate if at + // optimality, the objective is limited because of this variable lower + // bound or its upper bound. Note that there are other TODOs in the + // codebase about removing this FIXED_VALUE status. + if (value != ub && value != lb) { LogVariableStatusError(col, value, status, lb, ub); return false; } diff --git a/ortools/glop/preprocessor.h b/ortools/glop/preprocessor.h index 0ae5661dde..d94daf2a01 100644 --- a/ortools/glop/preprocessor.h +++ b/ortools/glop/preprocessor.h @@ -246,9 +246,6 @@ class EmptyColumnPreprocessor : public Preprocessor { // -------------------------------------------------------- // ProportionalColumnPreprocessor // -------------------------------------------------------- -// TODO(user): For now this preprocessor just logs the number of proportional -// columns. Do something with this information. -// // Removes the proportional columns from the problem when possible. Two columns // are proportional if one is a non-zero scalar multiple of the other. // diff --git a/ortools/graph/perfect_matching.h b/ortools/graph/perfect_matching.h index 1203ec0d3d..659d9948ed 100644 --- a/ortools/graph/perfect_matching.h +++ b/ortools/graph/perfect_matching.h @@ -170,6 +170,8 @@ class BlossomGraph { DEFINE_INT_TYPE(CostValue, int64); // Basic constants. + // NOTE(user): Those can't be constexpr because of the or-tools export, + // which complains for constexpr DEFINE_INT_TYPE. static const NodeIndex kNoNodeIndex; static const EdgeIndex kNoEdgeIndex; static const CostValue kMaxCostValue; diff --git a/ortools/linear_solver/cbc_interface.cc b/ortools/linear_solver/cbc_interface.cc index 95db61612e..24e321ce9e 100644 --- a/ortools/linear_solver/cbc_interface.cc +++ b/ortools/linear_solver/cbc_interface.cc @@ -420,9 +420,9 @@ MPSolver::ResultStatus CBCInterface::Solve(const MPSolverParameters& param) { result_status_ = MPSolver::UNBOUNDED; } else if (model.isProvenInfeasible()) { result_status_ = MPSolver::INFEASIBLE; + } else if (model.isAbandoned()) { + result_status_ = MPSolver::ABNORMAL; } else { - VLOG(1) << "Unknown solver status! Secondary status: " - << model.secondaryStatus(); result_status_ = MPSolver::ABNORMAL; } break; diff --git a/ortools/linear_solver/gurobi_proto_solver.cc b/ortools/linear_solver/gurobi_proto_solver.cc index 9a6a8802cb..33a69bf85a 100644 --- a/ortools/linear_solver/gurobi_proto_solver.cc +++ b/ortools/linear_solver/gurobi_proto_solver.cc @@ -294,7 +294,7 @@ util::StatusOr GurobiSolveProto( request.solver_specific_parameters(), GRBgetenv(gurobi_model)); if (!parameters_status.ok()) { response.set_status(MPSOLVER_MODEL_INVALID_SOLVER_PARAMETERS); - response.set_status_str(parameters_status.error_message()); + response.set_status_str(parameters_status.message()); return response; } } diff --git a/ortools/linear_solver/linear_solver.h b/ortools/linear_solver/linear_solver.h index a8f214f0d6..edaf876ea4 100644 --- a/ortools/linear_solver/linear_solver.h +++ b/ortools/linear_solver/linear_solver.h @@ -674,7 +674,7 @@ class MPSolver { */ int64 nodes() const; - /// Returns a std::string describing the underlying solver and its version. + /// Returns a string describing the underlying solver and its version. std::string SolverVersion() const; /** @@ -1622,7 +1622,7 @@ class MPSolverInterface { return result_status_; } - // Returns a std::string describing the underlying solver and its version. + // Returns a string describing the underlying solver and its version. virtual std::string SolverVersion() const = 0; // Returns the underlying solver. diff --git a/ortools/linear_solver/samples/SimpleMipProgram.cs b/ortools/linear_solver/samples/SimpleMipProgram.cs index ff28e7758a..895b31f34e 100644 --- a/ortools/linear_solver/samples/SimpleMipProgram.cs +++ b/ortools/linear_solver/samples/SimpleMipProgram.cs @@ -51,15 +51,15 @@ public class SimpleMipProgram // [START solve] Solver.ResultStatus resultStatus = solver.Solve(); + // [END solve] + + // [START print_solution] // Check that the problem has an optimal solution. if (resultStatus != Solver.ResultStatus.OPTIMAL) { Console.WriteLine("The problem does not have an optimal solution!"); return; } - // [END solve] - - // [START print_solution] Console.WriteLine("Solution:"); Console.WriteLine("Objective value = " + solver.Objective().Value()); Console.WriteLine("x = " + x.SolutionValue()); diff --git a/ortools/linear_solver/samples/mip_var_array.cc b/ortools/linear_solver/samples/mip_var_array.cc index c2acda30b3..70f5ced1fd 100644 --- a/ortools/linear_solver/samples/mip_var_array.cc +++ b/ortools/linear_solver/samples/mip_var_array.cc @@ -11,7 +11,6 @@ // See the License for the specific language governing permissions and // limitations under the License. -// Mixed Integer programming example that shows how to use the API. // [START program] // [START import] #include "ortools/linear_solver/linear_solver.h" @@ -39,12 +38,9 @@ void IntegerProgrammingExample() { // [END data] // [START solver] - // MOE:begin_strip // Create the mip solver with the CBC backend. MPSolver solver("simple_mip_program", MPSolver::CBC_MIXED_INTEGER_PROGRAMMING); - MPSolver solver("simple_mip_program", - MPSolver::CBC_MIXED_INTEGER_PROGRAMMING); */ // [END solver] // [START variables] @@ -54,6 +50,7 @@ void IntegerProgrammingExample() { for (int j = 0; j < data.num_vars; ++j) { x[j] = solver.MakeIntVar(0.0, infinity, ""); } + LOG(INFO) << "Number of variables = " << solver.NumVariables(); // [END variables] // [START constraints] @@ -64,28 +61,27 @@ void IntegerProgrammingExample() { constraint->SetCoefficient(x[j], data.constraint_coeffs[i][j]); } } + LOG(INFO) << "Number of constraints = " << solver.NumConstraints(); // [END constraints] // [START objective] // Create the objective function. MPObjective* const objective = solver.MutableObjective(); - for (int j = 0; j < data.num_vars; ++j) { objective->SetCoefficient(x[j], data.obj_coeffs[j]); } objective->SetMaximization(); // [END objective] - // [START print_solution] - LOG(INFO) << "Number of variables = " << solver.NumVariables(); - LOG(INFO) << "Number of constraints = " << solver.NumConstraints(); - + // [START solve] const MPSolver::ResultStatus result_status = solver.Solve(); + // [END solve] + + // [START print_solution] // Check that the problem has an optimal solution. if (result_status != MPSolver::OPTIMAL) { LOG(FATAL) << "The problem does not have an optimal solution."; } - LOG(INFO) << "Solution:"; LOG(INFO) << "Optimal objective value = " << objective->Value(); diff --git a/ortools/linear_solver/scip_proto_solver.cc b/ortools/linear_solver/scip_proto_solver.cc index cf0d289746..a26caf5fca 100644 --- a/ortools/linear_solver/scip_proto_solver.cc +++ b/ortools/linear_solver/scip_proto_solver.cc @@ -773,7 +773,7 @@ util::StatusOr ScipSolveProto( request.solver_specific_parameters(), scip); if (!parameters_status.ok()) { response.set_status(MPSOLVER_MODEL_INVALID_SOLVER_PARAMETERS); - response.set_status_str(parameters_status.error_message()); + response.set_status_str(parameters_status.message()); return response; } // Default clock type. We use wall clock time because getting CPU user seconds diff --git a/ortools/sat/cp_model_expand.cc b/ortools/sat/cp_model_expand.cc index 7c2b22e7b3..eec9b91ec2 100644 --- a/ortools/sat/cp_model_expand.cc +++ b/ortools/sat/cp_model_expand.cc @@ -657,8 +657,8 @@ void ExpandAutomaton(ConstraintProto* ct, PresolveContext* context) { const std::vector vars = {proto.vars().begin(), proto.vars().end()}; // Compute the set of reachable state at each time point. - const absl::flat_hash_set final_states({proto.final_states().begin(), - proto.final_states().end()}); + const absl::flat_hash_set final_states( + {proto.final_states().begin(), proto.final_states().end()}); std::vector> reachable_states(n + 1); reachable_states[0].insert(proto.starting_state()); diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index 30d4fec7e2..f99a6ab943 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -14,6 +14,7 @@ #include "ortools/sat/cp_model_lns.h" #include +#include #include #include "ortools/sat/cp_model.pb.h" @@ -252,6 +253,7 @@ void NeighborhoodGenerator::Synchronize() { int num_not_fully_solved_in_batch = 0; for (const SolveData& data : solve_data_) { + AdditionalProcessingOnSynchronize(data); ++num_calls_; // INFEASIBLE or OPTIMAL means that we "fully solved" the local problem. @@ -272,7 +274,9 @@ void NeighborhoodGenerator::Synchronize() { // // This might not be a final solution, but it does work ok for now. const IntegerValue best_objective_improvement = - data.initial_best_objective - data.new_objective; + IsRelaxationGenerator() + ? (data.new_objective_bound - data.initial_best_objective_bound) + : (data.initial_best_objective - data.new_objective); if (best_objective_improvement > 0) { num_consecutive_non_improving_calls_ = 0; } else { @@ -331,7 +335,7 @@ void GetRandomSubset(double relative_size, std::vector* base, Neighborhood SimpleNeighborhoodGenerator::Generate( const CpSolverResponse& initial_solution, double difficulty, - random_engine_t* random) const { + random_engine_t* random) { std::vector fixed_variables = helper_.ActiveVariables(); GetRandomSubset(1.0 - difficulty, &fixed_variables, random); return helper_.FixGivenVariables(initial_solution, fixed_variables); @@ -339,7 +343,7 @@ Neighborhood SimpleNeighborhoodGenerator::Generate( Neighborhood VariableGraphNeighborhoodGenerator::Generate( const CpSolverResponse& initial_solution, double difficulty, - random_engine_t* random) const { + random_engine_t* random) { const int num_active_vars = helper_.ActiveVariables().size(); const int num_model_vars = helper_.ModelProto().variables_size(); const int target_size = std::ceil(difficulty * num_active_vars); @@ -390,7 +394,7 @@ Neighborhood VariableGraphNeighborhoodGenerator::Generate( Neighborhood ConstraintGraphNeighborhoodGenerator::Generate( const CpSolverResponse& initial_solution, double difficulty, - random_engine_t* random) const { + random_engine_t* random) { const int num_active_vars = helper_.ActiveVariables().size(); const int num_model_vars = helper_.ModelProto().variables_size(); const int target_size = std::ceil(difficulty * num_active_vars); @@ -543,7 +547,7 @@ Neighborhood GenerateSchedulingNeighborhoodForRelaxation( Neighborhood SchedulingNeighborhoodGenerator::Generate( const CpSolverResponse& initial_solution, double difficulty, - random_engine_t* random) const { + random_engine_t* random) { const auto span = helper_.TypeToConstraints(ConstraintProto::kInterval); std::vector intervals_to_relax(span.begin(), span.end()); GetRandomSubset(difficulty, &intervals_to_relax, random); @@ -554,7 +558,7 @@ Neighborhood SchedulingNeighborhoodGenerator::Generate( Neighborhood SchedulingTimeWindowNeighborhoodGenerator::Generate( const CpSolverResponse& initial_solution, double difficulty, - random_engine_t* random) const { + random_engine_t* random) { std::vector> start_interval_pairs; for (const int i : helper_.TypeToConstraints(ConstraintProto::kInterval)) { const ConstraintProto& interval_ct = helper_.ModelProto().constraints(i); @@ -588,7 +592,7 @@ bool RelaxationInducedNeighborhoodGenerator::ReadyToGenerate() const { Neighborhood RelaxationInducedNeighborhoodGenerator::Generate( const CpSolverResponse& initial_solution, double difficulty, - random_engine_t* random) const { + random_engine_t* random) { Neighborhood neighborhood = helper_.FullNeighborhood(); neighborhood.is_generated = false; @@ -645,12 +649,12 @@ Neighborhood RelaxationInducedNeighborhoodGenerator::Generate( return neighborhood; } -Neighborhood RandomRelaxationNeighborhoodGenerator::Generate( +Neighborhood ConsecutiveConstraintsRelaxationNeighborhoodGenerator::Generate( const CpSolverResponse& initial_solution, double difficulty, - random_engine_t* random) const { - std::vector removed_constraints; + random_engine_t* random) { + std::vector removable_constraints; const int num_constraints = helper_.ModelProto().constraints_size(); - removed_constraints.reserve(num_constraints); + removable_constraints.reserve(num_constraints); for (int c = 0; c < num_constraints; ++c) { // Removing intervals is not easy because other constraint might require // them, so for now, we don't remove them. @@ -658,12 +662,153 @@ Neighborhood RandomRelaxationNeighborhoodGenerator::Generate( ConstraintProto::kInterval) { continue; } - removed_constraints.push_back(c); + removable_constraints.push_back(c); + } + + const int target_size = + std::round((1.0 - difficulty) * removable_constraints.size()); + + std::uniform_int_distribution random_var( + 0, removable_constraints.size() - 1); + const int random_start_index = random_var(*random); + std::vector removed_constraints; + removed_constraints.reserve(target_size); + int c = random_start_index; + while (removed_constraints.size() < target_size) { + removed_constraints.push_back(removable_constraints[c]); + ++c; + if (c == removable_constraints.size()) { + c = 0; + } } - GetRandomSubset(1.0 - difficulty, &removed_constraints, random); return helper_.RemoveMarkedConstraints(removed_constraints); } +WeightedRandomRelaxationNeighborhoodGenerator:: + WeightedRandomRelaxationNeighborhoodGenerator( + NeighborhoodGeneratorHelper const* helper, const std::string& name) + : NeighborhoodGenerator(name, helper) { + std::vector removable_constraints; + const int num_constraints = helper_.ModelProto().constraints_size(); + constraint_weights_.reserve(num_constraints); + // TODO(user): Experiment with different starting weights. + for (int c = 0; c < num_constraints; ++c) { + switch (helper_.ModelProto().constraints(c).constraint_case()) { + case ConstraintProto::kCumulative: + case ConstraintProto::kAllDiff: + case ConstraintProto::kElement: + case ConstraintProto::kRoutes: + case ConstraintProto::kCircuit: + case ConstraintProto::kCircuitCovering: + constraint_weights_.push_back(3.0); + num_removable_constraints_++; + break; + case ConstraintProto::kBoolOr: + case ConstraintProto::kBoolAnd: + case ConstraintProto::kBoolXor: + case ConstraintProto::kIntProd: + case ConstraintProto::kIntDiv: + case ConstraintProto::kIntMod: + case ConstraintProto::kIntMax: + case ConstraintProto::kIntMin: + case ConstraintProto::kNoOverlap: + case ConstraintProto::kNoOverlap2D: + constraint_weights_.push_back(2.0); + num_removable_constraints_++; + break; + case ConstraintProto::kLinear: + case ConstraintProto::kTable: + case ConstraintProto::kAutomaton: + case ConstraintProto::kInverse: + case ConstraintProto::kReservoir: + case ConstraintProto::kAtMostOne: + constraint_weights_.push_back(1.0); + num_removable_constraints_++; + break; + case ConstraintProto::CONSTRAINT_NOT_SET: + case ConstraintProto::kInterval: + // Removing intervals is not easy because other constraint might require + // them, so for now, we don't remove them. + constraint_weights_.push_back(0.0); + break; + } + } +} + +void WeightedRandomRelaxationNeighborhoodGenerator:: + AdditionalProcessingOnSynchronize(const SolveData& solve_data) { + const IntegerValue best_objective_improvement = + solve_data.new_objective_bound - solve_data.initial_best_objective_bound; + + const std::vector& removed_constraints = + removed_constraints_[solve_data.neighborhood_id]; + + // Heuristic: We change the weights of the removed constraints if the + // neighborhood is solved (status is OPTIMAL or INFEASIBLE) or we observe an + // improvement in objective bounds. Otherwise we assume that the + // difficulty/time wasn't right for us to record feedbacks. + // + // If the objective bounds are improved, we bump up the weights. If the + // objective bounds are worse and the problem status is OPTIMAL, we bump down + // the weights. Otherwise if the new objective bounds are same as current + // bounds (which happens a lot on some instances), we do not update the + // weights as we do not have a clear signal wheather the constraints removed + // were good choices or not. + // TODO(user): We can improve this heuristic with more experiments. + if (best_objective_improvement > 0) { + // Bump up the weights of all removed constraints. + for (int c : removed_constraints) { + if (constraint_weights_[c] <= 90.0) { + constraint_weights_[c] += 10.0; + } else { + constraint_weights_[c] = 100.0; + } + } + } else if (solve_data.status == CpSolverStatus::OPTIMAL && + best_objective_improvement < 0) { + // Bump down the weights of all removed constraints. + for (int c : removed_constraints) { + if (constraint_weights_[c] > 0.5) { + constraint_weights_[c] -= 0.5; + } + } + } + removed_constraints_.erase(solve_data.neighborhood_id); +} + +Neighborhood WeightedRandomRelaxationNeighborhoodGenerator::Generate( + const CpSolverResponse& initial_solution, double difficulty, + random_engine_t* random) { + const int target_size = + std::round((1.0 - difficulty) * num_removable_constraints_); + + std::vector removed_constraints; + + // Generate a random number between (0,1) = u[i] and use score[i] = + // u[i]^(1/w[i]) and then select top k items with largest scores. + // Reference: https://utopia.duth.gr/~pefraimi/research/data/2007EncOfAlg.pdf + std::vector> constraint_removal_scores; + std::uniform_real_distribution random_var(0.0, 1.0); + for (int c = 0; c < constraint_weights_.size(); ++c) { + if (constraint_weights_[c] <= 0) continue; + const double u = random_var(*random); + const double score = std::pow(u, (1 / constraint_weights_[c])); + constraint_removal_scores.push_back({score, c}); + } + std::sort(constraint_removal_scores.rbegin(), + constraint_removal_scores.rend()); + for (int i = 0; i < target_size; ++i) { + removed_constraints.push_back(constraint_removal_scores[i].second); + } + + Neighborhood result = helper_.RemoveMarkedConstraints(removed_constraints); + absl::MutexLock mutex_lock(&mutex_); + result.id = next_available_id_; + next_available_id_++; + removed_constraints_.insert({result.id, removed_constraints}); + return result; +} + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/cp_model_lns.h b/ortools/sat/cp_model_lns.h index 7793f40f32..2a0ad7eda0 100644 --- a/ortools/sat/cp_model_lns.h +++ b/ortools/sat/cp_model_lns.h @@ -16,6 +16,7 @@ #include +#include "absl/container/flat_hash_map.h" #include "absl/synchronization/mutex.h" #include "absl/types/span.h" #include "ortools/base/integral_types.h" @@ -42,6 +43,12 @@ struct Neighborhood { // Relaxed model. Any feasible solution to this "local" model should be a // feasible solution to the base model too. CpModelProto cp_model; + + // Neighborhood Id. Used to identify the neighborhood by a generator. + // Currently only used by WeightedRandomRelaxationNeighborhoodGenerator. + // TODO(user): Make sure that the id is unique for each generated + // neighborhood for each generator. + int64 id = 0; }; // Contains pre-computed information about a given CpModelProto that is meant @@ -114,7 +121,7 @@ class NeighborhoodGeneratorHelper : public SubSolver { const CpModelProto& ModelProto() const { return model_proto_; } const SatParameters& Parameters() const { return parameters_; } - // This mutex must be aquired before calling any of the function that access + // This mutex must be acquired before calling any of the function that access // data that can be updated by Synchronize(). // // TODO(user): Refactor the class to be thread-safe instead, it should be @@ -181,14 +188,13 @@ class NeighborhoodGenerator { // will be dynamically adjusted depending on whether or not we can solve the // subproblem in a given time limit. // - // The given initial_solution should contains a feasible solution to the + // The given initial_solution should contain a feasible solution to the // initial CpModelProto given to this class. Any solution to the returned // CPModelProto should also be valid solution to the same initial model. // // This function should be thread-safe. virtual Neighborhood Generate(const CpSolverResponse& initial_solution, - double difficulty, - random_engine_t* random) const = 0; + double difficulty, random_engine_t* random) = 0; // Returns true if the neighborhood generator can generate a neighborhood. virtual bool ReadyToGenerate() const; @@ -209,6 +215,10 @@ class NeighborhoodGenerator { // Adds solve data about one "solved" neighborhood. struct SolveData { + // Neighborhood Id. Used to identify the neighborhood by a generator. + // Currently only used by WeightedRandomRelaxationNeighborhoodGenerator. + int64 neighborhood_id = 0; + // The status of the sub-solve. CpSolverStatus status = CpSolverStatus::UNKNOWN; @@ -234,14 +244,22 @@ class NeighborhoodGenerator { IntegerValue base_objective = IntegerValue(0); IntegerValue new_objective = IntegerValue(0); + // Bounds data is only used by relaxation neighborhoods. + IntegerValue initial_best_objective_bound = IntegerValue(0); + IntegerValue new_objective_bound = IntegerValue(0); + // This is just used to construct a deterministic order for the updates. bool operator<(const SolveData& o) const { return std::tie(status, difficulty, deterministic_limit, deterministic_time, initial_best_objective, - base_objective, new_objective) < + base_objective, new_objective, + initial_best_objective_bound, new_objective_bound, + neighborhood_id) < std::tie(o.status, o.difficulty, o.deterministic_limit, o.deterministic_time, o.initial_best_objective, - base_objective, new_objective); + o.base_objective, o.new_objective, + o.initial_best_objective_bound, o.new_objective_bound, + o.neighborhood_id); } }; void AddSolveData(SolveData data) { @@ -287,11 +305,16 @@ class NeighborhoodGenerator { } protected: + // Triggered with each call to Synchronize() for each recently added + // SolveData. This is meant to be used for processing feedbacks by specific + // neighborhood generators to adjust the neighborhood generation process. + virtual void AdditionalProcessingOnSynchronize(const SolveData& solve_data) {} + const std::string name_; const NeighborhoodGeneratorHelper& helper_; + mutable absl::Mutex mutex_; private: - mutable absl::Mutex mutex_; std::vector solve_data_; // Current parameters to be used when generating/solving a neighborhood with @@ -315,7 +338,7 @@ class SimpleNeighborhoodGenerator : public NeighborhoodGenerator { NeighborhoodGeneratorHelper const* helper, const std::string& name) : NeighborhoodGenerator(name, helper) {} Neighborhood Generate(const CpSolverResponse& initial_solution, - double difficulty, random_engine_t* random) const final; + double difficulty, random_engine_t* random) final; }; // Pick a random subset of variables that are constructed by a BFS in the @@ -328,7 +351,7 @@ class VariableGraphNeighborhoodGenerator : public NeighborhoodGenerator { NeighborhoodGeneratorHelper const* helper, const std::string& name) : NeighborhoodGenerator(name, helper) {} Neighborhood Generate(const CpSolverResponse& initial_solution, - double difficulty, random_engine_t* random) const final; + double difficulty, random_engine_t* random) final; }; // Pick a random subset of constraint and relax all of their variables. We are a @@ -341,7 +364,7 @@ class ConstraintGraphNeighborhoodGenerator : public NeighborhoodGenerator { NeighborhoodGeneratorHelper const* helper, const std::string& name) : NeighborhoodGenerator(name, helper) {} Neighborhood Generate(const CpSolverResponse& initial_solution, - double difficulty, random_engine_t* random) const final; + double difficulty, random_engine_t* random) final; }; // Helper method for the scheduling neighborhood generators. Returns the model @@ -364,7 +387,7 @@ class SchedulingNeighborhoodGenerator : public NeighborhoodGenerator { : NeighborhoodGenerator(name, helper) {} Neighborhood Generate(const CpSolverResponse& initial_solution, - double difficulty, random_engine_t* random) const final; + double difficulty, random_engine_t* random) final; }; // Similar to SchedulingNeighborhoodGenerator except the set of intervals that @@ -376,7 +399,7 @@ class SchedulingTimeWindowNeighborhoodGenerator : public NeighborhoodGenerator { : NeighborhoodGenerator(name, helper) {} Neighborhood Generate(const CpSolverResponse& initial_solution, - double difficulty, random_engine_t* random) const final; + double difficulty, random_engine_t* random) final; }; // Generates a neighborhood by fixing the variables who have same solution value @@ -397,7 +420,7 @@ class RelaxationInducedNeighborhoodGenerator : public NeighborhoodGenerator { : NeighborhoodGenerator(name, helper), model_(model) {} Neighborhood Generate(const CpSolverResponse& initial_solution, - double difficulty, random_engine_t* random) const final; + double difficulty, random_engine_t* random) final; // Returns true if SharedRINSNeighborhoodManager has unexplored neighborhoods. bool ReadyToGenerate() const override; @@ -405,21 +428,60 @@ class RelaxationInducedNeighborhoodGenerator : public NeighborhoodGenerator { const Model* model_; }; -// Generates a relaxation of the original model by randomly removing some -// constraints. The number of constraints removed is in sync with the difficulty -// passed to the generator. -class RandomRelaxationNeighborhoodGenerator : public NeighborhoodGenerator { +// Generates a relaxation of the original model by removing a consecutive span +// of constraints starting at a random index. The number of constraints removed +// is in sync with the difficulty passed to the generator. +class ConsecutiveConstraintsRelaxationNeighborhoodGenerator + : public NeighborhoodGenerator { public: - explicit RandomRelaxationNeighborhoodGenerator( + explicit ConsecutiveConstraintsRelaxationNeighborhoodGenerator( NeighborhoodGeneratorHelper const* helper, const std::string& name) : NeighborhoodGenerator(name, helper) {} Neighborhood Generate(const CpSolverResponse& initial_solution, - double difficulty, random_engine_t* random) const final; + double difficulty, random_engine_t* random) final; bool IsRelaxationGenerator() const override { return true; } bool ReadyToGenerate() const override { return true; } }; +// Generates a relaxation of the original model by removing some constraints +// randomly with a given weight for each constraint that controls the +// probability of constraint getting removed. The number of constraints removed +// is in sync with the difficulty passed to the generator. Higher weighted +// constraints are more likely to get removed. +class WeightedRandomRelaxationNeighborhoodGenerator + : public NeighborhoodGenerator { + public: + WeightedRandomRelaxationNeighborhoodGenerator( + NeighborhoodGeneratorHelper const* helper, const std::string& name); + + // Generates the neighborhood as described above. Also stores the removed + // constraints indices for adjusting the weights. + Neighborhood Generate(const CpSolverResponse& initial_solution, + double difficulty, random_engine_t* random) final; + + bool IsRelaxationGenerator() const override { return true; } + bool ReadyToGenerate() const override { return true; } + + private: + // Adjusts the weights of the constraints removed to get the neighborhood + // based on the solve_data. + void AdditionalProcessingOnSynchronize(const SolveData& solve_data) override + EXCLUSIVE_LOCKS_REQUIRED(mutex_); + + // Higher weighted constraints are more likely to get removed. + std::vector constraint_weights_; + int num_removable_constraints_ = 0; + + // Indices of the removed constraints per generated neighborhood. + absl::flat_hash_map> removed_constraints_ + GUARDED_BY(mutex_); + + // TODO(user): Move this to parent class if other generators start using + // feedbacks. + int64 next_available_id_ GUARDED_BY(mutex_) = 0; +}; + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index abd1271a12..bfd0968782 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -2017,19 +2017,14 @@ class LnsSolver : public SubSolver { postsolve_mapping, shared_->wall_timer, &local_response); if (generator_->IsRelaxationGenerator()) { + data.neighborhood_id = neighborhood.id; data.status = local_response.status(); data.deterministic_time = local_response.deterministic_time(); - data.new_objective = data.base_objective; - // TODO(user): The objective value might not be a good signal to - // adjust difficulty. Use bounds instead. bool has_feasible_solution = false; if (local_response.status() == CpSolverStatus::OPTIMAL || local_response.status() == CpSolverStatus::FEASIBLE) { - data.new_objective = IntegerValue(ComputeInnerObjective( - shared_->model_proto->objective(), local_response)); has_feasible_solution = true; } - generator_->AddSolveData(data); if (local_response.status() == CpSolverStatus::INFEASIBLE) { shared_->response->NotifyThatImprovingProblemIsInfeasible( @@ -2043,27 +2038,17 @@ class LnsSolver : public SubSolver { const IntegerValue local_obj_lb = local_response_manager.GetInnerObjectiveLowerBound(); - const bool is_maximization = - (shared_->model_proto->objective().scaling_factor() < 0); - - const double scaled_current_obj_bound = ScaleObjectiveValue( - shared_->model_proto->objective(), current_obj_lb.value()); const double scaled_local_obj_bound = ScaleObjectiveValue( neighborhood.cp_model.objective(), local_obj_lb.value()); - // If the objective bounds are not improving, abort early. - if ((is_maximization && - scaled_local_obj_bound > scaled_current_obj_bound) || - (!is_maximization && - scaled_local_obj_bound < scaled_current_obj_bound)) { - return; - } - // Update the bound. const IntegerValue new_inner_obj_lb = IntegerValue( std::ceil(UnscaleObjectiveValue(shared_->model_proto->objective(), scaled_local_obj_bound) - 1e-6)); + data.new_objective_bound = new_inner_obj_lb; + data.initial_best_objective_bound = current_obj_lb; + if (new_inner_obj_lb > current_obj_lb) { const IntegerValue current_obj_ub = shared_->response->GetInnerObjectiveUpperBound(); @@ -2079,13 +2064,13 @@ class LnsSolver : public SubSolver { absl::StrCat(local_response.solution_info(), " ", solution_info)); } + // If we have a solution of the relaxed problem, we check if it is also + // a valid solution of the non-relaxed one. if (has_feasible_solution && SolutionIsFeasible( *shared_->model_proto, std::vector(local_response.solution().begin(), local_response.solution().end()))) { - // If we have a solution of the relaxed problem, we check if it is - // also a valid solution of the non-relaxed one. shared_->response->NewSolution(local_response, /*model=*/nullptr); @@ -2096,34 +2081,47 @@ class LnsSolver : public SubSolver { shared_->time_limit->Stop(); } } - return; - } - - if (!local_response.solution().empty()) { - CHECK(SolutionIsFeasible( - *shared_->model_proto, - std::vector(local_response.solution().begin(), - local_response.solution().end()))) - << solution_info; - } - - local_response_manager.BestSolutionInnerObjectiveValue(); - if (local_response.solution_info().empty()) { - local_response.set_solution_info(solution_info); } else { - local_response.set_solution_info( - absl::StrCat(local_response.solution_info(), " ", solution_info)); + if (!local_response.solution().empty()) { + CHECK(SolutionIsFeasible( + *shared_->model_proto, + std::vector(local_response.solution().begin(), + local_response.solution().end()))) + << solution_info; + } + + if (local_response.solution_info().empty()) { + local_response.set_solution_info(solution_info); + } else { + local_response.set_solution_info( + absl::StrCat(local_response.solution_info(), " ", solution_info)); + } + + // Finish to fill the SolveData now that the local solve is done. + data.status = local_response.status(); + data.deterministic_time = local_response.deterministic_time(); + data.new_objective = data.base_objective; + if (local_response.status() == CpSolverStatus::OPTIMAL || + local_response.status() == CpSolverStatus::FEASIBLE) { + data.new_objective = IntegerValue(ComputeInnerObjective( + shared_->model_proto->objective(), local_response)); + } + + // Report any feasible solution we have. + if (local_response.status() == CpSolverStatus::OPTIMAL || + local_response.status() == CpSolverStatus::FEASIBLE) { + shared_->response->NewSolution(local_response, + /*model=*/nullptr); + } + if (!neighborhood.is_reduced && + (local_response.status() == CpSolverStatus::OPTIMAL || + local_response.status() == CpSolverStatus::INFEASIBLE)) { + shared_->response->NotifyThatImprovingProblemIsInfeasible( + local_response.solution_info()); + shared_->time_limit->Stop(); + } } - // Finish to fill the SolveData now that the local solve is done. - data.status = local_response.status(); - data.deterministic_time = local_response.deterministic_time(); - data.new_objective = data.base_objective; - if (local_response.status() == CpSolverStatus::OPTIMAL || - local_response.status() == CpSolverStatus::FEASIBLE) { - data.new_objective = IntegerValue(ComputeInnerObjective( - shared_->model_proto->objective(), local_response)); - } generator_->AddSolveData(data); // The total number of call when this was called is the same as task_id. @@ -2136,20 +2134,6 @@ class LnsSolver : public SubSolver { << ", num calls: " << generator_->num_calls() << ", UCB1 Score: " << generator_->GetUCBScore(total_num_calls) << ", p: " << fully_solved_proportion << "]"; - - // Report any feasible solution we have. - if (local_response.status() == CpSolverStatus::OPTIMAL || - local_response.status() == CpSolverStatus::FEASIBLE) { - shared_->response->NewSolution(local_response, - /*model=*/nullptr); - } - if (!neighborhood.is_reduced && - (local_response.status() == CpSolverStatus::OPTIMAL || - local_response.status() == CpSolverStatus::INFEASIBLE)) { - shared_->response->NotifyThatImprovingProblemIsInfeasible( - local_response.solution_info()); - shared_->time_limit->Stop(); - } }; } @@ -2279,9 +2263,16 @@ void SolveCpModelParallel(const CpModelProto& model_proto, if (parameters.use_relaxation_lns()) { subsolvers.push_back(absl::make_unique( /*id=*/subsolvers.size(), - absl::make_unique( + absl::make_unique< + ConsecutiveConstraintsRelaxationNeighborhoodGenerator>( helper, absl::StrCat("rnd_rel_lns_", strategy_name)), local_params, helper, &shared)); + + subsolvers.push_back(absl::make_unique( + /*id=*/subsolvers.size(), + absl::make_unique( + helper, absl::StrCat("wgt_rel_lns_", strategy_name)), + local_params, helper, &shared)); } if (!helper->TypeToConstraints(ConstraintProto::kNoOverlap).empty()) { diff --git a/ortools/sat/cuts.cc b/ortools/sat/cuts.cc index db0b451acd..a4cb6c2bc0 100644 --- a/ortools/sat/cuts.cc +++ b/ortools/sat/cuts.cc @@ -22,6 +22,7 @@ #include "ortools/algorithms/knapsack_solver_for_cuts.h" #include "ortools/base/integral_types.h" +#include "ortools/base/stl_util.h" #include "ortools/sat/integer.h" #include "ortools/sat/linear_constraint.h" #include "ortools/util/time_limit.h" @@ -113,7 +114,7 @@ bool AllVarsTakeIntegerValue( // bound. // 2. Sort all terms (coefficient * shifted upper bound) in non decreasing // order. -// 3. Add terms in cover untill term sum is smaller or equal to upper bound. +// 3. Add terms in cover until term sum is smaller or equal to upper bound. // 4. Add the last item which violates the upper bound. This forms the smallest // cover. Return the size of this cover. int GetSmallestCoverSize(const LinearConstraint& constraint, @@ -611,7 +612,8 @@ std::function GetSuperAdditiveRoundingFunction( const IntegerValue size = divisor - rhs_remainder; if (max_scaling == 1) { // TODO(user): Use everywhere a two step computation to avoid overflow? - // First divide by divisor, then multiply by t. + // First divide by divisor, then multiply by t. For now, we limit the + // max_scaling so that we never have an overflow instead. return [t, divisor](IntegerValue coeff) { return FloorRatio(t * coeff, divisor); }; @@ -657,11 +659,13 @@ std::function GetSuperAdditiveRoundingFunction( } } -// TODO(user): This can currently take a significant portion of the run time on -// problem like opm2-z7-s8.mps.gz. Fix or just call less often. -void IntegerRoundingCut(RoundingOptions options, std::vector lp_values, - std::vector lower_bounds, - std::vector upper_bounds, +// TODO(user): This has been optimized a bit, but we can probably do even better +// as it still takes around 25% percent of the run time when all the cuts are on +// for the opm*mps.gz problems and others. +void IntegerRoundingCut(RoundingOptions options, + const std::vector& lp_values, + const std::vector& lower_bounds, + const std::vector& upper_bounds, LinearConstraint* cut) { const int size = lp_values.size(); if (size == 0) return; @@ -671,87 +675,151 @@ void IntegerRoundingCut(RoundingOptions options, std::vector lp_values, CHECK_EQ(cut->coeffs.size(), size); CHECK_EQ(cut->lb, kMinIntegerValue); - // Optimization: upper bound is no longer used once this is filled. - std::vector& bound_diffs = upper_bounds; + // To optimize the computation of the best divisor below, we only need to + // look at the indices with a shifted lp value that is not close to zero. + // + // TODO(user): use a class to reuse this memory. Note however that currently + // this do not appear in the cpu profile. + // + // TODO(user): sort by decreasing lp_values so that our early abort test in + // the critical loop below has more chance of returning early? I tried but it + // didn't seems to change much though. + std::vector relevant_indices; + std::vector relevant_lp_values; + std::vector relevant_coeffs; + std::vector relevant_bound_diffs; + std::vector divisors; + std::vector> adjusted_coeffs; // Shift each variable using its lower/upper bound so that no variable can // change sign. We eventually do a change of variable to its negation so // that all variable are non-negative. bool overflow = false; std::vector change_sign_at_postprocessing(size, false); - IntegerValue max_initial_magnitude(1); + IntegerValue max_magnitude(0); for (int i = 0; i < size; ++i) { if (cut->coeffs[i] == 0) continue; + // We might change them below. + IntegerValue lb = lower_bounds[i]; + double lp_value = lp_values[i]; + + const IntegerValue ub = upper_bounds[i]; + const IntegerValue bound_diff = + IntegerValue(CapSub(ub.value(), lb.value())); + // Note that since we use ToDouble() this code works fine with lb/ub at // min/max integer value. { - const double value = lp_values[i]; - const IntegerValue lb = lower_bounds[i]; - const IntegerValue ub = upper_bounds[i]; - bound_diffs[i] = IntegerValue(CapSub(ub.value(), lb.value())); - - if (std::abs(value - ToDouble(lb)) > std::abs(value - ToDouble(ub))) { + if (std::abs(lp_value - ToDouble(lb)) > + std::abs(lp_value - ToDouble(ub))) { // Change the variable sign. change_sign_at_postprocessing[i] = true; cut->coeffs[i] = -cut->coeffs[i]; - lp_values[i] = -lp_values[i]; - lower_bounds[i] = -ub; + lp_value = -lp_value; + lb = -ub; } } // Always shift to lb. // coeff * X = coeff * (X - shift) + coeff * shift. - lp_values[i] -= ToDouble(lower_bounds[i]); - if (!AddProductTo(-cut->coeffs[i], lower_bounds[i], &cut->ub)) { + lp_value -= ToDouble(lb); + if (!AddProductTo(-cut->coeffs[i], lb, &cut->ub)) { overflow = true; break; } // Deal with fixed variable, no need to shift back in this case, we can // just remove the term. - if (bound_diffs[i] == 0) { + if (bound_diff == 0) { cut->coeffs[i] = IntegerValue(0); - lp_values[i] = 0.0; + lp_value = 0.0; } - max_initial_magnitude = - std::max(max_initial_magnitude, IntTypeAbs(cut->coeffs[i])); + const IntegerValue magnitude = IntTypeAbs(cut->coeffs[i]); + if (std::abs(lp_value) > 1e-2) { + relevant_coeffs.push_back(cut->coeffs[i]); + relevant_indices.push_back(i); + relevant_lp_values.push_back(lp_value); + relevant_bound_diffs.push_back(bound_diff); + + divisors.push_back(magnitude); + } + max_magnitude = std::max(max_magnitude, magnitude); } - // Make sure that when we multiply the rhs or the coefficient by max_scaling, - // we do not have an integer overflow. - const IntegerValue max_scaling = std::min( - options.max_scaling, - kMaxIntegerValue / std::max(IntTypeAbs(cut->ub), max_initial_magnitude)); - if (overflow || max_scaling == 0) { - VLOG(1) << "Issue, overflow."; + // TODO(user): Maybe this shouldn't be called on such constraint. + if (relevant_coeffs.empty()) { + VLOG(2) << "Issue, nothing to cut."; *cut = LinearConstraint(IntegerValue(0), IntegerValue(0)); return; } + CHECK_NE(max_magnitude, 0); // Our heuristic will try to generate a few different cuts, and we will keep - // the most violated one. + // the most violated one scaled by the l2 norm of the relevant position. + // + // TODO(user): Experiment for the best value of this initial violation + // threshold. Note also that we use the l2 norm on the restricted position + // here. Maybe we should change that? On that note, the L2 norm usage seems a + // bit weird to me since it grows with the number of term in the cut. And + // often, we already have a good cut, and we make it stronger by adding extra + // terms that do not change its activity. + // + // The discussion above only concern the best_scaled_violation initial value. + // The remainder_threshold allows to not consider cuts for which the final + // efficacity is clearly lower than 1e-3 (it is a bound, so we could generate + // cuts with a lower efficacity than this). double best_scaled_violation = 0.01; - LinearConstraint best_cut(IntegerValue(0), IntegerValue(0)); - LinearConstraint temp_cut; + const IntegerValue remainder_threshold(max_magnitude / 1000); - // TODO(user): Avoid quadratic algorithm. - for (int index = 0; index < size; ++index) { - // Skip shifted variable almost at their lower bound. - if (lp_values[index] <= 1e-2) continue; - const IntegerValue divisor = IntTypeAbs(cut->coeffs[index]); + // The cut->ub might have grown quite a bit with the bound substitution, so + // we need to include it too since we will apply the rounding function on it. + max_magnitude = std::max(max_magnitude, IntTypeAbs(cut->ub)); + // Make sure that when we multiply the rhs or the coefficient by max_scaling, + // we do not have an integer overflow. Actually, we need a bit more room + // because we might round down a value to the next multiple of + // max_magnitude. + const IntegerValue threshold = kMaxIntegerValue / 2; + if (overflow || max_magnitude >= threshold) { + VLOG(2) << "Issue, overflow."; + *cut = LinearConstraint(IntegerValue(0), IntegerValue(0)); + return; + } + const IntegerValue max_scaling = + std::min(options.max_scaling, threshold / max_magnitude); + + // There is no point trying twice the same divisor or s divisor that is too + // small. Note that we use a higher threshold than the remainder_threshold + // because we can boost the remainder thanks to our adjusting heuristic below + // and also because this allows to have cuts with a small range of + // coefficients. + // + // TODO(user): Note that the std::sort() is visible in some cpu profile. + { + int new_size = 0; + const IntegerValue divisor_threshold = max_magnitude / 10; + for (int i = 0; i < divisors.size(); ++i) { + if (divisors[i] <= divisor_threshold) continue; + divisors[new_size++] = divisors[i]; + } + divisors.resize(new_size); + } + gtl::STLSortAndRemoveDuplicates(&divisors, std::greater()); + + // TODO(user): Avoid quadratic algorithm? Note that we are quadratic in + // relevant_indices not the full cut->coeffs.size(), but this is still too + // much on some problems. + IntegerValue best_divisor(0); + for (const IntegerValue divisor : divisors) { // Skip if we don't have the potential to generate a good enough cut. const IntegerValue initial_rhs_remainder = cut->ub - FloorRatio(cut->ub, divisor) * divisor; - if (ToDouble(initial_rhs_remainder) / ToDouble(max_initial_magnitude) <= - best_scaled_violation) { - continue; - } + if (initial_rhs_remainder <= remainder_threshold) continue; - // TODO(user): by recomputing, we should avoid the need to make this copy. - temp_cut = *cut; + IntegerValue temp_ub = cut->ub; + adjusted_coeffs.clear(); // We will adjust coefficient that are just under an exact multiple of // divisor to an exact multiple. This is meant to get rid of small errors @@ -771,64 +839,156 @@ void IntegerRoundingCut(RoundingOptions options, std::vector lp_values, const IntegerValue adjust_threshold = (divisor - initial_rhs_remainder - 1) / IntegerValue(size); if (adjust_threshold > 0) { - for (int i = 0; i < size; ++i) { - const IntegerValue diff = bound_diffs[i]; - if (diff > adjust_threshold) continue; + // Even before we finish the adjust, we can have a lower bound on the + // activily loss using this divisor, and so we can abort early. This is + // similar to what is done below in the function. + bool early_abort = false; + double loss_lb = 0.0; + const double threshold = ToDouble(initial_rhs_remainder); - // Adjust coeff of the form k * divisor - epsilon. - const IntegerValue coeff = temp_cut.coeffs[i]; + for (int i = 0; i < relevant_coeffs.size(); ++i) { + // Compute the difference of coeff with the next multiple of divisor. + const IntegerValue coeff = relevant_coeffs[i]; const IntegerValue remainder = CeilRatio(coeff, divisor) * divisor - coeff; - if (CapProd(diff.value(), remainder.value()) > adjust_threshold) { - continue; + + if (divisor - remainder <= initial_rhs_remainder) { + // We do not know exactly f() yet, but it will always round to the + // floor of the division by divisor in this case. + loss_lb += ToDouble(divisor - remainder) * relevant_lp_values[i]; + if (loss_lb >= threshold) { + early_abort = true; + break; + } + } + + // Adjust coeff of the form k * divisor - epsilon. + const IntegerValue diff = relevant_bound_diffs[i]; + if (remainder > 0 && remainder <= adjust_threshold && + CapProd(diff.value(), remainder.value()) <= adjust_threshold) { + temp_ub += remainder * diff; + adjusted_coeffs.push_back({i, coeff + remainder}); } - temp_cut.ub += remainder * diff; - temp_cut.coeffs[i] += remainder; } + + if (early_abort) continue; } // Create the super-additive function f(). const IntegerValue rhs_remainder = - temp_cut.ub - FloorRatio(temp_cut.ub, divisor) * divisor; + temp_ub - FloorRatio(temp_ub, divisor) * divisor; if (rhs_remainder == 0) continue; const auto f = GetSuperAdditiveRoundingFunction( !options.use_mir, rhs_remainder, divisor, max_scaling); - // Apply f() to the cut and compute the cut violation. - temp_cut.ub = f(temp_cut.ub); - double violation = -ToDouble(temp_cut.ub); - double max_magnitude = 1.0; - for (int i = 0; i < temp_cut.coeffs.size(); ++i) { - const IntegerValue coeff = temp_cut.coeffs[i]; + // As we round coefficients, we will compute the loss compared to the + // current scaled constraint activity. As soon as this loss crosses the + // slack, then we known that there is no violation and we can abort early. + // + // TODO(user): modulo the scaling, we could compute the exact threshold + // using our current best cut. Note that we also have to account the change + // in slack due to the adjust code above. + const double scaling = ToDouble(f(divisor)) / ToDouble(divisor); + const double threshold = scaling * ToDouble(rhs_remainder); + double loss = 0.0; + + // Apply f() to the cut and compute the cut violation. Note that it is + // okay to just look at the relevant indices since the other have a lp + // value which is almost zero. Doing it like this is faster, and even if + // the max_magnitude might be off it should still be relevant enough. + double violation = -ToDouble(f(temp_ub)); + double l2_norm = 0.0; + bool early_abort = false; + int adjusted_coeffs_index = 0; + for (int i = 0; i < relevant_coeffs.size(); ++i) { + IntegerValue coeff = relevant_coeffs[i]; + + // Adjust coeff according to our previous computation if needed. + if (adjusted_coeffs_index < adjusted_coeffs.size() && + adjusted_coeffs[adjusted_coeffs_index].first == i) { + coeff = adjusted_coeffs[adjusted_coeffs_index].second; + adjusted_coeffs_index++; + } + if (coeff == 0) continue; const IntegerValue new_coeff = f(coeff); - temp_cut.coeffs[i] = new_coeff; - max_magnitude = std::max(max_magnitude, std::abs(ToDouble(new_coeff))); - violation += ToDouble(new_coeff) * lp_values[i]; - } - violation /= max_magnitude; + const double new_coeff_double = ToDouble(new_coeff); + const double lp_value = relevant_lp_values[i]; - if (violation > 0.0) { - VLOG(2) << "lp_value: " << lp_values[index] << " divisor: " << divisor - << " cut_violation: " << violation; + l2_norm += new_coeff_double * new_coeff_double; + violation += new_coeff_double * lp_value; + loss += (scaling * ToDouble(coeff) - new_coeff_double) * lp_value; + if (loss >= threshold) { + early_abort = true; + break; + } } + if (early_abort) continue; + + // Here we scale by the L2 norm over the "relevant" positions. This seems + // to work slighly better in practice. + violation /= sqrt(l2_norm); if (violation > best_scaled_violation) { best_scaled_violation = violation; - best_cut = temp_cut; + best_divisor = divisor; } } + if (best_divisor == 0) { + *cut = LinearConstraint(IntegerValue(0), IntegerValue(0)); + return; + } + + // Adjust coefficients. + // + // TODO(user): It might make sense to also adjust the one with a small LP + // value, but then the cut will be slighlty different than the one we computed + // above. Try with and without maybe? + const IntegerValue initial_rhs_remainder = + cut->ub - FloorRatio(cut->ub, best_divisor) * best_divisor; + const IntegerValue adjust_threshold = + (best_divisor - initial_rhs_remainder - 1) / IntegerValue(size); + if (adjust_threshold > 0) { + for (int i = 0; i < relevant_indices.size(); ++i) { + const int index = relevant_indices[i]; + const IntegerValue diff = relevant_bound_diffs[i]; + if (diff > adjust_threshold) continue; + + // Adjust coeff of the form k * best_divisor - epsilon. + const IntegerValue coeff = cut->coeffs[index]; + const IntegerValue remainder = + CeilRatio(coeff, best_divisor) * best_divisor - coeff; + if (CapProd(diff.value(), remainder.value()) <= adjust_threshold) { + cut->ub += remainder * diff; + cut->coeffs[index] += remainder; + } + } + } + + // Create the super-additive function f(). + const IntegerValue rhs_remainder = + cut->ub - FloorRatio(cut->ub, best_divisor) * best_divisor; + const auto f = GetSuperAdditiveRoundingFunction( + !options.use_mir, rhs_remainder, best_divisor, max_scaling); + + // Apply f() to the cut. + // // Remove the bound shifts so the constraint is expressed in the original // variables and do some basic post-processing. - *cut = best_cut; - for (int i = 0; i < cut->coeffs.size(); ++i) { - const IntegerValue coeff = cut->coeffs[i]; + cut->ub = f(cut->ub); + for (int i = 0; i < size; ++i) { + IntegerValue coeff = cut->coeffs[i]; + if (coeff == 0) continue; + cut->coeffs[i] = coeff = f(coeff); if (coeff == 0) continue; - cut->ub = IntegerValue( - CapAdd((coeff * lower_bounds[i]).value(), cut->ub.value())); if (change_sign_at_postprocessing[i]) { + cut->ub = IntegerValue( + CapAdd((coeff * -upper_bounds[i]).value(), cut->ub.value())); cut->coeffs[i] = -coeff; + } else { + cut->ub = IntegerValue( + CapAdd((coeff * lower_bounds[i]).value(), cut->ub.value())); } } RemoveZeroTerms(cut); @@ -1052,7 +1212,7 @@ void ImpliedBoundsProcessor::ProcessUpperBoundedConstraint( } if (CapProd(std::abs(coeff.value()), diff.value()) >= kMaxIntegerValue) { - VLOG(1) << "Overflow"; + VLOG(2) << "Overflow"; return; } @@ -1060,7 +1220,7 @@ void ImpliedBoundsProcessor::ProcessUpperBoundedConstraint( // X >= Indicator * (bound - lb) + lb tmp_terms_.push_back({entry.literal_view, coeff * diff}); if (!AddProductTo(-coeff, lb, &new_ub)) { - VLOG(1) << "Overflow"; + VLOG(2) << "Overflow"; return; } } else { @@ -1068,7 +1228,7 @@ void ImpliedBoundsProcessor::ProcessUpperBoundedConstraint( // X >= -Indicator * (bound - lb) + bound tmp_terms_.push_back({entry.literal_view, -coeff * diff}); if (!AddProductTo(-coeff, entry.lower_bound, &new_ub)) { - VLOG(1) << "Overflow"; + VLOG(2) << "Overflow"; return; } } @@ -1091,7 +1251,7 @@ void ImpliedBoundsProcessor::ProcessUpperBoundedConstraint( } if (overflow_detection >= kMaxIntegerValue) { - VLOG(1) << "Overflow"; + VLOG(2) << "Overflow"; return; } if (!changed) return; diff --git a/ortools/sat/cuts.h b/ortools/sat/cuts.h index 83a313e6ab..7b3089e1e0 100644 --- a/ortools/sat/cuts.h +++ b/ortools/sat/cuts.h @@ -145,9 +145,10 @@ struct RoundingOptions { bool use_mir = false; IntegerValue max_scaling = IntegerValue(60); }; -void IntegerRoundingCut(RoundingOptions options, std::vector lp_values, - std::vector lower_bounds, - std::vector upper_bounds, +void IntegerRoundingCut(RoundingOptions options, + const std::vector& lp_values, + const std::vector& lower_bounds, + const std::vector& upper_bounds, LinearConstraint* cut); // If a variable is away from its upper bound by more than value 1.0, then it diff --git a/ortools/sat/integer.h b/ortools/sat/integer.h index ae46ed8eec..5b4bce94c7 100644 --- a/ortools/sat/integer.h +++ b/ortools/sat/integer.h @@ -91,7 +91,7 @@ inline IntegerValue Subtract(IntegerValue a, IntegerValue b) { inline IntegerValue CeilRatio(IntegerValue dividend, IntegerValue positive_divisor) { - CHECK_GT(positive_divisor, 0); + DCHECK_GT(positive_divisor, 0); const IntegerValue result = dividend / positive_divisor; const IntegerValue adjust = static_cast(result * positive_divisor < dividend); @@ -100,7 +100,7 @@ inline IntegerValue CeilRatio(IntegerValue dividend, inline IntegerValue FloorRatio(IntegerValue dividend, IntegerValue positive_divisor) { - CHECK_GT(positive_divisor, 0); + DCHECK_GT(positive_divisor, 0); const IntegerValue result = dividend / positive_divisor; const IntegerValue adjust = static_cast(result * positive_divisor > dividend); diff --git a/ortools/sat/linear_constraint_manager.cc b/ortools/sat/linear_constraint_manager.cc index 7910c27ca2..ada7f8e752 100644 --- a/ortools/sat/linear_constraint_manager.cc +++ b/ortools/sat/linear_constraint_manager.cc @@ -402,6 +402,9 @@ bool LinearConstraintManager::ChangeLp( // Inprocessing of the constraint. if (simplify_constraints && SimplifyConstraint(&constraint_infos_[i].constraint)) { + // Note that the canonicalization shouldn't be needed since the order + // of the variable is not changed by the simplification, and we only + // reduce the coefficients at both end of the spectrum. DivideByGCD(&constraint_infos_[i].constraint); DCHECK(DebugCheckConstraint(constraint_infos_[i].constraint)); @@ -413,6 +416,9 @@ bool LinearConstraintManager::ChangeLp( equiv_constraints_.erase(constraint_infos_[i].hash); constraint_infos_[i].hash = ComputeHashOfTerms(constraint_infos_[i].constraint); + + // TODO(user): Because we simplified this constraint, it is possible that + // it is now a duplicate of another one. Merge them. equiv_constraints_[constraint_infos_[i].hash] = i; } diff --git a/ortools/sat/linear_programming_constraint.cc b/ortools/sat/linear_programming_constraint.cc index 54438f7ab1..5e4a94fce4 100644 --- a/ortools/sat/linear_programming_constraint.cc +++ b/ortools/sat/linear_programming_constraint.cc @@ -28,6 +28,8 @@ #include "ortools/base/integral_types.h" #include "ortools/base/logging.h" #include "ortools/base/map_util.h" +#include "ortools/base/mathutil.h" +#include "ortools/base/stl_util.h" #include "ortools/glop/parameters.pb.h" #include "ortools/glop/preprocessor.h" #include "ortools/glop/status.h" @@ -58,6 +60,7 @@ LinearProgrammingConstraint::LinearProgrammingConstraint(Model* model) trail_(model->GetOrCreate()), model_heuristics_(model->GetOrCreate()), integer_encoder_(model->GetOrCreate()), + random_(model->GetOrCreate()), implied_bounds_processor_({}, integer_trail_, model->GetOrCreate()), dispatcher_(model->GetOrCreate()), @@ -230,7 +233,7 @@ LPSolveInfo LinearProgrammingConstraint::SolveLpForBranching() { LPSolveInfo info; glop::BasisState basis_state = simplex_.GetState(); - const auto status = simplex_.Solve(lp_data_, time_limit_); + const glop::Status status = simplex_.Solve(lp_data_, time_limit_); total_num_simplex_iterations_ += simplex_.GetNumberOfIterations(); simplex_.LoadStateForNextSolve(basis_state); if (!status.ok()) { @@ -596,10 +599,11 @@ void LinearProgrammingConstraint::AddCutFromConstraints( cut = ConvertToLinearConstraint(dense_cut, cut_ub); } - // This should be tight! + // This should be tight. + const double norm = ToDouble(ComputeInfinityNorm(cut)); if (std::abs(ComputeActivity(cut, expanded_lp_solution_) - ToDouble(cut.ub)) / - std::max(1.0, std::abs(ToDouble(cut.ub))) > - 1e-2) { + norm > + 1e-4) { VLOG(1) << "Cut not tight " << ComputeActivity(cut, expanded_lp_solution_) << " " << ToDouble(cut.ub); return; @@ -668,8 +672,7 @@ void LinearProgrammingConstraint::AddCutFromConstraints( RoundingOptions options; options.use_mir = sat_parameters_.use_mir_rounding(); options.max_scaling = sat_parameters_.max_integer_rounding_scaling(); - IntegerRoundingCut(options, std::move(lp_values), std::move(var_lbs), - std::move(var_ubs), &cut); + IntegerRoundingCut(options, lp_values, var_lbs, var_ubs, &cut); // Compute the activity. Warning: the cut no longer have the same size so we // cannot use lp_values. Note that the substitution below shouldn't change @@ -805,7 +808,34 @@ void LinearProgrammingConstraint::AddCGCuts() { } } +namespace { + +// For each element of a, adds a random one in b and append the pair to output. +void RandomPick(const std::vector& a, const std::vector& b, + ModelRandomGenerator* random, + std::vector>* output) { + if (a.empty() || b.empty()) return; + for (const RowIndex row : a) { + const RowIndex other = b[absl::Uniform(*random, 0, b.size() - 1)]; + if (other != row) { + output->push_back({row, other}); + } + } +} + +template +IntegerValue GetCoeff(ColIndex col, const ListOfTerms& terms) { + for (const auto& term : terms) { + if (term.first == col) return term.second; + } + return IntegerValue(0); +} + +} // namespace + void LinearProgrammingConstraint::AddMirCuts() { + // We first try to derive MIR cuts for just one constraints. We call these + // MIR1 cuts. CHECK_EQ(trail_->CurrentDecisionLevel(), 0); const RowIndex num_rows = lp_data_.num_constraints(); for (RowIndex row(0); row < num_rows; ++row) { @@ -813,9 +843,6 @@ void LinearProgrammingConstraint::AddMirCuts() { if (status == glop::ConstraintStatus::BASIC) continue; if (status == glop::ConstraintStatus::FREE) continue; - // TODO(user): Do not consider just one constraint, but take linear - // combination of a small number of constraints. There is a lot of - // literature on the possible heuristics here. if (status == glop::ConstraintStatus::AT_UPPER_BOUND || status == glop::ConstraintStatus::FIXED_VALUE) { std::vector> integer_multipliers; @@ -829,6 +856,118 @@ void LinearProgrammingConstraint::AddMirCuts() { AddCutFromConstraints("MIR1", integer_multipliers); } } + + // Now, try linear combination of 2 different rows (MIR2). + // + // To limit the combinations we try. We first pick the top 50 variables, i.e. + // the ones with a LP relaxation farther from their bounds. And we try for + // each row to combine it with a random other row in order to eliminate the + // considered variable from the base constraint used to derive the cut. + // + // TODO(user): Maybe an heuristic is better than random? consider constraint + // sizes, small coefficients, number of BASIC variable in both constraints? + // + // TODO(user): generalize this to a linear combination of a bit more rows. + // The literature seems to go up to 6. The code architecture need to change a + // bit though because the usual heuristic is to incrementality try to add one + // more row to an existing linear combination of some rows. We also have to + // be more careful about integer overflow. + const int num_cols_to_consider = 50; + + // First pick the top BASIC variables to consider. + std::vector> col_candidates; + const int num_cols = integer_variables_.size(); + for (ColIndex col(0); col < num_cols; ++col) { + const int col_degree = lp_data_.GetSparseColumn(col).num_entries().value(); + if (col_degree <= 1) continue; + if (simplex_.GetVariableStatus(col) != glop::VariableStatus::BASIC) { + continue; + } + + const IntegerVariable var = integer_variables_[col.value()]; + const double lp_value = expanded_lp_solution_[var]; + const double lb = ToDouble(integer_trail_->LowerBound(var)); + const double ub = ToDouble(integer_trail_->UpperBound(var)); + const double bound_distance = std::min(ub - lp_value, lp_value - lb); + + col_candidates.push_back({bound_distance, col}); + } + + std::sort(col_candidates.begin(), col_candidates.end(), + std::greater>()); + if (col_candidates.size() > num_cols_to_consider) { + col_candidates.resize(num_cols_to_consider); + } + + // For each columns, we split the rows according to their type and coefficient + // sign. + std::vector ub_positive_rows; + std::vector ub_negative_rows; + std::vector lb_positive_rows; + std::vector lb_negative_rows; + std::vector> to_try; + + for (const auto& entry : col_candidates) { + const ColIndex col = entry.second; + ub_positive_rows.clear(); + ub_negative_rows.clear(); + lb_positive_rows.clear(); + lb_negative_rows.clear(); + for (const auto entry : lp_data_.GetSparseColumn(col)) { + const RowIndex row = entry.row(); + const auto status = simplex_.GetConstraintStatus(row); + if (status == glop::ConstraintStatus::BASIC) continue; + + // TODO(user): Instead of using FIXED_VALUE consider also both direction + // when we almost have an equality? that is if the LP constraints bounds + // are close from each others (<1e-6 ?). Initial experiments shows it + // doesn't change much, so I kept this version for now. Note that it might + // just be better to use the side that constrain the current lp optimal + // solution (that we get from the status). + if (status == glop::ConstraintStatus::FIXED_VALUE || + status == glop::ConstraintStatus::AT_UPPER_BOUND) { + if (entry.coefficient() > 0.0) { + ub_positive_rows.push_back(row); + } else { + ub_negative_rows.push_back(row); + } + } + if (status == glop::ConstraintStatus::FIXED_VALUE || + status == glop::ConstraintStatus::AT_LOWER_BOUND) { + if (entry.coefficient() > 0.0) { + lb_positive_rows.push_back(row); + } else { + lb_negative_rows.push_back(row); + } + } + } + + // We combine the row in such a way that -coeff2 * row1 + coeff1 * row2 is + // an upper bounded constraint that is also tight under the current LP + // solution. + to_try.clear(); + RandomPick(ub_positive_rows, ub_negative_rows, random_, &to_try); + RandomPick(lb_negative_rows, lb_positive_rows, random_, &to_try); + + // Note that here, fixed rows will appear on both side. We will not add + // anything if we randomly pick twice the same row. + RandomPick(lb_positive_rows, ub_positive_rows, random_, &to_try); + RandomPick(ub_negative_rows, lb_negative_rows, random_, &to_try); + gtl::STLSortAndRemoveDuplicates(&to_try); + + for (const auto& entry : to_try) { + const IntegerValue coeff1 = GetCoeff(col, integer_lp_[entry.first].terms); + const IntegerValue coeff2 = + GetCoeff(col, integer_lp_[entry.second].terms); + if (coeff1 == 0 || coeff2 == 0) continue; + const IntegerValue gcd = IntegerValue( + MathUtil::GCD64(std::abs(coeff1.value()), std::abs(coeff2.value()))); + std::vector> integer_multipliers; + integer_multipliers.push_back({entry.first, -coeff2 / gcd}); + integer_multipliers.push_back({entry.second, coeff1 / gcd}); + AddCutFromConstraints("MIR2", integer_multipliers); + } + } } void LinearProgrammingConstraint::UpdateSimplexIterationLimit( @@ -989,7 +1128,8 @@ bool LinearProgrammingConstraint::Propagate() { VLOG(2) << "LP objective [ " << ToDouble(propagated_lb) << ", " << ToDouble(integer_trail_->UpperBound(objective_cp_)) << " ] approx_lb += " - << ToDouble(approximate_new_lb - propagated_lb); + << ToDouble(approximate_new_lb - propagated_lb) << " gap: " + << integer_trail_->UpperBound(objective_cp_) - propagated_lb; } } else { FillReducedCostReasonIn(simplex_.GetReducedCosts(), &integer_reason_); @@ -1341,7 +1481,10 @@ LinearProgrammingConstraint::ScaleLpMultiplier( } // We want max_sum * scaling to be <= 2 ^ max_pow and fit on an int64. - *scaling = std::ldexp(1, max_pow) / max_sum; + // We use a power of 2 as this seems to work better. + const double threshold = std::ldexp(1, max_pow) / max_sum; + *scaling = 1.0; + while (2 * *scaling <= threshold) *scaling *= 2; // Scale the multipliers by *scaling. // diff --git a/ortools/sat/linear_programming_constraint.h b/ortools/sat/linear_programming_constraint.h index 13dc97846d..5feddbc354 100644 --- a/ortools/sat/linear_programming_constraint.h +++ b/ortools/sat/linear_programming_constraint.h @@ -358,6 +358,7 @@ class LinearProgrammingConstraint : public PropagatorInterface, Trail* trail_; SearchHeuristicsVector* model_heuristics_; IntegerEncoder* integer_encoder_; + ModelRandomGenerator* random_; // Used while deriving cuts. ImpliedBoundsProcessor implied_bounds_processor_; diff --git a/ortools/sat/python/cp_model.py b/ortools/sat/python/cp_model.py index 2aeb1ba3f0..ba7b800725 100644 --- a/ortools/sat/python/cp_model.py +++ b/ortools/sat/python/cp_model.py @@ -157,6 +157,7 @@ class LinearExpr(object): model.Minimize(cp_model.LinearExpr.Sum(expressions)) model.Add(cp_model.LinearExpr.ScalProd(expressions, coefficients) >= 0) """ + @classmethod def Sum(cls, expressions): """Creates the expression sum(expressions).""" @@ -320,6 +321,7 @@ class LinearExpr(object): class _ProductCst(LinearExpr): """Represents the product of a LinearExpr by a constant.""" + def __init__(self, expr, coef): cp_model_helper.AssertIsInt64(coef) if isinstance(expr, _ProductCst): @@ -348,6 +350,7 @@ class _ProductCst(LinearExpr): class _SumArray(LinearExpr): """Represents the sum of a list of LinearExpr and a constant.""" + def __init__(self, expressions): self.__expressions = [] self.__constant = 0 @@ -380,6 +383,7 @@ class _SumArray(LinearExpr): class _ScalProd(LinearExpr): """Represents the scalar product of expressions with constants and a constant.""" + def __init__(self, expressions, coefficients): self.__expressions = [] self.__coefficients = [] @@ -452,6 +456,7 @@ class IntVar(LinearExpr): from the set of initial values (called the initial domain), such that the model is feasible, or optimal if you provided an objective function. """ + def __init__(self, model, domain, name): """See CpModel.NewIntVar below.""" self.__model = model @@ -472,7 +477,7 @@ class IntVar(LinearExpr): def __str__(self): if not self.__var.name: if len(self.__var.domain - ) == 2 and self.__var.domain[0] == self.__var.domain[1]: + ) == 2 and self.__var.domain[0] == self.__var.domain[1]: # Special case for constants. return str(self.__var.domain[0]) else: @@ -505,6 +510,7 @@ class IntVar(LinearExpr): class _NotBooleanVariable(LinearExpr): """Negation of a boolean variable.""" + def __init__(self, boolvar): self.__boolvar = boolvar @@ -526,6 +532,7 @@ class BoundedLinearExpression(object): model.Add(x + 2 * y -1 >= z) """ + def __init__(self, expr, bounds): self.__expr = expr self.__bounds = bounds @@ -571,6 +578,7 @@ class Constraint(object): model.Add(x + 2 * y == 5).OnlyEnforceIf(b.Not()) """ + def __init__(self, constraints): self.__index = len(constraints) self.__constraint = constraints.add() @@ -632,6 +640,7 @@ class IntervalVar(object): also set these enforcement literals to false if they cannot fit these intervals into the schedule. """ + def __init__(self, model, start_index, size_index, end_index, is_present_index, name): self.__model = model @@ -661,14 +670,14 @@ class IntervalVar(object): if self.__ct.enforcement_literal: return '%s(start = %s, size = %s, end = %s, is_present = %s)' % ( self.__ct.name, ShortName(self.__model, interval.start), - ShortName(self.__model, interval.size), - ShortName(self.__model, interval.end), + ShortName(self.__model, + interval.size), ShortName(self.__model, interval.end), ShortName(self.__model, self.__ct.enforcement_literal[0])) else: return '%s(start = %s, size = %s, end = %s)' % ( self.__ct.name, ShortName(self.__model, interval.start), - ShortName(self.__model, interval.size), - ShortName(self.__model, interval.end)) + ShortName(self.__model, + interval.size), ShortName(self.__model, interval.end)) def Name(self): return self.__ct.name @@ -682,6 +691,7 @@ class CpModel(object): * ```New``` create integer, boolean, or interval variables. * ```Add``` create new constraints and add them to the model. """ + def __init__(self): self.__model = cp_model_pb2.CpModelProto() self.__constant_map = {} @@ -1064,8 +1074,7 @@ class CpModel(object): ct = Constraint(self.__model.constraints) model_ct = self.__model.constraints[ct.Index()] - model_ct.reservoir.times.extend( - [self.GetOrMakeIndex(x) for x in times]) + model_ct.reservoir.times.extend([self.GetOrMakeIndex(x) for x in times]) model_ct.reservoir.demands.extend(demands) model_ct.reservoir.min_level = min_level model_ct.reservoir.max_level = max_level @@ -1117,8 +1126,7 @@ class CpModel(object): ct = Constraint(self.__model.constraints) model_ct = self.__model.constraints[ct.Index()] - model_ct.reservoir.times.extend( - [self.GetOrMakeIndex(x) for x in times]) + model_ct.reservoir.times.extend([self.GetOrMakeIndex(x) for x in times]) model_ct.reservoir.demands.extend(demands) model_ct.reservoir.actives.extend( [self.GetOrMakeIndex(x) for x in actives]) @@ -1383,9 +1391,8 @@ class CpModel(object): """Returns the index of a variable, its negation, or a number.""" if isinstance(arg, IntVar): return arg.Index() - elif (isinstance(arg, _ProductCst) - and isinstance(arg.Expression(), IntVar) - and arg.Coefficient() == -1): + elif (isinstance(arg, _ProductCst) and + isinstance(arg.Expression(), IntVar) and arg.Coefficient() == -1): return -arg.Expression().Index() - 1 elif isinstance(arg, numbers.Integral): cp_model_helper.AssertIsInt64(arg) @@ -1569,6 +1576,7 @@ class CpSolver(object): with the Value() and BooleanValue() methods, as well as general statistics about the solve procedure. """ + def __init__(self): self.__model = None self.__solution = None @@ -1695,6 +1703,7 @@ class CpSolverSolutionCallback(pywrapsat.SolutionCallback): These methods returns the same information as their counterpart in the `CpSolver` class. """ + def __init__(self): pywrapsat.SolutionCallback.__init__(self) @@ -1770,6 +1779,7 @@ class CpSolverSolutionCallback(pywrapsat.SolutionCallback): class ObjectiveSolutionPrinter(CpSolverSolutionCallback): """Display the objective value and time of intermediate solutions.""" + def __init__(self): CpSolverSolutionCallback.__init__(self) self.__solution_count = 0 @@ -1790,6 +1800,7 @@ class ObjectiveSolutionPrinter(CpSolverSolutionCallback): class VarArrayAndObjectiveSolutionPrinter(CpSolverSolutionCallback): """Print intermediate solutions (objective, variable values, time).""" + def __init__(self, variables): CpSolverSolutionCallback.__init__(self) self.__variables = variables @@ -1814,6 +1825,7 @@ class VarArrayAndObjectiveSolutionPrinter(CpSolverSolutionCallback): class VarArraySolutionPrinter(CpSolverSolutionCallback): """Print intermediate solutions (variable values, time).""" + def __init__(self, variables): CpSolverSolutionCallback.__init__(self) self.__variables = variables diff --git a/ortools/sat/samples/default.profraw b/ortools/sat/samples/default.profraw deleted file mode 100644 index b679fd1088a0543d3e23d58b3657ce073a651cdd..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 2056 zcmZoHO3N=Q$obF000E*l&*2L;iMANdub_1 zFaCen1?JxgQ1uF|5PdLln16Y;+dXhO_I~=+xBp%?%n#R`0af3CrXCiavudB)vCQO{ z3o~a2RQ&@s?B-vvle$yiWW6^~MRVd34=#oSQ1uJgQPsoD7ZkXA(f3n<^WK2yM@`S1 z8XrK_GjO1)7lt}O#`o|+gRAC&5>S4HIYiwBZdCO!|8~!4>S}%7c^am^0jmB34)qc) zm7BLZPW}K>-vL!Gz=PdjGT_21O(+A_jXeKBFCXddCiKEjn{pd7I9lH2X;iJ1B#)r``f1vYW;-u0r zbI@s+IWTb;pAe0%9>zzf$u$S24`v>WPl$%4D?;j)F`r?qxHYx!v>&q}kLy(VRRW^Y zWsCJRy1P~|zt&hXeP+fai%oYY{-|EIbY)0H>GSHEKO63zGkUs8p2=+6?Y+uN!=leO z^|4#Fif62sIsIYP3GodJjkfPNVds1>eRumjY4&AVmgkawL~QdA-E)Uqz_9Vv_tpy8 zUi()aiwrA{zGt?}@5sKtf9mdclY^v`7`3oca`2MPVG7n vGJo!q*$0oYNJ+G>7HZmMlakw@bxNtzMXR^frK!HaW^MSaI)*zJ=cxezDySf8 diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index bd88c284d7..eebafd43e6 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -541,7 +541,7 @@ message SatParameters { // // TODO(user): We should probably remove this parameters, and just always // generate cuts but only keep the best n or something. - optional int32 max_num_cuts = 91 [default = 5000]; + optional int32 max_num_cuts = 91 [default = 10000]; // For the cut that can be generated at any level, this control if we only // try to generate them at the root node. diff --git a/ortools/sat/swig_helper.h b/ortools/sat/swig_helper.h index 9b61d71c40..526607ea53 100644 --- a/ortools/sat/swig_helper.h +++ b/ortools/sat/swig_helper.h @@ -76,6 +76,7 @@ class SolutionCallback { void StopSearch() { stopped_ = true; } + // Reset the shared atomic Boolean used to stop the search. void ResetSharedBoolean() const { stopped_ = false; } std::atomic* stopped() const { return &stopped_; }