OR-Tools  8.2
lp_data/lp_utils.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// Basic utility functions on Fractional or row/column of Fractional.
15
16#ifndef OR_TOOLS_LP_DATA_LP_UTILS_H_
17#define OR_TOOLS_LP_DATA_LP_UTILS_H_
18
23
24namespace operations_research {
25namespace glop {
26
27// TODO(user): For some Fractional types, it may not gain much (or even nothing
28// if we are in infinite precision) to use this sum. A solution is to templatize
29// this class and specialize it to a normal sum for the Fractional type we want
30// so in this case the PreciseXXX() functions below will become equivalent to
31// their normal version.
33
34// Returns the square of a Fractional.
35// Useful to shorten the code when f is an expression or a long name.
36inline Fractional Square(Fractional f) { return f * f; }
37
38// Returns distance from a given fractional number to the closest integer. It
39// means that the result is always contained in range of [0.0, 0.5].
41 return std::abs(f - std::round(f));
42}
43
44// Returns the scalar product between u and v.
45// The precise versions use KahanSum and are about two times slower.
46template <class DenseRowOrColumn1, class DenseRowOrColumn2>
47Fractional ScalarProduct(const DenseRowOrColumn1& u,
48 const DenseRowOrColumn2& v) {
49 DCHECK_EQ(u.size().value(), v.size().value());
50 Fractional sum(0.0);
51 typename DenseRowOrColumn1::IndexType i(0);
52 typename DenseRowOrColumn2::IndexType j(0);
53 const size_t num_blocks = u.size().value() / 4;
54 for (size_t block = 0; block < num_blocks; ++block) {
55 // Computing the sum of 4 elements at once may allow the compiler to
56 // generate more efficient code, e.g. using SIMD and checking the loop
57 // condition much less frequently.
58 //
59 // This produces different results from the case where each multiplication
60 // is added to sum separately. An extreme example of this can be derived
61 // using the fact that 1e11 + 2e-6 == 1e11, but 1e11 + 8e-6 > 1e11.
62 //
63 // While the results are different, they aren't necessarily better or worse.
64 // Typically, sum will be of larger magnitude than any individual
65 // multiplication, so one might expect, in practice, this method to yield
66 // more accurate results. However, if accuracy is vital, use the precise
67 // version.
68 sum += (u[i++] * v[j++]) + (u[i++] * v[j++]) + (u[i++] * v[j++]) +
69 (u[i++] * v[j++]);
70 }
71 while (i < u.size()) {
72 sum += u[i++] * v[j++];
73 }
74 return sum;
75}
76
77// Note: This version is heavily used in the pricing.
78// TODO(user): Optimize this more (SSE or unroll with two sums). Another
79// option is to skip the u[col] that are 0.0 rather than fetching the coeff
80// and doing a Fractional multiplication.
81template <class DenseRowOrColumn>
82Fractional ScalarProduct(const DenseRowOrColumn& u, const SparseColumn& v) {
83 Fractional sum(0.0);
84 for (const SparseColumn::Entry e : v) {
85 sum += u[typename DenseRowOrColumn::IndexType(e.row().value())] *
86 e.coefficient();
87 }
88 return sum;
89}
90
91template <class DenseRowOrColumn, class DenseRowOrColumn2>
92Fractional PreciseScalarProduct(const DenseRowOrColumn& u,
93 const DenseRowOrColumn2& v) {
94 DCHECK_EQ(u.size().value(), v.size().value());
95 KahanSum sum;
96 for (typename DenseRowOrColumn::IndexType i(0); i < u.size(); ++i) {
97 sum.Add(u[i] * v[typename DenseRowOrColumn2::IndexType(i.value())]);
98 }
99 return sum.Value();
100}
101
102template <class DenseRowOrColumn>
103Fractional PreciseScalarProduct(const DenseRowOrColumn& u,
104 const SparseColumn& v) {
105 KahanSum sum;
106 for (const SparseColumn::Entry e : v) {
107 sum.Add(u[typename DenseRowOrColumn::IndexType(e.row().value())] *
108 e.coefficient());
109 }
110 return sum.Value();
111}
112
113template <class DenseRowOrColumn>
114Fractional PreciseScalarProduct(const DenseRowOrColumn& u,
115 const ScatteredColumn& v) {
116 DCHECK_EQ(u.size().value(), v.values.size().value());
117 if (v.ShouldUseDenseIteration()) {
118 return PreciseScalarProduct(u, v.values);
119 }
120 KahanSum sum;
121 for (const auto e : v) {
122 sum.Add(u[typename DenseRowOrColumn::IndexType(e.row().value())] *
123 e.coefficient());
124 }
125 return sum.Value();
126}
127
128// Computes a scalar product for entries with index not greater than max_index.
129template <class DenseRowOrColumn>
130Fractional PartialScalarProduct(const DenseRowOrColumn& u,
131 const SparseColumn& v, int max_index) {
132 Fractional sum(0.0);
133 for (const SparseColumn::Entry e : v) {
134 if (e.row().value() >= max_index) {
135 return sum;
136 }
137 sum += u[typename DenseRowOrColumn::IndexType(e.row().value())] *
138 e.coefficient();
139 }
140 return sum;
141}
142
143// Returns the norm^2 (sum of the square of the entries) of the given column.
144// The precise version uses KahanSum and are about two times slower.
145Fractional SquaredNorm(const SparseColumn& v);
146Fractional SquaredNorm(const DenseColumn& column);
147Fractional SquaredNorm(const ColumnView& v);
148Fractional PreciseSquaredNorm(const SparseColumn& v);
150Fractional PreciseSquaredNorm(const ScatteredColumn& v);
151
152// Returns the maximum of the |coefficients| of 'v'.
154Fractional InfinityNorm(const SparseColumn& v);
155Fractional InfinityNorm(const ColumnView& v);
156
157// Returns the fraction of non-zero entries of the given row.
158//
159// TODO(user): Take a Scattered row/col instead. This is only used to report
160// stats, but we should still have a sparse version to do it faster.
161double Density(const DenseRow& row);
162
163// Sets to 0.0 all entries of the given row whose fabs() is lower than the given
164// threshold.
166void RemoveNearZeroEntries(Fractional threshold, DenseColumn* column);
167
168// Transposition functions implemented below with a cast so it should actually
169// have no complexity cost.
170const DenseRow& Transpose(const DenseColumn& col);
171const DenseColumn& Transpose(const DenseRow& row);
172
173// Returns the maximum of the |coefficients| of the given column restricted
174// to the rows_to_consider. Also returns the first RowIndex 'row' that attains
175// this maximum. If the maximum is 0.0, then row_index is left untouched.
176Fractional RestrictedInfinityNorm(const ColumnView& column,
177 const DenseBooleanColumn& rows_to_consider,
178 RowIndex* row_index);
179
180// Sets to false the entry b[row] if column[row] is non null.
181// Note that if 'b' was true only on the non-zero position of column, this can
182// be used as a fast way to clear 'b'.
183void SetSupportToFalse(const ColumnView& column, DenseBooleanColumn* b);
184
185// Returns true iff for all 'row' we have '|column[row]| <= radius[row]'.
186bool IsDominated(const ColumnView& column, const DenseColumn& radius);
187
188// This cast based implementation should be safe, as long as DenseRow and
189// DenseColumn are implemented by the same underlying type.
190// We still do some DCHECK to be sure it works as expected in addition to the
191// unit tests.
192inline const DenseRow& Transpose(const DenseColumn& col) {
193 const DenseRow& row = reinterpret_cast<const DenseRow&>(col);
194 DCHECK_EQ(col.size(), ColToRowIndex(row.size()));
195 DCHECK(col.empty() || (&(col[RowIndex(0)]) == &(row[ColIndex(0)])));
196 return row;
197}
198
199// Similar comment as the other Transpose() implementation above.
200inline const DenseColumn& Transpose(const DenseRow& row) {
201 const DenseColumn& col = reinterpret_cast<const DenseColumn&>(row);
202 DCHECK_EQ(col.size(), ColToRowIndex(row.size()));
203 DCHECK(col.empty() || (&(col[RowIndex(0)]) == &(row[ColIndex(0)])));
204 return col;
205}
206
207// Computes the positions of the non-zeros of a dense vector.
208template <typename IndexType>
210 std::vector<IndexType>* non_zeros) {
211 non_zeros->clear();
212 const IndexType end = input.size();
213 for (IndexType index(0); index < end; ++index) {
214 if (input[index] != 0.0) {
215 non_zeros->push_back(index);
216 }
217 }
218}
219
220// Returns true if the given Fractional container is all zeros.
221template <typename Container>
222inline bool IsAllZero(const Container& input) {
223 for (Fractional value : input) {
224 if (value != 0.0) return false;
225 }
226 return true;
227}
228
229// Returns true if the given vector of bool is all false.
230template <typename BoolVector>
231bool IsAllFalse(const BoolVector& v) {
232 return std::all_of(v.begin(), v.end(), [](bool value) { return !value; });
233}
234
235// Permutes the given dense vector. It uses for this an all zero scratchpad.
236template <typename IndexType, typename PermutationIndexType>
238 const Permutation<PermutationIndexType>& permutation,
241 DCHECK(IsAllZero(*zero_scratchpad));
242 const IndexType size = input_output->size();
243 zero_scratchpad->swap(*input_output);
244 input_output->resize(size, 0.0);
245 for (IndexType index(0); index < size; ++index) {
246 const Fractional value = (*zero_scratchpad)[index];
247 if (value != 0.0) {
248 const IndexType permuted_index(
249 permutation[PermutationIndexType(index.value())].value());
250 (*input_output)[permuted_index] = value;
251 }
252 }
253 zero_scratchpad->assign(size, 0.0);
254}
255
256// Same as PermuteAndComputeNonZeros() except that we assume that the given
257// non-zeros are the initial non-zeros positions of output.
258template <typename IndexType>
260 const Permutation<IndexType>& permutation,
263 std::vector<IndexType>* non_zeros) {
264 DCHECK(IsAllZero(*zero_scratchpad));
265 zero_scratchpad->swap(*output);
266 output->resize(zero_scratchpad->size(), 0.0);
267 for (IndexType& index_ref : *non_zeros) {
268 const Fractional value = (*zero_scratchpad)[index_ref];
269 (*zero_scratchpad)[index_ref] = 0.0;
270 const IndexType permuted_index(permutation[index_ref]);
271 (*output)[permuted_index] = value;
272 index_ref = permuted_index;
273 }
274}
275
276// Sets a dense vector for which the non zeros are known to be non_zeros.
277template <typename IndexType, typename ScatteredRowOrCol>
278inline void ClearAndResizeVectorWithNonZeros(IndexType size,
279 ScatteredRowOrCol* v) {
280 // Only use the sparse version if there is less than 5% non-zeros positions
281 // compared to the wanted size. Note that in most cases the vector will
282 // already be of the correct size.
283 const double kSparseThreshold = 0.05;
284 if (!v->non_zeros.empty() &&
285 v->non_zeros.size() < kSparseThreshold * size.value()) {
286 for (const IndexType index : v->non_zeros) {
287 DCHECK_LT(index, v->values.size());
288 (*v)[index] = 0.0;
289 }
290 v->values.resize(size, 0.0);
291 DCHECK(IsAllZero(v->values));
292 } else {
293 v->values.AssignToZero(size);
294 }
295 v->non_zeros.clear();
296}
297
298// Changes the sign of all the entries in the given vector.
299template <typename IndexType>
301 const IndexType end = data->size();
302 for (IndexType i(0); i < end; ++i) {
303 (*data)[i] = -(*data)[i];
304 }
305}
306
307// Given N Fractional elements, this class maintains their sum and can
308// provide, for each element X, the sum of all elements except X.
309// The subtelty is that it works well with infinities: for example, if there is
310// exactly one infinite element X, then SumWithout(X) will be finite.
311//
312// Two flavors of this class are provided: SumWithPositiveInfiniteAndOneMissing
313// supports calling Add() with normal numbers and positive infinities (and will
314// DCHECK() that), and SumWithNegativeInfiniteAndOneMissing does the same with
315// negative infinities.
316//
317// The numerical accuracy suffers however. If X is 1e100 and SumWithout(X)
318// should be 1e-100, then the value actually returned by SumWithout(X) is likely
319// to be wrong (by up to std::numeric_limits<Fractional>::epsilon() ^ 2).
320template <bool supported_infinity_is_positive>
322 public:
323 SumWithOneMissing() : num_infinities_(0), sum_() {}
324
325 void Add(Fractional x) {
326 if (num_infinities_ > 1) return;
327 if (IsFinite(x)) {
328 sum_.Add(x);
329 return;
330 }
331 DCHECK_EQ(Infinity(), x);
332 ++num_infinities_;
333 }
334
335 Fractional Sum() const {
336 if (num_infinities_ > 0) return Infinity();
337 return sum_.Value();
338 }
339
341 if (IsFinite(x)) {
342 if (num_infinities_ > 0) return Infinity();
343 return sum_.Value() - x;
344 }
345 DCHECK_EQ(Infinity(), x);
346 if (num_infinities_ > 1) return Infinity();
347 return sum_.Value();
348 }
349
350 private:
351 Fractional Infinity() const {
352 return supported_infinity_is_positive ? kInfinity : -kInfinity;
353 }
354
355 // Count how many times Add() was called with an infinite value. The count is
356 // stopped at 2 to be a bit faster.
357 int num_infinities_;
358 KahanSum sum_; // stripped of all the infinite values.
359};
362
363} // namespace glop
364} // namespace operations_research
365
366#endif // OR_TOOLS_LP_DATA_LP_UTILS_H_
#define DCHECK_LT(val1, val2)
Definition: base/logging.h:888
#define DCHECK(condition)
Definition: base/logging.h:884
#define DCHECK_EQ(val1, val2)
Definition: base/logging.h:885
void swap(StrongVector &x)
void Add(const FpNumber &value)
Definition: accurate_sum.h:29
void assign(IntType size, const T &v)
Definition: lp_types.h:274
Fractional SumWithout(Fractional x) const
int64 value
ColIndex col
Definition: markowitz.cc:176
RowIndex row
Definition: markowitz.cc:175
RowIndex ColToRowIndex(ColIndex col)
Definition: lp_types.h:51
bool IsFinite(Fractional value)
Definition: lp_types.h:90
Fractional PreciseSquaredNorm(const SparseColumn &v)
static Fractional Fractionality(Fractional f)
Fractional InfinityNorm(const DenseColumn &v)
Fractional SquaredNorm(const SparseColumn &v)
AccurateSum< Fractional > KahanSum
bool IsAllZero(const Container &input)
Fractional ScalarProduct(const DenseRowOrColumn1 &u, const DenseRowOrColumn2 &v)
void ComputeNonZeros(const StrictITIVector< IndexType, Fractional > &input, std::vector< IndexType > *non_zeros)
Fractional Square(Fractional f)
StrictITIVector< ColIndex, Fractional > DenseRow
Definition: lp_types.h:299
void RemoveNearZeroEntries(Fractional threshold, DenseRow *row)
SumWithOneMissing< false > SumWithNegativeInfiniteAndOneMissing
void PermuteWithKnownNonZeros(const Permutation< IndexType > &permutation, StrictITIVector< IndexType, Fractional > *zero_scratchpad, StrictITIVector< IndexType, Fractional > *output, std::vector< IndexType > *non_zeros)
double Density(const DenseRow &row)
void SetSupportToFalse(const ColumnView &column, DenseBooleanColumn *b)
bool IsDominated(const ColumnView &column, const DenseColumn &radius)
const DenseRow & Transpose(const DenseColumn &col)
Fractional PreciseScalarProduct(const DenseRowOrColumn &u, const DenseRowOrColumn2 &v)
Fractional PartialScalarProduct(const DenseRowOrColumn &u, const SparseColumn &v, int max_index)
Fractional RestrictedInfinityNorm(const ColumnView &column, const DenseBooleanColumn &rows_to_consider, RowIndex *row_index)
void ClearAndResizeVectorWithNonZeros(IndexType size, ScatteredRowOrCol *v)
bool IsAllFalse(const BoolVector &v)
StrictITIVector< RowIndex, Fractional > DenseColumn
Definition: lp_types.h:328
void ChangeSign(StrictITIVector< IndexType, Fractional > *data)
SumWithOneMissing< true > SumWithPositiveInfiniteAndOneMissing
StrictITIVector< RowIndex, bool > DenseBooleanColumn
Definition: lp_types.h:331
const double kInfinity
Definition: lp_types.h:83
void PermuteWithScratchpad(const Permutation< PermutationIndexType > &permutation, StrictITIVector< IndexType, Fractional > *zero_scratchpad, StrictITIVector< IndexType, Fractional > *input_output)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
int index
Definition: pack.cc:508
static int input(yyscan_t yyscanner)
bool ShouldUseDenseIteration(double ratio_for_using_dense_representation) const
StrictITIVector< Index, Fractional > values