28 using ::operations_research::glop::ColIndex;
32 using ::operations_research::glop::LinearProgram;
33 using ::operations_research::glop::LPDecomposer;
34 using ::operations_research::glop::RowIndex;
35 using ::operations_research::glop::SparseColumn;
36 using ::operations_research::glop::SparseMatrix;
37 using ::operations_research::sat::LinearBooleanConstraint;
38 using ::operations_research::sat::LinearBooleanProblem;
39 using ::operations_research::sat::LinearObjective;
44 const double kTolerance = 1e-10;
45 return std::abs(x - round(x)) <= kTolerance;
51 bool ProblemIsBooleanAndHasOnlyIntegralConstraints(
52 const LinearProgram& linear_problem) {
53 const glop::SparseMatrix& matrix = linear_problem.GetSparseMatrix();
55 for (ColIndex
col(0);
col < linear_problem.num_variables(); ++
col) {
56 const Fractional lower_bound = linear_problem.variable_lower_bounds()[
col];
57 const Fractional upper_bound = linear_problem.variable_upper_bounds()[
col];
59 if (lower_bound <= -1.0 || upper_bound >= 2.0) {
64 for (
const SparseColumn::Entry e : matrix.column(
col)) {
77 void BuildBooleanProblemWithIntegralConstraints(
78 const LinearProgram& linear_problem,
const DenseRow& initial_solution,
79 LinearBooleanProblem* boolean_problem,
80 std::vector<bool>* boolean_initial_solution) {
81 CHECK(boolean_problem !=
nullptr);
82 boolean_problem->Clear();
84 const glop::SparseMatrix& matrix = linear_problem.GetSparseMatrix();
86 for (ColIndex
col(0);
col < matrix.num_cols(); ++
col) {
87 boolean_problem->add_var_names(linear_problem.GetVariableName(
col));
89 boolean_problem->set_num_variables(matrix.num_cols().value());
90 boolean_problem->set_name(linear_problem.name());
93 for (RowIndex
row(0);
row < matrix.num_rows(); ++
row) {
94 LinearBooleanConstraint*
const constraint =
95 boolean_problem->add_constraints();
96 constraint->set_name(linear_problem.GetConstraintName(
row));
97 if (linear_problem.constraint_lower_bounds()[
row] != -
kInfinity) {
98 constraint->set_lower_bound(
99 linear_problem.constraint_lower_bounds()[
row]);
101 if (linear_problem.constraint_upper_bounds()[
row] !=
kInfinity) {
102 constraint->set_upper_bound(
103 linear_problem.constraint_upper_bounds()[
row]);
108 for (ColIndex
col(0);
col < matrix.num_cols(); ++
col) {
109 for (
const SparseColumn::Entry e : matrix.column(
col)) {
110 LinearBooleanConstraint*
const constraint =
111 boolean_problem->mutable_constraints(e.row().value());
112 constraint->add_literals(
col.value() + 1);
113 constraint->add_coefficients(e.coefficient());
119 for (ColIndex
col(0);
col < matrix.num_cols(); ++
col) {
121 const int lb = std::round(linear_problem.variable_lower_bounds()[
col]);
122 const int ub = std::round(linear_problem.variable_upper_bounds()[
col]);
124 LinearBooleanConstraint*
ct = boolean_problem->add_constraints();
125 ct->set_lower_bound(ub);
126 ct->set_upper_bound(ub);
127 ct->add_literals(
col.value() + 1);
128 ct->add_coefficients(1.0);
134 for (ColIndex
col(0);
col < linear_problem.num_variables(); ++
col) {
135 const Fractional coeff = linear_problem.objective_coefficients()[
col];
138 double scaling_factor = 0.0;
139 double relative_error = 0.0;
142 &scaling_factor, &relative_error);
144 LinearObjective*
const objective = boolean_problem->mutable_objective();
145 objective->set_offset(linear_problem.objective_offset() * scaling_factor /
150 objective->set_scaling_factor(1.0 / scaling_factor * gcd);
151 for (ColIndex
col(0);
col < linear_problem.num_variables(); ++
col) {
152 const Fractional coeff = linear_problem.objective_coefficients()[
col];
153 const int64_t
value =
154 static_cast<int64_t
>(round(coeff * scaling_factor)) / gcd;
156 objective->add_literals(
col.value() + 1);
157 objective->add_coefficients(
value);
162 if (linear_problem.IsMaximizationProblem()) {
167 if (!initial_solution.empty()) {
168 CHECK(boolean_initial_solution !=
nullptr);
169 CHECK_EQ(boolean_problem->num_variables(), initial_solution.size());
170 boolean_initial_solution->assign(boolean_problem->num_variables(),
false);
171 for (
int i = 0; i < initial_solution.size(); ++i) {
172 (*boolean_initial_solution)[i] = (initial_solution[ColIndex(i)] != 0);
185 class IntegralVariable {
194 void BuildFromRange(
int start_var_index,
Fractional lower_bound,
198 void set_offset(int64_t offset) {
offset_ = offset; }
199 void set_weight(VariableIndex
var, int64_t
weight);
201 int GetNumberOfBooleanVariables()
const {
return bits_.size(); }
203 const std::vector<VariableIndex>& bits()
const {
return bits_; }
204 const std::vector<int64_t>& weights()
const {
return weights_; }
205 int64_t offset()
const {
return offset_; }
209 int64_t GetSolutionValue(
const BopSolution& solution)
const;
215 std::vector<bool> GetBooleanSolutionValues(int64_t integral_value)
const;
217 std::string DebugString()
const;
223 std::vector<VariableIndex> bits_;
224 std::vector<int64_t> weights_;
230 bool can_be_reversed_;
233 IntegralVariable::IntegralVariable()
234 : bits_(), weights_(),
offset_(0), can_be_reversed_(true) {}
236 void IntegralVariable::BuildFromRange(
int start_var_index,
246 const int64_t integral_lower_bound =
static_cast<int64_t
>(ceil(lower_bound));
247 const int64_t integral_upper_bound =
static_cast<int64_t
>(floor(upper_bound));
248 offset_ = integral_lower_bound;
249 const int64_t
delta = integral_upper_bound - integral_lower_bound;
251 for (
int i = 0; i < num_used_bits; ++i) {
252 bits_.push_back(VariableIndex(start_var_index + i));
253 weights_.push_back(1ULL << i);
257 void IntegralVariable::Clear() {
261 can_be_reversed_ =
true;
264 void IntegralVariable::set_weight(VariableIndex
var, int64_t
weight) {
265 bits_.push_back(
var);
266 weights_.push_back(
weight);
267 can_be_reversed_ =
false;
270 int64_t IntegralVariable::GetSolutionValue(
const BopSolution& solution)
const {
272 for (
int i = 0; i < bits_.size(); ++i) {
273 value += weights_[i] * solution.Value(bits_[i]);
278 std::vector<bool> IntegralVariable::GetBooleanSolutionValues(
279 int64_t integral_value)
const {
280 if (can_be_reversed_) {
281 DCHECK(std::is_sorted(weights_.begin(), weights_.end()));
282 std::vector<bool> boolean_values(weights_.size(),
false);
283 int64_t remaining_value = integral_value -
offset_;
284 for (
int i = weights_.size() - 1; i >= 0; --i) {
285 if (remaining_value >= weights_[i]) {
286 boolean_values[i] =
true;
287 remaining_value -= weights_[i];
291 <<
"Couldn't map integral value to boolean variables.";
292 return boolean_values;
294 return std::vector<bool>();
297 std::string IntegralVariable::DebugString()
const {
299 CHECK_EQ(bits_.size(), weights_.size());
300 for (
int i = 0; i < bits_.size(); ++i) {
301 str += absl::StrFormat(
"%d [%d] ", weights_[i], bits_[i].
value());
303 str += absl::StrFormat(
" Offset: %d",
offset_);
324 class IntegralProblemConverter {
326 IntegralProblemConverter();
332 bool ConvertToBooleanProblem(
const LinearProgram& linear_problem,
334 LinearBooleanProblem* boolean_problem,
335 std::vector<bool>* boolean_initial_solution);
339 int64_t GetSolutionValue(ColIndex global_col,
340 const BopSolution& solution)
const;
346 bool CheckProblem(
const LinearProgram& linear_problem)
const;
349 void InitVariableTypes(
const LinearProgram& linear_problem,
350 LinearBooleanProblem* boolean_problem);
353 void ConvertAllVariables(
const LinearProgram& linear_problem,
354 LinearBooleanProblem* boolean_problem);
357 void AddVariableConstraints(
const LinearProgram& linear_problem,
358 LinearBooleanProblem* boolean_problem);
361 void ConvertAllConstraints(
const LinearProgram& linear_problem,
362 LinearBooleanProblem* boolean_problem);
365 void ConvertObjective(
const LinearProgram& linear_problem,
366 LinearBooleanProblem* boolean_problem);
372 bool ConvertUsingExistingBooleans(
const LinearProgram& linear_problem,
374 IntegralVariable* integral_var);
385 bool CreateVariableUsingConstraint(
const LinearProgram& linear_problem,
387 IntegralVariable* integral_var);
401 double ScaleAndSparsifyWeights(
402 double scaling_factor, int64_t gcd,
406 bool HasNonZeroWeights(
409 bool problem_is_boolean_and_has_only_integral_constraints_;
416 std::vector<IntegralVariable> integral_variables_;
417 std::vector<ColIndex> integral_indices_;
418 int num_boolean_variables_;
420 enum VariableType { BOOLEAN, INTEGRAL, INTEGRAL_EXPRESSED_AS_BOOLEAN };
424 IntegralProblemConverter::IntegralProblemConverter()
425 : global_to_boolean_(),
426 integral_variables_(),
428 num_boolean_variables_(0),
431 bool IntegralProblemConverter::ConvertToBooleanProblem(
432 const LinearProgram& linear_problem,
const DenseRow& initial_solution,
433 LinearBooleanProblem* boolean_problem,
434 std::vector<bool>* boolean_initial_solution) {
435 bool use_initial_solution = (initial_solution.size() > 0);
436 if (use_initial_solution) {
437 CHECK_EQ(initial_solution.size(), linear_problem.num_variables())
438 <<
"The initial solution should have the same number of variables as "
439 "the LinearProgram.";
440 CHECK(boolean_initial_solution !=
nullptr);
442 if (!CheckProblem(linear_problem)) {
446 problem_is_boolean_and_has_only_integral_constraints_ =
447 ProblemIsBooleanAndHasOnlyIntegralConstraints(linear_problem);
448 if (problem_is_boolean_and_has_only_integral_constraints_) {
449 BuildBooleanProblemWithIntegralConstraints(linear_problem, initial_solution,
451 boolean_initial_solution);
455 InitVariableTypes(linear_problem, boolean_problem);
456 ConvertAllVariables(linear_problem, boolean_problem);
457 boolean_problem->set_num_variables(num_boolean_variables_);
458 boolean_problem->set_name(linear_problem.name());
460 AddVariableConstraints(linear_problem, boolean_problem);
461 ConvertAllConstraints(linear_problem, boolean_problem);
462 ConvertObjective(linear_problem, boolean_problem);
465 if (linear_problem.IsMaximizationProblem()) {
469 if (use_initial_solution) {
470 boolean_initial_solution->assign(boolean_problem->num_variables(),
false);
471 for (ColIndex global_col(0); global_col < global_to_boolean_.size();
473 const int col = global_to_boolean_[global_col];
475 (*boolean_initial_solution)[
col] = (initial_solution[global_col] != 0);
477 const IntegralVariable& integral_variable =
478 integral_variables_[-
col - 1];
479 const std::vector<VariableIndex>& boolean_cols =
480 integral_variable.bits();
481 const std::vector<bool>& boolean_values =
482 integral_variable.GetBooleanSolutionValues(
483 round(initial_solution[global_col]));
484 if (!boolean_values.empty()) {
485 CHECK_EQ(boolean_cols.size(), boolean_values.size());
486 for (
int i = 0; i < boolean_values.size(); ++i) {
487 const int boolean_col = boolean_cols[i].value();
488 (*boolean_initial_solution)[boolean_col] = boolean_values[i];
498 int64_t IntegralProblemConverter::GetSolutionValue(
499 ColIndex global_col,
const BopSolution& solution)
const {
500 if (problem_is_boolean_and_has_only_integral_constraints_) {
501 return solution.Value(VariableIndex(global_col.value()));
504 const int pos = global_to_boolean_[global_col];
505 return pos >= 0 ? solution.Value(VariableIndex(pos))
506 : integral_variables_[-pos - 1].GetSolutionValue(solution);
509 bool IntegralProblemConverter::CheckProblem(
510 const LinearProgram& linear_problem)
const {
511 for (ColIndex
col(0);
col < linear_problem.num_variables(); ++
col) {
512 if (!linear_problem.IsVariableInteger(
col)) {
513 LOG(
ERROR) <<
"Variable " << linear_problem.GetVariableName(
col)
514 <<
" is continuous. This is not supported by BOP.";
517 if (linear_problem.variable_lower_bounds()[
col] == -
kInfinity) {
518 LOG(
ERROR) <<
"Variable " << linear_problem.GetVariableName(
col)
519 <<
" has no lower bound. This is not supported by BOP.";
522 if (linear_problem.variable_upper_bounds()[
col] ==
kInfinity) {
523 LOG(
ERROR) <<
"Variable " << linear_problem.GetVariableName(
col)
524 <<
" has no upper bound. This is not supported by BOP.";
531 void IntegralProblemConverter::InitVariableTypes(
532 const LinearProgram& linear_problem,
533 LinearBooleanProblem* boolean_problem) {
534 global_to_boolean_.assign(linear_problem.num_variables().value(), 0);
535 variable_types_.assign(linear_problem.num_variables().value(), INTEGRAL);
536 for (ColIndex
col(0);
col < linear_problem.num_variables(); ++
col) {
537 const Fractional lower_bound = linear_problem.variable_lower_bounds()[
col];
538 const Fractional upper_bound = linear_problem.variable_upper_bounds()[
col];
540 if (lower_bound > -1.0 && upper_bound < 2.0) {
542 variable_types_[
col] = BOOLEAN;
543 global_to_boolean_[
col] = num_boolean_variables_;
544 ++num_boolean_variables_;
545 boolean_problem->add_var_names(linear_problem.GetVariableName(
col));
548 variable_types_[
col] = INTEGRAL;
549 integral_indices_.push_back(
col);
554 void IntegralProblemConverter::ConvertAllVariables(
555 const LinearProgram& linear_problem,
556 LinearBooleanProblem* boolean_problem) {
557 for (
const ColIndex
col : integral_indices_) {
559 IntegralVariable integral_var;
560 if (!ConvertUsingExistingBooleans(linear_problem,
col, &integral_var)) {
562 linear_problem.variable_lower_bounds()[
col];
564 linear_problem.variable_upper_bounds()[
col];
565 integral_var.BuildFromRange(num_boolean_variables_, lower_bound,
567 num_boolean_variables_ += integral_var.GetNumberOfBooleanVariables();
568 const std::string var_name = linear_problem.GetVariableName(
col);
569 for (
int i = 0; i < integral_var.bits().size(); ++i) {
570 boolean_problem->add_var_names(var_name + absl::StrFormat(
"_%d", i));
573 integral_variables_.push_back(integral_var);
574 global_to_boolean_[
col] = -integral_variables_.size();
575 variable_types_[
col] = INTEGRAL_EXPRESSED_AS_BOOLEAN;
579 void IntegralProblemConverter::ConvertAllConstraints(
580 const LinearProgram& linear_problem,
581 LinearBooleanProblem* boolean_problem) {
584 glop::SparseMatrix transpose;
585 transpose.PopulateFromTranspose(linear_problem.GetSparseMatrix());
587 double max_relative_error = 0.0;
588 double max_bound_error = 0.0;
590 double relative_error = 0.0;
591 double scaling_factor = 0.0;
593 for (RowIndex
row(0);
row < linear_problem.num_constraints(); ++
row) {
596 num_boolean_variables_, 0.0);
599 offset += AddWeightedIntegralVariable(
RowToColIndex(e.row()),
600 e.coefficient(), &dense_weights);
602 if (!HasNonZeroWeights(dense_weights)) {
608 for (VariableIndex
var(0);
var < num_boolean_variables_; ++
var) {
609 if (dense_weights[
var] != 0.0) {
615 &scaling_factor, &relative_error);
618 max_relative_error =
std::max(relative_error, max_relative_error);
621 LinearBooleanConstraint* constraint = boolean_problem->add_constraints();
622 constraint->set_name(linear_problem.GetConstraintName(
row));
623 const double bound_error =
624 ScaleAndSparsifyWeights(scaling_factor, gcd, dense_weights, constraint);
625 max_bound_error =
std::max(max_bound_error, bound_error);
628 linear_problem.constraint_lower_bounds()[
row];
630 const Fractional offset_lower_bound = lower_bound - offset;
631 const double offset_scaled_lower_bound =
632 round(offset_lower_bound * scaling_factor - bound_error);
633 if (offset_scaled_lower_bound >=
635 LOG(
WARNING) <<
"A constraint is trivially unsatisfiable.";
638 if (offset_scaled_lower_bound >
641 constraint->set_lower_bound(
642 static_cast<int64_t
>(offset_scaled_lower_bound) / gcd);
646 linear_problem.constraint_upper_bounds()[
row];
648 const Fractional offset_upper_bound = upper_bound - offset;
649 const double offset_scaled_upper_bound =
650 round(offset_upper_bound * scaling_factor + bound_error);
651 if (offset_scaled_upper_bound <=
653 LOG(
WARNING) <<
"A constraint is trivially unsatisfiable.";
656 if (offset_scaled_upper_bound <
659 constraint->set_upper_bound(
660 static_cast<int64_t
>(offset_scaled_upper_bound) / gcd);
666 void IntegralProblemConverter::ConvertObjective(
667 const LinearProgram& linear_problem,
668 LinearBooleanProblem* boolean_problem) {
669 LinearObjective* objective = boolean_problem->mutable_objective();
672 num_boolean_variables_, 0.0);
674 for (ColIndex
col(0);
col < linear_problem.num_variables(); ++
col) {
675 offset += AddWeightedIntegralVariable(
676 col, linear_problem.objective_coefficients()[
col], &dense_weights);
681 for (VariableIndex
var(0);
var < num_boolean_variables_; ++
var) {
682 if (dense_weights[
var] != 0.0) {
686 double scaling_factor = 0.0;
687 double max_relative_error = 0.0;
688 double relative_error = 0.0;
691 &scaling_factor, &relative_error);
693 max_relative_error =
std::max(relative_error, max_relative_error);
694 VLOG(1) <<
"objective relative error: " << relative_error;
695 VLOG(1) <<
"objective scaling factor: " << scaling_factor / gcd;
697 ScaleAndSparsifyWeights(scaling_factor, gcd, dense_weights, objective);
701 objective->set_scaling_factor(1.0 / scaling_factor * gcd);
702 objective->set_offset((linear_problem.objective_offset() + offset) *
703 scaling_factor / gcd);
706 void IntegralProblemConverter::AddVariableConstraints(
707 const LinearProgram& linear_problem,
708 LinearBooleanProblem* boolean_problem) {
709 for (ColIndex
col(0);
col < linear_problem.num_variables(); ++
col) {
710 const Fractional lower_bound = linear_problem.variable_lower_bounds()[
col];
711 const Fractional upper_bound = linear_problem.variable_upper_bounds()[
col];
712 const int pos = global_to_boolean_[
col];
716 const bool is_fixed = (lower_bound > -1.0 && upper_bound < 1.0) ||
717 (lower_bound > 0.0 && upper_bound < 2.0);
720 const int fixed_value = lower_bound > -1.0 && upper_bound < 1.0 ? 0 : 1;
721 LinearBooleanConstraint* constraint =
722 boolean_problem->add_constraints();
723 constraint->set_lower_bound(fixed_value);
724 constraint->set_upper_bound(fixed_value);
725 constraint->add_literals(pos + 1);
726 constraint->add_coefficients(1);
729 CHECK_EQ(INTEGRAL_EXPRESSED_AS_BOOLEAN, variable_types_[
col]);
732 const IntegralVariable& integral_var = integral_variables_[-pos - 1];
733 LinearBooleanConstraint* constraint =
734 boolean_problem->add_constraints();
735 for (
int i = 0; i < integral_var.bits().size(); ++i) {
736 constraint->add_literals(integral_var.bits()[i].value() + 1);
737 constraint->add_coefficients(integral_var.weights()[i]);
740 constraint->set_lower_bound(
static_cast<int64_t
>(ceil(lower_bound)) -
741 integral_var.offset());
744 constraint->set_upper_bound(
static_cast<int64_t
>(floor(upper_bound)) -
745 integral_var.offset());
752 bool IntegralProblemConverter::ConvertUsingExistingBooleans(
753 const LinearProgram& linear_problem, ColIndex
col,
754 IntegralVariable* integral_var) {
755 CHECK(
nullptr != integral_var);
758 const SparseMatrix& matrix = linear_problem.GetSparseMatrix();
759 const SparseMatrix& transpose = linear_problem.GetTransposeSparseMatrix();
760 for (
const SparseColumn::Entry var_entry : matrix.column(
col)) {
761 const RowIndex constraint = var_entry.row();
762 const Fractional lb = linear_problem.constraint_lower_bounds()[constraint];
763 const Fractional ub = linear_problem.constraint_upper_bounds()[constraint];
779 bool only_one_integral_variable =
true;
780 for (
const SparseColumn::Entry constraint_entry :
782 const ColIndex var_index =
RowToColIndex(constraint_entry.row());
783 if (var_index !=
col && variable_types_[var_index] == INTEGRAL) {
784 only_one_integral_variable =
false;
788 if (only_one_integral_variable &&
789 CreateVariableUsingConstraint(linear_problem, constraint,
795 integral_var->Clear();
799 bool IntegralProblemConverter::CreateVariableUsingConstraint(
800 const LinearProgram& linear_problem, RowIndex constraint,
801 IntegralVariable* integral_var) {
802 CHECK(
nullptr != integral_var);
803 integral_var->Clear();
805 const SparseMatrix& transpose = linear_problem.GetTransposeSparseMatrix();
807 num_boolean_variables_, 0.0);
809 int64_t variable_offset = 0;
810 for (
const SparseColumn::Entry constraint_entry :
813 if (variable_types_[
col] == INTEGRAL) {
814 scale = constraint_entry.coefficient();
815 }
else if (variable_types_[
col] == BOOLEAN) {
816 const int pos = global_to_boolean_[
col];
818 dense_weights[VariableIndex(pos)] -= constraint_entry.coefficient();
820 CHECK_EQ(INTEGRAL_EXPRESSED_AS_BOOLEAN, variable_types_[
col]);
821 const int pos = global_to_boolean_[
col];
823 const IntegralVariable& local_integral_var =
824 integral_variables_[-pos - 1];
826 constraint_entry.coefficient() * local_integral_var.offset();
827 for (
int i = 0; i < local_integral_var.bits().size(); ++i) {
828 dense_weights[local_integral_var.bits()[i]] -=
829 constraint_entry.coefficient() * local_integral_var.weights()[i];
835 const Fractional lb = linear_problem.constraint_lower_bounds()[constraint];
836 const Fractional offset = (lb + variable_offset) / scale;
840 integral_var->set_offset(
static_cast<int64_t
>(offset));
842 for (VariableIndex
var(0);
var < dense_weights.size(); ++
var) {
843 if (dense_weights[
var] != 0.0) {
848 integral_var->set_weight(
var,
static_cast<int64_t
>(
weight));
855 Fractional IntegralProblemConverter::AddWeightedIntegralVariable(
858 CHECK(
nullptr != dense_weights);
865 const int pos = global_to_boolean_[
col];
868 (*dense_weights)[VariableIndex(pos)] +=
weight;
871 const IntegralVariable& integral_var = integral_variables_[-pos - 1];
872 for (
int i = 0; i < integral_var.bits().size(); ++i) {
873 (*dense_weights)[integral_var.bits()[i]] +=
874 integral_var.weights()[i] *
weight;
876 offset +=
weight * integral_var.offset();
882 double IntegralProblemConverter::ScaleAndSparsifyWeights(
883 double scaling_factor, int64_t gcd,
887 double bound_error = 0.0;
888 for (VariableIndex
var(0);
var < dense_weights.
size(); ++
var) {
889 if (dense_weights[
var] != 0.0) {
890 const double scaled_weight = dense_weights[
var] * scaling_factor;
891 bound_error += fabs(round(scaled_weight) - scaled_weight);
892 t->add_literals(
var.value() + 1);
893 t->add_coefficients(
static_cast<int64_t
>(round(scaled_weight)) / gcd);
899 bool IntegralProblemConverter::HasNonZeroWeights(
913 const SparseMatrix& matrix = linear_problem.GetSparseMatrix();
914 for (ColIndex
col(0);
col < linear_problem.num_variables(); ++
col) {
915 const Fractional lower_bound = linear_problem.variable_lower_bounds()[
col];
916 const Fractional upper_bound = linear_problem.variable_upper_bounds()[
col];
918 if (lower_bound >
value || upper_bound <
value) {
920 <<
" should be in " << lower_bound <<
" .. " << upper_bound;
924 for (
const SparseColumn::Entry entry : matrix.column(
col)) {
925 constraint_values[entry.row()] += entry.coefficient() *
value;
929 for (RowIndex
row(0);
row < linear_problem.num_constraints(); ++
row) {
930 const Fractional lb = linear_problem.constraint_lower_bounds()[
row];
931 const Fractional ub = linear_problem.constraint_upper_bounds()[
row];
935 <<
" should be in " << lb <<
" .. " << ub;
950 CHECK(variable_values !=
nullptr);
951 CHECK(objective_value !=
nullptr);
952 CHECK(best_bound !=
nullptr);
953 const bool use_initial_solution = (initial_solution.size() > 0);
954 if (use_initial_solution) {
955 CHECK_EQ(initial_solution.size(), linear_problem.num_variables());
961 variable_values->resize(linear_problem.num_variables(), 0);
963 LinearBooleanProblem boolean_problem;
964 std::vector<bool> boolean_initial_solution;
965 IntegralProblemConverter converter;
966 if (!converter.ConvertToBooleanProblem(linear_problem, initial_solution,
968 &boolean_initial_solution)) {
969 return BopSolveStatus::INVALID_PROBLEM;
972 BopSolver bop_solver(boolean_problem);
975 if (use_initial_solution) {
976 BopSolution bop_solution(boolean_problem,
"InitialSolution");
977 CHECK_EQ(boolean_initial_solution.size(), boolean_problem.num_variables());
978 for (
int i = 0; i < boolean_initial_solution.size(); ++i) {
979 bop_solution.SetValue(VariableIndex(i), boolean_initial_solution[i]);
981 status = bop_solver.SolveWithTimeLimit(bop_solution,
time_limit);
983 status = bop_solver.SolveWithTimeLimit(
time_limit);
985 if (status == BopSolveStatus::OPTIMAL_SOLUTION_FOUND ||
986 status == BopSolveStatus::FEASIBLE_SOLUTION_FOUND) {
988 const BopSolution& solution = bop_solver.best_solution();
989 CHECK(solution.IsFeasible());
991 *objective_value = linear_problem.objective_offset();
992 for (ColIndex
col(0);
col < linear_problem.num_variables(); ++
col) {
993 const int64_t
value = converter.GetSolutionValue(
col, solution);
995 *objective_value +=
value * linear_problem.objective_coefficients()[
col];
1002 *best_bound = status == BopSolveStatus::OPTIMAL_SOLUTION_FOUND
1004 : bop_solver.GetScaledBestBound();
1009 void RunOneBop(
const BopParameters&
parameters,
int problem_index,
1011 LPDecomposer* decomposer,
DenseRow* variable_values,
1014 CHECK(decomposer !=
nullptr);
1015 CHECK(variable_values !=
nullptr);
1016 CHECK(objective_value !=
nullptr);
1017 CHECK(best_bound !=
nullptr);
1018 CHECK(status !=
nullptr);
1020 LinearProgram problem;
1021 decomposer->ExtractLocalProblem(problem_index, &problem);
1023 if (initial_solution.size() > 0) {
1024 local_initial_solution =
1025 decomposer->ExtractLocalAssignment(problem_index, initial_solution);
1029 const double total_num_variables =
std::max(
1030 1.0,
static_cast<double>(
1031 decomposer->original_problem().num_variables().value()));
1032 const double time_per_variable =
1033 parameters.max_time_in_seconds() / total_num_variables;
1034 const double deterministic_time_per_variable =
1035 parameters.max_deterministic_time() / total_num_variables;
1036 const int local_num_variables =
std::max(1, problem.num_variables().value());
1038 NestedTimeLimit subproblem_time_limit(
1040 std::max(time_per_variable * local_num_variables,
1041 parameters.decomposed_problem_min_time_in_seconds()),
1042 deterministic_time_per_variable * local_num_variables);
1044 *status = InternalSolve(problem,
parameters, local_initial_solution,
1045 subproblem_time_limit.GetTimeLimit(), variable_values,
1046 objective_value, best_bound);
1050 IntegralSolver::IntegralSolver()
1051 : parameters_(), variable_values_(), objective_value_(0.0) {}
1063 const LinearProgram& linear_problem,
1064 const DenseRow& user_provided_initial_solution) {
1072 const LinearProgram& linear_problem,
1075 DenseRow initial_solution = user_provided_initial_solution;
1076 if (initial_solution.size() > 0) {
1077 CHECK_EQ(initial_solution.size(), linear_problem.num_variables())
1078 <<
"The initial solution should have the same number of variables as "
1079 "the LinearProgram.";
1084 LinearProgram
const* lp = &linear_problem;
1087 if (lp->num_variables() >= parameters_.decomposer_num_variables_threshold()) {
1088 LPDecomposer decomposer;
1089 decomposer.Decompose(lp);
1090 const int num_sub_problems = decomposer.GetNumberOfProblems();
1091 VLOG(1) <<
"Problem is decomposable into " << num_sub_problems
1093 if (num_sub_problems > 1) {
1097 std::vector<Fractional> objective_values(num_sub_problems,
1099 std::vector<Fractional> best_bounds(num_sub_problems,
Fractional(0.0));
1100 std::vector<BopSolveStatus> statuses(num_sub_problems,
1103 for (
int i = 0; i < num_sub_problems; ++i) {
1104 RunOneBop(parameters_, i, initial_solution,
time_limit, &decomposer,
1106 &(best_bounds[i]), &(statuses[i]));
1111 objective_value_ = lp->objective_offset();
1113 for (
int i = 0; i < num_sub_problems; ++i) {
1114 objective_value_ += objective_values[i];
1115 best_bound_ += best_bounds[i];
1130 InternalSolve(*lp, parameters_, initial_solution,
time_limit,
1131 &variable_values_, &objective_value_, &best_bound_);
1134 status = InternalSolve(*lp, parameters_, initial_solution,
time_limit,
1135 &variable_values_, &objective_value_, &best_bound_);
#define CHECK_EQ(val1, val2)
#define CHECK_GT(val1, val2)
#define CHECK_NE(val1, val2)
#define DCHECK(condition)
#define CHECK_LE(val1, val2)
#define VLOG(verboselevel)
A simple class to enforce both an elapsed time limit and a deterministic time limit in the same threa...
static std::unique_ptr< TimeLimit > FromParameters(const Parameters ¶meters)
Creates a time limit object initialized from an object that provides methods max_time_in_seconds() an...
ABSL_MUST_USE_RESULT BopSolveStatus SolveWithTimeLimit(const glop::LinearProgram &linear_problem, TimeLimit *time_limit)
const glop::DenseRow & variable_values() const
ABSL_MUST_USE_RESULT BopSolveStatus Solve(const glop::LinearProgram &linear_problem)
SharedTimeLimit * time_limit
@ FEASIBLE_SOLUTION_FOUND
bool CheckSolution(const Model &model, const std::function< int64(IntegerVariable *)> &evaluator)
StrictITIVector< ColIndex, Fractional > DenseRow
ColIndex RowToColIndex(RowIndex row)
StrictITIVector< RowIndex, Fractional > DenseColumn
void ChangeOptimizationDirection(LinearBooleanProblem *problem)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
bool IsIntegerWithinTolerance(FloatType x, FloatType tolerance)
int MostSignificantBitPosition64(uint64 n)
int64 ComputeGcdOfRoundedDoubles(const std::vector< double > &x, double scaling_factor)
double GetBestScalingOfDoublesToInt64(const std::vector< double > &input, const std::vector< double > &lb, const std::vector< double > &ub, int64 max_absolute_sum)
double max_scaling_factor
std::vector< double > coefficients