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