OR-Tools  8.2
reduced_costs.h
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
14#ifndef OR_TOOLS_GLOP_REDUCED_COSTS_H_
15#define OR_TOOLS_GLOP_REDUCED_COSTS_H_
16
18#include "ortools/glop/parameters.pb.h"
20#include "ortools/glop/status.h"
27#include "ortools/util/stats.h"
28
29namespace operations_research {
30namespace glop {
31
32// Maintains the reduced costs of the non-basic variables and some related
33// quantities.
34//
35// Terminology:
36// - To each non-basic column 'a' of A, we can associate an "edge" in the
37// kernel of A equal to 1.0 on the index of 'a' and '-B^{-1}.a' on the basic
38// variables.
39// - 'B^{-1}.a' is called the "right inverse" of 'a'.
40// - The reduced cost of a column is equal to the scalar product of this
41// column's edge with the cost vector (objective_), and corresponds to the
42// variation in the objective function when we add this edge to the current
43// solution.
44// - The dual values are the "left inverse" of the basic objective by B.
45// That is 'basic_objective_.B^{-1}'
46// - The reduced cost of a column is also equal to the scalar product of this
47// column with the vector of the dual values.
49 public:
50 // Takes references to the linear program data we need.
51 ReducedCosts(const CompactSparseMatrix& matrix_, const DenseRow& objective,
52 const RowToColMapping& basis,
53 const VariablesInfo& variables_info,
54 const BasisFactorization& basis_factorization,
55 random_engine_t* random);
56
57 // If this is true, then the caller must re-factorize the basis before the
58 // next call to GetReducedCosts().
59 bool NeedsBasisRefactorization() const;
60
61 // Checks the precision of the entering variable choice now that the direction
62 // is computed, and return true if we can continue with this entering column,
63 // or false if this column is actually not good and ChooseEnteringColumn()
64 // need to be called again.
65 bool TestEnteringReducedCostPrecision(ColIndex entering_col,
66 const ScatteredColumn& direction,
67 Fractional* reduced_cost);
68
69 // Computes the current dual residual and infeasibility. Note that these
70 // functions are not really fast (many scalar products will be computed) and
71 // shouldn't be called at each iteration. They will return 0.0 if the reduced
72 // costs need to be recomputed first and fail in debug mode.
76
77 // Updates any internal data BEFORE the given simplex pivot is applied to B.
78 // Note that no updates are needed in case of a bound flip.
79 // The arguments are in order:
80 // - The index of the entering non-basic column of A.
81 // - The index in B of the leaving basic variable.
82 // - The 'direction', i.e. the right inverse of the entering column.
83 void UpdateBeforeBasisPivot(ColIndex entering_col, RowIndex leaving_row,
84 const ScatteredColumn& direction,
85 UpdateRow* update_row);
86
87 // Once a pivot has been done, this need to be called on the column that just
88 // left the basis before the next ChooseEnteringColumn(). On a bound flip,
89 // this must also be called with the variable that flipped.
90 //
91 // In both cases, the variable should be dual-feasible by construction. This
92 // sets is_dual_infeasible_[col] to false and checks in debug mode that the
93 // variable is indeed not an entering candidate.
95
96 // Sets the cost of the given non-basic variable to zero and updates its
97 // reduced cost. Note that changing the cost of a non-basic variable only
98 // impacts its reduced cost and not the one of any other variables.
99 // The current_cost pointer must be equal to the address of objective[col]
100 // where objective is the DenseRow passed at construction.
101 void SetNonBasicVariableCostToZero(ColIndex col, Fractional* current_cost);
102
103 // Sets the pricing parameters. This does not change the pricing rule.
104 void SetParameters(const GlopParameters& parameters);
105
106 // Returns true if the current reduced costs are computed with maximum
107 // precision.
108 bool AreReducedCostsPrecise() { return are_reduced_costs_precise_; }
109
110 // Returns true if the current reduced costs where just recomputed or will be
111 // on the next call to GetReducedCosts().
113 return recompute_reduced_costs_ || are_reduced_costs_recomputed_;
114 }
115
116 // Makes sure the next time the reduced cost are needed, they will be
117 // recomputed with maximum precision (i.e. from scratch with a basis
118 // refactorization first).
120
121 // Randomly perturb the costs. Both Koberstein and Huangfu recommend doing
122 // that before the dual simplex starts in their Phd thesis.
123 //
124 // The perturbation follows what is explained in Huangfu Q (2013) "High
125 // performance simplex solver", Ph.D, dissertation, University of Edinburgh,
126 // section 3.2.3, page 58.
127 void PerturbCosts();
128
129 // Shifts the cost of the given non-basic column such that its current reduced
130 // cost becomes 0.0. Actually, this shifts the cost a bit more, so that the
131 // reduced cost becomes slightly of the other sign.
132 //
133 // This is explained in Koberstein's thesis (section 6.2.2.3) and helps on
134 // degenerate problems. As of july 2013, this allowed to pass dano3mip and
135 // dbic1 without cycling forever. Note that contrary to what is explained in
136 // the thesis, we do not shift any other variable costs. If any becomes
137 // infeasible, it will be selected and shifted in subsequent iterations.
138 void ShiftCost(ColIndex col);
139
140 // Removes any cost shift and cost perturbation. This also lazily forces a
141 // recomputation of all the derived quantities. This effectively resets the
142 // class to its initial state.
144
145 // Invalidates all internal structure that depends on the objective function.
147
148 // Sets whether or not the bitset of the dual infeasible positions is updated.
149 void MaintainDualInfeasiblePositions(bool maintain);
150
151 // Invalidates the data that depends on the order of the column in basis_.
153
154 // Returns the current reduced costs. If AreReducedCostsPrecise() is true,
155 // then for basic columns, this gives the error between 'c_B' and 'y.B' and
156 // for non-basic columns, this is the classic reduced cost. If it is false,
157 // then this is defined only for the columns in
158 // variables_info_.GetIsRelevantBitRow().
159 const DenseRow& GetReducedCosts();
160
161 // Returns the non-basic columns that are dual-infeasible. These are also
162 // the primal simplex possible entering columns.
164 // TODO(user): recompute if needed?
165 DCHECK(are_dual_infeasible_positions_maintained_);
166 return is_dual_infeasible_;
167 }
168
169 // Returns the dual values associated to the current basis.
170 const DenseColumn& GetDualValues();
171
172 // Stats related functions.
173 std::string StatString() const { return stats_.StatString(); }
174
175 // Returns the current dual feasibility tolerance.
177 return dual_feasibility_tolerance_;
178 }
179
180 // Does basic checking of an entering candidate. To be used in DCHECK().
181 bool IsValidPrimalEnteringCandidate(ColIndex col) const;
182
183 // Visible for testing.
184 const DenseRow& GetCostPerturbations() const { return cost_perturbations_; }
185
186 private:
187 // Statistics about this class.
188 struct Stats : public StatsGroup {
189 Stats()
190 : StatsGroup("ReducedCosts"),
191 basic_objective_left_inverse_density(
192 "basic_objective_left_inverse_density", this),
193 reduced_costs_accuracy("reduced_costs_accuracy", this),
194 cost_shift("cost_shift", this) {}
195 RatioDistribution basic_objective_left_inverse_density;
196 DoubleDistribution reduced_costs_accuracy;
197 DoubleDistribution cost_shift;
198 };
199
200 // Small utility function to be called before reduced_costs_ and/or
201 // is_dual_infeasible_ are accessed.
202 void RecomputeReducedCostsAndPrimalEnteringCandidatesIfNeeded();
203
204 // All these Compute() functions fill the corresponding DenseRow using
205 // the current problem data.
206 void ComputeBasicObjective();
207 void ComputeReducedCosts();
208 void ComputeBasicObjectiveLeftInverse();
209
210 // Updates reduced_costs_ according to the given pivot. This adds a multiple
211 // of the vector equal to 1.0 on the leaving column and given by
212 // ComputeUpdateRow() on the non-basic columns. The multiple is such that the
213 // new leaving reduced cost is zero. If update_column_status is false, then
214 // is_dual_infeasible_ will not be updated.
215 void UpdateReducedCosts(ColIndex entering_col, ColIndex leaving_col,
216 RowIndex leaving_row, Fractional pivot,
217 UpdateRow* update_row);
218
219 // Recomputes from scratch the is_dual_infeasible_ bit row. Note that an
220 // entering candidate is by definition a dual-infeasible variable.
221 void ResetDualInfeasibilityBitSet();
222
223 // Recomputes is_dual_infeasible_ but only for the given column indices.
224 template <typename ColumnsToUpdate>
225 void UpdateEnteringCandidates(const ColumnsToUpdate& cols);
226
227 // Updates basic_objective_ according to the given pivot.
228 void UpdateBasicObjective(ColIndex entering_col, RowIndex leaving_row);
229
230 // Problem data that should be updated from outside.
231 const CompactSparseMatrix& matrix_;
232 const DenseRow& objective_;
233 const RowToColMapping& basis_;
234 const VariablesInfo& variables_info_;
235 const BasisFactorization& basis_factorization_;
236 random_engine_t* random_;
237
238 // Internal data.
239 GlopParameters parameters_;
240 mutable Stats stats_;
241
242 // Booleans to control what happens on the next ChooseEnteringColumn() call.
243 bool must_refactorize_basis_;
244 bool recompute_basic_objective_left_inverse_;
245 bool recompute_basic_objective_;
246 bool recompute_reduced_costs_;
247
248 // Indicates if we have computed the reduced costs with a good precision.
249 bool are_reduced_costs_precise_;
250 bool are_reduced_costs_recomputed_;
251
252 // Values of the objective on the columns of the basis. The order is given by
253 // the basis_ mapping. It is usually denoted as 'c_B' in the literature .
254 DenseRow basic_objective_;
255
256 // Perturbations to the objective function. This may be introduced to
257 // counter degenerecency. It will be removed at the end of the algorithm.
258 DenseRow cost_perturbations_;
259
260 // Reduced costs of the relevant columns of A.
261 DenseRow reduced_costs_;
262
263 // Left inverse by B of the basic_objective_. This is known as 'y' or 'pi' in
264 // the literature. Its scalar product with a column 'a' of A gives the value
265 // of the scalar product of the basic objective with the right inverse of 'a'.
266 //
267 // TODO(user): using the unit_row_left_inverse_, we can update the
268 // basic_objective_left_inverse_ at each iteration, this is not needed for the
269 // algorithm, but may gives us a good idea of the current precision of our
270 // estimates. It is also faster to compute the unit_row_left_inverse_ because
271 // of sparsity.
272 ScatteredRow basic_objective_left_inverse_;
273
274 // This is usually parameters_.dual_feasibility_tolerance() except when the
275 // dual residual error |y.B - c_B| is higher than it and we have to increase
276 // the tolerance.
277 Fractional dual_feasibility_tolerance_;
278
279 // True for the columns that can enter the basis.
280 // TODO(user): Investigate using a list of ColIndex instead (but we need
281 // dynamic update and may lose the ordering of the indices).
282 DenseBitRow is_dual_infeasible_;
283
284 // Indicates if the dual-infeasible positions are maintained or not.
285 bool are_dual_infeasible_positions_maintained_;
286
287 DISALLOW_COPY_AND_ASSIGN(ReducedCosts);
288};
289
290} // namespace glop
291} // namespace operations_research
292
293#endif // OR_TOOLS_GLOP_REDUCED_COSTS_H_
#define DCHECK(condition)
Definition: base/logging.h:884
ReducedCosts(const CompactSparseMatrix &matrix_, const DenseRow &objective, const RowToColMapping &basis, const VariablesInfo &variables_info, const BasisFactorization &basis_factorization, random_engine_t *random)
bool TestEnteringReducedCostPrecision(ColIndex entering_col, const ScatteredColumn &direction, Fractional *reduced_cost)
const DenseBitRow & GetDualInfeasiblePositions() const
bool IsValidPrimalEnteringCandidate(ColIndex col) const
void SetNonBasicVariableCostToZero(ColIndex col, Fractional *current_cost)
void SetAndDebugCheckThatColumnIsDualFeasible(ColIndex col)
void UpdateBeforeBasisPivot(ColIndex entering_col, RowIndex leaving_row, const ScatteredColumn &direction, UpdateRow *update_row)
Fractional GetDualFeasibilityTolerance() const
Fractional ComputeMaximumDualInfeasibility() const
void MaintainDualInfeasiblePositions(bool maintain)
Fractional ComputeSumOfDualInfeasibilities() const
const DenseRow & GetCostPerturbations() const
void SetParameters(const GlopParameters &parameters)
SatParameters parameters
ColIndex col
Definition: markowitz.cc:176
StrictITIVector< ColIndex, Fractional > DenseRow
Definition: lp_types.h:299
Bitset64< ColIndex > DenseBitRow
Definition: lp_types.h:323
StrictITIVector< RowIndex, ColIndex > RowToColMapping
Definition: lp_types.h:342
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