1 // Copyright 2010-2021 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 // Common types and constants used by the Linear Programming solver.
15
16 #ifndef OR_TOOLS_LP_DATA_LP_TYPES_H_
17 #define OR_TOOLS_LP_DATA_LP_TYPES_H_
18
19 #include <cmath>
20 #include <cstdint>
21 #include <limits>
22
23 #include "ortools/base/basictypes.h"
DeferredMessage()24 #include "ortools/base/int_type.h"
25 #include "ortools/base/logging.h"
26 #include "ortools/base/strong_vector.h"
27 #include "ortools/util/bitset.h"
28
29 // We use typedefs as much as possible to later permit the usage of
30 // types such as quad-doubles or rationals.
31
32 namespace operations_research {
33 namespace glop {
34
35 // This type is defined to avoid cast issues during index conversions,
36 // e.g. converting ColIndex into RowIndex.
37 // All types should use 'Index' instead of int32_t.
38 typedef int32_t Index;
39
40 // ColIndex is the type for integers representing column/variable indices.
41 // int32s are enough for handling even the largest problems.
42 DEFINE_INT_TYPE(ColIndex, Index);
43
44 // RowIndex is the type for integers representing row/constraint indices.
45 // int32s are enough for handling even the largest problems.
46 DEFINE_INT_TYPE(RowIndex, Index);
47
48 // Get the ColIndex corresponding to the column # row.
49 inline ColIndex RowToColIndex(RowIndex row) { return ColIndex(row.value()); }
50
51 // Get the RowIndex corresponding to the row # col.
52 inline RowIndex ColToRowIndex(ColIndex col) { return RowIndex(col.value()); }
53
54 // Get the integer index corresponding to the col.
55 inline Index ColToIntIndex(ColIndex col) { return col.value(); }
56
57 // Get the integer index corresponding to the row.
58 inline Index RowToIntIndex(RowIndex row) { return row.value(); }
59
60 // EntryIndex is the type for integers representing entry indices.
61 // An entry in a sparse matrix is a pair (row, value) for a given known column.
62 // See classes SparseColumn and SparseMatrix.
63 #if defined(__ANDROID__)
64 DEFINE_INT_TYPE(EntryIndex, int32_t);
65 #else
66 DEFINE_INT_TYPE(EntryIndex, int64_t);
67 #endif
68
69 static inline double ToDouble(double f) { return f; }
70
71 static inline double ToDouble(long double f) { return static_cast<double>(f); }
72
73 // The type Fractional denotes the type of numbers on which the computations are
74 // performed. This is defined as double here, but it could as well be float,
75 // DoubleDouble, QuadDouble, or infinite-precision rationals.
76 // Floating-point representations are binary fractional numbers, thus the name.
77 // (See http://en.wikipedia.org/wiki/Fraction_(mathematics) .)
78 typedef double Fractional;
79
80 // Range max for type Fractional. DBL_MAX for double for example.
81 const double kRangeMax = std::numeric_limits<double>::max();
82
83 // Infinity for type Fractional.
84 const double kInfinity = std::numeric_limits<double>::infinity();
85
86 // Epsilon for type Fractional, i.e. the smallest e such that 1.0 + e != 1.0 .
87 const double kEpsilon = std::numeric_limits<double>::epsilon();
88
89 // Returns true if the given value is finite, that means for a double:
90 // not a NaN and not +/- infinity.
91 inline bool IsFinite(Fractional value) {
92 return value > -kInfinity && value < kInfinity;
93 }
94
95 // Constants to represent invalid row or column index.
96 // It is important that their values be the same because during transposition,
97 // one needs to be converted into the other.
98 const RowIndex kInvalidRow(-1);
99 const ColIndex kInvalidCol(-1);
100
101 // Different statuses for a given problem.
102 enum class ProblemStatus : int8_t {
103 // The problem has been solved to optimality. Both the primal and dual have
104 // a feasible solution.
105 OPTIMAL,
106
107 // The problem has been proven primal-infeasible. Note that the problem is not
108 // necessarily DUAL_UNBOUNDED (See Chvatal p.60). The solver does not have a
109 // dual unbounded ray in this case.
110 PRIMAL_INFEASIBLE,
111
112 // The problem has been proven dual-infeasible. Note that the problem is not
113 // necessarily PRIMAL_UNBOUNDED (See Chvatal p.60). The solver does
114 // note have a primal unbounded ray in this case,
115 DUAL_INFEASIBLE,
116
117 // The problem is either INFEASIBLE or UNBOUNDED (this applies to both the
118 // primal and dual algorithms). This status is only returned by the presolve
119 // step and means that a primal or dual unbounded ray was found during
120 // presolve. Note that because some presolve techniques assume that a feasible
121 // solution exists to simplify the problem further, it is difficult to
122 // distinguish between infeasibility and unboundedness.
123 //
124 // If a client needs to distinguish, it is possible to run the primal
125 // algorithm on the same problem with a 0 objective function to know if the
126 // problem was PRIMAL_INFEASIBLE.
127 INFEASIBLE_OR_UNBOUNDED,
128
129 // The problem has been proven feasible and unbounded. That means that the
130 // problem is DUAL_INFEASIBLE and that the solver has a primal unbounded ray.
131 PRIMAL_UNBOUNDED,
132
133 // The problem has been proven dual-feasible and dual-unbounded. That means
134 // the problem is PRIMAL_INFEASIBLE and that the solver has a dual unbounded
135 // ray to prove it.
136 DUAL_UNBOUNDED,
137
138 // All the statuses below correspond to a case where the solver was
139 // interrupted. This can happen because of a timeout, an iteration limit or an
140 // error.
141
142 // The solver didn't had a chance to prove anything.
143 INIT,
144
145 // The problem has been proven primal-feasible but may still be
146 // PRIMAL_UNBOUNDED.
147 PRIMAL_FEASIBLE,
148
149 // The problem has been proven dual-feasible, but may still be DUAL_UNBOUNDED.
150 // That means that if the primal is feasible, then it has a finite optimal
151 // solution.
152 DUAL_FEASIBLE,
153
154 // An error occurred during the solving process.
155 ABNORMAL,
156
157 // The input problem was invalid (see LinearProgram.IsValid()).
158 INVALID_PROBLEM,
159
160 // The problem was solved to a feasible status, but the solution checker found
161 // the primal and/or dual infeasibilities too important for the specified
162 // parameters.
163 IMPRECISE,
164 };
165
166 // Returns the string representation of the ProblemStatus enum.
167 std::string GetProblemStatusString(ProblemStatus problem_status);
168
169 inline std::ostream& operator<<(std::ostream& os, ProblemStatus status) {
170 os << GetProblemStatusString(status);
171 return os;
172 }
173
174 // Different types of variables.
175 enum class VariableType : int8_t {
176 UNCONSTRAINED,
177 LOWER_BOUNDED,
178 UPPER_BOUNDED,
179 UPPER_AND_LOWER_BOUNDED,
180 FIXED_VARIABLE
181 };
182
183 // Returns the string representation of the VariableType enum.
184 std::string GetVariableTypeString(VariableType variable_type);
185
186 inline std::ostream& operator<<(std::ostream& os, VariableType type) {
187 os << GetVariableTypeString(type);
188 return os;
189 }
190
191 // Different variables statuses.
192 // If a solution is OPTIMAL or FEASIBLE, then all the properties described here
193 // should be satisfied. These properties should also be true during the
194 // execution of the revised simplex algorithm, except that because of
195 // bound-shifting, the variable may not be at their exact bounds until the
196 // shifts are removed.
197 enum class VariableStatus : int8_t {
198 // The basic status is special and takes precedence over all the other
199 // statuses. It means that the variable is part of the basis.
200 BASIC,
201 // Only possible status of a FIXED_VARIABLE not in the basis. The variable
202 // value should be exactly equal to its bounds (which are the same).
203 FIXED_VALUE,
204 // Only possible statuses of a non-basic variable which is not UNCONSTRAINED
205 // or FIXED. The variable value should be at its exact specified bound (which
206 // must be finite).
207 AT_LOWER_BOUND,
208 AT_UPPER_BOUND,
209 // Only possible status of an UNCONSTRAINED non-basic variable.
210 // Its value should be zero.
211 //
212 // Note that during crossover, this status is relaxed, and any variable that
213 // can currently move in both directions can be marked as free.
214 FREE,
215 };
216
217 // Returns the string representation of the VariableStatus enum.
218 std::string GetVariableStatusString(VariableStatus status);
219
220 inline std::ostream& operator<<(std::ostream& os, VariableStatus status) {
221 os << GetVariableStatusString(status);
222 return os;
223 }
224
225 // Different constraints statuses.
226 // The meaning is the same for the constraint activity relative to its bounds as
227 // it is for a variable value relative to its bounds. Actually, this is the
228 // VariableStatus of the slack variable associated to a constraint modulo a
229 // change of sign. The difference is that because of precision error, a
230 // constraint activity cannot exactly be equal to one of its bounds or to zero.
231 enum class ConstraintStatus : int8_t {
232 BASIC,
233 FIXED_VALUE,
234 AT_LOWER_BOUND,
235 AT_UPPER_BOUND,
236 FREE,
237 };
238
239 // Returns the string representation of the ConstraintStatus enum.
240 std::string GetConstraintStatusString(ConstraintStatus status);
241
242 inline std::ostream& operator<<(std::ostream& os, ConstraintStatus status) {
243 os << GetConstraintStatusString(status);
244 return os;
245 }
246
247 // Returns the ConstraintStatus corresponding to a given VariableStatus.
248 ConstraintStatus VariableToConstraintStatus(VariableStatus status);
249
250 // Wrapper around an ITIVector to allow (and enforce) creation/resize/assign
251 // to use the index type for the size.
252 //
253 // TODO(user): This should probably move into ITIVector, but note that this
254 // version is more strict and does not allow any other size types.
255 template <typename IntType, typename T>
256 class StrictITIVector : public absl::StrongVector<IntType, T> {
257 public:
258 typedef IntType IndexType;
259 typedef absl::StrongVector<IntType, T> ParentType;
260 // This allows for brace initialization, which is really useful in tests.
261 // It is not 'explicit' by design, so one can do vector = {...};
262 #if !defined(__ANDROID__) && (!defined(_MSC_VER) || (_MSC_VER >= 1800))
263 StrictITIVector(std::initializer_list<T> init_list) // NOLINT
264 : ParentType(init_list.begin(), init_list.end()) {}
265 #endif
266 StrictITIVector() : ParentType() {}
267 explicit StrictITIVector(IntType size) : ParentType(size.value()) {}
268 StrictITIVector(IntType size, const T& v) : ParentType(size.value(), v) {}
269 template <typename InputIteratorType>
270 StrictITIVector(InputIteratorType first, InputIteratorType last)
271 : ParentType(first, last) {}
272
273 void resize(IntType size) { ParentType::resize(size.value()); }
274 void resize(IntType size, const T& v) { ParentType::resize(size.value(), v); }
275
276 void reserve(IntType size) { ParentType::reserve(size.value()); }
277
278 void assign(IntType size, const T& v) { ParentType::assign(size.value(), v); }
279
280 IntType size() const { return IntType(ParentType::size()); }
281
282 IntType capacity() const { return IntType(ParentType::capacity()); }
283
284 // Since calls to resize() must use a default value, we introduce a new
285 // function for convenience to reduce the size of a vector.
286 void resize_down(IntType size) {
287 DCHECK_GE(size, IntType(0));
288 DCHECK_LE(size, IntType(ParentType::size()));
289 ParentType::resize(size.value());
290 }
291
292 // This function can be up to 4 times faster than calling assign(size, 0).
293 // Note that it only works with StrictITIVector of basic types.
294 void AssignToZero(IntType size) {
295 resize(size, 0);
296 memset(ParentType::data(), 0, size.value() * sizeof(T));
297 }
298 };
299
300 // Row-vector types. Row-vector types are indexed by a column index.
301
302 // Row of fractional values.
303 typedef StrictITIVector<ColIndex, Fractional> DenseRow;
304
305 // Row of booleans.
306 typedef StrictITIVector<ColIndex, bool> DenseBooleanRow;
307
308 // Row of column indices. Used to represent mappings between columns.
309 typedef StrictITIVector<ColIndex, ColIndex> ColMapping;
310
311 // Vector of row or column indices. Useful to list the non-zero positions.
312 typedef std::vector<ColIndex> ColIndexVector;
313 typedef std::vector<RowIndex> RowIndexVector;
314
315 // Row of row indices.
316 // Useful for knowing which row corresponds to a particular column in the basis,
317 // or for storing the number of rows for a given column.
318 typedef StrictITIVector<ColIndex, RowIndex> ColToRowMapping;
319
320 // Row of variable types.
321 typedef StrictITIVector<ColIndex, VariableType> VariableTypeRow;
322
323 // Row of variable statuses.
324 typedef StrictITIVector<ColIndex, VariableStatus> VariableStatusRow;
325
326 // Row of bits.
327 typedef Bitset64<ColIndex> DenseBitRow;
328
329 // Column-vector types. Column-vector types are indexed by a row index.
330
331 // Column of fractional values.
332 typedef StrictITIVector<RowIndex, Fractional> DenseColumn;
333
334 // Column of booleans.
335 typedef StrictITIVector<RowIndex, bool> DenseBooleanColumn;
336
337 // Column of bits.
338 typedef Bitset64<RowIndex> DenseBitColumn;
339
340 // Column of row indices. Used to represent mappings between rows.
341 typedef StrictITIVector<RowIndex, RowIndex> RowMapping;
342
343 // Column of column indices.
344 // Used to represent which column corresponds to a particular row in the basis,
345 // or for storing the number of columns for a given row.
346 typedef StrictITIVector<RowIndex, ColIndex> RowToColMapping;
347
348 // Column of constraints (slack variables) statuses.
349 typedef StrictITIVector<RowIndex, ConstraintStatus> ConstraintStatusColumn;
350
351 // --------------------------------------------------------
352 // VectorIterator
353 // --------------------------------------------------------
354
355 // An iterator over the elements of a sparse data structure that stores the
356 // elements in arrays for indices and coefficients. The iterator is
357 // built as a wrapper over a sparse vector entry class; the concrete entry class
358 // is provided through the template argument EntryType.
359 template <typename EntryType>
360 class VectorIterator : EntryType {
361 public:
362 using Index = typename EntryType::Index;
363 using Entry = EntryType;
364
365 VectorIterator(const Index* indices, const Fractional* coefficients,
366 EntryIndex i)
367 : EntryType(indices, coefficients, i) {}
368
369 void operator++() { ++this->i_; }
370 bool operator!=(const VectorIterator& other) const {
371 // This operator is intended for use in natural range iteration ONLY.
372 // Therefore, we prefer to use '<' so that a buggy range iteration which
373 // start point is *after* its end point stops immediately, instead of
374 // iterating 2^(number of bits of EntryIndex) times.
375 return this->i_ < other.i_;
376 }
377 const Entry& operator*() const { return *this; }
378 };
379
380 // This is used during the deterministic time computation to convert a given
381 // number of floating-point operations to something in the same order of
382 // magnitude as a second (on a 2014 desktop).
383 static inline double DeterministicTimeForFpOperations(int64_t n) {
384 const double kConversionFactor = 2e-9;
385 return kConversionFactor * static_cast<double>(n);
386 }
387
388 } // namespace glop
389 } // namespace operations_research
390
391 #endif // OR_TOOLS_LP_DATA_LP_TYPES_H_
392