diff --git a/ortools/flatzinc/cp_model_fz_solver.cc b/ortools/flatzinc/cp_model_fz_solver.cc index 628a59d0bd..d927d5943a 100644 --- a/ortools/flatzinc/cp_model_fz_solver.cc +++ b/ortools/flatzinc/cp_model_fz_solver.cc @@ -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{ + 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{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 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 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. diff --git a/ortools/glop/lp_solver.cc b/ortools/glop/lp_solver.cc index 1b3690ff81..b657b53009 100644 --- a/ortools/glop/lp_solver.cc +++ b/ortools/glop/lp_solver.cc @@ -174,8 +174,7 @@ ProblemStatus LPSolver::SolveWithTimeLimit(const LinearProgram& lp, current_linear_program_.PopulateFromLinearProgram(lp); // Preprocess. - MainLpPreprocessor preprocessor; - preprocessor.SetParameters(parameters_); + MainLpPreprocessor preprocessor(¶meters_); preprocessor.SetTimeLimit(time_limit); const bool postsolve_is_needed = preprocessor.Run(¤t_linear_program_); diff --git a/ortools/glop/parameters.proto b/ortools/glop/parameters.proto index 486a31ec74..44fe14a86b 100644 --- a/ortools/glop/parameters.proto +++ b/ortools/glop/parameters.proto @@ -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 diff --git a/ortools/glop/preprocessor.cc b/ortools/glop/preprocessor.cc index 32208e0a89..33e218952f 100644 --- a/ortools/glop/preprocessor.cc +++ b/ortools/glop/preprocessor.cc @@ -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(new name()), #name, \ - time_limit_, lp) +#define RUN_PREPROCESSOR(name) \ + RunAndPushIfRelevant(std::unique_ptr(new name(¶meters_)), \ + #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 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 diff --git a/ortools/glop/preprocessor.h b/ortools/glop/preprocessor.h index d1d95911fd..e1b8a07625 100644 --- a/ortools/glop/preprocessor.h +++ b/ortools/glop/preprocessor.h @@ -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 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 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 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 column_transform_; - - // Contains the initial problem bounds. - DenseRow variable_initial_lbs_; - DenseRow variable_initial_ubs_; - DISALLOW_COPY_AND_ASSIGN(SolowHalimPreprocessor); }; } // namespace glop diff --git a/ortools/linear_solver/linear_solver.cc b/ortools/linear_solver/linear_solver.cc index 3586080f4a..d9826c7bcb 100644 --- a/ortools/linear_solver/linear_solver.cc +++ b/ortools/linear_solver/linear_solver.cc @@ -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(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 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; diff --git a/ortools/linear_solver/linear_solver.h b/ortools/linear_solver/linear_solver.h index 2b4fa9cf0e..4ae1dbbb77 100644 --- a/ortools/linear_solver/linear_solver.h +++ b/ortools/linear_solver/linear_solver.h @@ -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 diff --git a/ortools/sat/cp_model_presolve.cc b/ortools/sat/cp_model_presolve.cc index 91811134a3..0cdb1939c0 100644 --- a/ortools/sat/cp_model_presolve.cc +++ b/ortools/sat/cp_model_presolve.cc @@ -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 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& 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 tmp_literals; - std::vector tmp_domain; - std::vector> tmp_term_domains; - std::vector> tmp_left_domains; + std::vector tmp_term_domains; + std::vector tmp_left_domains; // Each time a domain is modified this is set to true. SparseBitset modified_domains; private: // The current domain of each variables. - std::vector> domains; + std::vector domains; }; // ============================================================================= @@ -529,12 +528,10 @@ bool PresolveIntMax(ConstraintProto* ct, PresolveContext* context) { // Update the target domain. bool domain_reduced = false; if (!HasEnforcementLiteral(*ct)) { - std::vector 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 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& implied_rhs = left_domains[num_vars]; + const Domain& implied_rhs = left_domains[num_vars]; // Abort if intersection is empty. - const std::vector 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 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 new_domain; - std::vector 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 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 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 infered_domain; - const std::vector target_domain = - context->GetRefDomain(target_ref); - const std::vector 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( + PositiveRef(ref), Domain::FromValues(std::vector( 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 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 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 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. diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 4c6a827819..0259346363 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -304,9 +304,9 @@ std::vector ValuesFromProto(const Values& values) { } // Returns the size of the given domain capped to int64max. -int64 DomainSize(const std::vector& 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& domain) { // never have any trivial inequalities. void ModelWithMapping::ExtractEncoding(const CpModelProto& model_proto) { IntegerEncoder* encoder = GetOrCreate(); + IntegerTrail* integer_trail = GetOrCreate(); // 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 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 union_of_non_constant_domains; + Domain union_of_non_constant_domains; std::map 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 constant_values; - std::map, int, ExactVectorOfDomainComparator> - num_vars_per_domains; + std::map 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(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(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()->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? diff --git a/ortools/sat/cp_model_utils.h b/ortools/sat/cp_model_utils.h index 9b4cf49574..f57e40194f 100644 --- a/ortools/sat/cp_model_utils.h +++ b/ortools/sat/cp_model_utils.h @@ -92,6 +92,10 @@ void FillDomain(const std::vector& domain, proto->add_domain(interval.end); } } +template +void FillDomainInProto(const Domain& domain, ProtoWithDomain* proto) { + FillDomain(domain.intervals(), proto); +} // Extract a sorted interval list from the domain field of a proto. template @@ -103,6 +107,10 @@ std::vector ReadDomain(const ProtoWithDomain& proto) { CHECK(IntervalsAreSortedAndDisjoint(result)); return result; } +template +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. diff --git a/ortools/sat/disjunctive.cc b/ortools/sat/disjunctive.cc index 7dd1fac72b..a895d43fd1 100644 --- a/ortools/sat/disjunctive.cc +++ b/ortools/sat/disjunctive.cc @@ -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); diff --git a/ortools/sat/integer.cc b/ortools/sat/integer.cc index 7f38e2252d..b81f442c8f 100644 --- a/ortools/sat/integer.cc +++ b/ortools/sat/integer.cc @@ -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& 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 IntegerTrail::InitialVariableDomain( - IntegerVariable var) const { - std::vector 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 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 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 temp = - NegationOfSortedDisjointIntervals(domain); + const std::vector temp = domain.intervals(); + (*domains_)[var].assign(temp.begin(), temp.end()); + } + { + const std::vector 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 coeffs, + std::vector* reason) const { + CHECK_GE(slack, 0); + if (slack == 0) return; + const int size = reason->size(); + std::vector 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_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); diff --git a/ortools/sat/integer.h b/ortools/sat/integer.h index a80f1d5deb..693ac12752 100644 --- a/ortools/sat/integer.h +++ b/ortools/sat/integer.h @@ -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& 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 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 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 coeffs, + std::vector* 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 NewIntegerVariable(int64 lb, } inline std::function NewIntegerVariable( - const std::vector& domain) { + const Domain& domain) { return [=](Model* model) { return model->GetOrCreate()->AddIntegerVariable(domain); }; diff --git a/ortools/sat/integer_expr.cc b/ortools/sat/integer_expr.cc index d86ae462fd..07969e27da 100644 --- a/ortools/sat/integer_expr.cc +++ b/ortools/sat/integer_expr.cc @@ -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 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 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; diff --git a/ortools/sat/intervals.cc b/ortools/sat/intervals.cc index d5595853ec..1844f05340 100644 --- a/ortools/sat/intervals.cc +++ b/ortools/sat/intervals.cc @@ -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 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_); } diff --git a/ortools/sat/intervals.h b/ortools/sat/intervals.h index 20f60966c9..69aa8dec25 100644 --- a/ortools/sat/intervals.h +++ b/ortools/sat/intervals.h @@ -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* MutableLiteralReason() { return &literal_reason_; } diff --git a/ortools/sat/precedences.cc b/ortools/sat/precedences.cc index 34346358db..0636d550e3 100644 --- a/ortools/sat/precedences.cc +++ b/ortools/sat/precedences.cc @@ -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 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_); diff --git a/ortools/sat/python/cp_model.py b/ortools/sat/python/cp_model.py index 793d8eedf1..cf9eeb5cd8 100644 --- a/ortools/sat/python/cp_model.py +++ b/ortools/sat/python/cp_model.py @@ -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. diff --git a/ortools/sat/table.cc b/ortools/sat/table.cc index c348f2e4df..cc047bce8f 100644 --- a/ortools/sat/table.cc +++ b/ortools/sat/table.cc @@ -61,13 +61,12 @@ std::unordered_map GetEncoding(IntegerVariable var, void FilterValues(IntegerVariable var, Model* model, std::unordered_set* values) { - std::vector domain = - model->Get()->InitialVariableDomain(var); + const Domain domain = model->Get()->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 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 TransitionConstraint( const auto domain = integer_trail->InitialVariableDomain(vars[time]); for (const std::vector& 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 TransitionConstraint( std::vector 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 { diff --git a/ortools/util/sorted_interval_list.cc b/ortools/util/sorted_interval_list.cc index a674d5cf67..d8cba450d1 100644 --- a/ortools/util/sorted_interval_list.cc +++ b/ortools/util/sorted_interval_list.cc @@ -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& d1, - const std::vector& 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& intervals) { std::string result; for (ClosedInterval interval : intervals) { @@ -54,13 +34,13 @@ std::string IntervalsAsString(const std::vector& 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& intervals) { - return out << IntervalsAsString(intervals); +// TODO(user): binary search if size is large? +bool SortedDisjointIntervalsContain(absl::Span intervals, + int64 value) { + for (const ClosedInterval& interval : intervals) { + if (interval.start <= value && interval.end >= value) return true; + } + return false; } std::vector SortedDisjointIntervalsFromValues( @@ -95,14 +75,27 @@ bool IntervalsAreSortedAndDisjoint( return true; } -bool SortedDisjointIntervalsContain(absl::Span 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* 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 IntersectionOfSortedDisjointIntervals( const std::vector& a, const std::vector& b) { @@ -163,27 +156,6 @@ std::vector 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* 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 UnionOfSortedDisjointIntervals( const std::vector& a, const std::vector& b) { @@ -297,6 +269,151 @@ std::vector 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& 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 values) { + Domain result; + result.intervals_ = SortedDisjointIntervalsFromValues( + std::vector(values.begin(), values.end())); + return result; +} + +Domain Domain::FromIntervals(absl::Span 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& 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& d1 = intervals_; + const std::vector& 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( diff --git a/ortools/util/sorted_interval_list.h b/ortools/util/sorted_interval_list.h index 9a20c2a9d5..53eb5464ab 100644 --- a/ortools/util/sorted_interval_list.h +++ b/ortools/util/sorted_interval_list.h @@ -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& d1, - const std::vector& d2) const; -}; - -// Returns a compact std::string of a vector of intervals like -// "[1,4][6][10,20]". -std::string IntervalsAsString(const std::vector& intervals); - std::ostream& operator<<(std::ostream& out, const ClosedInterval& interval); std::ostream& operator<<(std::ostream& out, const std::vector& 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 SortedDisjointIntervalsFromValues( - std::vector 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& 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 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 IntersectionOfSortedDisjointIntervals( - const std::vector& a, const std::vector& 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 values); + + // Creates a domain from the union of an unsorted list of intervals. + static Domain FromIntervals(absl::Span 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 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 intervals_; +}; + +std::ostream& operator<<(std::ostream& out, const Domain& domain); + +// TODO(user): Remove. These are deprecated. Use the Domain class instead. std::vector UnionOfSortedDisjointIntervals( const std::vector& a, const std::vector& b); - -// Returns the domain of x + y given that the domain of x is a and the one of y -// is b. -std::vector AdditionOfSortedDisjointIntervals( - const std::vector& a, const std::vector& b); - -// Returns [kint64min, kint64max] minus the given intervals. -std::vector ComplementOfSortedDisjointIntervals( - const std::vector& 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 NegationOfSortedDisjointIntervals( std::vector 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 MultiplicationOfSortedDisjointIntervals( - std::vector intervals, int64 coeff); -std::vector PreciseMultiplicationOfSortedDisjointIntervals( - std::vector 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 InverseMultiplicationOfSortedDisjointIntervals( - std::vector intervals, int64 coeff); - -// Given the domain of x, this returns the domain of x / coeff. -std::vector DivisionOfSortedDisjointIntervals( - std::vector 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].