From c5e9d022a704c22df655cbb3f0323afddee3baf7 Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Fri, 5 Jul 2019 09:32:17 +0200 Subject: [PATCH] cleaning --- ortools/glop/basis_representation.cc | 2 +- ortools/glop/basis_representation.h | 2 +- ortools/glop/initial_basis.cc | 37 ++++++++++++------------ ortools/glop/initial_basis.h | 7 +++-- ortools/glop/lu_factorization.cc | 2 +- ortools/glop/lu_factorization.h | 2 +- ortools/glop/primal_edge_norms.cc | 22 +++++++------- ortools/glop/primal_edge_norms.h | 4 +-- ortools/glop/revised_simplex.cc | 43 ++++++++++++---------------- ortools/glop/revised_simplex.h | 10 +------ ortools/lp_data/lp_utils.cc | 34 +++++++++++++++++----- ortools/lp_data/lp_utils.h | 8 ++++-- ortools/lp_data/sparse_column.h | 16 +++++++++++ 13 files changed, 105 insertions(+), 84 deletions(-) diff --git a/ortools/glop/basis_representation.cc b/ortools/glop/basis_representation.cc index 0f31b7991d..a54e0c43eb 100644 --- a/ortools/glop/basis_representation.cc +++ b/ortools/glop/basis_representation.cc @@ -466,7 +466,7 @@ void BasisFactorization::RightSolveForProblemColumn(ColIndex col, } Fractional BasisFactorization::RightSolveSquaredNorm( - const SparseColumn& a) const { + const ColumnView& a) const { SCOPED_TIME_STAT(&stats_); DCHECK(IsRefactorized()); BumpDeterministicTimeForSolve(a.num_entries().value()); diff --git a/ortools/glop/basis_representation.h b/ortools/glop/basis_representation.h index 284c13c0aa..fdd809de1a 100644 --- a/ortools/glop/basis_representation.h +++ b/ortools/glop/basis_representation.h @@ -235,7 +235,7 @@ class BasisFactorization { // Returns the norm of B^{-1}.a, this is a specific function because // it is a bit faster and it avoids polluting the stats of RightSolve(). // It can be called only when IsRefactorized() is true. - Fractional RightSolveSquaredNorm(const SparseColumn& a) const; + Fractional RightSolveSquaredNorm(const ColumnView& a) const; // Returns the norm of (B^T)^{-1}.e_row where e is an unit vector. // This is a bit faster and avoids polluting the stats of LeftSolve(). diff --git a/ortools/glop/initial_basis.cc b/ortools/glop/initial_basis.cc index ce7e7c3966..a1b80f4536 100644 --- a/ortools/glop/initial_basis.cc +++ b/ortools/glop/initial_basis.cc @@ -21,14 +21,15 @@ namespace operations_research { namespace glop { -InitialBasis::InitialBasis(const MatrixView& matrix, const DenseRow& objective, +InitialBasis::InitialBasis(const CompactSparseMatrix& compact_matrix, + const DenseRow& objective, const DenseRow& lower_bound, const DenseRow& upper_bound, const VariableTypeRow& variable_type) : max_scaled_abs_cost_(0.0), bixby_column_comparator_(*this), triangular_column_comparator_(*this), - matrix_(matrix), + compact_matrix_(compact_matrix), objective_(objective), lower_bound_(lower_bound), upper_bound_(upper_bound), @@ -38,7 +39,7 @@ void InitialBasis::CompleteBixbyBasis(ColIndex num_cols, RowToColMapping* basis) { // Initialize can_be_replaced ('I' in Bixby's paper) and has_zero_coefficient // ('r' in Bixby's paper). - const RowIndex num_rows = matrix_.num_rows(); + const RowIndex num_rows = compact_matrix_.num_rows(); DenseBooleanColumn can_be_replaced(num_rows, false); DenseBooleanColumn has_zero_coefficient(num_rows, false); DCHECK_EQ(num_rows, basis->size()); @@ -51,7 +52,7 @@ void InitialBasis::CompleteBixbyBasis(ColIndex num_cols, } // This is 'v' in Bixby's paper. - DenseColumn scaled_diagonal_abs(matrix_.num_rows(), kInfinity); + DenseColumn scaled_diagonal_abs(compact_matrix_.num_rows(), kInfinity); // Compute a list of candidate indices and sort them using the heuristic // described in Bixby's paper. @@ -63,7 +64,7 @@ void InitialBasis::CompleteBixbyBasis(ColIndex num_cols, for (int i = 0; i < candidates.size(); ++i) { bool enter_basis = false; const ColIndex candidate_col_index = candidates[i]; - const SparseColumn& candidate_col = matrix_.column(candidate_col_index); + const auto& candidate_col = compact_matrix_.column(candidate_col_index); // Bixby's heuristic only works with scaled columns. This should be the // case by default since we only use this when the matrix is scaled, but @@ -120,7 +121,7 @@ template void InitialBasis::CompleteTriangularBasis(ColIndex num_cols, RowToColMapping* basis) { // Initialize can_be_replaced. - const RowIndex num_rows = matrix_.num_rows(); + const RowIndex num_rows = compact_matrix_.num_rows(); DenseBooleanColumn can_be_replaced(num_rows, false); DCHECK_EQ(num_rows, basis->size()); basis->resize(num_rows, kInvalidCol); @@ -135,7 +136,7 @@ void InitialBasis::CompleteTriangularBasis(ColIndex num_cols, residual_pattern.Reset(num_rows, num_cols); for (ColIndex col(0); col < num_cols; ++col) { if (only_allow_zero_cost_column && objective_[col] != 0.0) continue; - for (const SparseColumn::Entry e : matrix_.column(col)) { + for (const SparseColumn::Entry e : compact_matrix_.column(col)) { if (can_be_replaced[e.row()]) { residual_pattern.AddEntry(e.row(), col); } @@ -173,7 +174,7 @@ void InitialBasis::CompleteTriangularBasis(ColIndex num_cols, RowIndex row(kInvalidRow); Fractional coeff = 0.0; Fractional max_magnitude = 0.0; - for (const SparseColumn::Entry e : matrix_.column(candidate)) { + for (const SparseColumn::Entry e : compact_matrix_.column(candidate)) { max_magnitude = std::max(max_magnitude, std::abs(e.coefficient())); if (can_be_replaced[e.row()]) { row = e.row(); @@ -218,8 +219,8 @@ int InitialBasis::GetMarosPriority(ColIndex col) const { int InitialBasis::GetMarosPriority(RowIndex row) const { // Priority values for rows are equal to // 3 - row priority values as defined in Maros's book - ColIndex slack_index(RowToColIndex(row) + matrix_.num_cols() - - RowToColIndex(matrix_.num_rows())); + ColIndex slack_index(RowToColIndex(row) + compact_matrix_.num_cols() - + RowToColIndex(compact_matrix_.num_rows())); return GetMarosPriority(slack_index); } @@ -229,7 +230,7 @@ void InitialBasis::GetMarosBasis(ColIndex num_cols, RowToColMapping* basis) { VLOG(1) << "Starting Maros crash procedure."; // Initialize basis to the all-slack basis. - const RowIndex num_rows = matrix_.num_rows(); + const RowIndex num_rows = compact_matrix_.num_rows(); const ColIndex first_slack = num_cols - RowToColIndex(num_rows); DCHECK_EQ(num_rows, basis->size()); basis->resize(num_rows); @@ -255,7 +256,7 @@ void InitialBasis::GetMarosBasis(ColIndex num_cols, RowToColMapping* basis) { MatrixNonZeroPattern residual_pattern; residual_pattern.Reset(num_rows, num_cols); for (ColIndex col(0); col < first_slack; ++col) { - for (const SparseColumn::Entry e : matrix_.column(col)) { + for (const SparseColumn::Entry e : compact_matrix_.column(col)) { if (available[RowToColIndex(e.row())] && available[col]) { residual_pattern.AddEntry(e.row(), col); } @@ -300,7 +301,7 @@ void InitialBasis::GetMarosBasis(ColIndex num_cols, RowToColMapping* basis) { // Make sure that the pivotal entry is not too small in magnitude. Fractional max_magnitude = 0; pivot_absolute_value = 0.0; - const SparseColumn& column_values = matrix_.column(col); + const auto& column_values = compact_matrix_.column(col); for (const SparseColumn::Entry e : column_values) { const Fractional absolute_value = std::fabs(e.coefficient()); if (e.row() == max_rpf_row) pivot_absolute_value = absolute_value; @@ -355,7 +356,7 @@ void InitialBasis::ComputeCandidates(ColIndex num_cols, max_scaled_abs_cost_ = 0.0; for (ColIndex col(0); col < num_cols; ++col) { if (variable_type_[col] != VariableType::FIXED_VARIABLE && - matrix_.column(col).num_entries() > 0) { + compact_matrix_.column(col).num_entries() > 0) { candidates->push_back(col); max_scaled_abs_cost_ = std::max(max_scaled_abs_cost_, std::abs(objective_[col])); @@ -427,10 +428,10 @@ bool InitialBasis::TriangularColumnComparator::operator()( // // TODO(user): Experiments more with this comparator or the // BixbyColumnComparator. - if (initial_basis_.matrix_.column(col_a).num_entries() != - initial_basis_.matrix_.column(col_b).num_entries()) { - return initial_basis_.matrix_.column(col_a).num_entries() > - initial_basis_.matrix_.column(col_b).num_entries(); + if (initial_basis_.compact_matrix_.column(col_a).num_entries() != + initial_basis_.compact_matrix_.column(col_b).num_entries()) { + return initial_basis_.compact_matrix_.column(col_a).num_entries() > + initial_basis_.compact_matrix_.column(col_b).num_entries(); } return initial_basis_.GetColumnPenalty(col_a) > initial_basis_.GetColumnPenalty(col_b); diff --git a/ortools/glop/initial_basis.h b/ortools/glop/initial_basis.h index 710bfc207a..a3587d7994 100644 --- a/ortools/glop/initial_basis.h +++ b/ortools/glop/initial_basis.h @@ -44,8 +44,9 @@ namespace glop { class InitialBasis { public: // Takes references to the linear program data we need. - InitialBasis(const MatrixView& matrix, const DenseRow& objective, - const DenseRow& lower_bound, const DenseRow& upper_bound, + InitialBasis(const CompactSparseMatrix& compact_matrix, + const DenseRow& objective, const DenseRow& lower_bound, + const DenseRow& upper_bound, const VariableTypeRow& variable_type); // Completes the entries of the given basis that are equal to kInvalidCol with @@ -119,7 +120,7 @@ class InitialBasis { const InitialBasis& initial_basis_; } triangular_column_comparator_; - const MatrixView& matrix_; + const CompactSparseMatrix& compact_matrix_; const DenseRow& objective_; const DenseRow& lower_bound_; const DenseRow& upper_bound_; diff --git a/ortools/glop/lu_factorization.cc b/ortools/glop/lu_factorization.cc index a44e90f807..d6e11d9dfc 100644 --- a/ortools/glop/lu_factorization.cc +++ b/ortools/glop/lu_factorization.cc @@ -105,7 +105,7 @@ Fractional ComputeSquaredNormAndResetToZero( } } // namespace -Fractional LuFactorization::RightSolveSquaredNorm(const SparseColumn& a) const { +Fractional LuFactorization::RightSolveSquaredNorm(const ColumnView& a) const { SCOPED_TIME_STAT(&stats_); if (is_identity_factorization_) return SquaredNorm(a); diff --git a/ortools/glop/lu_factorization.h b/ortools/glop/lu_factorization.h index a6176a49a3..b6401a3903 100644 --- a/ortools/glop/lu_factorization.h +++ b/ortools/glop/lu_factorization.h @@ -125,7 +125,7 @@ class LuFactorization { const SparseColumn& GetColumnOfU(ColIndex col) const; // Returns the norm of B^{-1}.a - Fractional RightSolveSquaredNorm(const SparseColumn& a) const; + Fractional RightSolveSquaredNorm(const ColumnView& a) const; // Returns the norm of (B^T)^{-1}.e_row where e is an unit vector. Fractional DualEdgeSquaredNorm(RowIndex row) const; diff --git a/ortools/glop/primal_edge_norms.cc b/ortools/glop/primal_edge_norms.cc index 299d8a4987..79d2d97970 100644 --- a/ortools/glop/primal_edge_norms.cc +++ b/ortools/glop/primal_edge_norms.cc @@ -19,12 +19,10 @@ namespace operations_research { namespace glop { -PrimalEdgeNorms::PrimalEdgeNorms(const MatrixView& matrix, - const CompactSparseMatrix& compact_matrix, +PrimalEdgeNorms::PrimalEdgeNorms(const CompactSparseMatrix& compact_matrix, const VariablesInfo& variables_info, const BasisFactorization& basis_factorization) - : matrix_(matrix), - compact_matrix_(compact_matrix), + : compact_matrix_(compact_matrix), variables_info_(variables_info), basis_factorization_(basis_factorization), stats_(), @@ -115,10 +113,10 @@ void PrimalEdgeNorms::UpdateBeforeBasisPivot(ColIndex entering_col, void PrimalEdgeNorms::ComputeMatrixColumnNorms() { SCOPED_TIME_STAT(&stats_); - 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(); + matrix_column_norms_.resize(compact_matrix_.num_cols(), 0.0); + for (ColIndex col(0); col < compact_matrix_.num_cols(); ++col) { + matrix_column_norms_[col] = sqrt(SquaredNorm(compact_matrix_.column(col))); + num_operations_ += compact_matrix_.column(col).num_entries().value(); } } @@ -128,12 +126,12 @@ void PrimalEdgeNorms::ComputeEdgeSquaredNorms() { // Since we will do a lot of inversions, it is better to be as efficient and // precise as possible by refactorizing the basis. DCHECK(basis_factorization_.IsRefactorized()); - edge_squared_norms_.resize(matrix_.num_cols(), 0.0); + edge_squared_norms_.resize(compact_matrix_.num_cols(), 0.0); for (const ColIndex col : variables_info_.GetIsRelevantBitRow()) { // Note the +1.0 in the squared norm for the component of the edge on the // 'entering_col'. - edge_squared_norms_[col] = - 1.0 + basis_factorization_.RightSolveSquaredNorm(matrix_.column(col)); + edge_squared_norms_[col] = 1.0 + basis_factorization_.RightSolveSquaredNorm( + compact_matrix_.column(col)); } recompute_edge_squared_norms_ = false; } @@ -256,7 +254,7 @@ void PrimalEdgeNorms::ResetDevexWeights() { if (parameters_.initialize_devex_with_column_norms()) { devex_weights_ = GetMatrixColumnNorms(); } else { - devex_weights_.assign(matrix_.num_cols(), 1.0); + devex_weights_.assign(compact_matrix_.num_cols(), 1.0); } num_devex_updates_since_reset_ = 0; reset_devex_weights_ = false; diff --git a/ortools/glop/primal_edge_norms.h b/ortools/glop/primal_edge_norms.h index a5dfad5e06..2197bb862d 100644 --- a/ortools/glop/primal_edge_norms.h +++ b/ortools/glop/primal_edge_norms.h @@ -55,8 +55,7 @@ class PrimalEdgeNorms { // Takes references to the linear program data we need. Note that we assume // that the matrix will never change in our back, but the other references are // supposed to reflect the correct state. - PrimalEdgeNorms(const MatrixView& matrix, - const CompactSparseMatrix& compact_matrix, + PrimalEdgeNorms(const CompactSparseMatrix& compact_matrix, const VariablesInfo& variables_info, const BasisFactorization& basis_factorization); @@ -162,7 +161,6 @@ class PrimalEdgeNorms { const UpdateRow& update_row); // Problem data that should be updated from outside. - const MatrixView& matrix_; const CompactSparseMatrix& compact_matrix_; const VariablesInfo& variables_info_; const BasisFactorization& basis_factorization_; diff --git a/ortools/glop/revised_simplex.cc b/ortools/glop/revised_simplex.cc index d66cc70eef..9e9f5971ed 100644 --- a/ortools/glop/revised_simplex.cc +++ b/ortools/glop/revised_simplex.cc @@ -90,7 +90,7 @@ RevisedSimplex::RevisedSimplex() variable_values_(parameters_, compact_matrix_, basis_, variables_info_, basis_factorization_), dual_edge_norms_(basis_factorization_), - primal_edge_norms_(matrix_with_slack_, compact_matrix_, variables_info_, + primal_edge_norms_(compact_matrix_, variables_info_, basis_factorization_), update_row_(compact_matrix_, transposed_matrix_, variables_info_, basis_, basis_factorization_), @@ -177,9 +177,9 @@ Status RevisedSimplex::Solve(const LinearProgram& lp, TimeLimit* time_limit) { const bool use_dual = parameters_.use_dual_simplex(); VLOG(1) << "------ " << (use_dual ? "Dual simplex." : "Primal simplex."); - VLOG(1) << "The matrix has " << matrix_with_slack_.num_rows() << " rows, " - << matrix_with_slack_.num_cols() << " columns, " - << matrix_with_slack_.num_entries() << " entries."; + VLOG(1) << "The matrix has " << compact_matrix_.num_rows() << " rows, " + << compact_matrix_.num_cols() << " columns, " + << compact_matrix_.num_entries() << " entries."; // TODO(user): Avoid doing the first phase checks when we know from the // incremental solve that the solution is already dual or primal feasible. @@ -589,10 +589,9 @@ void RevisedSimplex::UseSingletonColumnInInitialBasis(RowToColMapping* basis) { std::vector singleton_column; DenseRow cost_variation(num_cols_, 0.0); for (ColIndex col(0); col < num_cols_; ++col) { - if (matrix_with_slack_.column(col).num_entries() != 1) continue; + if (compact_matrix_.column(col).num_entries() != 1) continue; if (lower_bound_[col] == upper_bound_[col]) continue; - const Fractional slope = - matrix_with_slack_.column(col).GetFirstCoefficient(); + const Fractional slope = compact_matrix_.column(col).GetFirstCoefficient(); if (variable_values_.Get(col) == lower_bound_[col]) { cost_variation[col] = objective_[col] / std::abs(slope); } else { @@ -693,12 +692,6 @@ bool RevisedSimplex::InitializeMatrixAndTestIfUnchanged( // 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. - // - // TODO(user): use compact_matrix_ everywhere instead. - matrix_with_slack_.PopulateFromMatrix(lp.GetSparseMatrix()); return true; } @@ -716,8 +709,7 @@ bool RevisedSimplex::InitializeMatrixAndTestIfUnchanged( *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()); + // Initialize first_slack_. first_slack_col_ = lp.GetFirstSlackVariable(); // Initialize the new dimensions. @@ -727,7 +719,8 @@ bool RevisedSimplex::InitializeMatrixAndTestIfUnchanged( // 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. - compact_matrix_.PopulateFromMatrixView(matrix_with_slack_); + // TODO(user): This can be sped up by removing the MatrixView. + compact_matrix_.PopulateFromMatrixView(MatrixView(lp.GetSparseMatrix())); if (parameters_.use_transposed_matrix()) { transposed_matrix_.PopulateFromTranspose(compact_matrix_); } @@ -951,7 +944,7 @@ Status RevisedSimplex::CreateInitialBasis() { // the value of the boxed singleton column with a non-zero cost to the best // of their two bounds. for (ColIndex col(0); col < num_cols_; ++col) { - if (matrix_with_slack_.column(col).num_entries() != 1) continue; + if (compact_matrix_.column(col).num_entries() != 1) continue; const VariableStatus status = variables_info_.GetStatusRow()[col]; const Fractional objective = objective_[col]; if (objective > 0 && IsFinite(lower_bound_[col]) && @@ -991,7 +984,7 @@ Status RevisedSimplex::CreateInitialBasis() { return InitializeFirstBasis(basis); } if (parameters_.initial_basis() == GlopParameters::MAROS) { - InitialBasis initial_basis(matrix_with_slack_, objective_, lower_bound_, + InitialBasis initial_basis(compact_matrix_, objective_, lower_bound_, upper_bound_, variables_info_.GetTypeRow()); if (parameters_.use_dual_simplex()) { // This dual version only uses zero-cost columns to complete the @@ -1027,7 +1020,7 @@ Status RevisedSimplex::CreateInitialBasis() { // Then complete the basis with an advanced initial basis algorithm. VLOG(1) << "Trying to remove " << num_fixed_variables << " fixed variables from the initial basis."; - InitialBasis initial_basis(matrix_with_slack_, objective_, lower_bound_, + InitialBasis initial_basis(compact_matrix_, objective_, lower_bound_, upper_bound_, variables_info_.GetTypeRow()); if (parameters_.initial_basis() == GlopParameters::BIXBY) { @@ -1333,7 +1326,7 @@ void RevisedSimplex::SaveState() { RowIndex RevisedSimplex::ComputeNumberOfEmptyRows() { DenseBooleanColumn contains_data(num_rows_, false); for (ColIndex col(0); col < num_cols_; ++col) { - for (const SparseColumn::Entry e : matrix_with_slack_.column(col)) { + for (const SparseColumn::Entry e : compact_matrix_.column(col)) { contains_data[e.row()] = true; } } @@ -1350,7 +1343,7 @@ RowIndex RevisedSimplex::ComputeNumberOfEmptyRows() { ColIndex RevisedSimplex::ComputeNumberOfEmptyColumns() { ColIndex num_empty_cols(0); for (ColIndex col(0); col < num_cols_; ++col) { - if (matrix_with_slack_.column(col).IsEmpty()) { + if (compact_matrix_.column(col).IsEmpty()) { ++num_empty_cols; VLOG(1) << "Column " << col << " is empty."; } @@ -3094,10 +3087,10 @@ void RevisedSimplex::DisplayProblem() const { for (RowIndex row(0); row < num_rows_; ++row) { output = ""; for (ColIndex col(0); col < num_cols_; ++col) { - absl::StrAppend( - &output, StringifyMonomialWithFlags( - matrix_with_slack_.column(col).LookUpCoefficient(row), - variable_name_[col])); + absl::StrAppend(&output, + StringifyMonomialWithFlags( + compact_matrix_.column(col).LookUpCoefficient(row), + variable_name_[col])); } VLOG(3) << output << " = 0;"; } diff --git a/ortools/glop/revised_simplex.h b/ortools/glop/revised_simplex.h index 88a6eeec04..79a972454c 100644 --- a/ortools/glop/revised_simplex.h +++ b/ortools/glop/revised_simplex.h @@ -579,15 +579,7 @@ class RevisedSimplex { // it's as fast as std::unique_ptr as long as the size is properly reserved // beforehand. - // Temporary view of the matrix given to Solve(). Note that it is an error - // to access this view once Solve() is finished since there is no guarantee - // that the stored pointers are still valid. - // TODO(user): The matrix view is identical to the matrix of the linear - // program after pre-processing. Investigate if we could get rid of it and use - // compact_matrix_ in all places. - MatrixView matrix_with_slack_; - - // The compact version of matrix_with_slack_. + // Compact version of the matrix given to Solve(). CompactSparseMatrix compact_matrix_; // The transpose of compact_matrix_, it may be empty if it is not needed. diff --git a/ortools/lp_data/lp_utils.cc b/ortools/lp_data/lp_utils.cc index 4675879289..43e26c4c88 100644 --- a/ortools/lp_data/lp_utils.cc +++ b/ortools/lp_data/lp_utils.cc @@ -13,17 +13,28 @@ #include "ortools/lp_data/lp_utils.h" +#include "ortools/lp_data/sparse_column.h" + namespace operations_research { namespace glop { -Fractional SquaredNorm(const SparseColumn& v) { +template +Fractional SquaredNormTemplate(const SparseColumnLike& column) { Fractional sum(0.0); - for (const SparseColumn::Entry e : v) { + for (const SparseColumn::Entry e : column) { sum += Square(e.coefficient()); } return sum; } +Fractional SquaredNorm(const SparseColumn& v) { + return SquaredNormTemplate(v); +} + +Fractional SquaredNorm(const ColumnView& v) { + return SquaredNormTemplate(v); +} + Fractional PreciseSquaredNorm(const SparseColumn& v) { KahanSum sum; for (const SparseColumn::Entry e : v) { @@ -75,14 +86,23 @@ Fractional InfinityNorm(const DenseColumn& v) { return infinity_norm; } -Fractional InfinityNorm(const SparseColumn& v) { +template +Fractional InfinityNormTemplate(const SparseColumnLike& column) { Fractional infinity_norm = 0.0; - for (const SparseColumn::Entry e : v) { + for (const SparseColumn::Entry e : column) { infinity_norm = std::max(infinity_norm, fabs(e.coefficient())); } return infinity_norm; } +Fractional InfinityNorm(const SparseColumn& v) { + return InfinityNormTemplate(v); +} + +Fractional InfinityNorm(const ColumnView& v) { + return InfinityNormTemplate(v); +} + double Density(const DenseRow& row) { if (row.empty()) return 0.0; int sum = 0.0; @@ -110,7 +130,7 @@ void RemoveNearZeroEntries(Fractional threshold, DenseColumn* column) { } } -Fractional RestrictedInfinityNorm(const SparseColumn& column, +Fractional RestrictedInfinityNorm(const ColumnView& column, const DenseBooleanColumn& rows_to_consider, RowIndex* row_index) { Fractional infinity_norm = 0.0; @@ -123,7 +143,7 @@ Fractional RestrictedInfinityNorm(const SparseColumn& column, return infinity_norm; } -void SetSupportToFalse(const SparseColumn& column, DenseBooleanColumn* b) { +void SetSupportToFalse(const ColumnView& column, DenseBooleanColumn* b) { for (const SparseColumn::Entry e : column) { if (e.coefficient() != 0.0) { (*b)[e.row()] = false; @@ -131,7 +151,7 @@ void SetSupportToFalse(const SparseColumn& column, DenseBooleanColumn* b) { } } -bool IsDominated(const SparseColumn& column, const DenseColumn& radius) { +bool IsDominated(const ColumnView& column, const DenseColumn& radius) { for (const SparseColumn::Entry e : column) { DCHECK_GE(radius[e.row()], 0.0); if (fabs(e.coefficient()) > radius[e.row()]) return false; diff --git a/ortools/lp_data/lp_utils.h b/ortools/lp_data/lp_utils.h index 6ad45541a2..1534437252 100644 --- a/ortools/lp_data/lp_utils.h +++ b/ortools/lp_data/lp_utils.h @@ -142,6 +142,7 @@ Fractional PartialScalarProduct(const DenseRowOrColumn& u, // The precise version uses KahanSum and are about two times slower. Fractional SquaredNorm(const SparseColumn& v); Fractional SquaredNorm(const DenseColumn& column); +Fractional SquaredNorm(const ColumnView& v); Fractional PreciseSquaredNorm(const SparseColumn& v); Fractional PreciseSquaredNorm(const DenseColumn& column); Fractional PreciseSquaredNorm(const ScatteredColumn& v); @@ -149,6 +150,7 @@ Fractional PreciseSquaredNorm(const ScatteredColumn& v); // Returns the maximum of the |coefficients| of 'v'. Fractional InfinityNorm(const DenseColumn& v); Fractional InfinityNorm(const SparseColumn& v); +Fractional InfinityNorm(const ColumnView& v); // Returns the fraction of non-zero entries of the given row. // @@ -169,17 +171,17 @@ const DenseColumn& Transpose(const DenseRow& row); // Returns the maximum of the |coefficients| of the given column restricted // to the rows_to_consider. Also returns the first RowIndex 'row' that attains // this maximum. If the maximum is 0.0, then row_index is left untouched. -Fractional RestrictedInfinityNorm(const SparseColumn& column, +Fractional RestrictedInfinityNorm(const ColumnView& column, const DenseBooleanColumn& rows_to_consider, RowIndex* row_index); // Sets to false the entry b[row] if column[row] is non null. // Note that if 'b' was true only on the non-zero position of column, this can // be used as a fast way to clear 'b'. -void SetSupportToFalse(const SparseColumn& column, DenseBooleanColumn* b); +void SetSupportToFalse(const ColumnView& column, DenseBooleanColumn* b); // Returns true iff for all 'row' we have '|column[row]| <= radius[row]'. -bool IsDominated(const SparseColumn& column, const DenseColumn& radius); +bool IsDominated(const ColumnView& column, const DenseColumn& radius); // This cast based implementation should be safe, as long as DenseRow and // DenseColumn are implemented by the same underlying type. diff --git a/ortools/lp_data/sparse_column.h b/ortools/lp_data/sparse_column.h index e9bac5c9c7..902782dc43 100644 --- a/ortools/lp_data/sparse_column.h +++ b/ortools/lp_data/sparse_column.h @@ -97,6 +97,22 @@ class ColumnView { return Iterator(this->rows_, this->coefficients_, num_entries_); } + Fractional LookUpCoefficient(RowIndex index) const { + Fractional value(0.0); + for (const auto e : *this) { + if (e.row() == index) { + // Keep in mind the vector may contains several entries with the same + // index. In such a case the last one is returned. + // TODO(user): investigate whether an optimized version of + // LookUpCoefficient for "clean" columns yields speed-ups. + value = e.coefficient(); + } + } + return value; + } + + bool IsEmpty() const { return num_entries_ == EntryIndex(0); } + private: const EntryIndex num_entries_; const RowIndex* const rows_;