always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

internal improvemensts

Browse Source
This commit is contained in:
Laurent Perron
2019-12-05 16:36:11 +01:00
parent fd04db1add
commit dc3d9ccf84
24 changed files with 753 additions and 227 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
ortools/base/random.h
View File

@@ -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) {}

11
ortools/glop/lp_solver.cc
View File

@@ -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;
}

3
ortools/glop/preprocessor.h
View File

@@ -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.
//

2
ortools/graph/perfect_matching.h
View File

@@ -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;

4
ortools/linear_solver/cbc_interface.cc
View File

@@ -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;

2
ortools/linear_solver/gurobi_proto_solver.cc
View File

@@ -294,7 +294,7 @@ util::StatusOr<MPSolutionResponse> 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;
}
}

4
ortools/linear_solver/linear_solver.h
View File

@@ -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.

6
ortools/linear_solver/samples/SimpleMipProgram.cs
View File

@@ -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());

16
ortools/linear_solver/samples/mip_var_array.cc
View File

@@ -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();

2
ortools/linear_solver/scip_proto_solver.cc
View File

@@ -773,7 +773,7 @@ util::StatusOr<MPSolutionResponse> 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

4
ortools/sat/cp_model_expand.cc
View File

@@ -657,8 +657,8 @@ void ExpandAutomaton(ConstraintProto* ct, PresolveContext* context) {
const std::vector<int> vars = {proto.vars().begin(), proto.vars().end()};
// Compute the set of reachable state at each time point.
const absl::flat_hash_set<int64> final_states({proto.final_states().begin(),
proto.final_states().end()});
const absl::flat_hash_set<int64> final_states(
{proto.final_states().begin(), proto.final_states().end()});
std::vector<absl::flat_hash_set<int64>> reachable_states(n + 1);
reachable_states[0].insert(proto.starting_state());

171
ortools/sat/cp_model_lns.cc
View File

@@ -14,6 +14,7 @@
#include "ortools/sat/cp_model_lns.h"
#include <numeric>
#include <random>
#include <vector>
#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<int>* base,
Neighborhood SimpleNeighborhoodGenerator::Generate(
const CpSolverResponse& initial_solution, double difficulty,
random_engine_t* random) const {
random_engine_t* random) {
std::vector<int> 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<int> 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<std::pair<int64, int>> 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<int> removed_constraints;
random_engine_t* random) {
std::vector<int> 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<int> random_var(
0, removable_constraints.size() - 1);
const int random_start_index = random_var(*random);
std::vector<int> 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<int> 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<int>& 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<int> 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<std::pair<double, int>> constraint_removal_scores;
std::uniform_real_distribution<double> 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

100
ortools/sat/cp_model_lns.h
View File

@@ -16,6 +16,7 @@
#include <vector>
#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<SolveData> 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<double> constraint_weights_;
int num_removable_constraints_ = 0;
// Indices of the removed constraints per generated neighborhood.
absl::flat_hash_map<int64, std::vector<int>> 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

113
ortools/sat/cp_model_solver.cc
View File

@@ -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<int64>(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<int64>(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<int64>(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<LnsSolver>(
/*id=*/subsolvers.size(),
absl::make_unique<RandomRelaxationNeighborhoodGenerator>(
absl::make_unique<
ConsecutiveConstraintsRelaxationNeighborhoodGenerator>(
helper, absl::StrCat("rnd_rel_lns_", strategy_name)),
local_params, helper, &shared));
subsolvers.push_back(absl::make_unique<LnsSolver>(
/*id=*/subsolvers.size(),
absl::make_unique<WeightedRandomRelaxationNeighborhoodGenerator>(
helper, absl::StrCat("wgt_rel_lns_", strategy_name)),
local_params, helper, &shared));
}
if (!helper->TypeToConstraints(ConstraintProto::kNoOverlap).empty()) {

318
ortools/sat/cuts.cc
View File

@@ -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<IntegerValue(IntegerValue)> 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<IntegerValue(IntegerValue)> 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<double> lp_values,
std::vector<IntegerValue> lower_bounds,
std::vector<IntegerValue> 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<double>& lp_values,
const std::vector<IntegerValue>& lower_bounds,
const std::vector<IntegerValue>& upper_bounds,
LinearConstraint* cut) {
const int size = lp_values.size();
if (size == 0) return;
@@ -671,87 +675,151 @@ void IntegerRoundingCut(RoundingOptions options, std::vector<double> 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<IntegerValue>& 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<int> relevant_indices;
std::vector<double> relevant_lp_values;
std::vector<IntegerValue> relevant_coeffs;
std::vector<IntegerValue> relevant_bound_diffs;
std::vector<IntegerValue> divisors;
std::vector<std::pair<int, IntegerValue>> 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<bool> 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<IntegerValue>());
// 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<double> 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;

7
ortools/sat/cuts.h
View File

@@ -145,9 +145,10 @@ struct RoundingOptions {
bool use_mir = false;
IntegerValue max_scaling = IntegerValue(60);
};
void IntegerRoundingCut(RoundingOptions options, std::vector<double> lp_values,
std::vector<IntegerValue> lower_bounds,
std::vector<IntegerValue> upper_bounds,
void IntegerRoundingCut(RoundingOptions options,
const std::vector<double>& lp_values,
const std::vector<IntegerValue>& lower_bounds,
const std::vector<IntegerValue>& upper_bounds,
LinearConstraint* cut);
// If a variable is away from its upper bound by more than value 1.0, then it

4
ortools/sat/integer.h
View File

@@ -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<IntegerValue>(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<IntegerValue>(result * positive_divisor > dividend);

6
ortools/sat/linear_constraint_manager.cc
View File

@@ -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;
}

165
ortools/sat/linear_programming_constraint.cc
View File

@@ -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<Trail>()),
model_heuristics_(model->GetOrCreate<SearchHeuristicsVector>()),
integer_encoder_(model->GetOrCreate<IntegerEncoder>()),
random_(model->GetOrCreate<ModelRandomGenerator>()),
implied_bounds_processor_({}, integer_trail_,
model->GetOrCreate<ImpliedBounds>()),
dispatcher_(model->GetOrCreate<LinearProgrammingDispatcher>()),
@@ -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<RowIndex>& a, const std::vector<RowIndex>& b,
ModelRandomGenerator* random,
std::vector<std::pair<RowIndex, RowIndex>>* output) {
if (a.empty() || b.empty()) return;
for (const RowIndex row : a) {
const RowIndex other = b[absl::Uniform<int>(*random, 0, b.size() - 1)];
if (other != row) {
output->push_back({row, other});
}
}
}
template <class ListOfTerms>
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<std::pair<RowIndex, IntegerValue>> 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<std::pair<double, ColIndex>> 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<std::pair<double, ColIndex>>());
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<RowIndex> ub_positive_rows;
std::vector<RowIndex> ub_negative_rows;
std::vector<RowIndex> lb_positive_rows;
std::vector<RowIndex> lb_negative_rows;
std::vector<std::pair<RowIndex, RowIndex>> 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<std::pair<RowIndex, IntegerValue>> 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.
//

1
ortools/sat/linear_programming_constraint.h
View File

@@ -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_;

36
ortools/sat/python/cp_model.py
View File

@@ -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

BIN
ortools/sat/samples/default.profraw
View File

Binary file not shown.

2
ortools/sat/sat_parameters.proto
View File

@@ -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.

1
ortools/sat/swig_helper.h
View File

@@ -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<bool>* stopped() const { return &stopped_; }