diff --git a/src/glop/lp_solver.h b/src/glop/lp_solver.h index 5a715f1a92..946b897963 100644 --- a/src/glop/lp_solver.h +++ b/src/glop/lp_solver.h @@ -64,8 +64,8 @@ class LPSolver { void Clear(); // Advanced usage. This should be called before calling Solve(). It will - // configure the solver to try to start from the given point. Note that - // calling Clear() will invalidate this information. + // configure the solver to try to start from the given point for the next + // Solve() only. Note that calling Clear() will invalidate this information. // // If the set of variables/constraints with a BASIC status does not form a // basis a warning will be logged and the code will ignore it. Otherwise, the diff --git a/src/glop/lu_factorization.cc b/src/glop/lu_factorization.cc index a0f63f4eab..167cf2b570 100644 --- a/src/glop/lu_factorization.cc +++ b/src/glop/lu_factorization.cc @@ -602,20 +602,6 @@ Fractional LuFactorization::ComputeInfinityNormConditionNumber( return matrix.ComputeInfinityNorm() * ComputeInverseInfinityNorm(); } -namespace { -// Returns the density of the sparse column 'b' w.r.t. the given permutation. -double ComputeDensity(const SparseColumn& b, const RowPermutation& row_perm) { - double density = 0.0; - for (const SparseColumn::Entry e : b) { - if (row_perm[e.row()] != kNonPivotal && e.coefficient() != 0.0) { - ++density; - } - } - const RowIndex num_rows = row_perm.size(); - return density / num_rows.value(); -} -} // anonymous namespace - void LuFactorization::ComputeTransposeUpper() { SCOPED_TIME_STAT(&stats_); transpose_upper_.PopulateFromTranspose(upper_); diff --git a/src/glop/preprocessor.cc b/src/glop/preprocessor.cc index 748b79f3a1..6e131c07b5 100644 --- a/src/glop/preprocessor.cc +++ b/src/glop/preprocessor.cc @@ -24,10 +24,6 @@ namespace operations_research { namespace glop { namespace { -// Returns an interval as an human readable std::string for debugging. -std::string IntervalString(Fractional lb, Fractional ub) { - return StringPrintf("[%g, %g]", lb, ub); -} #if defined(_MSC_VER) double trunc(double d) { return d > 0 ? floor(d) : ceil(d); } @@ -104,14 +100,15 @@ bool MainLpPreprocessor::Run(LinearProgram* lp, TimeLimit* time_limit) { RUN_PREPROCESSOR(EmptyColumnPreprocessor); RUN_PREPROCESSOR(EmptyConstraintPreprocessor); } + + RUN_PREPROCESSOR(SingletonColumnSignPreprocessor); } - // These are implemented as preprocessors, but are not controlled by the - // use_preprocessing() parameter. - RUN_PREPROCESSOR(SingletonColumnSignPreprocessor); + // The scaling is controled by use_scaling, not use_preprocessing. RUN_PREPROCESSOR(ScalingPreprocessor); - RUN_PREPROCESSOR(AddSlackVariablesPreprocessor); + // This one must always run. It is needed by the revised simplex code. + RUN_PREPROCESSOR(AddSlackVariablesPreprocessor); return !preprocessors_.empty(); } diff --git a/src/glop/primal_edge_norms.h b/src/glop/primal_edge_norms.h index 000cd5d4d6..5faac8e679 100644 --- a/src/glop/primal_edge_norms.h +++ b/src/glop/primal_edge_norms.h @@ -172,7 +172,6 @@ class PrimalEdgeNorms { Stats stats_; // Booleans to control what happens on the next ChooseEnteringColumn() call. - bool must_refactorize_basis_; bool recompute_edge_squared_norms_; bool reset_devex_weights_; diff --git a/src/glop/rank_one_update.h b/src/glop/rank_one_update.h index a483d76458..d519d621ee 100644 --- a/src/glop/rank_one_update.h +++ b/src/glop/rank_one_update.h @@ -119,13 +119,6 @@ class RankOneUpdateElementaryMatrix { } private: - // This is only used in debug mode. - Fractional ComputeUScalarV() const { - DenseColumn dense_u; - storage_->ColumnCopyToDenseColumn(u_index_, &dense_u); - return storage_->ColumnScalarProduct(v_index_, Transpose(dense_u)); - } - // Note that we allow copy and assignment so we can store a // RankOneUpdateElementaryMatrix in an STL container. const CompactSparseMatrix* storage_; diff --git a/src/glop/revised_simplex.cc b/src/glop/revised_simplex.cc index e7735fd3cb..45de38ac41 100644 --- a/src/glop/revised_simplex.cc +++ b/src/glop/revised_simplex.cc @@ -672,19 +672,29 @@ void RevisedSimplex::UseSingletonColumnInInitialBasis(RowToColMapping* basis) { } } -bool RevisedSimplex::InitializeMatrixAndTestIfUnchanged(const LinearProgram& lp, - bool* only_new_rows) { +bool RevisedSimplex::InitializeMatrixAndTestIfUnchanged( + const LinearProgram& lp, bool* only_change_is_new_rows, + bool* only_change_is_new_cols, ColIndex* num_new_cols) { SCOPED_TIME_STAT(&function_stats_); - DCHECK_NE(only_new_rows, nullptr); + DCHECK_NE(only_change_is_new_rows, nullptr); + DCHECK_NE(only_change_is_new_cols, nullptr); + DCHECK_NE(num_new_cols, nullptr); DCHECK_NE(kInvalidCol, lp.GetFirstSlackVariable()); DCHECK_EQ(num_cols_, compact_matrix_.num_cols()); DCHECK_EQ(num_rows_, compact_matrix_.num_rows()); - // Test if the matrix is unchanged, and if yes, just returns true. Note that - // we compare also the columns corresponding to the slack variables. - if (lp.num_constraints() == num_rows_ && lp.num_variables() == num_cols_ && + DCHECK_EQ(lp.num_variables(), + lp.GetFirstSlackVariable() + RowToColIndex(lp.num_constraints())); + DCHECK(IsRightMostSquareMatrixIdentity(lp.GetSparseMatrix())); + const bool old_part_of_matrix_is_unchanged = AreFirstColumnsAndRowsExactlyEquals( - num_rows_, num_cols_, lp.GetSparseMatrix(), compact_matrix_)) { + num_rows_, first_slack_col_, lp.GetSparseMatrix(), compact_matrix_); + + // Test if the matrix is unchanged, and if yes, just returns true. Note that + // this doesn't check the columns corresponding to the slack variables, + // because they were checked by lp.IsInEquationForm() when Solve() was called. + if (old_part_of_matrix_is_unchanged && lp.num_constraints() == num_rows_ && + lp.num_variables() == num_cols_) { // IMPORTANT: we need to recreate matrix_with_slack_ because this matrix // view was refering to a previous lp.GetSparseMatrix(). The matrices are // the same, but we do need to update the pointers. @@ -695,17 +705,18 @@ bool RevisedSimplex::InitializeMatrixAndTestIfUnchanged(const LinearProgram& lp, } // Check if the new matrix can be derived from the old one just by adding - // new rows (i.e new constraints). In this case, we make two separate tests: - // first we compare the coefficients for the non-slack variables, and then we - // check that the coefficients for slack variables form an identity matrix. - *only_new_rows = - lp.num_constraints() > num_rows_ && - lp.GetFirstSlackVariable() == first_slack_col_ && - lp.num_variables() == - lp.GetFirstSlackVariable() + RowToColIndex(lp.num_constraints()) && - AreFirstColumnsAndRowsExactlyEquals( - num_rows_, first_slack_col_, lp.GetSparseMatrix(), compact_matrix_) && - IsRightMostSquareMatrixIdentity(lp.GetSparseMatrix()); + // new rows (i.e new constraints). + *only_change_is_new_rows = old_part_of_matrix_is_unchanged && + lp.num_constraints() > num_rows_ && + lp.GetFirstSlackVariable() == first_slack_col_; + + // Check if the new matrix can be derived from the old one just by adding + // new columns (i.e new variables). + *only_change_is_new_cols = old_part_of_matrix_is_unchanged && + lp.num_constraints() == num_rows_ && + lp.GetFirstSlackVariable() > first_slack_col_; + *num_new_cols = + *only_change_is_new_cols ? lp.num_variables() - num_cols_ : ColIndex(0); // Initialize matrix_with_slack_. matrix_with_slack_.PopulateFromMatrix(lp.GetSparseMatrix()); @@ -725,6 +736,37 @@ bool RevisedSimplex::InitializeMatrixAndTestIfUnchanged(const LinearProgram& lp, return false; } +bool RevisedSimplex::OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero( + const LinearProgram& lp, ColIndex num_new_cols) { + SCOPED_TIME_STAT(&function_stats_); + DCHECK_EQ(lp.num_variables(), num_cols_); + DCHECK_LE(num_new_cols, first_slack_col_); + const ColIndex first_new_col(first_slack_col_ - num_new_cols); + + // Check the original variable bounds. + for (ColIndex col(0); col < first_new_col; ++col) { + if (lower_bound_[col] != lp.variable_lower_bounds()[col] || + upper_bound_[col] != lp.variable_upper_bounds()[col]) { + return false; + } + } + // Check that each new variable has a bound of zero. + for (ColIndex col(first_new_col); col < first_slack_col_; ++col) { + if (lp.variable_lower_bounds()[col] != 0.0 && + lp.variable_upper_bounds()[col] != 0.0) { + return false; + } + } + // Check that the slack bounds are unchanged. + for (ColIndex col(first_slack_col_); col < num_cols_; ++col) { + if (lower_bound_[col - num_new_cols] != lp.variable_lower_bounds()[col] || + upper_bound_[col - num_new_cols] != lp.variable_upper_bounds()[col]) { + return false; + } + } + return true; +} + bool RevisedSimplex::InitializeBoundsAndTestIfUnchanged( const LinearProgram& lp) { SCOPED_TIME_STAT(&function_stats_); @@ -821,10 +863,11 @@ void RevisedSimplex::InitializeObjectiveLimit(const LinearProgram& lp) { } void RevisedSimplex::InitializeVariableStatusesForWarmStart( - const BasisState& state) { + const BasisState& state, ColIndex num_new_cols) { variables_info_.Initialize(); RowIndex num_basic_variables(0); - + DCHECK_LE(num_new_cols, first_slack_col_); + const ColIndex first_new_col(first_slack_col_ - num_new_cols); // Compute the status for all the columns (note that the slack variables are // already added at the end of the matrix at this stage). for (ColIndex col(0); col < num_cols_; ++col) { @@ -832,8 +875,11 @@ void RevisedSimplex::InitializeVariableStatusesForWarmStart( // Start with the given "warm" status from the BasisState if it exists. VariableStatus status = default_status; - if (col < state.statuses.size()) { + if (col < first_new_col && col < state.statuses.size()) { status = state.statuses[col]; + } else if (col >= first_slack_col_ && + col - num_new_cols < state.statuses.size()) { + status = state.statuses[col - num_new_cols]; } // Remove incompatibilities between the warm status and the variable bounds. @@ -1007,9 +1053,14 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) { // // Note that these functions can't depend on use_dual_simplex() since we may // change it below. - bool only_new_rows = false; - const bool is_matrix_unchanged = - InitializeMatrixAndTestIfUnchanged(lp, &only_new_rows); + bool only_change_is_new_rows = false; + bool only_change_is_new_cols = false; + ColIndex num_new_cols(0); + const bool is_matrix_unchanged = InitializeMatrixAndTestIfUnchanged( + lp, &only_change_is_new_rows, &only_change_is_new_cols, &num_new_cols); + const bool only_new_bounds = + only_change_is_new_cols && num_new_cols > 0 && + OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero(lp, num_new_cols); const bool objective_is_unchanged = InitializeObjectiveAndTestIfUnchanged(lp); const bool bounds_are_unchanged = InitializeBoundsAndTestIfUnchanged(lp); @@ -1048,7 +1099,7 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) { if (solution_state_has_been_set_externally_) { // If an external basis has been provided we need to perform more work, // e.g., factorize and validate it. - InitializeVariableStatusesForWarmStart(solution_state_); + InitializeVariableStatusesForWarmStart(solution_state_, ColIndex(0)); basis_.assign(num_rows_, kInvalidCol); RowIndex row(0); for (ColIndex col : variables_info_.GetIsBasicBitRow()) { @@ -1065,13 +1116,14 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) { reduced_costs_.ClearAndRemoveCostShifts(); solve_from_scratch = false; } else { - VLOG(1) << "RevisedSimplex is not using the externally provided basis " - "because it is not factorizable."; + LOG(WARNING) << "RevisedSimplex is not using the externally provided " + "basis because it is not factorizable."; } } else if (!parameters_.use_dual_simplex()) { // With primal simplex, always clear dual norms and dual pricing. - // Incrementality is supported only if both the matrix and the bounds - // remain the same (objective may change). + // Incrementality is supported only if only change to the matrix and + // bounds is adding new columns (objective may change), and that all + // new columns have a bound equal to zero. dual_edge_norms_.Clear(); dual_pricing_vector_.clear(); if (is_matrix_unchanged && bounds_are_unchanged) { @@ -1079,6 +1131,20 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) { // this seems to break something. Investigate. reduced_costs_.ClearAndRemoveCostShifts(); solve_from_scratch = false; + } else if (only_change_is_new_cols && only_new_bounds) { + InitializeVariableStatusesForWarmStart(solution_state_, num_new_cols); + const ColIndex first_new_col(first_slack_col_ - num_new_cols); + for (ColIndex& col_ref : basis_) { + if (col_ref >= first_new_col) { + col_ref += num_new_cols; + } + } + // Make sure the primal edge norm are recomputed from scratch. + // TODO(user): only the norms of the new columns actually need to be + // computed. + primal_edge_norms_.Clear(); + reduced_costs_.ClearAndRemoveCostShifts(); + solve_from_scratch = false; } } else { // With dual simplex, always clear primal norms. Incrementality is @@ -1088,14 +1154,15 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) { if (objective_is_unchanged) { if (is_matrix_unchanged) { if (!bounds_are_unchanged) { - InitializeVariableStatusesForWarmStart(solution_state_); + InitializeVariableStatusesForWarmStart(solution_state_, + ColIndex(0)); variable_values_.RecomputeBasicVariableValues(); } solve_from_scratch = false; - } else if (only_new_rows) { + } else if (only_change_is_new_rows) { // For the dual-simplex, we also perform a warm start if a couple of // new rows where added. - InitializeVariableStatusesForWarmStart(solution_state_); + InitializeVariableStatusesForWarmStart(solution_state_, ColIndex(0)); // TODO(user): Both the edge norms and the reduced costs do not really // need to be recomputed. We just need to initialize the ones of the @@ -1265,8 +1332,9 @@ void RevisedSimplex::ComputeDirection(ColIndex col) { std::max(direction_infinity_norm_, fabs(direction_[row])); } IF_STATS_ENABLED(ratio_test_stats_.direction_density.Add( - num_rows_ == 0 ? 0.0 : static_cast(direction_non_zero_.size()) / - static_cast(num_rows_.value()))); + num_rows_ == 0 ? 0.0 + : static_cast(direction_non_zero_.size()) / + static_cast(num_rows_.value()))); } Fractional RevisedSimplex::ComputeDirectionError(ColIndex col) { diff --git a/src/glop/revised_simplex.h b/src/glop/revised_simplex.h index 40afff2ee7..d043a7b6b5 100644 --- a/src/glop/revised_simplex.h +++ b/src/glop/revised_simplex.h @@ -102,14 +102,14 @@ #include "base/macros.h" #include "glop/basis_representation.h" #include "glop/dual_edge_norms.h" +#include "glop/entering_variable.h" #include "glop/parameters.pb.h" #include "glop/primal_edge_norms.h" -#include "glop/entering_variable.h" #include "glop/reduced_costs.h" #include "glop/status.h" #include "glop/update_row.h" -#include "glop/variables_info.h" #include "glop/variable_values.h" +#include "glop/variables_info.h" #include "lp_data/lp_data.h" #include "lp_data/lp_print_utils.h" #include "lp_data/lp_types.h" @@ -120,7 +120,7 @@ namespace operations_research { namespace glop { -// Holds the statuses of the all the variables, including slack variables. There +// Holds the statuses of all the variables, including slack variables. There // is no point storing constraint statuses since internally all constraints are // always fixed to zero. // @@ -305,15 +305,25 @@ class RevisedSimplex { VariableStatus leaving_variable_status); // Initializes matrix-related internal data. Returns true if this data was - // unchanged. If not, also sets only_new_rows to true if compared to the - // current matrix, the only difference is that new rows have been added (with - // their corresponding extra slack variables). + // unchanged. If not, also sets only_change_is_new_rows to true if compared + // to the current matrix, the only difference is that new rows have been + // added (with their corresponding extra slack variables). Similarly, sets + // only_change_is_new_cols to true if the only difference is that new columns + // have been added, in which case also sets num_new_cols to the number of + // new columns. bool InitializeMatrixAndTestIfUnchanged(const LinearProgram& lp, - bool* only_new_rows); + bool* only_change_is_new_rows, + bool* only_change_is_new_cols, + ColIndex* num_new_cols); // Initializes bound-related internal data. Returns true if unchanged. bool InitializeBoundsAndTestIfUnchanged(const LinearProgram& lp); + // Checks if the only change to the bounds is the addition of new columns, + // and that the new columns have at least one bound equal to zero. + bool OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero( + const LinearProgram& lp, ColIndex num_new_cols); + // Initializes objective-related internal data. Returns true if unchanged. bool InitializeObjectiveAndTestIfUnchanged(const LinearProgram& lp); @@ -321,7 +331,8 @@ class RevisedSimplex { void InitializeObjectiveLimit(const LinearProgram& lp); // Initializes the variable statuses using a warm-start basis. - void InitializeVariableStatusesForWarmStart(const BasisState& state); + void InitializeVariableStatusesForWarmStart(const BasisState& state, + ColIndex num_new_cols); // Initializes the starting basis. In most cases it starts by the all slack // basis and tries to apply some heuristics to replace fixed variables. diff --git a/src/glop/update_row.cc b/src/glop/update_row.cc index 8419f9dd04..3c89c645a5 100644 --- a/src/glop/update_row.cc +++ b/src/glop/update_row.cc @@ -30,7 +30,6 @@ UpdateRow::UpdateRow(const CompactSparseMatrix& matrix, : matrix_(matrix), transposed_matrix_(transposed_matrix), variables_info_(variables_info), - basis_(basis), basis_factorization_(basis_factorization), unit_row_left_inverse_(), non_zero_position_list_(), diff --git a/src/glop/update_row.h b/src/glop/update_row.h index 9928f3c104..6c901c6835 100644 --- a/src/glop/update_row.h +++ b/src/glop/update_row.h @@ -103,7 +103,6 @@ class UpdateRow { const CompactSparseMatrix& matrix_; const CompactSparseMatrix& transposed_matrix_; const VariablesInfo& variables_info_; - const RowToColMapping& basis_; const BasisFactorization& basis_factorization_; // Left inverse by B of a unit row. Its scalar product with a column 'a' of A @@ -118,7 +117,6 @@ class UpdateRow { DenseRow coefficient_; // Boolean used to avoid recomputing many times the same thing. - bool compute_unit_row_left_inverse_; bool compute_update_row_; // Statistics about this class. diff --git a/src/linear_solver/glop_interface.cc b/src/linear_solver/glop_interface.cc index 6210d35e18..3eb63b4bbb 100644 --- a/src/linear_solver/glop_interface.cc +++ b/src/linear_solver/glop_interface.cc @@ -35,6 +35,7 @@ #include "base/hash.h" #include "glop/lp_solver.h" #include "glop/parameters.pb.h" +#include "linear_solver/glop_utils.h" #include "linear_solver/linear_solver.h" #include "lp_data/lp_data.h" #include "lp_data/lp_types.h" @@ -44,69 +45,6 @@ namespace operations_research { namespace { -MPSolver::ResultStatus TranslateProblemStatus(glop::ProblemStatus status) { - switch (status) { - case glop::ProblemStatus::OPTIMAL: - return MPSolver::OPTIMAL; - case glop::ProblemStatus::PRIMAL_FEASIBLE: - return MPSolver::FEASIBLE; - case glop::ProblemStatus::PRIMAL_INFEASIBLE: // PASS_THROUGH_INTENDED - case glop::ProblemStatus::DUAL_UNBOUNDED: - return MPSolver::INFEASIBLE; - case glop::ProblemStatus::PRIMAL_UNBOUNDED: - return MPSolver::UNBOUNDED; - case glop::ProblemStatus::DUAL_FEASIBLE: // PASS_THROUGH_INTENDED - case glop::ProblemStatus::INIT: - return MPSolver::NOT_SOLVED; - // TODO(user): Glop may return ProblemStatus::DUAL_INFEASIBLE or - // ProblemStatus::INFEASIBLE_OR_UNBOUNDED. - // Unfortunatley, the wrapper does not support this return status at this - // point (even though Cplex and Gurobi have the equivalent). So we convert - // it to MPSolver::ABNORMAL instead. - case glop::ProblemStatus::DUAL_INFEASIBLE: // PASS_THROUGH_INTENDED - case glop::ProblemStatus::INFEASIBLE_OR_UNBOUNDED: // PASS_THROUGH_INTENDED - case glop::ProblemStatus::ABNORMAL: // PASS_THROUGH_INTENDED - case glop::ProblemStatus::IMPRECISE: // PASS_THROUGH_INTENDED - case glop::ProblemStatus::INVALID_PROBLEM: - return MPSolver::ABNORMAL; - } - LOG(DFATAL) << "Invalid glop::ProblemStatus " << status; - return MPSolver::ABNORMAL; -} - -MPSolver::BasisStatus TranslateVariableStatus(glop::VariableStatus status) { - switch (status) { - case glop::VariableStatus::FREE: - return MPSolver::FREE; - case glop::VariableStatus::AT_LOWER_BOUND: - return MPSolver::AT_LOWER_BOUND; - case glop::VariableStatus::AT_UPPER_BOUND: - return MPSolver::AT_UPPER_BOUND; - case glop::VariableStatus::FIXED_VALUE: - return MPSolver::FIXED_VALUE; - case glop::VariableStatus::BASIC: - return MPSolver::BASIC; - } - LOG(DFATAL) << "Unknown variable status: " << status; - return MPSolver::FREE; -} - -MPSolver::BasisStatus TranslateConstraintStatus(glop::ConstraintStatus status) { - switch (status) { - case glop::ConstraintStatus::FREE: - return MPSolver::FREE; - case glop::ConstraintStatus::AT_LOWER_BOUND: - return MPSolver::AT_LOWER_BOUND; - case glop::ConstraintStatus::AT_UPPER_BOUND: - return MPSolver::AT_UPPER_BOUND; - case glop::ConstraintStatus::FIXED_VALUE: - return MPSolver::FIXED_VALUE; - case glop::ConstraintStatus::BASIC: - return MPSolver::BASIC; - } - LOG(DFATAL) << "Unknown constraint status: " << status; - return MPSolver::FREE; -} } // Anonymous namespace class GLOPInterface : public MPSolverInterface { @@ -154,6 +92,10 @@ class GLOPInterface : public MPSolverInterface { void ExtractNewConstraints() override; void ExtractObjective() override; + void SetStartingLpBasis( + const std::vector& variable_statuses, + const std::vector& constraint_statuses) override; + void SetParameters(const MPSolverParameters& param) override; void SetRelativeMipGap(double value) override; void SetPrimalTolerance(double value) override; @@ -186,8 +128,11 @@ GLOPInterface::GLOPInterface(MPSolver* const solver) GLOPInterface::~GLOPInterface() {} MPSolver::ResultStatus GLOPInterface::Solve(const MPSolverParameters& param) { - // Reset extraction as Glop is not incremental yet. - Reset(); + // Re-extract the problem from scratch. We don't support modifying the + // LinearProgram in sync with changes done in the MPSolver. + ResetExtractionInformation(); + linear_program_.Clear(); + interrupt_solver_ = false; ExtractModel(); SetParameters(param); @@ -212,7 +157,7 @@ MPSolver::ResultStatus GLOPInterface::Solve(const MPSolverParameters& param) { // The solution must be marked as synchronized even when no solution exists. sync_status_ = SOLUTION_SYNCHRONIZED; - result_status_ = TranslateProblemStatus(status); + result_status_ = GlopToMPSolverResultStatus(status); objective_value_ = lp_solver_.GetObjectiveValue(); const size_t num_vars = solver_->variables_.size(); @@ -231,7 +176,7 @@ MPSolver::ResultStatus GLOPInterface::Solve(const MPSolverParameters& param) { const glop::VariableStatus variable_status = lp_solver_.variable_statuses()[lp_solver_var_id]; - column_status_.at(var_id) = TranslateVariableStatus(variable_status); + column_status_.at(var_id) = GlopToMPSolverVariableStatus(variable_status); } const size_t num_constraints = solver_->constraints_.size(); @@ -246,7 +191,7 @@ MPSolver::ResultStatus GLOPInterface::Solve(const MPSolverParameters& param) { const glop::ConstraintStatus constraint_status = lp_solver_.constraint_statuses()[lp_solver_ct_id]; - row_status_.at(ct_id) = TranslateConstraintStatus(constraint_status); + row_status_.at(ct_id) = GlopToMPSolverConstraintStatus(constraint_status); } return result_status_; @@ -258,9 +203,9 @@ bool GLOPInterface::InterruptSolve() { } void GLOPInterface::Reset() { - ResetExtractionInformation(); - linear_program_.Clear(); - interrupt_solver_ = false; + // Ignore any incremental info for the next solve. Note that the parameters + // will not be reset as we re-read them on each Solve(). + lp_solver_.Clear(); } void GLOPInterface::SetOptimizationDirection(bool maximize) { @@ -391,6 +336,20 @@ void GLOPInterface::ExtractObjective() { } } +void GLOPInterface::SetStartingLpBasis( + const std::vector& variable_statuses, + const std::vector& constraint_statuses) { + glop::VariableStatusRow glop_variable_statuses; + glop::ConstraintStatusColumn glop_constraint_statuses; + for (const MPSolver::BasisStatus& status : variable_statuses) { + glop_variable_statuses.push_back(MPSolverToGlopVariableStatus(status)); + } + for (const MPSolver::BasisStatus& status : constraint_statuses) { + glop_constraint_statuses.push_back(MPSolverToGlopConstraintStatus(status)); + } + lp_solver_.SetInitialBasis(glop_variable_statuses, glop_constraint_statuses); +} + void GLOPInterface::SetParameters(const MPSolverParameters& param) { parameters_.Clear(); SetCommonParameters(param); diff --git a/src/linear_solver/glop_utils.cc b/src/linear_solver/glop_utils.cc new file mode 100644 index 0000000000..7cc21ca731 --- /dev/null +++ b/src/linear_solver/glop_utils.cc @@ -0,0 +1,116 @@ +// 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 "linear_solver/glop_utils.h" + +namespace operations_research { + +MPSolver::ResultStatus GlopToMPSolverResultStatus(glop::ProblemStatus s) { + switch (s) { + case glop::ProblemStatus::OPTIMAL: + return MPSolver::OPTIMAL; + case glop::ProblemStatus::PRIMAL_FEASIBLE: + return MPSolver::FEASIBLE; + case glop::ProblemStatus::PRIMAL_INFEASIBLE: // PASS_THROUGH_INTENDED + case glop::ProblemStatus::DUAL_UNBOUNDED: + return MPSolver::INFEASIBLE; + case glop::ProblemStatus::PRIMAL_UNBOUNDED: + return MPSolver::UNBOUNDED; + case glop::ProblemStatus::DUAL_FEASIBLE: // PASS_THROUGH_INTENDED + case glop::ProblemStatus::INIT: + return MPSolver::NOT_SOLVED; + // TODO(user): Glop may return ProblemStatus::DUAL_INFEASIBLE or + // ProblemStatus::INFEASIBLE_OR_UNBOUNDED. + // Unfortunatley, the wrapper does not support this return status at this + // point (even though Cplex and Gurobi have the equivalent). So we convert + // it to MPSolver::ABNORMAL instead. + case glop::ProblemStatus::DUAL_INFEASIBLE: // PASS_THROUGH_INTENDED + case glop::ProblemStatus::INFEASIBLE_OR_UNBOUNDED: // PASS_THROUGH_INTENDED + case glop::ProblemStatus::ABNORMAL: // PASS_THROUGH_INTENDED + case glop::ProblemStatus::IMPRECISE: // PASS_THROUGH_INTENDED + case glop::ProblemStatus::INVALID_PROBLEM: + return MPSolver::ABNORMAL; + } + LOG(DFATAL) << "Invalid glop::ProblemStatus " << s; + return MPSolver::ABNORMAL; +} + +MPSolver::BasisStatus GlopToMPSolverVariableStatus(glop::VariableStatus s) { + switch (s) { + case glop::VariableStatus::FREE: + return MPSolver::FREE; + case glop::VariableStatus::AT_LOWER_BOUND: + return MPSolver::AT_LOWER_BOUND; + case glop::VariableStatus::AT_UPPER_BOUND: + return MPSolver::AT_UPPER_BOUND; + case glop::VariableStatus::FIXED_VALUE: + return MPSolver::FIXED_VALUE; + case glop::VariableStatus::BASIC: + return MPSolver::BASIC; + } + LOG(DFATAL) << "Unknown variable status: " << s; + return MPSolver::FREE; +} + +glop::VariableStatus MPSolverToGlopVariableStatus(MPSolver::BasisStatus s) { + switch (s) { + case MPSolver::FREE: + return glop::VariableStatus::FREE; + case MPSolver::AT_LOWER_BOUND: + return glop::VariableStatus::AT_LOWER_BOUND; + case MPSolver::AT_UPPER_BOUND: + return glop::VariableStatus::AT_UPPER_BOUND; + case MPSolver::FIXED_VALUE: + return glop::VariableStatus::FIXED_VALUE; + case MPSolver::BASIC: + return glop::VariableStatus::BASIC; + } + LOG(DFATAL) << "Unknown variable status: " << s; + return glop::VariableStatus::FREE; +} + +MPSolver::BasisStatus GlopToMPSolverConstraintStatus(glop::ConstraintStatus s) { + switch (s) { + case glop::ConstraintStatus::FREE: + return MPSolver::FREE; + case glop::ConstraintStatus::AT_LOWER_BOUND: + return MPSolver::AT_LOWER_BOUND; + case glop::ConstraintStatus::AT_UPPER_BOUND: + return MPSolver::AT_UPPER_BOUND; + case glop::ConstraintStatus::FIXED_VALUE: + return MPSolver::FIXED_VALUE; + case glop::ConstraintStatus::BASIC: + return MPSolver::BASIC; + } + LOG(DFATAL) << "Unknown constraint status: " << s; + return MPSolver::FREE; +} + +glop::ConstraintStatus MPSolverToGlopConstraintStatus(MPSolver::BasisStatus s) { + switch (s) { + case MPSolver::FREE: + return glop::ConstraintStatus::FREE; + case MPSolver::AT_LOWER_BOUND: + return glop::ConstraintStatus::AT_LOWER_BOUND; + case MPSolver::AT_UPPER_BOUND: + return glop::ConstraintStatus::AT_UPPER_BOUND; + case MPSolver::FIXED_VALUE: + return glop::ConstraintStatus::FIXED_VALUE; + case MPSolver::BASIC: + return glop::ConstraintStatus::BASIC; + } + LOG(DFATAL) << "Unknown constraint status: " << s; + return glop::ConstraintStatus::FREE; +} + +} // namespace operations_research diff --git a/src/linear_solver/glop_utils.h b/src/linear_solver/glop_utils.h new file mode 100644 index 0000000000..6027bab789 --- /dev/null +++ b/src/linear_solver/glop_utils.h @@ -0,0 +1,32 @@ +// 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_LINEAR_SOLVER_GLOP_UTILS_H_ +#define OR_TOOLS_LINEAR_SOLVER_GLOP_UTILS_H_ + +#include "linear_solver/linear_solver.h" +#include "lp_data/lp_types.h" + +namespace operations_research { + +MPSolver::ResultStatus GlopToMPSolverResultStatus(glop::ProblemStatus s); + +MPSolver::BasisStatus GlopToMPSolverVariableStatus(glop::VariableStatus s); +glop::VariableStatus MPSolverToGlopVariableStatus(MPSolver::BasisStatus s); + +MPSolver::BasisStatus GlopToMPSolverConstraintStatus(glop::ConstraintStatus s); +glop::ConstraintStatus MPSolverToGlopConstraintStatus(MPSolver::BasisStatus s); + +} // namespace operations_research + +#endif // OR_TOOLS_LINEAR_SOLVER_GLOP_UTILS_H_ diff --git a/src/linear_solver/linear_solver.cc b/src/linear_solver/linear_solver.cc index 6ffb5c4f12..f85bae4b63 100644 --- a/src/linear_solver/linear_solver.cc +++ b/src/linear_solver/linear_solver.cc @@ -752,6 +752,12 @@ void MPSolver::Reset() { interface_->Reset(); } bool MPSolver::InterruptSolve() { return interface_->InterruptSolve(); } +void MPSolver::SetStartingLpBasis( + const std::vector& variable_statuses, + const std::vector& constraint_statuses) { + interface_->SetStartingLpBasis(variable_statuses, constraint_statuses); +} + MPVariable* MPSolver::MakeVar(double lb, double ub, bool integer, const std::string& name) { const int var_index = NumVariables(); diff --git a/src/linear_solver/linear_solver.h b/src/linear_solver/linear_solver.h index e63343e74f..75fb5fe0ef 100644 --- a/src/linear_solver/linear_solver.h +++ b/src/linear_solver/linear_solver.h @@ -459,6 +459,18 @@ class MPSolver { BASIC }; + // Advanced usage: Incrementality. This function takes a starting basis to be + // used in the next LP Solve() call. The statuses of a current solution can be + // retrieved via the basis_status() function of a MPVariable or a + // MPConstraint. + // + // WARNING: With Glop, you should disable presolve when using this because + // this information will not be modified in sync with the presolve and will + // likely not mean much on the presolved problem. + void SetStartingLpBasis( + const std::vector& variable_statuses, + const std::vector& constraint_statuses); + // Infinity. You can use -MPSolver::infinity() for negative infinity. static double infinity() { return std::numeric_limits::infinity(); } @@ -847,9 +859,13 @@ class MPConstraint { // Advanced usage: returns the dual value of the constraint in the // current solution (only available for continuous problems). double dual_value() const; - // Advanced usage: returns the basis status of the slack variable - // associated with the constraint (only available for continuous - // problems). + + // Advanced usage: returns the basis status of the constraint (only available + // for continuous problems). Note that if a constraint "linear_expression in + // [lb, ub]" is transformed into "linear_expression + slack = 0" with slack in + // [-ub, -lb], then this status is the same as the status of the slack + // variable with AT_UPPER_BOUND and AT_LOWER_BOUND swapped. + // // @see MPSolver::BasisStatus. MPSolver::BasisStatus basis_status() const; @@ -1232,6 +1248,13 @@ class MPSolverInterface { // problems and only implemented in GLPK. virtual double ComputeExactConditionNumber() const; + // See MPSolver::SetStartingLpBasis(). + virtual void SetStartingLpBasis( + const std::vector& variable_statuses, + const std::vector& constraint_statuses) { + LOG(FATAL) << "Not supported by this solver."; + } + virtual bool InterruptSolve() { return false; } friend class MPSolver;