always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

split lp_data out of glop; port rest of code

Browse Source
This commit is contained in:
lperron@google.com
2014-12-16 11:35:09 +00:00
parent dd997176d6
commit aeb6434297
70 changed files with 875 additions and 527 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
examples/cpp/mps_driver.cc
View File

@@ -27,11 +27,11 @@
#include "google/protobuf/message.h"
#include "google/protobuf/text_format.h"
#include "base/strutil.h"
#include "glop/lp_print_utils.h"
#include "glop/lp_solver.h"
#include "glop/mps_reader.h"
#include "glop/parameters.pb.h"
#include "glop/proto_utils.h"
#include "lp_data/lp_print_utils.h"
#include "lp_data/mps_reader.h"
#include "linear_solver/proto_tools.h"
#include "base/status.h"

2
examples/cpp/opb_reader.h
View File

@@ -39,7 +39,6 @@ class OpbReader {
bool Load(const std::string& filename, LinearBooleanProblem* problem) {
problem->Clear();
problem->set_name(ExtractProblemName(filename));
problem->set_type(LinearBooleanProblem::SATISFIABILITY);
num_variables_ = 0;
int num_lines = 0;
@@ -72,7 +71,6 @@ class OpbReader {
}
if (words[0] == "min:") {
problem->set_type(LinearBooleanProblem::MINIMIZATION);
LinearObjective* objective = problem->mutable_objective();
for (int i = 1; i < words.size(); ++i) {
const std::string& word = words[i];

5
examples/cpp/sat_cnf_reader.h
View File

@@ -127,11 +127,6 @@ class SatCnfReader {
if (words_[1] == "wcnf") {
is_wcnf_ = true;
hard_weight_ = (words_.size() > 4) ? StringPieceAtoi(words_[4]) : 0;
problem->set_type(LinearBooleanProblem::MINIMIZATION);
} else {
problem->set_type(interpret_cnf_as_max_sat_
? LinearBooleanProblem::MINIMIZATION
: LinearBooleanProblem::SATISFIABILITY);
}
} else {
LOG(FATAL) << "Unknow file type: " << words_[1];

2
examples/cpp/sat_runner.cc
View File

@@ -223,7 +223,7 @@ int Run() {
problem, !FLAGS_lower_bound.empty(),
Coefficient(atoi64(FLAGS_lower_bound)), !FLAGS_upper_bound.empty(),
Coefficient(atoi64(FLAGS_upper_bound)), solver.get())) {
LOG(FATAL) << "Issue when setting the objective bounds.";
LOG(INFO) << "UNSAT when setting the objective constraint.";
}
// Symmetries!

96
makefiles/Makefile.cpp.mk
View File

@@ -886,29 +886,66 @@ $(LIB_DIR)/$(LIBPREFIX)base.$(STATIC_LIB_SUFFIX): $(BASE_LIB_OBJS)
$(STATIC_LINK_CMD) $(STATIC_LINK_PREFIX)$(LIB_DIR)$S$(LIBPREFIX)base.$(STATIC_LIB_SUFFIX) $(BASE_LIB_OBJS)
endif
GLOP_LIB_OBJS= \
$(OBJ_DIR)/glop/parameters.pb.$O \
$(OBJ_DIR)/glop/entering_variable.$O \
# Glop library.
LP_DATA_OBJS= \
$(OBJ_DIR)/lp_data/lp_data.$O \
$(OBJ_DIR)/lp_data/lp_print_utils.$O \
$(OBJ_DIR)/lp_data/lp_types.$O \
$(OBJ_DIR)/lp_data/lp_utils.$O \
$(OBJ_DIR)/lp_data/matrix_scaler.$O \
$(OBJ_DIR)/lp_data/matrix_utils.$O \
$(OBJ_DIR)/lp_data/mps_reader.$O \
$(OBJ_DIR)/lp_data/sparse.$O \
$(OBJ_DIR)/lp_data/sparse_column.$O \
$(OBJ_DIR)/lp_data/lp_data.$O:$(SRC_DIR)/lp_data/lp_data.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Slp_data$Slp_data.cc $(OBJ_OUT)$(OBJ_DIR)$Slp_data$Slp_data.$O
$(OBJ_DIR)/lp_data/lp_print_utils.$O:$(SRC_DIR)/lp_data/lp_print_utils.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Slp_data$Slp_print_utils.cc $(OBJ_OUT)$(OBJ_DIR)$Slp_data$Slp_print_utils.$O
$(OBJ_DIR)/lp_data/lp_types.$O:$(SRC_DIR)/lp_data/lp_types.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Slp_data$Slp_types.cc $(OBJ_OUT)$(OBJ_DIR)$Slp_data$Slp_types.$O
$(OBJ_DIR)/lp_data/lp_utils.$O:$(SRC_DIR)/lp_data/lp_utils.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Slp_data$Slp_utils.cc $(OBJ_OUT)$(OBJ_DIR)$Slp_data$Slp_utils.$O
$(OBJ_DIR)/lp_data/matrix_scaler.$O:$(SRC_DIR)/lp_data/matrix_scaler.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Slp_data$Smatrix_scaler.cc $(OBJ_OUT)$(OBJ_DIR)$Slp_data$Smatrix_scaler.$O
$(OBJ_DIR)/lp_data/matrix_utils.$O:$(SRC_DIR)/lp_data/matrix_utils.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Slp_data$Smatrix_utils.cc $(OBJ_OUT)$(OBJ_DIR)$Slp_data$Smatrix_utils.$O
$(OBJ_DIR)/lp_data/mps_reader.$O:$(SRC_DIR)/lp_data/mps_reader.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Slp_data$Smps_reader.cc $(OBJ_OUT)$(OBJ_DIR)$Slp_data$Smps_reader.$O
$(OBJ_DIR)/lp_data/mps_to_png.$O:$(SRC_DIR)/lp_data/mps_to_png.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Slp_data$Smps_to_png.cc $(OBJ_OUT)$(OBJ_DIR)$Slp_data$Smps_to_png.$O
$(OBJ_DIR)/lp_data/png_dump.$O:$(SRC_DIR)/lp_data/png_dump.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Slp_data$Spng_dump.cc $(OBJ_OUT)$(OBJ_DIR)$Slp_data$Spng_dump.$O
$(OBJ_DIR)/lp_data/sparse.$O:$(SRC_DIR)/lp_data/sparse.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Slp_data$Ssparse.cc $(OBJ_OUT)$(OBJ_DIR)$Slp_data$Ssparse.$O
$(OBJ_DIR)/lp_data/sparse_column.$O:$(SRC_DIR)/lp_data/sparse_column.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Slp_data$Ssparse_column.cc $(OBJ_OUT)$(OBJ_DIR)$Slp_data$Ssparse_column.$O
GLOP_LIB_OBJS= $(LP_DATA_OBJS) \
$(OBJ_DIR)/glop/basis_representation.$O \
$(OBJ_DIR)/glop/dual_edge_norms.$O \
$(OBJ_DIR)/glop/entering_variable.$O \
$(OBJ_DIR)/glop/initial_basis.$O \
$(OBJ_DIR)/glop/lp_data.$O \
$(OBJ_DIR)/glop/lp_print_utils.$O \
$(OBJ_DIR)/glop/lp_solver.$O \
$(OBJ_DIR)/glop/lp_types.$O \
$(OBJ_DIR)/glop/lp_utils.$O \
$(OBJ_DIR)/glop/lu_factorization.$O \
$(OBJ_DIR)/glop/markowitz.$O \
$(OBJ_DIR)/glop/matrix_scaler.$O \
$(OBJ_DIR)/glop/matrix_utils.$O \
$(OBJ_DIR)/glop/mps_reader.$O \
$(OBJ_DIR)/glop/parameters.pb.$O \
$(OBJ_DIR)/glop/preprocessor.$O \
$(OBJ_DIR)/glop/primal_edge_norms.$O \
$(OBJ_DIR)/glop/proto_utils.$O \
$(OBJ_DIR)/glop/reduced_costs.$O \
$(OBJ_DIR)/glop/revised_simplex.$O \
$(OBJ_DIR)/glop/sparse.$O \
$(OBJ_DIR)/glop/sparse_column.$O \
$(OBJ_DIR)/glop/status.$O \
$(OBJ_DIR)/glop/update_row.$O \
$(OBJ_DIR)/glop/variables_info.$O \
@@ -934,42 +971,15 @@ $(OBJ_DIR)/glop/entering_variable.$O:$(SRC_DIR)/glop/entering_variable.cc
$(OBJ_DIR)/glop/initial_basis.$O:$(SRC_DIR)/glop/initial_basis.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Sinitial_basis.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Sinitial_basis.$O
$(OBJ_DIR)/glop/lp_data.$O:$(SRC_DIR)/glop/lp_data.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Slp_data.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Slp_data.$O
$(OBJ_DIR)/glop/lp_print_utils.$O:$(SRC_DIR)/glop/lp_print_utils.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Slp_print_utils.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Slp_print_utils.$O
$(OBJ_DIR)/glop/lp_solver.$O:$(SRC_DIR)/glop/lp_solver.cc $(GEN_DIR)/linear_solver/linear_solver2.pb.h $(GEN_DIR)/glop/parameters.pb.h
$(OBJ_DIR)/glop/lp_solver.$O:$(SRC_DIR)/glop/lp_solver.cc $(GEN_DIR)/linear_solver/linear_solver2.pb.h
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Slp_solver.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Slp_solver.$O
$(OBJ_DIR)/glop/lp_utils.$O:$(SRC_DIR)/glop/lp_utils.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Slp_utils.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Slp_utils.$O
$(OBJ_DIR)/glop/lp_types.$O:$(SRC_DIR)/glop/lp_types.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Slp_types.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Slp_types.$O
$(OBJ_DIR)/glop/lu_factorization.$O:$(SRC_DIR)/glop/lu_factorization.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Slu_factorization.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Slu_factorization.$O
$(OBJ_DIR)/glop/markowitz.$O:$(SRC_DIR)/glop/markowitz.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Smarkowitz.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Smarkowitz.$O
$(OBJ_DIR)/glop/matrix_scaler.$O:$(SRC_DIR)/glop/matrix_scaler.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Smatrix_scaler.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Smatrix_scaler.$O
$(OBJ_DIR)/glop/matrix_utils.$O:$(SRC_DIR)/glop/matrix_utils.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Smatrix_utils.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Smatrix_utils.$O
$(OBJ_DIR)/glop/mps_reader.$O:$(SRC_DIR)/glop/mps_reader.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Smps_reader.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Smps_reader.$O
$(OBJ_DIR)/glop/mps_to_png.$O:$(SRC_DIR)/glop/mps_to_png.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Smps_to_png.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Smps_to_png.$O
$(OBJ_DIR)/glop/png_dump.$O:$(SRC_DIR)/glop/png_dump.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Spng_dump.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Spng_dump.$O
$(OBJ_DIR)/glop/preprocessor.$O:$(SRC_DIR)/glop/preprocessor.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Spreprocessor.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Spreprocessor.$O
@@ -991,12 +1001,6 @@ $(OBJ_DIR)/glop/reduced_costs.$O:$(SRC_DIR)/glop/reduced_costs.cc
$(OBJ_DIR)/glop/revised_simplex.$O:$(SRC_DIR)/glop/revised_simplex.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Srevised_simplex.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Srevised_simplex.$O
$(OBJ_DIR)/glop/sparse.$O:$(SRC_DIR)/glop/sparse.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Ssparse.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Ssparse.$O
$(OBJ_DIR)/glop/sparse_column.$O:$(SRC_DIR)/glop/sparse_column.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Ssparse_column.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Ssparse_column.$O
$(OBJ_DIR)/glop/status.$O:$(SRC_DIR)/glop/status.cc
$(CCC) $(CFLAGS) -c $(SRC_DIR)$Sglop$Sstatus.cc $(OBJ_OUT)$(OBJ_DIR)$Sglop$Sstatus.$O

2
src/base/status.h
View File

@@ -45,6 +45,8 @@ struct Status {
return StrCat("ERROR #", error_code_, ": '", error_message_, "'");
}
std::string error_message() const { return error_message_; }
void IgnoreError() const {}
private:

38
src/glop/basis_representation.cc
View File

@@ -15,8 +15,8 @@
#include "glop/basis_representation.h"
#include "base/stl_util.h"
#include "glop/lp_utils.h"
#include "glop/status.h"
#include "lp_data/lp_utils.h"
namespace operations_research {
namespace glop {
@@ -184,7 +184,8 @@ BasisFactorization::BasisFactorization(const MatrixView& matrix,
max_num_updates_(0),
num_updates_(0),
eta_factorization_(),
lu_factorization_() {}
lu_factorization_(),
deterministic_time_(0.0) {}
BasisFactorization::~BasisFactorization() {}
@@ -222,7 +223,13 @@ Status BasisFactorization::ForceRefactorization() {
Clear();
MatrixView basis_matrix;
basis_matrix.PopulateFromBasis(matrix_, basis_);
return lu_factorization_.ComputeFactorization(basis_matrix);
const Status status = lu_factorization_.ComputeFactorization(basis_matrix);
const double kLuComplexityFactor = 10;
deterministic_time_ +=
kLuComplexityFactor * DeterministicTimeForFpOperations(
lu_factorization_.NumberOfEntries().value());
return status;
}
// This update formula can be derived by:
@@ -296,6 +303,7 @@ Status BasisFactorization::Update(ColIndex entering_col,
void BasisFactorization::LeftSolve(DenseRow* y) const {
SCOPED_TIME_STAT(&stats_);
RETURN_IF_NULL(y);
BumpDeterministicTimeForSolve(matrix_.num_rows().value());
if (use_middle_product_form_update_) {
lu_factorization_.LeftSolveU(y);
rank_one_factorization_.LeftSolve(y);
@@ -310,6 +318,7 @@ void BasisFactorization::LeftSolveWithNonZeros(
DenseRow* y, ColIndexVector* non_zeros) const {
SCOPED_TIME_STAT(&stats_);
RETURN_IF_NULL(y);
BumpDeterministicTimeForSolve(matrix_.num_rows().value());
if (use_middle_product_form_update_) {
lu_factorization_.LeftSolveU(y);
rank_one_factorization_.LeftSolve(y);
@@ -324,6 +333,7 @@ void BasisFactorization::LeftSolveWithNonZeros(
void BasisFactorization::RightSolve(DenseColumn* d) const {
SCOPED_TIME_STAT(&stats_);
RETURN_IF_NULL(d);
BumpDeterministicTimeForSolve(matrix_.num_rows().value());
if (use_middle_product_form_update_) {
lu_factorization_.RightSolveL(d);
rank_one_factorization_.RightSolve(d);
@@ -338,6 +348,7 @@ void BasisFactorization::RightSolveWithNonZeros(
DenseColumn* d, std::vector<RowIndex>* non_zeros) const {
SCOPED_TIME_STAT(&stats_);
RETURN_IF_NULL(d);
BumpDeterministicTimeForSolve(non_zeros->size());
if (use_middle_product_form_update_) {
lu_factorization_.RightSolveL(d);
rank_one_factorization_.RightSolve(d);
@@ -353,6 +364,7 @@ void BasisFactorization::RightSolveWithNonZeros(
DenseColumn* BasisFactorization::RightSolveForTau(ScatteredColumnReference a)
const {
SCOPED_TIME_STAT(&stats_);
BumpDeterministicTimeForSolve(matrix_.num_rows().value());
if (use_middle_product_form_update_) {
lu_factorization_.RightSolveLWithPermutedInput(a.dense_column, &tau_);
rank_one_factorization_.RightSolve(&tau_);
@@ -369,6 +381,7 @@ void BasisFactorization::LeftSolveForUnitRow(ColIndex j, DenseRow* y,
ColIndexVector* non_zeros) const {
SCOPED_TIME_STAT(&stats_);
RETURN_IF_NULL(y);
BumpDeterministicTimeForSolve(1);
ClearAndResizeVectorWithNonZeros(RowToColIndex(matrix_.num_rows()), y,
non_zeros);
@@ -414,6 +427,7 @@ void BasisFactorization::RightSolveForProblemColumn(
ColIndex col, DenseColumn* d, std::vector<RowIndex>* non_zeros) const {
SCOPED_TIME_STAT(&stats_);
RETURN_IF_NULL(d);
BumpDeterministicTimeForSolve(matrix_.column(col).num_entries().value());
if (!use_middle_product_form_update_) {
lu_factorization_.SparseRightSolve(matrix_.column(col), matrix_.num_rows(),
d);
@@ -435,12 +449,14 @@ Fractional BasisFactorization::RightSolveSquaredNorm(const SparseColumn& a)
const {
SCOPED_TIME_STAT(&stats_);
DCHECK(IsRefactorized());
BumpDeterministicTimeForSolve(a.num_entries().value());
return lu_factorization_.RightSolveSquaredNorm(a);
}
Fractional BasisFactorization::DualEdgeSquaredNorm(RowIndex row) const {
SCOPED_TIME_STAT(&stats_);
DCHECK(IsRefactorized());
BumpDeterministicTimeForSolve(1);
return lu_factorization_.DualEdgeSquaredNorm(row);
}
@@ -527,5 +543,21 @@ Fractional BasisFactorization::ComputeInfinityNormConditionNumber() const {
return ComputeInfinityNorm() * ComputeInverseInfinityNorm();
}
double BasisFactorization::DeterministicTime() const {
return deterministic_time_;
}
void BasisFactorization::BumpDeterministicTimeForSolve(int num_entries) const {
// TODO(user): Spend more time finding a good approximation here.
if (matrix_.num_rows().value() == 0) return;
const double density = static_cast<double>(num_entries) /
static_cast<double>(matrix_.num_rows().value());
deterministic_time_ +=
(1.0 + density) * DeterministicTimeForFpOperations(
lu_factorization_.NumberOfEntries().value()) +
DeterministicTimeForFpOperations(
rank_one_factorization_.num_entries().value());
}
} // namespace glop
} // namespace operations_research

17
src/glop/basis_representation.h
View File

@@ -16,12 +16,12 @@
#define OR_TOOLS_GLOP_BASIS_REPRESENTATION_H_
#include "base/logging.h"
#include "glop/lp_types.h"
#include "glop/lu_factorization.h"
#include "glop/parameters.pb.h"
#include "glop/sparse.h"
#include "glop/status.h"
#include "glop/rank_one_update.h"
#include "glop/status.h"
#include "lp_data/lp_types.h"
#include "lp_data/sparse.h"
#include "util/stats.h"
namespace operations_research {
@@ -277,6 +277,10 @@ class BasisFactorization {
}
void ResetStats() { stats_.Reset(); }
// The deterministic time used by this class. It is incremented for each
// solve and each factorization.
double DeterministicTime() const;
private:
// Return true if the submatrix of matrix_ given by basis_ is exactly the
// identity (without permutation).
@@ -288,6 +292,10 @@ class BasisFactorization {
Status MiddleProductFormUpdate(ColIndex entering_col,
RowIndex leaving_variable_row) MUST_USE_RESULT;
// Increases the deterministic time for a solve operation with a vector having
// this number of non-zero entries (it can be an approximation).
void BumpDeterministicTimeForSolve(int num_entries) const;
// Stats about this class.
struct Stats : public StatsGroup {
Stats()
@@ -331,6 +339,9 @@ class BasisFactorization {
EtaFactorization eta_factorization_;
LuFactorization lu_factorization_;
// mutable because the Solve() functions are const but need to update this.
mutable double deterministic_time_;
DISALLOW_COPY_AND_ASSIGN(BasisFactorization);
};

2
src/glop/dual_edge_norms.cc
View File

@@ -13,7 +13,7 @@
#include "glop/dual_edge_norms.h"
#include "glop/lp_utils.h"
#include "lp_data/lp_utils.h"
namespace operations_research {
namespace glop {

4
src/glop/dual_edge_norms.h
View File

@@ -15,9 +15,9 @@
#define OR_TOOLS_GLOP_DUAL_EDGE_NORMS_H_
#include "glop/basis_representation.h"
#include "glop/lp_data.h"
#include "glop/lp_types.h"
#include "glop/parameters.pb.h"
#include "lp_data/lp_data.h"
#include "lp_data/lp_types.h"
#include "util/stats.h"
namespace operations_research {

2
src/glop/entering_variable.cc
View File

@@ -17,7 +17,7 @@
#include <queue>
#include "base/timer.h"
#include "glop/lp_utils.h"
#include "lp_data/lp_utils.h"
namespace operations_research {
namespace glop {

6
src/glop/entering_variable.h
View File

@@ -15,14 +15,14 @@
#define OR_TOOLS_GLOP_ENTERING_VARIABLE_H_
#include "glop/basis_representation.h"
#include "glop/lp_data.h"
#include "glop/lp_types.h"
#include "glop/parameters.pb.h"
#include "glop/primal_edge_norms.h"
#include "glop/reduced_costs.h"
#include "glop/status.h"
#include "glop/update_row.h"
#include "glop/variables_info.h"
#include "glop/status.h"
#include "lp_data/lp_data.h"
#include "lp_data/lp_types.h"
#include "util/bitset.h"
#include "base/random.h"
#include "util/stats.h"

2
src/glop/initial_basis.cc
View File

@@ -15,8 +15,8 @@
#include <queue>
#include "glop/initial_basis.h"
#include "glop/lp_utils.h"
#include "glop/markowitz.h"
#include "lp_data/lp_utils.h"
namespace operations_research {
namespace glop {

6
src/glop/initial_basis.h
View File

@@ -14,9 +14,9 @@
#ifndef OR_TOOLS_GLOP_INITIAL_BASIS_H_
#define OR_TOOLS_GLOP_INITIAL_BASIS_H_
#include "glop/lp_data.h"
#include "glop/lp_types.h"
#include "glop/sparse.h"
#include "lp_data/lp_data.h"
#include "lp_data/lp_types.h"
#include "lp_data/sparse.h"
namespace operations_research {
namespace glop {

194
src/glop/lp_solver.cc
View File

@@ -23,12 +23,12 @@
#include "base/join.h"
#include "base/strutil.h"
#include "glop/lp_types.h"
#include "glop/lp_utils.h"
#include "glop/preprocessor.h"
#include "glop/proto_utils.h"
#include "glop/status.h"
#include "linear_solver/proto_tools.h"
#include "lp_data/lp_types.h"
#include "lp_data/lp_utils.h"
#include "util/fp_utils.h"
DEFINE_bool(lp_solver_enable_fp_exceptions, true,
@@ -126,7 +126,7 @@ ProblemStatus LPSolver::Solve(const LinearProgram& lp) {
ScopedFloatingPointEnv scoped_fenv;
if (FLAGS_lp_solver_enable_fp_exceptions) {
#ifdef _MSC_VER
scoped_fenv.EnableExceptions(_EM_INVALID | _EM_ZERODIVIDE);
scoped_fenv.EnableExceptions(_EM_INVALID | EM_ZERODIVIDE);
#else
scoped_fenv.EnableExceptions(FE_DIVBYZERO | FE_INVALID);
#endif
@@ -174,7 +174,6 @@ ProblemStatus LPSolver::LoadAndVerifySolution(const LinearProgram& lp,
constraint_statuses_ = solution.constraint_statuses;
status_ = solution.status;
// Verify solution.
if (status_ == ProblemStatus::OPTIMAL &&
parameters_.provide_strong_optimal_guarantee()) {
MovePrimalValuesWithinBounds(lp);
@@ -182,31 +181,49 @@ ProblemStatus LPSolver::LoadAndVerifySolution(const LinearProgram& lp,
}
ComputeReducedCosts(lp);
ComputeConstraintActivities(lp);
max_absolute_primal_infeasibility_ = 0.0;
max_absolute_dual_infeasibility_ = 0.0;
ComputeColumnInfeasibilities(lp);
ComputeRowInfeasibilities(lp);
ComputeObjective(lp);
ComputeDualObjective(lp);
// Note that this needs primal_values_ and reduced_costs_ to be up to date.
if (status_ == ProblemStatus::OPTIMAL &&
parameters_.provide_strong_optimal_guarantee()) {
// Note that this call may modify the solution values.
Fractional cost_perturbation =
GetMaxCostPerturbationToEnforceOptimality(lp);
VLOG(1) << "Cost perturbation = " << cost_perturbation
<< " (Difference with dual infeasibility = "
<< cost_perturbation - max_absolute_dual_infeasibility_ << ")";
VLOG(1) << "Cost perturbation = " << cost_perturbation;
if (cost_perturbation > parameters_.solution_feasibility_tolerance()) {
VLOG(1) << "The needed cost perturbation is too high !!";
status_ = ProblemStatus::IMPRECISE;
}
VLOG(1) << "Rhs perturbation = " << max_absolute_primal_infeasibility_
<< " (Same as primal infeasibility)";
}
ComputeConstraintActivities(lp);
const double primal_infeasibility = ComputePrimalValueInfeasibility(lp);
const double dual_infeasibility = ComputeDualValueInfeasibility(lp);
const double primal_residual = ComputeActivityInfeasibility(lp);
const double dual_residual = ComputeReducedCostInfeasibility(lp);
const double objective_value = ComputeObjective(lp);
const double dual_objective_value = ComputeDualObjective(lp);
objective_value_with_offset_ = objective_value + lp.objective_offset();
// If the primal/dual values were moved to the bounds, then the primal/dual
// infeasibilities should be exactly zero (but not the residuals).
if (status_ == ProblemStatus::OPTIMAL &&
parameters_.provide_strong_optimal_guarantee()) {
if (primal_infeasibility != 0.0 || dual_infeasibility != 0.0) {
VLOG(1) << "If you see this, there is a bug in "
"MovePrimalValuesWithinBounds() or in "
"MoveDualValuesWithinBounds().";
}
}
// TODO(user): the name is not really consistent since in practice those are
// the "residual" since the primal/dual infeasibility are zero when
// parameters_.provide_strong_optimal_guarantee() is true.
max_absolute_primal_infeasibility_ =
std::max(primal_infeasibility, primal_residual);
max_absolute_dual_infeasibility_ =
std::max(dual_infeasibility, dual_residual);
// Check precision and optimality of the result. See Chvatal pp. 61-62.
//
// For that we check that:
@@ -223,7 +240,8 @@ ProblemStatus LPSolver::LoadAndVerifySolution(const LinearProgram& lp,
VLOG(1) << "Max. primal infeasibility = "
<< max_absolute_primal_infeasibility_;
VLOG(1) << "Max. dual infeasibility = " << max_absolute_dual_infeasibility_;
VLOG(1) << "Gap = " << objective_value_ - dual_objective_value_;
VLOG(1) << "Objective (with offset) = " << objective_value_with_offset_;
VLOG(1) << "Objective gap = " << objective_value - dual_objective_value;
if ((status_ == ProblemStatus::OPTIMAL ||
status_ == ProblemStatus::PRIMAL_FEASIBLE) &&
@@ -241,7 +259,7 @@ ProblemStatus LPSolver::LoadAndVerifySolution(const LinearProgram& lp,
}
if (status_ == ProblemStatus::OPTIMAL) {
if (!AreWithinAbsoluteOrRelativeTolerances(
objective_value_, dual_objective_value_,
objective_value, dual_objective_value,
parameters_.solution_objective_gap_tolerance(),
parameters_.solution_objective_gap_tolerance())) {
VLOG(1) << "The final solution objective gap is too high !!";
@@ -311,6 +329,11 @@ int LPSolver::GetNumberOfSimplexIterations() const {
return num_revised_simplex_iterations_;
}
double LPSolver::DeterministicTime() const {
return revised_simplex_ == nullptr ? 0.0
: revised_simplex_->DeterministicTime();
}
void LPSolver::MovePrimalValuesWithinBounds(const LinearProgram& lp) {
const ColIndex num_cols = lp.num_variables();
DCHECK_EQ(num_cols, primal_values_.size());
@@ -636,12 +659,12 @@ Fractional LPSolver::GetMaxCostPerturbationToEnforceOptimality(
DCHECK_LE(variable_value, upper_bound);
// We correct the reduced cost, so we have a minimization problem and
// thus the dual_objective_value_ will be a lower bound of the primal
// thus the dual objective value will be a lower bound of the primal
// objective.
const Fractional reduced_cost = optimization_sign * reduced_costs_[col];
// Perturbation needed to move the cost to 0.0
const Fractional cost_delta = fabs(reduced_cost);
const Fractional cost_delta = std::abs(reduced_cost);
// We are enforcing the complementary slackness conditions, see the comment
// in ComputeDualObjective() for more explanations. By default, we move the
@@ -710,18 +733,17 @@ void LPSolver::ComputeReducedCosts(const LinearProgram& lp) {
}
}
void LPSolver::ComputeObjective(const LinearProgram& lp) {
double LPSolver::ComputeObjective(const LinearProgram& lp) {
const ColIndex num_cols = lp.num_variables();
DCHECK_EQ(num_cols, primal_values_.size());
KahanSum sum;
for (ColIndex col(0); col < num_cols; ++col) {
sum.Add(lp.objective_coefficients()[col] * primal_values_[col]);
}
objective_value_ = sum.Value();
objective_value_with_offset_ = objective_value_ + lp.objective_offset();
return sum.Value();
}
void LPSolver::ComputeDualObjective(const LinearProgram& lp) {
double LPSolver::ComputeDualObjective(const LinearProgram& lp) {
// By the duality theorem, the dual objective is a bound on the primal
// objective obtained by taking the linear combinaison of the constraints
// given by dual_values_.
@@ -794,36 +816,53 @@ void LPSolver::ComputeDualObjective(const LinearProgram& lp) {
// We now need to apply the correction in the good direction!
dual_objective.Add(optimization_sign * correction);
}
dual_objective_value_ = dual_objective.Value();
return dual_objective.Value();
}
void LPSolver::UpdateMaxPrimalInfeasibility(Fractional update) {
max_absolute_primal_infeasibility_ =
std::max(max_absolute_primal_infeasibility_, update);
double LPSolver::ComputePrimalValueInfeasibility(const LinearProgram& lp) {
double infeasibility = 0.0;
const ColIndex num_cols = lp.num_variables();
for (ColIndex col(0); col < num_cols; ++col) {
const Fractional lower_bound = lp.variable_lower_bounds()[col];
const Fractional upper_bound = lp.variable_upper_bounds()[col];
DCHECK(IsFinite(primal_values_[col]));
if (lower_bound == upper_bound) {
infeasibility =
std::max(infeasibility, std::abs(primal_values_[col] - upper_bound));
continue;
}
if (primal_values_[col] > upper_bound) {
infeasibility =
std::max(infeasibility, primal_values_[col] - upper_bound);
}
if (primal_values_[col] < lower_bound) {
infeasibility =
std::max(infeasibility, lower_bound - primal_values_[col]);
}
}
return infeasibility;
}
void LPSolver::UpdateMaxDualInfeasibility(Fractional update) {
max_absolute_dual_infeasibility_ =
std::max(max_absolute_dual_infeasibility_, update);
}
void LPSolver::ComputeRowInfeasibilities(const LinearProgram& lp) {
RowIndex num_problematic_rows(0);
double LPSolver::ComputeActivityInfeasibility(const LinearProgram& lp) {
double infeasibility = 0.0;
int num_problematic_rows(0);
const RowIndex num_rows = lp.num_constraints();
const Fractional optimization_sign = lp.IsMaximizationProblem() ? -1.0 : 1.0;
for (RowIndex row(0); row < num_rows; ++row) {
const Fractional activity = constraint_activities_[row];
const Fractional lower_bound = lp.constraint_lower_bounds()[row];
const Fractional upper_bound = lp.constraint_upper_bounds()[row];
const Fractional tolerance = parameters_.solution_feasibility_tolerance();
DCHECK(IsFinite(activity));
if (lower_bound == upper_bound) {
if (fabs(activity - upper_bound) > tolerance) {
if (std::abs(activity - upper_bound) > tolerance) {
VLOG(2) << "Row " << row.value() << " has activity " << activity
<< " which is different from " << upper_bound << " by "
<< activity - upper_bound;
++num_problematic_rows;
}
UpdateMaxPrimalInfeasibility(fabs(activity - upper_bound));
infeasibility = std::max(infeasibility, std::abs(activity - upper_bound));
continue;
}
if (activity > upper_bound) {
@@ -834,7 +873,7 @@ void LPSolver::ComputeRowInfeasibilities(const LinearProgram& lp) {
<< row_excess;
++num_problematic_rows;
}
UpdateMaxPrimalInfeasibility(row_excess);
infeasibility = std::max(infeasibility, row_excess);
}
if (activity < lower_bound) {
const Fractional row_deficit = lower_bound - activity;
@@ -844,53 +883,52 @@ void LPSolver::ComputeRowInfeasibilities(const LinearProgram& lp) {
<< row_deficit;
++num_problematic_rows;
}
UpdateMaxPrimalInfeasibility(row_deficit);
}
// For a minimization problem, we want a lower bound.
const Fractional minimization_dual_value =
optimization_sign * dual_values_[row];
if (lower_bound == -kInfinity) {
UpdateMaxDualInfeasibility(minimization_dual_value);
}
if (upper_bound == kInfinity) {
UpdateMaxDualInfeasibility(-minimization_dual_value);
infeasibility = std::max(infeasibility, row_deficit);
}
}
VLOG(1) << "Number of problematic rows " << num_problematic_rows.value();
VLOG(1) << "Number of problematic rows: " << num_problematic_rows;
return infeasibility;
}
void LPSolver::ComputeColumnInfeasibilities(const LinearProgram& lp) {
const ColIndex num_cols = lp.num_variables();
double LPSolver::ComputeDualValueInfeasibility(const LinearProgram& lp) {
const Fractional optimization_sign = lp.IsMaximizationProblem() ? -1.0 : 1.0;
for (ColIndex col(0); col < num_cols; ++col) {
const Fractional lower_bound = lp.variable_lower_bounds()[col];
const Fractional upper_bound = lp.variable_upper_bounds()[col];
DCHECK(IsFinite(primal_values_[col]));
DCHECK(IsFinite(reduced_costs_[col]));
// Compute the infeasibility on the bound.
if (lower_bound == upper_bound) {
UpdateMaxPrimalInfeasibility(fabs(primal_values_[col] - upper_bound));
continue;
}
if (primal_values_[col] > upper_bound) {
UpdateMaxPrimalInfeasibility(primal_values_[col] - upper_bound);
}
if (primal_values_[col] < lower_bound) {
UpdateMaxPrimalInfeasibility(lower_bound - primal_values_[col]);
}
// For a minimization problem, we want a lower bound.
const Fractional minimization_reduced_cost =
optimization_sign * reduced_costs_[col];
double infeasibility = 0.0;
const RowIndex num_rows = lp.num_constraints();
for (RowIndex row(0); row < num_rows; ++row) {
const Fractional dual_value = dual_values_[row];
const Fractional lower_bound = lp.constraint_lower_bounds()[row];
const Fractional upper_bound = lp.constraint_upper_bounds()[row];
DCHECK(IsFinite(dual_value));
const Fractional minimization_dual_value = optimization_sign * dual_value;
if (lower_bound == -kInfinity) {
UpdateMaxDualInfeasibility(minimization_reduced_cost);
infeasibility = std::max(infeasibility, minimization_dual_value);
}
if (upper_bound == kInfinity) {
UpdateMaxDualInfeasibility(-minimization_reduced_cost);
infeasibility = std::max(infeasibility, -minimization_dual_value);
}
}
return infeasibility;
}
double LPSolver::ComputeReducedCostInfeasibility(const LinearProgram& lp) {
const Fractional optimization_sign = lp.IsMaximizationProblem() ? -1.0 : 1.0;
double infeasibility = 0.0;
const ColIndex num_cols = lp.num_variables();
for (ColIndex col(0); col < num_cols; ++col) {
const Fractional reduced_cost = reduced_costs_[col];
const Fractional lower_bound = lp.variable_lower_bounds()[col];
const Fractional upper_bound = lp.variable_upper_bounds()[col];
DCHECK(IsFinite(reduced_cost));
const Fractional minimization_reduced_cost =
optimization_sign * reduced_cost;
if (lower_bound == -kInfinity) {
infeasibility = std::max(infeasibility, minimization_reduced_cost);
}
if (upper_bound == kInfinity) {
infeasibility = std::max(infeasibility, -minimization_reduced_cost);
}
}
return infeasibility;
}
} // namespace glop

37
src/glop/lp_solver.h
View File

@@ -16,10 +16,10 @@
#include "base/unique_ptr.h"
#include "glop/lp_data.h"
#include "glop/lp_types.h"
#include "glop/parameters.pb.h"
#include "glop/preprocessor.h"
#include "lp_data/lp_data.h"
#include "lp_data/lp_types.h"
#include "util/time_limit.h"
namespace operations_research {
@@ -113,6 +113,16 @@ class LPSolver {
// Returns the number of simplex iterations used by the last Solve().
int GetNumberOfSimplexIterations() const;
// Returns the "deterministic time" since the creation of the solver. Note
// That this time is only increased when some operations take place in this
// class.
//
// TODO(user): Currently, this is only modified when the simplex code is
// executed.
//
// TODO(user): Improve the correlation with the running time.
double DeterministicTime() const;
private:
// Resizes all the solution vectors to the given sizes.
// This is used in case of error to make sure all the getter functions will
@@ -175,20 +185,21 @@ class LPSolver {
Fractional GetMaxCostPerturbationToEnforceOptimality(const LinearProgram& lp);
// Computes derived quantities from the solution.
void ComputeObjective(const LinearProgram& lp);
void ComputeDualObjective(const LinearProgram& lp);
void ComputeReducedCosts(const LinearProgram& lp);
void ComputeConstraintActivities(const LinearProgram& lp);
// Shortcuts that just take the max of the current value and the given update.
// Since the initial infeasibility value is 0.0, a negative update is
// not considered an infeasibility.
void UpdateMaxPrimalInfeasibility(Fractional update);
void UpdateMaxDualInfeasibility(Fractional update);
// Computes the primal/dual objectives (without the offset). Note that the
// dual objective needs the reduced costs in addition to the dual values.
double ComputeObjective(const LinearProgram& lp);
double ComputeDualObjective(const LinearProgram& lp);
// Compute the primal/dual infeasibilities.
void ComputeRowInfeasibilities(const LinearProgram& lp);
void ComputeColumnInfeasibilities(const LinearProgram& lp);
// Computes the maximum of the infeasibilities associated with each values.
// The returned infeasibilities are the maximum of the "absolute" errors of
// each vector coefficients.
double ComputePrimalValueInfeasibility(const LinearProgram& lp);
double ComputeActivityInfeasibility(const LinearProgram& lp);
double ComputeDualValueInfeasibility(const LinearProgram& lp);
double ComputeReducedCostInfeasibility(const LinearProgram& lp);
// Dimension of the linear program given to the last Solve().
// This is used for displaying purpose only.
@@ -221,9 +232,7 @@ class LPSolver {
// Quantities computed from the solution and the linear program.
DenseRow reduced_costs_;
DenseColumn constraint_activities_;
Fractional objective_value_;
Fractional objective_value_with_offset_;
Fractional dual_objective_value_;
bool may_have_multiple_solutions_;
Fractional max_absolute_primal_infeasibility_;
Fractional max_absolute_dual_infeasibility_;

7
src/glop/lu_factorization.cc
View File

@@ -13,7 +13,7 @@
#include "glop/lu_factorization.h"
#include "glop/lp_utils.h"
#include "lp_data/lp_utils.h"
namespace operations_research {
namespace glop {
@@ -420,6 +420,11 @@ double LuFactorization::GetFillInPercentage(const MatrixView& matrix) const {
static_cast<double>(initial_num_entries);
}
EntryIndex LuFactorization::NumberOfEntries() const {
return is_identity_factorization_ ? EntryIndex(0) : lower_.num_entries() +
upper_.num_entries();
}
Fractional LuFactorization::ComputeDeterminant() const {
if (is_identity_factorization_) return 1.0;
DCHECK_EQ(upper_.num_rows().value(), upper_.num_cols().value());

16
src/glop/lu_factorization.h
View File

@@ -11,15 +11,15 @@
// See the License for the specific language governing permissions and
// limitations under the License.
#include "glop/markowitz.h"
#include "glop/parameters.pb.h"
#include "glop/sparse.h"
#include "glop/status.h"
#include "util/stats.h"
#ifndef OR_TOOLS_GLOP_LU_FACTORIZATION_H_
#define OR_TOOLS_GLOP_LU_FACTORIZATION_H_
#include "glop/markowitz.h"
#include "glop/parameters.pb.h"
#include "glop/status.h"
#include "lp_data/sparse.h"
#include "util/stats.h"
namespace operations_research {
namespace glop {
@@ -154,6 +154,10 @@ class LuFactorization {
// the number of entries in B.
double GetFillInPercentage(const MatrixView& matrix) const;
// Returns the number of entries in L + U.
// If the factorization is the identity, this returns 0.
EntryIndex NumberOfEntries() const;
// Computes the determinant of the input matrix B.
// Since P.B.Q^{-1} = L.U, det(P) * det(B) * det(Q^{-1}) = det(L) * det(U).
// det(L) = 1 since L is a lower-triangular matrix with 1 on the diagonal.

2
src/glop/markowitz.cc
View File

@@ -15,7 +15,7 @@
#include <limits>
#include "base/stringprintf.h"
#include "glop/lp_utils.h"
#include "lp_data/lp_utils.h"
namespace operations_research {
namespace glop {

6
src/glop/markowitz.h
View File

@@ -76,10 +76,10 @@
#include <queue>
#include "base/logging.h"
#include "glop/lp_types.h"
#include "glop/sparse.h"
#include "glop/status.h"
#include "glop/parameters.pb.h"
#include "glop/status.h"
#include "lp_data/lp_types.h"
#include "lp_data/sparse.h"
#include "util/stats.h"
namespace operations_research {

161
src/glop/mps_driver.cc
View File

@@ -1,161 +0,0 @@
// Copyright 2010-2014 Google
// 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.
// Driver for reading and solving files in the MPS format and in
// the linear_solver2.proto format.
#include <stdio.h>
#include <string>
#include "base/commandlineflags.h"
#include "base/commandlineflags.h"
#include "base/logging.h"
#include "base/timer.h"
#include "base/file.h"
#include "google/protobuf/descriptor.h"
#include "google/protobuf/message.h"
#include "google/protobuf/text_format.h"
#include "base/strutil.h"
#include "glop/lp_print_utils.h"
#include "glop/lp_solver.h"
#include "glop/mps_reader.h"
#include "glop/parameters.pb.h"
#include "glop/proto_utils.h"
#include "linear_solver/proto_tools.h"
#include "base/status.h"
DEFINE_bool(mps_dump_problem, false, "Dumps problem in readable form.");
DEFINE_bool(mps_solve, true, "Solves problem.");
DEFINE_bool(mps_terse_result, false,
"Displays the result in form of a single CSV line.");
DEFINE_bool(mps_verbose_result, true, "Displays the result in verbose form.");
DEFINE_bool(mps_display_full_path, true,
"Displays the full path of the input file in the result line.");
DEFINE_string(input, "", "File pattern for problems to be optimized.");
DEFINE_string(params_file, "", "Path to a GlopParameters file in text format.");
DEFINE_string(params, "",
"GlopParameters in text format. If --params_file was "
"also specified, the --params will be merged onto "
"them (i.e. in case of conflicts, --params wins)");
using operations_research::FullProtocolMessageAsString;
using operations_research::ReadFileToProto;
using operations_research::glop::GetProblemStatusString;
using operations_research::glop::GlopParameters;
using operations_research::glop::LinearProgram;
using operations_research::glop::LPSolver;
using operations_research::glop::MPSReader;
using operations_research::glop::MPModelProtoToLinearProgram;
using operations_research::glop::ProblemStatus;
using operations_research::glop::ToDouble;
using google::protobuf::Descriptor;
using google::protobuf::FieldDescriptor;
using google::protobuf::Message;
using google::protobuf::Reflection;
using google::protobuf::TextFormat;
using operations_research::HasSuffixString;
using operations_research::ScopedWallTime;
// Parse glop parameters from the flags --params_file and --params.
void ReadGlopParameters(GlopParameters* parameters) {
if (!FLAGS_params_file.empty()) {
std::string params;
CHECK(TextFormat::ParseFromString(params, parameters)) << params;
}
if (!FLAGS_params.empty()) {
CHECK(TextFormat::MergeFromString(FLAGS_params, parameters))
<< FLAGS_params;
}
if (FLAGS_mps_verbose_result) {
printf("GlopParameters {\n%s}\n",
FullProtocolMessageAsString(*parameters, 1).c_str());
}
}
int main(int argc, char* argv[]) {
google::ParseCommandLineFlags(&argc, &argv, true);
GlopParameters parameters;
ReadGlopParameters(&parameters);
LinearProgram linear_program;
std::vector<std::string> file_list;
// Replace this with your favorite match function.
file_list.push_back(FLAGS_input);
for (int i = 0; i < file_list.size(); ++i) {
const std::string& file_name = file_list[i];
MPSReader mps_reader;
operations_research::new_proto::MPModelProto model_proto;
if (HasSuffixString(file_name, ".mps") ||
HasSuffixString(file_name, ".mps.gz")) {
if (!mps_reader.LoadFileAndTryFreeFormOnFail(file_name,
&linear_program)) {
LOG(INFO) << "Parse error for " << file_name;
continue;
}
} else {
ReadFileToProto(file_name, &model_proto);
MPModelProtoToLinearProgram(model_proto, &linear_program);
}
if (FLAGS_mps_dump_problem) {
printf("%s", linear_program.Dump().c_str());
}
// Create the solver with the correct parameters.
LPSolver solver;
solver.SetParameters(parameters);
ProblemStatus solve_status = ProblemStatus::INIT;
const char* status_string;
double objective_value;
double solving_time_in_sec = 0;
if (FLAGS_mps_solve) {
ScopedWallTime timer(&solving_time_in_sec);
solve_status = solver.Solve(linear_program);
status_string = GetProblemStatusString(solve_status).c_str();
objective_value = ToDouble(solver.GetObjectiveValue());
}
if (FLAGS_mps_terse_result) {
if (FLAGS_mps_display_full_path) {
printf("%s,", file_name.c_str());
}
printf("%s,", mps_reader.GetProblemName().c_str());
if (FLAGS_mps_solve) {
printf("%15.15e,%s,%-6.4g,", objective_value, status_string,
solving_time_in_sec);
}
printf("%s,%s\n", linear_program.GetProblemStats().c_str(),
linear_program.GetNonZeroStats().c_str());
}
if (FLAGS_mps_verbose_result) {
if (FLAGS_mps_display_full_path) {
printf("%-45s: %s\n", "File path", file_name.c_str());
}
printf("%-45s: %s\n", "Problem name",
mps_reader.GetProblemName().c_str());
if (FLAGS_mps_solve) {
printf("%-45s: %15.15e\n", "Objective value", objective_value);
printf("%-45s: %s\n", "Problem status", status_string);
printf("%-45s: %-6.4g\n", "Solving time", solving_time_in_sec);
}
printf("%s%s", linear_program.GetPrettyProblemStats().c_str(),
linear_program.GetPrettyNonZeroStats().c_str());
}
}
return EXIT_SUCCESS;
}

7
src/glop/parameters.proto
View File

@@ -224,6 +224,13 @@ message GlopParameters {
// Maximum time allowed in seconds to solve a problem.
optional double max_time_in_seconds = 26 [default = inf];
// Maximum deterministic time allowed to solve a problem. The deterministic
// time is more or less correlated to the running time, and its unit should
// be around the second (at least on a Xeon(R) CPU E5-1650 v2 @ 3.50GHz).
//
// TODO(user): Improve the correlation.
optional double max_deterministic_time = 45 [default = inf];
// Maximum number of simplex iterations to solve a problem.
// A value of -1 means no limit.
optional int64 max_number_of_iterations = 27 [default = -1];

73
src/glop/preprocessor.cc
View File

@@ -15,10 +15,10 @@
#include "glop/preprocessor.h"
#include "base/stringprintf.h"
#include "glop/lp_utils.h"
#include "glop/matrix_utils.h"
#include "glop/revised_simplex.h"
#include "glop/status.h"
#include "lp_data/lp_utils.h"
#include "lp_data/matrix_utils.h"
namespace operations_research {
namespace glop {
@@ -1817,11 +1817,9 @@ void EmptyConstraintPreprocessor::StoreSolution(
// --------------------------------------------------------
SingletonUndo::SingletonUndo(OperationType type, const LinearProgram& lp,
const SparseColumn& sparse_column_or_row,
MatrixEntry e, ConstraintStatus status)
: type_(type),
is_maximization_(lp.IsMaximizationProblem()),
sparse_column_or_row_(sparse_column_or_row),
e_(e),
cost_(lp.objective_coefficients()[e.col]),
variable_lower_bound_(lp.variable_lower_bounds()[e.col]),
@@ -1830,16 +1828,18 @@ SingletonUndo::SingletonUndo(OperationType type, const LinearProgram& lp,
constraint_upper_bound_(lp.constraint_upper_bounds()[e.row]),
constraint_status_(status) {}
void SingletonUndo::Undo(ProblemSolution* solution) const {
void SingletonUndo::Undo(const SparseMatrix& deleted_columns,
const SparseMatrix& deleted_rows,
ProblemSolution* solution) const {
switch (type_) {
case SINGLETON_ROW:
SingletonRowUndo(solution);
SingletonRowUndo(deleted_columns, solution);
break;
case ZERO_COST_SINGLETON_COLUMN:
ZeroCostSingletonColumnUndo(solution);
ZeroCostSingletonColumnUndo(deleted_rows, solution);
break;
case SINGLETON_COLUMN_IN_EQUALITY:
SingletonColumnInEqualityUndo(solution);
SingletonColumnInEqualityUndo(deleted_rows, solution);
break;
case MAKE_CONSTRAINT_AN_EQUALITY:
MakeConstraintAnEqualityUndo(solution);
@@ -1878,14 +1878,18 @@ void SingletonPreprocessor::DeleteSingletonRow(MatrixEntry e,
DCHECK_EQ(new_lower_bound, new_upper_bound);
}
row_deletion_helper_.MarkRowForDeletion(e.row);
undo_stack_.push_back(SingletonUndo(SingletonUndo::SINGLETON_ROW, *lp,
lp->GetSparseColumn(e.col), e,
undo_stack_.push_back(SingletonUndo(SingletonUndo::SINGLETON_ROW, *lp, e,
ConstraintStatus::FREE));
if (deleted_columns_.column(e.col).IsEmpty()) {
deleted_columns_.mutable_column(e.col)
->PopulateFromSparseVector(lp->GetSparseColumn(e.col));
}
lp->SetVariableBounds(e.col, new_lower_bound, new_upper_bound);
}
// The dual value of the row needs to be corrected to stay at the optimal.
void SingletonUndo::SingletonRowUndo(ProblemSolution* solution) const {
void SingletonUndo::SingletonRowUndo(const SparseMatrix& deleted_columns,
ProblemSolution* solution) const {
DCHECK_EQ(0, solution->dual_values[e_.row]);
// If the variable is basic or free, we can just keep the constraint
@@ -1910,7 +1914,8 @@ void SingletonUndo::SingletonRowUndo(ProblemSolution* solution) const {
// This is the reduced cost of the variable before the singleton constraint is
// added back.
const Fractional reduced_cost =
cost_ - ScalarProduct(solution->dual_values, sparse_column_or_row_);
cost_ -
ScalarProduct(solution->dual_values, deleted_columns.column(e_.col));
const Fractional reduced_cost_for_minimization =
is_maximization_ ? -reduced_cost : reduced_cost;
@@ -1964,16 +1969,21 @@ void SingletonPreprocessor::UpdateConstraintBoundsWithVariableBounds(
void SingletonPreprocessor::DeleteZeroCostSingletonColumn(
const SparseMatrix& transpose, MatrixEntry e, LinearProgram* lp) {
undo_stack_.push_back(SingletonUndo(
SingletonUndo::ZERO_COST_SINGLETON_COLUMN, *lp,
transpose.column(RowToColIndex(e.row)), e, ConstraintStatus::FREE));
const ColIndex transpose_col = RowToColIndex(e.row);
const SparseColumn& column = transpose.column(transpose_col);
undo_stack_.push_back(SingletonUndo(SingletonUndo::ZERO_COST_SINGLETON_COLUMN,
*lp, e, ConstraintStatus::FREE));
if (deleted_rows_.column(transpose_col).IsEmpty()) {
deleted_rows_.mutable_column(transpose_col)
->PopulateFromSparseVector(column);
}
UpdateConstraintBoundsWithVariableBounds(e, lp);
column_deletion_helper_.MarkColumnForDeletion(e.col);
}
// We need to restore the variable value in order to satisfy the constraint.
void SingletonUndo::ZeroCostSingletonColumnUndo(
ProblemSolution* solution) const {
const SparseMatrix& deleted_rows, ProblemSolution* solution) const {
// If the variable was fixed, this is easy. Note that this is the only
// possible case if the current constraint status is FIXED.
if (variable_upper_bound_ == variable_lower_bound_) {
@@ -2004,8 +2014,9 @@ void SingletonUndo::ZeroCostSingletonColumnUndo(
// This is the activity of the constraint before the singleton variable is
// added back to it.
const ColIndex row_as_col = RowToColIndex(e_.row);
const Fractional activity =
ScalarProduct(solution->primal_values, sparse_column_or_row_);
ScalarProduct(solution->primal_values, deleted_rows.column(row_as_col));
// TODO(user): use parameters_ here (it needs to be passed to SingletonUndo)
const Fractional tolerance = 1e-10;
@@ -2080,12 +2091,13 @@ void SingletonUndo::ZeroCostSingletonColumnUndo(
void SingletonPreprocessor::DeleteSingletonColumnInEquality(
const SparseMatrix& transpose, MatrixEntry e, LinearProgram* lp) {
// Save information for the undo.
// TODO(user): This can be shared between the singleton columns on the same
// row.
const SparseColumn& column = transpose.column(RowToColIndex(e.row));
const ColIndex transpose_col = RowToColIndex(e.row);
const SparseColumn& row_as_column = transpose.column(transpose_col);
undo_stack_.push_back(
SingletonUndo(SingletonUndo::SINGLETON_COLUMN_IN_EQUALITY, *lp, column, e,
SingletonUndo(SingletonUndo::SINGLETON_COLUMN_IN_EQUALITY, *lp, e,
ConstraintStatus::FREE));
deleted_rows_.mutable_column(transpose_col)
->PopulateFromSparseVector(row_as_column);
// Update the objective function using the equality constraint. We have
// v_col*coeff + expression = rhs,
@@ -2096,7 +2108,7 @@ void SingletonPreprocessor::DeleteSingletonColumnInEquality(
const Fractional cost = lp->objective_coefficients()[e.col];
const Fractional multiplier = cost / e.coeff;
lp->SetObjectiveOffset(lp->objective_offset() + rhs * multiplier);
for (const SparseColumn::Entry e : column) {
for (const SparseColumn::Entry e : row_as_column) {
const ColIndex col = RowToColIndex(e.row());
if (!column_deletion_helper_.IsColumnMarked(col)) {
Fractional new_cost =
@@ -2121,9 +2133,9 @@ void SingletonPreprocessor::DeleteSingletonColumnInEquality(
}
void SingletonUndo::SingletonColumnInEqualityUndo(
ProblemSolution* solution) const {
const SparseMatrix& deleted_rows, ProblemSolution* solution) const {
// First do the same as a zero-cost singleton column.
ZeroCostSingletonColumnUndo(solution);
ZeroCostSingletonColumnUndo(deleted_rows, solution);
// Then, restore the dual optimal value taking into account the cost
// modification.
@@ -2253,10 +2265,8 @@ bool SingletonPreprocessor::MakeConstraintAnEqualityIfPossible(
if (lp->constraint_lower_bounds()[e.row] ==
lp->constraint_upper_bounds()[e.row]) {
SparseColumn unused;
undo_stack_.push_back(
SingletonUndo(SingletonUndo::MAKE_CONSTRAINT_AN_EQUALITY, *lp, unused,
e, relaxed_status));
undo_stack_.push_back(SingletonUndo(
SingletonUndo::MAKE_CONSTRAINT_AN_EQUALITY, *lp, e, relaxed_status));
return true;
}
return false;
@@ -2269,6 +2279,10 @@ bool SingletonPreprocessor::Run(LinearProgram* lp) {
// Initialize column_to_process with the current singleton columns.
ColIndex num_cols(matrix.num_cols());
RowIndex num_rows(matrix.num_rows());
deleted_columns_.PopulateFromZero(num_rows, num_cols);
deleted_rows_.PopulateFromZero(ColToRowIndex(num_cols),
RowToColIndex(num_rows));
StrictITIVector<ColIndex, EntryIndex> column_degree(num_cols, EntryIndex(0));
std::vector<ColIndex> column_to_process;
for (ColIndex col(0); col < num_cols; ++col) {
@@ -2279,7 +2293,6 @@ bool SingletonPreprocessor::Run(LinearProgram* lp) {
}
// Initialize row_to_process with the current singleton rows.
RowIndex num_rows(matrix.num_rows());
StrictITIVector<RowIndex, EntryIndex> row_degree(num_rows, EntryIndex(0));
std::vector<RowIndex> row_to_process;
for (RowIndex row(0); row < num_rows; ++row) {
@@ -2347,7 +2360,7 @@ void SingletonPreprocessor::StoreSolution(ProblemSolution* solution) const {
// It is important to indo the operations in the correct order, i.e. in the
// reverse order in which they were done.
for (int i = undo_stack_.size() - 1; i >= 0; --i) {
undo_stack_[i].Undo(solution);
undo_stack_[i].Undo(deleted_columns_, deleted_rows_, solution);
}
}

40
src/glop/preprocessor.h
View File

@@ -21,10 +21,11 @@
#ifndef OR_TOOLS_GLOP_PREPROCESSOR_H_
#define OR_TOOLS_GLOP_PREPROCESSOR_H_
#include "glop/lp_types.h"
#include "glop/matrix_scaler.h"
#include "glop/parameters.pb.h"
#include "glop/revised_simplex.h"
#include "lp_data/lp_types.h"
#include "lp_data/matrix_scaler.h"
namespace operations_research {
namespace glop {
@@ -328,29 +329,34 @@ class SingletonUndo {
MAKE_CONSTRAINT_AN_EQUALITY,
} OperationType;
// Stores all the information needed to undo an operation with the given type.
// All the arguments must refer to the linear program BEFORE the operation is
// applied.
// Stores the information, which together with the field deleted_columns_ and
// deleted_rows_ of SingletonPreprocessor, are needed to undo an operation
// with the given type. Note that all the arguments must refer to the linear
// program BEFORE the operation is applied.
SingletonUndo(OperationType type, const LinearProgram& linear_program,
const SparseColumn& sparse_column_or_row, MatrixEntry e,
ConstraintStatus status);
MatrixEntry e, ConstraintStatus status);
// Undo the operation saved in this class. Note that the operations must be
// undone in the reverse order of the one in which they were applied.
void Undo(ProblemSolution* solution) const;
// Undo the operation saved in this class, taking into account the deleted
// columns and rows passed by the calling instance of SingletonPreprocessor.
// Note that the operations must be undone in the reverse order of the one
// in which they were applied.
void Undo(const SparseMatrix& deleted_columns,
const SparseMatrix& deleted_rows, ProblemSolution* solution) const;
private:
// Actual undo functions for each OperationType.
// Undo() just calls the correct one.
void SingletonRowUndo(ProblemSolution* solution) const;
void ZeroCostSingletonColumnUndo(ProblemSolution* solution) const;
void SingletonColumnInEqualityUndo(ProblemSolution* solution) const;
void SingletonRowUndo(const SparseMatrix& deleted_columns,
ProblemSolution* solution) const;
void ZeroCostSingletonColumnUndo(const SparseMatrix& deleted_rows,
ProblemSolution* solution) const;
void SingletonColumnInEqualityUndo(const SparseMatrix& deleted_rows,
ProblemSolution* solution) const;
void MakeConstraintAnEqualityUndo(ProblemSolution* solution) const;
// All the information needed during undo.
OperationType type_;
bool is_maximization_;
SparseColumn sparse_column_or_row_;
MatrixEntry e_;
Fractional cost_;
@@ -429,6 +435,12 @@ class SingletonPreprocessor : public Preprocessor {
RowDeletionHelper row_deletion_helper_;
std::vector<SingletonUndo> undo_stack_;
// The columns that are deleted by this preprocessor.
SparseMatrix deleted_columns_;
// The transpose of the rows that are deleted by this preprocessor.
// TODO(user): implement a RowMajorSparseMatrix class to simplify the code.
SparseMatrix deleted_rows_;
DISALLOW_COPY_AND_ASSIGN(SingletonPreprocessor);
};

9
src/glop/primal_edge_norms.cc
View File

@@ -18,7 +18,7 @@
#endif
#include "base/timer.h"
#include "glop/lp_utils.h"
#include "lp_data/lp_utils.h"
namespace operations_research {
namespace glop {
@@ -37,7 +37,8 @@ PrimalEdgeNorms::PrimalEdgeNorms(const MatrixView& matrix,
edge_squared_norms_(),
matrix_column_norms_(),
devex_weights_(),
direction_left_inverse_() {}
direction_left_inverse_(),
num_operations_(0) {}
void PrimalEdgeNorms::Clear() {
SCOPED_TIME_STAT(&stats_);
@@ -121,6 +122,7 @@ void PrimalEdgeNorms::ComputeMatrixColumnNorms() {
matrix_column_norms_.resize(matrix_.num_cols(), 0.0);
for (ColIndex col(0); col < matrix_.num_cols(); ++col) {
matrix_column_norms_[col] = sqrt(SquaredNorm(matrix_.column(col)));
num_operations_ += matrix_.column(col).num_entries().value();
}
}
@@ -217,6 +219,7 @@ void PrimalEdgeNorms::UpdateEdgeSquaredNorms(ColIndex entering_col,
const Fractional coeff = update_row.GetCoefficient(col);
const Fractional scalar_product =
compact_matrix_.ColumnScalarProduct(col, direction_left_inverse_);
num_operations_ += compact_matrix_.column(col).num_entries().value();
// Update the edge squared norm of this column. Note that the update
// formula used is important to maximize the precision. See an explanation
@@ -241,6 +244,8 @@ void PrimalEdgeNorms::UpdateEdgeSquaredNorms(ColIndex entering_col,
#ifdef OMP
// In the multi-threaded case, perform the same computation as in the
// single-threaded case above.
//
// TODO(user): also update num_operations_.
std::vector<int> thread_local_stat_lower_bounded_norms(num_omp_threads, 0);
const std::vector<ColIndex>& relevant_rows = update_row.GetNonZeroPositions();
const int parallel_loop_size = relevant_rows.size();

12
src/glop/primal_edge_norms.h
View File

@@ -15,11 +15,11 @@
#define OR_TOOLS_GLOP_PRIMAL_EDGE_NORMS_H_
#include "glop/basis_representation.h"
#include "glop/lp_data.h"
#include "glop/lp_types.h"
#include "glop/parameters.pb.h"
#include "glop/update_row.h"
#include "glop/variables_info.h"
#include "lp_data/lp_data.h"
#include "lp_data/lp_types.h"
#include "util/stats.h"
namespace operations_research {
@@ -114,6 +114,11 @@ class PrimalEdgeNorms {
// Returns a std::string with statistics about this class.
std::string StatString() const { return stats_.StatString(); }
// Deterministic time used by the scalar product computation of this class.
double DeterministicTime() const {
return DeterministicTimeForFpOperations(num_operations_);
}
private:
// Statistics about this class.
struct Stats : public StatsGroup {
@@ -193,6 +198,9 @@ class PrimalEdgeNorms {
DenseRow direction_left_inverse_;
ColIndexVector direction_left_inverse_non_zeros_;
// Used by DeterministicTime().
int64 num_operations_;
DISALLOW_COPY_AND_ASSIGN(PrimalEdgeNorms);
};

2
src/glop/proto_utils.h
View File

@@ -15,7 +15,7 @@
#define OR_TOOLS_GLOP_PROTO_UTILS_H_
#include "linear_solver/linear_solver2.pb.h"
#include "glop/lp_data.h"
#include "lp_data/lp_data.h"
namespace operations_research {
namespace glop {

20
src/glop/rank_one_update.h
View File

@@ -15,10 +15,10 @@
#define OR_TOOLS_GLOP_RANK_ONE_UPDATE_H_
#include "base/logging.h"
#include "glop/lp_types.h"
#include "glop/lp_utils.h"
#include "glop/sparse.h"
#include "glop/status.h"
#include "lp_data/lp_types.h"
#include "lp_data/lp_utils.h"
#include "lp_data/sparse.h"
namespace operations_research {
namespace glop {
@@ -87,6 +87,11 @@ class RankOneUpdateElementaryMatrix {
reinterpret_cast<DenseColumn*>(y));
}
EntryIndex num_entries() const {
return storage_->column(u_index_).num_entries() +
storage_->column(v_index_).num_entries();
}
private:
// This is only used in debug mode.
Fractional ComputeUScalarV() const {
@@ -110,11 +115,15 @@ class RankOneUpdateFactorization {
RankOneUpdateFactorization() : elementary_matrices_() {}
// Deletes all elementary matrices of this factorization.
void Clear() { elementary_matrices_.clear(); }
void Clear() {
elementary_matrices_.clear();
num_entries_ = 0;
}
// Updates the factorization.
void Update(const RankOneUpdateElementaryMatrix& update_matrix) {
elementary_matrices_.push_back(update_matrix);
num_entries_ += update_matrix.num_entries();
}
// Left-solves all systems from right to left, i.e. y_i = y_{i+1}.(T_i)^{-1}
@@ -134,7 +143,10 @@ class RankOneUpdateFactorization {
}
}
EntryIndex num_entries() const { return num_entries_; }
private:
EntryIndex num_entries_;
std::vector<RankOneUpdateElementaryMatrix> elementary_matrices_;
DISALLOW_COPY_AND_ASSIGN(RankOneUpdateFactorization);
};

2
src/glop/reduced_costs.cc
View File

@@ -17,7 +17,7 @@
#include <omp.h>
#endif
#include "glop/lp_utils.h"
#include "lp_data/lp_utils.h"
namespace operations_research {
namespace glop {

4
src/glop/reduced_costs.h
View File

@@ -15,13 +15,13 @@
#define OR_TOOLS_GLOP_REDUCED_COSTS_H_
#include "glop/basis_representation.h"
#include "glop/lp_data.h"
#include "glop/lp_types.h"
#include "glop/parameters.pb.h"
#include "glop/primal_edge_norms.h"
#include "glop/status.h"
#include "glop/update_row.h"
#include "glop/variables_info.h"
#include "lp_data/lp_data.h"
#include "lp_data/lp_types.h"
#include "util/stats.h"
namespace operations_research {

19
src/glop/revised_simplex.cc
View File

@@ -27,11 +27,11 @@
#include "base/stringprintf.h"
#include "base/timer.h"
#include "glop/initial_basis.h"
#include "glop/lp_data.h"
#include "glop/lp_print_utils.h"
#include "glop/lp_utils.h"
#include "glop/matrix_utils.h"
#include "glop/parameters.pb.h"
#include "lp_data/lp_data.h"
#include "lp_data/lp_print_utils.h"
#include "lp_data/lp_utils.h"
#include "lp_data/matrix_utils.h"
#include "util/fp_utils.h"
DEFINE_bool(simplex_display_numbers_as_fractions, false,
@@ -338,6 +338,14 @@ std::string RevisedSimplex::GetPrettySolverStats() const {
FLAGS_simplex_stop_after_first_basis);
}
double RevisedSimplex::DeterministicTime() const {
// TODO(user): Also take into account the dual edge norms and the reduced cost
// updates.
return basis_factorization_.DeterministicTime() +
update_row_.DeterministicTime() +
primal_edge_norms_.DeterministicTime();
}
void RevisedSimplex::SetVariableNames() {
variable_name_.resize(num_cols_, "");
for (ColIndex col(0); col < first_slack_col_; ++col) {
@@ -2106,7 +2114,8 @@ Status RevisedSimplex::Minimize(TimeLimit* time_limit) {
// ProblemStatus::PRIMAL_FEASIBLE if it is the case at the beginning of the
// algorithm.
if (num_iterations_ == parameters_.max_number_of_iterations() ||
time_limit->LimitReached()) {
time_limit->LimitReached() ||
DeterministicTime() > parameters_.max_deterministic_time()) {
break;
}

13
src/glop/revised_simplex.h
View File

@@ -102,11 +102,6 @@
#include "base/macros.h"
#include "glop/basis_representation.h"
#include "glop/dual_edge_norms.h"
#include "glop/lp_data.h"
#include "glop/lp_print_utils.h"
#include "glop/lp_types.h"
#include "glop/matrix_scaler.h"
#include "glop/mps_reader.h"
#include "glop/parameters.pb.h"
#include "glop/primal_edge_norms.h"
#include "glop/entering_variable.h"
@@ -115,6 +110,10 @@
#include "glop/update_row.h"
#include "glop/variables_info.h"
#include "glop/variable_values.h"
#include "lp_data/lp_data.h"
#include "lp_data/lp_print_utils.h"
#include "lp_data/lp_types.h"
#include "lp_data/matrix_scaler.h"
#include "base/random.h"
#include "util/time_limit.h"
@@ -190,6 +189,7 @@ class RevisedSimplex {
VariableStatus GetVariableStatus(ColIndex col) const;
ConstraintStatus GetConstraintStatus(RowIndex row) const;
const BasisState& GetState() const;
double DeterministicTime() const;
// If the problem status is PRIMAL_UNBOUNDED (respectively DUAL_UNBOUNDED),
// then the solver has a corresponding primal (respectively dual) ray to show
@@ -531,7 +531,7 @@ class RevisedSimplex {
ColIndex first_slack_col_;
// We're using vectors after profiling and looking at the generated assembly
// it's as fast as std::unique_ptr as long as the size is properly reserved
// it's as fast as scoped_ptr as long as the size is properly reserved
// beforehand.
// Temporary view of the matrix given to Solve() with extra slack columns.
@@ -730,6 +730,7 @@ class RevisedSimplex {
// A random number generator.
MTRandom random_;
DISALLOW_COPY_AND_ASSIGN(RevisedSimplex);
};

15
src/glop/status.h
View File

@@ -87,21 +87,6 @@ std::string GetErrorCodeString(Status::ErrorCode error_code);
return Status(Status::ERROR_NULL, error_message); \
}
// Macros to replace CHECK_NOTNULL() so we don't crash in production.
// Logs a FATAL message in debug mode, and an ERROR message in production.
// It is not perfect, but more robust than crashing right away.
#define RETURN_IF_NULL(x) \
if (x == nullptr) { \
LOG(DFATAL) << #x << " == NULL"; \
return; \
}
#define RETURN_VALUE_IF_NULL(x, v) \
if (x == nullptr) { \
LOG(DFATAL) << #x << " == NULL"; \
return v; \
}
} // namespace glop
} // namespace operations_research

11
src/glop/update_row.cc
View File

@@ -17,7 +17,7 @@
#include <omp.h>
#endif
#include "glop/lp_utils.h"
#include "lp_data/lp_utils.h"
namespace operations_research {
namespace glop {
@@ -37,6 +37,7 @@ UpdateRow::UpdateRow(const CompactSparseMatrix& matrix,
non_zero_position_set_(),
coefficient_(),
compute_update_row_(true),
num_operations_(0),
parameters_(),
stats_() {}
@@ -96,14 +97,22 @@ void UpdateRow::ComputeUpdateRow(RowIndex leaving_row) {
if (lhs < 0.5 * static_cast<double>(num_col_wise_entries.value())) {
if (lhs < 1.1 * static_cast<double>(matrix_.num_cols().value())) {
ComputeUpdatesRowWiseHypersparse();
num_operations_ += num_row_wise_entries.value();
} else {
ComputeUpdatesRowWise();
num_operations_ +=
num_row_wise_entries.value() + matrix_.num_rows().value();
}
} else {
ComputeUpdatesColumnWise();
num_operations_ +=
num_col_wise_entries.value() + matrix_.num_cols().value();
}
} else {
ComputeUpdatesColumnWise();
num_operations_ +=
variables_info_.GetNumEntriesInRelevantColumns().value() +
matrix_.num_cols().value();
}
IF_STATS_ENABLED(stats_.update_row_density.Add(
static_cast<double>(non_zero_position_list_.size()) /

11
src/glop/update_row.h
View File

@@ -15,9 +15,9 @@
#define OR_TOOLS_GLOP_UPDATE_ROW_H_
#include "glop/basis_representation.h"
#include "glop/lp_types.h"
#include "glop/parameters.pb.h"
#include "glop/variables_info.h"
#include "lp_data/lp_types.h"
#include "util/stats.h"
namespace operations_research {
@@ -83,6 +83,11 @@ class UpdateRow {
void ComputeUpdateRowForBenchmark(const DenseRow& lhs,
const std::string& algorithm);
// Deterministic time used by the scalar product computation of this class.
double DeterministicTime() const {
return DeterministicTimeForFpOperations(num_operations_);
}
private:
// Computes the left inverse of the given unit row, and stores it in
// unit_row_left_inverse_ and unit_row_left_inverse_non_zeros_.
@@ -129,6 +134,10 @@ class UpdateRow {
RatioDistribution update_row_density;
};
// Track the number of basic floating point multiplication.
// Used by DeterministicTime().
int64 num_operations_;
// Glop standard classes.
GlopParameters parameters_;
Stats stats_;

2
src/glop/variable_values.cc
View File

@@ -13,7 +13,7 @@
#include "glop/variable_values.h"
#include "glop/lp_utils.h"
#include "lp_data/lp_utils.h"
#include "util/iterators.h"
namespace operations_research {

2
src/glop/variable_values.h
View File

@@ -15,8 +15,8 @@
#define OR_TOOLS_GLOP_VARIABLE_VALUES_H_
#include "glop/basis_representation.h"
#include "glop/lp_types.h"
#include "glop/variables_info.h"
#include "lp_data/lp_types.h"
#include "util/stats.h"
namespace operations_research {

4
src/glop/variables_info.h
View File

@@ -14,8 +14,8 @@
#ifndef OR_TOOLS_GLOP_VARIABLES_INFO_H_
#define OR_TOOLS_GLOP_VARIABLES_INFO_H_
#include "glop/lp_types.h"
#include "glop/sparse.h"
#include "lp_data/lp_types.h"
#include "lp_data/sparse.h"
namespace operations_research {
namespace glop {

6
src/linear_solver/glop_interface.cc
View File

@@ -25,11 +25,11 @@
#include "base/file.h"
#include "google/protobuf/text_format.h"
#include "base/hash.h"
#include "glop/lp_data.h"
#include "glop/lp_solver.h"
#include "glop/lp_types.h"
#include "glop/parameters.pb.h"
#include "linear_solver/linear_solver.h"
#include "lp_data/lp_data.h"
#include "lp_data/lp_types.h"
DECLARE_double(solver_timeout_in_seconds);
DECLARE_string(solver_write_model);
@@ -189,7 +189,7 @@ MPSolver::ResultStatus GLOPInterface::Solve(const MPSolverParameters& param) {
if (solver_->time_limit()) {
VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms.";
parameters_.set_max_time_in_seconds(
1000.0 * static_cast<double>(solver_->time_limit()));
static_cast<double>(solver_->time_limit()) / 1000.0);
}
solver_->SetSolverSpecificParametersAsString(

5
src/linear_solver/gurobi_interface.cc
View File

@@ -36,6 +36,8 @@ extern "C" {
#define CHECKED_GUROBI_CALL(x) CHECK_EQ(0, x)
DEFINE_int32(num_gurobi_threads, 4, "Number of threads available for Gurobi.");
namespace operations_research {
class GurobiInterface : public MPSolverInterface {
@@ -174,6 +176,9 @@ GurobiInterface::GurobiInterface(MPSolver* const solver, bool mip)
NULL)); // varnanes
CHECKED_GUROBI_CALL(
GRBsetintattr(model_, GRB_INT_ATTR_MODELSENSE, maximize_ ? -1 : 1));
CHECKED_GUROBI_CALL(
GRBsetintparam(env_, GRB_INT_PAR_THREADS, FLAGS_num_gurobi_threads));
}
GurobiInterface::~GurobiInterface() {

19
src/linear_solver/linear_solver.cc
View File

@@ -503,6 +503,7 @@ MPSolver::LoadStatus MPSolver::LoadModelFromProto(
const new_proto::MPConstraintProto& ct_proto = input_model.constraint(i);
MPConstraint* const ct = MakeRowConstraint(
ct_proto.lower_bound(), ct_proto.upper_bound(), ct_proto.name());
ct->set_is_lazy(ct_proto.is_lazy());
if (ct_proto.var_index_size() != ct_proto.coefficient_size()) {
LOG(ERROR) << "In constraint #" << i << " (name: '" << ct_proto.name()
<< "'):"
@@ -562,6 +563,10 @@ void MPSolver::FillSolutionResponseProto(
for (int i = 0; i < variables_.size(); ++i) {
response->add_variable_value(variables_[i]->solution_value());
}
if (interface_->IsMIP()) {
response->set_best_objective_bound(interface_->best_objective_bound());
}
}
}
@@ -727,6 +732,8 @@ void MPSolver::EnableOutput() { interface_->set_quiet(false); }
void MPSolver::SuppressOutput() { interface_->set_quiet(true); }
bool MPSolver::InterruptSolve() { return interface_->InterruptSolve(); }
MPVariable* MPSolver::MakeVar(double lb, double ub, bool integer,
const std::string& name) {
const int var_index = NumVariables();
@@ -1219,8 +1226,11 @@ double MPSolverInterface::ComputeExactConditionNumber() const {
}
void MPSolverInterface::SetCommonParameters(const MPSolverParameters& param) {
// TODO(user): Overhaul the code that sets parameters to enable changing
// GLOP parameters without issuing warnings.
// By default, we let GLOP keep its own default tolerance, much more accurate
// than for the rest of the solvers.
//
#if defined(USE_GLOP)
if (solver_->ProblemType() != MPSolver::GLOP_LINEAR_PROGRAMMING) {
#endif
@@ -1241,7 +1251,14 @@ void MPSolverInterface::SetCommonParameters(const MPSolverParameters& param) {
}
void MPSolverInterface::SetMIPParameters(const MPSolverParameters& param) {
SetRelativeMipGap(param.GetDoubleParam(MPSolverParameters::RELATIVE_MIP_GAP));
#if defined(USE_GLOP)
if (solver_->ProblemType() != MPSolver::GLOP_LINEAR_PROGRAMMING) {
#endif
SetRelativeMipGap(param.GetDoubleParam(
MPSolverParameters::RELATIVE_MIP_GAP));
#if defined(USE_GLOP)
}
#endif
}
void MPSolverInterface::SetUnsupportedDoubleParam(

9
src/linear_solver/linear_solver.h
View File

@@ -341,6 +341,13 @@ class MPSolver {
// everything was reconstructed from scratch.
void Reset();
// Interrupts the Solve() execution to terminate processing early if possible.
// If the underlying interface supports interruption; it does that and returns
// true regardless of whether there's an ongoing Solve() or not. The Solve()
// call may still linger for a while depending on the conditions. If
// interruption is not supported; returns false and does nothing.
bool InterruptSolve();
// ----- Methods using protocol buffers -----
// The status of loading the problem from a protocol buffer.
@@ -1182,6 +1189,8 @@ class MPSolverInterface {
// problems and only implemented in GLPK.
virtual double ComputeExactConditionNumber() const;
virtual bool InterruptSolve() { return false; }
friend class MPSolver;
// To access the maximize_ bool and the MPSolver.

14
src/linear_solver/linear_solver2.proto
View File

@@ -132,15 +132,16 @@ message MPModelRequest {
// program?).
// This must remain consistent with MPSolver::OptimizationProblemType.
enum SolverType {
GLOP_LINEAR_PROGRAMMING = 2; // Recommended default for LP models.
GLOP_LINEAR_PROGRAMMING = 2; // Recommended default for LP models.
CLP_LINEAR_PROGRAMMING = 0;
GLPK_LINEAR_PROGRAMMING = 1;
GUROBI_LINEAR_PROGRAMMING = 6; // Commercial, needs a valid license.
GUROBI_LINEAR_PROGRAMMING = 6; // Commercial, needs a valid license.
SCIP_MIXED_INTEGER_PROGRAMMING = 3; // Recommended default for MIP models.
SCIP_MIXED_INTEGER_PROGRAMMING = 3; // Recommended default for MIP models.
GLPK_MIXED_INTEGER_PROGRAMMING = 4;
CBC_MIXED_INTEGER_PROGRAMMING = 5;
GUROBI_MIXED_INTEGER_PROGRAMMING = 7; // Commercial, needs a valid license.
}
optional SolverType solver_type = 2;
@@ -186,6 +187,13 @@ message MPSolutionResponse {
// This is set iff 'status' is OPTIMAL or FEASIBLE.
optional double objective_value = 2;
// This field is only filled for MIP problems. For a minimization problem,
// this is a lower bound on the optimal objective value. For a maximization
// problem, it is an upper bound. It is only filled if the status is OPTIMAL
// or FEASIBLE. In the former case, best_objective_bound should be equal to
// objective_value (modulo numerical errors).
optional double best_objective_bound = 5;
// Variable values in the same order as the MPModelProto::variable field.
// This is a dense representation. These are set iff 'status' is OPTIMAL or
// FEASIBLE.

27
src/linear_solver/model_exporter.cc
View File

@@ -34,6 +34,9 @@ DEFINE_int32(lp_max_line_length, 10000,
"Maximum line length in exported .lp files. The default was chosen"
" so that SCIP can read the files.");
DEFINE_bool(lp_log_invalid_name, false,
"Whether to log invalid variable and contraint names.");
namespace operations_research {
using new_proto::MPConstraintProto;
@@ -56,30 +59,32 @@ MPModelProtoExporter::MPModelProtoExporter(const MPModelProto& proto)
// Note(user): This method is static. It is also used by MPSolver.
bool MPModelProtoExporter::CheckNameValidity(const std::string& name) {
if (name.empty()) {
LOG(WARNING) << "CheckNameValidity() should not be passed an empty name.";
LOG_IF(WARNING, FLAGS_lp_log_invalid_name)
<< "CheckNameValidity() should not be passed an empty name.";
return false;
}
// Allow names that conform to the LP and MPS format.
const int kMaxNameLength = 255;
if (name.size() > kMaxNameLength) {
LOG(WARNING) << "Invalid name " << name << ": length > " << kMaxNameLength
<< "."
<< " Will be unable to write model to file.";
LOG_IF(WARNING, FLAGS_lp_log_invalid_name)
<< "Invalid name " << name << ": length > " << kMaxNameLength << "."
<< " Will be unable to write model to file.";
return false;
}
const std::string kForbiddenChars = " +-*/<>=:\\";
if (name.find_first_of(kForbiddenChars) != std::string::npos) {
LOG(WARNING) << "Invalid name " << name
<< " contains forbidden character: " << kForbiddenChars
<< " or space."
<< " Will be unable to write model to file.";
LOG_IF(WARNING, FLAGS_lp_log_invalid_name)
<< "Invalid name " << name
<< " contains forbidden character: " << kForbiddenChars << " or space."
<< " Will be unable to write model to file.";
return false;
}
const std::string kForbiddenFirstChars = "$.0123456789";
if (kForbiddenFirstChars.find(name[0]) != std::string::npos) {
LOG(WARNING) << "Invalid name " << name
<< ". First character is one of: " << kForbiddenFirstChars
<< " Will be unable to write model to file.";
LOG_IF(WARNING, FLAGS_lp_log_invalid_name)
<< "Invalid name " << name
<< ". First character is one of: " << kForbiddenFirstChars
<< " Will be unable to write model to file.";
return false;
}
return true;

13
src/linear_solver/scip_interface.cc
View File

@@ -121,6 +121,11 @@ class SCIPInterface : public MPSolverInterface {
SCIPlpiGetSolverName());
}
bool InterruptSolve() override {
if (scip_ != NULL) SCIPinterruptSolve(scip_);
return true;
}
virtual void* underlying_solver() { return reinterpret_cast<void*>(scip_); }
private:
@@ -163,6 +168,10 @@ void SCIPInterface::Reset() {
void SCIPInterface::CreateSCIP() {
ORTOOLS_SCIP_CALL(SCIPcreate(&scip_));
ORTOOLS_SCIP_CALL(SCIPincludeDefaultPlugins(scip_));
// Default clock type (1: CPU user seconds, 2: wall clock time). We use
// wall clock time because getting CPU user seconds involves calling
// times() which is very expensive.
ORTOOLS_SCIP_CALL(SCIPsetIntParam(scip_, "timing/clocktype", 2));
ORTOOLS_SCIP_CALL(SCIPcreateProb(scip_, solver_->name_.c_str(), NULL, NULL,
NULL, NULL, NULL, NULL, NULL));
ORTOOLS_SCIP_CALL(SCIPsetObjsense(
@@ -481,7 +490,9 @@ MPSolver::ResultStatus SCIPInterface::Solve(const MPSolverParameters& param) {
// Solve.
timer.Restart();
ORTOOLS_SCIP_CALL(SCIPsolve(scip_));
if (SCIPsolve(scip_) != SCIP_OKAY) {
return MPSolver::ABNORMAL;
}
VLOG(1) << StringPrintf("Solved in %.3f seconds.", timer.Get());
// Get the results.

7
src/glop/lp_data.cc → src/lp_data/lp_data.cc
View File

@@ -12,7 +12,7 @@
// limitations under the License.
#include "glop/lp_data.h"
#include "lp_data/lp_data.h"
#include <algorithm>
#include <cmath>
@@ -23,9 +23,8 @@
#include "base/logging.h"
#include "base/stringprintf.h"
#include "base/join.h"
#include "glop/lp_print_utils.h"
#include "glop/lp_utils.h"
#include "glop/status.h"
#include "lp_data/lp_print_utils.h"
#include "lp_data/lp_utils.h"
namespace operations_research {
namespace glop {

12
src/glop/lp_data.h → src/lp_data/lp_data.h
View File

@@ -20,8 +20,8 @@
// - bounds for each variable,
// - bounds for each constraint.
#ifndef OR_TOOLS_GLOP_LP_DATA_H_
#define OR_TOOLS_GLOP_LP_DATA_H_
#ifndef OR_TOOLS_LP_DATA_LP_DATA_H_
#define OR_TOOLS_LP_DATA_LP_DATA_H_
#include <algorithm> // for max
#include <map>
@@ -34,9 +34,9 @@
#include "base/int_type.h"
#include "base/int_type_indexed_vector.h"
#include "base/hash.h"
#include "glop/lp_types.h"
#include "glop/matrix_scaler.h"
#include "glop/sparse.h"
#include "lp_data/lp_types.h"
#include "lp_data/matrix_scaler.h"
#include "lp_data/sparse.h"
#include "util/fp_utils.h"
namespace operations_research {
@@ -451,4 +451,4 @@ inline bool AreBoundsValid(Fractional lower_bound, Fractional upper_bound) {
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_GLOP_LP_DATA_H_
#endif // OR_TOOLS_LP_DATA_LP_DATA_H_

4
src/glop/lp_print_utils.cc → src/lp_data/lp_print_utils.cc
View File

@@ -19,8 +19,8 @@
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/join.h"
#include "glop/lp_print_utils.h"
#include "glop/lp_types.h"
#include "lp_data/lp_print_utils.h"
#include "lp_data/lp_types.h"
#include "util/rational_approximation.h"
namespace operations_research {

8
src/glop/lp_print_utils.h → src/lp_data/lp_print_utils.h
View File

@@ -14,14 +14,14 @@
// Utilities to display linear expression in a human-readable way.
#ifndef OR_TOOLS_GLOP_LP_PRINT_UTILS_H_
#define OR_TOOLS_GLOP_LP_PRINT_UTILS_H_
#ifndef OR_TOOLS_LP_DATA_LP_PRINT_UTILS_H_
#define OR_TOOLS_LP_DATA_LP_PRINT_UTILS_H_
#include <string>
#include "base/integral_types.h"
#include "base/stringprintf.h"
#include "glop/lp_types.h"
#include "lp_data/lp_types.h"
namespace operations_research {
namespace glop {
@@ -56,4 +56,4 @@ std::string StringifyMonomial(const Fractional a, const std::string& x, bool fra
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_GLOP_LP_PRINT_UTILS_H_
#endif // OR_TOOLS_LP_DATA_LP_PRINT_UTILS_H_

2
src/glop/lp_types.cc → src/lp_data/lp_types.cc
View File

@@ -11,7 +11,7 @@
// See the License for the specific language governing permissions and
// limitations under the License.
#include "glop/lp_types.h"
#include "lp_data/lp_types.h"
namespace operations_research {
namespace glop {

14
src/glop/lp_types.h → src/lp_data/lp_types.h
View File

@@ -14,8 +14,8 @@
// Common types and constants used by the Linear Programming solver.
#ifndef OR_TOOLS_GLOP_LP_TYPES_H_
#define OR_TOOLS_GLOP_LP_TYPES_H_
#ifndef OR_TOOLS_LP_DATA_LP_TYPES_H_
#define OR_TOOLS_LP_DATA_LP_TYPES_H_
#include <math.h>
#include <limits>
@@ -367,7 +367,15 @@ struct ScatteredColumnReference {
static const double kDenseThresholdForPreciseSum;
};
// This is used during the deterministic time computation to convert a given
// number of floating-point operations to something in the same order of
// magnitude as a second (on a 2014 desktop).
static inline double DeterministicTimeForFpOperations(int64 n) {
const double kConvertionFactor = 2e-9;
return kConvertionFactor * static_cast<double>(n);
}
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_GLOP_LP_TYPES_H_
#endif // OR_TOOLS_LP_DATA_LP_TYPES_H_

2
src/glop/lp_utils.cc → src/lp_data/lp_utils.cc
View File

@@ -12,7 +12,7 @@
// limitations under the License.
#include "glop/lp_utils.h"
#include "lp_data/lp_utils.h"
namespace operations_research {
namespace glop {

10
src/glop/lp_utils.h → src/lp_data/lp_utils.h
View File

@@ -14,12 +14,12 @@
// Basic utility functions on Fractional or row/column of Fractional.
#ifndef OR_TOOLS_GLOP_LP_UTILS_H_
#define OR_TOOLS_GLOP_LP_UTILS_H_
#ifndef OR_TOOLS_LP_DATA_LP_UTILS_H_
#define OR_TOOLS_LP_DATA_LP_UTILS_H_
#include "base/accurate_sum.h"
#include "glop/lp_types.h"
#include "glop/sparse_column.h"
#include "lp_data/lp_types.h"
#include "lp_data/sparse_column.h"
namespace operations_research {
namespace glop {
@@ -315,4 +315,4 @@ typedef SumWithOneMissing<false> SumWithNegativeInfiniteAndOneMissing;
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_GLOP_LP_UTILS_H_
#endif // OR_TOOLS_LP_DATA_LP_UTILS_H_

7
src/glop/matrix_scaler.cc → src/lp_data/matrix_scaler.cc
View File

@@ -12,7 +12,7 @@
// limitations under the License.
#include "glop/matrix_scaler.h"
#include "lp_data/matrix_scaler.h"
#include <algorithm>
#include <cmath>
@@ -20,9 +20,8 @@
#include "base/logging.h"
#include "base/stringprintf.h"
#include "glop/lp_utils.h"
#include "glop/sparse.h"
#include "glop/status.h"
#include "lp_data/lp_utils.h"
#include "lp_data/sparse.h"
namespace operations_research {
namespace glop {

8
src/glop/matrix_scaler.h → src/lp_data/matrix_scaler.h
View File

@@ -58,15 +58,15 @@
// and:
// A'.x' = R.A.C.C^-1.x = R.A.x = R.b = b'.
#ifndef OR_TOOLS_GLOP_MATRIX_SCALER_H_
#define OR_TOOLS_GLOP_MATRIX_SCALER_H_
#ifndef OR_TOOLS_LP_DATA_MATRIX_SCALER_H_
#define OR_TOOLS_LP_DATA_MATRIX_SCALER_H_
#include <vector>
#include "base/integral_types.h"
#include "base/macros.h"
#include "base/int_type_indexed_vector.h"
#include "glop/lp_types.h"
#include "lp_data/lp_types.h"
namespace operations_research {
namespace glop {
@@ -159,4 +159,4 @@ class SparseMatrixScaler {
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_GLOP_MATRIX_SCALER_H_
#endif // OR_TOOLS_LP_DATA_MATRIX_SCALER_H_

220
src/lp_data/matrix_utils.cc Normal file
View File

@@ -0,0 +1,220 @@
// Copyright 2010-2014 Google
// 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 "lp_data/matrix_utils.h"
#include <algorithm>
#include "base/hash.h"
namespace operations_research {
namespace glop {
namespace {
// Returns true iff the two given sparse columns are proportional. The two
// sparse columns must be ordered by row and must not contain any zero entry.
//
// See the header comment on FindProportionalColumns() for the exact definition
// of two proportional columns with a given tolerance.
bool AreColumnsProportional(const SparseColumn& a, const SparseColumn& b,
Fractional tolerance) {
DCHECK(a.IsCleanedUp());
DCHECK(b.IsCleanedUp());
if (a.num_entries() != b.num_entries()) return false;
Fractional multiple = 0.0;
bool a_is_larger = true;
for (const EntryIndex i : a.AllEntryIndices()) {
if (a.EntryRow(i) != b.EntryRow(i)) return false;
const Fractional coeff_a = a.EntryCoefficient(i);
const Fractional coeff_b = b.EntryCoefficient(i);
if (multiple == 0.0) {
a_is_larger = fabs(coeff_a) > fabs(coeff_b);
multiple = a_is_larger ? coeff_a / coeff_b : coeff_b / coeff_a;
} else {
if (a_is_larger) {
if (fabs(coeff_a / coeff_b - multiple) > tolerance) return false;
} else {
if (fabs(coeff_b / coeff_a - multiple) > tolerance) return false;
}
}
}
return true;
}
// A column index together with its fingerprint. See ComputeFingerprint().
struct ColumnFingerprint {
ColumnFingerprint(ColIndex _col, int32 _hash, double _value)
: col(_col), hash(_hash), value(_value) {}
ColIndex col;
int32 hash;
double value;
// This order has the property that if AreProportionalCandidates() is true for
// two given columns, then in a sorted list of columns
// AreProportionalCandidates() will be true for all the pairs of columns
// between the two given ones (included).
bool operator<(const ColumnFingerprint& other) const {
if (hash == other.hash) {
return value < other.value;
}
return hash < other.hash;
}
};
// Two columns can be proportional only if:
// - Their non-zero pattern hashes are the same.
// - Their double fingerprints are close to each other.
bool AreProportionalCandidates(ColumnFingerprint a, ColumnFingerprint b,
Fractional tolerance) {
if (a.hash != b.hash) return false;
return fabs(a.value - b.value) < tolerance;
}
// The fingerprint of a column has two parts:
// - A hash value of the column non-zero pattern.
// - A double value which should be the same for two proportional columns
// modulo numerical errors.
ColumnFingerprint ComputeFingerprint(ColIndex col, const SparseColumn& column) {
int32 non_zero_pattern_hash = 0;
Fractional min_abs = std::numeric_limits<Fractional>::max();
Fractional max_abs = 0.0;
Fractional sum = 0.0;
for (const SparseColumn::Entry e : column) {
non_zero_pattern_hash =
Hash32NumWithSeed(e.row().value(), non_zero_pattern_hash);
sum += e.coefficient();
min_abs = std::min(min_abs, fabs(e.coefficient()));
max_abs = std::max(max_abs, fabs(e.coefficient()));
}
// The two scaled values are in [0, 1].
// TODO(user): A better way to discriminate columns would be to take the
// scalar product with a constant but random vector scaled by max_abs.
DCHECK_NE(0.0, max_abs);
const double inverse_dynamic_range = min_abs / max_abs;
const double scaled_average =
fabs(sum) / (static_cast<double>(column.num_entries().value()) * max_abs);
return ColumnFingerprint(col, non_zero_pattern_hash,
inverse_dynamic_range + scaled_average);
}
} // namespace
ColMapping FindProportionalColumns(const SparseMatrix& matrix,
Fractional tolerance) {
const ColIndex num_cols = matrix.num_cols();
ColMapping mapping(num_cols, kInvalidCol);
// Compute the fingerprint of each columns and sort them.
std::vector<ColumnFingerprint> fingerprints;
for (ColIndex col(0); col < num_cols; ++col) {
if (!matrix.column(col).IsEmpty()) {
fingerprints.push_back(ComputeFingerprint(col, matrix.column(col)));
}
}
std::sort(fingerprints.begin(), fingerprints.end());
// Find a representative of each proportional columns class. This only
// compares columns with a close-enough fingerprint.
for (int i = 0; i < fingerprints.size(); ++i) {
const ColIndex col_a = fingerprints[i].col;
if (mapping[col_a] != kInvalidCol) continue;
for (int j = i + 1; j < fingerprints.size(); ++j) {
const ColIndex col_b = fingerprints[j].col;
if (mapping[col_b] != kInvalidCol) continue;
// Note that we use the same tolerance for the fingerprints.
// TODO(user): Derive precise bounds on what this tolerance should be so
// that no proportional columns are missed.
if (!AreProportionalCandidates(fingerprints[i], fingerprints[j],
tolerance))
break;
if (AreColumnsProportional(matrix.column(col_a), matrix.column(col_b),
tolerance)) {
mapping[col_b] = col_a;
}
}
}
// Sort the mapping so that the representative of each class is the smallest
// column. To achieve this, the current representative is used as a pointer
// to the new one, a bit like in an union find algorithm.
for (ColIndex col(0); col < num_cols; ++col) {
if (mapping[col] == kInvalidCol) continue;
const ColIndex new_representative = mapping[mapping[col]];
if (new_representative != kInvalidCol) {
mapping[col] = new_representative;
} else {
if (mapping[col] > col) {
mapping[mapping[col]] = col;
mapping[col] = kInvalidCol;
}
}
}
return mapping;
}
ColMapping FindProportionalColumnsUsingSimpleAlgorithm(
const SparseMatrix& matrix, Fractional tolerance) {
const ColIndex num_cols = matrix.num_cols();
ColMapping mapping(num_cols, kInvalidCol);
for (ColIndex col_a(0); col_a < num_cols; ++col_a) {
if (mapping[col_a] != kInvalidCol) continue;
for (ColIndex col_b(col_a + 1); col_b < num_cols; ++col_b) {
if (mapping[col_b] != kInvalidCol) continue;
if (AreColumnsProportional(matrix.column(col_a), matrix.column(col_b),
tolerance)) {
mapping[col_b] = col_a;
}
}
}
return mapping;
}
bool AreFirstColumnsAndRowsExactlyEquals(RowIndex num_rows, ColIndex num_cols,
const SparseMatrix& matrix_a,
const CompactSparseMatrix& matrix_b) {
// TODO(user): Also DCHECK() that matrix_b is ordered by rows.
DCHECK(matrix_a.IsCleanedUp());
if (num_rows > matrix_a.num_rows() || num_rows > matrix_b.num_rows() ||
num_cols > matrix_a.num_cols() || num_cols > matrix_b.num_cols()) {
return false;
}
for (ColIndex col(0); col < num_cols; ++col) {
const SparseColumn& col_a = matrix_a.column(col);
const CompactSparseMatrix::ColumnView& col_b = matrix_b.column(col);
const EntryIndex end = std::min(col_a.num_entries(), col_b.num_entries());
for (EntryIndex i(0); i < end; ++i) {
if (col_a.EntryRow(i) != col_b.EntryRow(i)) {
if (col_a.EntryRow(i) < num_rows || col_b.EntryRow(i) < num_rows) {
return false;
} else {
break;
}
}
if (col_a.EntryCoefficient(i) != col_b.EntryCoefficient(i)) {
return false;
}
if (col_a.num_entries() > end && col_a.EntryRow(end) < num_rows) {
return false;
}
if (col_b.num_entries() > end && col_b.EntryRow(end) < num_rows) {
return false;
}
}
}
return true;
}
} // namespace glop
} // namespace operations_research

55
src/lp_data/matrix_utils.h Normal file
View File

@@ -0,0 +1,55 @@
// Copyright 2010-2014 Google
// 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.
#ifndef OR_TOOLS_LP_DATA_MATRIX_UTILS_H_
#define OR_TOOLS_LP_DATA_MATRIX_UTILS_H_
#include "lp_data/lp_types.h"
#include "lp_data/sparse.h"
namespace operations_research {
namespace glop {
// Finds the proportional columns of the given matrix: all the pairs of columns
// such that one is a non-zero scalar multiple of the other. The returned
// mapping for a given column will either be:
// - The index of the first column which is proportional with this one.
// - Or kInvalidCol if this column is not proportional to any other.
//
// Because of precision errors, two columns are declared proportional if the
// multiples (>= 1.0) defined by all the couple of coefficients of the same row
// are equal up to the given absolute tolerance.
//
// The complexity is in most cases O(num entries of the matrix). However,
// compared to the less efficient algorithm below, it is highly unlikely but
// possible that some pairs of proportional columns are not detected.
ColMapping FindProportionalColumns(const SparseMatrix& matrix,
Fractional tolerance);
// A simple version of FindProportionalColumns() that compares all the columns
// pairs one by one. This is slow, but here for reference. The complexity is
// O(num_cols * num_entries).
ColMapping FindProportionalColumnsUsingSimpleAlgorithm(
const SparseMatrix& matrix, Fractional tolerance);
// Returns true iff the two given matrices have exactly the same first num_rows
// entries on the first num_cols columns. The two given matrices must be ordered
// by rows (this is DCHECKed, but only for the first one at this point).
bool AreFirstColumnsAndRowsExactlyEquals(RowIndex num_rows, ColIndex num_cols,
const SparseMatrix& matrix_a,
const CompactSparseMatrix& matrix_b);
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_LP_DATA_MATRIX_UTILS_H_

5
src/glop/mps_reader.cc → src/lp_data/mps_reader.cc
View File

@@ -12,7 +12,7 @@
// limitations under the License.
#include "glop/mps_reader.h"
#include "lp_data/mps_reader.h"
#include <math.h>
#include <algorithm>
@@ -30,7 +30,8 @@
#include "base/numbers.h" // for safe_strtod
#include "base/split.h"
#include "base/strutil.h"
#include "glop/lp_print_utils.h"
#include "lp_data/lp_print_utils.h"
#include "base/status.h"
DEFINE_bool(mps_free_form, false, "Read MPS files in free form.");
DEFINE_bool(mps_stop_after_first_error, true, "Stop after the first error.");

10
src/glop/mps_reader.h → src/lp_data/mps_reader.h
View File

@@ -20,8 +20,8 @@
// standard. We have developed this reader to be able to read benchmark data
// files. Using the MPS file format for new models is discouraged.
#ifndef OR_TOOLS_GLOP_MPS_READER_H_
#define OR_TOOLS_GLOP_MPS_READER_H_
#ifndef OR_TOOLS_LP_DATA_MPS_READER_H_
#define OR_TOOLS_LP_DATA_MPS_READER_H_
#include <algorithm> // for max
#include <map>
@@ -35,8 +35,8 @@
#include "base/int_type_indexed_vector.h"
#include "base/map_util.h" // for FindOrNull, FindWithDefault
#include "base/hash.h"
#include "glop/lp_data.h"
#include "glop/lp_types.h"
#include "lp_data/lp_data.h"
#include "lp_data/lp_types.h"
namespace operations_research {
namespace glop {
@@ -239,4 +239,4 @@ class MPSReader {
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_GLOP_MPS_READER_H_
#endif // OR_TOOLS_LP_DATA_MPS_READER_H_

6
src/glop/mps_to_png.cc → src/lp_data/mps_to_png.cc
View File

@@ -22,9 +22,9 @@
#include "base/commandlineflags.h"
#include "base/logging.h"
#include "base/file.h"
#include "glop/mps_reader.h"
#include "glop/lp_data.h"
#include "glop/png_dump.h"
#include "lp_data/mps_reader.h"
#include "lp_data/lp_data.h"
#include "lp_data/png_dump.h"
#include "base/status.h"
DEFINE_string(mps_file, "", "MPS input file.");

10
src/glop/permutation.h → src/lp_data/permutation.h
View File

@@ -12,11 +12,11 @@
// limitations under the License.
#ifndef OR_TOOLS_GLOP_PERMUTATION_H_
#define OR_TOOLS_GLOP_PERMUTATION_H_
#ifndef OR_TOOLS_LP_DATA_PERMUTATION_H_
#define OR_TOOLS_LP_DATA_PERMUTATION_H_
#include "glop/lp_types.h"
#include "glop/status.h"
#include "lp_data/lp_types.h"
#include "util/return_macros.h"
namespace operations_research {
namespace glop {
@@ -220,4 +220,4 @@ void ApplyInversePermutation(const Permutation<IndexType>& perm,
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_GLOP_PERMUTATION_H_
#endif // OR_TOOLS_LP_DATA_PERMUTATION_H_

8
src/glop/png_dump.cc → src/lp_data/png_dump.cc
View File

@@ -12,13 +12,13 @@
// limitations under the License.
#include "glop/png_dump.h"
#include "lp_data/png_dump.h"
#include "image/base/rawimage.h"
#include "image/codec/pngencoder.h"
#include "glop/lp_data.h"
#include "glop/lp_types.h"
#include "glop/sparse.h"
#include "lp_data/lp_data.h"
#include "lp_data/lp_types.h"
#include "lp_data/sparse.h"
namespace operations_research {
namespace glop {

8
src/glop/png_dump.h → src/lp_data/png_dump.h
View File

@@ -12,12 +12,12 @@
// limitations under the License.
#ifndef OR_TOOLS_GLOP_PNG_DUMP_H_
#define OR_TOOLS_GLOP_PNG_DUMP_H_
#ifndef OR_TOOLS_LP_DATA_PNG_DUMP_H_
#define OR_TOOLS_LP_DATA_PNG_DUMP_H_
#include <string>
#include "glop/lp_data.h"
#include "lp_data/lp_data.h"
namespace operations_research {
namespace glop {
@@ -28,4 +28,4 @@ std::string DumpConstraintMatrixToPng(const LinearProgram& linear_program);
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_GLOP_PNG_DUMP_H_
#endif // OR_TOOLS_LP_DATA_PNG_DUMP_H_

6
src/glop/sparse.cc → src/lp_data/sparse.cc
View File

@@ -16,9 +16,9 @@
#include "base/stringprintf.h"
#include "base/join.h"
#include "glop/lp_data.h"
#include "glop/sparse.h"
#include "glop/status.h"
#include "lp_data/lp_data.h"
#include "lp_data/sparse.h"
#include "util/return_macros.h"
namespace operations_research {
namespace glop {

22
src/glop/sparse.h → src/lp_data/sparse.h
View File

@@ -24,16 +24,16 @@
//
// Both books also contain a wealth of references.
#ifndef OR_TOOLS_GLOP_SPARSE_H_
#define OR_TOOLS_GLOP_SPARSE_H_
#ifndef OR_TOOLS_LP_DATA_SPARSE_H_
#define OR_TOOLS_LP_DATA_SPARSE_H_
#include <string>
#include "base/integral_types.h"
#include "glop/lp_types.h"
#include "glop/permutation.h"
#include "glop/sparse_column.h"
#include "glop/status.h"
#include "lp_data/lp_types.h"
#include "lp_data/permutation.h"
#include "lp_data/sparse_column.h"
#include "util/return_macros.h"
namespace operations_research {
namespace glop {
@@ -360,9 +360,13 @@ class CompactSparseMatrix {
};
ColumnView column(ColIndex col) const {
DCHECK_LT(col, num_cols_);
// Note that the start may be equals to row.size() if the last columns
// are empty, it is why we don't use &row[start].
const EntryIndex start = starts_[col];
return ColumnView(starts_[col + 1] - start, &(rows_[start]),
&(coefficients_[start]));
return ColumnView(starts_[col + 1] - start, rows_.data() + start.value(),
coefficients_.data() + start.value());
}
// Returns true if the given column is empty. Note that for triangular matrix
@@ -723,4 +727,4 @@ class TriangularMatrix : private CompactSparseMatrix {
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_GLOP_SPARSE_H_
#endif // OR_TOOLS_LP_DATA_SPARSE_H_

5
src/glop/sparse_column.cc → src/lp_data/sparse_column.cc
View File

@@ -14,9 +14,8 @@
#include <algorithm>
#include "base/stringprintf.h"
#include "glop/lp_types.h"
#include "glop/sparse_column.h"
#include "glop/status.h"
#include "lp_data/lp_types.h"
#include "lp_data/sparse_column.h"
namespace operations_research {
namespace glop {

8
src/glop/sparse_column.h → src/lp_data/sparse_column.h
View File

@@ -11,10 +11,10 @@
// See the License for the specific language governing permissions and
// limitations under the License.
#ifndef OR_TOOLS_GLOP_SPARSE_COLUMN_H_
#define OR_TOOLS_GLOP_SPARSE_COLUMN_H_
#ifndef OR_TOOLS_LP_DATA_SPARSE_COLUMN_H_
#define OR_TOOLS_LP_DATA_SPARSE_COLUMN_H_
#include "glop/sparse_vector.h"
#include "lp_data/sparse_vector.h"
namespace operations_research {
namespace glop {
@@ -134,4 +134,4 @@ class RandomAccessSparseColumn {
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_GLOP_SPARSE_COLUMN_H_
#endif // OR_TOOLS_LP_DATA_SPARSE_COLUMN_H_

12
src/glop/sparse_vector.h → src/lp_data/sparse_vector.h
View File

@@ -25,18 +25,18 @@
//
// Both books also contain a wealth of references.
#ifndef OR_TOOLS_GLOP_SPARSE_VECTOR_H_
#define OR_TOOLS_GLOP_SPARSE_VECTOR_H_
#ifndef OR_TOOLS_LP_DATA_SPARSE_VECTOR_H_
#define OR_TOOLS_LP_DATA_SPARSE_VECTOR_H_
#include <algorithm>
#include <string>
#include "base/logging.h" // for CHECK*
#include "base/integral_types.h"
#include "glop/lp_types.h"
#include "glop/permutation.h"
#include "glop/status.h"
#include "lp_data/lp_types.h"
#include "lp_data/permutation.h"
#include "util/iterators.h"
#include "util/return_macros.h"
namespace operations_research {
namespace glop {
@@ -796,4 +796,4 @@ std::string SparseVector<IndexType>::DebugString() const {
} // namespace glop
} // namespace operations_research
#endif // OR_TOOLS_GLOP_SPARSE_VECTOR_H_
#endif // OR_TOOLS_LP_DATA_SPARSE_VECTOR_H_

1
src/sat/boolean_problem.h
View File

@@ -100,6 +100,7 @@ void MakeAllLiteralsPositive(LinearBooleanProblem* problem);
void FindLinearBooleanProblemSymmetries(
const LinearBooleanProblem& problem,
std::vector<std::unique_ptr<SparsePermutation>>* generators);
// Maps all the literals of the problem. Note that this converts the cost of a
// variable correctly, that is if a variable with cost is mapped to another, the
// cost of the later is updated.

2
src/sat/sat_solver.h
View File

@@ -141,7 +141,7 @@ class SatSolver {
// will also be the unique index associated to the next constraint that will
// be added. This unique index is used by UnsatCore() to indicates what
// constraints are part of the core.
int NumAddedConstraints() { return num_constraints_; }
int NumAddedConstraints() const { return num_constraints_; }
// Gives a hint so the solver tries to find a solution with the given literal
// set to true. Currently this take precedence over the phase saving heuristic