1146 lines
46 KiB
C++
1146 lines
46 KiB
C++
// Copyright 2010-2022 Google LLC
|
|
// Licensed under the Apache License, Version 2.0 (the "License");
|
|
// you may not use this file except in compliance with the License.
|
|
// You may obtain a copy of the License at
|
|
//
|
|
// http://www.apache.org/licenses/LICENSE-2.0
|
|
//
|
|
// Unless required by applicable law or agreed to in writing, software
|
|
// distributed under the License is distributed on an "AS IS" BASIS,
|
|
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
|
// See the License for the specific language governing permissions and
|
|
// limitations under the License.
|
|
|
|
#include "ortools/bop/integral_solver.h"
|
|
|
|
#include <math.h>
|
|
|
|
#include <algorithm>
|
|
#include <cstdint>
|
|
#include <limits>
|
|
#include <memory>
|
|
#include <string>
|
|
#include <vector>
|
|
|
|
#include "ortools/bop/bop_solver.h"
|
|
#include "ortools/lp_data/lp_decomposer.h"
|
|
|
|
namespace operations_research {
|
|
namespace bop {
|
|
|
|
using ::operations_research::glop::ColIndex;
|
|
using ::operations_research::glop::DenseRow;
|
|
using ::operations_research::glop::Fractional;
|
|
using ::operations_research::glop::kInfinity;
|
|
using ::operations_research::glop::LinearProgram;
|
|
using ::operations_research::glop::LPDecomposer;
|
|
using ::operations_research::glop::RowIndex;
|
|
using ::operations_research::glop::SparseColumn;
|
|
using ::operations_research::glop::SparseMatrix;
|
|
using ::operations_research::sat::LinearBooleanConstraint;
|
|
using ::operations_research::sat::LinearBooleanProblem;
|
|
using ::operations_research::sat::LinearObjective;
|
|
|
|
namespace {
|
|
// TODO(user): Use an existing one or move it to util.
|
|
bool IsIntegerWithinTolerance(Fractional x) {
|
|
const double kTolerance = 1e-10;
|
|
return std::abs(x - round(x)) <= kTolerance;
|
|
}
|
|
|
|
// Returns true when all the variables of the problem are Boolean, and all the
|
|
// constraints have integer coefficients.
|
|
// TODO(user): Move to SAT util.
|
|
bool ProblemIsBooleanAndHasOnlyIntegralConstraints(
|
|
const LinearProgram& linear_problem) {
|
|
const glop::SparseMatrix& matrix = linear_problem.GetSparseMatrix();
|
|
|
|
for (ColIndex col(0); col < linear_problem.num_variables(); ++col) {
|
|
const Fractional lower_bound = linear_problem.variable_lower_bounds()[col];
|
|
const Fractional upper_bound = linear_problem.variable_upper_bounds()[col];
|
|
|
|
if (lower_bound <= -1.0 || upper_bound >= 2.0) {
|
|
// Integral variable.
|
|
return false;
|
|
}
|
|
|
|
for (const SparseColumn::Entry e : matrix.column(col)) {
|
|
if (!IsIntegerWithinTolerance(e.coefficient())) {
|
|
// Floating coefficient.
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
// Builds a LinearBooleanProblem based on a LinearProgram with all the variables
|
|
// being booleans and all the constraints having only integral coefficients.
|
|
// TODO(user): Move to SAT util.
|
|
void BuildBooleanProblemWithIntegralConstraints(
|
|
const LinearProgram& linear_problem, const DenseRow& initial_solution,
|
|
LinearBooleanProblem* boolean_problem,
|
|
std::vector<bool>* boolean_initial_solution) {
|
|
CHECK(boolean_problem != nullptr);
|
|
boolean_problem->Clear();
|
|
|
|
const glop::SparseMatrix& matrix = linear_problem.GetSparseMatrix();
|
|
// Create Boolean variables.
|
|
for (ColIndex col(0); col < matrix.num_cols(); ++col) {
|
|
boolean_problem->add_var_names(linear_problem.GetVariableName(col));
|
|
}
|
|
boolean_problem->set_num_variables(matrix.num_cols().value());
|
|
boolean_problem->set_name(linear_problem.name());
|
|
|
|
// Create constraints.
|
|
for (RowIndex row(0); row < matrix.num_rows(); ++row) {
|
|
LinearBooleanConstraint* const constraint =
|
|
boolean_problem->add_constraints();
|
|
constraint->set_name(linear_problem.GetConstraintName(row));
|
|
if (linear_problem.constraint_lower_bounds()[row] != -kInfinity) {
|
|
constraint->set_lower_bound(
|
|
linear_problem.constraint_lower_bounds()[row]);
|
|
}
|
|
if (linear_problem.constraint_upper_bounds()[row] != kInfinity) {
|
|
constraint->set_upper_bound(
|
|
linear_problem.constraint_upper_bounds()[row]);
|
|
}
|
|
}
|
|
|
|
// Store the constraint coefficients.
|
|
for (ColIndex col(0); col < matrix.num_cols(); ++col) {
|
|
for (const SparseColumn::Entry e : matrix.column(col)) {
|
|
LinearBooleanConstraint* const constraint =
|
|
boolean_problem->mutable_constraints(e.row().value());
|
|
constraint->add_literals(col.value() + 1);
|
|
constraint->add_coefficients(e.coefficient());
|
|
}
|
|
}
|
|
|
|
// Add the unit constraints to fix the variables since the variable bounds
|
|
// are always [0, 1] in a BooleanLinearProblem.
|
|
for (ColIndex col(0); col < matrix.num_cols(); ++col) {
|
|
// TODO(user): double check the rounding, and add unit test for this.
|
|
const int lb = std::round(linear_problem.variable_lower_bounds()[col]);
|
|
const int ub = std::round(linear_problem.variable_upper_bounds()[col]);
|
|
if (lb == ub) {
|
|
LinearBooleanConstraint* ct = boolean_problem->add_constraints();
|
|
ct->set_lower_bound(ub);
|
|
ct->set_upper_bound(ub);
|
|
ct->add_literals(col.value() + 1);
|
|
ct->add_coefficients(1.0);
|
|
}
|
|
}
|
|
|
|
// Create the minimization objective.
|
|
std::vector<double> coefficients;
|
|
for (ColIndex col(0); col < linear_problem.num_variables(); ++col) {
|
|
const Fractional coeff = linear_problem.objective_coefficients()[col];
|
|
if (coeff != 0.0) coefficients.push_back(coeff);
|
|
}
|
|
double scaling_factor = 0.0;
|
|
double relative_error = 0.0;
|
|
GetBestScalingOfDoublesToInt64(coefficients,
|
|
std::numeric_limits<int64_t>::max(),
|
|
&scaling_factor, &relative_error);
|
|
const int64_t gcd = ComputeGcdOfRoundedDoubles(coefficients, scaling_factor);
|
|
LinearObjective* const objective = boolean_problem->mutable_objective();
|
|
objective->set_offset(linear_problem.objective_offset() * scaling_factor /
|
|
gcd);
|
|
|
|
// Note that here we set the scaling factor for the inverse operation of
|
|
// getting the "true" objective value from the scaled one. Hence the inverse.
|
|
objective->set_scaling_factor(1.0 / scaling_factor * gcd);
|
|
for (ColIndex col(0); col < linear_problem.num_variables(); ++col) {
|
|
const Fractional coeff = linear_problem.objective_coefficients()[col];
|
|
const int64_t value =
|
|
static_cast<int64_t>(round(coeff * scaling_factor)) / gcd;
|
|
if (value != 0) {
|
|
objective->add_literals(col.value() + 1);
|
|
objective->add_coefficients(value);
|
|
}
|
|
}
|
|
|
|
// If the problem was a maximization one, we need to modify the objective.
|
|
if (linear_problem.IsMaximizationProblem()) {
|
|
sat::ChangeOptimizationDirection(boolean_problem);
|
|
}
|
|
|
|
// Fill the Boolean initial solution.
|
|
if (!initial_solution.empty()) {
|
|
CHECK(boolean_initial_solution != nullptr);
|
|
CHECK_EQ(boolean_problem->num_variables(), initial_solution.size());
|
|
boolean_initial_solution->assign(boolean_problem->num_variables(), false);
|
|
for (int i = 0; i < initial_solution.size(); ++i) {
|
|
(*boolean_initial_solution)[i] = (initial_solution[ColIndex(i)] != 0);
|
|
}
|
|
}
|
|
}
|
|
|
|
//------------------------------------------------------------------------------
|
|
// IntegralVariable
|
|
//------------------------------------------------------------------------------
|
|
// Model an integral variable using Boolean variables.
|
|
// TODO(user): Enable discrete representation by value, i.e. use three Boolean
|
|
// variables when only possible values are 10, 12, 32.
|
|
// In the same way, when only two consecutive values are possible
|
|
// use only one Boolean variable with an offset.
|
|
class IntegralVariable {
|
|
public:
|
|
IntegralVariable();
|
|
|
|
// Creates the minimal number of Boolean variables to represent an integral
|
|
// variable with range [lower_bound, upper_bound]. start_var_index corresponds
|
|
// to the next available Boolean variable index. If three Boolean variables
|
|
// are needed to model the integral variable, the used variables will have
|
|
// indices start_var_index, start_var_index +1, and start_var_index +2.
|
|
void BuildFromRange(int start_var_index, Fractional lower_bound,
|
|
Fractional upper_bound);
|
|
|
|
void Clear();
|
|
void set_offset(int64_t offset) { offset_ = offset; }
|
|
void set_weight(VariableIndex var, int64_t weight);
|
|
|
|
int GetNumberOfBooleanVariables() const { return bits_.size(); }
|
|
|
|
const std::vector<VariableIndex>& bits() const { return bits_; }
|
|
const std::vector<int64_t>& weights() const { return weights_; }
|
|
int64_t offset() const { return offset_; }
|
|
|
|
// Returns the value of the integral variable based on the Boolean conversion
|
|
// and the Boolean solution to the problem.
|
|
int64_t GetSolutionValue(const BopSolution& solution) const;
|
|
|
|
// Returns the values of the Boolean variables based on the Boolean conversion
|
|
// and the integral value of this variable. This only works for variables that
|
|
// were constructed using BuildFromRange() (for which can_be_reversed_ is
|
|
// true).
|
|
std::vector<bool> GetBooleanSolutionValues(int64_t integral_value) const;
|
|
|
|
std::string DebugString() const;
|
|
|
|
private:
|
|
// The value of the integral variable is expressed as
|
|
// sum_i(weights[i] * Value(bits[i])) + offset.
|
|
// Note that weights can be negative to represent negative values.
|
|
std::vector<VariableIndex> bits_;
|
|
std::vector<int64_t> weights_;
|
|
int64_t offset_;
|
|
// True if the values of the boolean variables representing this integral
|
|
// variable can be deduced from the integral variable's value. Namely, this is
|
|
// true for variables built using BuildFromRange() but usually false for
|
|
// variables built using set_weight().
|
|
bool can_be_reversed_;
|
|
};
|
|
|
|
IntegralVariable::IntegralVariable()
|
|
: bits_(), weights_(), offset_(0), can_be_reversed_(true) {}
|
|
|
|
void IntegralVariable::BuildFromRange(int start_var_index,
|
|
Fractional lower_bound,
|
|
Fractional upper_bound) {
|
|
Clear();
|
|
|
|
// Integral variable. Split the variable into the minimum number of bits
|
|
// required to model the upper bound.
|
|
CHECK_NE(-kInfinity, lower_bound);
|
|
CHECK_NE(kInfinity, upper_bound);
|
|
|
|
const int64_t integral_lower_bound = static_cast<int64_t>(ceil(lower_bound));
|
|
const int64_t integral_upper_bound = static_cast<int64_t>(floor(upper_bound));
|
|
offset_ = integral_lower_bound;
|
|
const int64_t delta = integral_upper_bound - integral_lower_bound;
|
|
const int num_used_bits = MostSignificantBitPosition64(delta) + 1;
|
|
for (int i = 0; i < num_used_bits; ++i) {
|
|
bits_.push_back(VariableIndex(start_var_index + i));
|
|
weights_.push_back(1ULL << i);
|
|
}
|
|
}
|
|
|
|
void IntegralVariable::Clear() {
|
|
bits_.clear();
|
|
weights_.clear();
|
|
offset_ = 0;
|
|
can_be_reversed_ = true;
|
|
}
|
|
|
|
void IntegralVariable::set_weight(VariableIndex var, int64_t weight) {
|
|
bits_.push_back(var);
|
|
weights_.push_back(weight);
|
|
can_be_reversed_ = false;
|
|
}
|
|
|
|
int64_t IntegralVariable::GetSolutionValue(const BopSolution& solution) const {
|
|
int64_t value = offset_;
|
|
for (int i = 0; i < bits_.size(); ++i) {
|
|
value += weights_[i] * solution.Value(bits_[i]);
|
|
}
|
|
return value;
|
|
}
|
|
|
|
std::vector<bool> IntegralVariable::GetBooleanSolutionValues(
|
|
int64_t integral_value) const {
|
|
if (can_be_reversed_) {
|
|
DCHECK(std::is_sorted(weights_.begin(), weights_.end()));
|
|
std::vector<bool> boolean_values(weights_.size(), false);
|
|
int64_t remaining_value = integral_value - offset_;
|
|
for (int i = weights_.size() - 1; i >= 0; --i) {
|
|
if (remaining_value >= weights_[i]) {
|
|
boolean_values[i] = true;
|
|
remaining_value -= weights_[i];
|
|
}
|
|
}
|
|
CHECK_EQ(0, remaining_value)
|
|
<< "Couldn't map integral value to boolean variables.";
|
|
return boolean_values;
|
|
}
|
|
return std::vector<bool>();
|
|
}
|
|
|
|
std::string IntegralVariable::DebugString() const {
|
|
std::string str;
|
|
CHECK_EQ(bits_.size(), weights_.size());
|
|
for (int i = 0; i < bits_.size(); ++i) {
|
|
str += absl::StrFormat("%d [%d] ", weights_[i], bits_[i].value());
|
|
}
|
|
str += absl::StrFormat(" Offset: %d", offset_);
|
|
return str;
|
|
}
|
|
|
|
//------------------------------------------------------------------------------
|
|
// IntegralProblemConverter
|
|
//------------------------------------------------------------------------------
|
|
// This class is used to convert a LinearProblem containing integral variables
|
|
// into a LinearBooleanProblem that Bop can consume.
|
|
// The converter tries to reuse existing Boolean variables as much as possible,
|
|
// but there are no guarantees to model all integral variables using the total
|
|
// minimal number of Boolean variables.
|
|
// Consider for instance the constraint "x - 2 * y = 0".
|
|
// Depending on the declaration order, two different outcomes are possible:
|
|
// - When x is considered first, the converter will generate new variables
|
|
// for both x and y as we only consider integral weights, i.e. y = x / 2.
|
|
// - When y is considered first, the converter will reuse Boolean variables
|
|
// from y to model x as x = 2 * y (integral weight).
|
|
//
|
|
// Note that the converter only deals with integral variables, i.e. no
|
|
// continuous variables.
|
|
class IntegralProblemConverter {
|
|
public:
|
|
IntegralProblemConverter();
|
|
|
|
// Converts the LinearProgram into a LinearBooleanProblem. If an initial
|
|
// solution is given (i.e. if its size is not zero), converts it into a
|
|
// Boolean solution.
|
|
// Returns false when the conversion fails.
|
|
bool ConvertToBooleanProblem(const LinearProgram& linear_problem,
|
|
const DenseRow& initial_solution,
|
|
LinearBooleanProblem* boolean_problem,
|
|
std::vector<bool>* boolean_initial_solution);
|
|
|
|
// Returns the value of a variable of the original problem based on the
|
|
// Boolean conversion and the Boolean solution to the problem.
|
|
int64_t GetSolutionValue(ColIndex global_col,
|
|
const BopSolution& solution) const;
|
|
|
|
private:
|
|
// Returns true when the linear_problem_ can be converted into a Boolean
|
|
// problem. Note that floating weights and continuous variables are not
|
|
// supported.
|
|
bool CheckProblem(const LinearProgram& linear_problem) const;
|
|
|
|
// Initializes the type of each variable of the linear_problem_.
|
|
void InitVariableTypes(const LinearProgram& linear_problem,
|
|
LinearBooleanProblem* boolean_problem);
|
|
|
|
// Converts all variables of the problem.
|
|
void ConvertAllVariables(const LinearProgram& linear_problem,
|
|
LinearBooleanProblem* boolean_problem);
|
|
|
|
// Adds all variables constraints, i.e. lower and upper bounds of variables.
|
|
void AddVariableConstraints(const LinearProgram& linear_problem,
|
|
LinearBooleanProblem* boolean_problem);
|
|
|
|
// Converts all constraints from LinearProgram to LinearBooleanProblem.
|
|
void ConvertAllConstraints(const LinearProgram& linear_problem,
|
|
LinearBooleanProblem* boolean_problem);
|
|
|
|
// Converts the objective from LinearProgram to LinearBooleanProblem.
|
|
void ConvertObjective(const LinearProgram& linear_problem,
|
|
LinearBooleanProblem* boolean_problem);
|
|
|
|
// Converts the integral variable represented by col in the linear_problem_
|
|
// into an IntegralVariable using existing Boolean variables.
|
|
// Returns false when existing Boolean variables are not enough to model
|
|
// the integral variable.
|
|
bool ConvertUsingExistingBooleans(const LinearProgram& linear_problem,
|
|
ColIndex col,
|
|
IntegralVariable* integral_var);
|
|
|
|
// Creates the integral_var using the given linear_problem_ constraint.
|
|
// The constraint is an equality constraint and contains only one integral
|
|
// variable (already the case in the model or thanks to previous
|
|
// booleanization of other integral variables), i.e.
|
|
// bound <= w * integral_var + sum(w_i * b_i) <= bound
|
|
// The remaining integral variable can then be expressed:
|
|
// integral_var == (bound + sum(-w_i * b_i)) / w
|
|
// Note that all divisions by w have to be integral as Bop only deals with
|
|
// integral coefficients.
|
|
bool CreateVariableUsingConstraint(const LinearProgram& linear_problem,
|
|
RowIndex constraint,
|
|
IntegralVariable* integral_var);
|
|
|
|
// Adds weighted integral variable represented by col to the current dense
|
|
// constraint.
|
|
Fractional AddWeightedIntegralVariable(
|
|
ColIndex col, Fractional weight,
|
|
absl::StrongVector<VariableIndex, Fractional>* dense_weights);
|
|
|
|
// Scales weights and adds all non-zero scaled weights and literals to t.
|
|
// t is a constraint or the objective.
|
|
// Returns the bound error due to the scaling.
|
|
// The weight is scaled using:
|
|
// static_cast<int64_t>(round(weight * scaling_factor)) / gcd;
|
|
template <class T>
|
|
double ScaleAndSparsifyWeights(
|
|
double scaling_factor, int64_t gcd,
|
|
const absl::StrongVector<VariableIndex, Fractional>& dense_weights, T* t);
|
|
|
|
// Returns true when at least one element is non-zero.
|
|
bool HasNonZeroWeights(
|
|
const absl::StrongVector<VariableIndex, Fractional>& dense_weights) const;
|
|
|
|
bool problem_is_boolean_and_has_only_integral_constraints_;
|
|
|
|
// global_to_boolean_[i] represents the Boolean variable index in Bop; when
|
|
// negative -global_to_boolean_[i] - 1 represents the index of the
|
|
// integral variable in integral_variables_.
|
|
absl::StrongVector</*global_col*/ glop::ColIndex, /*boolean_col*/ int>
|
|
global_to_boolean_;
|
|
std::vector<IntegralVariable> integral_variables_;
|
|
std::vector<ColIndex> integral_indices_;
|
|
int num_boolean_variables_;
|
|
|
|
enum VariableType { BOOLEAN, INTEGRAL, INTEGRAL_EXPRESSED_AS_BOOLEAN };
|
|
absl::StrongVector<glop::ColIndex, VariableType> variable_types_;
|
|
};
|
|
|
|
IntegralProblemConverter::IntegralProblemConverter()
|
|
: global_to_boolean_(),
|
|
integral_variables_(),
|
|
integral_indices_(),
|
|
num_boolean_variables_(0),
|
|
variable_types_() {}
|
|
|
|
bool IntegralProblemConverter::ConvertToBooleanProblem(
|
|
const LinearProgram& linear_problem, const DenseRow& initial_solution,
|
|
LinearBooleanProblem* boolean_problem,
|
|
std::vector<bool>* boolean_initial_solution) {
|
|
bool use_initial_solution = (initial_solution.size() > 0);
|
|
if (use_initial_solution) {
|
|
CHECK_EQ(initial_solution.size(), linear_problem.num_variables())
|
|
<< "The initial solution should have the same number of variables as "
|
|
"the LinearProgram.";
|
|
CHECK(boolean_initial_solution != nullptr);
|
|
}
|
|
if (!CheckProblem(linear_problem)) {
|
|
return false;
|
|
}
|
|
|
|
problem_is_boolean_and_has_only_integral_constraints_ =
|
|
ProblemIsBooleanAndHasOnlyIntegralConstraints(linear_problem);
|
|
if (problem_is_boolean_and_has_only_integral_constraints_) {
|
|
BuildBooleanProblemWithIntegralConstraints(linear_problem, initial_solution,
|
|
boolean_problem,
|
|
boolean_initial_solution);
|
|
return true;
|
|
}
|
|
|
|
InitVariableTypes(linear_problem, boolean_problem);
|
|
ConvertAllVariables(linear_problem, boolean_problem);
|
|
boolean_problem->set_num_variables(num_boolean_variables_);
|
|
boolean_problem->set_name(linear_problem.name());
|
|
|
|
AddVariableConstraints(linear_problem, boolean_problem);
|
|
ConvertAllConstraints(linear_problem, boolean_problem);
|
|
ConvertObjective(linear_problem, boolean_problem);
|
|
|
|
// A BooleanLinearProblem is always in the minimization form.
|
|
if (linear_problem.IsMaximizationProblem()) {
|
|
sat::ChangeOptimizationDirection(boolean_problem);
|
|
}
|
|
|
|
if (use_initial_solution) {
|
|
boolean_initial_solution->assign(boolean_problem->num_variables(), false);
|
|
for (ColIndex global_col(0); global_col < global_to_boolean_.size();
|
|
++global_col) {
|
|
const int col = global_to_boolean_[global_col];
|
|
if (col >= 0) {
|
|
(*boolean_initial_solution)[col] = (initial_solution[global_col] != 0);
|
|
} else {
|
|
const IntegralVariable& integral_variable =
|
|
integral_variables_[-col - 1];
|
|
const std::vector<VariableIndex>& boolean_cols =
|
|
integral_variable.bits();
|
|
const std::vector<bool>& boolean_values =
|
|
integral_variable.GetBooleanSolutionValues(
|
|
round(initial_solution[global_col]));
|
|
if (!boolean_values.empty()) {
|
|
CHECK_EQ(boolean_cols.size(), boolean_values.size());
|
|
for (int i = 0; i < boolean_values.size(); ++i) {
|
|
const int boolean_col = boolean_cols[i].value();
|
|
(*boolean_initial_solution)[boolean_col] = boolean_values[i];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
int64_t IntegralProblemConverter::GetSolutionValue(
|
|
ColIndex global_col, const BopSolution& solution) const {
|
|
if (problem_is_boolean_and_has_only_integral_constraints_) {
|
|
return solution.Value(VariableIndex(global_col.value()));
|
|
}
|
|
|
|
const int pos = global_to_boolean_[global_col];
|
|
return pos >= 0 ? solution.Value(VariableIndex(pos))
|
|
: integral_variables_[-pos - 1].GetSolutionValue(solution);
|
|
}
|
|
|
|
bool IntegralProblemConverter::CheckProblem(
|
|
const LinearProgram& linear_problem) const {
|
|
for (ColIndex col(0); col < linear_problem.num_variables(); ++col) {
|
|
if (!linear_problem.IsVariableInteger(col)) {
|
|
LOG(ERROR) << "Variable " << linear_problem.GetVariableName(col)
|
|
<< " is continuous. This is not supported by BOP.";
|
|
return false;
|
|
}
|
|
if (linear_problem.variable_lower_bounds()[col] == -kInfinity) {
|
|
LOG(ERROR) << "Variable " << linear_problem.GetVariableName(col)
|
|
<< " has no lower bound. This is not supported by BOP.";
|
|
return false;
|
|
}
|
|
if (linear_problem.variable_upper_bounds()[col] == kInfinity) {
|
|
LOG(ERROR) << "Variable " << linear_problem.GetVariableName(col)
|
|
<< " has no upper bound. This is not supported by BOP.";
|
|
return false;
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
void IntegralProblemConverter::InitVariableTypes(
|
|
const LinearProgram& linear_problem,
|
|
LinearBooleanProblem* boolean_problem) {
|
|
global_to_boolean_.assign(linear_problem.num_variables().value(), 0);
|
|
variable_types_.assign(linear_problem.num_variables().value(), INTEGRAL);
|
|
for (ColIndex col(0); col < linear_problem.num_variables(); ++col) {
|
|
const Fractional lower_bound = linear_problem.variable_lower_bounds()[col];
|
|
const Fractional upper_bound = linear_problem.variable_upper_bounds()[col];
|
|
|
|
if (lower_bound > -1.0 && upper_bound < 2.0) {
|
|
// Boolean variable.
|
|
variable_types_[col] = BOOLEAN;
|
|
global_to_boolean_[col] = num_boolean_variables_;
|
|
++num_boolean_variables_;
|
|
boolean_problem->add_var_names(linear_problem.GetVariableName(col));
|
|
} else {
|
|
// Integral variable.
|
|
variable_types_[col] = INTEGRAL;
|
|
integral_indices_.push_back(col);
|
|
}
|
|
}
|
|
}
|
|
|
|
void IntegralProblemConverter::ConvertAllVariables(
|
|
const LinearProgram& linear_problem,
|
|
LinearBooleanProblem* boolean_problem) {
|
|
for (const ColIndex col : integral_indices_) {
|
|
CHECK_EQ(INTEGRAL, variable_types_[col]);
|
|
IntegralVariable integral_var;
|
|
if (!ConvertUsingExistingBooleans(linear_problem, col, &integral_var)) {
|
|
const Fractional lower_bound =
|
|
linear_problem.variable_lower_bounds()[col];
|
|
const Fractional upper_bound =
|
|
linear_problem.variable_upper_bounds()[col];
|
|
integral_var.BuildFromRange(num_boolean_variables_, lower_bound,
|
|
upper_bound);
|
|
num_boolean_variables_ += integral_var.GetNumberOfBooleanVariables();
|
|
const std::string var_name = linear_problem.GetVariableName(col);
|
|
for (int i = 0; i < integral_var.bits().size(); ++i) {
|
|
boolean_problem->add_var_names(var_name + absl::StrFormat("_%d", i));
|
|
}
|
|
}
|
|
integral_variables_.push_back(integral_var);
|
|
global_to_boolean_[col] = -integral_variables_.size();
|
|
variable_types_[col] = INTEGRAL_EXPRESSED_AS_BOOLEAN;
|
|
}
|
|
}
|
|
|
|
void IntegralProblemConverter::ConvertAllConstraints(
|
|
const LinearProgram& linear_problem,
|
|
LinearBooleanProblem* boolean_problem) {
|
|
// TODO(user): This is the way it's done in glop/proto_utils.cc but having
|
|
// to transpose looks unnecessary costly.
|
|
glop::SparseMatrix transpose;
|
|
transpose.PopulateFromTranspose(linear_problem.GetSparseMatrix());
|
|
|
|
double max_relative_error = 0.0;
|
|
double max_bound_error = 0.0;
|
|
double max_scaling_factor = 0.0;
|
|
double relative_error = 0.0;
|
|
double scaling_factor = 0.0;
|
|
std::vector<double> coefficients;
|
|
for (RowIndex row(0); row < linear_problem.num_constraints(); ++row) {
|
|
Fractional offset = 0.0;
|
|
absl::StrongVector<VariableIndex, Fractional> dense_weights(
|
|
num_boolean_variables_, 0.0);
|
|
for (const SparseColumn::Entry e : transpose.column(RowToColIndex(row))) {
|
|
// Cast in ColIndex due to the transpose.
|
|
offset += AddWeightedIntegralVariable(RowToColIndex(e.row()),
|
|
e.coefficient(), &dense_weights);
|
|
}
|
|
if (!HasNonZeroWeights(dense_weights)) {
|
|
continue;
|
|
}
|
|
|
|
// Compute the scaling for non-integral weights.
|
|
coefficients.clear();
|
|
for (VariableIndex var(0); var < num_boolean_variables_; ++var) {
|
|
if (dense_weights[var] != 0.0) {
|
|
coefficients.push_back(dense_weights[var]);
|
|
}
|
|
}
|
|
GetBestScalingOfDoublesToInt64(coefficients,
|
|
std::numeric_limits<int64_t>::max(),
|
|
&scaling_factor, &relative_error);
|
|
const int64_t gcd =
|
|
ComputeGcdOfRoundedDoubles(coefficients, scaling_factor);
|
|
max_relative_error = std::max(relative_error, max_relative_error);
|
|
max_scaling_factor = std::max(scaling_factor / gcd, max_scaling_factor);
|
|
|
|
LinearBooleanConstraint* constraint = boolean_problem->add_constraints();
|
|
constraint->set_name(linear_problem.GetConstraintName(row));
|
|
const double bound_error =
|
|
ScaleAndSparsifyWeights(scaling_factor, gcd, dense_weights, constraint);
|
|
max_bound_error = std::max(max_bound_error, bound_error);
|
|
|
|
const Fractional lower_bound =
|
|
linear_problem.constraint_lower_bounds()[row];
|
|
if (lower_bound != -kInfinity) {
|
|
const Fractional offset_lower_bound = lower_bound - offset;
|
|
const double offset_scaled_lower_bound =
|
|
round(offset_lower_bound * scaling_factor - bound_error);
|
|
if (offset_scaled_lower_bound >=
|
|
static_cast<double>(std::numeric_limits<int64_t>::max())) {
|
|
LOG(WARNING) << "A constraint is trivially unsatisfiable.";
|
|
return;
|
|
}
|
|
if (offset_scaled_lower_bound >
|
|
-static_cast<double>(std::numeric_limits<int64_t>::max())) {
|
|
// Otherwise, the constraint is not needed.
|
|
constraint->set_lower_bound(
|
|
static_cast<int64_t>(offset_scaled_lower_bound) / gcd);
|
|
}
|
|
}
|
|
const Fractional upper_bound =
|
|
linear_problem.constraint_upper_bounds()[row];
|
|
if (upper_bound != kInfinity) {
|
|
const Fractional offset_upper_bound = upper_bound - offset;
|
|
const double offset_scaled_upper_bound =
|
|
round(offset_upper_bound * scaling_factor + bound_error);
|
|
if (offset_scaled_upper_bound <=
|
|
-static_cast<double>(std::numeric_limits<int64_t>::max())) {
|
|
LOG(WARNING) << "A constraint is trivially unsatisfiable.";
|
|
return;
|
|
}
|
|
if (offset_scaled_upper_bound <
|
|
static_cast<double>(std::numeric_limits<int64_t>::max())) {
|
|
// Otherwise, the constraint is not needed.
|
|
constraint->set_upper_bound(
|
|
static_cast<int64_t>(offset_scaled_upper_bound) / gcd);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void IntegralProblemConverter::ConvertObjective(
|
|
const LinearProgram& linear_problem,
|
|
LinearBooleanProblem* boolean_problem) {
|
|
LinearObjective* objective = boolean_problem->mutable_objective();
|
|
Fractional offset = 0.0;
|
|
absl::StrongVector<VariableIndex, Fractional> dense_weights(
|
|
num_boolean_variables_, 0.0);
|
|
// Compute the objective weights for the binary variable model.
|
|
for (ColIndex col(0); col < linear_problem.num_variables(); ++col) {
|
|
offset += AddWeightedIntegralVariable(
|
|
col, linear_problem.objective_coefficients()[col], &dense_weights);
|
|
}
|
|
|
|
// Compute the scaling for non-integral weights.
|
|
std::vector<double> coefficients;
|
|
for (VariableIndex var(0); var < num_boolean_variables_; ++var) {
|
|
if (dense_weights[var] != 0.0) {
|
|
coefficients.push_back(dense_weights[var]);
|
|
}
|
|
}
|
|
double scaling_factor = 0.0;
|
|
double max_relative_error = 0.0;
|
|
double relative_error = 0.0;
|
|
GetBestScalingOfDoublesToInt64(coefficients,
|
|
std::numeric_limits<int64_t>::max(),
|
|
&scaling_factor, &relative_error);
|
|
const int64_t gcd = ComputeGcdOfRoundedDoubles(coefficients, scaling_factor);
|
|
max_relative_error = std::max(relative_error, max_relative_error);
|
|
VLOG(1) << "objective relative error: " << relative_error;
|
|
VLOG(1) << "objective scaling factor: " << scaling_factor / gcd;
|
|
|
|
ScaleAndSparsifyWeights(scaling_factor, gcd, dense_weights, objective);
|
|
|
|
// Note that here we set the scaling factor for the inverse operation of
|
|
// getting the "true" objective value from the scaled one. Hence the inverse.
|
|
objective->set_scaling_factor(1.0 / scaling_factor * gcd);
|
|
objective->set_offset((linear_problem.objective_offset() + offset) *
|
|
scaling_factor / gcd);
|
|
}
|
|
|
|
void IntegralProblemConverter::AddVariableConstraints(
|
|
const LinearProgram& linear_problem,
|
|
LinearBooleanProblem* boolean_problem) {
|
|
for (ColIndex col(0); col < linear_problem.num_variables(); ++col) {
|
|
const Fractional lower_bound = linear_problem.variable_lower_bounds()[col];
|
|
const Fractional upper_bound = linear_problem.variable_upper_bounds()[col];
|
|
const int pos = global_to_boolean_[col];
|
|
if (pos >= 0) {
|
|
// Boolean variable.
|
|
CHECK_EQ(BOOLEAN, variable_types_[col]);
|
|
const bool is_fixed = (lower_bound > -1.0 && upper_bound < 1.0) ||
|
|
(lower_bound > 0.0 && upper_bound < 2.0);
|
|
if (is_fixed) {
|
|
// Set the variable.
|
|
const int fixed_value = lower_bound > -1.0 && upper_bound < 1.0 ? 0 : 1;
|
|
LinearBooleanConstraint* constraint =
|
|
boolean_problem->add_constraints();
|
|
constraint->set_lower_bound(fixed_value);
|
|
constraint->set_upper_bound(fixed_value);
|
|
constraint->add_literals(pos + 1);
|
|
constraint->add_coefficients(1);
|
|
}
|
|
} else {
|
|
CHECK_EQ(INTEGRAL_EXPRESSED_AS_BOOLEAN, variable_types_[col]);
|
|
// Integral variable.
|
|
if (lower_bound != -kInfinity || upper_bound != kInfinity) {
|
|
const IntegralVariable& integral_var = integral_variables_[-pos - 1];
|
|
LinearBooleanConstraint* constraint =
|
|
boolean_problem->add_constraints();
|
|
for (int i = 0; i < integral_var.bits().size(); ++i) {
|
|
constraint->add_literals(integral_var.bits()[i].value() + 1);
|
|
constraint->add_coefficients(integral_var.weights()[i]);
|
|
}
|
|
if (lower_bound != -kInfinity) {
|
|
constraint->set_lower_bound(static_cast<int64_t>(ceil(lower_bound)) -
|
|
integral_var.offset());
|
|
}
|
|
if (upper_bound != kInfinity) {
|
|
constraint->set_upper_bound(static_cast<int64_t>(floor(upper_bound)) -
|
|
integral_var.offset());
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
bool IntegralProblemConverter::ConvertUsingExistingBooleans(
|
|
const LinearProgram& linear_problem, ColIndex col,
|
|
IntegralVariable* integral_var) {
|
|
CHECK(nullptr != integral_var);
|
|
CHECK_EQ(INTEGRAL, variable_types_[col]);
|
|
|
|
const SparseMatrix& matrix = linear_problem.GetSparseMatrix();
|
|
const SparseMatrix& transpose = linear_problem.GetTransposeSparseMatrix();
|
|
for (const SparseColumn::Entry var_entry : matrix.column(col)) {
|
|
const RowIndex constraint = var_entry.row();
|
|
const Fractional lb = linear_problem.constraint_lower_bounds()[constraint];
|
|
const Fractional ub = linear_problem.constraint_upper_bounds()[constraint];
|
|
if (lb != ub) {
|
|
// To replace an integral variable by a weighted sum of Boolean variables,
|
|
// the constraint has to be an equality.
|
|
continue;
|
|
}
|
|
|
|
if (transpose.column(RowToColIndex(constraint)).num_entries() <= 1) {
|
|
// Can't replace the integer variable by Boolean variables when there are
|
|
// no Boolean variables.
|
|
// TODO(user): We could actually simplify the problem when the variable
|
|
// is constant, but this should be done by the preprocessor,
|
|
// not here. Consider activating the MIP preprocessing.
|
|
continue;
|
|
}
|
|
|
|
bool only_one_integral_variable = true;
|
|
for (const SparseColumn::Entry constraint_entry :
|
|
transpose.column(RowToColIndex(constraint))) {
|
|
const ColIndex var_index = RowToColIndex(constraint_entry.row());
|
|
if (var_index != col && variable_types_[var_index] == INTEGRAL) {
|
|
only_one_integral_variable = false;
|
|
break;
|
|
}
|
|
}
|
|
if (only_one_integral_variable &&
|
|
CreateVariableUsingConstraint(linear_problem, constraint,
|
|
integral_var)) {
|
|
return true;
|
|
}
|
|
}
|
|
|
|
integral_var->Clear();
|
|
return false;
|
|
}
|
|
|
|
bool IntegralProblemConverter::CreateVariableUsingConstraint(
|
|
const LinearProgram& linear_problem, RowIndex constraint,
|
|
IntegralVariable* integral_var) {
|
|
CHECK(nullptr != integral_var);
|
|
integral_var->Clear();
|
|
|
|
const SparseMatrix& transpose = linear_problem.GetTransposeSparseMatrix();
|
|
absl::StrongVector<VariableIndex, Fractional> dense_weights(
|
|
num_boolean_variables_, 0.0);
|
|
Fractional scale = 1.0;
|
|
int64_t variable_offset = 0;
|
|
for (const SparseColumn::Entry constraint_entry :
|
|
transpose.column(RowToColIndex(constraint))) {
|
|
const ColIndex col = RowToColIndex(constraint_entry.row());
|
|
if (variable_types_[col] == INTEGRAL) {
|
|
scale = constraint_entry.coefficient();
|
|
} else if (variable_types_[col] == BOOLEAN) {
|
|
const int pos = global_to_boolean_[col];
|
|
CHECK_LE(0, pos);
|
|
dense_weights[VariableIndex(pos)] -= constraint_entry.coefficient();
|
|
} else {
|
|
CHECK_EQ(INTEGRAL_EXPRESSED_AS_BOOLEAN, variable_types_[col]);
|
|
const int pos = global_to_boolean_[col];
|
|
CHECK_GT(0, pos);
|
|
const IntegralVariable& local_integral_var =
|
|
integral_variables_[-pos - 1];
|
|
variable_offset -=
|
|
constraint_entry.coefficient() * local_integral_var.offset();
|
|
for (int i = 0; i < local_integral_var.bits().size(); ++i) {
|
|
dense_weights[local_integral_var.bits()[i]] -=
|
|
constraint_entry.coefficient() * local_integral_var.weights()[i];
|
|
}
|
|
}
|
|
}
|
|
|
|
// Rescale using the weight of the integral variable.
|
|
const Fractional lb = linear_problem.constraint_lower_bounds()[constraint];
|
|
const Fractional offset = (lb + variable_offset) / scale;
|
|
if (!IsIntegerWithinTolerance(offset)) {
|
|
return false;
|
|
}
|
|
integral_var->set_offset(static_cast<int64_t>(offset));
|
|
|
|
for (VariableIndex var(0); var < dense_weights.size(); ++var) {
|
|
if (dense_weights[var] != 0.0) {
|
|
const Fractional weight = dense_weights[var] / scale;
|
|
if (!IsIntegerWithinTolerance(weight)) {
|
|
return false;
|
|
}
|
|
integral_var->set_weight(var, static_cast<int64_t>(weight));
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
Fractional IntegralProblemConverter::AddWeightedIntegralVariable(
|
|
ColIndex col, Fractional weight,
|
|
absl::StrongVector<VariableIndex, Fractional>* dense_weights) {
|
|
CHECK(nullptr != dense_weights);
|
|
|
|
if (weight == 0.0) {
|
|
return 0;
|
|
}
|
|
|
|
Fractional offset = 0;
|
|
const int pos = global_to_boolean_[col];
|
|
if (pos >= 0) {
|
|
// Boolean variable.
|
|
(*dense_weights)[VariableIndex(pos)] += weight;
|
|
} else {
|
|
// Integral variable.
|
|
const IntegralVariable& integral_var = integral_variables_[-pos - 1];
|
|
for (int i = 0; i < integral_var.bits().size(); ++i) {
|
|
(*dense_weights)[integral_var.bits()[i]] +=
|
|
integral_var.weights()[i] * weight;
|
|
}
|
|
offset += weight * integral_var.offset();
|
|
}
|
|
return offset;
|
|
}
|
|
|
|
template <class T>
|
|
double IntegralProblemConverter::ScaleAndSparsifyWeights(
|
|
double scaling_factor, int64_t gcd,
|
|
const absl::StrongVector<VariableIndex, Fractional>& dense_weights, T* t) {
|
|
CHECK(nullptr != t);
|
|
|
|
double bound_error = 0.0;
|
|
for (VariableIndex var(0); var < dense_weights.size(); ++var) {
|
|
if (dense_weights[var] != 0.0) {
|
|
const double scaled_weight = dense_weights[var] * scaling_factor;
|
|
bound_error += fabs(round(scaled_weight) - scaled_weight);
|
|
t->add_literals(var.value() + 1);
|
|
t->add_coefficients(static_cast<int64_t>(round(scaled_weight)) / gcd);
|
|
}
|
|
}
|
|
|
|
return bound_error;
|
|
}
|
|
bool IntegralProblemConverter::HasNonZeroWeights(
|
|
const absl::StrongVector<VariableIndex, Fractional>& dense_weights) const {
|
|
for (const Fractional weight : dense_weights) {
|
|
if (weight != 0.0) {
|
|
return true;
|
|
}
|
|
}
|
|
return false;
|
|
}
|
|
|
|
bool CheckSolution(const LinearProgram& linear_problem,
|
|
const glop::DenseRow& variable_values) {
|
|
glop::DenseColumn constraint_values(linear_problem.num_constraints(), 0);
|
|
|
|
const SparseMatrix& matrix = linear_problem.GetSparseMatrix();
|
|
for (ColIndex col(0); col < linear_problem.num_variables(); ++col) {
|
|
const Fractional lower_bound = linear_problem.variable_lower_bounds()[col];
|
|
const Fractional upper_bound = linear_problem.variable_upper_bounds()[col];
|
|
const Fractional value = variable_values[col];
|
|
if (lower_bound > value || upper_bound < value) {
|
|
LOG(ERROR) << "Variable " << col << " out of bound: " << value
|
|
<< " should be in " << lower_bound << " .. " << upper_bound;
|
|
return false;
|
|
}
|
|
|
|
for (const SparseColumn::Entry entry : matrix.column(col)) {
|
|
constraint_values[entry.row()] += entry.coefficient() * value;
|
|
}
|
|
}
|
|
|
|
for (RowIndex row(0); row < linear_problem.num_constraints(); ++row) {
|
|
const Fractional lb = linear_problem.constraint_lower_bounds()[row];
|
|
const Fractional ub = linear_problem.constraint_upper_bounds()[row];
|
|
const Fractional value = constraint_values[row];
|
|
if (lb > value || ub < value) {
|
|
LOG(ERROR) << "Constraint " << row << " out of bound: " << value
|
|
<< " should be in " << lb << " .. " << ub;
|
|
return false;
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
// Solves the given linear program and returns the solve status.
|
|
BopSolveStatus InternalSolve(const LinearProgram& linear_problem,
|
|
const BopParameters& parameters,
|
|
const DenseRow& initial_solution,
|
|
TimeLimit* time_limit, DenseRow* variable_values,
|
|
Fractional* objective_value,
|
|
Fractional* best_bound) {
|
|
CHECK(variable_values != nullptr);
|
|
CHECK(objective_value != nullptr);
|
|
CHECK(best_bound != nullptr);
|
|
const bool use_initial_solution = (initial_solution.size() > 0);
|
|
if (use_initial_solution) {
|
|
CHECK_EQ(initial_solution.size(), linear_problem.num_variables());
|
|
}
|
|
|
|
// Those values will only make sense when a solution is found, however we
|
|
// resize here such that one can access the values even if they don't mean
|
|
// anything.
|
|
variable_values->resize(linear_problem.num_variables(), 0);
|
|
|
|
LinearBooleanProblem boolean_problem;
|
|
std::vector<bool> boolean_initial_solution;
|
|
IntegralProblemConverter converter;
|
|
if (!converter.ConvertToBooleanProblem(linear_problem, initial_solution,
|
|
&boolean_problem,
|
|
&boolean_initial_solution)) {
|
|
return BopSolveStatus::INVALID_PROBLEM;
|
|
}
|
|
|
|
BopSolver bop_solver(boolean_problem);
|
|
bop_solver.SetParameters(parameters);
|
|
BopSolveStatus status = BopSolveStatus::NO_SOLUTION_FOUND;
|
|
if (use_initial_solution) {
|
|
BopSolution bop_solution(boolean_problem, "InitialSolution");
|
|
CHECK_EQ(boolean_initial_solution.size(), boolean_problem.num_variables());
|
|
for (int i = 0; i < boolean_initial_solution.size(); ++i) {
|
|
bop_solution.SetValue(VariableIndex(i), boolean_initial_solution[i]);
|
|
}
|
|
status = bop_solver.SolveWithTimeLimit(bop_solution, time_limit);
|
|
} else {
|
|
status = bop_solver.SolveWithTimeLimit(time_limit);
|
|
}
|
|
if (status == BopSolveStatus::OPTIMAL_SOLUTION_FOUND ||
|
|
status == BopSolveStatus::FEASIBLE_SOLUTION_FOUND) {
|
|
// Compute objective value.
|
|
const BopSolution& solution = bop_solver.best_solution();
|
|
CHECK(solution.IsFeasible());
|
|
|
|
*objective_value = linear_problem.objective_offset();
|
|
for (ColIndex col(0); col < linear_problem.num_variables(); ++col) {
|
|
const int64_t value = converter.GetSolutionValue(col, solution);
|
|
(*variable_values)[col] = value;
|
|
*objective_value += value * linear_problem.objective_coefficients()[col];
|
|
}
|
|
|
|
CheckSolution(linear_problem, *variable_values);
|
|
|
|
// TODO(user): Check that the scaled best bound from Bop is a valid one
|
|
// even after conversion. If yes, remove the optimality test.
|
|
*best_bound = status == BopSolveStatus::OPTIMAL_SOLUTION_FOUND
|
|
? *objective_value
|
|
: bop_solver.GetScaledBestBound();
|
|
}
|
|
return status;
|
|
}
|
|
|
|
void RunOneBop(const BopParameters& parameters, int problem_index,
|
|
const DenseRow& initial_solution, TimeLimit* time_limit,
|
|
LPDecomposer* decomposer, DenseRow* variable_values,
|
|
Fractional* objective_value, Fractional* best_bound,
|
|
BopSolveStatus* status) {
|
|
CHECK(decomposer != nullptr);
|
|
CHECK(variable_values != nullptr);
|
|
CHECK(objective_value != nullptr);
|
|
CHECK(best_bound != nullptr);
|
|
CHECK(status != nullptr);
|
|
|
|
LinearProgram problem;
|
|
decomposer->ExtractLocalProblem(problem_index, &problem);
|
|
DenseRow local_initial_solution;
|
|
if (initial_solution.size() > 0) {
|
|
local_initial_solution =
|
|
decomposer->ExtractLocalAssignment(problem_index, initial_solution);
|
|
}
|
|
// TODO(user): Investigate a better approximation of the time needed to
|
|
// solve the problem than just the number of variables.
|
|
const double total_num_variables = std::max(
|
|
1.0, static_cast<double>(
|
|
decomposer->original_problem().num_variables().value()));
|
|
const double time_per_variable =
|
|
parameters.max_time_in_seconds() / total_num_variables;
|
|
const double deterministic_time_per_variable =
|
|
parameters.max_deterministic_time() / total_num_variables;
|
|
const int local_num_variables = std::max(1, problem.num_variables().value());
|
|
|
|
NestedTimeLimit subproblem_time_limit(
|
|
time_limit,
|
|
std::max(time_per_variable * local_num_variables,
|
|
parameters.decomposed_problem_min_time_in_seconds()),
|
|
deterministic_time_per_variable * local_num_variables);
|
|
|
|
*status = InternalSolve(problem, parameters, local_initial_solution,
|
|
subproblem_time_limit.GetTimeLimit(), variable_values,
|
|
objective_value, best_bound);
|
|
}
|
|
} // anonymous namespace
|
|
|
|
IntegralSolver::IntegralSolver()
|
|
: parameters_(), variable_values_(), objective_value_(0.0) {}
|
|
|
|
BopSolveStatus IntegralSolver::Solve(const LinearProgram& linear_problem) {
|
|
return Solve(linear_problem, DenseRow());
|
|
}
|
|
|
|
BopSolveStatus IntegralSolver::SolveWithTimeLimit(
|
|
const LinearProgram& linear_problem, TimeLimit* time_limit) {
|
|
return SolveWithTimeLimit(linear_problem, DenseRow(), time_limit);
|
|
}
|
|
|
|
BopSolveStatus IntegralSolver::Solve(
|
|
const LinearProgram& linear_problem,
|
|
const DenseRow& user_provided_initial_solution) {
|
|
std::unique_ptr<TimeLimit> time_limit =
|
|
TimeLimit::FromParameters(parameters_);
|
|
return SolveWithTimeLimit(linear_problem, user_provided_initial_solution,
|
|
time_limit.get());
|
|
}
|
|
|
|
BopSolveStatus IntegralSolver::SolveWithTimeLimit(
|
|
const LinearProgram& linear_problem,
|
|
const DenseRow& user_provided_initial_solution, TimeLimit* time_limit) {
|
|
// We make a copy so that we can clear it if the presolve is active.
|
|
DenseRow initial_solution = user_provided_initial_solution;
|
|
if (initial_solution.size() > 0) {
|
|
CHECK_EQ(initial_solution.size(), linear_problem.num_variables())
|
|
<< "The initial solution should have the same number of variables as "
|
|
"the LinearProgram.";
|
|
}
|
|
|
|
// Some code path requires to copy the given linear_problem. When this
|
|
// happens, we will simply change the target of this pointer.
|
|
LinearProgram const* lp = &linear_problem;
|
|
|
|
BopSolveStatus status;
|
|
if (lp->num_variables() >= parameters_.decomposer_num_variables_threshold()) {
|
|
LPDecomposer decomposer;
|
|
decomposer.Decompose(lp);
|
|
const int num_sub_problems = decomposer.GetNumberOfProblems();
|
|
VLOG(1) << "Problem is decomposable into " << num_sub_problems
|
|
<< " components!";
|
|
if (num_sub_problems > 1) {
|
|
// The problem can be decomposed. Solve each sub-problem and aggregate the
|
|
// result.
|
|
std::vector<DenseRow> variable_values(num_sub_problems);
|
|
std::vector<Fractional> objective_values(num_sub_problems,
|
|
Fractional(0.0));
|
|
std::vector<Fractional> best_bounds(num_sub_problems, Fractional(0.0));
|
|
std::vector<BopSolveStatus> statuses(num_sub_problems,
|
|
BopSolveStatus::INVALID_PROBLEM);
|
|
|
|
for (int i = 0; i < num_sub_problems; ++i) {
|
|
RunOneBop(parameters_, i, initial_solution, time_limit, &decomposer,
|
|
&(variable_values[i]), &(objective_values[i]),
|
|
&(best_bounds[i]), &(statuses[i]));
|
|
}
|
|
|
|
// Aggregate results.
|
|
status = BopSolveStatus::OPTIMAL_SOLUTION_FOUND;
|
|
objective_value_ = lp->objective_offset();
|
|
best_bound_ = 0.0;
|
|
for (int i = 0; i < num_sub_problems; ++i) {
|
|
objective_value_ += objective_values[i];
|
|
best_bound_ += best_bounds[i];
|
|
if (statuses[i] == BopSolveStatus::NO_SOLUTION_FOUND ||
|
|
statuses[i] == BopSolveStatus::INFEASIBLE_PROBLEM ||
|
|
statuses[i] == BopSolveStatus::INVALID_PROBLEM) {
|
|
return statuses[i];
|
|
}
|
|
|
|
if (statuses[i] == BopSolveStatus::FEASIBLE_SOLUTION_FOUND) {
|
|
status = BopSolveStatus::FEASIBLE_SOLUTION_FOUND;
|
|
}
|
|
}
|
|
variable_values_ = decomposer.AggregateAssignments(variable_values);
|
|
CheckSolution(*lp, variable_values_);
|
|
} else {
|
|
status =
|
|
InternalSolve(*lp, parameters_, initial_solution, time_limit,
|
|
&variable_values_, &objective_value_, &best_bound_);
|
|
}
|
|
} else {
|
|
status = InternalSolve(*lp, parameters_, initial_solution, time_limit,
|
|
&variable_values_, &objective_value_, &best_bound_);
|
|
}
|
|
|
|
return status;
|
|
}
|
|
|
|
} // namespace bop
|
|
} // namespace operations_research
|