20#include "absl/log/check.h"
21#include "absl/log/log.h"
22#include "absl/random/bit_gen_ref.h"
43 absl::BitGenRef random)
45 objective_(objective),
47 variables_info_(variables_info),
48 basis_factorization_(basis_factorization),
52 must_refactorize_basis_(false),
53 recompute_basic_objective_left_inverse_(true),
54 recompute_basic_objective_(true),
55 recompute_reduced_costs_(true),
56 are_reduced_costs_precise_(false),
57 are_reduced_costs_recomputed_(false),
60 basic_objective_left_inverse_(),
61 dual_feasibility_tolerance_() {}
64 return must_refactorize_basis_;
70 if (recompute_basic_objective_) {
71 ComputeBasicObjective();
73 const Fractional old_reduced_cost = reduced_costs_[entering_col];
75 objective_[entering_col] + cost_perturbations_[entering_col] -
79 reduced_costs_[entering_col] = precise_reduced_cost;
87 if (!recompute_reduced_costs_) {
88 const Fractional estimated_reduced_costs_accuracy =
89 old_reduced_cost - precise_reduced_cost;
91 (std::abs(precise_reduced_cost) <= 1.0) ? 1.0 : precise_reduced_cost;
92 stats_.reduced_costs_accuracy.Add(estimated_reduced_costs_accuracy / scale);
93 if (std::abs(estimated_reduced_costs_accuracy) / scale >
94 parameters_.recompute_reduced_costs_threshold()) {
95 VLOG(1) <<
"Recomputing reduced costs, value = " << precise_reduced_cost
97 << std::abs(precise_reduced_cost - old_reduced_cost);
102 return precise_reduced_cost;
108 const RowIndex num_rows = matrix_.num_rows();
110 for (RowIndex row(0); row < num_rows; ++row) {
111 const ColIndex basic_col = basis_[row];
113 objective_[basic_col] + cost_perturbations_[basic_col] -
114 matrix_.ColumnScalarProduct(basic_col, dual_values);
115 dual_residual_error = std::max(dual_residual_error, std::abs(residual));
117 return dual_residual_error;
127 const DenseBitRow& can_decrease = variables_info_.GetCanDecreaseBitRow();
128 const DenseBitRow& can_increase = variables_info_.GetCanIncreaseBitRow();
129 for (
const ColIndex col : variables_info_.GetIsRelevantBitRow()) {
131 if ((can_increase.
IsSet(col) && rc < 0.0) ||
132 (can_decrease.
IsSet(col) && rc > 0.0)) {
133 maximum_dual_infeasibility =
134 std::max(maximum_dual_infeasibility, std::abs(rc));
137 return maximum_dual_infeasibility;
147 const DenseBitRow& can_decrease = variables_info_.GetCanDecreaseBitRow();
148 const DenseBitRow& can_increase = variables_info_.GetCanIncreaseBitRow();
149 const DenseBitRow& is_boxed = variables_info_.GetNonBasicBoxedVariables();
150 for (
const ColIndex col : variables_info_.GetNotBasicBitRow()) {
151 if (is_boxed[col])
continue;
153 if ((can_increase.
IsSet(col) && rc < 0.0) ||
154 (can_decrease.
IsSet(col) && rc > 0.0)) {
155 maximum_dual_infeasibility =
156 std::max(maximum_dual_infeasibility, std::abs(rc));
159 return maximum_dual_infeasibility;
169 const DenseBitRow& can_decrease = variables_info_.GetCanDecreaseBitRow();
170 const DenseBitRow& can_increase = variables_info_.GetCanIncreaseBitRow();
171 for (
const ColIndex col : variables_info_.GetIsRelevantBitRow()) {
173 if ((can_increase.
IsSet(col) && rc < 0.0) ||
174 (can_decrease.
IsSet(col) && rc > 0.0)) {
175 dual_infeasibility_sum += std::abs(std::abs(rc));
178 return dual_infeasibility_sum;
182 RowIndex leaving_row,
186 const ColIndex leaving_col = basis_[leaving_row];
187 DCHECK(!variables_info_.GetIsBasicBitRow().IsSet(entering_col));
188 DCHECK(variables_info_.GetIsBasicBitRow().IsSet(leaving_col));
191 if (!recompute_reduced_costs_) {
192 UpdateReducedCosts(entering_col, leaving_col, leaving_row,
193 direction[leaving_row], update_row);
198 UpdateBasicObjective(entering_col, leaving_row);
204 DCHECK_EQ(current_cost, &objective_[col]);
205 reduced_costs_[col] -= objective_[col];
210 parameters_ = parameters;
215 recompute_basic_objective_ =
true;
216 recompute_basic_objective_left_inverse_ =
true;
217 are_reduced_costs_precise_ =
false;
218 SetRecomputeReducedCostsAndNotifyWatchers();
223 recompute_basic_objective_ =
true;
224 recompute_basic_objective_left_inverse_ =
true;
229 if (are_reduced_costs_precise_)
return;
230 must_refactorize_basis_ =
true;
231 recompute_basic_objective_left_inverse_ =
true;
232 SetRecomputeReducedCostsAndNotifyWatchers();
237 VLOG(1) <<
"Perturbing the costs ... ";
240 const ColIndex structural_size =
242 for (ColIndex col(0); col < structural_size; ++col) {
244 std::max(max_cost_magnitude, std::abs(objective_[col]));
247 cost_perturbations_.AssignToZero(matrix_.num_cols());
248 for (ColIndex col(0); col < structural_size; ++col) {
251 (1.0 + std::uniform_real_distribution<double>()(random_)) *
252 (parameters_.relative_cost_perturbation() * std::abs(objective) +
253 parameters_.relative_max_cost_perturbation() * max_cost_magnitude);
254 DCHECK_GE(magnitude, 0.0);
258 const VariableType type = variables_info_.GetTypeRow()[col];
265 cost_perturbations_[col] = magnitude;
268 cost_perturbations_[col] = -magnitude;
276 if (objective > 0.0) {
277 cost_perturbations_[col] = magnitude;
278 }
else if (objective < 0.0) {
279 cost_perturbations_[col] = -magnitude;
293 parameters_.degenerate_ministep_factor() * dual_feasibility_tolerance_;
294 if (increasing_rc_is_needed && reduced_costs_[col] <= -minimum_delta)
return;
295 if (!increasing_rc_is_needed && reduced_costs_[col] >= minimum_delta)
return;
298 increasing_rc_is_needed ? minimum_delta : -minimum_delta;
300 cost_perturbations_[col] -= reduced_costs_[col] + delta;
301 reduced_costs_[col] = -delta;
302 has_cost_shift_ =
true;
307 if (increasing_rc_is_needed && reduced_costs_[col] >= 0.0)
return true;
308 if (!increasing_rc_is_needed && reduced_costs_[col] <= 0.0)
return true;
314 has_cost_shift_ =
false;
315 cost_perturbations_.AssignToZero(matrix_.num_cols());
316 recompute_basic_objective_ =
true;
317 recompute_basic_objective_left_inverse_ =
true;
318 are_reduced_costs_precise_ =
false;
319 SetRecomputeReducedCostsAndNotifyWatchers();
324 if (!are_reduced_costs_recomputed_) {
325 SetRecomputeReducedCostsAndNotifyWatchers();
332 if (basis_factorization_.IsRefactorized()) {
333 must_refactorize_basis_ =
false;
335 if (recompute_reduced_costs_) {
336 ComputeReducedCosts();
338 return reduced_costs_.const_view();
343 ComputeBasicObjectiveLeftInverse();
344 return Transpose(basic_objective_left_inverse_.values);
347void ReducedCosts::ComputeBasicObjective() {
351 basic_objective_.
resize(num_cols_in_basis, 0.0);
352 for (ColIndex col(0); col < num_cols_in_basis; ++col) {
354 basic_objective_[col] =
355 objective_[basis_col] + cost_perturbations_[basis_col];
357 recompute_basic_objective_ =
false;
358 recompute_basic_objective_left_inverse_ =
true;
361void ReducedCosts::ComputeReducedCosts() {
363 if (recompute_basic_objective_left_inverse_) {
364 ComputeBasicObjectiveLeftInverse();
367 const ColIndex num_cols = matrix_.num_cols();
368 const ColIndex first_slack = num_cols -
RowToColIndex(matrix_.num_rows());
370 reduced_costs_.resize(num_cols, 0.0);
371 const DenseBitRow& is_basic = variables_info_.GetIsBasicBitRow();
372 for (ColIndex col(0); col < first_slack; ++col) {
373 reduced_costs_[col] =
374 objective_[col] + cost_perturbations_[col] -
375 matrix_.ColumnScalarProduct(col, basic_objective_left_inverse_.values);
378 if (is_basic.IsSet(col)) {
379 dual_residual_error =
380 std::max(dual_residual_error, std::abs(reduced_costs_[col]));
383 for (ColIndex col(first_slack); col < num_cols; ++col) {
384 reduced_costs_[col] = objective_[col] + cost_perturbations_[col] -
385 basic_objective_left_inverse_[col - first_slack];
388 if (is_basic.IsSet(col)) {
389 dual_residual_error =
390 std::max(dual_residual_error, std::abs(reduced_costs_[col]));
394 deterministic_time_ +=
396 recompute_reduced_costs_ =
false;
397 are_reduced_costs_recomputed_ =
true;
398 are_reduced_costs_precise_ = basis_factorization_.IsRefactorized();
405 dual_feasibility_tolerance_ = parameters_.dual_feasibility_tolerance();
406 if (dual_residual_error > dual_feasibility_tolerance_) {
407 VLOG(2) <<
"Changing dual_feasibility_tolerance to " << dual_residual_error;
408 dual_feasibility_tolerance_ = dual_residual_error;
412void ReducedCosts::ComputeBasicObjectiveLeftInverse() {
414 if (recompute_basic_objective_) {
415 ComputeBasicObjective();
417 basic_objective_left_inverse_.values = basic_objective_;
418 basic_objective_left_inverse_.non_zeros.clear();
419 basis_factorization_.LeftSolve(&basic_objective_left_inverse_);
420 recompute_basic_objective_left_inverse_ =
false;
422 Density(basic_objective_left_inverse_.values)));
431void ReducedCosts::UpdateReducedCosts(ColIndex entering_col,
432 ColIndex leaving_col,
435 DCHECK_NE(entering_col, leaving_col);
436 DCHECK_NE(pivot, 0.0);
437 if (recompute_reduced_costs_)
return;
440 const Fractional entering_reduced_cost = reduced_costs_[entering_col];
444 if (entering_reduced_cost == 0.0) {
445 VLOG(2) <<
"Reduced costs didn't change.";
451 are_reduced_costs_precise_ =
false;
455 are_reduced_costs_recomputed_ =
false;
456 are_reduced_costs_precise_ =
false;
457 update_row->ComputeUpdateRow(leaving_row);
464 const Fractional new_leaving_reduced_cost = entering_reduced_cost / -pivot;
465 auto rc = reduced_costs_.view();
466 auto update_coeffs = update_row->GetCoefficients().const_view();
467 for (
const ColIndex col : update_row->GetNonZeroPositions()) {
468 rc[col] += new_leaving_reduced_cost * update_coeffs[col];
470 rc[leaving_col] = new_leaving_reduced_cost;
475 rc[entering_col] = 0.0;
479 const Fractional reduced_cost = reduced_costs_[col];
480 const DenseBitRow& can_decrease = variables_info_.GetCanDecreaseBitRow();
481 const DenseBitRow& can_increase = variables_info_.GetCanIncreaseBitRow();
482 const Fractional tolerance = dual_feasibility_tolerance_;
483 return (can_increase.
IsSet(col) && (reduced_cost < -tolerance)) ||
484 (can_decrease.
IsSet(col) && (reduced_cost > tolerance));
487void ReducedCosts::UpdateBasicObjective(ColIndex entering_col,
488 RowIndex leaving_row) {
491 objective_[entering_col] + cost_perturbations_[entering_col];
492 recompute_basic_objective_left_inverse_ =
true;
495void ReducedCosts::SetRecomputeReducedCostsAndNotifyWatchers() {
496 recompute_reduced_costs_ =
true;
497 for (
bool* watcher : watchers_) *watcher =
true;
505 variables_info_(variables_info),
506 primal_edge_norms_(primal_edge_norms),
507 reduced_costs_(reduced_costs) {
508 reduced_costs_->AddRecomputationWatcher(&recompute_);
515 if (recompute_)
return;
520 UpdateEnteringCandidates<
false>(
525 prices_.Remove(entering_col);
529 if (recompute_)
return;
530 if (reduced_costs_->IsValidPrimalEnteringCandidate(col)) {
532 primal_edge_norms_->GetSquaredNorms();
534 DCHECK_NE(0.0, squared_norms[col]);
535 prices_.AddOrUpdate(col,
Square(reduced_costs[col]) / squared_norms[col]);
544 if (recompute_)
return;
546 DCHECK(!reduced_costs_->IsValidPrimalEnteringCandidate(col));
553 prices_.ClearAndResize(reduced_costs.
size());
554 UpdateEnteringCandidates<
true>(
555 variables_info_.GetIsRelevantBitRow());
558 return prices_.GetMaximum();
566template <
bool from_clean_state,
typename ColumnsToUpdate>
567void PrimalPrices::UpdateEnteringCandidates(
const ColumnsToUpdate& cols) {
569 const DenseBitRow::ConstView can_decrease =
571 const DenseBitRow::ConstView can_increase =
576 for (
const ColIndex col : cols) {
577 const Fractional reduced_cost = reduced_costs[col];
585 col, reduced_cost > tolerance, can_decrease, reduced_cost < -tolerance,
587 if (is_dual_infeasible) {
593 if (!from_clean_state) prices_.
Remove(col);
bool IsSet(IndexType i) const
Returns true if the bit at position i is set.
static uint64_t ConditionalXorOfTwoBits(IndexType i, uint64_t use1, Bitset64< IndexType >::ConstView set1, uint64_t use2, Bitset64< IndexType >::ConstView set2)
ConstView const_view() const
RowIndex num_rows() const
ColIndex num_cols() const
void Remove(Index position)
Removes the given index from the set of candidates.
void AddOrUpdate(Index position, Fractional value)
DenseRow::ConstView GetSquaredNorms()
void AddRecomputationWatcher(bool *watcher)
void RecomputePriceAt(ColIndex col)
Triggers a recomputation of the price at the given column only.
void SetAndDebugCheckThatColumnIsDualFeasible(ColIndex col)
void UpdateBeforeBasisPivot(ColIndex entering_col, UpdateRow *update_row)
PrimalPrices(absl::BitGenRef random, const VariablesInfo &variables_info, PrimalEdgeNorms *primal_edge_norms, ReducedCosts *reduced_costs)
ColIndex GetBestEnteringColumn()
void SetNonBasicVariableCostToZero(ColIndex col, Fractional *current_cost)
Fractional ComputeMaximumDualInfeasibilityOnNonBoxedVariables()
Fractional ComputeMaximumDualInfeasibility()
Fractional GetDualFeasibilityTolerance() const
Returns the current dual feasibility tolerance.
void SetParameters(const GlopParameters ¶meters)
Sets the pricing parameters. This does not change the pricing rule.
bool StepIsDualDegenerate(bool increasing_rc_is_needed, ColIndex col)
bool NeedsBasisRefactorization() const
Fractional ComputeMaximumDualResidual()
Fractional TestEnteringReducedCostPrecision(ColIndex entering_col, const ScatteredColumn &direction)
bool IsValidPrimalEnteringCandidate(ColIndex col) const
Does basic checking of an entering candidate.
void UpdateDataOnBasisPermutation()
Invalidates the data that depends on the order of the column in basis_.
void ClearAndRemoveCostShifts()
ReducedCosts(const CompactSparseMatrix &matrix_, const DenseRow &objective, const RowToColMapping &basis, const VariablesInfo &variables_info, const BasisFactorization &basis_factorization, absl::BitGenRef random)
Takes references to the linear program data we need.
const DenseColumn & GetDualValues()
Returns the dual values associated to the current basis.
void ShiftCostIfNeeded(bool increasing_rc_is_needed, ColIndex col)
void UpdateBeforeBasisPivot(ColIndex entering_col, RowIndex leaving_row, const ScatteredColumn &direction, UpdateRow *update_row)
DenseRow::ConstView GetFullReducedCosts()
DenseRow::ConstView GetReducedCosts()
void MakeReducedCostsPrecise()
Fractional ComputeSumOfDualInfeasibilities()
void ResetForNewObjective()
Invalidates all internal structure that depends on the objective function.
StrictITISpan< ColIndex, const Fractional > ConstView
void resize(IntType size)
absl::Span< const ColIndex > GetNonZeroPositions() const
const DenseBitRow & GetCanDecreaseBitRow() const
const DenseBitRow & GetCanIncreaseBitRow() const
double Density(const DenseRow &row)
StrictITIVector< RowIndex, ColIndex > RowToColMapping
Fractional Square(Fractional f)
Fractional ScalarProduct(const DenseRowOrColumn1 &u, const DenseRowOrColumn2 &v)
Bitset64< ColIndex > DenseBitRow
Row of bits.
const DenseRow & Transpose(const DenseColumn &col)
Similar comment as the other Transpose() implementation above.
ColIndex RowToColIndex(RowIndex row)
Get the ColIndex corresponding to the column # row.
VariableType
Different types of variables.
@ UPPER_AND_LOWER_BOUNDED
RowIndex ColToRowIndex(ColIndex col)
Get the RowIndex corresponding to the row # col.
StrictITIVector< RowIndex, Fractional > DenseColumn
Column-vector types. Column-vector types are indexed by a row index.
StrictITIVector< ColIndex, Fractional > DenseRow
Row-vector types. Row-vector types are indexed by a column index.
static double DeterministicTimeForFpOperations(int64_t n)
In SWIG mode, we don't want anything besides these top-level includes.
#define IF_STATS_ENABLED(instructions)
#define SCOPED_TIME_STAT(stats)