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

refactor code

Browse Source
This commit is contained in:
Laurent Perron
2017-03-28 16:07:49 +02:00
parent 6a7fef42ed
commit fe66d48cee
16 changed files with 395 additions and 281 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
src/glop/basis_representation.cc
View File

@@ -477,6 +477,12 @@ void BasisFactorization::RightSolveForProblemColumn(
lu_factorization_.RightSolveLForSparseColumn(matrix_.column(col), d,
non_zeros);
rank_one_factorization_.RightSolveWithNonZeros(d, non_zeros);
if (col >= right_pool_mapping_.size()) {
// This is needed because when we do an incremental solve with only new
// columns, we still reuse the current factorization without calling
// Refactorize() which would have resized this vector.
right_pool_mapping_.resize(col + 1, kInvalidCol);
}
if (non_zeros->empty()) {
right_pool_mapping_[col] = right_storage_.AddDenseColumn(*d);
} else {
@@ -546,7 +552,7 @@ Fractional BasisFactorization::ComputeInverseOneNorm() const {
Fractional column_norm = 0.0;
// Compute sum_i |inverse_ij|.
for (RowIndex row(0); row < num_rows; ++row) {
column_norm += fabs(right_hand_side[row]);
column_norm += std::abs(right_hand_side[row]);
}
// Compute max_j sum_i |inverse_ij|
norm = std::max(norm, column_norm);
@@ -566,7 +572,7 @@ Fractional BasisFactorization::ComputeInverseInfinityNorm() const {
RightSolve(&right_hand_side);
// Compute sum_j |inverse_ij|.
for (RowIndex row(0); row < num_rows; ++row) {
row_sum[row] += fabs(right_hand_side[row]);
row_sum[row] += std::abs(right_hand_side[row]);
}
}
// Compute max_i sum_j |inverse_ij|

2
src/glop/dual_edge_norms.cc
View File

@@ -61,7 +61,7 @@ void DualEdgeNorms::UpdateBeforeBasisPivot(
(sqrt(leaving_squared_norm) - sqrt(old_squared_norm)) /
sqrt(leaving_squared_norm);
stats_.edge_norms_accuracy.Add(estimated_edge_norms_accuracy);
if (fabs(estimated_edge_norms_accuracy) >
if (std::abs(estimated_edge_norms_accuracy) >
parameters_.recompute_edges_norm_threshold()) {
VLOG(1) << "Recomputing edge norms: " << sqrt(leaving_squared_norm)
<< " vs " << sqrt(old_squared_norm);

16
src/glop/entering_variable.cc
View File

@@ -201,7 +201,7 @@ Status EnteringVariable::DualChooseEnteringColumn(
*entering_col = kInvalidCol;
bound_flip_candidates->clear();
Fractional best_coeff = -1.0;
Fractional variation_magnitude = fabs(cost_variation);
Fractional variation_magnitude = std::abs(cost_variation);
equivalent_entering_choices_.clear();
while (!breakpoints.empty()) {
const ColWithRatio top = breakpoints.front();
@@ -327,7 +327,7 @@ Status EnteringVariable::DualPhaseIChooseEnteringColumn(
DCHECK_NE(variable_type[col], VariableType::FIXED_VARIABLE);
// Skip if the coeff is too small to be a numerically stable pivot.
if (fabs(update_coefficient[col]) < threshold) continue;
if (std::abs(update_coefficient[col]) < threshold) continue;
// We will add ratio * coeff to this column. cost_variation makes sure
// the leaving variable will be dual-feasible (its update coeff is
@@ -340,7 +340,7 @@ Status EnteringVariable::DualPhaseIChooseEnteringColumn(
// Only proceed if there is a transition, note that if reduced_costs[col]
// is close to zero, then the variable is supposed to be dual-feasible.
if (fabs(reduced_costs[col]) <= dual_feasibility_tolerance) {
if (std::abs(reduced_costs[col]) <= dual_feasibility_tolerance) {
// Continue if the variation goes in the dual-feasible direction.
if (coeff > 0 && !can_decrease.IsSet(col)) continue;
if (coeff < 0 && !can_increase.IsSet(col)) continue;
@@ -357,7 +357,7 @@ Status EnteringVariable::DualPhaseIChooseEnteringColumn(
// We are sure there is a transition, add it to the set of breakpoints.
breakpoints.push_back(
ColWithRatio(col, fabs(reduced_costs[col]), fabs(coeff)));
ColWithRatio(col, std::abs(reduced_costs[col]), std::abs(coeff)));
}
// Process the breakpoints in priority order.
@@ -372,7 +372,7 @@ Status EnteringVariable::DualPhaseIChooseEnteringColumn(
// numerically stable pivot.
*entering_col = kInvalidCol;
*step = -1.0;
Fractional improvement = fabs(cost_variation);
Fractional improvement = std::abs(cost_variation);
while (!breakpoints.empty()) {
const ColWithRatio top = breakpoints.front();
@@ -390,7 +390,7 @@ Status EnteringVariable::DualPhaseIChooseEnteringColumn(
// improvment, we also render it worse if we keep going in the same
// direction.
if (can_decrease.IsSet(top.col) && can_increase.IsSet(top.col) &&
fabs(reduced_costs[top.col]) > threshold) {
std::abs(reduced_costs[top.col]) > threshold) {
improvement -= top.coeff_magnitude;
}
@@ -444,7 +444,7 @@ void EnteringVariable::DantzigChooseEnteringColumn(ColIndex* entering_col) {
*entering_col = kInvalidCol;
for (const ColIndex col : reduced_costs_->GetDualInfeasiblePositions()) {
if (nested_pricing && !unused_columns_.IsSet(col)) continue;
const Fractional unormalized_price = fabs(reduced_costs[col]);
const Fractional unormalized_price = std::abs(reduced_costs[col]);
if (normalize) {
if (unormalized_price > best_price * matrix_column_norms[col]) {
best_price = unormalized_price / matrix_column_norms[col];
@@ -489,7 +489,7 @@ void EnteringVariable::NormalizedChooseEnteringColumn(ColIndex* entering_col) {
*entering_col = col;
}
} else {
const Fractional positive_reduced_cost = fabs(reduced_costs[col]);
const Fractional positive_reduced_cost = std::abs(reduced_costs[col]);
if (positive_reduced_cost >= best_price * weights[col]) {
if (positive_reduced_cost == best_price * weights[col]) {
equivalent_entering_choices_.push_back(col);

22
src/glop/initial_basis.cc
View File

@@ -51,7 +51,7 @@ void InitialBasis::CompleteBixbyBasis(ColIndex num_cols,
}
// This is 'v' in Bixby's paper.
DenseColumn scaled_diagonal_fabs(matrix_.num_rows(), kInfinity);
DenseColumn scaled_diagonal_abs(matrix_.num_rows(), kInfinity);
// Compute a list of candidate indices and sort them using the heuristic
// described in Bixby's paper.
@@ -77,7 +77,7 @@ void InitialBasis::CompleteBixbyBasis(ColIndex num_cols,
const Fractional kBixbyHighThreshold = 0.99;
if (candidate_coeff > kBixbyHighThreshold) {
enter_basis = true;
} else if (IsDominated(candidate_col, scaled_diagonal_fabs)) {
} else if (IsDominated(candidate_col, scaled_diagonal_abs)) {
candidate_coeff = RestrictedInfinityNorm(candidate_col, can_be_replaced,
&candidate_row);
if (candidate_coeff != 0.0) {
@@ -89,8 +89,8 @@ void InitialBasis::CompleteBixbyBasis(ColIndex num_cols,
can_be_replaced[candidate_row] = false;
SetSupportToFalse(candidate_col, &has_zero_coefficient);
const Fractional kBixbyLowThreshold = 0.01;
scaled_diagonal_fabs[candidate_row] =
kBixbyLowThreshold * fabs(candidate_coeff);
scaled_diagonal_abs[candidate_row] =
kBixbyLowThreshold * std::abs(candidate_coeff);
(*basis)[candidate_row] = candidate_col_index;
}
}
@@ -138,7 +138,7 @@ bool InitialBasis::CompleteTriangularBasis(ColIndex num_cols,
max_scaled_abs_cost_ = 0.0;
for (ColIndex col(0); col < num_cols; ++col) {
max_scaled_abs_cost_ =
std::max(max_scaled_abs_cost_, fabs(objective_[col]));
std::max(max_scaled_abs_cost_, std::abs(objective_[col]));
if (residual_pattern.ColDegree(col) == 1) {
residual_singleton_column.push_back(col);
}
@@ -172,7 +172,7 @@ bool InitialBasis::CompleteTriangularBasis(ColIndex num_cols,
Fractional coeff = 0.0;
Fractional max_magnitude = 0.0;
for (const SparseColumn::Entry e : matrix_.column(candidate)) {
max_magnitude = std::max(max_magnitude, fabs(e.coefficient()));
max_magnitude = std::max(max_magnitude, std::abs(e.coefficient()));
if (can_be_replaced[e.row()]) {
row = e.row();
coeff = e.coefficient();
@@ -180,11 +180,11 @@ bool InitialBasis::CompleteTriangularBasis(ColIndex num_cols,
}
}
const Fractional kStabilityThreshold = 0.01;
if (fabs(coeff) < kStabilityThreshold * max_magnitude) continue;
if (std::abs(coeff) < kStabilityThreshold * max_magnitude) continue;
DCHECK_NE(kInvalidRow, row);
partial_diagonal_product *= coeff;
if (fabs(partial_diagonal_product) < kMinimumProductMagnitude) {
if (std::abs(partial_diagonal_product) < kMinimumProductMagnitude) {
LOG(INFO) << "Numerical difficulties detected. The product of the "
<< "diagonal coefficients is currently equal to "
<< partial_diagonal_product;
@@ -204,7 +204,7 @@ bool InitialBasis::CompleteTriangularBasis(ColIndex num_cols,
}
}
return fabs(partial_diagonal_product) >= kMinimumProductMagnitude;
return std::abs(partial_diagonal_product) >= kMinimumProductMagnitude;
}
void InitialBasis::ComputeCandidates(ColIndex num_cols,
@@ -216,7 +216,7 @@ void InitialBasis::ComputeCandidates(ColIndex num_cols,
matrix_.column(col).num_entries() > 0) {
candidates->push_back(col);
max_scaled_abs_cost_ =
std::max(max_scaled_abs_cost_, fabs(objective_[col]));
std::max(max_scaled_abs_cost_, std::abs(objective_[col]));
}
}
const Fractional kBixbyWeight = 1000.0;
@@ -258,7 +258,7 @@ Fractional InitialBasis::GetColumnPenalty(ColIndex col) const {
if (type == VariableType::UPPER_AND_LOWER_BOUNDED) {
penalty = lower_bound_[col] - upper_bound_[col];
}
return penalty + fabs(objective_[col]) / max_scaled_abs_cost_;
return penalty + std::abs(objective_[col]) / max_scaled_abs_cost_;
}
bool InitialBasis::BixbyColumnComparator::operator()(ColIndex col_a,

6
src/glop/lu_factorization.cc
View File

@@ -559,7 +559,7 @@ Fractional LuFactorization::ComputeInverseOneNorm() const {
Fractional column_norm = 0.0;
// Compute sum_i |basis_matrix_ij|.
for (RowIndex row(0); row < num_rows; ++row) {
column_norm += fabs(right_hand_side[row]);
column_norm += std::abs(right_hand_side[row]);
}
// Compute max_j sum_i |basis_matrix_ij|
norm = std::max(norm, column_norm);
@@ -579,7 +579,7 @@ Fractional LuFactorization::ComputeInverseInfinityNorm() const {
RightSolve(&right_hand_side);
// Compute sum_j |basis_matrix_ij|.
for (RowIndex row(0); row < num_rows; ++row) {
row_sum[row] += fabs(right_hand_side[row]);
row_sum[row] += std::abs(right_hand_side[row]);
}
}
// Compute max_i sum_j |basis_matrix_ij|
@@ -646,7 +646,7 @@ bool LuFactorization::CheckFactorization(const MatrixView& matrix,
for (ColIndex col(0); col < should_be_zero.num_cols(); ++col) {
for (const SparseColumn::Entry e : should_be_zero.column(col)) {
const Fractional magnitude = fabs(e.coefficient());
const Fractional magnitude = std::abs(e.coefficient());
if (magnitude > tolerance) {
VLOG(2) << magnitude << " != 0, at column " << col;
return false;

8
src/glop/markowitz.cc
View File

@@ -76,7 +76,7 @@ Status Markowitz::ComputeRowAndColumnPermutation(const MatrixView& basis_matrix,
// but its magnitude can be really close to zero. In both cases, we
// report the singularity of the matrix.
if (pivot_row == kInvalidRow || pivot_col == kInvalidCol ||
fabs(pivot_coefficient) <= singularity_threshold) {
std::abs(pivot_coefficient) <= singularity_threshold) {
RETURN_AND_LOG_ERROR(Status::ERROR_LU,
StringPrintf("The matrix is singular! pivot = %E",
pivot_coefficient));
@@ -427,7 +427,7 @@ int64 Markowitz::FindPivot(const RowPermutation& row_perm,
Fractional max_magnitude = 0.0;
for (const SparseColumn::Entry e : column) {
max_magnitude = std::max(max_magnitude, fabs(e.coefficient()));
max_magnitude = std::max(max_magnitude, std::abs(e.coefficient()));
}
if (max_magnitude == 0.0) {
// All symbolic non-zero entries have been cancelled!
@@ -438,7 +438,7 @@ int64 Markowitz::FindPivot(const RowPermutation& row_perm,
const Fractional skip_threshold = threshold * max_magnitude;
for (const SparseColumn::Entry e : column) {
const Fractional magnitude = fabs(e.coefficient());
const Fractional magnitude = std::abs(e.coefficient());
if (magnitude < skip_threshold) continue;
const int row_degree = residual_matrix_non_zero_.RowDegree(e.row());
@@ -446,7 +446,7 @@ int64 Markowitz::FindPivot(const RowPermutation& row_perm,
DCHECK_NE(markowitz_number, 0);
if (markowitz_number < min_markowitz_number ||
((markowitz_number == min_markowitz_number) &&
magnitude > fabs(*pivot_coefficient))) {
magnitude > std::abs(*pivot_coefficient))) {
min_markowitz_number = markowitz_number;
*pivot_col = col;
*pivot_row = e.row();

55
src/glop/preprocessor.cc
View File

@@ -306,12 +306,12 @@ VariableStatus ComputeVariableStatus(Fractional value, Fractional lower_bound,
// Returns the input with the smallest magnitude or zero if both are infinite.
Fractional MinInMagnitudeOrZeroIfInfinite(Fractional a, Fractional b) {
const Fractional value = fabs(a) < fabs(b) ? a : b;
const Fractional value = std::abs(a) < std::abs(b) ? a : b;
return IsFinite(value) ? value : 0.0;
}
Fractional MagnitudeOrZeroIfInfinite(Fractional value) {
return IsFinite(value) ? fabs(value) : 0.0;
return IsFinite(value) ? std::abs(value) : 0.0;
}
// Returns the maximum magnitude of the finite variable bounds of the given
@@ -598,7 +598,7 @@ bool ProportionalColumnPreprocessor::Run(LinearProgram* lp,
int num_merged = 0;
for (++i; i < sorted_columns.size(); ++i) {
if (sorted_columns[i].representative != target_representative) break;
if (fabs(sorted_columns[i].scaled_cost - target_scaled_cost) >=
if (std::abs(sorted_columns[i].scaled_cost - target_scaled_cost) >=
parameters_.preprocessor_zero_tolerance()) {
break;
}
@@ -710,7 +710,7 @@ void ProportionalColumnPreprocessor::RecoverSolution(
const Fractional bound_factor =
column_factors_[col] / column_factors_[representative];
const Fractional scaled_distance =
distance_to_bound[representative] / fabs(bound_factor);
distance_to_bound[representative] / std::abs(bound_factor);
const Fractional width = upper_bounds_[col] - lower_bounds_[col];
const bool to_upper_bound =
(bound_factor > 0.0) == is_distance_to_upper_bound[representative];
@@ -720,7 +720,7 @@ void ProportionalColumnPreprocessor::RecoverSolution(
solution->variable_statuses[col] =
ComputeVariableStatus(solution->primal_values[col],
lower_bounds_[col], upper_bounds_[col]);
distance_to_bound[representative] -= width * fabs(bound_factor);
distance_to_bound[representative] -= width * std::abs(bound_factor);
} else {
solution->primal_values[col] =
to_upper_bound ? upper_bounds_[col] - scaled_distance
@@ -890,7 +890,8 @@ bool ProportionalRowPreprocessor::Run(LinearProgram* lp,
// proportionality factor.
RowIndex new_representative = lower_source;
RowIndex other = upper_source;
if (fabs(row_factors_[new_representative]) < fabs(row_factors_[other])) {
if (std::abs(row_factors_[new_representative]) <
std::abs(row_factors_[other])) {
std::swap(new_representative, other);
}
@@ -1533,7 +1534,7 @@ bool DoubletonFreeColumnPreprocessor::Run(LinearProgram* lp,
// We prefer deleting the row with the larger coefficient magnitude because
// we will divide by this magnitude. TODO(user): Impact?
if (fabs(r.coeff[DELETED]) < fabs(r.coeff[MODIFIED])) {
if (std::abs(r.coeff[DELETED]) < std::abs(r.coeff[MODIFIED])) {
std::swap(r.coeff[DELETED], r.coeff[MODIFIED]);
std::swap(r.row[DELETED], r.row[MODIFIED]);
}
@@ -1577,9 +1578,9 @@ bool DoubletonFreeColumnPreprocessor::Run(LinearProgram* lp,
// the numerical error in the formula above, we have a really low
// objective instead. The logic is the same as in
// AddMultipleToSparseVectorAndIgnoreCommonIndex().
if (fabs(new_objective) > std::numeric_limits<Fractional>::epsilon() *
2.0 *
fabs(lp->objective_coefficients()[col])) {
if (std::abs(new_objective) >
std::numeric_limits<Fractional>::epsilon() * 2.0 *
std::abs(lp->objective_coefficients()[col])) {
lp->SetObjectiveCoefficient(col, new_objective);
} else {
lp->SetObjectiveCoefficient(col, 0.0);
@@ -1791,10 +1792,10 @@ bool UnconstrainedVariablePreprocessor::Run(LinearProgram* lp,
}
}
// Change the rhs to reflect the fixed variables. Note that is is important to
// do that after all the calls to RemoveZeroCostUnconstrainedVariable()
// because RemoveZeroCostUnconstrainedVariable() needs to store the rhs before
// this modification!
// Change the rhs to reflect the fixed variables. Note that is important to do
// that after all the calls to RemoveZeroCostUnconstrainedVariable() because
// RemoveZeroCostUnconstrainedVariable() needs to store the rhs before this
// modification!
const ColIndex end = column_deletion_helper_.GetMarkedColumns().size();
for (ColIndex col(0); col < end; ++col) {
if (column_deletion_helper_.IsColumnMarked(col)) {
@@ -1855,7 +1856,7 @@ void UnconstrainedVariablePreprocessor::RecoverSolution(
// unbounded direction.
if (activity * activity_sign_correction_[row] < 0.0) {
const Fractional bound = activity / e.coefficient();
if (fabs(bound) > fabs(primal_value_shift)) {
if (std::abs(bound) > std::abs(primal_value_shift)) {
primal_value_shift = bound;
row_at_bound = row;
}
@@ -2266,7 +2267,7 @@ void SingletonPreprocessor::DeleteSingletonColumnInEquality(
// tolerances in a few preprocessors. Like an empty column with a cost of
// 1e-17 and unbounded towards infinity is currently implying that the
// problem is unbounded. This will need fixing.
if (fabs(new_cost) < parameters_.preprocessor_zero_tolerance()) {
if (std::abs(new_cost) < parameters_.preprocessor_zero_tolerance()) {
new_cost = 0.0;
}
lp->SetObjectiveCoefficient(col, new_cost);
@@ -2613,7 +2614,7 @@ bool RemoveNearZeroEntriesPreprocessor::Run(LinearProgram* lp,
// and column bounds, and "propagate" the bounds as much as possible so we
// can use this better estimate here and remove more near-zero entries.
const Fractional max_magnitude =
std::max(fabs(lower_bound), fabs(upper_bound));
std::max(std::abs(lower_bound), std::abs(upper_bound));
if (max_magnitude == kInfinity || max_magnitude == 0) continue;
const Fractional threshold = allowed_impact / max_magnitude;
lp->GetMutableSparseColumn(col)->RemoveNearZeroEntriesWithWeights(
@@ -2621,7 +2622,7 @@ bool RemoveNearZeroEntriesPreprocessor::Run(LinearProgram* lp,
if (lp->objective_coefficients()[col] != 0.0 &&
num_non_zero_objective_coefficients *
fabs(lp->objective_coefficients()[col]) <
std::abs(lp->objective_coefficients()[col]) <
threshold) {
lp->SetObjectiveCoefficient(col, 0.0);
++num_zeroed_objective_coefficients;
@@ -2756,6 +2757,7 @@ bool DoubletonEqualityRowPreprocessor::Run(LinearProgram* lp,
continue;
}
// Look at the bounds of both variables and exit early if we can delegate
// to another pre-processor; otherwise adjust the bounds of the remaining
// variable as necessary.
@@ -2777,6 +2779,8 @@ bool DoubletonEqualityRowPreprocessor::Run(LinearProgram* lp,
status_ = ProblemStatus::ABNORMAL;
break;
}
Fractional lb = r.lb[MODIFIED];
Fractional ub = r.ub[MODIFIED];
Fractional carried_over_lb =
@@ -2934,6 +2938,17 @@ void DoubletonEqualityRowPreprocessor::RecoverSolution(
}
}
void DoubletonEqualityRowPreprocessor::
SwapDeletedAndModifiedVariableRestoreInfo(RestoreInfo* r) {
using std::swap;
swap(r->col[DELETED], r->col[MODIFIED]);
swap(r->coeff[DELETED], r->coeff[MODIFIED]);
swap(r->lb[DELETED], r->lb[MODIFIED]);
swap(r->ub[DELETED], r->ub[MODIFIED]);
swap(r->column[DELETED], r->column[MODIFIED]);
swap(r->objective_coefficient[DELETED], r->objective_coefficient[MODIFIED]);
}
// --------------------------------------------------------
// DualizerPreprocessor
// --------------------------------------------------------
@@ -2949,7 +2964,7 @@ bool DualizerPreprocessor::Run(LinearProgram* lp, TimeLimit* time_limit) {
primal_num_rows_ = lp->num_constraints();
primal_is_maximization_problem_ = lp->IsMaximizationProblem();
// If we need to decide wether or not to take the dual, we only take it when
// If we need to decide whether or not to take the dual, we only take it when
// the matrix has more rows than columns. The number of rows of a linear
// program gives the size of the square matrices we need to invert and the
// order of iterations of the simplex method. So solving a program with less
@@ -3219,7 +3234,7 @@ bool ShiftVariableBoundsPreprocessor::Run(LinearProgram* lp,
if (0.0 < variable_initial_lbs_[col] || 0.0 > variable_initial_ubs_[col]) {
Fractional offset = MinInMagnitudeOrZeroIfInfinite(
variable_initial_lbs_[col], variable_initial_ubs_[col]);
if (in_mip_context_ && lp->is_variable_integer()[col]) {
if (in_mip_context_ && lp->IsVariableInteger(col)) {
// In the integer case, we truncate the number because if for instance
// the lower bound is a positive integer + epsilon, we only want to
// shift by the integer and leave the lower bound at epsilon.

5
src/glop/preprocessor.h
View File

@@ -394,7 +394,6 @@ class SingletonPreprocessor : public Preprocessor {
~SingletonPreprocessor() final {}
bool Run(LinearProgram* linear_program, TimeLimit* time_limit) final;
void RecoverSolution(ProblemSolution* solution) const final;
void UseInMipContext() final { LOG(FATAL) << "Not implemented."; }
private:
// Returns the MatrixEntry of the given singleton row or column, taking into
@@ -550,7 +549,6 @@ class ImpliedFreePreprocessor : public Preprocessor {
~ImpliedFreePreprocessor() final {}
bool Run(LinearProgram* linear_program, TimeLimit* time_limit) final;
void RecoverSolution(ProblemSolution* solution) const final;
void UseInMipContext() final { LOG(FATAL) << "Not implemented."; }
private:
// This preprocessor adds fixed offsets to some variables. We remember those
@@ -753,7 +751,6 @@ class DoubletonEqualityRowPreprocessor : public Preprocessor {
~DoubletonEqualityRowPreprocessor() final {}
bool Run(LinearProgram* linear_program, TimeLimit* time_limit) final;
void RecoverSolution(ProblemSolution* solution) const final;
void UseInMipContext() final { LOG(FATAL) << "Not implemented."; }
private:
enum ColChoice {
@@ -799,6 +796,8 @@ class DoubletonEqualityRowPreprocessor : public Preprocessor {
ColChoiceAndStatus bound_backtracking_at_lower_bound;
ColChoiceAndStatus bound_backtracking_at_upper_bound;
};
void SwapDeletedAndModifiedVariableRestoreInfo(RestoreInfo* r);
std::vector<RestoreInfo> restore_stack_;
DISALLOW_COPY_AND_ASSIGN(DoubletonEqualityRowPreprocessor);

6
src/glop/primal_edge_norms.cc
View File

@@ -80,7 +80,7 @@ void PrimalEdgeNorms::TestEnteringEdgeNormPrecision(
const Fractional estimated_edges_norm_accuracy =
(precise_norm - sqrt(old_squared_norm)) / precise_norm;
stats_.edges_norm_accuracy.Add(estimated_edges_norm_accuracy);
if (fabs(estimated_edges_norm_accuracy) >
if (std::abs(estimated_edges_norm_accuracy) >
parameters_.recompute_edges_norm_threshold()) {
VLOG(1) << "Recomputing edge norms: " << sqrt(precise_squared_norm)
<< " vs " << sqrt(old_squared_norm);
@@ -285,12 +285,12 @@ void PrimalEdgeNorms::UpdateDevexWeights(
// Compared to steepest edge update, the DEVEX weight uses the largest of the
// norms of two vectors to approximate the norm of the sum.
const Fractional entering_norm = sqrt(PreciseSquaredNorm(direction));
const Fractional pivot_magnitude = fabs(direction[leaving_row]);
const Fractional pivot_magnitude = std::abs(direction[leaving_row]);
const Fractional leaving_norm =
std::max(1.0, entering_norm / pivot_magnitude);
for (const ColIndex col : update_row.GetNonZeroPositions()) {
const Fractional coeff = update_row.GetCoefficient(col);
const Fractional update_vector_norm = fabs(coeff) * leaving_norm;
const Fractional update_vector_norm = std::abs(coeff) * leaving_norm;
devex_weights_[col] = std::max(devex_weights_[col], update_vector_norm);
}
devex_weights_[leaving_col] = leaving_norm;

14
src/glop/reduced_costs.cc
View File

@@ -95,9 +95,9 @@ bool ReducedCosts::TestEnteringReducedCostPrecision(
const Fractional estimated_reduced_costs_accuracy =
old_reduced_cost - precise_reduced_cost;
const Fractional scale =
(fabs(precise_reduced_cost) <= 1.0) ? 1.0 : precise_reduced_cost;
(std::abs(precise_reduced_cost) <= 1.0) ? 1.0 : precise_reduced_cost;
stats_.reduced_costs_accuracy.Add(estimated_reduced_costs_accuracy / scale);
if (fabs(estimated_reduced_costs_accuracy) / scale >
if (std::abs(estimated_reduced_costs_accuracy) / scale >
parameters_.recompute_reduced_costs_threshold()) {
VLOG(1) << "Recomputing reduced costs: " << precise_reduced_cost << " vs "
<< old_reduced_cost;
@@ -127,7 +127,7 @@ Fractional ReducedCosts::ComputeMaximumDualResidual() const {
const Fractional residual =
objective_[basic_col] + objective_perturbation_[basic_col] -
matrix_.ColumnScalarProduct(basic_col, dual_values);
dual_residual_error = std::max(dual_residual_error, fabs(residual));
dual_residual_error = std::max(dual_residual_error, std::abs(residual));
}
return dual_residual_error;
}
@@ -144,7 +144,7 @@ Fractional ReducedCosts::ComputeMaximumDualInfeasibility() const {
if ((can_increase.IsSet(col) && rc < 0.0) ||
(can_decrease.IsSet(col) && rc > 0.0)) {
maximum_dual_infeasibility =
std::max(maximum_dual_infeasibility, fabs(rc));
std::max(maximum_dual_infeasibility, std::abs(rc));
}
}
return maximum_dual_infeasibility;
@@ -161,7 +161,7 @@ Fractional ReducedCosts::ComputeSumOfDualInfeasibilities() const {
const Fractional rc = reduced_costs_[col];
if ((can_increase.IsSet(col) && rc < 0.0) ||
(can_decrease.IsSet(col) && rc > 0.0)) {
dual_infeasibility_sum += fabs(fabs(rc));
dual_infeasibility_sum += std::abs(std::abs(rc));
}
}
return dual_infeasibility_sum;
@@ -372,7 +372,7 @@ void ReducedCosts::ComputeReducedCosts() {
// We also compute the dual residual error y.B - c_B.
if (is_basic.IsSet(col)) {
dual_residual_error =
std::max(dual_residual_error, fabs(reduced_costs_[col]));
std::max(dual_residual_error, std::abs(reduced_costs_[col]));
}
}
} else {
@@ -392,7 +392,7 @@ void ReducedCosts::ComputeReducedCosts() {
if (is_basic.IsSet(col)) {
thread_local_dual_residual_error[omp_get_thread_num()] =
std::max(thread_local_dual_residual_error[omp_get_thread_num()],
fabs(reduced_costs_[col]));
std::abs(reduced_costs_[col]));
}
}
// end of omp parallel for

368
src/glop/revised_simplex.cc
View File

@@ -26,6 +26,7 @@
#include "base/integral_types.h"
#include "base/logging.h"
#include "base/stringprintf.h"
#include "base/join.h"
#include "glop/initial_basis.h"
#include "glop/parameters.pb.h"
#include "lp_data/lp_data.h"
@@ -113,7 +114,7 @@ RevisedSimplex::RevisedSimplex()
test_lu_(),
feasibility_phase_(true),
random_("This is a deterministic seed.") {
PropagateParameters();
SetParameters(parameters_);
}
void RevisedSimplex::ClearStateForNextSolve() {
@@ -141,7 +142,6 @@ Status RevisedSimplex::Solve(const LinearProgram& lp, TimeLimit* time_limit) {
// Initialization. Note That Initialize() must be called first since it
// analyzes the current solver state.
const double start_time = time_limit->GetElapsedTime();
random_.Reset(parameters_.random_seed());
RETURN_IF_ERROR(Initialize(lp));
dual_infeasibility_improvement_direction_.clear();
update_row_.Invalidate();
@@ -155,7 +155,7 @@ Status RevisedSimplex::Solve(const LinearProgram& lp, TimeLimit* time_limit) {
optimization_time_ = 0.0;
total_time_ = 0.0;
if (FLAGS_v > 0) {
if (VLOG_IS_ON(1)) {
ComputeNumberOfEmptyRows();
ComputeNumberOfEmptyColumns();
DisplayBasicVariableStatistics();
@@ -174,6 +174,8 @@ Status RevisedSimplex::Solve(const LinearProgram& lp, TimeLimit* time_limit) {
current_objective_ = objective_;
// TODO(user): Avoid doing the first phase checks when we know from the
// incremental solve that the solution is already dual or primal feasible.
VLOG(1) << "------ First phase: feasibility.";
entering_variable_.SetPricingRule(parameters_.feasibility_rule());
if (use_dual) {
@@ -421,6 +423,11 @@ const DenseColumn& RevisedSimplex::GetDualRay() const {
return solution_dual_ray_;
}
const DenseRow& RevisedSimplex::GetDualRayRowCombination() const {
DCHECK_EQ(problem_status_, ProblemStatus::DUAL_UNBOUNDED);
return solution_dual_ray_row_combination_;
}
ColIndex RevisedSimplex::GetBasis(RowIndex row) const { return basis_[row]; }
const BasisFactorization& RevisedSimplex::GetBasisFactorization() const {
@@ -491,14 +498,14 @@ VariableStatus RevisedSimplex::ComputeDefaultVariableStatus(
// Returns the bound with the lowest magnitude. Note that it must be finite
// because the VariableStatus::FREE case was tested earlier.
DCHECK(IsFinite(lower_bound_[col]) || IsFinite(upper_bound_[col]));
return fabs(lower_bound_[col]) <= fabs(upper_bound_[col])
return std::abs(lower_bound_[col]) <= std::abs(upper_bound_[col])
? VariableStatus::AT_LOWER_BOUND
: VariableStatus::AT_UPPER_BOUND;
}
void RevisedSimplex::SetNonBasicVariableStatusAndDeriveValue(
ColIndex col, VariableStatus status) {
variables_info_.Update(col, status);
variables_info_.UpdateToNonBasicStatus(col, status);
variable_values_.SetNonBasicVariableValueFromStatus(col);
}
@@ -594,9 +601,9 @@ void RevisedSimplex::UseSingletonColumnInInitialBasis(RowToColMapping* basis) {
const Fractional slope =
matrix_with_slack_.column(col).GetFirstCoefficient();
if (variable_values_.Get(col) == lower_bound_[col]) {
cost_variation[col] = objective_[col] / fabs(slope);
cost_variation[col] = objective_[col] / std::abs(slope);
} else {
cost_variation[col] = -objective_[col] / fabs(slope);
cost_variation[col] = -objective_[col] / std::abs(slope);
}
singleton_column.push_back(col);
}
@@ -882,40 +889,31 @@ void RevisedSimplex::InitializeVariableStatusesForWarmStart(
status = state.statuses[col - num_new_cols];
}
// Remove incompatibilities between the warm status and the variable bounds.
// We use the default status as an indication of the bounds type.
if ((status == VariableStatus::FREE &&
default_status != VariableStatus::FREE) ||
(status == VariableStatus::FIXED_VALUE &&
default_status != VariableStatus::FIXED_VALUE) ||
(status == VariableStatus::AT_LOWER_BOUND &&
lower_bound_[col] == -kInfinity) ||
(status == VariableStatus::AT_UPPER_BOUND &&
upper_bound_[col] == kInfinity)) {
status = default_status;
}
if (default_status == VariableStatus::FIXED_VALUE &&
(status == VariableStatus::AT_UPPER_BOUND ||
status == VariableStatus::AT_LOWER_BOUND)) {
status = VariableStatus::FIXED_VALUE;
}
// Do not allow more than num_rows_ VariableStatus::BASIC variables.
if (status == VariableStatus::BASIC) {
// Do not allow more than num_rows_ VariableStatus::BASIC variables.
if (num_basic_variables == num_rows_) {
VLOG(1) << "Too many basic variables in the warm-start basis."
<< "Only keeping the first ones as VariableStatus::BASIC.";
status = default_status;
SetNonBasicVariableStatusAndDeriveValue(col, default_status);
} else {
++num_basic_variables;
variables_info_.UpdateToBasicStatus(col);
}
}
// Set the status.
if (status != VariableStatus::BASIC) {
SetNonBasicVariableStatusAndDeriveValue(col, status);
} else {
variables_info_.Update(col, VariableStatus::BASIC);
// Remove incompatibilities between the warm status and the variable
// bounds. We use the default status as an indication of the bounds
// type.
if ((status != default_status) &&
((default_status == VariableStatus::FIXED_VALUE) ||
(status == VariableStatus::FREE) ||
(status == VariableStatus::FIXED_VALUE) ||
(status == VariableStatus::AT_LOWER_BOUND &&
lower_bound_[col] == -kInfinity) ||
(status == VariableStatus::AT_UPPER_BOUND &&
upper_bound_[col] == kInfinity))) {
status = default_status;
}
SetNonBasicVariableStatusAndDeriveValue(col, status);
}
}
}
@@ -1083,7 +1081,7 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) {
// Computes the variable name as soon as possible for logging.
// TODO(user): do we really need to store them? we could just compute them
// on the fly since we do not need the speed.
if (FLAGS_v > 0) {
if (VLOG_IS_ON(1)) {
SetVariableNames();
}
@@ -1329,7 +1327,7 @@ void RevisedSimplex::ComputeDirection(ColIndex col) {
direction_infinity_norm_ = 0.0;
for (const RowIndex row : direction_non_zero_) {
direction_infinity_norm_ =
std::max(direction_infinity_norm_, fabs(direction_[row]));
std::max(direction_infinity_norm_, std::abs(direction_[row]));
}
IF_STATS_ENABLED(ratio_test_stats_.direction_density.Add(
num_rows_ == 0 ? 0.0
@@ -1352,7 +1350,7 @@ void RevisedSimplex::SkipVariableForRatioTest(RowIndex row) {
// during the ratio test. The direction vector will be restored later from
// the information in direction_ignored_position_.
IF_STATS_ENABLED(
ratio_test_stats_.abs_skipped_pivot.Add(fabs(direction_[row])));
ratio_test_stats_.abs_skipped_pivot.Add(std::abs(direction_[row])));
VLOG(1) << "Skipping leaving variable with coefficient " << direction_[row];
direction_ignored_position_.SetCoefficient(row, direction_[row]);
direction_[row] = 0.0;
@@ -1397,7 +1395,7 @@ Fractional RevisedSimplex::ComputeHarrisRatioAndLeavingCandidates(
leaving_candidates->Clear();
const Fractional threshold = parameters_.ratio_test_zero_threshold();
for (const RowIndex row : direction_non_zero_) {
const Fractional magnitude = fabs(direction_[row]);
const Fractional magnitude = std::abs(direction_[row]);
if (magnitude < threshold) continue;
Fractional ratio = GetRatio<is_entering_reduced_cost_positive>(row);
// TODO(user): The perturbation is currently disabled, so no need to test
@@ -1405,7 +1403,7 @@ Fractional RevisedSimplex::ComputeHarrisRatioAndLeavingCandidates(
if (false && ratio < 0.0) {
// If the variable is already pass its bound, we use the perturbed version
// of the bound (if bound_perturbation_[basis_[row]] is not zero).
ratio += fabs(bound_perturbation_[basis_[row]] / direction_[row]);
ratio += std::abs(bound_perturbation_[basis_[row]] / direction_[row]);
}
if (ratio <= harris_ratio) {
leaving_candidates->SetCoefficient(row, ratio);
@@ -1510,7 +1508,7 @@ Status RevisedSimplex::ChooseLeavingVariableRow(
// If the magnitudes are the same, we choose the leaving variable with
// what is probably the more stable ratio, see
// IsRatioMoreOrEquallyStable().
const Fractional candidate_magnitude = fabs(direction_[row]);
const Fractional candidate_magnitude = std::abs(direction_[row]);
if (candidate_magnitude < pivot_magnitude) continue;
if (candidate_magnitude == pivot_magnitude) {
if (!IsRatioMoreOrEquallyStable(ratio, current_ratio)) continue;
@@ -1590,7 +1588,7 @@ Status RevisedSimplex::ChooseLeavingVariableRow(
// Note(user): This reduces quite a bit the number of iterations.
// What is its impact? Is it dangerous?
if (fabs(direction_[*leaving_row]) <
if (std::abs(direction_[*leaving_row]) <
parameters_.minimum_acceptable_pivot()) {
SkipVariableForRatioTest(*leaving_row);
continue;
@@ -1633,7 +1631,7 @@ Status RevisedSimplex::ChooseLeavingVariableRow(
equivalent_leaving_choices_.size());
}
if (*leaving_row != kInvalidRow) {
ratio_test_stats_.abs_used_pivot.Add(fabs(direction_[*leaving_row]));
ratio_test_stats_.abs_used_pivot.Add(std::abs(direction_[*leaving_row]));
}
});
return Status::OK;
@@ -1725,7 +1723,7 @@ void RevisedSimplex::PrimalPhaseIChooseLeavingVariableRow(
for (RowIndex row : direction_non_zero_) {
const Fractional direction =
reduced_cost > 0.0 ? direction_[row] : -direction_[row];
const Fractional magnitude = fabs(direction);
const Fractional magnitude = std::abs(direction);
if (magnitude < tolerance) continue;
// Computes by how much we can add 'direction' to the basic variable value
@@ -1766,7 +1764,7 @@ void RevisedSimplex::PrimalPhaseIChooseLeavingVariableRow(
// Select the last breakpoint that still improves the infeasibility and has
// the largest coefficient magnitude.
Fractional improvement = fabs(reduced_cost);
Fractional improvement = std::abs(reduced_cost);
Fractional best_magnitude = 0.0;
*leaving_row = kInvalidRow;
while (!breakpoints.empty()) {
@@ -2635,15 +2633,22 @@ Status RevisedSimplex::DualMinimize(TimeLimit* time_limit) {
} else {
problem_status_ = ProblemStatus::DUAL_UNBOUNDED;
solution_dual_ray_ = update_row_.GetUnitRowLeftInverse().dense_column;
update_row_.RecomputeFullUpdateRow(leaving_row);
solution_dual_ray_row_combination_.assign(num_cols_, 0.0);
for (const ColIndex col : update_row_.GetNonZeroPositions()) {
solution_dual_ray_row_combination_[col] =
update_row_.GetCoefficient(col);
}
if (cost_variation < 0) {
ChangeSign(&solution_dual_ray_);
ChangeSign(&solution_dual_ray_row_combination_);
}
}
return Status::OK;
}
// If the coefficient is too small, we recompute the reduced costs.
if (fabs(entering_coeff) < parameters_.dual_small_pivot_threshold() &&
if (std::abs(entering_coeff) < parameters_.dual_small_pivot_threshold() &&
!reduced_costs_.AreReducedCostsPrecise()) {
VLOG(1) << "Trying not to pivot by " << entering_coeff;
refactorize = true;
@@ -2655,7 +2660,7 @@ Status RevisedSimplex::DualMinimize(TimeLimit* time_limit) {
// direction_[leaving_row] is actually 0 which causes a floating
// point exception below.
ComputeDirection(entering_col);
if (fabs(direction_[leaving_row]) <
if (std::abs(direction_[leaving_row]) <
parameters_.minimum_acceptable_pivot()) {
VLOG(1) << "Do not pivot by " << entering_coeff
<< " because the direction is " << direction_[leaving_row];
@@ -2719,7 +2724,8 @@ Status RevisedSimplex::DualMinimize(TimeLimit* time_limit) {
// This is slow, but otherwise we have a really bad precision on the
// variable values ...
if (fabs(primal_step) * parameters_.primal_feasibility_tolerance() > 1.0) {
if (std::abs(primal_step) * parameters_.primal_feasibility_tolerance() >
1.0) {
refactorize = true;
}
++num_iterations_;
@@ -2770,6 +2776,7 @@ Fractional RevisedSimplex::ComputeInitialProblemObjectiveValue() const {
void RevisedSimplex::SetParameters(const GlopParameters& parameters) {
SCOPED_TIME_STAT(&function_stats_);
random_.Reset(parameters_.random_seed());
initial_parameters_ = parameters;
parameters_ = parameters;
PropagateParameters();
@@ -2788,35 +2795,37 @@ void RevisedSimplex::PropagateParameters() {
}
void RevisedSimplex::DisplayIterationInfo() const {
if (FLAGS_v < 1) return;
const int iter = feasibility_phase_
? num_iterations_
: num_iterations_ - num_feasibility_iterations_;
// Note that in the dual phase II, ComputeObjectiveValue() is also computing
// the dual objective even if it uses the variable values. This is because if
// we modify the bounds to make the problem primal-feasible, we are at the
// optimal and hence the two objectives are the same.
const Fractional objective =
!feasibility_phase_
? ComputeInitialProblemObjectiveValue()
: (parameters_.use_dual_simplex()
? reduced_costs_.ComputeSumOfDualInfeasibilities()
: variable_values_.ComputeSumOfPrimalInfeasibilities());
VLOG(1) << (feasibility_phase_ ? "Feasibility" : "Optimization")
<< " phase, iteration # " << iter
<< ", objective = " << StringPrintf("%.15E", objective);
if (VLOG_IS_ON(1)) {
const int iter = feasibility_phase_
? num_iterations_
: num_iterations_ - num_feasibility_iterations_;
// Note that in the dual phase II, ComputeObjectiveValue() is also computing
// the dual objective even if it uses the variable values. This is because
// if we modify the bounds to make the problem primal-feasible, we are at
// the optimal and hence the two objectives are the same.
const Fractional objective =
!feasibility_phase_
? ComputeInitialProblemObjectiveValue()
: (parameters_.use_dual_simplex()
? reduced_costs_.ComputeSumOfDualInfeasibilities()
: variable_values_.ComputeSumOfPrimalInfeasibilities());
VLOG(1) << (feasibility_phase_ ? "Feasibility" : "Optimization")
<< " phase, iteration # " << iter
<< ", objective = " << StringPrintf("%.15E", objective);
}
}
void RevisedSimplex::DisplayErrors() const {
if (FLAGS_v < 1) return;
VLOG(1) << "Primal infeasibility (bounds) = "
<< variable_values_.ComputeMaximumPrimalInfeasibility();
VLOG(1) << "Primal residual |A.x - b| = "
<< variable_values_.ComputeMaximumPrimalResidual();
VLOG(1) << "Dual infeasibility (reduced costs) = "
<< reduced_costs_.ComputeMaximumDualInfeasibility();
VLOG(1) << "Dual residual |c_B - y.B| = "
<< reduced_costs_.ComputeMaximumDualResidual();
if (VLOG_IS_ON(1)) {
VLOG(1) << "Primal infeasibility (bounds) = "
<< variable_values_.ComputeMaximumPrimalInfeasibility();
VLOG(1) << "Primal residual |A.x - b| = "
<< variable_values_.ComputeMaximumPrimalResidual();
VLOG(1) << "Dual infeasibility (reduced costs) = "
<< reduced_costs_.ComputeMaximumDualInfeasibility();
VLOG(1) << "Dual residual |c_B - y.B| = "
<< reduced_costs_.ComputeMaximumDualResidual();
}
}
namespace {
@@ -2848,144 +2857,129 @@ std::string RevisedSimplex::SimpleVariableInfo(ColIndex col) const {
}
void RevisedSimplex::DisplayInfoOnVariables() const {
if (FLAGS_v < 3) return;
for (ColIndex col(0); col < num_cols_; ++col) {
const Fractional variable_value = variable_values_.Get(col);
const Fractional objective_coefficient = current_objective_[col];
const Fractional objective_contribution =
objective_coefficient * variable_value;
VLOG(3) << SimpleVariableInfo(col) << ". " << variable_name_[col] << " = "
<< StringifyWithFlags(variable_value) << " * "
<< StringifyWithFlags(objective_coefficient)
<< "(obj) = " << StringifyWithFlags(objective_contribution);
if (VLOG_IS_ON(3)) {
for (ColIndex col(0); col < num_cols_; ++col) {
const Fractional variable_value = variable_values_.Get(col);
const Fractional objective_coefficient = current_objective_[col];
const Fractional objective_contribution =
objective_coefficient * variable_value;
VLOG(3) << SimpleVariableInfo(col) << ". " << variable_name_[col] << " = "
<< StringifyWithFlags(variable_value) << " * "
<< StringifyWithFlags(objective_coefficient)
<< "(obj) = " << StringifyWithFlags(objective_contribution);
}
VLOG(3) << "------";
}
VLOG(3) << "------";
}
void RevisedSimplex::DisplayVariableBounds() {
if (FLAGS_v < 3) return;
const VariableTypeRow& variable_type = variables_info_.GetTypeRow();
for (ColIndex col(0); col < num_cols_; ++col) {
switch (variable_type[col]) {
case VariableType::UNCONSTRAINED:
break;
case VariableType::LOWER_BOUNDED:
VLOG(3) << variable_name_[col]
<< " >= " << StringifyWithFlags(lower_bound_[col]) << ";";
break;
case VariableType::UPPER_BOUNDED:
VLOG(3) << variable_name_[col]
<< " <= " << StringifyWithFlags(upper_bound_[col]) << ";";
break;
case VariableType::UPPER_AND_LOWER_BOUNDED:
VLOG(3) << StringifyWithFlags(lower_bound_[col])
<< " <= " << variable_name_[col]
<< " <= " << StringifyWithFlags(upper_bound_[col]) << ";";
break;
case VariableType::FIXED_VARIABLE:
VLOG(3) << variable_name_[col] << " = "
<< StringifyWithFlags(lower_bound_[col]) << ";";
break;
default: // This should never happen.
LOG(DFATAL) << "Column " << col
<< " has no meaningful status.";
break;
if (VLOG_IS_ON(3)) {
const VariableTypeRow& variable_type = variables_info_.GetTypeRow();
for (ColIndex col(0); col < num_cols_; ++col) {
switch (variable_type[col]) {
case VariableType::UNCONSTRAINED:
break;
case VariableType::LOWER_BOUNDED:
VLOG(3) << variable_name_[col]
<< " >= " << StringifyWithFlags(lower_bound_[col]) << ";";
break;
case VariableType::UPPER_BOUNDED:
VLOG(3) << variable_name_[col]
<< " <= " << StringifyWithFlags(upper_bound_[col]) << ";";
break;
case VariableType::UPPER_AND_LOWER_BOUNDED:
VLOG(3) << StringifyWithFlags(lower_bound_[col])
<< " <= " << variable_name_[col]
<< " <= " << StringifyWithFlags(upper_bound_[col]) << ";";
break;
case VariableType::FIXED_VARIABLE:
VLOG(3) << variable_name_[col] << " = "
<< StringifyWithFlags(lower_bound_[col]) << ";";
break;
default: // This should never happen.
LOG(DFATAL) << "Column " << col
<< " has no meaningful status.";
break;
}
}
}
}
namespace {
struct MatrixElementForDisplayDictionary {
MatrixElementForDisplayDictionary(ColIndex _col, RowIndex _row,
Fractional _value)
: col(_col), row(_row), value(_value) {}
ColIndex col;
RowIndex row;
Fractional value;
};
Fractional DictionaryLookUpValue(
const std::vector<MatrixElementForDisplayDictionary>& matrix, ColIndex col,
RowIndex row) {
const size_t num_elements = matrix.size();
for (int i = 0; i < num_elements; ++i) {
if (matrix[i].col == col && matrix[i].row == row) {
return matrix[i].value;
}
}
return Fractional(0.0);
}
} // anonymous namespace
void RevisedSimplex::DisplayDictionary() {
// This function has a complexity in
// O(num_cols_ * num_rows_ * non zeros in column).
// It uses LookUpValue(row, col) which has a time complexity
// in O(non zeros in column).
// This choice was made because it is not expected that people use
// this function on large problems and/or very dense problems.
// TODO(user): re-evaluate this design decision if need be.
if (FLAGS_v < 4) return;
std::vector<MatrixElementForDisplayDictionary> dictionary;
DisplayInfoOnVariables();
ITIVector<RowIndex, SparseRow> RevisedSimplex::ComputeDictionary() {
ITIVector<RowIndex, SparseRow> dictionary(num_rows_.value());
for (ColIndex col(0); col < num_cols_; ++col) {
ComputeDirection(col);
for (const RowIndex row : direction_non_zero_) {
dictionary.push_back(
MatrixElementForDisplayDictionary(col, row, direction_[row]));
dictionary[row].SetCoefficient(col, direction_[row]);
}
}
return dictionary;
}
std::string output = "z = " + StringifyWithFlags(ComputeObjectiveValue());
const DenseRow& reduced_costs = reduced_costs_.GetReducedCosts();
for (const ColIndex col : variables_info_.GetNotBasicBitRow()) {
output +=
StringifyMonomialWithFlags(reduced_costs[col], variable_name_[col]);
}
VLOG(3) << output << ";";
for (RowIndex row(0); row < num_rows_; ++row) {
output = "";
ColIndex basic_col = basis_[row];
output += variable_name_[basic_col] + " = " +
StringifyWithFlags(variable_values_.Get(basic_col));
for (ColIndex col(0); col < num_cols_; ++col) {
if (col != basic_col) {
const Fractional value = DictionaryLookUpValue(dictionary, col, row);
output += StringifyMonomialWithFlags(value, variable_name_[col]);
}
void RevisedSimplex::DisplayRevisedSimplexDebugInfo() {
if (VLOG_IS_ON(3)) {
// This function has a complexity in O(num_non_zeros_in_matrix).
DisplayInfoOnVariables();
std::string output = "z = " + StringifyWithFlags(ComputeObjectiveValue());
const DenseRow& reduced_costs = reduced_costs_.GetReducedCosts();
for (const ColIndex col : variables_info_.GetNotBasicBitRow()) {
StrAppend(&output, StringifyMonomialWithFlags(reduced_costs[col],
variable_name_[col]));
}
VLOG(3) << output << ";";
const RevisedSimplexDictionary dictionary(this);
RowIndex r(0);
for (const SparseRow& row : dictionary) {
output.clear();
ColIndex basic_col = basis_[r];
StrAppend(&output, variable_name_[basic_col], " = ",
StringifyWithFlags(variable_values_.Get(basic_col)));
for (const SparseRowEntry e : row) {
if (e.col() != basic_col) {
StrAppend(&output,
StringifyMonomialWithFlags(e.coefficient(),
variable_name_[e.col()]));
}
}
VLOG(3) << output << ";";
}
VLOG(3) << "------";
DisplayVariableBounds();
++r;
}
VLOG(3) << "------";
DisplayVariableBounds();
}
void RevisedSimplex::DisplayProblem() const {
// TODO(user): see comment for DisplayDictionary()
if (FLAGS_v < 3) return;
DisplayInfoOnVariables();
std::string output = "min: ";
bool has_objective = false;
for (ColIndex col(0); col < num_cols_; ++col) {
const Fractional coeff = current_objective_[col];
has_objective |= (coeff != 0.0);
output += StringifyMonomialWithFlags(coeff, variable_name_[col]);
}
if (!has_objective) {
output += " 0";
}
VLOG(3) << output << ";";
for (RowIndex row(0); row < num_rows_; ++row) {
output = "";
// This function has a complexity in O(num_rows * num_cols *
// num_non_zeros_in_row).
if (VLOG_IS_ON(3)) {
DisplayInfoOnVariables();
std::string output = "min: ";
bool has_objective = false;
for (ColIndex col(0); col < num_cols_; ++col) {
output += StringifyMonomialWithFlags(
matrix_with_slack_.column(col).LookUpCoefficient(row),
variable_name_[col]);
const Fractional coeff = current_objective_[col];
has_objective |= (coeff != 0.0);
StrAppend(&output,
StringifyMonomialWithFlags(coeff, variable_name_[col]));
}
VLOG(3) << output << " = 0;";
if (!has_objective) {
StrAppend(&output, " 0");
}
VLOG(3) << output << ";";
for (RowIndex row(0); row < num_rows_; ++row) {
output = "";
for (ColIndex col(0); col < num_cols_; ++col) {
StrAppend(
&output, StringifyMonomialWithFlags(
matrix_with_slack_.column(col).LookUpCoefficient(row),
variable_name_[col]));
}
VLOG(3) << output << " = 0;";
}
VLOG(3) << "------";
}
VLOG(3) << "------";
}
void RevisedSimplex::AdvanceDeterministicTime(TimeLimit* time_limit) {

47
src/glop/revised_simplex.h
View File

@@ -114,6 +114,7 @@
#include "lp_data/lp_print_utils.h"
#include "lp_data/lp_types.h"
#include "lp_data/matrix_scaler.h"
#include "lp_data/sparse_row.h"
#include "base/random.h"
#include "util/time_limit.h"
@@ -206,6 +207,9 @@ class RevisedSimplex {
const DenseRow& GetPrimalRay() const;
const DenseColumn& GetDualRay() const;
// This is the "dual ray" linear combination of the matrix rows.
const DenseRow& GetDualRayRowCombination() const;
// Returns the index of the column in the basis and the basis factorization.
// Note that the order of the column in the basis is important since it is the
// one used by the various solve functions provided by the BasisFactorization
@@ -216,6 +220,11 @@ class RevisedSimplex {
// Returns statistics about this class as a std::string.
std::string StatString();
// Computes the dictionary B^-1*N on-the-fly row by row. Returns the resulting
// matrix as a vector of sparse rows so that it is easy to use it on the left
// side in the matrix multiplication. Runs in O(num_non_zeros_in_matrix).
RowMajorSparseMatrix ComputeDictionary();
private:
// Propagates parameters_ to all the other classes that need it.
//
@@ -254,8 +263,12 @@ class RevisedSimplex {
// Displays the bounds of the variables.
void DisplayVariableBounds();
// Displays the Linear Programming problem as a dictionary,
// taking into account the iterations that have been made.
// Displays the following information:
// * Linear Programming problem as a dictionary, taking into
// account the iterations that have been made;
// * Variable info;
// * Reduced costs;
// * Variable bounds.
// A dictionary is in the form:
// xB = value + sum_{j in N} pa_ij x_j
// z = objective_value + sum_{i in N} rc_i x_i
@@ -264,7 +277,7 @@ class RevisedSimplex {
// after the pivotings.
// Dictionaries are the modern way of presenting the result of an iteration
// of the Simplex algorithm in the literature.
void DisplayDictionary();
void DisplayRevisedSimplexDebugInfo();
// Displays the Linear Programming problem as it was input.
void DisplayProblem() const;
@@ -385,7 +398,7 @@ class RevisedSimplex {
// See Chvatal's book for more detail (Chapter 7).
void ComputeDirection(ColIndex col);
// Computes a - B.d in error_ and return the maximum fabs() of its coeffs.
// Computes a - B.d in error_ and return the maximum std::abs() of its coeffs.
// TODO(user): Use this to trigger a refactorization of B? Or to track the
// error created by calls to SkipVariableForRatioTest()?
Fractional ComputeDirectionError(ColIndex col);
@@ -633,6 +646,7 @@ class RevisedSimplex {
DenseRow solution_reduced_costs_;
DenseRow solution_primal_ray_;
DenseColumn solution_dual_ray_;
DenseRow solution_dual_ray_row_combination_;
BasisState solution_state_;
bool solution_state_has_been_set_externally_;
@@ -771,6 +785,31 @@ class RevisedSimplex {
DISALLOW_COPY_AND_ASSIGN(RevisedSimplex);
};
// Hides the details of the dictionary matrix implementation. In the future,
// GLOP will support generating the dictionary one row at a time without having
// to store the whole matrix in memory.
class RevisedSimplexDictionary {
public:
typedef RowMajorSparseMatrix::const_iterator ConstIterator;
// RevisedSimplex cannot be passed const because we have to call a non-const
// method ComputeDictionary.
explicit RevisedSimplexDictionary(RevisedSimplex* revised_simplex)
: dictionary_(CHECK_NOTNULL(revised_simplex)->ComputeDictionary()) {}
ConstIterator begin() const { return dictionary_.begin(); }
ConstIterator end() const { return dictionary_.end(); }
size_t size() const { return dictionary_.size(); }
private:
const RowMajorSparseMatrix dictionary_;
DISALLOW_COPY_AND_ASSIGN(RevisedSimplexDictionary);
};
// TODO(user): When a row-by-row generation of the dictionary is supported,
// implement DictionaryIterator class that would call it inside operator*().
} // namespace glop
} // namespace operations_research

33
src/glop/update_row.cc
View File

@@ -162,7 +162,7 @@ void UpdateRow::ComputeUpdatesRowWise() {
non_zero_position_list_.clear();
const Fractional drop_tolerance = parameters_.drop_tolerance();
for (const ColIndex col : variables_info_.GetIsRelevantBitRow()) {
if (fabs(coefficient_[col]) > drop_tolerance) {
if (std::abs(coefficient_[col]) > drop_tolerance) {
non_zero_position_list_.push_back(col);
}
}
@@ -202,12 +202,37 @@ void UpdateRow::ComputeUpdatesRowWiseHypersparse() {
// non-zero coefficients contiguously in a vector is better than keeping
// them as they are. Note however that we will iterate only twice on the
// update row coefficients during an iteration.
if (fabs(coefficient_[col]) > drop_tolerance) {
if (std::abs(coefficient_[col]) > drop_tolerance) {
non_zero_position_list_.push_back(col);
}
}
}
// Note that we use the same algo as ComputeUpdatesColumnWise() here. The
// others version might be faster, but this is called only once per solve, so
// it shouldn't be too bad.
void UpdateRow::RecomputeFullUpdateRow(RowIndex leaving_row) {
CHECK(!compute_update_row_);
const ColIndex num_cols = matrix_.num_cols();
const Fractional drop_tolerance = parameters_.drop_tolerance();
coefficient_.resize(num_cols, 0.0);
non_zero_position_list_.clear();
// Fills the only position at one in the basic columns.
coefficient_[basis_[leaving_row]] = 1.0;
non_zero_position_list_.push_back(basis_[leaving_row]);
// Fills the non-basic column.
for (const ColIndex col : variables_info_.GetNotBasicBitRow()) {
const Fractional coeff =
matrix_.ColumnScalarProduct(col, unit_row_left_inverse_);
if (std::abs(coeff) > drop_tolerance) {
non_zero_position_list_.push_back(col);
coefficient_[col] = coeff;
}
}
}
void UpdateRow::ComputeUpdatesColumnWise() {
SCOPED_TIME_STAT(&stats_);
@@ -227,7 +252,7 @@ void UpdateRow::ComputeUpdatesColumnWise() {
// sparsity. Note that it shouldn't be too bad to use a non-zero drop
// tolerance here because even if we introduce some precision issues, the
// quantities updated by this update row will eventually be recomputed.
if (fabs(coeff) > drop_tolerance) {
if (std::abs(coeff) > drop_tolerance) {
non_zero_position_list_.push_back(col);
coefficient_[col] = coeff;
}
@@ -248,7 +273,7 @@ void UpdateRow::ComputeUpdatesColumnWise() {
}
// End of omp parallel for.
for (const ColIndex col : variables_info_.GetIsRelevantBitRow()) {
if (fabs(coefficient_[col]) > drop_tolerance) {
if (std::abs(coefficient_[col]) > drop_tolerance) {
non_zero_position_list_.push_back(col);
}
}

4
src/glop/update_row.h
View File

@@ -66,6 +66,10 @@ class UpdateRow {
return coefficient_[col];
}
// This must be called after a call to ComputeUpdateRow(). It will fill
// all the non-relevant positions that where not filled by ComputeUpdateRow().
void RecomputeFullUpdateRow(RowIndex leaving_row);
// Sets to zero the coefficient for column col.
void IgnoreUpdatePosition(ColIndex col);

70
src/glop/variables_info.cc
View File

@@ -32,7 +32,7 @@ void VariablesInfo::Initialize() {
variable_status_.resize(num_cols, VariableStatus::FREE);
can_increase_.ClearAndResize(num_cols);
can_decrease_.ClearAndResize(num_cols);
is_relevant_.ClearAndResize(num_cols);
relevance_.ClearAndResize(num_cols);
is_basic_.ClearAndResize(num_cols);
not_basic_.ClearAndResize(num_cols);
non_basic_boxed_variables_.ClearAndResize(num_cols);
@@ -40,24 +40,56 @@ void VariablesInfo::Initialize() {
void VariablesInfo::MakeBoxedVariableRelevant(bool value) {
if (value == boxed_variables_are_relevant_) return;
boxed_variables_are_relevant_ = value;
for (ColIndex col(0); col < is_relevant_.size(); ++col) {
SetRelevance(col, variable_status_[col]);
if (value) {
boxed_variables_are_relevant_ = true;
for (const ColIndex col : non_basic_boxed_variables_) {
SetRelevance(col, variable_status_[col] != VariableStatus::FIXED_VALUE);
}
} else {
boxed_variables_are_relevant_ = false;
for (const ColIndex col : non_basic_boxed_variables_) {
SetRelevance(col, false);
}
}
}
void VariablesInfo::Update(ColIndex col, VariableStatus status) {
if (status == VariableStatus::BASIC) {
UpdateToBasicStatus(col);
} else {
UpdateToNonBasicStatus(col, status);
}
}
void VariablesInfo::UpdateToBasicStatus(ColIndex col) {
variable_type_[col] = ComputeVariableType(col);
variable_status_[col] = VariableStatus::BASIC;
is_basic_.Set(col, true);
not_basic_.Set(col, false);
can_increase_.Set(col, false);
can_decrease_.Set(col, false);
non_basic_boxed_variables_.Set(col, false);
SetRelevance(col, false);
}
void VariablesInfo::UpdateToNonBasicStatus(ColIndex col,
VariableStatus status) {
DCHECK_NE(status, VariableStatus::BASIC);
const VariableType type = ComputeVariableType(col);
variable_type_[col] = type;
variable_status_[col] = status;
is_basic_.Set(col, status == VariableStatus::BASIC);
not_basic_.Set(col, status != VariableStatus::BASIC);
is_basic_.Set(col, false);
not_basic_.Set(col, true);
can_increase_.Set(col, status == VariableStatus::AT_LOWER_BOUND ||
status == VariableStatus::FREE);
status == VariableStatus::FREE);
can_decrease_.Set(col, status == VariableStatus::AT_UPPER_BOUND ||
status == VariableStatus::FREE);
SetRelevance(col, status);
non_basic_boxed_variables_.Set(col, status != VariableStatus::BASIC &&
variable_type_[col] == VariableType::UPPER_AND_LOWER_BOUNDED);
status == VariableStatus::FREE);
non_basic_boxed_variables_.Set(col,
type == VariableType::UPPER_AND_LOWER_BOUNDED);
const bool relevance = status != VariableStatus::FIXED_VALUE &&
(boxed_variables_are_relevant_ ||
type != VariableType::UPPER_AND_LOWER_BOUNDED);
SetRelevance(col, relevance);
}
const VariableTypeRow& VariablesInfo::GetTypeRow() const {
@@ -77,7 +109,7 @@ const DenseBitRow& VariablesInfo::GetCanDecreaseBitRow() const {
}
const DenseBitRow& VariablesInfo::GetIsRelevantBitRow() const {
return is_relevant_;
return relevance_;
}
const DenseBitRow& VariablesInfo::GetIsBasicBitRow() const { return is_basic_; }
@@ -109,17 +141,13 @@ VariableType VariablesInfo::ComputeVariableType(ColIndex col) const {
}
}
void VariablesInfo::SetRelevance(ColIndex col, VariableStatus status) {
const bool is_relevant = (status != VariableStatus::BASIC &&
status != VariableStatus::FIXED_VALUE &&
(boxed_variables_are_relevant_ ||
variable_type_[col] != VariableType::UPPER_AND_LOWER_BOUNDED));
if (is_relevant_.IsSet(col) == is_relevant) return;
if (is_relevant) {
is_relevant_.Set(col);
void VariablesInfo::SetRelevance(ColIndex col, bool relevance) {
if (relevance_.IsSet(col) == relevance) return;
if (relevance) {
relevance_.Set(col);
num_entries_in_relevant_columns_ += matrix_.ColumnNumEntries(col);
} else {
is_relevant_.Clear(col);
relevance_.Clear(col);
num_entries_in_relevant_columns_ -= matrix_.ColumnNumEntries(col);
}
}

10
src/glop/variables_info.h
View File

@@ -41,6 +41,10 @@ class VariablesInfo {
// to call this if the status or the bound of a variable didn't change.
void Update(ColIndex col, VariableStatus status);
// Slighlty optimized version of Update() above for the two separate cases.
void UpdateToBasicStatus(ColIndex col);
void UpdateToNonBasicStatus(ColIndex col, VariableStatus status);
// Various getter, see the corresponding member declaration below for more
// information.
const VariableTypeRow& GetTypeRow() const;
@@ -76,7 +80,7 @@ class VariablesInfo {
VariableType ComputeVariableType(ColIndex col) const;
// Sets the column relevance and updates num_entries_in_relevant_columns_.
void SetRelevance(ColIndex col, VariableStatus status);
void SetRelevance(ColIndex col, bool relevance);
// Problem data that should be updated from outside.
const CompactSparseMatrix& matrix_;
@@ -98,7 +102,7 @@ class VariablesInfo {
// Indicates if we should consider this variable for entering the basis during
// the simplex algorithm. Only non-fixed and non-basic columns are relevant.
DenseBitRow is_relevant_;
DenseBitRow relevance_;
// Indicates if a variable is BASIC or not. There are currently two members
// because the DenseBitRow class only supports a nice range-based iteration on
@@ -109,7 +113,7 @@ class VariablesInfo {
// Set of boxed variables that are non-basic.
DenseBitRow non_basic_boxed_variables_;
// Number of entries for the relevant matrix columns (see is_relevant_).
// Number of entries for the relevant matrix columns (see relevance_).
EntryIndex num_entries_in_relevant_columns_;
// Whether or not a boxed variable should be considered relevant.