From 790c94745344894cd2aa59f5c5f7369e9824530e Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Fri, 23 Jul 2021 21:54:59 +0200 Subject: [PATCH] remove forced slacks in glop --- ortools/glop/lp_solver.cc | 5 +- ortools/glop/preprocessor.cc | 2 - ortools/glop/revised_simplex.cc | 117 ++++++++++++++++++++++---------- ortools/glop/revised_simplex.h | 16 +++-- ortools/glop/variables_info.cc | 41 +++++++++++ ortools/glop/variables_info.h | 7 ++ ortools/lp_data/sparse.cc | 27 ++++++++ ortools/lp_data/sparse.h | 4 ++ 8 files changed, 172 insertions(+), 47 deletions(-) diff --git a/ortools/glop/lp_solver.cc b/ortools/glop/lp_solver.cc index 563d4e342d..2837a1ecc5 100644 --- a/ortools/glop/lp_solver.cc +++ b/ortools/glop/lp_solver.cc @@ -581,8 +581,9 @@ void LPSolver::RunRevisedSimplexIfNeeded(ProblemSolution* solution, num_revised_simplex_iterations_ = revised_simplex_->GetNumberOfIterations(); solution->status = revised_simplex_->GetProblemStatus(); - const ColIndex num_cols = revised_simplex_->GetProblemNumCols(); - DCHECK_EQ(solution->primal_values.size(), num_cols); + // Make sure we do not copy the slacks added by revised_simplex_. + const ColIndex num_cols = solution->primal_values.size(); + DCHECK_LE(num_cols, revised_simplex_->GetProblemNumCols()); for (ColIndex col(0); col < num_cols; ++col) { solution->primal_values[col] = revised_simplex_->GetVariableValue(col); solution->variable_statuses[col] = diff --git a/ortools/glop/preprocessor.cc b/ortools/glop/preprocessor.cc index c954806080..9b7a979d1f 100644 --- a/ortools/glop/preprocessor.cc +++ b/ortools/glop/preprocessor.cc @@ -123,8 +123,6 @@ bool MainLpPreprocessor::Run(LinearProgram* lp) { // The scaling is controled by use_scaling, not use_preprocessing. RUN_PREPROCESSOR(ScalingPreprocessor); - // This one must always run. It is needed by the revised simplex code. - RUN_PREPROCESSOR(AddSlackVariablesPreprocessor); return !preprocessors_.empty(); } diff --git a/ortools/glop/revised_simplex.cc b/ortools/glop/revised_simplex.cc index 015e28a213..e33dbfc82a 100644 --- a/ortools/glop/revised_simplex.cc +++ b/ortools/glop/revised_simplex.cc @@ -138,10 +138,6 @@ Status RevisedSimplex::Solve(const LinearProgram& lp, TimeLimit* time_limit) { SCOPED_TIME_STAT(&function_stats_); DCHECK(lp.IsCleanedUp()); GLOP_RETURN_ERROR_IF_NULL(time_limit); - if (!lp.IsInEquationForm()) { - return Status(Status::ERROR_INVALID_PROBLEM, - "The problem is not in the equations form."); - } Cleanup update_deterministic_time_on_return( [this, time_limit]() { AdvanceDeterministicTime(time_limit); }); @@ -806,28 +802,30 @@ void RevisedSimplex::UseSingletonColumnInInitialBasis(RowToColMapping* basis) { } bool RevisedSimplex::InitializeMatrixAndTestIfUnchanged( - const LinearProgram& lp, bool* only_change_is_new_rows, - bool* only_change_is_new_cols, ColIndex* num_new_cols) { + const LinearProgram& lp, bool lp_is_in_equation_form, + bool* only_change_is_new_rows, bool* only_change_is_new_cols, + ColIndex* num_new_cols) { SCOPED_TIME_STAT(&function_stats_); DCHECK(only_change_is_new_rows != nullptr); DCHECK(only_change_is_new_cols != nullptr); DCHECK(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()); - DCHECK_EQ(lp.num_variables(), - lp.GetFirstSlackVariable() + RowToColIndex(lp.num_constraints())); - DCHECK(IsRightMostSquareMatrixIdentity(lp.GetSparseMatrix())); + // This works whether the lp is in equation form (with slack) or not. const bool old_part_of_matrix_is_unchanged = AreFirstColumnsAndRowsExactlyEquals( num_rows_, first_slack_col_, lp.GetSparseMatrix(), compact_matrix_); + // This is the only adaptation we need for the test below. + const ColIndex lp_first_slack = + lp_is_in_equation_form ? lp.GetFirstSlackVariable() : lp.num_variables(); + // 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_) { + lp_first_slack == first_slack_col_) { return true; } @@ -835,38 +833,45 @@ bool RevisedSimplex::InitializeMatrixAndTestIfUnchanged( // 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_; + lp_first_slack == 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); + lp_first_slack > first_slack_col_; + *num_new_cols = *only_change_is_new_cols ? lp_first_slack - first_slack_col_ + : ColIndex(0); // Initialize first_slack_. - first_slack_col_ = lp.GetFirstSlackVariable(); + first_slack_col_ = lp_first_slack; // Initialize the new dimensions. num_rows_ = lp.num_constraints(); - num_cols_ = lp.num_variables(); + num_cols_ = lp_first_slack + RowToColIndex(lp.num_constraints()); // Populate compact_matrix_ and transposed_matrix_ if needed. Note that we // already added all the slack variables at this point, so matrix_ will not // change anymore. - // TODO(user): This can be sped up by removing the MatrixView. - compact_matrix_.PopulateFromMatrixView(MatrixView(lp.GetSparseMatrix())); + if (lp_is_in_equation_form) { + // TODO(user): This can be sped up by removing the MatrixView, but then + // this path will likely go away. + compact_matrix_.PopulateFromMatrixView(MatrixView(lp.GetSparseMatrix())); + } else { + compact_matrix_.PopulateFromSparseMatrixAndAddSlacks(lp.GetSparseMatrix()); + } if (parameters_.use_transposed_matrix()) { transposed_matrix_.PopulateFromTranspose(compact_matrix_); } return false; } +// Preconditions: This should only be called if there are only new variable +// in the lp. bool RevisedSimplex::OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero( - const LinearProgram& lp, ColIndex num_new_cols) { + const LinearProgram& lp, bool lp_is_in_equation_form, + 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); @@ -879,6 +884,7 @@ bool RevisedSimplex::OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero( 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 && @@ -886,11 +892,25 @@ bool RevisedSimplex::OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero( return false; } } + // Check that the slack bounds are unchanged. - for (ColIndex col(first_slack_col_); col < num_cols_; ++col) { - if (lower_bounds[col - num_new_cols] != lp.variable_lower_bounds()[col] || - upper_bounds[col - num_new_cols] != lp.variable_upper_bounds()[col]) { - return false; + if (lp_is_in_equation_form) { + for (ColIndex col(first_slack_col_); col < num_cols_; ++col) { + if (lower_bounds[col - num_new_cols] != lp.variable_lower_bounds()[col] || + upper_bounds[col - num_new_cols] != lp.variable_upper_bounds()[col]) { + return false; + } + } + } else { + DCHECK_EQ(num_rows_, lp.num_constraints()); + for (RowIndex row(0); row < num_rows_; ++row) { + const ColIndex col = first_slack_col_ + RowToColIndex(row); + if (lower_bounds[col - num_new_cols] != + -lp.constraint_upper_bounds()[row] || + upper_bounds[col - num_new_cols] != + -lp.constraint_lower_bounds()[row]) { + return false; + } } } return true; @@ -902,31 +922,40 @@ bool RevisedSimplex::InitializeObjectiveAndTestIfUnchanged( bool objective_is_unchanged = true; objective_.resize(num_cols_, 0.0); - DCHECK_EQ(num_cols_, lp.num_variables()); + + // This function work whether the lp is in equation form (with slack) or + // without, since the objective of the slacks are always zero. + DCHECK_GE(num_cols_, lp.num_variables()); + for (ColIndex col(lp.num_variables()); col < num_cols_; ++col) { + if (objective_[col] != 0.0) { + objective_is_unchanged = false; + objective_[col] = 0.0; + } + } + if (lp.IsMaximizationProblem()) { // Note that we use the minimization version of the objective internally. for (ColIndex col(0); col < lp.num_variables(); ++col) { const Fractional coeff = -lp.objective_coefficients()[col]; if (objective_[col] != coeff) { objective_is_unchanged = false; + objective_[col] = coeff; } - objective_[col] = coeff; } objective_offset_ = -lp.objective_offset(); objective_scaling_factor_ = -lp.objective_scaling_factor(); } else { for (ColIndex col(0); col < lp.num_variables(); ++col) { - if (objective_[col] != lp.objective_coefficients()[col]) { + const Fractional coeff = lp.objective_coefficients()[col]; + if (objective_[col] != coeff) { objective_is_unchanged = false; - break; + objective_[col] = coeff; } } - if (!objective_is_unchanged) { - objective_ = lp.objective_coefficients(); - } objective_offset_ = lp.objective_offset(); objective_scaling_factor_ = lp.objective_scaling_factor(); } + return objective_is_unchanged; } @@ -1176,6 +1205,14 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) { parameters_ = initial_parameters_; PropagateParameters(); + // We accept both kind of input. + // + // TODO(user): Ideally there should be no need to ever put the slack in the + // LinearProgram. That take extra memory (one big SparseColumn per slack) and + // just add visible overhead in incremental solve when one wants to add/remove + // constraints. But for historical reason, we handle both for now. + const bool lp_is_in_equation_form = lp.IsInEquationForm(); + // Calling InitializeMatrixAndTestIfUnchanged() first is important because // this is where num_rows_ and num_cols_ are computed. // @@ -1188,13 +1225,15 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) { bool only_new_bounds = false; if (solution_state_.IsEmpty() || !notify_that_matrix_is_unchanged_) { matrix_is_unchanged = InitializeMatrixAndTestIfUnchanged( - lp, &only_change_is_new_rows, &only_change_is_new_cols, &num_new_cols); + lp, lp_is_in_equation_form, &only_change_is_new_rows, + &only_change_is_new_cols, &num_new_cols); only_new_bounds = only_change_is_new_cols && num_new_cols > 0 && OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero( - lp, num_new_cols); + lp, lp_is_in_equation_form, num_new_cols); } else if (DEBUG_MODE) { CHECK(InitializeMatrixAndTestIfUnchanged( - lp, &only_change_is_new_rows, &only_change_is_new_cols, &num_new_cols)); + lp, lp_is_in_equation_form, &only_change_is_new_rows, + &only_change_is_new_cols, &num_new_cols)); } notify_that_matrix_is_unchanged_ = false; @@ -1202,8 +1241,12 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) { const bool objective_is_unchanged = InitializeObjectiveAndTestIfUnchanged(lp); const bool bounds_are_unchanged = - variables_info_.LoadBoundsAndReturnTrueIfUnchanged( - lp.variable_lower_bounds(), lp.variable_upper_bounds()); + lp_is_in_equation_form + ? variables_info_.LoadBoundsAndReturnTrueIfUnchanged( + lp.variable_lower_bounds(), lp.variable_upper_bounds()) + : variables_info_.LoadBoundsAndReturnTrueIfUnchanged( + lp.variable_lower_bounds(), lp.variable_upper_bounds(), + lp.constraint_lower_bounds(), lp.constraint_upper_bounds()); // If parameters_.allow_simplex_algorithm_change() is true and we already have // a primal (resp. dual) feasible solution, then we use the primal (resp. diff --git a/ortools/glop/revised_simplex.h b/ortools/glop/revised_simplex.h index 2258e36d32..0042a55f6f 100644 --- a/ortools/glop/revised_simplex.h +++ b/ortools/glop/revised_simplex.h @@ -131,11 +131,13 @@ class RevisedSimplex { // Solves the given linear program. // - // Expects that the linear program is in the equations form Ax = 0 created by - // LinearProgram::AddSlackVariablesForAllRows, i.e. the rightmost square - // submatrix of A is an identity matrix, all its columns have been marked as - // slack variables, and the bounds of all constraints have been set to [0, 0]. - // Returns ERROR_INVALID_PROBLEM, if these assumptions are violated. + // We accept two forms of LinearProgram: + // - The lp can be in the equations form Ax = 0 created by + // LinearProgram::AddSlackVariablesForAllRows(), i.e. the rightmost square + // submatrix of A is an identity matrix, all its columns have been marked as + // slack variables, and the bounds of all constraints have been set to 0. + // - If not, we will convert it internally while copying it to the internal + // structure used. // // By default, the algorithm tries to exploit the computation done during the // last Solve() call. It will analyze the difference of the new linear program @@ -323,6 +325,7 @@ class RevisedSimplex { // have been added, in which case also sets num_new_cols to the number of // new columns. bool InitializeMatrixAndTestIfUnchanged(const LinearProgram& lp, + bool lp_is_in_equation_form, bool* only_change_is_new_rows, bool* only_change_is_new_cols, ColIndex* num_new_cols); @@ -330,7 +333,8 @@ class RevisedSimplex { // 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); + const LinearProgram& lp, bool lp_is_in_equation_form, + ColIndex num_new_cols); // Initializes objective-related internal data. Returns true if unchanged. bool InitializeObjectiveAndTestIfUnchanged(const LinearProgram& lp); diff --git a/ortools/glop/variables_info.cc b/ortools/glop/variables_info.cc index f9626d1a97..0a64224c5d 100644 --- a/ortools/glop/variables_info.cc +++ b/ortools/glop/variables_info.cc @@ -39,6 +39,47 @@ bool VariablesInfo::LoadBoundsAndReturnTrueIfUnchanged( return false; } +bool VariablesInfo::LoadBoundsAndReturnTrueIfUnchanged( + const DenseRow& variable_lower_bounds, + const DenseRow& variable_upper_bounds, + const DenseColumn& constraint_lower_bounds, + const DenseColumn& constraint_upper_bounds) { + const ColIndex num_cols = matrix_.num_cols(); + const ColIndex num_variables = variable_upper_bounds.size(); + const RowIndex num_rows = constraint_lower_bounds.size(); + + bool is_unchanged = (num_cols == lower_bounds_.size()); + DCHECK_EQ(num_cols, num_variables + RowToColIndex(num_rows)); + lower_bounds_.resize(num_cols, 0.0); + upper_bounds_.resize(num_cols, 0.0); + variable_type_.resize(num_cols, VariableType::FIXED_VARIABLE); + + // Copy bounds of the variables. + for (ColIndex col(0); col < num_variables; ++col) { + if (lower_bounds_[col] != variable_lower_bounds[col] || + upper_bounds_[col] != variable_upper_bounds[col]) { + lower_bounds_[col] = variable_lower_bounds[col]; + upper_bounds_[col] = variable_upper_bounds[col]; + is_unchanged = false; + variable_type_[col] = ComputeVariableType(col); + } + } + + // Copy bounds of the slack. + for (RowIndex row(0); row < num_rows; ++row) { + const ColIndex col = num_variables + RowToColIndex(row); + if (lower_bounds_[col] != -constraint_upper_bounds[row] || + upper_bounds_[col] != -constraint_lower_bounds[row]) { + lower_bounds_[col] = -constraint_upper_bounds[row]; + upper_bounds_[col] = -constraint_lower_bounds[row]; + is_unchanged = false; + variable_type_[col] = ComputeVariableType(col); + } + } + + return is_unchanged; +} + void VariablesInfo::ResetStatusInfo() { const ColIndex num_cols = matrix_.num_cols(); DCHECK_EQ(num_cols, lower_bounds_.size()); diff --git a/ortools/glop/variables_info.h b/ortools/glop/variables_info.h index 411671c214..5bef42c465 100644 --- a/ortools/glop/variables_info.h +++ b/ortools/glop/variables_info.h @@ -66,6 +66,13 @@ class VariablesInfo { bool LoadBoundsAndReturnTrueIfUnchanged(const DenseRow& new_lower_bounds, const DenseRow& new_upper_bounds); + // Same for an LP not in equation form. + bool LoadBoundsAndReturnTrueIfUnchanged( + const DenseRow& variable_lower_bounds, + const DenseRow& variable_upper_bounds, + const DenseColumn& constraint_lower_bounds, + const DenseColumn& constraint_upper_bounds); + // Initializes the status according to the given BasisState. Incompatible // statuses will be corrected, and we transfrom the state correctly if new // columns / rows were added. Note however that one will need to update the diff --git a/ortools/lp_data/sparse.cc b/ortools/lp_data/sparse.cc index 670e048abf..ef9e17d633 100644 --- a/ortools/lp_data/sparse.cc +++ b/ortools/lp_data/sparse.cc @@ -453,6 +453,33 @@ void CompactSparseMatrix::PopulateFromMatrixView(const MatrixView& input) { starts_[input.num_cols()] = index; } +void CompactSparseMatrix::PopulateFromSparseMatrixAndAddSlacks( + const SparseMatrix& input) { + num_cols_ = input.num_cols() + RowToColIndex(input.num_rows()); + num_rows_ = input.num_rows(); + const EntryIndex num_entries = + input.num_entries() + EntryIndex(num_rows_.value()); + starts_.assign(num_cols_ + 1, EntryIndex(0)); + coefficients_.assign(num_entries, 0.0); + rows_.assign(num_entries, RowIndex(0)); + EntryIndex index(0); + for (ColIndex col(0); col < input.num_cols(); ++col) { + starts_[col] = index; + for (const SparseColumn::Entry e : input.column(col)) { + coefficients_[index] = e.coefficient(); + rows_[index] = e.row(); + ++index; + } + } + for (RowIndex row(0); row < num_rows_; ++row) { + starts_[input.num_cols() + RowToColIndex(row)] = index; + coefficients_[index] = 1.0; + rows_[index] = row; + ++index; + } + starts_[num_cols_] = index; +} + void CompactSparseMatrix::PopulateFromTranspose( const CompactSparseMatrix& input) { num_cols_ = RowToColIndex(input.num_rows()); diff --git a/ortools/lp_data/sparse.h b/ortools/lp_data/sparse.h index 6e9847cf8e..55abb72af5 100644 --- a/ortools/lp_data/sparse.h +++ b/ortools/lp_data/sparse.h @@ -301,6 +301,10 @@ class CompactSparseMatrix { // each column is preserved. void PopulateFromMatrixView(const MatrixView& input); + // Creates a CompactSparseMatrix by copying the input and adding an identity + // matrix to the left of it. + void PopulateFromSparseMatrixAndAddSlacks(const SparseMatrix& input); + // Creates a CompactSparseMatrix from the transpose of the given // CompactSparseMatrix. Note that the entries in each columns will be ordered // by row indices.