remove solow-halim glop preprocessor as the interface has changed a lot; add a proper domain class for CP-SAT integer variables, rewrite preprocessor and other parts of the solver to use it

This commit is contained in:
Laurent Perron
2018-10-02 10:35:52 +02:00
parent bb86602f7c
commit 4e893e5080
21 changed files with 857 additions and 759 deletions

View File

@@ -270,10 +270,10 @@ void CpModelProtoWithMapping::FillConstraint(const fz::Constraint& fz_ct,
arg->add_vars(LookupVar(fz_ct.arguments[0]));
arg->add_coeffs(1);
if (fz_ct.arguments[1].type == fz::Argument::INT_LIST) {
FillDomain(
SortedDisjointIntervalsFromValues({fz_ct.arguments[1].values.begin(),
fz_ct.arguments[1].values.end()}),
arg);
FillDomainInProto(Domain::FromValues(std::vector<int64>{
fz_ct.arguments[1].values.begin(),
fz_ct.arguments[1].values.end()}),
arg);
} else if (fz_ct.arguments[1].type == fz::Argument::INT_INTERVAL) {
FillDomain({{fz_ct.arguments[1].values[0], fz_ct.arguments[1].values[1]}},
arg);
@@ -285,15 +285,16 @@ void CpModelProtoWithMapping::FillConstraint(const fz::Constraint& fz_ct,
arg->add_vars(LookupVar(fz_ct.arguments[0]));
arg->add_coeffs(1);
if (fz_ct.arguments[1].type == fz::Argument::INT_LIST) {
FillDomain(
ComplementOfSortedDisjointIntervals(SortedDisjointIntervalsFromValues(
{fz_ct.arguments[1].values.begin(),
fz_ct.arguments[1].values.end()})),
FillDomainInProto(
Domain::FromValues(
std::vector<int64>{fz_ct.arguments[1].values.begin(),
fz_ct.arguments[1].values.end()})
.Complement(),
arg);
} else if (fz_ct.arguments[1].type == fz::Argument::INT_INTERVAL) {
FillDomain(
ComplementOfSortedDisjointIntervals(
{{fz_ct.arguments[1].values[0], fz_ct.arguments[1].values[1]}}),
FillDomainInProto(
Domain(fz_ct.arguments[1].values[0], fz_ct.arguments[1].values[1])
.Complement(),
arg);
} else {
LOG(FATAL) << "Wrong format";
@@ -456,21 +457,18 @@ void CpModelProtoWithMapping::FillConstraint(const fz::Constraint& fz_ct,
int index = min_index;
const bool is_circuit = (fz_ct.type == "circuit");
for (const int var : LookupVars(fz_ct.arguments[0])) {
Domain domain = ReadDomainFromProto(proto.variables(var));
// Restrict the domain of var to [min_index, max_index]
FillDomain(
IntersectionOfSortedDisjointIntervals(
ReadDomain(proto.variables(var)), {{min_index, max_index}}),
proto.mutable_variables(var));
domain = domain.IntersectionWith(Domain(min_index, max_index));
if (is_circuit) {
// We simply make sure that the variable cannot take the value index.
FillDomain(IntersectionOfSortedDisjointIntervals(
ReadDomain(proto.variables(var)),
{{kint64min, index - 1}, {index + 1, kint64max}}),
proto.mutable_variables(var));
domain = domain.IntersectionWith(Domain::FromIntervals(
{{kint64min, index - 1}, {index + 1, kint64max}}));
}
FillDomainInProto(domain, proto.mutable_variables(var));
const auto domain = ReadDomain(proto.variables(var));
for (const auto interval : domain) {
for (const ClosedInterval interval : domain.intervals()) {
for (int64 value = interval.start; value <= interval.end; ++value) {
// Create one Boolean variable for this arc.
const int literal = proto.variables_size();
@@ -540,20 +538,20 @@ void CpModelProtoWithMapping::FillConstraint(const fz::Constraint& fz_ct,
for (const int var : direct_variables) {
arg->add_f_direct(var);
// Intersect domains with offset + [0, num_variables).
FillDomain(IntersectionOfSortedDisjointIntervals(
ReadDomain(proto.variables(var)),
{{offset, num_variables - 1 + offset}}),
proto.mutable_variables(var));
FillDomainInProto(
ReadDomainFromProto(proto.variables(var))
.IntersectionWith(Domain(offset, num_variables - 1 + offset)),
proto.mutable_variables(var));
}
if (is_one_based) arg->add_f_inverse(LookupConstant(0));
for (const int var : inverse_variables) {
arg->add_f_inverse(var);
// Intersect domains with offset + [0, num_variables).
FillDomain(IntersectionOfSortedDisjointIntervals(
ReadDomain(proto.variables(var)),
{{offset, num_variables - 1 + offset}}),
proto.mutable_variables(var));
FillDomainInProto(
ReadDomainFromProto(proto.variables(var))
.IntersectionWith(Domain(offset, num_variables - 1 + offset)),
proto.mutable_variables(var));
}
} else if (fz_ct.type == "cumulative") {
const std::vector<int> starts = LookupVars(fz_ct.arguments[0]);
@@ -863,12 +861,7 @@ void SolveFzWithCpModelProto(const fz::Model& fz_model,
var->add_domain(fz_var->domain.values[1]);
}
} else {
const std::vector<ClosedInterval> domain =
SortedDisjointIntervalsFromValues(fz_var->domain.values);
for (const ClosedInterval& interval : domain) {
var->add_domain(interval.start);
var->add_domain(interval.end);
}
FillDomainInProto(Domain::FromValues(fz_var->domain.values), var);
}
// Some variables in flatzinc have large domain and we don't really support
@@ -877,9 +870,9 @@ void SolveFzWithCpModelProto(const fz::Model& fz_model,
// variable domains with [kint32min, kint32max].
//
// TODO(user): Warn when we reduce the domain.
FillDomain(IntersectionOfSortedDisjointIntervals(ReadDomain(*var),
{{kint32min, kint32max}}),
var);
FillDomainInProto(ReadDomainFromProto(*var).IntersectionWith(
Domain(kint32min, kint32max)),
var);
}
// Translate the constraints.

View File

@@ -174,8 +174,7 @@ ProblemStatus LPSolver::SolveWithTimeLimit(const LinearProgram& lp,
current_linear_program_.PopulateFromLinearProgram(lp);
// Preprocess.
MainLpPreprocessor preprocessor;
preprocessor.SetParameters(parameters_);
MainLpPreprocessor preprocessor(&parameters_);
preprocessor.SetTimeLimit(time_limit);
const bool postsolve_is_needed = preprocessor.Run(&current_linear_program_);

View File

@@ -377,15 +377,6 @@ message GlopParameters {
// not create any OMP threads and will remain single-threaded.
optional int32 num_omp_threads = 44 [default = 1];
// Whether or not use Solow-Halim preprocessor. It tries to begin the phase 1
// with a closer point to the optimal solution.
//
// Reference: Improving the Efficiency of the Simplex Algorithm Based on a
// Geometric Explanation of Phase 1
// https://weatherhead.case.edu/faculty/research/library/detail?id=12128593921
// DOI: 10.1504/IJOR.2009.025701
optional bool use_solowhalim = 46 [default = false];
// When this is true, then the costs are randomly perturbed before the dual
// simplex is even started. This has been shown to improve the dual simplex
// performance. For a good reference, see Huangfu Q (2013) "High performance

View File

@@ -17,6 +17,7 @@
#include "ortools/glop/revised_simplex.h"
#include "ortools/glop/status.h"
#include "ortools/lp_data/lp_data_utils.h"
#include "ortools/lp_data/lp_types.h"
#include "ortools/lp_data/lp_utils.h"
#include "ortools/lp_data/matrix_utils.h"
@@ -39,9 +40,9 @@ double trunc(double d) { return d > 0 ? floor(d) : ceil(d); }
// --------------------------------------------------------
// Preprocessor
// --------------------------------------------------------
Preprocessor::Preprocessor()
Preprocessor::Preprocessor(const GlopParameters* parameters)
: status_(ProblemStatus::INIT),
parameters_(),
parameters_(*parameters),
in_mip_context_(false),
time_limit_(TimeLimit::Infinite().get()) {}
Preprocessor::~Preprocessor() {}
@@ -50,9 +51,9 @@ Preprocessor::~Preprocessor() {}
// MainLpPreprocessor
// --------------------------------------------------------
#define RUN_PREPROCESSOR(name) \
RunAndPushIfRelevant(std::unique_ptr<Preprocessor>(new name()), #name, \
time_limit_, lp)
#define RUN_PREPROCESSOR(name) \
RunAndPushIfRelevant(std::unique_ptr<Preprocessor>(new name(&parameters_)), \
#name, time_limit_, lp)
bool MainLpPreprocessor::Run(LinearProgram* lp) {
RETURN_VALUE_IF_NULL(lp, false);
@@ -101,7 +102,6 @@ bool MainLpPreprocessor::Run(LinearProgram* lp) {
// preprocessor so that the rhs/objective of the dual are with a good
// magnitude.
RUN_PREPROCESSOR(DualizerPreprocessor);
RUN_PREPROCESSOR(SolowHalimPreprocessor);
if (old_stack_size != preprocessors_.size()) {
RUN_PREPROCESSOR(SingletonPreprocessor);
RUN_PREPROCESSOR(FreeConstraintPreprocessor);
@@ -131,7 +131,6 @@ void MainLpPreprocessor::RunAndPushIfRelevant(
if (status_ != ProblemStatus::INIT || time_limit->LimitReached()) return;
const double start_time = time_limit->GetElapsedTime();
preprocessor->SetParameters(parameters_);
preprocessor->SetTimeLimit(time_limit);
// No need to run the preprocessor if the lp is empty.
@@ -145,7 +144,7 @@ void MainLpPreprocessor::RunAndPushIfRelevant(
const EntryIndex new_num_entries = lp->num_entries();
const double preprocess_time = time_limit->GetElapsedTime() - start_time;
VLOG(1) << absl::StrFormat(
"%s(%fs): %d(%d) rows, %d(%d) columns, %d(%d) entries.", name.c_str(),
"%s(%fs): %d(%d) rows, %d(%d) columns, %d(%d) entries.", name,
preprocess_time, lp->num_constraints().value(),
(lp->num_constraints() - initial_num_rows_).value(),
lp->num_variables().value(),
@@ -3625,149 +3624,5 @@ void AddSlackVariablesPreprocessor::RecoverSolution(
solution->variable_statuses.resize(first_slack_col_, VariableStatus::FREE);
}
// --------------------------------------------------------
// SolowHalimPreprocessor
// --------------------------------------------------------
bool SolowHalimPreprocessor::Run(LinearProgram* lp) {
RETURN_VALUE_IF_NULL(lp, false);
if (!parameters_.use_solowhalim()) {
return false;
}
// mip context not implemented
// TODO : in order to manage mip context we must take care
// of truncated offsets
if (in_mip_context_) {
return false;
}
const bool lp_is_maximization_problem = lp->IsMaximizationProblem();
const ColIndex num_cols = lp->num_variables();
const RowIndex num_rows = lp->num_constraints();
ColIndex num_shifted_cols(0);
ColIndex num_shifted_opposite_cols(0);
variable_initial_lbs_.assign(num_cols, 0.0);
variable_initial_ubs_.assign(num_cols, 0.0);
for (ColIndex col(0); col < num_cols; ++col) {
variable_initial_lbs_[col] = lp->variable_lower_bounds()[col];
variable_initial_ubs_[col] = lp->variable_upper_bounds()[col];
}
KahanSum objective_offset;
gtl::ITIVector<RowIndex, KahanSum> row_offsets(num_rows.value());
column_transform_.resize(num_cols.value(), NOT_MODIFIED);
for (ColIndex col(0); col < num_cols; ++col) {
const Fractional coeff = lp->objective_coefficients()[col];
if (coeff != 0.0) {
bool coeff_opposite_direction =
((coeff < 0.0 && !lp_is_maximization_problem) ||
(coeff > 0.0 && lp_is_maximization_problem));
const Fractional ub = lp->variable_upper_bounds()[col];
const Fractional lb = lp->variable_lower_bounds()[col];
if (IsFinite(ub) && IsFinite(lb)) {
ColumnTransformType column_transform = NOT_MODIFIED;
double shift_value = 0.0;
if (coeff_opposite_direction) {
SparseColumn* mutable_sparse_column = lp->GetMutableSparseColumn(col);
for (const SparseColumn::Entry e : (*mutable_sparse_column)) {
row_offsets[e.row()].Add(e.coefficient() * ub);
}
mutable_sparse_column->MultiplyByConstant(-1);
lp->SetObjectiveCoefficient(col, -coeff);
shift_value = ub;
column_transform = SHIFTED_OPPOSITE_DIRECTION;
num_shifted_opposite_cols++;
} else {
const SparseColumn& sparse_column = lp->GetSparseColumn(col);
for (const SparseColumn::Entry e : sparse_column) {
row_offsets[e.row()].Add(e.coefficient() * lb);
}
shift_value = lb;
column_transform = SHIFTED;
num_shifted_cols++;
}
if (column_transform != NOT_MODIFIED) {
column_transform_[col] = column_transform;
objective_offset.Add(coeff * shift_value);
lp->SetVariableBounds(col, 0, ub - lb);
}
}
}
}
for (RowIndex row(0); row < num_rows; ++row) {
lp->SetConstraintBounds(
row, lp->constraint_lower_bounds()[row] - row_offsets[row].Value(),
lp->constraint_upper_bounds()[row] - row_offsets[row].Value());
}
lp->SetObjectiveOffset(lp->objective_offset() + objective_offset.Value());
VLOG(1) << "Shifted " << num_shifted_cols << " variables.";
VLOG(1) << "Shifted opposite " << num_shifted_opposite_cols << " variables.";
VLOG(1) << "Objective offset : " << objective_offset.Value();
return true;
}
void SolowHalimPreprocessor::RecoverSolution(ProblemSolution* solution) const {
RETURN_IF_NULL(solution);
const ColIndex num_cols = solution->variable_statuses.size();
for (ColIndex col(0); col < num_cols; ++col) {
VLOG(2) << "col = " << col << "\t" << column_transform_[col];
VLOG(2) << "\tinitial range : \t [" << variable_initial_lbs_[col] << " ; "
<< variable_initial_ubs_[col] << "]";
VLOG(2) << "\tstatus : " << solution->variable_statuses[col]
<< "\t raw value : " << solution->primal_values[col];
switch (column_transform_[col]) {
case NOT_MODIFIED:
break;
case SHIFTED:
switch (solution->variable_statuses[col]) {
case VariableStatus::AT_LOWER_BOUND:
solution->primal_values[col] = variable_initial_lbs_[col];
break;
case VariableStatus::AT_UPPER_BOUND:
solution->primal_values[col] = variable_initial_ubs_[col];
break;
case VariableStatus::BASIC:
solution->primal_values[col] =
variable_initial_lbs_[col] + solution->primal_values[col];
break;
case VariableStatus::FIXED_VALUE:
FALLTHROUGH_INTENDED;
case VariableStatus::FREE:
break;
}
break;
case SHIFTED_OPPOSITE_DIRECTION:
switch (solution->variable_statuses[col]) {
case VariableStatus::AT_LOWER_BOUND:
solution->primal_values[col] = variable_initial_ubs_[col];
solution->variable_statuses[col] = VariableStatus::AT_UPPER_BOUND;
break;
case VariableStatus::AT_UPPER_BOUND:
solution->primal_values[col] = variable_initial_lbs_[col];
solution->variable_statuses[col] = VariableStatus::AT_LOWER_BOUND;
break;
case VariableStatus::BASIC:
solution->primal_values[col] =
variable_initial_ubs_[col] - solution->primal_values[col];
break;
case VariableStatus::FIXED_VALUE:
FALLTHROUGH_INTENDED;
case VariableStatus::FREE:
break;
}
break;
}
VLOG(2) << " recover value : " << solution->primal_values[col];
}
}
} // namespace glop
} // namespace operations_research

View File

@@ -41,7 +41,9 @@ namespace glop {
// as expected. Fix? or document and crash in debug if this happens.
class Preprocessor {
public:
Preprocessor();
explicit Preprocessor(const GlopParameters* parameters);
Preprocessor(const Preprocessor&) = delete;
Preprocessor& operator=(const Preprocessor&) = delete;
virtual ~Preprocessor();
// Runs the preprocessor by modifying the given linear program. Returns true
@@ -61,11 +63,6 @@ class Preprocessor {
// solved and there is not need to call subsequent preprocessors.
ProblemStatus status() const { return status_; }
// Stores the parameters for use by the different preprocessors.
void SetParameters(const GlopParameters& parameters) {
parameters_ = parameters;
}
// Some preprocessors only need minimal changes when used with integer
// variables in a MIP context. Setting this to true allows to consider integer
// variables as integer in these preprocessors.
@@ -91,12 +88,9 @@ class Preprocessor {
}
ProblemStatus status_;
GlopParameters parameters_;
const GlopParameters& parameters_;
bool in_mip_context_;
TimeLimit* time_limit_;
private:
DISALLOW_COPY_AND_ASSIGN(Preprocessor);
};
// --------------------------------------------------------
@@ -106,7 +100,10 @@ class Preprocessor {
// preprocessors in this file, possibly more than once.
class MainLpPreprocessor : public Preprocessor {
public:
MainLpPreprocessor() {}
explicit MainLpPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
MainLpPreprocessor(const MainLpPreprocessor&) = delete;
MainLpPreprocessor& operator=(const MainLpPreprocessor&) = delete;
~MainLpPreprocessor() override {}
bool Run(LinearProgram* lp) final;
@@ -130,8 +127,6 @@ class MainLpPreprocessor : public Preprocessor {
EntryIndex initial_num_entries_;
RowIndex initial_num_rows_;
ColIndex initial_num_cols_;
DISALLOW_COPY_AND_ASSIGN(MainLpPreprocessor);
};
// --------------------------------------------------------
@@ -141,6 +136,8 @@ class MainLpPreprocessor : public Preprocessor {
class ColumnDeletionHelper {
public:
ColumnDeletionHelper() {}
ColumnDeletionHelper(const ColumnDeletionHelper&) = delete;
ColumnDeletionHelper& operator=(const ColumnDeletionHelper&) = delete;
// Remember the given column as "deleted" so that it can later be restored
// by RestoreDeletedColumns(). Optionally, the caller may indicate the
@@ -186,7 +183,6 @@ class ColumnDeletionHelper {
// data structure so columns can be deleted in any order if needed.
DenseRow stored_value_;
VariableStatusRow stored_status_;
DISALLOW_COPY_AND_ASSIGN(ColumnDeletionHelper);
};
// --------------------------------------------------------
@@ -196,6 +192,8 @@ class ColumnDeletionHelper {
class RowDeletionHelper {
public:
RowDeletionHelper() {}
RowDeletionHelper(const RowDeletionHelper&) = delete;
RowDeletionHelper& operator=(const RowDeletionHelper&) = delete;
// Returns true if no rows have been marked for deletion.
bool IsEmpty() const { return is_row_deleted_.empty(); }
@@ -224,8 +222,6 @@ class RowDeletionHelper {
private:
DenseBooleanColumn is_row_deleted_;
DISALLOW_COPY_AND_ASSIGN(RowDeletionHelper);
};
// --------------------------------------------------------
@@ -234,14 +230,16 @@ class RowDeletionHelper {
// Removes the empty columns from the problem.
class EmptyColumnPreprocessor : public Preprocessor {
public:
EmptyColumnPreprocessor() {}
explicit EmptyColumnPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
EmptyColumnPreprocessor(const EmptyColumnPreprocessor&) = delete;
EmptyColumnPreprocessor& operator=(const EmptyColumnPreprocessor&) = delete;
~EmptyColumnPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
private:
ColumnDeletionHelper column_deletion_helper_;
DISALLOW_COPY_AND_ASSIGN(EmptyColumnPreprocessor);
};
// --------------------------------------------------------
@@ -259,7 +257,12 @@ class EmptyColumnPreprocessor : public Preprocessor {
// so it makes sense to use the more general notion of proportional columns.
class ProportionalColumnPreprocessor : public Preprocessor {
public:
ProportionalColumnPreprocessor() {}
explicit ProportionalColumnPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
ProportionalColumnPreprocessor(const ProportionalColumnPreprocessor&) =
delete;
ProportionalColumnPreprocessor& operator=(
const ProportionalColumnPreprocessor&) = delete;
~ProportionalColumnPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
@@ -285,8 +288,6 @@ class ProportionalColumnPreprocessor : public Preprocessor {
DenseRow new_upper_bounds_;
ColumnDeletionHelper column_deletion_helper_;
DISALLOW_COPY_AND_ASSIGN(ProportionalColumnPreprocessor);
};
// --------------------------------------------------------
@@ -297,7 +298,11 @@ class ProportionalColumnPreprocessor : public Preprocessor {
// same remark above for columns in ProportionalColumnPreprocessor.
class ProportionalRowPreprocessor : public Preprocessor {
public:
ProportionalRowPreprocessor() {}
explicit ProportionalRowPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
ProportionalRowPreprocessor(const ProportionalRowPreprocessor&) = delete;
ProportionalRowPreprocessor& operator=(const ProportionalRowPreprocessor&) =
delete;
~ProportionalRowPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
@@ -310,7 +315,6 @@ class ProportionalRowPreprocessor : public Preprocessor {
bool lp_is_maximization_problem_;
RowDeletionHelper row_deletion_helper_;
DISALLOW_COPY_AND_ASSIGN(ProportionalRowPreprocessor);
};
// --------------------------------------------------------
@@ -393,7 +397,10 @@ class SingletonUndo {
// each time we delete a row or a column, new singletons may be created.
class SingletonPreprocessor : public Preprocessor {
public:
SingletonPreprocessor() {}
explicit SingletonPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
SingletonPreprocessor(const SingletonPreprocessor&) = delete;
SingletonPreprocessor& operator=(const SingletonPreprocessor&) = delete;
~SingletonPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
@@ -463,8 +470,6 @@ class SingletonPreprocessor : public Preprocessor {
// The transpose of the rows that are deleted by this preprocessor.
// TODO(user): implement a RowMajorSparseMatrix class to simplify the code.
SparseMatrix deleted_rows_;
DISALLOW_COPY_AND_ASSIGN(SingletonPreprocessor);
};
// --------------------------------------------------------
@@ -473,14 +478,17 @@ class SingletonPreprocessor : public Preprocessor {
// Removes the fixed variables from the problem.
class FixedVariablePreprocessor : public Preprocessor {
public:
FixedVariablePreprocessor() {}
explicit FixedVariablePreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
FixedVariablePreprocessor(const FixedVariablePreprocessor&) = delete;
FixedVariablePreprocessor& operator=(const FixedVariablePreprocessor&) =
delete;
~FixedVariablePreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
private:
ColumnDeletionHelper column_deletion_helper_;
DISALLOW_COPY_AND_ASSIGN(FixedVariablePreprocessor);
};
// --------------------------------------------------------
@@ -505,7 +513,13 @@ class FixedVariablePreprocessor : public Preprocessor {
// * Otherwise, wo do nothing.
class ForcingAndImpliedFreeConstraintPreprocessor : public Preprocessor {
public:
ForcingAndImpliedFreeConstraintPreprocessor() {}
explicit ForcingAndImpliedFreeConstraintPreprocessor(
const GlopParameters* parameters)
: Preprocessor(parameters) {}
ForcingAndImpliedFreeConstraintPreprocessor(
const ForcingAndImpliedFreeConstraintPreprocessor&) = delete;
ForcingAndImpliedFreeConstraintPreprocessor& operator=(
const ForcingAndImpliedFreeConstraintPreprocessor&) = delete;
~ForcingAndImpliedFreeConstraintPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
@@ -517,8 +531,6 @@ class ForcingAndImpliedFreeConstraintPreprocessor : public Preprocessor {
DenseBooleanColumn is_forcing_up_;
ColumnDeletionHelper column_deletion_helper_;
RowDeletionHelper row_deletion_helper_;
DISALLOW_COPY_AND_ASSIGN(ForcingAndImpliedFreeConstraintPreprocessor);
};
// --------------------------------------------------------
@@ -547,7 +559,10 @@ class ForcingAndImpliedFreeConstraintPreprocessor : public Preprocessor {
// problem thanks to the DoubletonFreeColumnPreprocessor.
class ImpliedFreePreprocessor : public Preprocessor {
public:
ImpliedFreePreprocessor() {}
explicit ImpliedFreePreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
ImpliedFreePreprocessor(const ImpliedFreePreprocessor&) = delete;
ImpliedFreePreprocessor& operator=(const ImpliedFreePreprocessor&) = delete;
~ImpliedFreePreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
@@ -562,8 +577,6 @@ class ImpliedFreePreprocessor : public Preprocessor {
// value of these variables; which will only be used (eg. restored) if the
// variable actually turns out to be VariableStatus::FREE.
VariableStatusRow postsolve_status_of_free_variables_;
DISALLOW_COPY_AND_ASSIGN(ImpliedFreePreprocessor);
};
// --------------------------------------------------------
@@ -592,7 +605,12 @@ class ImpliedFreePreprocessor : public Preprocessor {
// required. Most probably, commercial solvers do use it though.
class DoubletonFreeColumnPreprocessor : public Preprocessor {
public:
DoubletonFreeColumnPreprocessor() {}
explicit DoubletonFreeColumnPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
DoubletonFreeColumnPreprocessor(const DoubletonFreeColumnPreprocessor&) =
delete;
DoubletonFreeColumnPreprocessor& operator=(
const DoubletonFreeColumnPreprocessor&) = delete;
~DoubletonFreeColumnPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
@@ -621,7 +639,6 @@ class DoubletonFreeColumnPreprocessor : public Preprocessor {
std::vector<RestoreInfo> restore_stack_;
RowDeletionHelper row_deletion_helper_;
DISALLOW_COPY_AND_ASSIGN(DoubletonFreeColumnPreprocessor);
};
// --------------------------------------------------------
@@ -639,7 +656,12 @@ class DoubletonFreeColumnPreprocessor : public Preprocessor {
// the Andersen & Andersen paper.
class UnconstrainedVariablePreprocessor : public Preprocessor {
public:
UnconstrainedVariablePreprocessor() {}
explicit UnconstrainedVariablePreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
UnconstrainedVariablePreprocessor(const UnconstrainedVariablePreprocessor&) =
delete;
UnconstrainedVariablePreprocessor& operator=(
const UnconstrainedVariablePreprocessor&) = delete;
~UnconstrainedVariablePreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
@@ -680,8 +702,6 @@ class UnconstrainedVariablePreprocessor : public Preprocessor {
SparseMatrix deleted_columns_;
SparseMatrix deleted_rows_as_column_;
DISALLOW_COPY_AND_ASSIGN(UnconstrainedVariablePreprocessor);
};
// --------------------------------------------------------
@@ -690,14 +710,17 @@ class UnconstrainedVariablePreprocessor : public Preprocessor {
// Removes the constraints with no bounds from the problem.
class FreeConstraintPreprocessor : public Preprocessor {
public:
FreeConstraintPreprocessor() {}
explicit FreeConstraintPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
FreeConstraintPreprocessor(const FreeConstraintPreprocessor&) = delete;
FreeConstraintPreprocessor& operator=(const FreeConstraintPreprocessor&) =
delete;
~FreeConstraintPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
private:
RowDeletionHelper row_deletion_helper_;
DISALLOW_COPY_AND_ASSIGN(FreeConstraintPreprocessor);
};
// --------------------------------------------------------
@@ -706,14 +729,17 @@ class FreeConstraintPreprocessor : public Preprocessor {
// Removes the constraints with no coefficients from the problem.
class EmptyConstraintPreprocessor : public Preprocessor {
public:
EmptyConstraintPreprocessor() {}
explicit EmptyConstraintPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
EmptyConstraintPreprocessor(const EmptyConstraintPreprocessor&) = delete;
EmptyConstraintPreprocessor& operator=(const EmptyConstraintPreprocessor&) =
delete;
~EmptyConstraintPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
private:
RowDeletionHelper row_deletion_helper_;
DISALLOW_COPY_AND_ASSIGN(EmptyConstraintPreprocessor);
};
// --------------------------------------------------------
@@ -728,13 +754,17 @@ class EmptyConstraintPreprocessor : public Preprocessor {
// coefficients are small! Run it after ScalingPreprocessor or fix the code.
class RemoveNearZeroEntriesPreprocessor : public Preprocessor {
public:
RemoveNearZeroEntriesPreprocessor() {}
explicit RemoveNearZeroEntriesPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
RemoveNearZeroEntriesPreprocessor(const RemoveNearZeroEntriesPreprocessor&) =
delete;
RemoveNearZeroEntriesPreprocessor& operator=(
const RemoveNearZeroEntriesPreprocessor&) = delete;
~RemoveNearZeroEntriesPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
private:
DISALLOW_COPY_AND_ASSIGN(RemoveNearZeroEntriesPreprocessor);
};
// --------------------------------------------------------
@@ -746,14 +776,18 @@ class RemoveNearZeroEntriesPreprocessor : public Preprocessor {
// efficient solve when this column is involved.
class SingletonColumnSignPreprocessor : public Preprocessor {
public:
SingletonColumnSignPreprocessor() {}
explicit SingletonColumnSignPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
SingletonColumnSignPreprocessor(const SingletonColumnSignPreprocessor&) =
delete;
SingletonColumnSignPreprocessor& operator=(
const SingletonColumnSignPreprocessor&) = delete;
~SingletonColumnSignPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
private:
std::vector<ColIndex> changed_columns_;
DISALLOW_COPY_AND_ASSIGN(SingletonColumnSignPreprocessor);
};
// --------------------------------------------------------
@@ -764,7 +798,12 @@ class SingletonColumnSignPreprocessor : public Preprocessor {
// in all the constraints that it is involved in.
class DoubletonEqualityRowPreprocessor : public Preprocessor {
public:
DoubletonEqualityRowPreprocessor() {}
explicit DoubletonEqualityRowPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
DoubletonEqualityRowPreprocessor(const DoubletonEqualityRowPreprocessor&) =
delete;
DoubletonEqualityRowPreprocessor& operator=(
const DoubletonEqualityRowPreprocessor&) = delete;
~DoubletonEqualityRowPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
@@ -818,8 +857,6 @@ class DoubletonEqualityRowPreprocessor : public Preprocessor {
std::vector<RestoreInfo> restore_stack_;
DenseColumn saved_row_lower_bounds_;
DenseColumn saved_row_upper_bounds_;
DISALLOW_COPY_AND_ASSIGN(DoubletonEqualityRowPreprocessor);
};
// Because of numerical imprecision, a preprocessor like
@@ -848,7 +885,10 @@ void FixConstraintWithFixedStatuses(const DenseColumn& row_lower_bounds,
// preprocessor does not deal correctly with free constraints.
class DualizerPreprocessor : public Preprocessor {
public:
DualizerPreprocessor() {}
explicit DualizerPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
DualizerPreprocessor(const DualizerPreprocessor&) = delete;
DualizerPreprocessor& operator=(const DualizerPreprocessor&) = delete;
~DualizerPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
@@ -872,8 +912,6 @@ class DualizerPreprocessor : public Preprocessor {
// For postsolving the variable/constraint statuses.
VariableStatusRow dual_status_correspondence_;
ColMapping slack_or_surplus_mapping_;
DISALLOW_COPY_AND_ASSIGN(DualizerPreprocessor);
};
// --------------------------------------------------------
@@ -905,7 +943,12 @@ class DualizerPreprocessor : public Preprocessor {
// a free variable so that only having a domain containing 0.0 is enough?
class ShiftVariableBoundsPreprocessor : public Preprocessor {
public:
ShiftVariableBoundsPreprocessor() {}
explicit ShiftVariableBoundsPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
ShiftVariableBoundsPreprocessor(const ShiftVariableBoundsPreprocessor&) =
delete;
ShiftVariableBoundsPreprocessor& operator=(
const ShiftVariableBoundsPreprocessor&) = delete;
~ShiftVariableBoundsPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
@@ -919,7 +962,6 @@ class ShiftVariableBoundsPreprocessor : public Preprocessor {
// numerical accuracy for variables at their bound after postsolve.
DenseRow variable_initial_lbs_;
DenseRow variable_initial_ubs_;
DISALLOW_COPY_AND_ASSIGN(ShiftVariableBoundsPreprocessor);
};
// --------------------------------------------------------
@@ -929,7 +971,10 @@ class ShiftVariableBoundsPreprocessor : public Preprocessor {
// This is only applied if the parameter use_scaling is true.
class ScalingPreprocessor : public Preprocessor {
public:
ScalingPreprocessor() {}
explicit ScalingPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
ScalingPreprocessor(const ScalingPreprocessor&) = delete;
ScalingPreprocessor& operator=(const ScalingPreprocessor&) = delete;
~ScalingPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
@@ -941,8 +986,6 @@ class ScalingPreprocessor : public Preprocessor {
Fractional cost_scaling_factor_;
Fractional bound_scaling_factor_;
SparseMatrixScaler scaler_;
DISALLOW_COPY_AND_ASSIGN(ScalingPreprocessor);
};
// --------------------------------------------------------
@@ -953,13 +996,14 @@ class ScalingPreprocessor : public Preprocessor {
// The preprocessor is kept here, because it could be used by Glop too.
class ToMinimizationPreprocessor : public Preprocessor {
public:
ToMinimizationPreprocessor() {}
explicit ToMinimizationPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
ToMinimizationPreprocessor(const ToMinimizationPreprocessor&) = delete;
ToMinimizationPreprocessor& operator=(const ToMinimizationPreprocessor&) =
delete;
~ToMinimizationPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
private:
DISALLOW_COPY_AND_ASSIGN(ToMinimizationPreprocessor);
};
// --------------------------------------------------------
@@ -977,54 +1021,17 @@ class ToMinimizationPreprocessor : public Preprocessor {
// so that the rightmost square sub-matrix is always the identity matrix.
class AddSlackVariablesPreprocessor : public Preprocessor {
public:
AddSlackVariablesPreprocessor() {}
explicit AddSlackVariablesPreprocessor(const GlopParameters* parameters)
: Preprocessor(parameters) {}
AddSlackVariablesPreprocessor(const AddSlackVariablesPreprocessor&) = delete;
AddSlackVariablesPreprocessor& operator=(
const AddSlackVariablesPreprocessor&) = delete;
~AddSlackVariablesPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
private:
ColIndex first_slack_col_;
DISALLOW_COPY_AND_ASSIGN(AddSlackVariablesPreprocessor);
};
// --------------------------------------------------------
// SolowHalimPreprocessor
// --------------------------------------------------------
// Modifies the problem by changing variables to begin
// geometrically near to the optimal solution according to the suggestion of
// Solow & Halim.
// Variable changes have a very simple form :
// x_j' = upper_bound(j) - x_j if cost(j) > 0
// x_j' = -lower_bound(j) + x_j if cost(j) < 0
// for a maximization problem
//
// reference:
// Improving the Efficiency of the Simplex Algorithm Based on a
// Geometric Explanation of Phase 1
// https://weatherhead.case.edu/faculty/research/library/detail?id=12128593921
// DOI: 10.1504/IJOR.2009.025701
class SolowHalimPreprocessor : public Preprocessor {
public:
SolowHalimPreprocessor() {}
~SolowHalimPreprocessor() final {}
bool Run(LinearProgram* lp) final;
void RecoverSolution(ProblemSolution* solution) const final;
private:
typedef enum {
NOT_MODIFIED = 0,
SHIFTED = 1,
SHIFTED_OPPOSITE_DIRECTION = 2
} ColumnTransformType;
// Contains the coordinate change information for each column
gtl::ITIVector<ColIndex, ColumnTransformType> column_transform_;
// Contains the initial problem bounds.
DenseRow variable_initial_lbs_;
DenseRow variable_initial_ubs_;
DISALLOW_COPY_AND_ASSIGN(SolowHalimPreprocessor);
};
} // namespace glop

View File

@@ -749,8 +749,8 @@ void MPSolver::ExportModelToProto(MPModelProto* output_model) const {
}
}
util::Status MPSolver::LoadSolutionFromProto(
const MPSolutionResponse& response) {
util::Status MPSolver::LoadSolutionFromProto(const MPSolutionResponse& response,
double tolerance) {
interface_->result_status_ = static_cast<ResultStatus>(response.status());
if (response.status() != MPSOLVER_OPTIMAL &&
response.status() != MPSOLVER_FEASIBLE) {
@@ -768,34 +768,35 @@ util::Status MPSolver::LoadSolutionFromProto(
response.variable_value_size(),
") does not correspond to the Solver's (", variables_.size(), ")"));
}
double largest_error = 0;
interface_->ExtractModel();
// Look further: verify that the variable values are within the bounds.
int num_vars_out_of_bounds = 0;
const double tolerance = MPSolverParameters::kDefaultPrimalTolerance;
int last_offending_var = -1;
for (int i = 0; i < response.variable_value_size(); ++i) {
const double var_value = response.variable_value(i);
MPVariable* var = variables_[i];
// TODO(user): Use parameter when they become available in this class.
const double lb_error = var->lb() - var_value;
const double ub_error = var_value - var->ub();
if (lb_error > tolerance || ub_error > tolerance) {
++num_vars_out_of_bounds;
largest_error = std::max(largest_error, std::max(lb_error, ub_error));
last_offending_var = i;
if (tolerance != infinity()) {
// Look further: verify that the variable values are within the bounds.
double largest_error = 0;
int num_vars_out_of_bounds = 0;
int last_offending_var = -1;
for (int i = 0; i < response.variable_value_size(); ++i) {
const double var_value = response.variable_value(i);
MPVariable* var = variables_[i];
// TODO(user): Use parameter when they become available in this class.
const double lb_error = var->lb() - var_value;
const double ub_error = var_value - var->ub();
if (lb_error > tolerance || ub_error > tolerance) {
++num_vars_out_of_bounds;
largest_error = std::max(largest_error, std::max(lb_error, ub_error));
last_offending_var = i;
}
}
if (num_vars_out_of_bounds > 0) {
return util::InvalidArgumentError(absl::StrCat(
"Loaded a solution whose variables matched the solver's, but ",
num_vars_out_of_bounds, " of ", variables_.size(),
" variables were out of their bounds, by more than the primal"
" tolerance which is: ",
tolerance, ". Max error: ", largest_error, ", last offender var is #",
last_offending_var, ": '", variables_[last_offending_var]->name(),
"'"));
}
}
if (num_vars_out_of_bounds > 0) {
return util::InvalidArgumentError(absl::StrCat(
"Loaded a solution whose variables matched the solver's, but ",
num_vars_out_of_bounds, " of ", variables_.size(),
" variables were out of their bounds, by more than the primal"
" tolerance which is: ",
tolerance, ". Max error: ", largest_error, ", last offendir var is #",
last_offending_var, ": '", variables_[last_offending_var]->name(),
"'"));
}
for (int i = 0; i < response.variable_value_size(); ++i) {
variables_[i]->set_solution_value(response.variable_value(i));
@@ -1075,6 +1076,24 @@ std::string PrettyPrintConstraint(const MPConstraint& constraint) {
}
} // namespace
util::Status MPSolver::ClampSolutionWithinBounds() {
interface_->ExtractModel();
for (MPVariable* const variable : variables_) {
const double value = variable->solution_value();
if (std::isnan(value)) {
return util::InvalidArgumentError(
absl::StrCat("NaN value for ", PrettyPrintVar(*variable)));
}
if (value < variable->lb()) {
variable->set_solution_value(variable->lb());
} else if (value > variable->ub()) {
variable->set_solution_value(variable->ub());
}
}
interface_->sync_status_ = MPSolverInterface::SOLUTION_SYNCHRONIZED;
return util::OkStatus();
}
std::vector<double> MPSolver::ComputeConstraintActivities() const {
// TODO(user): test this failure case.
if (!interface_->CheckSolutionIsSynchronizedAndExists()) return {};
@@ -1354,7 +1373,7 @@ bool MPSolverInterface::CheckSolutionIsSynchronized() const {
if (sync_status_ != SOLUTION_SYNCHRONIZED) {
LOG(DFATAL)
<< "The model has been changed since the solution was last computed."
<< " MPSolverInterface::status_ = " << sync_status_;
<< " MPSolverInterface::sync_status_ = " << sync_status_;
return false;
}
return true;
@@ -1527,7 +1546,8 @@ std::string MPSolverInterface::ValidFileExtensionForParameterFile() const {
const double MPSolverParameters::kDefaultRelativeMipGap = 1e-4;
// For the primal and dual tolerances, choose the same default as CLP and GLPK.
const double MPSolverParameters::kDefaultPrimalTolerance = 1e-7;
const double MPSolverParameters::kDefaultPrimalTolerance =
operations_research::kDefaultPrimalTolerance;
const double MPSolverParameters::kDefaultDualTolerance = 1e-7;
const MPSolverParameters::PresolveValues MPSolverParameters::kDefaultPresolve =
MPSolverParameters::PRESOLVE_ON;

View File

@@ -163,6 +163,8 @@
namespace operations_research {
constexpr double kDefaultPrimalTolerance = 1e-07;
class MPConstraint;
class MPObjective;
class MPSolverInterface;
@@ -449,7 +451,15 @@ class MPSolver {
// - loading a solution with a status other than OPTIMAL / FEASIBLE.
// Note: the objective value isnn't checked. You can use VerifySolution()
// for that.
util::Status LoadSolutionFromProto(const MPSolutionResponse& response);
// TODO(b/116117536) split this into two separate functions: Load...() without
// checking for tolerance and SolutionIsFeasibleWithTolerance().
util::Status LoadSolutionFromProto(
const MPSolutionResponse& response,
double tolerance = kDefaultPrimalTolerance);
// Resets values of out of bound variables to the corresponding bound
// and returns an error if any of the variables have NaN value.
util::Status ClampSolutionWithinBounds();
// ----- Export model to files or strings -----
// Shortcuts to the homonymous MPModelProtoExporter methods, via

View File

@@ -52,43 +52,42 @@ namespace {
// in-memory domain of each variables and the constraint variable graph.
struct PresolveContext {
bool DomainIsEmpty(int ref) const {
return domains[PositiveRef(ref)].empty();
return domains[PositiveRef(ref)].IsEmpty();
}
bool IsFixed(int ref) const {
CHECK(!DomainIsEmpty(ref));
return domains[PositiveRef(ref)].front().start ==
domains[PositiveRef(ref)].back().end;
return domains[PositiveRef(ref)].Min() == domains[PositiveRef(ref)].Max();
}
bool LiteralIsTrue(int lit) const {
if (!IsFixed(lit)) return false;
if (RefIsPositive(lit)) {
return domains[lit].front().start == 1ll;
return domains[lit].Min() == 1ll;
} else {
return domains[PositiveRef(lit)].front().start == 0ll;
return domains[PositiveRef(lit)].Min() == 0ll;
}
}
bool LiteralIsFalse(int lit) const {
if (!IsFixed(lit)) return false;
if (RefIsPositive(lit)) {
return domains[lit].front().start == 0ll;
return domains[lit].Min() == 0ll;
} else {
return domains[PositiveRef(lit)].front().start == 1ll;
return domains[PositiveRef(lit)].Min() == 1ll;
}
}
int64 MinOf(int ref) const {
CHECK(!DomainIsEmpty(ref));
return RefIsPositive(ref) ? domains[PositiveRef(ref)].front().start
: -domains[PositiveRef(ref)].back().end;
return RefIsPositive(ref) ? domains[PositiveRef(ref)].Min()
: -domains[PositiveRef(ref)].Max();
}
int64 MaxOf(int ref) const {
CHECK(!DomainIsEmpty(ref));
return RefIsPositive(ref) ? domains[PositiveRef(ref)].back().end
: -domains[PositiveRef(ref)].front().start;
return RefIsPositive(ref) ? domains[PositiveRef(ref)].Max()
: -domains[PositiveRef(ref)].Min();
}
// Returns true if this ref only appear in one constraint.
@@ -96,27 +95,28 @@ struct PresolveContext {
return var_to_constraints[PositiveRef(ref)].size() == 1;
}
std::vector<ClosedInterval> GetRefDomain(int ref) const {
Domain DomainOf(int ref) const {
if (RefIsPositive(ref)) return domains[ref];
return NegationOfSortedDisjointIntervals(domains[PositiveRef(ref)]);
return domains[PositiveRef(ref)].Negation();
}
// Returns false if the domain changed.
bool IntersectDomainWith(int ref, const std::vector<ClosedInterval>& domain) {
// Returns true iff the domain changed.
bool IntersectDomainWith(int ref, const Domain& domain) {
CHECK(!DomainIsEmpty(ref));
const int var = PositiveRef(ref);
tmp_domain = IntersectionOfSortedDisjointIntervals(
domains[var], RefIsPositive(ref)
? domain
: NegationOfSortedDisjointIntervals(domain));
if (tmp_domain != domains[var]) {
modified_domains.Set(var);
domains[var] = tmp_domain;
if (tmp_domain.empty()) {
is_unsat = true;
}
return true;
if (RefIsPositive(ref)) {
if (domains[var].IsIncludedIn(domain)) return false;
domains[var] = domains[var].IntersectionWith(domain);
} else {
const Domain temp = domain.Negation();
if (domains[var].IsIncludedIn(temp)) return false;
domains[var] = domains[var].IntersectionWith(temp);
}
return false;
modified_domains.Set(var);
if (domains[var].IsEmpty()) is_unsat = true;
return true;
}
void SetLiteralToFalse(int lit) {
@@ -128,7 +128,7 @@ struct PresolveContext {
is_unsat = true;
}
} else {
IntersectDomainWith(var, {{value, value}});
IntersectDomainWith(var, Domain(value));
}
}
@@ -240,7 +240,7 @@ struct PresolveContext {
// Create the internal structure for any new variables in working_model.
void InitializeNewDomains() {
for (int i = domains.size(); i < working_model->variables_size(); ++i) {
domains.push_back(ReadDomain(working_model->variables(i)));
domains.push_back(ReadDomainFromProto(working_model->variables(i)));
if (IsFixed(i)) ExploitFixedDomain(i);
}
modified_domains.Resize(domains.size());
@@ -290,16 +290,15 @@ struct PresolveContext {
// Temporary storage.
std::vector<int> tmp_literals;
std::vector<ClosedInterval> tmp_domain;
std::vector<std::vector<ClosedInterval>> tmp_term_domains;
std::vector<std::vector<ClosedInterval>> tmp_left_domains;
std::vector<Domain> tmp_term_domains;
std::vector<Domain> tmp_left_domains;
// Each time a domain is modified this is set to true.
SparseBitset<int64> modified_domains;
private:
// The current domain of each variables.
std::vector<std::vector<ClosedInterval>> domains;
std::vector<Domain> domains;
};
// =============================================================================
@@ -529,12 +528,10 @@ bool PresolveIntMax(ConstraintProto* ct, PresolveContext* context) {
// Update the target domain.
bool domain_reduced = false;
if (!HasEnforcementLiteral(*ct)) {
std::vector<ClosedInterval> infered_domain;
Domain infered_domain;
for (const int ref : ct->int_max().vars()) {
infered_domain = UnionOfSortedDisjointIntervals(
infered_domain,
IntersectionOfSortedDisjointIntervals(context->GetRefDomain(ref),
{{target_min, target_max}}));
infered_domain = infered_domain.UnionWith(
context->DomainOf(ref).IntersectionWith({target_min, target_max}));
}
domain_reduced |= context->IntersectDomainWith(target_ref, infered_domain);
}
@@ -546,7 +543,7 @@ bool PresolveIntMax(ConstraintProto* ct, PresolveContext* context) {
for (const int ref : ct->int_max().vars()) {
if (!HasEnforcementLiteral(*ct)) {
domain_reduced |=
context->IntersectDomainWith(ref, {{kint64min, target_max}});
context->IntersectDomainWith(ref, Domain(kint64min, target_max));
}
if (context->MaxOf(ref) >= target_min) {
ct->mutable_int_max()->set_vars(new_size++, ref);
@@ -626,7 +623,7 @@ bool PresolveIntProd(ConstraintProto* ct, PresolveContext* context) {
}
// This is a bool constraint!
context->IntersectDomainWith(target_ref, {{0, 1}});
context->IntersectDomainWith(target_ref, Domain(0, 1));
context->UpdateRuleStats("int_prod: all Boolean.");
{
ConstraintProto* new_ct = context->working_model->add_constraints();
@@ -662,8 +659,7 @@ bool PresolveIntDiv(ConstraintProto* ct, PresolveContext* context) {
context->UpdateRuleStats("TODO int_div: rewrite to equality");
}
if (context->IntersectDomainWith(
target, DivisionOfSortedDisjointIntervals(
context->GetRefDomain(ref_x), divisor))) {
target, context->DomainOf(ref_x).DivisionBy(divisor))) {
context->UpdateRuleStats(
"int_div: updated domain of target in target = X / cte");
}
@@ -721,7 +717,7 @@ bool ExploitEquivalenceRelations(ConstraintProto* ct,
bool PresolveLinear(ConstraintProto* ct, PresolveContext* context) {
bool var_constraint_graph_changed = false;
std::vector<ClosedInterval> rhs = ReadDomain(ct->linear());
Domain rhs = ReadDomainFromProto(ct->linear());
// First regroup the terms on the same variables and sum the fixed ones.
// Note that we use a map to sort the variables and because we expect most
@@ -781,15 +777,15 @@ bool PresolveLinear(ConstraintProto* ct, PresolveContext* context) {
if (context->IsUnique(var) &&
context->affine_relations.ClassSize(var) == 1) {
bool success;
const auto term_domain = PreciseMultiplicationOfSortedDisjointIntervals(
context->GetRefDomain(var), -coeff, &success);
const auto term_domain =
context->DomainOf(var).MultiplicationBy(-coeff, &success);
if (success) {
// Note that we can't do that if we loose information in the
// multiplication above because the new domain might not be as strict
// as the initial constraint otherwise. TODO(user): because of the
// addition, it might be possible to cover more cases though.
var_to_erase.push_back(var);
rhs = AdditionOfSortedDisjointIntervals(rhs, term_domain);
rhs = rhs.AdditionWith(term_domain);
continue;
}
}
@@ -837,11 +833,10 @@ bool PresolveLinear(ConstraintProto* ct, PresolveContext* context) {
// Rewrite the constraint in canonical form and update rhs (it will be copied
// to the constraint later).
if (sum_of_fixed_terms != 0) {
rhs = AdditionOfSortedDisjointIntervals(
rhs, {{-sum_of_fixed_terms, -sum_of_fixed_terms}});
rhs = rhs.AdditionWith({-sum_of_fixed_terms, -sum_of_fixed_terms});
}
if (gcd > 1) {
rhs = InverseMultiplicationOfSortedDisjointIntervals(rhs, gcd);
rhs = rhs.InverseMultiplicationBy(gcd);
}
ct->mutable_linear()->clear_vars();
ct->mutable_linear()->clear_coeffs();
@@ -854,7 +849,7 @@ bool PresolveLinear(ConstraintProto* ct, PresolveContext* context) {
// Empty constraint?
if (ct->linear().vars().empty()) {
context->UpdateRuleStats("linear: empty");
if (SortedDisjointIntervalsContain(rhs, 0)) {
if (rhs.Contains(0)) {
return RemoveConstraint(ct, context);
} else {
return MarkConstraintAsFalse(ct, context);
@@ -871,7 +866,7 @@ bool PresolveLinear(ConstraintProto* ct, PresolveContext* context) {
context->IntersectDomainWith(var, rhs);
} else {
DCHECK_EQ(coeff, -1); // Because of the GCD above.
context->IntersectDomainWith(var, NegationOfSortedDisjointIntervals(rhs));
context->IntersectDomainWith(var, rhs.Negation());
}
return RemoveConstraint(ct, context);
}
@@ -883,69 +878,69 @@ bool PresolveLinear(ConstraintProto* ct, PresolveContext* context) {
const int num_vars = arg.vars_size();
term_domains.resize(num_vars + 1);
left_domains.resize(num_vars + 1);
left_domains[0] = {{0, 0}};
left_domains[0] = Domain(0);
for (int i = 0; i < num_vars; ++i) {
const int var = PositiveRef(arg.vars(i));
const int64 coeff = arg.coeffs(i);
const auto& domain = context->GetRefDomain(var);
// TODO(user): Try PreciseMultiplicationOfSortedDisjointIntervals() if
// the size is reasonable.
term_domains[i] = MultiplicationOfSortedDisjointIntervals(domain, coeff);
left_domains[i + 1] =
AdditionOfSortedDisjointIntervals(left_domains[i], term_domains[i]);
if (left_domains[i + 1].size() > kDomainComplexityLimit) {
term_domains[i] = context->DomainOf(var).ContinuousMultiplicationBy(coeff);
left_domains[i + 1] = left_domains[i].AdditionWith(term_domains[i]);
if (left_domains[i + 1].intervals().size() > kDomainComplexityLimit) {
// We take a super-set, otherwise it will be too slow.
//
// TODO(user): We could be smarter in how we compute this if we allow for
// more than one intervals.
left_domains[i + 1] = {
{left_domains[i + 1].front().start, left_domains[i + 1].back().end}};
left_domains[i + 1] =
Domain(left_domains[i + 1].Min(), left_domains[i + 1].Max());
}
}
const std::vector<ClosedInterval>& implied_rhs = left_domains[num_vars];
const Domain& implied_rhs = left_domains[num_vars];
// Abort if intersection is empty.
const std::vector<ClosedInterval> restricted_rhs =
IntersectionOfSortedDisjointIntervals(rhs, implied_rhs);
if (restricted_rhs.empty()) {
const Domain restricted_rhs = rhs.IntersectionWith(implied_rhs);
if (restricted_rhs.IsEmpty()) {
context->UpdateRuleStats("linear: infeasible");
return MarkConstraintAsFalse(ct, context);
}
// Relax the constraint rhs for faster propagation.
// TODO(user): add an IntersectionIsEmpty() function.
rhs.clear();
for (const ClosedInterval i : UnionOfSortedDisjointIntervals(
restricted_rhs, ComplementOfSortedDisjointIntervals(implied_rhs))) {
if (!IntersectionOfSortedDisjointIntervals({i}, restricted_rhs).empty()) {
rhs.push_back(i);
std::vector<ClosedInterval> rhs_intervals;
for (const ClosedInterval i :
restricted_rhs.UnionWith(implied_rhs.Complement()).intervals()) {
if (!Domain::FromIntervals({i})
.IntersectionWith(restricted_rhs)
.IsEmpty()) {
rhs_intervals.push_back(i);
}
}
if (rhs.size() == 1 && rhs[0].start == kint64min && rhs[0].end == kint64max) {
rhs = Domain::FromIntervals(rhs_intervals);
if (rhs == Domain::AllValues()) {
context->UpdateRuleStats("linear: always true");
return RemoveConstraint(ct, context);
}
if (rhs != ReadDomain(ct->linear())) {
if (rhs != ReadDomainFromProto(ct->linear())) {
context->UpdateRuleStats("linear: simplified rhs");
}
FillDomain(rhs, ct->mutable_linear());
FillDomainInProto(rhs, ct->mutable_linear());
// Propagate the variable bounds.
if (!HasEnforcementLiteral(*ct)) {
bool new_bounds = false;
std::vector<ClosedInterval> new_domain;
std::vector<ClosedInterval> right_domain = {{0, 0}};
term_domains[num_vars] = NegationOfSortedDisjointIntervals(rhs);
Domain new_domain;
Domain right_domain(0, 0);
term_domains[num_vars] = rhs.Negation();
for (int i = num_vars - 1; i >= 0; --i) {
right_domain =
AdditionOfSortedDisjointIntervals(right_domain, term_domains[i + 1]);
if (right_domain.size() > kDomainComplexityLimit) {
right_domain = right_domain.AdditionWith(term_domains[i + 1]);
if (right_domain.intervals().size() > kDomainComplexityLimit) {
// We take a super-set, otherwise it will be too slow.
right_domain = {{right_domain.front().start, right_domain.back().end}};
right_domain = Domain(right_domain.Min(), right_domain.Max());
}
new_domain = InverseMultiplicationOfSortedDisjointIntervals(
AdditionOfSortedDisjointIntervals(left_domains[i], right_domain),
-arg.coeffs(i));
new_domain = left_domains[i]
.AdditionWith(right_domain)
.InverseMultiplicationBy(-arg.coeffs(i));
if (context->IntersectDomainWith(arg.vars(i), new_domain)) {
new_bounds = true;
}
@@ -961,8 +956,8 @@ bool PresolveLinear(ConstraintProto* ct, PresolveContext* context) {
// a coefficient of magnitude 1, and later the one with larger coeffs.
if (!was_affine && !HasEnforcementLiteral(*ct)) {
const LinearConstraintProto& arg = ct->linear();
const int64 rhs_min = rhs.front().start;
const int64 rhs_max = rhs.back().end;
const int64 rhs_min = rhs.Min();
const int64 rhs_max = rhs.Max();
if (rhs_min == rhs_max && arg.vars_size() == 2) {
const int v1 = arg.vars(0);
const int v2 = arg.vars(1);
@@ -1008,9 +1003,9 @@ bool PresolveLinearIntoClauses(ConstraintProto* ct, PresolveContext* context) {
// Detect clauses and reified ands.
// TODO(user): split an == 1 constraint or similar into a clause and a <= 1
// constraint?
const std::vector<ClosedInterval> domain = ReadDomain(arg);
DCHECK(!domain.empty());
if (offset + min_coeff > domain.back().end) {
const Domain domain = ReadDomainFromProto(arg);
DCHECK(!domain.IsEmpty());
if (offset + min_coeff > domain.Max()) {
// All Boolean are false if the reified literal is true.
context->UpdateRuleStats("linear: reified and");
const auto copy = arg;
@@ -1020,8 +1015,8 @@ bool PresolveLinearIntoClauses(ConstraintProto* ct, PresolveContext* context) {
copy.coeffs(i) > 0 ? NegatedRef(copy.vars(i)) : copy.vars(i));
}
return PresolveBoolAnd(ct, context);
} else if (offset + min_coeff >= domain[0].start &&
domain[0].end == kint64max) {
} else if (offset + min_coeff >= domain.Min() &&
domain.intervals()[0].end == kint64max) {
// At least one Boolean is true.
context->UpdateRuleStats("linear: clause");
const auto copy = arg;
@@ -1045,7 +1040,7 @@ bool PresolveLinearIntoClauses(ConstraintProto* ct, PresolveContext* context) {
for (int i = 0; i < num_vars; ++i) {
if ((mask >> i) & 1) value += arg.coeffs(i);
}
if (SortedDisjointIntervalsContain(domain, value)) continue;
if (domain.Contains(value)) continue;
// Add a new clause to exclude this bad assignment.
ConstraintProto* new_ct = context->working_model->add_constraints();
@@ -1069,16 +1064,13 @@ bool PresolveInterval(ConstraintProto* ct, PresolveContext* context) {
const int size = ct->interval().size();
bool changed = false;
changed |= context->IntersectDomainWith(
end, AdditionOfSortedDisjointIntervals(context->GetRefDomain(start),
context->GetRefDomain(size)));
end, context->DomainOf(start).AdditionWith(context->DomainOf(size)));
changed |= context->IntersectDomainWith(
start, AdditionOfSortedDisjointIntervals(
context->GetRefDomain(end), NegationOfSortedDisjointIntervals(
context->GetRefDomain(size))));
start,
context->DomainOf(end).AdditionWith(context->DomainOf(size).Negation()));
changed |= context->IntersectDomainWith(
size, AdditionOfSortedDisjointIntervals(
context->GetRefDomain(end), NegationOfSortedDisjointIntervals(
context->GetRefDomain(start))));
size,
context->DomainOf(end).AdditionWith(context->DomainOf(start).Negation()));
if (changed) {
context->UpdateRuleStats("interval: reduced domains");
}
@@ -1110,45 +1102,34 @@ bool PresolveElement(ConstraintProto* ct, PresolveContext* context) {
bool all_included_in_target_domain = true;
bool reduced_index_domain = false;
const std::vector<ClosedInterval> index_domain =
context->GetRefDomain(index_ref);
if (index_domain.front().start < 0 ||
index_domain.back().end >= ct->element().vars_size()) {
if (context->IntersectDomainWith(index_ref,
Domain(0, ct->element().vars_size() - 1))) {
reduced_index_domain = true;
context->IntersectDomainWith(index_ref,
{{0, ct->element().vars_size() - 1}});
}
std::vector<ClosedInterval> infered_domain;
const std::vector<ClosedInterval> target_domain =
context->GetRefDomain(target_ref);
const std::vector<ClosedInterval> complement_of_target_domain =
ComplementOfSortedDisjointIntervals(context->GetRefDomain(target_ref));
for (const ClosedInterval interval : context->GetRefDomain(index_ref)) {
for (int i = interval.start; i <= interval.end; ++i) {
CHECK_GE(i, 0);
CHECK_LT(i, ct->element().vars_size());
const int ref = ct->element().vars(i);
const auto& domain = context->GetRefDomain(ref);
if (IntersectionOfSortedDisjointIntervals(target_domain, domain)
.empty()) {
context->IntersectDomainWith(index_ref,
{{kint64min, i - 1}, {i + 1, kint64max}});
Domain infered_domain;
const Domain target_domain = context->DomainOf(target_ref);
for (const ClosedInterval interval :
context->DomainOf(index_ref).intervals()) {
for (int value = interval.start; value <= interval.end; ++value) {
CHECK_GE(value, 0);
CHECK_LT(value, ct->element().vars_size());
const int ref = ct->element().vars(value);
const Domain domain = context->DomainOf(ref);
if (domain.IntersectionWith(target_domain).IsEmpty()) {
context->IntersectDomainWith(index_ref, Domain(value).Complement());
reduced_index_domain = true;
} else {
++num_vars;
if (domain.front().start == domain.back().end) {
constant_set.insert(domain.front().start);
if (domain.Min() == domain.Max()) {
constant_set.insert(domain.Min());
} else {
all_constants = false;
}
if (!IntersectionOfSortedDisjointIntervals(domain,
complement_of_target_domain)
.empty()) {
if (!domain.IsIncludedIn(target_domain)) {
all_included_in_target_domain = false;
}
infered_domain = UnionOfSortedDisjointIntervals(infered_domain, domain);
infered_domain = infered_domain.UnionWith(domain);
}
}
}
@@ -1156,6 +1137,7 @@ bool PresolveElement(ConstraintProto* ct, PresolveContext* context) {
context->UpdateRuleStats("element: reduced index domain");
}
if (context->IntersectDomainWith(target_ref, infered_domain)) {
if (context->DomainOf(target_ref).IsEmpty()) return true;
context->UpdateRuleStats("element: reduced target domain");
}
@@ -1169,6 +1151,7 @@ bool PresolveElement(ConstraintProto* ct, PresolveContext* context) {
*(context->mapping_model->add_constraints()) = *ct;
return RemoveConstraint(ct, context);
}
const bool unique_target =
context->IsUnique(target_ref) || context->IsFixed(target_ref);
if (all_included_in_target_domain && unique_target) {
@@ -1217,7 +1200,7 @@ bool PresolveTable(ConstraintProto* ct, PresolveContext* context) {
const int ref = ct->table().vars(j);
const int64 v = ct->table().values(i * num_vars + j);
tuple[j] = v;
if (!SortedDisjointIntervalsContain(context->GetRefDomain(ref), v)) {
if (!context->DomainOf(ref).Contains(v)) {
delete_row = true;
break;
}
@@ -1248,7 +1231,7 @@ bool PresolveTable(ConstraintProto* ct, PresolveContext* context) {
for (int j = 0; j < num_vars; ++j) {
const int ref = ct->table().vars(j);
changed |= context->IntersectDomainWith(
PositiveRef(ref), SortedDisjointIntervalsFromValues(std::vector<int64>(
PositiveRef(ref), Domain::FromValues(std::vector<int64>(
new_domains[j].begin(), new_domains[j].end())));
}
if (changed) {
@@ -1390,7 +1373,7 @@ bool PresolveCumulative(ConstraintProto* ct, PresolveContext* context) {
} else if (demand_max > capacity) {
if (ct.enforcement_literal().empty()) {
context->UpdateRuleStats("cumulative: demand_max exceeds capacity.");
context->IntersectDomainWith(demand_ref, {{kint64min, capacity}});
context->IntersectDomainWith(demand_ref, Domain(kint64min, capacity));
} else {
// TODO(user): we abort because we cannot convert this to a no_overlap
// for instance.
@@ -1989,8 +1972,9 @@ void PresolveCpModel(bool log_info, CpModelProto* presolved_model,
// it allows us to update the proto objective domain too.
CHECK_EQ(context.working_model->objective().vars_size(), 1);
CHECK_EQ(context.working_model->objective().coeffs(0), 1);
FillDomain(context.GetRefDomain(context.working_model->objective().vars(0)),
context.working_model->mutable_objective());
FillDomainInProto(
context.DomainOf(context.working_model->objective().vars(0)),
context.working_model->mutable_objective());
// We use a general code even if we currently have only one term (see checks
// above).
@@ -2108,25 +2092,22 @@ void PresolveCpModel(bool log_info, CpModelProto* presolved_model,
// If the objective variable wasn't used in other constraints and it can
// be reconstructed whatever the value of the other variables, we can
// remove the constraint. Note that we can't do that if the magnitude of
// objective_coeff_in_expanded_constraint is not one because of integer
// division.
// remove the constraint.
//
// TODO(user): It should be possible to refactor the code so this is
// automatically done by the linear constraint singleton presolve rule.
if (context.var_to_constraints[objective_var].size() == 1 &&
std::abs(objective_coeff_in_expanded_constraint) == 1) {
if (context.var_to_constraints[objective_var].size() == 1) {
// Compute implied domain on objective_var.
std::vector<ClosedInterval> rhs = {
{ct.linear().domain(0), ct.linear().domain(0)}};
Domain implied_domain = ReadDomainFromProto(ct.linear());
for (int i = 0; i < num_terms; ++i) {
const int ref = ct.linear().vars(i);
if (PositiveRef(ref) == objective_var) continue;
std::vector<ClosedInterval> domain = context.GetRefDomain(ref);
domain = MultiplicationOfSortedDisjointIntervals(
domain, -ct.linear().coeffs(i));
rhs = AdditionOfSortedDisjointIntervals(rhs, domain);
implied_domain = implied_domain.AdditionWith(
context.DomainOf(ref).ContinuousMultiplicationBy(
-ct.linear().coeffs(i)));
}
implied_domain = implied_domain.InverseMultiplicationBy(
objective_coeff_in_expanded_constraint);
// Remove the constraint if the implied domain is included in the
// domain of the objective_var term.
@@ -2134,15 +2115,8 @@ void PresolveCpModel(bool log_info, CpModelProto* presolved_model,
// Note the special case for the first expansion where any domain
// restriction will be handled by the objective domain because we
// called EncodeObjectiveAsSingleVariable() above.
std::vector<ClosedInterval> objective_var_domain =
context.GetRefDomain(objective_var);
objective_var_domain = MultiplicationOfSortedDisjointIntervals(
objective_var_domain, objective_coeff_in_expanded_constraint);
if (num_expansions == 0 ||
IntersectionOfSortedDisjointIntervals(
ComplementOfSortedDisjointIntervals(objective_var_domain),
rhs)
.empty()) {
implied_domain.IsIncludedIn(context.DomainOf(objective_var))) {
context.UpdateRuleStats("objective: removed objective constraint.");
*(context.mapping_model->add_constraints()) = ct;
context.working_model->mutable_constraints(expanded_linear_index)
@@ -2164,10 +2138,9 @@ void PresolveCpModel(bool log_info, CpModelProto* presolved_model,
}
mutable_objective->set_offset(mutable_objective->offset() +
objective_offset);
FillDomain(AdditionOfSortedDisjointIntervals(
ReadDomain(*mutable_objective),
{{-objective_offset, -objective_offset}}),
mutable_objective);
FillDomainInProto(ReadDomainFromProto(*mutable_objective)
.AdditionWith(Domain(-objective_offset)),
mutable_objective);
}
// Remove all empty or affine constraints (they will be re-added later if
@@ -2225,10 +2198,9 @@ void PresolveCpModel(bool log_info, CpModelProto* presolved_model,
CpModelProto* proto;
if (context.var_to_constraints[var].empty()) {
// Make sure that domain(representative) is tight.
const auto implied = InverseMultiplicationOfSortedDisjointIntervals(
AdditionOfSortedDisjointIntervals({{-r.offset, -r.offset}},
context.GetRefDomain(var)),
r.coeff);
const Domain implied = context.DomainOf(var)
.AdditionWith({-r.offset, -r.offset})
.InverseMultiplicationBy(r.coeff);
if (context.IntersectDomainWith(r.representative, implied)) {
LOG(WARNING) << "Domain of " << r.representative
<< " was not fully propagated using the affine relation "
@@ -2309,7 +2281,8 @@ void PresolveCpModel(bool log_info, CpModelProto* presolved_model,
// Update the variables domain of the presolved_model.
for (int i = 0; i < presolved_model->variables_size(); ++i) {
FillDomain(context.GetRefDomain(i), presolved_model->mutable_variables(i));
FillDomainInProto(context.DomainOf(i),
presolved_model->mutable_variables(i));
}
// Set the variables of the mapping_model.

View File

@@ -304,9 +304,9 @@ std::vector<int64> ValuesFromProto(const Values& values) {
}
// Returns the size of the given domain capped to int64max.
int64 DomainSize(const std::vector<ClosedInterval>& domain) {
int64 DomainSize(const Domain& domain) {
int64 size = 0;
for (const ClosedInterval interval : domain) {
for (const ClosedInterval interval : domain.intervals()) {
size += operations_research::CapAdd(
1, operations_research::CapSub(interval.end, interval.start));
}
@@ -318,6 +318,7 @@ int64 DomainSize(const std::vector<ClosedInterval>& domain) {
// never have any trivial inequalities.
void ModelWithMapping::ExtractEncoding(const CpModelProto& model_proto) {
IntegerEncoder* encoder = GetOrCreate<IntegerEncoder>();
IntegerTrail* integer_trail = GetOrCreate<IntegerTrail>();
// Detection of literal equivalent to (i_var == value). We collect all the
// half-reified constraint lit => equality or lit => inequality for a given
@@ -369,24 +370,27 @@ void ModelWithMapping::ExtractEncoding(const CpModelProto& model_proto) {
const sat::Literal enforcement_literal = Literal(ct.enforcement_literal(0));
const int ref = ct.linear().vars(0);
const int var = PositiveRef(ref);
const auto domain = ReadDomain(model_proto.variables(var));
const auto rhs = InverseMultiplicationOfSortedDisjointIntervals(
ReadDomain(ct.linear()),
ct.linear().coeffs(0) * (RefIsPositive(ref) ? 1 : -1));
const Domain domain = ReadDomainFromProto(model_proto.variables(var));
const Domain domain_if_enforced =
ReadDomainFromProto(ct.linear())
.InverseMultiplicationBy(ct.linear().coeffs(0) *
(RefIsPositive(ref) ? 1 : -1));
// Detect enforcement_literal => (var >= value or var <= value).
if (rhs.size() == 1) {
// We relax by 1 because we may take the negation of the rhs above.
if (rhs[0].end >= domain.back().end &&
rhs[0].start > domain.front().start) {
inequalities.push_back({&ct, enforcement_literal,
IntegerLiteral::GreaterOrEqual(
Integer(var), IntegerValue(rhs[0].start))});
} else if (rhs[0].start <= domain.front().start &&
rhs[0].end < domain.back().end) {
inequalities.push_back({&ct, enforcement_literal,
IntegerLiteral::LowerOrEqual(
Integer(var), IntegerValue(rhs[0].end))});
if (domain_if_enforced.intervals().size() == 1) {
if (domain_if_enforced.Max() >= domain.Max() &&
domain_if_enforced.Min() > domain.Min()) {
inequalities.push_back(
{&ct, enforcement_literal,
IntegerLiteral::GreaterOrEqual(
Integer(var), IntegerValue(domain_if_enforced.Min()))});
} else if (domain_if_enforced.Min() <= domain.Min() &&
domain_if_enforced.Max() < domain.Max()) {
inequalities.push_back(
{&ct, enforcement_literal,
IntegerLiteral::LowerOrEqual(
Integer(var), IntegerValue(domain_if_enforced.Max()))});
}
}
@@ -396,18 +400,18 @@ void ModelWithMapping::ExtractEncoding(const CpModelProto& model_proto) {
// and != 1. Similarly, for a domain in [min, max], we should both detect
// (== min) and (<= min), and both detect (== max) and (>= max).
{
const auto inter = IntersectionOfSortedDisjointIntervals(domain, rhs);
if (inter.size() == 1 && inter[0].start == inter[0].end) {
const Domain inter = domain.IntersectionWith(domain_if_enforced);
if (!inter.IsEmpty() && inter.Min() == inter.Max()) {
var_to_equalities[var].push_back(
{&ct, enforcement_literal, inter[0].start, true});
{&ct, enforcement_literal, inter.Min(), true});
}
}
{
const auto inter = IntersectionOfSortedDisjointIntervals(
domain, ComplementOfSortedDisjointIntervals(rhs));
if (inter.size() == 1 && inter[0].start == inter[0].end) {
const Domain inter =
domain.IntersectionWith(domain_if_enforced.Complement());
if (!inter.IsEmpty() && inter.Min() == inter.Max()) {
var_to_equalities[var].push_back(
{&ct, enforcement_literal, inter[0].start, false});
{&ct, enforcement_literal, inter.Min(), false});
}
}
}
@@ -419,6 +423,14 @@ void ModelWithMapping::ExtractEncoding(const CpModelProto& model_proto) {
if (inequalities[i].literal != inequalities[i + 1].literal.Negated()) {
continue;
}
if (inequalities[i].i_lit.bound <=
integer_trail->LowerBound(inequalities[i].i_lit.var)) {
continue;
}
if (inequalities[i + 1].i_lit.bound <=
integer_trail->LowerBound(inequalities[i + 1].i_lit.var)) {
continue;
}
const auto pair_a = encoder->Canonicalize(inequalities[i].i_lit);
const auto pair_b = encoder->Canonicalize(inequalities[i + 1].i_lit);
if (pair_a.first == pair_b.second) {
@@ -465,8 +477,7 @@ void ModelWithMapping::ExtractEncoding(const CpModelProto& model_proto) {
// Detect fully encoded variables and mark them as such.
//
// TODO(user): Also fully encode variable that are almost fully encoded.
const std::vector<ClosedInterval> domain =
ReadDomain(model_proto.variables(i));
const Domain domain = ReadDomainFromProto(model_proto.variables(i));
if (DomainSize(domain) == values.size()) {
++num_fully_encoded;
if (!encoder->VariableIsFullyEncoded(integers_[i])) {
@@ -580,7 +591,7 @@ ModelWithMapping::ModelWithMapping(const CpModelProto& model_proto,
integers_.resize(num_proto_variables, kNoIntegerVariable);
for (const int i : var_to_instantiate_as_integer) {
const auto& var_proto = model_proto.variables(i);
integers_[i] = Add(NewIntegerVariable(ReadDomain(var_proto)));
integers_[i] = Add(NewIntegerVariable(ReadDomainFromProto(var_proto)));
if (integers_[i] >= reverse_integer_map_.size()) {
reverse_integer_map_.resize(integers_[i].value() + 1, -1);
}
@@ -1108,7 +1119,7 @@ void DetectEquivalencesInElementConstraint(const ConstraintProto& ct,
if (m->Get(IsFixed(index))) return;
std::vector<ClosedInterval> union_of_non_constant_domains;
Domain union_of_non_constant_domains;
std::map<IntegerValue, int> constant_to_num;
for (const auto literal_value : m->Add(FullyEncodeVariable(index))) {
const int i = literal_value.value.value();
@@ -1116,16 +1127,14 @@ void DetectEquivalencesInElementConstraint(const ConstraintProto& ct,
const IntegerValue value(m->Get(Value(vars[i])));
constant_to_num[value]++;
} else {
union_of_non_constant_domains = UnionOfSortedDisjointIntervals(
union_of_non_constant_domains,
union_of_non_constant_domains = union_of_non_constant_domains.UnionWith(
integer_trail->InitialVariableDomain(vars[i]));
}
}
// Bump the number if the constant appear in union_of_non_constant_domains.
for (const auto entry : constant_to_num) {
if (SortedDisjointIntervalsContain(union_of_non_constant_domains,
entry.first.value())) {
if (union_of_non_constant_domains.Contains(entry.first.value())) {
constant_to_num[entry.first]++;
}
}
@@ -1496,14 +1505,13 @@ std::string CpModelStats(const CpModelProto& model_proto) {
int num_constants = 0;
std::set<int64> constant_values;
std::map<std::vector<ClosedInterval>, int, ExactVectorOfDomainComparator>
num_vars_per_domains;
std::map<Domain, int> num_vars_per_domains;
for (const IntegerVariableProto& var : model_proto.variables()) {
if (var.domain_size() == 2 && var.domain(0) == var.domain(1)) {
++num_constants;
constant_values.insert(var.domain(0));
} else {
num_vars_per_domains[ReadDomain(var)]++;
num_vars_per_domains[ReadDomainFromProto(var)]++;
}
}
@@ -1537,8 +1545,8 @@ std::string CpModelStats(const CpModelProto& model_proto) {
objective_string.c_str(), "\n");
if (num_vars_per_domains.size() < 50) {
for (const auto& entry : num_vars_per_domains) {
const std::string temp = absl::StrCat(
" - ", entry.second, " in ", IntervalsAsString(entry.first), "\n");
const std::string temp = absl::StrCat(" - ", entry.second, " in ",
entry.first.ToString().c_str(), "\n");
absl::StrAppend(&result, Summarize(temp));
}
} else {
@@ -1546,10 +1554,10 @@ std::string CpModelStats(const CpModelProto& model_proto) {
int64 min = kint64max;
int64 max = kint64min;
for (const auto& entry : num_vars_per_domains) {
min = std::min(min, entry.first.front().start);
max = std::max(max, entry.first.back().end);
max_complexity =
std::max(max_complexity, static_cast<int64>(entry.first.size()));
min = std::min(min, entry.first.Min());
max = std::max(max, entry.first.Max());
max_complexity = std::max(
max_complexity, static_cast<int64>(entry.first.intervals().size()));
}
absl::StrAppend(&result, " - ", num_vars_per_domains.size(),
" different domains in [", min, ",", max,
@@ -2517,8 +2525,8 @@ CpSolverResponse SolveCpModelInternal(
// Intersect the objective domain with the given one if any.
if (!model_proto.objective().domain().empty()) {
const auto user_domain = ReadDomain(model_proto.objective());
const auto automatic_domain =
const Domain user_domain = ReadDomainFromProto(model_proto.objective());
const Domain automatic_domain =
model->GetOrCreate<IntegerTrail>()->InitialVariableDomain(
objective_var);
VLOG(2) << "Objective offset:" << model_proto.objective().offset()
@@ -2741,10 +2749,10 @@ void PostsolveResponse(const CpModelProto& model_proto,
}
for (int i = 0; i < response->solution_lower_bounds_size(); ++i) {
auto* var_proto = mapping_proto.mutable_variables(postsolve_mapping[i]);
FillDomain(
IntersectionOfSortedDisjointIntervals(
ReadDomain(*var_proto), {{response->solution_lower_bounds(i),
response->solution_upper_bounds(i)}}),
FillDomainInProto(
ReadDomainFromProto(*var_proto)
.IntersectionWith({response->solution_lower_bounds(i),
response->solution_upper_bounds(i)}),
var_proto);
}
@@ -3297,8 +3305,8 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) {
if (!FLAGS_cp_model_dump_file.empty()) {
LOG(INFO) << "Dumping cp model proto to '" << FLAGS_cp_model_dump_file
<< "'.";
CHECK_OK(file::SetBinaryProto(FLAGS_cp_model_dump_file, model_proto,
file::Defaults()));
CHECK_OK(file::SetTextProto(FLAGS_cp_model_dump_file, model_proto,
file::Defaults()));
}
// Override parameters?

View File

@@ -92,6 +92,10 @@ void FillDomain(const std::vector<ClosedInterval>& domain,
proto->add_domain(interval.end);
}
}
template <typename ProtoWithDomain>
void FillDomainInProto(const Domain& domain, ProtoWithDomain* proto) {
FillDomain(domain.intervals(), proto);
}
// Extract a sorted interval list from the domain field of a proto.
template <typename ProtoWithDomain>
@@ -103,6 +107,10 @@ std::vector<ClosedInterval> ReadDomain(const ProtoWithDomain& proto) {
CHECK(IntervalsAreSortedAndDisjoint(result));
return result;
}
template <typename ProtoWithDomain>
Domain ReadDomainFromProto(const ProtoWithDomain& proto) {
return Domain::FromIntervals(ReadDomain(proto));
}
// Returns the list of values in a given domain.
// This will fail if the domain contains more than one millions values.

View File

@@ -225,8 +225,7 @@ bool DisjunctiveWithTwoItems::Propagate() {
helper_->StartMin(task_after) < helper_->EndMin(task_before)) {
// Reason for precedences if both present.
helper_->ClearReason();
helper_->AddStartMaxReason(task_before, helper_->EndMin(task_after) - 1);
helper_->AddEndMinReason(task_after, helper_->EndMin(task_after));
helper_->AddReasonForBeingBefore(task_before, task_after);
// Reason for the bound push.
helper_->AddPresenceReason(task_before);
@@ -240,8 +239,7 @@ bool DisjunctiveWithTwoItems::Propagate() {
helper_->EndMax(task_before) > helper_->StartMax(task_after)) {
// Reason for precedences if both present.
helper_->ClearReason();
helper_->AddStartMaxReason(task_before, helper_->EndMin(task_after) - 1);
helper_->AddEndMinReason(task_after, helper_->EndMin(task_after));
helper_->AddReasonForBeingBefore(task_before, task_after);
// Reason for the bound push.
helper_->AddPresenceReason(task_after);

View File

@@ -314,8 +314,8 @@ void IntegerEncoder::AssociateToIntegerEqualValue(Literal literal,
// Detect literal view. Note that the same literal can be associated to more
// than one variable, and thus already have a view. We don't change it in
// this case.
const auto& domain = (*domains_)[var];
if (value == 1 && domain.front().start >= 0 && domain.back().end <= 1) {
const Domain domain = Domain::FromIntervals((*domains_)[var]);
if (value == 1 && domain.Min() >= 0 && domain.Max() <= 1) {
if (literal.Index() >= literal_view_.size()) {
literal_view_.resize(literal.Index().value() + 1, kNoIntegerVariable);
literal_view_[literal.Index()] = var;
@@ -323,7 +323,7 @@ void IntegerEncoder::AssociateToIntegerEqualValue(Literal literal,
literal_view_[literal.Index()] = var;
}
}
if (value == -1 && domain.front().start >= -1 && domain.back().end <= 0) {
if (value == -1 && domain.Min() >= -1 && domain.Max() <= 0) {
if (literal.Index() >= literal_view_.size()) {
literal_view_.resize(literal.Index().value() + 1, kNoIntegerVariable);
literal_view_[literal.Index()] = NegationOf(var);
@@ -333,17 +333,17 @@ void IntegerEncoder::AssociateToIntegerEqualValue(Literal literal,
}
// Fix literal for value outside the domain or for singleton domain.
if (!SortedDisjointIntervalsContain(domain, value.value())) {
if (!domain.Contains(value.value())) {
sat_solver_->AddUnitClause(literal.Negated());
return;
}
if (value == domain.front().start && value == domain.back().end) {
if (value == domain.Min() && value == domain.Max()) {
sat_solver_->AddUnitClause(literal);
return;
}
// Special case for the first and last value.
if (value == domain.front().start) {
if (value == domain.Min()) {
// Note that this will recursively call AssociateToIntegerEqualValue() but
// since equality_to_associated_literal_[] is now set, the recursion will
// stop there. When a domain has just 2 values, this allows to call just
@@ -353,7 +353,7 @@ void IntegerEncoder::AssociateToIntegerEqualValue(Literal literal,
IntegerLiteral::LowerOrEqual(var, value));
return;
}
if (value == domain.back().end) {
if (value == domain.Max()) {
AssociateToIntegerLiteral(literal,
IntegerLiteral::GreaterOrEqual(var, value));
return;
@@ -537,12 +537,10 @@ IntegerVariable IntegerTrail::AddIntegerVariable(IntegerValue lower_bound,
return i;
}
IntegerVariable IntegerTrail::AddIntegerVariable(
const std::vector<ClosedInterval>& domain) {
CHECK(!domain.empty());
CHECK(IntervalsAreSortedAndDisjoint(domain));
const IntegerVariable var = AddIntegerVariable(
IntegerValue(domain.front().start), IntegerValue(domain.back().end));
IntegerVariable IntegerTrail::AddIntegerVariable(const Domain& domain) {
CHECK(!domain.IsEmpty());
const IntegerVariable var = AddIntegerVariable(IntegerValue(domain.Min()),
IntegerValue(domain.Max()));
CHECK(UpdateInitialDomain(var, domain));
return var;
}
@@ -551,32 +549,30 @@ IntegerVariable IntegerTrail::AddIntegerVariable(
//
// TODO(user): we could return a reference, but this is likely not in any
// critical code path.
std::vector<ClosedInterval> IntegerTrail::InitialVariableDomain(
IntegerVariable var) const {
std::vector<ClosedInterval> domain;
for (const ClosedInterval& i : (*domains_)[var]) domain.push_back(i);
return domain;
Domain IntegerTrail::InitialVariableDomain(IntegerVariable var) const {
return Domain::FromIntervals((*domains_)[var]);
}
bool IntegerTrail::UpdateInitialDomain(IntegerVariable var,
std::vector<ClosedInterval> domain) {
bool IntegerTrail::UpdateInitialDomain(IntegerVariable var, Domain domain) {
CHECK_EQ(trail_->CurrentDecisionLevel(), 0);
// TODO(user): A bit inefficient as this recreate a vector for no reason.
// The IntersectionOfSortedDisjointIntervals() should take a Span<> instead.
std::vector<ClosedInterval> old_domain = InitialVariableDomain(var);
domain = IntersectionOfSortedDisjointIntervals(domain, old_domain);
const Domain old_domain = InitialVariableDomain(var);
domain = domain.IntersectionWith(old_domain);
if (old_domain == domain) return true;
if (domain.empty()) return false;
if (domain.IsEmpty()) return false;
(*domains_)[var].assign(domain.begin(), domain.end());
{
std::vector<ClosedInterval> temp =
NegationOfSortedDisjointIntervals(domain);
const std::vector<ClosedInterval> temp = domain.intervals();
(*domains_)[var].assign(temp.begin(), temp.end());
}
{
const std::vector<ClosedInterval> temp = domain.Negation().intervals();
(*domains_)[NegationOf(var)].assign(temp.begin(), temp.end());
}
if (domain.size() > 1) {
if (domain.intervals().size() > 1) {
var_to_current_lb_interval_index_.Set(var, 0);
var_to_current_lb_interval_index_.Set(NegationOf(var), 0);
}
@@ -585,21 +581,20 @@ bool IntegerTrail::UpdateInitialDomain(IntegerVariable var,
// bounds here directly. This is because these function might call again
// UpdateInitialDomain(), and we will abort after realizing that the domain
// didn't change this time.
CHECK(Enqueue(
IntegerLiteral::GreaterOrEqual(var, IntegerValue(domain.front().start)),
{}, {}));
CHECK(Enqueue(
IntegerLiteral::LowerOrEqual(var, IntegerValue(domain.back().end)), {},
{}));
CHECK(Enqueue(IntegerLiteral::GreaterOrEqual(var, IntegerValue(domain.Min())),
{}, {}));
CHECK(Enqueue(IntegerLiteral::LowerOrEqual(var, IntegerValue(domain.Max())),
{}, {}));
// If the variable is fully encoded, set to false excluded literals.
if (encoder_->VariableIsFullyEncoded(var)) {
int i = 0;
int num_fixed = 0;
const auto encoding = encoder_->FullDomainEncoding(var);
const auto intervals = domain.intervals();
for (const auto pair : encoding) {
while (i < domain.size() && pair.value > domain[i].end) ++i;
if (i == domain.size() || pair.value < domain[i].start) {
while (i < intervals.size() && pair.value > intervals[i].end) ++i;
if (i == intervals.size() || pair.value < intervals[i].start) {
// Set the literal to false;
++num_fixed;
if (trail_->Assignment().LiteralIsTrue(pair.literal)) return false;
@@ -674,6 +669,44 @@ int IntegerTrail::FindLowestTrailIndexThatExplainBound(
}
}
// We try to relax the reason in a smart way here by minimizing the maximum
// trail indices of the literals appearing in reason.
//
// TODO(user): use priority queue instead of O(n^2) algo.
void IntegerTrail::RelaxLinearReason(
IntegerValue slack, absl::Span<IntegerValue> coeffs,
std::vector<IntegerLiteral>* reason) const {
CHECK_GE(slack, 0);
if (slack == 0) return;
const int size = reason->size();
std::vector<int> indices(size);
for (int i = 0; i < size; ++i) {
CHECK_EQ((*reason)[i].bound, LowerBound((*reason)[i].var));
CHECK_GT(coeffs[i], 0);
indices[i] = vars_[(*reason)[i].var].current_trail_index;
}
const int num_vars = vars_.size();
while (slack != 0) {
int best_i = -1;
for (int i = 0; i < size; ++i) {
if (indices[i] < num_vars) continue; // level zero.
if (best_i != -1 && indices[i] < indices[best_i]) continue;
const TrailEntry& entry = integer_trail_[indices[i]];
const TrailEntry& previous_entry = integer_trail_[entry.prev_trail_index];
if (coeffs[i] * (entry.bound - previous_entry.bound) > slack) continue;
best_i = i;
}
if (best_i == -1) return;
const TrailEntry& entry = integer_trail_[indices[best_i]];
const TrailEntry& previous_entry = integer_trail_[entry.prev_trail_index];
indices[best_i] = entry.prev_trail_index;
(*reason)[best_i].bound = previous_entry.bound;
slack -= coeffs[best_i] * (entry.bound - previous_entry.bound);
}
}
bool IntegerTrail::EnqueueAssociatedLiteral(
Literal literal, int trail_index_with_same_reason,
absl::Span<Literal> literal_reason,
@@ -923,8 +956,9 @@ bool IntegerTrail::Enqueue(IntegerLiteral i_lit,
integer_trail_[i_lit.var.value()].bound = i_lit.bound;
// We also update the initial domain.
CHECK(UpdateInitialDomain(i_lit.var, {{LowerBound(i_lit.var).value(),
UpperBound(i_lit.var).value()}}));
CHECK(UpdateInitialDomain(
i_lit.var,
Domain(LowerBound(i_lit.var).value(), UpperBound(i_lit.var).value())));
return true;
}
DCHECK_GT(trail_->CurrentDecisionLevel(), 0);

View File

@@ -447,29 +447,19 @@ class IntegerTrail : public SatPropagator {
IntegerValue upper_bound);
// Same as above but for a more complex domain specified as a sorted list of
// disjoint intervals. Note that the ClosedInterval struct use int64 instead
// of integer values (but we will convert them internally).
//
// Precondition: we check that IntervalsAreSortedAndDisjoint(domain) is true.
IntegerVariable AddIntegerVariable(const std::vector<ClosedInterval>& domain);
// disjoint intervals. See the Domain class.
IntegerVariable AddIntegerVariable(const Domain& domain);
// Returns the initial domain of the given variable. Note that for variables
// whose domain is a single interval, this is updated with level zero
// propagations, but not if the domain is more complex.
std::vector<ClosedInterval> InitialVariableDomain(IntegerVariable var) const;
// Returns the initial domain of the given variable. Note that the min/max
// are updated with level zero propagation, but not holes.
Domain InitialVariableDomain(IntegerVariable var) const;
// Takes the intersection with the current initial variable domain.
//
// TODO(user): There is some memory inefficiency if this is called many time
// because of the underlying data structure we use. In practice, when used
// with a presolve, this is not often used, so that is fine though.
//
// TODO(user): The Enqueue() done at level zero on a variable are not
// reflected on its initial domain. That can causes issue if the variable
// is fully encoded afterwards because literals will be created for the values
// no longer relevant, and these will not be propagated right away.
bool UpdateInitialDomain(IntegerVariable var,
std::vector<ClosedInterval> domain);
bool UpdateInitialDomain(IntegerVariable var, Domain domain);
// Same as AddIntegerVariable(value, value), but this is a bit more efficient
// because it reuses another constant with the same value if its exist.
@@ -532,6 +522,25 @@ class IntegerTrail : public SatPropagator {
IntegerLiteral LowerBoundAsLiteral(IntegerVariable i) const;
IntegerLiteral UpperBoundAsLiteral(IntegerVariable i) const;
// Advanced usage. Given the reason for
// (Sum_i coeffs[i] * reason[i].var >= current_lb) initially in reason,
// this function relaxes the reason given that we only need the explanation of
// (Sum_i coeffs[i] * reason[i].var >= current_lb - slack).
//
// Preconditions:
// - coeffs must be of same size as reason, and all entry must be positive.
// - *reason must initially contains the trivial initial reason, that is
// the current lower-bound of each variables.
//
// TODO(user): Requiring all initial literal to be at their current bound is
// not really clean. Maybe we can change the API to only take IntegerVariable
// and produce the reason directly.
//
// TODO(user): change API so that this work is performed during the conflict
// analysis. Note that we could be smarter there.
void RelaxLinearReason(IntegerValue slack, absl::Span<IntegerValue> coeffs,
std::vector<IntegerLiteral>* reason) const;
// Enqueue new information about a variable bound. Calling this with a less
// restrictive bound than the current one will have no effect.
//
@@ -1023,7 +1032,7 @@ inline std::function<IntegerVariable(Model*)> NewIntegerVariable(int64 lb,
}
inline std::function<IntegerVariable(Model*)> NewIntegerVariable(
const std::vector<ClosedInterval>& domain) {
const Domain& domain) {
return [=](Model* model) {
return model->GetOrCreate<IntegerTrail>()->AddIntegerVariable(domain);
};

View File

@@ -111,27 +111,15 @@ bool IntegerSumLE::Propagate() {
const IntegerValue new_lb = rev_lb_fixed_vars_ + lb_unfixed_vars;
// Conflict?
IntegerValue slack = upper_bound_ - new_lb;
const IntegerValue slack = upper_bound_ - new_lb;
if (slack < 0) {
// Like FillIntegerReason() but try to relax the reason a bit.
//
// TODO(user): if not all the slack is consumed, we could relax it even
// more. It might also be advantageous to relax first the variable with the
// highest "trail index".
integer_reason_.clear();
FillIntegerReason();
std::vector<IntegerValue> coeffs;
for (int i = 0; i < vars_.size(); ++i) {
const IntegerVariable var = vars_[i];
const IntegerValue lb = integer_trail_->LowerBound(var);
const IntegerValue prev_lb = integer_trail_->PreviousLowerBound(var);
if (lb == prev_lb) continue; // level zero.
const IntegerValue diff = (lb - prev_lb) * coeffs_[i];
if (slack + diff < 0) {
integer_reason_.push_back(IntegerLiteral::GreaterOrEqual(var, prev_lb));
slack += diff;
} else {
integer_reason_.push_back(IntegerLiteral::GreaterOrEqual(var, lb));
}
if (index_in_integer_reason_[i] == -1) continue;
coeffs.push_back(coeffs_[i]);
}
integer_trail_->RelaxLinearReason(-slack - 1, coeffs, &integer_reason_);
if (num_unassigned_enforcement_literal == 1) {
// Propagate the only non-true literal to false.
@@ -550,8 +538,7 @@ std::function<void(Model*)> IsOneOf(IntegerVariable var,
}
gtl::STLSortAndRemoveDuplicates(&unique_values);
integer_trail->UpdateInitialDomain(
var, SortedDisjointIntervalsFromValues(unique_values));
integer_trail->UpdateInitialDomain(var, Domain::FromValues(unique_values));
if (unique_values.size() == 1) {
model->Add(ClauseConstraint(selectors));
return;

View File

@@ -154,6 +154,20 @@ SchedulingConstraintHelper::TaskByDecreasingEndMax() {
return task_by_decreasing_max_end_;
}
// Produces a relaxed reason for StartMax(before) < EndMin(after).
void SchedulingConstraintHelper::AddReasonForBeingBefore(int before,
int after) {
DCHECK_LT(StartMax(before), EndMin(after));
const IntegerValue slack = EndMin(after) - StartMax(before) - 1;
std::vector<IntegerLiteral> temp;
temp.push_back(integer_trail_->LowerBoundAsLiteral(end_vars_[after]));
temp.push_back(
integer_trail_->LowerBoundAsLiteral(NegationOf(start_vars_[before])));
integer_trail_->RelaxLinearReason(slack, {IntegerValue(1), IntegerValue(1)},
&temp);
integer_reason_.insert(integer_reason_.end(), temp.begin(), temp.end());
}
bool SchedulingConstraintHelper::PushIntegerLiteral(IntegerLiteral bound) {
return integer_trail_->Enqueue(bound, literal_reason_, integer_reason_);
}

View File

@@ -179,6 +179,10 @@ class SchedulingConstraintHelper {
void AddEndMinReason(int t, IntegerValue lower_bound);
void AddEndMaxReason(int t, IntegerValue upper_bound);
// Adds the reason why task "before" must be before task "after".
// That is StartMax(before) < EndMin(after).
void AddReasonForBeingBefore(int before, int after);
// It is also possible to directly manipulates the underlying reason vectors
// that will be used when pushing something.
std::vector<Literal>* MutableLiteralReason() { return &literal_reason_; }

View File

@@ -452,6 +452,37 @@ bool PrecedencesPropagator::EnqueueAndCheck(const ArcInfo& arc,
AppendLowerBoundReasonIfValid(arc.offset_var, *integer_trail_,
&integer_reason_);
// The code works without this block since Enqueue() below can already take
// care of conflicts. However, it is better to deal with the conflict
// ourselves because we can be smarter about the reason this way.
//
// The reason for a "precedence" conflict is always a linear reason
// involving the tail lower_bound, the head upper bound and eventually the
// size lower bound. Because of that, we can use the RelaxLinearReason()
// code.
if (new_head_lb > integer_trail_->UpperBound(arc.head_var)) {
const IntegerValue slack =
new_head_lb - integer_trail_->UpperBound(arc.head_var) - 1;
integer_reason_.push_back(
integer_trail_->UpperBoundAsLiteral(arc.head_var));
std::vector<IntegerValue> coeffs(integer_reason_.size(), IntegerValue(1));
integer_trail_->RelaxLinearReason(slack, coeffs, &integer_reason_);
if (!integer_trail_->IsOptional(arc.head_var)) {
return integer_trail_->ReportConflict(literal_reason_, integer_reason_);
} else {
CHECK(!integer_trail_->IsCurrentlyIgnored(arc.head_var));
const Literal l = integer_trail_->IsIgnoredLiteral(arc.head_var);
if (trail->Assignment().LiteralIsFalse(l)) {
literal_reason_.push_back(l);
return integer_trail_->ReportConflict(literal_reason_, integer_reason_);
} else {
integer_trail_->EnqueueLiteral(l, literal_reason_, integer_reason_);
return true;
}
}
}
return integer_trail_->Enqueue(
IntegerLiteral::GreaterOrEqual(arc.head_var, new_head_lb),
literal_reason_, integer_reason_);

View File

@@ -463,7 +463,6 @@ class Constraint(object):
class IntervalVar(object):
"""Represents a Interval variable.
An interval variable is both a constraint and a variable. It is defined by
three integer variables: start, size, and end.

View File

@@ -61,13 +61,12 @@ std::unordered_map<IntegerValue, Literal> GetEncoding(IntegerVariable var,
void FilterValues(IntegerVariable var, Model* model,
std::unordered_set<int64>* values) {
std::vector<ClosedInterval> domain =
model->Get<IntegerTrail>()->InitialVariableDomain(var);
const Domain domain = model->Get<IntegerTrail>()->InitialVariableDomain(var);
for (auto it = values->begin(); it != values->end();) {
const int64 v = *it;
auto copy = it++;
// TODO(user): quadratic! improve.
if (!SortedDisjointIntervalsContain(domain, v)) {
if (!domain.Contains(v)) {
values->erase(copy);
}
}
@@ -171,8 +170,8 @@ std::function<void(Model*)> TableConstraint(
[first](int64 v) { return v == first; })) {
model->Add(Equality(vars[i], first));
} else {
integer_trail->UpdateInitialDomain(
vars[i], SortedDisjointIntervalsFromValues(tr_tuples[i]));
integer_trail->UpdateInitialDomain(vars[i],
Domain::FromValues(tr_tuples[i]));
model->Add(FullyEncodeVariable(vars[i]));
ProcessOneColumn(
tuple_literals,
@@ -322,7 +321,7 @@ std::function<void(Model*)> TransitionConstraint(
const auto domain = integer_trail->InitialVariableDomain(vars[time]);
for (const std::vector<int64>& transition : automata) {
// TODO(user): quadratic algo, improve!
if (SortedDisjointIntervalsContain(domain, transition[1])) {
if (domain.Contains(transition[1])) {
possible_values[time].insert(transition[1]);
}
}
@@ -411,8 +410,8 @@ std::function<void(Model*)> TransitionConstraint(
std::vector<int64> values;
values.reserve(s.size());
for (IntegerValue v : s) values.push_back(v.value());
integer_trail->UpdateInitialDomain(
vars[time], SortedDisjointIntervalsFromValues(values));
integer_trail->UpdateInitialDomain(vars[time],
Domain::FromValues(values));
model->Add(FullyEncodeVariable(vars[time]));
encoding = GetEncoding(vars[time], model);
} else {

View File

@@ -26,26 +26,6 @@ std::string ClosedInterval::DebugString() const {
return absl::StrFormat("[%" GG_LL_FORMAT "d,%" GG_LL_FORMAT "d]", start, end);
}
bool ExactDomainComparator::operator()(const ClosedInterval& i1,
const ClosedInterval& i2) const {
return i1.start < i2.start || (i1.start == i2.start && i1.end < i2.end);
}
bool ExactVectorOfDomainComparator::operator()(
const std::vector<ClosedInterval>& d1,
const std::vector<ClosedInterval>& d2) const {
const int common_size = std::min(d1.size(), d2.size());
for (int i = 0; i < common_size; ++i) {
const ClosedInterval& i1 = d1[i];
const ClosedInterval& i2 = d2[i];
if (i1.start < i2.start) return true;
if (i1.start > i2.start) return false;
if (i1.end < i2.end) return true;
if (i1.end > i2.end) return false;
}
return d1.size() < d2.size();
}
std::string IntervalsAsString(const std::vector<ClosedInterval>& intervals) {
std::string result;
for (ClosedInterval interval : intervals) {
@@ -54,13 +34,13 @@ std::string IntervalsAsString(const std::vector<ClosedInterval>& intervals) {
return result;
}
std::ostream& operator<<(std::ostream& out, const ClosedInterval& interval) {
return out << interval.DebugString();
}
std::ostream& operator<<(std::ostream& out,
const std::vector<ClosedInterval>& intervals) {
return out << IntervalsAsString(intervals);
// TODO(user): binary search if size is large?
bool SortedDisjointIntervalsContain(absl::Span<ClosedInterval> intervals,
int64 value) {
for (const ClosedInterval& interval : intervals) {
if (interval.start <= value && interval.end >= value) return true;
}
return false;
}
std::vector<ClosedInterval> SortedDisjointIntervalsFromValues(
@@ -95,14 +75,27 @@ bool IntervalsAreSortedAndDisjoint(
return true;
}
bool SortedDisjointIntervalsContain(absl::Span<ClosedInterval> intervals,
int64 value) {
for (const ClosedInterval& interval : intervals) {
if (interval.start <= value && interval.end >= value) return true;
namespace {
// Transforms a sorted list of intervals in a sorted DISJOINT list for which
// IntervalsAreSortedAndDisjoint() would return true.
void UnionOfSortedIntervals(std::vector<ClosedInterval>* intervals) {
DCHECK(std::is_sorted(intervals->begin(), intervals->end()));
int new_size = 0;
for (const ClosedInterval& i : *intervals) {
if (new_size > 0 && i.start <= CapAdd((*intervals)[new_size - 1].end, 1)) {
(*intervals)[new_size - 1].end =
std::max(i.end, (*intervals)[new_size - 1].end);
} else {
(*intervals)[new_size++] = i;
}
}
return false;
intervals->resize(new_size);
DCHECK(IntervalsAreSortedAndDisjoint(*intervals)) << *intervals;
}
} // namespace
std::vector<ClosedInterval> IntersectionOfSortedDisjointIntervals(
const std::vector<ClosedInterval>& a,
const std::vector<ClosedInterval>& b) {
@@ -163,27 +156,6 @@ std::vector<ClosedInterval> NegationOfSortedDisjointIntervals(
return intervals;
}
namespace {
// Transforms a sorted list of intervals in a sorted DISJOINT list for which
// IntervalsAreSortedAndDisjoint() would return true.
void UnionOfSortedIntervals(std::vector<ClosedInterval>* intervals) {
DCHECK(std::is_sorted(intervals->begin(), intervals->end()));
int new_size = 0;
for (const ClosedInterval& i : *intervals) {
if (new_size > 0 && i.start <= CapAdd((*intervals)[new_size - 1].end, 1)) {
(*intervals)[new_size - 1].end =
std::max(i.end, (*intervals)[new_size - 1].end);
} else {
(*intervals)[new_size++] = i;
}
}
intervals->resize(new_size);
DCHECK(IntervalsAreSortedAndDisjoint(*intervals)) << *intervals;
}
} // namespace
std::vector<ClosedInterval> UnionOfSortedDisjointIntervals(
const std::vector<ClosedInterval>& a,
const std::vector<ClosedInterval>& b) {
@@ -297,6 +269,151 @@ std::vector<ClosedInterval> DivisionOfSortedDisjointIntervals(
return coeff > 0 ? intervals : NegationOfSortedDisjointIntervals(intervals);
}
std::ostream& operator<<(std::ostream& out, const ClosedInterval& interval) {
return out << interval.DebugString();
}
std::ostream& operator<<(std::ostream& out,
const std::vector<ClosedInterval>& intervals) {
return out << IntervalsAsString(intervals);
}
std::ostream& operator<<(std::ostream& out, const Domain& domain) {
return out << IntervalsAsString(domain.intervals());
}
Domain::Domain(int64 value) : intervals_({{value, value}}) {}
Domain::Domain(int64 left, int64 right) : intervals_({{left, right}}) {
if (left > right) intervals_.clear();
}
Domain Domain::AllValues() { return Domain(kint64min, kint64max); }
Domain Domain::FromValues(absl::Span<int64> values) {
Domain result;
result.intervals_ = SortedDisjointIntervalsFromValues(
std::vector<int64>(values.begin(), values.end()));
return result;
}
Domain Domain::FromIntervals(absl::Span<ClosedInterval> intervals) {
Domain result;
result.intervals_.assign(intervals.begin(), intervals.end());
std::sort(result.intervals_.begin(), result.intervals_.end());
UnionOfSortedIntervals(&result.intervals_);
return result;
}
bool Domain::IsEmpty() const { return intervals_.empty(); }
int64 Domain::Min() const {
CHECK(!IsEmpty());
return intervals_.front().start;
}
int64 Domain::Max() const {
CHECK(!IsEmpty());
return intervals_.back().end;
}
// TODO(user): In all the Domain::Function() remove the indirection and
// optimize.
bool Domain::Contains(int64 value) const {
return SortedDisjointIntervalsContain(intervals_, value);
}
bool Domain::IsIncludedIn(const Domain& domain) const {
int i = 0;
const std::vector<ClosedInterval>& others = domain.intervals_;
for (const ClosedInterval interval : intervals_) {
// Find the unique interval in others that contains interval if any.
for (; i < others.size() && interval.end > others[i].end; ++i) {
}
if (i == others.size()) return false;
if (interval.start < others[i].start) return false;
}
return true;
}
Domain Domain::Complement() const {
Domain result;
result.intervals_ = ComplementOfSortedDisjointIntervals(intervals_);
return result;
}
Domain Domain::Negation() const {
Domain result;
result.intervals_ = NegationOfSortedDisjointIntervals(intervals_);
return result;
}
Domain Domain::IntersectionWith(const Domain& domain) const {
Domain result;
result.intervals_ =
IntersectionOfSortedDisjointIntervals(intervals_, domain.intervals_);
return result;
}
Domain Domain::UnionWith(const Domain& domain) const {
Domain result;
result.intervals_ =
UnionOfSortedDisjointIntervals(intervals_, domain.intervals_);
return result;
}
Domain Domain::AdditionWith(const Domain& domain) const {
Domain result;
result.intervals_ =
AdditionOfSortedDisjointIntervals(intervals_, domain.intervals_);
return result;
}
Domain Domain::MultiplicationBy(int64 coeff, bool* success) const {
Domain result;
result.intervals_ = PreciseMultiplicationOfSortedDisjointIntervals(
intervals_, coeff, success);
return result;
}
Domain Domain::ContinuousMultiplicationBy(int64 coeff) const {
Domain result;
result.intervals_ =
MultiplicationOfSortedDisjointIntervals(intervals_, coeff);
return result;
}
Domain Domain::DivisionBy(int64 coeff) const {
Domain result;
result.intervals_ = DivisionOfSortedDisjointIntervals(intervals_, coeff);
return result;
}
Domain Domain::InverseMultiplicationBy(const int64 coeff) const {
Domain result;
result.intervals_ =
InverseMultiplicationOfSortedDisjointIntervals(intervals_, coeff);
return result;
}
bool Domain::operator<(const Domain& other) const {
const std::vector<ClosedInterval>& d1 = intervals_;
const std::vector<ClosedInterval>& d2 = other.intervals_;
const int common_size = std::min(d1.size(), d2.size());
for (int i = 0; i < common_size; ++i) {
const ClosedInterval& i1 = d1[i];
const ClosedInterval& i2 = d2[i];
if (i1.start < i2.start) return true;
if (i1.start > i2.start) return false;
if (i1.end < i2.end) return true;
if (i1.end > i2.end) return false;
}
return d1.size() < d2.size();
}
std::string Domain::ToString() const { return IntervalsAsString(intervals_); }
SortedDisjointIntervalList::SortedDisjointIntervalList() {}
SortedDisjointIntervalList::SortedDisjointIntervalList(

View File

@@ -42,101 +42,143 @@ struct ClosedInterval {
}
};
// Custom exact comparators.
class ExactDomainComparator {
public:
bool operator()(const ClosedInterval& i1, const ClosedInterval& i2) const;
};
class ExactVectorOfDomainComparator {
public:
bool operator()(const std::vector<ClosedInterval>& d1,
const std::vector<ClosedInterval>& d2) const;
};
// Returns a compact std::string of a vector of intervals like
// "[1,4][6][10,20]".
std::string IntervalsAsString(const std::vector<ClosedInterval>& intervals);
std::ostream& operator<<(std::ostream& out, const ClosedInterval& interval);
std::ostream& operator<<(std::ostream& out,
const std::vector<ClosedInterval>& intervals);
// TODO(user): Regroup all the functions below in a SortedDisjointIntervalVector
// class, it will lead to shorter/easier to use names. This will also allow to
// use an inlined vector for the most common case of one or two intervals.
// Converts an unsorted list of integer values to the unique list of
// non-adjacent, disjoint ClosedInterval spanning exactly these values. Eg. for
// values (after sorting): {1, 2, 3, 5, 7, 8, 10, 11}, it returns the list of
// intervals: [1,3] [5] [7,8] [10,11]. Input values may be repeated, with no
// consequence on the output.
//
// The output will satisfy the criteria of IntervalsAreSortedAndDisjoint().
std::vector<ClosedInterval> SortedDisjointIntervalsFromValues(
std::vector<int64> values);
// Returns true iff we have:
// - The intervals appear in increasing order.
// - for all i: intervals[i].start <= intervals[i].end
// - for all i but the last: intervals[i].end + 1 < intervals[i+1].start
//
// TODO(user): rename to IntervalsAreSortedAndNonAdjacent().
bool IntervalsAreSortedAndDisjoint(
const std::vector<ClosedInterval>& intervals);
// Returns true iff the given intervals contain the given value.
// We call "domain" any subset of Int64 = [kint64min, kint64max].
//
// TODO(user): This works in O(n), but could be made to work in O(log n) for
// long list of intervals.
bool SortedDisjointIntervalsContain(absl::Span<ClosedInterval> intervals,
int64 value);
// This class can be used to represent such set efficiently as a sorted and
// non-adjacent list of intervals. This is efficient as long as the size of such
// list stays reasonable.
//
// In the comments below, the domain of *this will always be written 'D'.
// Note that all the functions are safe with respect to integer overflow.
class Domain {
public:
// By default, Domain will be empty.
Domain() {}
// Returns the intersection of two lists of sorted disjoint intervals in a
// sorted disjoint interval form.
std::vector<ClosedInterval> IntersectionOfSortedDisjointIntervals(
const std::vector<ClosedInterval>& a, const std::vector<ClosedInterval>& b);
// Constructor for the common case of a singleton domain.
explicit Domain(int64 value);
// Returns the union of two lists of sorted disjoint intervals in a
// sorted disjoint interval form.
// Constructor for the common case of a single interval [left, right].
// If left > right, this will result in the empty domain.
Domain(int64 left, int64 right);
// Returns the full domain Int64.
static Domain AllValues();
// Creates a domain from the union of an unsorted list of integer values.
// Input values may be repeated, with no consequence on the output
static Domain FromValues(absl::Span<int64> values);
// Creates a domain from the union of an unsorted list of intervals.
static Domain FromIntervals(absl::Span<ClosedInterval> intervals);
// Returns true if this is the empty set.
bool IsEmpty() const;
// Returns the domain min/max value.
// This Checks that the domain is not empty.
int64 Min() const;
int64 Max() const;
// Returns true iff value is in Domain.
bool Contains(int64 value) const;
// Returns true iff D is included in the given domain.
bool IsIncludedIn(const Domain& domain) const;
// Returns the set Int64 D.
Domain Complement() const;
// Returns {x ∈ Int64, ∃ e ∈ D, x = -e}.
//
// Note in particular that if the negation of Int64 is not Int64 but
// Int64 \ {kint64min} !!
Domain Negation() const;
// Returns the set D ∩ domain.
Domain IntersectionWith(const Domain& domain) const;
// Returns the set D domain.
Domain UnionWith(const Domain& domain) const;
// Returns {x ∈ Int64, ∃ a ∈ D, ∃ b ∈ domain, x = a + b}.
Domain AdditionWith(const Domain& domain) const;
// Returns {x ∈ Int64, ∃ e ∈ D, x = e * coeff}.
//
// Note that because the resulting domain will only contains multiple of
// coeff, the size of intervals.size() can become really large. If it is
// larger than a fixed constant, success will be set to false and an empty
// Domain will be returned.
Domain MultiplicationBy(int64 coeff, bool* success) const;
// Returns a super-set of MultiplicationBy() to avoid the explosion in the
// representation size. This behaves as if we replace the set D of
// non-adjacent integer intervals by the set of floating-point element in the
// same intervals.
//
// For instance, [1, 100] * 2 will be transformed in [2, 200] and not in
// [2][4][6]...[200] like in MultiplicationBy(). Note that this would be
// similar to a InverseDivisionBy(), but not quite the same because if we
// look for {x ∈ Int64, ∃ e ∈ D, x / coeff = e}, then we will get [2, 201] in
// the case above.
Domain ContinuousMultiplicationBy(int64 coeff) const;
// Returns {x ∈ Int64, ∃ e ∈ D, x = e / coeff}.
//
// For instance Domain(1, 7).DivisionBy(2) == Domain(0, 3).
Domain DivisionBy(int64 coeff) const;
// Returns {x ∈ Int64, ∃ e ∈ D, x * coeff = e}.
//
// For instance Domain(1, 7).InverseMultiplicationBy(2) == Domain(1, 3).
Domain InverseMultiplicationBy(const int64 coeff) const;
// Returns the representation of D as the unique sorted list of non-adjacent
// intervals. This will always satisfy IntervalsAreSortedAndDisjoint().
std::vector<ClosedInterval> intervals() const { return intervals_; }
// Returns a compact std::string of a vector of intervals like
// "[1,4][6][10,20]"
std::string ToString() const;
// Lexicographic order on the intervals() representation.
bool operator<(const Domain& other) const;
bool operator==(const Domain& other) const {
return intervals_ == other.intervals_;
}
bool operator!=(const Domain& other) const {
return intervals_ != other.intervals_;
}
private:
// Invariant: will always satisfy IntervalsAreSortedAndDisjoint().
std::vector<ClosedInterval> intervals_;
};
std::ostream& operator<<(std::ostream& out, const Domain& domain);
// TODO(user): Remove. These are deprecated. Use the Domain class instead.
std::vector<ClosedInterval> UnionOfSortedDisjointIntervals(
const std::vector<ClosedInterval>& a, const std::vector<ClosedInterval>& b);
// Returns the domain of x + y given that the domain of x is a and the one of y
// is b.
std::vector<ClosedInterval> AdditionOfSortedDisjointIntervals(
const std::vector<ClosedInterval>& a, const std::vector<ClosedInterval>& b);
// Returns [kint64min, kint64max] minus the given intervals.
std::vector<ClosedInterval> ComplementOfSortedDisjointIntervals(
const std::vector<ClosedInterval>& intervals);
// For an x in the given intervals, this returns the domain of -x.
//
// Tricky: because the negation of kint64min doesn't fit, we always remove
// kint64min from the given intervals.
std::vector<ClosedInterval> NegationOfSortedDisjointIntervals(
std::vector<ClosedInterval> intervals);
// Returns the domain of x * coeff given the domain of x.
// To avoid an explosion in the size of the returned vector, the first function
// will actually return a super-set of the domain. For instance [1, 100] * 2
// will be transformed in [2, 200] not [2][4][6]...[200]. The second version
// will try to be exact as long as the result is not too large, and will set
// success to true when this is the case.
std::vector<ClosedInterval> MultiplicationOfSortedDisjointIntervals(
std::vector<ClosedInterval> intervals, int64 coeff);
std::vector<ClosedInterval> PreciseMultiplicationOfSortedDisjointIntervals(
std::vector<ClosedInterval> intervals, int64 coeff, bool* success);
// If x * coeff is in the given intervals, this returns the domain of x. Note
// that it is not the same as given the domains of x, return the domain of x /
// coeff because of how the integer division work.
std::vector<ClosedInterval> InverseMultiplicationOfSortedDisjointIntervals(
std::vector<ClosedInterval> intervals, int64 coeff);
// Given the domain of x, this returns the domain of x / coeff.
std::vector<ClosedInterval> DivisionOfSortedDisjointIntervals(
std::vector<ClosedInterval> intervals, int64 coeff);
// This class represents a sorted list of disjoint, closed intervals. When an
// interval is inserted, all intervals that overlap it or that are even adjacent
// to it are merged into one. I.e. [0,14] and [15,30] will be merged to [0,30].