From aeb6434297e029e2c19855ce1cbfde06e811d98c Mon Sep 17 00:00:00 2001 From: "lperron@google.com" Date: Tue, 16 Dec 2014 11:35:09 +0000 Subject: [PATCH] split lp_data out of glop; port rest of code --- examples/cpp/mps_driver.cc | 4 +- examples/cpp/opb_reader.h | 2 - examples/cpp/sat_cnf_reader.h | 5 - examples/cpp/sat_runner.cc | 2 +- makefiles/Makefile.cpp.mk | 96 ++++++----- src/base/status.h | 2 + src/glop/basis_representation.cc | 38 +++- src/glop/basis_representation.h | 17 +- src/glop/dual_edge_norms.cc | 2 +- src/glop/dual_edge_norms.h | 4 +- src/glop/entering_variable.cc | 2 +- src/glop/entering_variable.h | 6 +- src/glop/initial_basis.cc | 2 +- src/glop/initial_basis.h | 6 +- src/glop/lp_solver.cc | 194 ++++++++++++--------- src/glop/lp_solver.h | 37 ++-- src/glop/lu_factorization.cc | 7 +- src/glop/lu_factorization.h | 16 +- src/glop/markowitz.cc | 2 +- src/glop/markowitz.h | 6 +- src/glop/mps_driver.cc | 161 ----------------- src/glop/parameters.proto | 7 + src/glop/preprocessor.cc | 73 ++++---- src/glop/preprocessor.h | 40 +++-- src/glop/primal_edge_norms.cc | 9 +- src/glop/primal_edge_norms.h | 12 +- src/glop/proto_utils.h | 2 +- src/glop/rank_one_update.h | 20 ++- src/glop/reduced_costs.cc | 2 +- src/glop/reduced_costs.h | 4 +- src/glop/revised_simplex.cc | 19 +- src/glop/revised_simplex.h | 13 +- src/glop/status.h | 15 -- src/glop/update_row.cc | 11 +- src/glop/update_row.h | 11 +- src/glop/variable_values.cc | 2 +- src/glop/variable_values.h | 2 +- src/glop/variables_info.h | 4 +- src/linear_solver/glop_interface.cc | 6 +- src/linear_solver/gurobi_interface.cc | 5 + src/linear_solver/linear_solver.cc | 19 +- src/linear_solver/linear_solver.h | 9 + src/linear_solver/linear_solver2.proto | 14 +- src/linear_solver/model_exporter.cc | 27 +-- src/linear_solver/scip_interface.cc | 13 +- src/{glop => lp_data}/lp_data.cc | 7 +- src/{glop => lp_data}/lp_data.h | 12 +- src/{glop => lp_data}/lp_print_utils.cc | 4 +- src/{glop => lp_data}/lp_print_utils.h | 8 +- src/{glop => lp_data}/lp_types.cc | 2 +- src/{glop => lp_data}/lp_types.h | 14 +- src/{glop => lp_data}/lp_utils.cc | 2 +- src/{glop => lp_data}/lp_utils.h | 10 +- src/{glop => lp_data}/matrix_scaler.cc | 7 +- src/{glop => lp_data}/matrix_scaler.h | 8 +- src/lp_data/matrix_utils.cc | 220 ++++++++++++++++++++++++ src/lp_data/matrix_utils.h | 55 ++++++ src/{glop => lp_data}/mps_reader.cc | 5 +- src/{glop => lp_data}/mps_reader.h | 10 +- src/{glop => lp_data}/mps_to_png.cc | 6 +- src/{glop => lp_data}/permutation.h | 10 +- src/{glop => lp_data}/png_dump.cc | 8 +- src/{glop => lp_data}/png_dump.h | 8 +- src/{glop => lp_data}/sparse.cc | 6 +- src/{glop => lp_data}/sparse.h | 22 ++- src/{glop => lp_data}/sparse_column.cc | 5 +- src/{glop => lp_data}/sparse_column.h | 8 +- src/{glop => lp_data}/sparse_vector.h | 12 +- src/sat/boolean_problem.h | 1 + src/sat/sat_solver.h | 2 +- 70 files changed, 875 insertions(+), 527 deletions(-) delete mode 100644 src/glop/mps_driver.cc rename src/{glop => lp_data}/lp_data.cc (99%) rename src/{glop => lp_data}/lp_data.h (98%) rename src/{glop => lp_data}/lp_print_utils.cc (96%) rename src/{glop => lp_data}/lp_print_utils.h (93%) rename src/{glop => lp_data}/lp_types.cc (99%) rename src/{glop => lp_data}/lp_types.h (96%) rename src/{glop => lp_data}/lp_utils.cc (99%) rename src/{glop => lp_data}/lp_utils.h (98%) rename src/{glop => lp_data}/matrix_scaler.cc (99%) rename src/{glop => lp_data}/matrix_scaler.h (97%) create mode 100644 src/lp_data/matrix_utils.cc create mode 100644 src/lp_data/matrix_utils.h rename src/{glop => lp_data}/mps_reader.cc (99%) rename src/{glop => lp_data}/mps_reader.h (97%) rename src/{glop => lp_data}/mps_to_png.cc (94%) rename src/{glop => lp_data}/permutation.h (97%) rename src/{glop => lp_data}/png_dump.cc (93%) rename src/{glop => lp_data}/png_dump.h (85%) rename src/{glop => lp_data}/sparse.cc (99%) rename src/{glop => lp_data}/sparse.h (98%) rename src/{glop => lp_data}/sparse_column.cc (96%) rename src/{glop => lp_data}/sparse_column.h (96%) rename src/{glop => lp_data}/sparse_vector.h (99%) diff --git a/examples/cpp/mps_driver.cc b/examples/cpp/mps_driver.cc index b184a5f270..685b4a6b64 100644 --- a/examples/cpp/mps_driver.cc +++ b/examples/cpp/mps_driver.cc @@ -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" diff --git a/examples/cpp/opb_reader.h b/examples/cpp/opb_reader.h index ab9929c52e..04ea8558a6 100644 --- a/examples/cpp/opb_reader.h +++ b/examples/cpp/opb_reader.h @@ -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]; diff --git a/examples/cpp/sat_cnf_reader.h b/examples/cpp/sat_cnf_reader.h index ea77f98e35..0174edcfd5 100644 --- a/examples/cpp/sat_cnf_reader.h +++ b/examples/cpp/sat_cnf_reader.h @@ -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]; diff --git a/examples/cpp/sat_runner.cc b/examples/cpp/sat_runner.cc index f45da6d921..3fc6053d80 100644 --- a/examples/cpp/sat_runner.cc +++ b/examples/cpp/sat_runner.cc @@ -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! diff --git a/makefiles/Makefile.cpp.mk b/makefiles/Makefile.cpp.mk index fe18dba5d9..f4f8fc3319 100644 --- a/makefiles/Makefile.cpp.mk +++ b/makefiles/Makefile.cpp.mk @@ -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 diff --git a/src/base/status.h b/src/base/status.h index 20e3ed1603..7546482d3e 100644 --- a/src/base/status.h +++ b/src/base/status.h @@ -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: diff --git a/src/glop/basis_representation.cc b/src/glop/basis_representation.cc index d83653b1e4..2dcb43eee6 100644 --- a/src/glop/basis_representation.cc +++ b/src/glop/basis_representation.cc @@ -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* 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* 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(num_entries) / + static_cast(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 diff --git a/src/glop/basis_representation.h b/src/glop/basis_representation.h index b09fd606a3..d67e001810 100644 --- a/src/glop/basis_representation.h +++ b/src/glop/basis_representation.h @@ -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); }; diff --git a/src/glop/dual_edge_norms.cc b/src/glop/dual_edge_norms.cc index bbe48ce767..ba6d7f0ed6 100644 --- a/src/glop/dual_edge_norms.cc +++ b/src/glop/dual_edge_norms.cc @@ -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 { diff --git a/src/glop/dual_edge_norms.h b/src/glop/dual_edge_norms.h index e53a6474d7..281547a853 100644 --- a/src/glop/dual_edge_norms.h +++ b/src/glop/dual_edge_norms.h @@ -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 { diff --git a/src/glop/entering_variable.cc b/src/glop/entering_variable.cc index 359f3ae3d5..22fc4ab1a2 100644 --- a/src/glop/entering_variable.cc +++ b/src/glop/entering_variable.cc @@ -17,7 +17,7 @@ #include #include "base/timer.h" -#include "glop/lp_utils.h" +#include "lp_data/lp_utils.h" namespace operations_research { namespace glop { diff --git a/src/glop/entering_variable.h b/src/glop/entering_variable.h index 172c565550..b00390ab69 100644 --- a/src/glop/entering_variable.h +++ b/src/glop/entering_variable.h @@ -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" diff --git a/src/glop/initial_basis.cc b/src/glop/initial_basis.cc index d9fe379465..71d1ec5c91 100644 --- a/src/glop/initial_basis.cc +++ b/src/glop/initial_basis.cc @@ -15,8 +15,8 @@ #include #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 { diff --git a/src/glop/initial_basis.h b/src/glop/initial_basis.h index 85749d19be..78f1da0dc4 100644 --- a/src/glop/initial_basis.h +++ b/src/glop/initial_basis.h @@ -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 { diff --git a/src/glop/lp_solver.cc b/src/glop/lp_solver.cc index a04e6c62f3..3dc0ef181f 100644 --- a/src/glop/lp_solver.cc +++ b/src/glop/lp_solver.cc @@ -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 diff --git a/src/glop/lp_solver.h b/src/glop/lp_solver.h index 09fbfaa7cf..3b6d8e8806 100644 --- a/src/glop/lp_solver.h +++ b/src/glop/lp_solver.h @@ -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_; diff --git a/src/glop/lu_factorization.cc b/src/glop/lu_factorization.cc index 041a8f4653..dc2004eb79 100644 --- a/src/glop/lu_factorization.cc +++ b/src/glop/lu_factorization.cc @@ -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(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()); diff --git a/src/glop/lu_factorization.h b/src/glop/lu_factorization.h index 857db38f15..f0d79be791 100644 --- a/src/glop/lu_factorization.h +++ b/src/glop/lu_factorization.h @@ -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. diff --git a/src/glop/markowitz.cc b/src/glop/markowitz.cc index 5fcc95ba80..2c4c2bf49b 100644 --- a/src/glop/markowitz.cc +++ b/src/glop/markowitz.cc @@ -15,7 +15,7 @@ #include #include "base/stringprintf.h" -#include "glop/lp_utils.h" +#include "lp_data/lp_utils.h" namespace operations_research { namespace glop { diff --git a/src/glop/markowitz.h b/src/glop/markowitz.h index ffea067be4..9c0b87f481 100644 --- a/src/glop/markowitz.h +++ b/src/glop/markowitz.h @@ -76,10 +76,10 @@ #include #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 { diff --git a/src/glop/mps_driver.cc b/src/glop/mps_driver.cc deleted file mode 100644 index 5c16563359..0000000000 --- a/src/glop/mps_driver.cc +++ /dev/null @@ -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 -#include - -#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(¶meters); - - LinearProgram linear_program; - std::vector 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; -} diff --git a/src/glop/parameters.proto b/src/glop/parameters.proto index 4a9faaa436..49bd775458 100644 --- a/src/glop/parameters.proto +++ b/src/glop/parameters.proto @@ -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]; diff --git a/src/glop/preprocessor.cc b/src/glop/preprocessor.cc index f70d65af80..2ae3988052 100644 --- a/src/glop/preprocessor.cc +++ b/src/glop/preprocessor.cc @@ -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 column_degree(num_cols, EntryIndex(0)); std::vector 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 row_degree(num_rows, EntryIndex(0)); std::vector 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); } } diff --git a/src/glop/preprocessor.h b/src/glop/preprocessor.h index cad3254d6a..9e00b7d5f1 100644 --- a/src/glop/preprocessor.h +++ b/src/glop/preprocessor.h @@ -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 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); }; diff --git a/src/glop/primal_edge_norms.cc b/src/glop/primal_edge_norms.cc index 21c4b8dec4..52a754006f 100644 --- a/src/glop/primal_edge_norms.cc +++ b/src/glop/primal_edge_norms.cc @@ -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 thread_local_stat_lower_bounded_norms(num_omp_threads, 0); const std::vector& relevant_rows = update_row.GetNonZeroPositions(); const int parallel_loop_size = relevant_rows.size(); diff --git a/src/glop/primal_edge_norms.h b/src/glop/primal_edge_norms.h index 6874637332..000cd5d4d6 100644 --- a/src/glop/primal_edge_norms.h +++ b/src/glop/primal_edge_norms.h @@ -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); }; diff --git a/src/glop/proto_utils.h b/src/glop/proto_utils.h index 8fb60213fc..7f26a714ef 100644 --- a/src/glop/proto_utils.h +++ b/src/glop/proto_utils.h @@ -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 { diff --git a/src/glop/rank_one_update.h b/src/glop/rank_one_update.h index 3df7a5b9c3..1b5a8b93cb 100644 --- a/src/glop/rank_one_update.h +++ b/src/glop/rank_one_update.h @@ -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(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 elementary_matrices_; DISALLOW_COPY_AND_ASSIGN(RankOneUpdateFactorization); }; diff --git a/src/glop/reduced_costs.cc b/src/glop/reduced_costs.cc index e1eb4087c4..42f495d294 100644 --- a/src/glop/reduced_costs.cc +++ b/src/glop/reduced_costs.cc @@ -17,7 +17,7 @@ #include #endif -#include "glop/lp_utils.h" +#include "lp_data/lp_utils.h" namespace operations_research { namespace glop { diff --git a/src/glop/reduced_costs.h b/src/glop/reduced_costs.h index 12ab2dba0b..25140be20c 100644 --- a/src/glop/reduced_costs.h +++ b/src/glop/reduced_costs.h @@ -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 { diff --git a/src/glop/revised_simplex.cc b/src/glop/revised_simplex.cc index 4f87e24bdd..1f6c77b03b 100644 --- a/src/glop/revised_simplex.cc +++ b/src/glop/revised_simplex.cc @@ -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; } diff --git a/src/glop/revised_simplex.h b/src/glop/revised_simplex.h index 2ccaccbf18..4a0c18da7c 100644 --- a/src/glop/revised_simplex.h +++ b/src/glop/revised_simplex.h @@ -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); }; diff --git a/src/glop/status.h b/src/glop/status.h index 93b8b946e0..dac99de138 100644 --- a/src/glop/status.h +++ b/src/glop/status.h @@ -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 diff --git a/src/glop/update_row.cc b/src/glop/update_row.cc index db95c2abe7..74f7f5fef5 100644 --- a/src/glop/update_row.cc +++ b/src/glop/update_row.cc @@ -17,7 +17,7 @@ #include #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(num_col_wise_entries.value())) { if (lhs < 1.1 * static_cast(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(non_zero_position_list_.size()) / diff --git a/src/glop/update_row.h b/src/glop/update_row.h index e3dad95185..9928f3c104 100644 --- a/src/glop/update_row.h +++ b/src/glop/update_row.h @@ -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_; diff --git a/src/glop/variable_values.cc b/src/glop/variable_values.cc index d1bbbc9b39..dc9bd5106e 100644 --- a/src/glop/variable_values.cc +++ b/src/glop/variable_values.cc @@ -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 { diff --git a/src/glop/variable_values.h b/src/glop/variable_values.h index e93456e0d4..0b1f9458d2 100644 --- a/src/glop/variable_values.h +++ b/src/glop/variable_values.h @@ -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 { diff --git a/src/glop/variables_info.h b/src/glop/variables_info.h index 38b2eac59e..da242342aa 100644 --- a/src/glop/variables_info.h +++ b/src/glop/variables_info.h @@ -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 { diff --git a/src/linear_solver/glop_interface.cc b/src/linear_solver/glop_interface.cc index a947198a94..beffb2a02c 100644 --- a/src/linear_solver/glop_interface.cc +++ b/src/linear_solver/glop_interface.cc @@ -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(solver_->time_limit())); + static_cast(solver_->time_limit()) / 1000.0); } solver_->SetSolverSpecificParametersAsString( diff --git a/src/linear_solver/gurobi_interface.cc b/src/linear_solver/gurobi_interface.cc index b6aca53303..6feb39fee0 100644 --- a/src/linear_solver/gurobi_interface.cc +++ b/src/linear_solver/gurobi_interface.cc @@ -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() { diff --git a/src/linear_solver/linear_solver.cc b/src/linear_solver/linear_solver.cc index 949fa9d1d3..1cea15167f 100644 --- a/src/linear_solver/linear_solver.cc +++ b/src/linear_solver/linear_solver.cc @@ -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( diff --git a/src/linear_solver/linear_solver.h b/src/linear_solver/linear_solver.h index 78901c9dbd..fa46f741b8 100644 --- a/src/linear_solver/linear_solver.h +++ b/src/linear_solver/linear_solver.h @@ -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. diff --git a/src/linear_solver/linear_solver2.proto b/src/linear_solver/linear_solver2.proto index 2901f76058..7c0418cdc1 100644 --- a/src/linear_solver/linear_solver2.proto +++ b/src/linear_solver/linear_solver2.proto @@ -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. diff --git a/src/linear_solver/model_exporter.cc b/src/linear_solver/model_exporter.cc index dfad4abe77..d1ac97f7f8 100644 --- a/src/linear_solver/model_exporter.cc +++ b/src/linear_solver/model_exporter.cc @@ -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; diff --git a/src/linear_solver/scip_interface.cc b/src/linear_solver/scip_interface.cc index c8a2949089..cdf560a497 100644 --- a/src/linear_solver/scip_interface.cc +++ b/src/linear_solver/scip_interface.cc @@ -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(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. diff --git a/src/glop/lp_data.cc b/src/lp_data/lp_data.cc similarity index 99% rename from src/glop/lp_data.cc rename to src/lp_data/lp_data.cc index 2b4beb467e..a1c0e72135 100644 --- a/src/glop/lp_data.cc +++ b/src/lp_data/lp_data.cc @@ -12,7 +12,7 @@ // limitations under the License. -#include "glop/lp_data.h" +#include "lp_data/lp_data.h" #include #include @@ -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 { diff --git a/src/glop/lp_data.h b/src/lp_data/lp_data.h similarity index 98% rename from src/glop/lp_data.h rename to src/lp_data/lp_data.h index 03e81a77ff..476d3fbb07 100644 --- a/src/glop/lp_data.h +++ b/src/lp_data/lp_data.h @@ -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 // for max #include @@ -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_ diff --git a/src/glop/lp_print_utils.cc b/src/lp_data/lp_print_utils.cc similarity index 96% rename from src/glop/lp_print_utils.cc rename to src/lp_data/lp_print_utils.cc index 1b64f0bc01..617348a30a 100644 --- a/src/glop/lp_print_utils.cc +++ b/src/lp_data/lp_print_utils.cc @@ -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 { diff --git a/src/glop/lp_print_utils.h b/src/lp_data/lp_print_utils.h similarity index 93% rename from src/glop/lp_print_utils.h rename to src/lp_data/lp_print_utils.h index 3735df9e72..043b16dc78 100644 --- a/src/glop/lp_print_utils.h +++ b/src/lp_data/lp_print_utils.h @@ -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 #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_ diff --git a/src/glop/lp_types.cc b/src/lp_data/lp_types.cc similarity index 99% rename from src/glop/lp_types.cc rename to src/lp_data/lp_types.cc index fd386e84e8..9f73707c55 100644 --- a/src/glop/lp_types.cc +++ b/src/lp_data/lp_types.cc @@ -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 { diff --git a/src/glop/lp_types.h b/src/lp_data/lp_types.h similarity index 96% rename from src/glop/lp_types.h rename to src/lp_data/lp_types.h index d98979a051..d5076d773e 100644 --- a/src/glop/lp_types.h +++ b/src/lp_data/lp_types.h @@ -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 #include @@ -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(n); +} + } // namespace glop } // namespace operations_research -#endif // OR_TOOLS_GLOP_LP_TYPES_H_ +#endif // OR_TOOLS_LP_DATA_LP_TYPES_H_ diff --git a/src/glop/lp_utils.cc b/src/lp_data/lp_utils.cc similarity index 99% rename from src/glop/lp_utils.cc rename to src/lp_data/lp_utils.cc index 4a642280f3..9520dcf5cd 100644 --- a/src/glop/lp_utils.cc +++ b/src/lp_data/lp_utils.cc @@ -12,7 +12,7 @@ // limitations under the License. -#include "glop/lp_utils.h" +#include "lp_data/lp_utils.h" namespace operations_research { namespace glop { diff --git a/src/glop/lp_utils.h b/src/lp_data/lp_utils.h similarity index 98% rename from src/glop/lp_utils.h rename to src/lp_data/lp_utils.h index 66bdc0e547..5691e3129f 100644 --- a/src/glop/lp_utils.h +++ b/src/lp_data/lp_utils.h @@ -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 SumWithNegativeInfiniteAndOneMissing; } // namespace glop } // namespace operations_research -#endif // OR_TOOLS_GLOP_LP_UTILS_H_ +#endif // OR_TOOLS_LP_DATA_LP_UTILS_H_ diff --git a/src/glop/matrix_scaler.cc b/src/lp_data/matrix_scaler.cc similarity index 99% rename from src/glop/matrix_scaler.cc rename to src/lp_data/matrix_scaler.cc index 44c903ee31..70a2917cee 100644 --- a/src/glop/matrix_scaler.cc +++ b/src/lp_data/matrix_scaler.cc @@ -12,7 +12,7 @@ // limitations under the License. -#include "glop/matrix_scaler.h" +#include "lp_data/matrix_scaler.h" #include #include @@ -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 { diff --git a/src/glop/matrix_scaler.h b/src/lp_data/matrix_scaler.h similarity index 97% rename from src/glop/matrix_scaler.h rename to src/lp_data/matrix_scaler.h index 1f3158bcfa..3f2491cacb 100644 --- a/src/glop/matrix_scaler.h +++ b/src/lp_data/matrix_scaler.h @@ -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 #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_ diff --git a/src/lp_data/matrix_utils.cc b/src/lp_data/matrix_utils.cc new file mode 100644 index 0000000000..b4c1440618 --- /dev/null +++ b/src/lp_data/matrix_utils.cc @@ -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 +#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::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(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 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 diff --git a/src/lp_data/matrix_utils.h b/src/lp_data/matrix_utils.h new file mode 100644 index 0000000000..f0193cd7de --- /dev/null +++ b/src/lp_data/matrix_utils.h @@ -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_ diff --git a/src/glop/mps_reader.cc b/src/lp_data/mps_reader.cc similarity index 99% rename from src/glop/mps_reader.cc rename to src/lp_data/mps_reader.cc index c1d81fc12e..87dd662cde 100644 --- a/src/glop/mps_reader.cc +++ b/src/lp_data/mps_reader.cc @@ -12,7 +12,7 @@ // limitations under the License. -#include "glop/mps_reader.h" +#include "lp_data/mps_reader.h" #include #include @@ -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."); diff --git a/src/glop/mps_reader.h b/src/lp_data/mps_reader.h similarity index 97% rename from src/glop/mps_reader.h rename to src/lp_data/mps_reader.h index 444429dfbd..67e7fef8d0 100644 --- a/src/glop/mps_reader.h +++ b/src/lp_data/mps_reader.h @@ -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 // for max #include @@ -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_ diff --git a/src/glop/mps_to_png.cc b/src/lp_data/mps_to_png.cc similarity index 94% rename from src/glop/mps_to_png.cc rename to src/lp_data/mps_to_png.cc index 89a22a5445..6dae876350 100644 --- a/src/glop/mps_to_png.cc +++ b/src/lp_data/mps_to_png.cc @@ -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."); diff --git a/src/glop/permutation.h b/src/lp_data/permutation.h similarity index 97% rename from src/glop/permutation.h rename to src/lp_data/permutation.h index c3c2d91e43..637b4d0434 100644 --- a/src/glop/permutation.h +++ b/src/lp_data/permutation.h @@ -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& perm, } // namespace glop } // namespace operations_research -#endif // OR_TOOLS_GLOP_PERMUTATION_H_ +#endif // OR_TOOLS_LP_DATA_PERMUTATION_H_ diff --git a/src/glop/png_dump.cc b/src/lp_data/png_dump.cc similarity index 93% rename from src/glop/png_dump.cc rename to src/lp_data/png_dump.cc index 391ecea2e8..4c4686d8e2 100644 --- a/src/glop/png_dump.cc +++ b/src/lp_data/png_dump.cc @@ -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 { diff --git a/src/glop/png_dump.h b/src/lp_data/png_dump.h similarity index 85% rename from src/glop/png_dump.h rename to src/lp_data/png_dump.h index e955dd545f..fdf1e7be53 100644 --- a/src/glop/png_dump.h +++ b/src/lp_data/png_dump.h @@ -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 -#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_ diff --git a/src/glop/sparse.cc b/src/lp_data/sparse.cc similarity index 99% rename from src/glop/sparse.cc rename to src/lp_data/sparse.cc index b62bb2089d..86e87c84ec 100644 --- a/src/glop/sparse.cc +++ b/src/lp_data/sparse.cc @@ -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 { diff --git a/src/glop/sparse.h b/src/lp_data/sparse.h similarity index 98% rename from src/glop/sparse.h rename to src/lp_data/sparse.h index 0d4f2309f5..6182d6fdfa 100644 --- a/src/glop/sparse.h +++ b/src/lp_data/sparse.h @@ -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 #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_ diff --git a/src/glop/sparse_column.cc b/src/lp_data/sparse_column.cc similarity index 96% rename from src/glop/sparse_column.cc rename to src/lp_data/sparse_column.cc index 7b089cf986..17f654135d 100644 --- a/src/glop/sparse_column.cc +++ b/src/lp_data/sparse_column.cc @@ -14,9 +14,8 @@ #include #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 { diff --git a/src/glop/sparse_column.h b/src/lp_data/sparse_column.h similarity index 96% rename from src/glop/sparse_column.h rename to src/lp_data/sparse_column.h index adaa6f3d44..8a5851b786 100644 --- a/src/glop/sparse_column.h +++ b/src/lp_data/sparse_column.h @@ -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_ diff --git a/src/glop/sparse_vector.h b/src/lp_data/sparse_vector.h similarity index 99% rename from src/glop/sparse_vector.h rename to src/lp_data/sparse_vector.h index 0551379281..3e73449ef3 100644 --- a/src/glop/sparse_vector.h +++ b/src/lp_data/sparse_vector.h @@ -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 #include #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::DebugString() const { } // namespace glop } // namespace operations_research -#endif // OR_TOOLS_GLOP_SPARSE_VECTOR_H_ +#endif // OR_TOOLS_LP_DATA_SPARSE_VECTOR_H_ diff --git a/src/sat/boolean_problem.h b/src/sat/boolean_problem.h index 16e1a50f62..236aae7ada 100644 --- a/src/sat/boolean_problem.h +++ b/src/sat/boolean_problem.h @@ -100,6 +100,7 @@ void MakeAllLiteralsPositive(LinearBooleanProblem* problem); void FindLinearBooleanProblemSymmetries( const LinearBooleanProblem& problem, std::vector>* 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. diff --git a/src/sat/sat_solver.h b/src/sat/sat_solver.h index ecf1bd0997..563fc74127 100644 --- a/src/sat/sat_solver.h +++ b/src/sat/sat_solver.h @@ -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