OR-Tools  8.2
entering_variable.cc
Go to the documentation of this file.
1// Copyright 2010-2018 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
15
16#include <queue>
17
18#include "ortools/base/timer.h"
21
22namespace operations_research {
23namespace glop {
24
26 random_engine_t* random,
27 ReducedCosts* reduced_costs,
28 PrimalEdgeNorms* primal_edge_norms)
29 : variables_info_(variables_info),
30 random_(random),
31 reduced_costs_(reduced_costs),
32 primal_edge_norms_(primal_edge_norms),
33 parameters_(),
34 rule_(GlopParameters::DANTZIG),
35 unused_columns_() {}
36
38 SCOPED_TIME_STAT(&stats_);
39 GLOP_RETURN_ERROR_IF_NULL(entering_col);
40
41 // For better redability of the templated function calls below.
42 const bool kNormalize = true;
43 const bool kNested = true;
44 const bool kSteepest = true;
45
46 switch (rule_) {
47 case GlopParameters::DANTZIG:
48 if (parameters_.use_nested_pricing()) {
49 if (unused_columns_.size() != variables_info_.GetNumberOfColumns()) {
51 }
52 if (parameters_.normalize_using_column_norm()) {
53 DantzigChooseEnteringColumn<kNormalize, kNested>(entering_col);
54 } else {
55 DantzigChooseEnteringColumn<!kNormalize, kNested>(entering_col);
56 }
57 if (*entering_col != kInvalidCol) {
58 unused_columns_.Clear(*entering_col);
59 return Status::OK();
60 }
62 if (parameters_.normalize_using_column_norm()) {
63 DantzigChooseEnteringColumn<kNormalize, kNested>(entering_col);
64 } else {
65 DantzigChooseEnteringColumn<!kNormalize, kNested>(entering_col);
66 }
67 } else {
68 if (parameters_.normalize_using_column_norm()) {
69 DantzigChooseEnteringColumn<kNormalize, !kNested>(entering_col);
70 } else {
71 DantzigChooseEnteringColumn<!kNormalize, !kNested>(entering_col);
72 }
73 }
74 return Status::OK();
75 case GlopParameters::STEEPEST_EDGE:
76 NormalizedChooseEnteringColumn<kSteepest>(entering_col);
77 return Status::OK();
78 case GlopParameters::DEVEX:
79 NormalizedChooseEnteringColumn<!kSteepest>(entering_col);
80 return Status::OK();
81 }
82 LOG(DFATAL) << "Unknown pricing rule: "
83 << ProtoEnumToString<GlopParameters::PricingRule>(rule_)
84 << ". Using steepest edge.";
85 NormalizedChooseEnteringColumn<kSteepest>(entering_col);
86 return Status::OK();
87}
88
90 const UpdateRow& update_row, Fractional cost_variation,
91 std::vector<ColIndex>* bound_flip_candidates, ColIndex* entering_col,
92 Fractional* step) {
93 GLOP_RETURN_ERROR_IF_NULL(entering_col);
95 const DenseRow& update_coefficient = update_row.GetCoefficients();
96 const DenseRow& reduced_costs = reduced_costs_->GetReducedCosts();
97 SCOPED_TIME_STAT(&stats_);
98
99 breakpoints_.clear();
100 breakpoints_.reserve(update_row.GetNonZeroPositions().size());
101 const Fractional threshold = parameters_.ratio_test_zero_threshold();
102 const DenseBitRow& can_decrease = variables_info_.GetCanDecreaseBitRow();
103 const DenseBitRow& can_increase = variables_info_.GetCanIncreaseBitRow();
104 const DenseBitRow& is_boxed = variables_info_.GetNonBasicBoxedVariables();
105
106 // Harris ratio test. See below for more explanation. Here this is used to
107 // prune the first pass by not enqueueing ColWithRatio for columns that have
108 // a ratio greater than the current harris_ratio.
109 const Fractional harris_tolerance =
110 parameters_.harris_tolerance_ratio() *
111 reduced_costs_->GetDualFeasibilityTolerance();
113
114 for (const ColIndex col : update_row.GetNonZeroPositions()) {
115 // We will add ratio * coeff to this column with a ratio positive or zero.
116 // cost_variation makes sure the leaving variable will be dual-feasible
117 // (its update coeff is sign(cost_variation) * 1.0).
118 const Fractional coeff = (cost_variation > 0.0) ? update_coefficient[col]
119 : -update_coefficient[col];
120
121 // In this case, at some point the reduced cost will be positive if not
122 // already, and the column will be dual-infeasible.
123 if (can_decrease.IsSet(col) && coeff > threshold) {
124 if (!is_boxed[col]) {
125 if (-reduced_costs[col] > harris_ratio * coeff) continue;
126 harris_ratio = std::min(
127 harris_ratio, (-reduced_costs[col] + harris_tolerance) / coeff);
128 harris_ratio = std::max(0.0, harris_ratio);
129 }
130 breakpoints_.push_back(ColWithRatio(col, -reduced_costs[col], coeff));
131 continue;
132 }
133
134 // In this case, at some point the reduced cost will be negative if not
135 // already, and the column will be dual-infeasible.
136 if (can_increase.IsSet(col) && coeff < -threshold) {
137 if (!is_boxed[col]) {
138 if (reduced_costs[col] > harris_ratio * -coeff) continue;
139 harris_ratio = std::min(
140 harris_ratio, (reduced_costs[col] + harris_tolerance) / -coeff);
141 harris_ratio = std::max(0.0, harris_ratio);
142 }
143 breakpoints_.push_back(ColWithRatio(col, reduced_costs[col], -coeff));
144 continue;
145 }
146 }
147
148 // Process the breakpoints in priority order as suggested by Maros in
149 // I. Maros, "A generalized dual phase-2 simplex algorithm", European Journal
150 // of Operational Research, 149(1):1-16, 2003.
151 // We use directly make_heap() to avoid a copy of breakpoints, benchmark shows
152 // that it is slightly faster.
153 std::make_heap(breakpoints_.begin(), breakpoints_.end());
154
155 // Harris ratio test. Since we process the breakpoints by increasing ratio, we
156 // do not need a two-pass algorithm as described in the literature. Each time
157 // we process a new breakpoint, we update the harris_ratio of all the
158 // processed breakpoints. For the first new breakpoint with a ratio greater
159 // than the current harris_ratio we know that:
160 // - All the unprocessed breakpoints will have a ratio greater too, so they
161 // will not contribute to the minimum Harris ratio.
162 // - We thus have the actual harris_ratio.
163 // - We have processed all breakpoints with a ratio smaller than it.
165
166 *entering_col = kInvalidCol;
167 bound_flip_candidates->clear();
168 Fractional best_coeff = -1.0;
169 Fractional variation_magnitude = std::abs(cost_variation);
170 equivalent_entering_choices_.clear();
171 while (!breakpoints_.empty()) {
172 const ColWithRatio top = breakpoints_.front();
173 if (top.ratio > harris_ratio) break;
174
175 // If the column is boxed, we can just switch its bounds and
176 // ignore the breakpoint! But we need to see if the entering row still
177 // improve the objective. This is called the bound flipping ratio test in
178 // the literature. See for instance:
179 // http://www.mpi-inf.mpg.de/conferences/adfocs-03/Slides/Bixby_2.pdf
180 //
181 // For each bound flip, |cost_variation| decreases by
182 // |upper_bound - lower_bound| times |coeff|.
183 //
184 // Note that the actual flipping will be done afterwards by
185 // MakeBoxedVariableDualFeasible() in revised_simplex.cc.
186 if (variation_magnitude > threshold) {
187 if (is_boxed[top.col]) {
188 variation_magnitude -=
189 variables_info_.GetBoundDifference(top.col) * top.coeff_magnitude;
190 if (variation_magnitude > threshold) {
191 bound_flip_candidates->push_back(top.col);
192 std::pop_heap(breakpoints_.begin(), breakpoints_.end());
193 breakpoints_.pop_back();
194 continue;
195 }
196 }
197 }
198
199 // TODO(user): We want to maximize both the ratio (objective improvement)
200 // and the coeff_magnitude (stable pivot), so we have to make some
201 // trade-offs. Investigate alternative strategies.
202 if (top.coeff_magnitude >= best_coeff) {
203 // Update harris_ratio. Note that because we process ratio in order, the
204 // harris ratio can only get smaller if the coeff_magnitude is bigger
205 // than the one of the best coefficient.
206 harris_ratio = std::min(
207 harris_ratio, top.ratio + harris_tolerance / top.coeff_magnitude);
208
209 // If the dual infeasibility is too high, the harris_ratio can be
210 // negative. In this case we set it to 0.0, allowing any infeasible
211 // position to enter the basis. This is quite important because its
212 // helps in the choice of a stable pivot.
213 harris_ratio = std::max(harris_ratio, 0.0);
214
215 if (top.coeff_magnitude == best_coeff && top.ratio == *step) {
216 DCHECK_NE(*entering_col, kInvalidCol);
217 equivalent_entering_choices_.push_back(top.col);
218 } else {
219 equivalent_entering_choices_.clear();
220 best_coeff = top.coeff_magnitude;
221 *entering_col = top.col;
222
223 // Note that the step is not directly used, so it is okay to leave it
224 // negative.
225 *step = top.ratio;
226 }
227 }
228
229 // Remove the top breakpoint and maintain the heap structure.
230 // This is the same as doing a pop() on a priority_queue.
231 std::pop_heap(breakpoints_.begin(), breakpoints_.end());
232 breakpoints_.pop_back();
233 }
234
235 // Break the ties randomly.
236 if (!equivalent_entering_choices_.empty()) {
237 equivalent_entering_choices_.push_back(*entering_col);
238 *entering_col =
239 equivalent_entering_choices_[std::uniform_int_distribution<int>(
240 0, equivalent_entering_choices_.size() - 1)(*random_)];
242 stats_.num_perfect_ties.Add(equivalent_entering_choices_.size()));
243 }
244
245 if (*entering_col == kInvalidCol) return Status::OK();
246
247 // If the step is 0.0, we make sure the reduced cost is 0.0 so
248 // UpdateReducedCosts() will not take a step that goes in the wrong way (a few
249 // experiments seems to indicate that this is not a good idea). See comment
250 // at the top of UpdateReducedCosts().
251 //
252 // Note that ShiftCost() actually shifts the cost a bit more in order to do a
253 // non-zero step. This helps on degenerate problems. See the comment of
254 // ShiftCost() for more detail.
255 //
256 // TODO(user): Do not do that if we do not end up using this pivot?
257 if (*step <= 0.0) {
258 // In order to be mathematically consistent, we shift the cost of the
259 // entering column in such a way that its reduced cost is indeed zero. This
260 // is called cost-shifting or perturbation in the literature and it does
261 // really help on degenerate problems. The pertubation will be removed once
262 // the pertubed problem is solved to the optimal.
263 reduced_costs_->ShiftCost(*entering_col);
264 }
265 return Status::OK();
266}
267
269 const UpdateRow& update_row, Fractional cost_variation,
270 ColIndex* entering_col, Fractional* step) {
271 GLOP_RETURN_ERROR_IF_NULL(entering_col);
273 const DenseRow& update_coefficient = update_row.GetCoefficients();
274 const DenseRow& reduced_costs = reduced_costs_->GetReducedCosts();
275 SCOPED_TIME_STAT(&stats_);
276
277 // List of breakpoints where a variable change from feasibility to
278 // infeasibility or the opposite.
279 breakpoints_.clear();
280 breakpoints_.reserve(update_row.GetNonZeroPositions().size());
281
282 // Ratio test.
283 const Fractional threshold = parameters_.ratio_test_zero_threshold();
284 const Fractional dual_feasibility_tolerance =
285 reduced_costs_->GetDualFeasibilityTolerance();
286 const DenseBitRow& can_decrease = variables_info_.GetCanDecreaseBitRow();
287 const DenseBitRow& can_increase = variables_info_.GetCanIncreaseBitRow();
288 const VariableTypeRow& variable_type = variables_info_.GetTypeRow();
289 for (const ColIndex col : update_row.GetNonZeroPositions()) {
290 // Boxed variables shouldn't be in the update position list because they
291 // will be dealt with afterwards by MakeBoxedVariableDualFeasible().
293
294 // Fixed variable shouldn't be in the update position list.
296
297 // Skip if the coeff is too small to be a numerically stable pivot.
298 if (std::abs(update_coefficient[col]) < threshold) continue;
299
300 // We will add ratio * coeff to this column. cost_variation makes sure
301 // the leaving variable will be dual-feasible (its update coeff is
302 // sign(cost_variation) * 1.0).
303 //
304 // TODO(user): This is the same in DualChooseEnteringColumn(), remove
305 // duplication?
306 const Fractional coeff = (cost_variation > 0.0) ? update_coefficient[col]
307 : -update_coefficient[col];
308
309 // Only proceed if there is a transition, note that if reduced_costs[col]
310 // is close to zero, then the variable is supposed to be dual-feasible.
311 if (std::abs(reduced_costs[col]) <= dual_feasibility_tolerance) {
312 // Continue if the variation goes in the dual-feasible direction.
313 if (coeff > 0 && !can_decrease.IsSet(col)) continue;
314 if (coeff < 0 && !can_increase.IsSet(col)) continue;
315
316 // Note that here, a variable which is already dual-infeasible will still
317 // have a positive ratio. This may sounds weird, but the idea is to put
318 // first in the sorted breakpoint list a variable which has a reduced
319 // costs close to zero in order to minimize the magnitude of a step in the
320 // wrong direction.
321 } else {
322 // If the two are of the same sign, there is no transition, skip.
323 if (coeff * reduced_costs[col] > 0) continue;
324 }
325
326 // We are sure there is a transition, add it to the set of breakpoints.
327 breakpoints_.push_back(
328 ColWithRatio(col, std::abs(reduced_costs[col]), std::abs(coeff)));
329 }
330
331 // Process the breakpoints in priority order.
332 std::make_heap(breakpoints_.begin(), breakpoints_.end());
333
334 // Because of our priority queue, it is easy to choose a sub-optimal step to
335 // have a stable pivot. The pivot with the highest magnitude and that reduces
336 // the infeasibility the most is chosen.
337 Fractional pivot_magnitude = 0.0;
338
339 // Select the last breakpoint that still improves the infeasibility and has a
340 // numerically stable pivot.
341 *entering_col = kInvalidCol;
342 *step = -1.0;
343 Fractional improvement = std::abs(cost_variation);
344 while (!breakpoints_.empty()) {
345 const ColWithRatio top = breakpoints_.front();
346
347 // We keep the greatest coeff_magnitude for the same ratio.
348 DCHECK(top.ratio > *step ||
349 (top.ratio == *step && top.coeff_magnitude <= pivot_magnitude));
350 if (top.ratio > *step && top.coeff_magnitude >= pivot_magnitude) {
351 *entering_col = top.col;
352 *step = top.ratio;
353 pivot_magnitude = top.coeff_magnitude;
354 }
355 improvement -= top.coeff_magnitude;
356
357 // If the variable is free, then not only do we loose the infeasibility
358 // improvment, we also render it worse if we keep going in the same
359 // direction.
360 if (can_decrease.IsSet(top.col) && can_increase.IsSet(top.col) &&
361 std::abs(reduced_costs[top.col]) > threshold) {
362 improvement -= top.coeff_magnitude;
363 }
364
365 if (improvement <= 0.0) break;
366 std::pop_heap(breakpoints_.begin(), breakpoints_.end());
367 breakpoints_.pop_back();
368 }
369 return Status::OK();
370}
371
372void EnteringVariable::SetParameters(const GlopParameters& parameters) {
373 parameters_ = parameters;
374}
375
376void EnteringVariable::SetPricingRule(GlopParameters::PricingRule rule) {
377 rule_ = rule;
378}
379
381 SCOPED_TIME_STAT(&stats_);
382 const ColIndex num_cols = variables_info_.GetNumberOfColumns();
383 if (unused_columns_.size() != num_cols) {
384 unused_columns_.ClearAndResize(num_cols);
385 }
386
387 // Invert the set of unused columns, minus the basis.
388 const DenseBitRow& is_basic = variables_info_.GetIsBasicBitRow();
389 for (ColIndex col(0); col < num_cols; ++col) {
390 if (unused_columns_.IsSet(col)) {
391 unused_columns_.Clear(col);
392 } else {
393 if (!is_basic.IsSet(col)) {
394 unused_columns_.Set(col);
395 }
396 }
397 }
398 return &unused_columns_;
399}
400
401template <bool normalize, bool nested_pricing>
402void EnteringVariable::DantzigChooseEnteringColumn(ColIndex* entering_col) {
403 DenseRow dummy;
404 const DenseRow& matrix_column_norms =
405 normalize ? primal_edge_norms_->GetMatrixColumnNorms() : dummy;
406 const DenseRow& reduced_costs = reduced_costs_->GetReducedCosts();
407 SCOPED_TIME_STAT(&stats_);
408
409 Fractional best_price(0.0);
410 *entering_col = kInvalidCol;
411 for (const ColIndex col : reduced_costs_->GetDualInfeasiblePositions()) {
412 if (nested_pricing && !unused_columns_.IsSet(col)) continue;
413 const Fractional unormalized_price = std::abs(reduced_costs[col]);
414 if (normalize) {
415 if (unormalized_price > best_price * matrix_column_norms[col]) {
416 best_price = unormalized_price / matrix_column_norms[col];
417 *entering_col = col;
418 }
419 } else {
420 if (unormalized_price > best_price) {
421 best_price = unormalized_price;
422 *entering_col = col;
423 }
424 }
425 }
426}
427
428// TODO(user): Here we could fill a priority queue with the normalized
429// reduced cost of the top n candidate columns. This makes it possible
430// - To respond right away after each bound flip iteration.
431// - To return the top-n choices if we want to consider multiple candidates in
432// the other parts of the simplex algorithm.
433template <bool use_steepest_edge>
434void EnteringVariable::NormalizedChooseEnteringColumn(ColIndex* entering_col) {
435 const DenseRow& weights = use_steepest_edge
436 ? primal_edge_norms_->GetEdgeSquaredNorms()
437 : primal_edge_norms_->GetDevexWeights();
438 const DenseRow& reduced_costs = reduced_costs_->GetReducedCosts();
439 SCOPED_TIME_STAT(&stats_);
440
441 Fractional best_price(0.0);
442 *entering_col = kInvalidCol;
443 equivalent_entering_choices_.clear();
444 for (const ColIndex col : reduced_costs_->GetDualInfeasiblePositions()) {
445 if (use_steepest_edge) {
446 // Note that here the weights are squared.
447 const Fractional squared_reduced_cost = Square(reduced_costs[col]);
448 if (squared_reduced_cost >= best_price * weights[col]) {
449 if (squared_reduced_cost == best_price * weights[col]) {
450 equivalent_entering_choices_.push_back(col);
451 continue;
452 }
453 equivalent_entering_choices_.clear();
454 best_price = squared_reduced_cost / weights[col];
455 *entering_col = col;
456 }
457 } else {
458 const Fractional positive_reduced_cost = std::abs(reduced_costs[col]);
459 if (positive_reduced_cost >= best_price * weights[col]) {
460 if (positive_reduced_cost == best_price * weights[col]) {
461 equivalent_entering_choices_.push_back(col);
462 continue;
463 }
464 equivalent_entering_choices_.clear();
465 best_price = positive_reduced_cost / weights[col];
466 *entering_col = col;
467 }
468 }
469 }
470 // Break the ties randomly.
471 if (!equivalent_entering_choices_.empty()) {
472 equivalent_entering_choices_.push_back(*entering_col);
473 *entering_col =
474 equivalent_entering_choices_[std::uniform_int_distribution<int>(
475 0, equivalent_entering_choices_.size() - 1)(*random_)];
477 stats_.num_perfect_ties.Add(equivalent_entering_choices_.size()));
478 }
479}
480
481} // namespace glop
482} // namespace operations_research
int64 min
Definition: alldiff_cst.cc:138
int64 max
Definition: alldiff_cst.cc:139
#define DCHECK_NE(val1, val2)
Definition: base/logging.h:886
#define LOG(severity)
Definition: base/logging.h:420
#define DCHECK(condition)
Definition: base/logging.h:884
void ClearAndResize(IndexType size)
Definition: bitset.h:436
IndexType size() const
Definition: bitset.h:419
void Clear(IndexType i)
Definition: bitset.h:453
void Set(IndexType i)
Definition: bitset.h:491
bool IsSet(IndexType i) const
Definition: bitset.h:481
ABSL_MUST_USE_RESULT Status PrimalChooseEnteringColumn(ColIndex *entering_col)
ABSL_MUST_USE_RESULT Status DualChooseEnteringColumn(const UpdateRow &update_row, Fractional cost_variation, std::vector< ColIndex > *bound_flip_candidates, ColIndex *entering_col, Fractional *step)
EnteringVariable(const VariablesInfo &variables_info, random_engine_t *random, ReducedCosts *reduced_costs, PrimalEdgeNorms *primal_edge_norms)
ABSL_MUST_USE_RESULT Status DualPhaseIChooseEnteringColumn(const UpdateRow &update_row, Fractional cost_variation, ColIndex *entering_col, Fractional *step)
void SetPricingRule(GlopParameters::PricingRule rule)
void SetParameters(const GlopParameters &parameters)
const DenseBitRow & GetDualInfeasiblePositions() const
Fractional GetDualFeasibilityTolerance() const
static const Status OK()
Definition: status.h:54
const DenseRow & GetCoefficients() const
Definition: update_row.cc:168
const ColIndexVector & GetNonZeroPositions() const
Definition: update_row.cc:170
const DenseBitRow & GetIsBasicBitRow() const
const DenseBitRow & GetNonBasicBoxedVariables() const
Fractional GetBoundDifference(ColIndex col) const
const DenseBitRow & GetCanIncreaseBitRow() const
const DenseBitRow & GetCanDecreaseBitRow() const
const VariableTypeRow & GetTypeRow() const
const ColIndex GetNumberOfColumns() const
SatParameters parameters
ColIndex col
Definition: markowitz.cc:176
Fractional Square(Fractional f)
const ColIndex kInvalidCol(-1)
StrictITIVector< ColIndex, Fractional > DenseRow
Definition: lp_types.h:299
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
std::mt19937 random_engine_t
Definition: random_engine.h:23
#define IF_STATS_ENABLED(instructions)
Definition: stats.h:437
#define SCOPED_TIME_STAT(stats)
Definition: stats.h:438
#define GLOP_RETURN_ERROR_IF_NULL(arg)
Definition: status.h:85