always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
146566e2be8f5fdd4725cb36a70719d92ebf32c7
ortools-clone/ortools/sat/feasibility_pump.cc

739 lines
28 KiB
C++
Raw Blame History

// Copyright 2010-2025 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/sat/feasibility_pump.h"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <limits>
#include <utility>
#include <vector>
#include "absl/container/flat_hash_map.h"
#include "absl/log/check.h"
#include "absl/meta/type_traits.h"
#include "ortools/base/logging.h"
#include "ortools/base/strong_vector.h"
#include "ortools/glop/parameters.pb.h"
#include "ortools/glop/revised_simplex.h"
#include "ortools/glop/status.h"
#include "ortools/lp_data/lp_data.h"
#include "ortools/lp_data/lp_data_utils.h"
#include "ortools/lp_data/lp_types.h"
#include "ortools/lp_data/sparse_column.h"
#include "ortools/sat/cp_model_mapping.h"
#include "ortools/sat/integer.h"
#include "ortools/sat/integer_base.h"
#include "ortools/sat/linear_constraint.h"
#include "ortools/sat/model.h"
#include "ortools/sat/sat_base.h"
#include "ortools/sat/sat_parameters.pb.h"
#include "ortools/sat/sat_solver.h"
#include "ortools/sat/synchronization.h"
#include "ortools/sat/util.h"
#include "ortools/util/saturated_arithmetic.h"
#include "ortools/util/sorted_interval_list.h"
#include "ortools/util/strong_integers.h"
#include "ortools/util/time_limit.h"
namespace operations_research {
namespace sat {
using glop::ColIndex;
using glop::ConstraintStatus;
using glop::Fractional;
using glop::RowIndex;
const double FeasibilityPump::kCpEpsilon = 1e-4;
FeasibilityPump::FeasibilityPump(Model* model)
: sat_parameters_(*(model->GetOrCreate<SatParameters>())),
time_limit_(model->GetOrCreate<TimeLimit>()),
integer_trail_(model->GetOrCreate<IntegerTrail>()),
trail_(model->GetOrCreate<Trail>()),
integer_encoder_(model->GetOrCreate<IntegerEncoder>()),
incomplete_solutions_(model->Mutable<SharedIncompleteSolutionManager>()),
sat_solver_(model->GetOrCreate<SatSolver>()),
domains_(model->GetOrCreate<IntegerDomains>()),
mapping_(model->Get<CpModelMapping>()) {
// Tweak the default parameters to make the solve incremental.
glop::GlopParameters parameters;
// Note(user): Primal simplex does better here since we have a limit on
// simplex iterations. So dual simplex sometimes fails to find a LP feasible
// solution.
parameters.set_use_dual_simplex(false);
parameters.set_max_number_of_iterations(2000);
simplex_.SetParameters(parameters);
lp_data_.Clear();
integer_lp_.clear();
}
FeasibilityPump::~FeasibilityPump() {
VLOG(1) << "Feasibility Pump Total number of simplex iterations: "
<< total_num_simplex_iterations_;
}
void FeasibilityPump::AddLinearConstraint(const LinearConstraint& ct) {
// We still create the mirror variable right away though.
for (const IntegerVariable var : ct.VarsAsSpan()) {
GetOrCreateMirrorVariable(PositiveVariable(var));
}
integer_lp_.push_back(LinearConstraintInternal());
LinearConstraintInternal& new_ct = integer_lp_.back();
new_ct.lb = ct.lb;
new_ct.ub = ct.ub;
CHECK_LE(ct.lb, ct.ub);
for (int i = 0; i < ct.num_terms; ++i) {
// We only use positive variable inside this class.
IntegerVariable var = ct.vars[i];
IntegerValue coeff = ct.coeffs[i];
if (!VariableIsPositive(var)) {
var = NegationOf(var);
coeff = -coeff;
}
new_ct.terms.push_back({GetOrCreateMirrorVariable(var), coeff});
}
// Important to keep lp_data_ "clean".
std::sort(new_ct.terms.begin(), new_ct.terms.end());
}
void FeasibilityPump::SetObjectiveCoefficient(IntegerVariable ivar,
IntegerValue coeff) {
objective_is_defined_ = true;
const IntegerVariable pos_var =
VariableIsPositive(ivar) ? ivar : NegationOf(ivar);
if (ivar != pos_var) coeff = -coeff;
const auto it = mirror_lp_variable_.find(pos_var);
if (it == mirror_lp_variable_.end()) return;
const ColIndex col = it->second;
integer_objective_.push_back({col, coeff});
objective_infinity_norm_ =
std::max(objective_infinity_norm_, IntTypeAbs(coeff));
}
ColIndex FeasibilityPump::GetOrCreateMirrorVariable(
IntegerVariable positive_variable) {
DCHECK(VariableIsPositive(positive_variable));
const auto it = mirror_lp_variable_.find(positive_variable);
if (it == mirror_lp_variable_.end()) {
const int model_var =
mapping_->GetProtoVariableFromIntegerVariable(positive_variable);
model_vars_size_ = std::max(model_vars_size_, model_var + 1);
const ColIndex col(integer_variables_.size());
mirror_lp_variable_[positive_variable] = col;
integer_variables_.push_back(positive_variable);
var_is_binary_.push_back(false);
lp_solution_.push_back(std::numeric_limits<double>::infinity());
integer_solution_.push_back(0);
return col;
}
return it->second;
}
void FeasibilityPump::PrintStats() {
if (lp_solution_is_set_) {
VLOG(2) << "Fractionality: " << lp_solution_fractionality_;
} else {
VLOG(2) << "Fractionality: NA";
VLOG(2) << "simplex status: " << simplex_.GetProblemStatus();
}
if (integer_solution_is_set_) {
VLOG(2) << "#Infeasible const: " << num_infeasible_constraints_;
VLOG(2) << "Infeasibility: " << integer_solution_infeasibility_;
} else {
VLOG(2) << "Infeasibility: NA";
}
}
bool FeasibilityPump::Solve() {
if (lp_data_.num_variables() == 0) {
InitializeWorkingLP();
}
UpdateBoundsOfLpVariables();
lp_solution_is_set_ = false;
integer_solution_is_set_ = false;
// Restore the original objective
for (ColIndex col(0); col < lp_data_.num_variables(); ++col) {
lp_data_.SetObjectiveCoefficient(col, 0.0);
}
for (const auto& term : integer_objective_) {
lp_data_.SetObjectiveCoefficient(term.first, ToDouble(term.second));
}
mixing_factor_ = 1.0;
for (int i = 0; i < max_fp_iterations_; ++i) {
if (time_limit_->LimitReached()) break;
L1DistanceMinimize();
if (!SolveLp()) break;
if (lp_solution_is_integer_) break;
if (!Round()) break;
// We don't end this loop if the integer solutions is feasible in hope to
// get better solution.
if (integer_solution_is_feasible_) MaybePushToRepo();
}
if (sat_solver_->ModelIsUnsat()) return false;
PrintStats();
MaybePushToRepo();
return true;
}
void FeasibilityPump::MaybePushToRepo() {
if (incomplete_solutions_ == nullptr) return;
const double kInfinity = std::numeric_limits<double>::infinity();
// TODO(user): Consider adding solutions that have low fractionality.
if (lp_solution_is_integer_) {
// Fill the solution using LP solution values.
tmp_solution_.assign(model_vars_size_, kInfinity);
for (const IntegerVariable positive_var : integer_variables_) {
const int model_var =
mapping_->GetProtoVariableFromIntegerVariable(positive_var);
if (model_var >= 0 && model_var < model_vars_size_) {
tmp_solution_[model_var] = GetLPSolutionValue(positive_var);
}
}
incomplete_solutions_->AddSolution(tmp_solution_);
}
if (integer_solution_is_feasible_) {
// Fill the solution using Integer solution values.
tmp_solution_.assign(model_vars_size_, kInfinity);
for (const IntegerVariable positive_var : integer_variables_) {
const int model_var =
mapping_->GetProtoVariableFromIntegerVariable(positive_var);
if (model_var >= 0 && model_var < model_vars_size_) {
tmp_solution_[model_var] = GetIntegerSolutionValue(positive_var);
}
}
incomplete_solutions_->AddSolution(tmp_solution_);
}
}
// ----------------------------------------------------------------
// -------------------LPSolving------------------------------------
// ----------------------------------------------------------------
void FeasibilityPump::InitializeWorkingLP() {
lp_data_.Clear();
// Create variables.
for (int i = 0; i < integer_variables_.size(); ++i) {
CHECK_EQ(ColIndex(i), lp_data_.CreateNewVariable());
lp_data_.SetVariableType(ColIndex(i),
glop::LinearProgram::VariableType::INTEGER);
}
// Add constraints.
for (const LinearConstraintInternal& ct : integer_lp_) {
const ConstraintIndex row = lp_data_.CreateNewConstraint();
lp_data_.SetConstraintBounds(row, ToDouble(ct.lb), ToDouble(ct.ub));
for (const auto& term : ct.terms) {
lp_data_.SetCoefficient(row, term.first, ToDouble(term.second));
}
}
// Add objective.
for (const auto& term : integer_objective_) {
lp_data_.SetObjectiveCoefficient(term.first, ToDouble(term.second));
}
const int num_vars = integer_variables_.size();
for (int i = 0; i < num_vars; i++) {
const IntegerVariable cp_var = integer_variables_[i];
const double lb = ToDouble(integer_trail_->LevelZeroLowerBound(cp_var));
const double ub = ToDouble(integer_trail_->LevelZeroUpperBound(cp_var));
lp_data_.SetVariableBounds(ColIndex(i), lb, ub);
}
objective_normalization_factor_ = 0.0;
glop::ColIndexVector integer_variables;
const ColIndex num_cols = lp_data_.num_variables();
for (ColIndex col : lp_data_.IntegerVariablesList()) {
var_is_binary_[col.value()] = lp_data_.IsVariableBinary(col);
if (!var_is_binary_[col.value()]) {
integer_variables.push_back(col);
}
// The aim of this normalization value is to compute a coefficient of the
// d_i variables that should be minimized.
objective_normalization_factor_ +=
std::abs(lp_data_.GetObjectiveCoefficientForMinimizationVersion(col));
}
CHECK_GT(lp_data_.IntegerVariablesList().size(), 0);
objective_normalization_factor_ =
objective_normalization_factor_ / lp_data_.IntegerVariablesList().size();
if (!integer_variables.empty()) {
// Update the LpProblem with norm variables and constraints.
norm_variables_.assign(num_cols, ColIndex(-1));
norm_lhs_constraints_.assign(num_cols, RowIndex(-1));
norm_rhs_constraints_.assign(num_cols, RowIndex(-1));
// For each integer non-binary variable x_i we introduce one new variable
// d_i subject to two new constraints:
// d_i - x_i >= -round(x'_i)
// d_i + x_i >= +round(x'_i)
// That's round(x'_i) - d_i <= x_i <= round(x'_i) + d_i, where d_i is an
// unbounded non-negative, and x'_i is the value of variable i from the
// previous solution obtained during the projection step. Consequently
// coefficients of the constraints are set here, but bounds of the
// constraints are updated at each iteration of the feasibility pump. Also
// coefficients of the objective are set here: x_i's are not present in the
// objective (i.e., coefficients set to 0.0), and d_i's are present in the
// objective with coefficients set to 1.0.
// Note that the treatment of integer non-binary variables is different
// from the treatment of binary variables. Binary variables do not impose
// any extra variables, nor extra constraints, but their objective
// coefficients are changed in the linear projection steps.
for (const ColIndex col : integer_variables) {
const ColIndex norm_variable = lp_data_.CreateNewVariable();
norm_variables_[col] = norm_variable;
lp_data_.SetVariableBounds(norm_variable, 0.0, glop::kInfinity);
const RowIndex row_a = lp_data_.CreateNewConstraint();
norm_lhs_constraints_[col] = row_a;
lp_data_.SetCoefficient(row_a, norm_variable, 1.0);
lp_data_.SetCoefficient(row_a, col, -1.0);
const RowIndex row_b = lp_data_.CreateNewConstraint();
norm_rhs_constraints_[col] = row_b;
lp_data_.SetCoefficient(row_b, norm_variable, 1.0);
lp_data_.SetCoefficient(row_b, col, 1.0);
}
}
scaler_.Scale(&lp_data_);
lp_data_.AddSlackVariablesWhereNecessary(
/*detect_integer_constraints=*/false);
}
void FeasibilityPump::L1DistanceMinimize() {
std::vector<double> new_obj_coeffs(lp_data_.num_variables().value(), 0.0);
// Set the original subobjective. The coefficients are scaled by mixing factor
// and the offset remains at 0 (because it does not affect the solution).
const ColIndex num_cols(lp_data_.objective_coefficients().size());
for (ColIndex col(0); col < num_cols; ++col) {
new_obj_coeffs[col.value()] =
mixing_factor_ * lp_data_.objective_coefficients()[col];
}
// Set the norm subobjective. The coefficients are scaled by 1 - mixing factor
// and the offset remains at 0 (because it does not affect the solution).
for (const ColIndex col : lp_data_.IntegerVariablesList()) {
if (var_is_binary_[col.value()]) {
const Fractional objective_coefficient =
mixing_factor_ * lp_data_.objective_coefficients()[col] +
(1 - mixing_factor_) * objective_normalization_factor_ *
(1 - 2 * integer_solution_[col.value()]);
new_obj_coeffs[col.value()] = objective_coefficient;
} else { // The variable is integer.
// Update the bounds of the constraints added in
// InitializeIntegerVariables() (see there for more details):
// d_i - x_i >= -round(x'_i)
// d_i + x_i >= +round(x'_i)
// TODO(user): We change both the objective and the bounds, thus
// breaking the incrementality. Handle integer variables differently,
// e.g., intensify rounding, or use soft fixing from: Fischetti, Lodi,
// "Local Branching", Math Program Ser B 98:23-47 (2003).
const Fractional objective_coefficient =
(1 - mixing_factor_) * objective_normalization_factor_;
new_obj_coeffs[norm_variables_[col].value()] = objective_coefficient;
// At this point, constraint bounds have already been transformed into
// bounds of slack variables. Instead of updating the constraints, we need
// to update the slack variables corresponding to them.
const ColIndex norm_lhs_slack_variable =
lp_data_.GetSlackVariable(norm_lhs_constraints_[col]);
const double lhs_scaling_factor =
scaler_.VariableScalingFactorWithSlack(norm_lhs_slack_variable);
lp_data_.SetVariableBounds(
norm_lhs_slack_variable, -glop::kInfinity,
lhs_scaling_factor * integer_solution_[col.value()]);
const ColIndex norm_rhs_slack_variable =
lp_data_.GetSlackVariable(norm_rhs_constraints_[col]);
const double rhs_scaling_factor =
scaler_.VariableScalingFactorWithSlack(norm_rhs_slack_variable);
lp_data_.SetVariableBounds(
norm_rhs_slack_variable, -glop::kInfinity,
-rhs_scaling_factor * integer_solution_[col.value()]);
}
}
for (ColIndex col(0); col < lp_data_.num_variables(); ++col) {
lp_data_.SetObjectiveCoefficient(col, new_obj_coeffs[col.value()]);
}
// TODO(user): Tune this or expose as parameter.
mixing_factor_ *= 0.8;
}
bool FeasibilityPump::SolveLp() {
const int num_vars = integer_variables_.size();
VLOG(3) << "LP relaxation: " << lp_data_.GetDimensionString() << ".";
const auto status = simplex_.Solve(lp_data_, time_limit_);
total_num_simplex_iterations_ += simplex_.GetNumberOfIterations();
if (!status.ok()) {
VLOG(1) << "The LP solver encountered an error: " << status.error_message();
simplex_.ClearStateForNextSolve();
return false;
}
// TODO(user): This shouldn't really happen except if the problem is UNSAT.
// But we can't just rely on a potentially imprecise LP to close the problem.
// The rest of the solver should do that with exact precision.
VLOG(3) << "simplex status: " << simplex_.GetProblemStatus();
if (simplex_.GetProblemStatus() == glop::ProblemStatus::PRIMAL_INFEASIBLE) {
return false;
}
lp_solution_fractionality_ = 0.0;
if (simplex_.GetProblemStatus() == glop::ProblemStatus::OPTIMAL ||
simplex_.GetProblemStatus() == glop::ProblemStatus::DUAL_FEASIBLE ||
simplex_.GetProblemStatus() == glop::ProblemStatus::PRIMAL_FEASIBLE ||
simplex_.GetProblemStatus() == glop::ProblemStatus::IMPRECISE) {
lp_solution_is_set_ = true;
for (int i = 0; i < num_vars; i++) {
const double value = GetVariableValueAtCpScale(ColIndex(i));
lp_solution_[i] = value;
lp_solution_fractionality_ = std::max(
lp_solution_fractionality_, std::abs(value - std::round(value)));
}
// Compute the objective value.
lp_objective_ = 0;
for (const auto& term : integer_objective_) {
lp_objective_ += lp_solution_[term.first.value()] * term.second.value();
}
lp_solution_is_integer_ = lp_solution_fractionality_ < kCpEpsilon;
}
return true;
}
void FeasibilityPump::UpdateBoundsOfLpVariables() {
const int num_vars = integer_variables_.size();
for (int i = 0; i < num_vars; i++) {
const IntegerVariable cp_var = integer_variables_[i];
const double lb = ToDouble(integer_trail_->LevelZeroLowerBound(cp_var));
const double ub = ToDouble(integer_trail_->LevelZeroUpperBound(cp_var));
const double factor = scaler_.VariableScalingFactor(ColIndex(i));
lp_data_.SetVariableBounds(ColIndex(i), lb * factor, ub * factor);
}
}
double FeasibilityPump::GetLPSolutionValue(IntegerVariable variable) const {
return lp_solution_[mirror_lp_variable_.at(variable).value()];
}
double FeasibilityPump::GetVariableValueAtCpScale(ColIndex var) {
return scaler_.UnscaleVariableValue(var, simplex_.GetVariableValue(var));
}
// ----------------------------------------------------------------
// -------------------Rounding-------------------------------------
// ----------------------------------------------------------------
int64_t FeasibilityPump::GetIntegerSolutionValue(
IntegerVariable variable) const {
return integer_solution_[mirror_lp_variable_.at(variable).value()];
}
bool FeasibilityPump::Round() {
bool rounding_successful = true;
if (sat_parameters_.fp_rounding() == SatParameters::NEAREST_INTEGER) {
rounding_successful = NearestIntegerRounding();
} else if (sat_parameters_.fp_rounding() == SatParameters::LOCK_BASED) {
rounding_successful = LockBasedRounding();
} else if (sat_parameters_.fp_rounding() ==
SatParameters::ACTIVE_LOCK_BASED) {
rounding_successful = ActiveLockBasedRounding();
} else if (sat_parameters_.fp_rounding() ==
SatParameters::PROPAGATION_ASSISTED) {
rounding_successful = PropagationRounding();
}
if (!rounding_successful) return false;
FillIntegerSolutionStats();
return true;
}
bool FeasibilityPump::NearestIntegerRounding() {
if (!lp_solution_is_set_) return false;
for (int i = 0; i < lp_solution_.size(); ++i) {
integer_solution_[i] = static_cast<int64_t>(std::round(lp_solution_[i]));
}
integer_solution_is_set_ = true;
return true;
}
bool FeasibilityPump::LockBasedRounding() {
if (!lp_solution_is_set_) return false;
const int num_vars = integer_variables_.size();
// We compute the number of locks based on variable coefficient in constraints
// and constraint bounds. This doesn't change over time so we cache it.
if (var_up_locks_.empty()) {
var_up_locks_.resize(num_vars, 0);
var_down_locks_.resize(num_vars, 0);
for (int i = 0; i < num_vars; ++i) {
for (const auto entry : lp_data_.GetSparseColumn(ColIndex(i))) {
ColIndex slack = lp_data_.GetSlackVariable(entry.row());
const bool constraint_upper_bounded =
lp_data_.variable_lower_bounds()[slack] > -glop::kInfinity;
const bool constraint_lower_bounded =
lp_data_.variable_upper_bounds()[slack] < glop::kInfinity;
if (entry.coefficient() > 0) {
var_up_locks_[i] += constraint_upper_bounded;
var_down_locks_[i] += constraint_lower_bounded;
} else {
var_up_locks_[i] += constraint_lower_bounded;
var_down_locks_[i] += constraint_upper_bounded;
}
}
}
}
for (int i = 0; i < lp_solution_.size(); ++i) {
if (std::abs(lp_solution_[i] - std::round(lp_solution_[i])) < 0.1 ||
var_up_locks_[i] == var_down_locks_[i]) {
integer_solution_[i] = static_cast<int64_t>(std::round(lp_solution_[i]));
} else if (var_up_locks_[i] > var_down_locks_[i]) {
integer_solution_[i] = static_cast<int64_t>(std::floor(lp_solution_[i]));
} else {
integer_solution_[i] = static_cast<int64_t>(std::ceil(lp_solution_[i]));
}
}
integer_solution_is_set_ = true;
return true;
}
bool FeasibilityPump::ActiveLockBasedRounding() {
if (!lp_solution_is_set_) return false;
const int num_vars = integer_variables_.size();
// We compute the number of locks based on variable coefficient in constraints
// and constraint bounds of active constraints. We consider the bound of the
// constraint that is tight for the current lp solution.
for (int i = 0; i < num_vars; ++i) {
if (std::abs(lp_solution_[i] - std::round(lp_solution_[i])) < 0.1) {
integer_solution_[i] = static_cast<int64_t>(std::round(lp_solution_[i]));
}
int up_locks = 0;
int down_locks = 0;
for (const auto entry : lp_data_.GetSparseColumn(ColIndex(i))) {
const ConstraintStatus row_status =
simplex_.GetConstraintStatus(entry.row());
if (row_status == ConstraintStatus::AT_LOWER_BOUND) {
if (entry.coefficient() > 0) {
down_locks++;
} else {
up_locks++;
}
} else if (row_status == ConstraintStatus::AT_UPPER_BOUND) {
if (entry.coefficient() > 0) {
up_locks++;
} else {
down_locks++;
}
}
}
if (up_locks == down_locks) {
integer_solution_[i] = static_cast<int64_t>(std::round(lp_solution_[i]));
} else if (up_locks > down_locks) {
integer_solution_[i] = static_cast<int64_t>(std::floor(lp_solution_[i]));
} else {
integer_solution_[i] = static_cast<int64_t>(std::ceil(lp_solution_[i]));
}
}
integer_solution_is_set_ = true;
return true;
}
bool FeasibilityPump::PropagationRounding() {
if (!lp_solution_is_set_) return false;
if (!sat_solver_->ResetToLevelZero()) return false;
// Compute an order in which we will fix variables and do the propagation.
std::vector<int> rounding_order;
{
std::vector<std::pair<double, int>> binary_fractionality_vars;
std::vector<std::pair<double, int>> general_fractionality_vars;
for (int i = 0; i < lp_solution_.size(); ++i) {
const double fractionality =
std::abs(std::round(lp_solution_[i]) - lp_solution_[i]);
if (var_is_binary_[i]) {
binary_fractionality_vars.push_back({fractionality, i});
} else {
general_fractionality_vars.push_back({fractionality, i});
}
}
std::sort(binary_fractionality_vars.begin(),
binary_fractionality_vars.end());
std::sort(general_fractionality_vars.begin(),
general_fractionality_vars.end());
for (int i = 0; i < binary_fractionality_vars.size(); ++i) {
rounding_order.push_back(binary_fractionality_vars[i].second);
}
for (int i = 0; i < general_fractionality_vars.size(); ++i) {
rounding_order.push_back(general_fractionality_vars[i].second);
}
}
for (const int var_index : rounding_order) {
if (time_limit_->LimitReached()) return false;
// Get the bounds of the variable.
const IntegerVariable var = integer_variables_[var_index];
CHECK(VariableIsPositive(var));
const Domain& domain = (*domains_)[GetPositiveOnlyIndex(var)];
const IntegerValue lb = integer_trail_->LowerBound(var);
const IntegerValue ub = integer_trail_->UpperBound(var);
if (lb == ub) {
integer_solution_[var_index] = lb.value();
continue;
}
const int64_t rounded_value =
SafeDoubleToInt64(std::round(lp_solution_[var_index]));
const int64_t floor_value =
SafeDoubleToInt64(std::floor(lp_solution_[var_index]));
const int64_t ceil_value =
SafeDoubleToInt64(std::ceil(lp_solution_[var_index]));
const bool floor_is_in_domain =
(domain.Contains(floor_value) && lb.value() <= floor_value);
const bool ceil_is_in_domain =
(domain.Contains(ceil_value) && ub.value() >= ceil_value);
if (domain.IsEmpty()) {
integer_solution_[var_index] = rounded_value;
sat_solver_->NotifyThatModelIsUnsat();
return false;
}
if (ceil_value < lb.value()) {
integer_solution_[var_index] = lb.value();
} else if (floor_value > ub.value()) {
integer_solution_[var_index] = ub.value();
} else if (ceil_is_in_domain && floor_is_in_domain) {
DCHECK(domain.Contains(rounded_value));
integer_solution_[var_index] = rounded_value;
} else if (ceil_is_in_domain) {
integer_solution_[var_index] = ceil_value;
} else if (floor_is_in_domain) {
integer_solution_[var_index] = floor_value;
} else {
const std::pair<IntegerLiteral, IntegerLiteral> values_in_domain =
integer_encoder_->Canonicalize(
IntegerLiteral::GreaterOrEqual(var, IntegerValue(rounded_value)));
const int64_t lower_value = values_in_domain.first.bound.value();
const int64_t higher_value = -values_in_domain.second.bound.value();
const int64_t distance_from_lower_value =
std::abs(lower_value - rounded_value);
const int64_t distance_from_higher_value =
std::abs(higher_value - rounded_value);
integer_solution_[var_index] =
(distance_from_lower_value < distance_from_higher_value)
? lower_value
: higher_value;
}
CHECK(domain.Contains(integer_solution_[var_index]));
CHECK_GE(integer_solution_[var_index], lb);
CHECK_LE(integer_solution_[var_index], ub);
// Propagate the value.
//
// When we want to fix the variable at its lb or ub, we do not create an
// equality literal to minimize the number of new literal we create. This
// is because creating an "== value" literal will implicitly also create
// a ">= value" and a "<= value" literals.
Literal to_enqueue;
const IntegerValue value(integer_solution_[var_index]);
if (value == lb) {
to_enqueue = integer_encoder_->GetOrCreateAssociatedLiteral(
IntegerLiteral::LowerOrEqual(var, value));
} else if (value == ub) {
to_enqueue = integer_encoder_->GetOrCreateAssociatedLiteral(
IntegerLiteral::GreaterOrEqual(var, value));
} else {
to_enqueue =
integer_encoder_->GetOrCreateLiteralAssociatedToEquality(var, value);
}
if (!sat_solver_->FinishPropagation()) return false;
const SatSolver::Status decision_status =
sat_solver_->EnqueueDecisionAndBacktrackOnConflict(to_enqueue);
if (decision_status != SatSolver::Status::FEASIBLE) return false;
if (sat_solver_->ModelIsUnsat()) return false;
}
integer_solution_is_set_ = true;
return sat_solver_->ResetToLevelZero();
}
void FeasibilityPump::FillIntegerSolutionStats() {
// Compute the objective value.
integer_solution_objective_ = 0;
for (const auto& term : integer_objective_) {
integer_solution_objective_ +=
integer_solution_[term.first.value()] * term.second.value();
}
integer_solution_is_feasible_ = true;
num_infeasible_constraints_ = 0;
integer_solution_infeasibility_ = 0;
for (RowIndex i(0); i < integer_lp_.size(); ++i) {
int64_t activity = 0;
for (const auto& term : integer_lp_[i].terms) {
const int64_t prod =
CapProd(integer_solution_[term.first.value()], term.second.value());
if (prod <= std::numeric_limits<int64_t>::min() ||
prod >= std::numeric_limits<int64_t>::max()) {
activity = prod;
break;
}
activity = CapAdd(activity, prod);
if (activity <= std::numeric_limits<int64_t>::min() ||
activity >= std::numeric_limits<int64_t>::max())
break;
}
if (activity > integer_lp_[i].ub || activity < integer_lp_[i].lb) {
integer_solution_is_feasible_ = false;
num_infeasible_constraints_++;
const int64_t ub_infeasibility =
activity > integer_lp_[i].ub.value()
? activity - integer_lp_[i].ub.value()
: 0;
const int64_t lb_infeasibility =
activity < integer_lp_[i].lb.value()
? integer_lp_[i].lb.value() - activity
: 0;
integer_solution_infeasibility_ =
std::max(integer_solution_infeasibility_,
std::max(ub_infeasibility, lb_infeasibility));
}
}
}
} // namespace sat
} // namespace operations_research
Reference in New Issue View Git Blame Copy Permalink