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 
24 namespace operations_research {
25 namespace 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.
36 inline 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.
46 template <class DenseRowOrColumn1, class DenseRowOrColumn2>
47 Fractional 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.
81 template <class DenseRowOrColumn>
82 Fractional 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 
91 template <class DenseRowOrColumn, class DenseRowOrColumn2>
92 Fractional 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 
102 template <class DenseRowOrColumn>
103 Fractional 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 
113 template <class DenseRowOrColumn>
114 Fractional 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.
129 template <class DenseRowOrColumn>
130 Fractional 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.
145 Fractional SquaredNorm(const SparseColumn& v);
146 Fractional SquaredNorm(const DenseColumn& column);
147 Fractional SquaredNorm(const ColumnView& v);
148 Fractional PreciseSquaredNorm(const SparseColumn& v);
150 Fractional PreciseSquaredNorm(const ScatteredColumn& v);
151 
152 // Returns the maximum of the |coefficients| of 'v'.
154 Fractional InfinityNorm(const SparseColumn& v);
155 Fractional 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.
161 double 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.
165 void RemoveNearZeroEntries(Fractional threshold, DenseRow* row);
166 void RemoveNearZeroEntries(Fractional threshold, DenseColumn* column);
167 
168 // Transposition functions implemented below with a cast so it should actually
169 // have no complexity cost.
170 const DenseRow& Transpose(const DenseColumn& col);
171 const 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.
176 Fractional 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'.
183 void SetSupportToFalse(const ColumnView& column, DenseBooleanColumn* b);
184 
185 // Returns true iff for all 'row' we have '|column[row]| <= radius[row]'.
186 bool 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.
192 inline 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.
200 inline 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.
208 template <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.
221 template <typename Container>
222 inline 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.
230 template <typename BoolVector>
231 bool 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.
236 template <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.
258 template <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.
277 template <typename IndexType, typename ScatteredRowOrCol>
278 inline 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.
299 template <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).
320 template <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
void PermuteWithScratchpad(const Permutation< PermutationIndexType > &permutation, StrictITIVector< IndexType, Fractional > *zero_scratchpad, StrictITIVector< IndexType, Fractional > *input_output)
Fractional Square(Fractional f)
Fractional PreciseSquaredNorm(const SparseColumn &v)
Fractional InfinityNorm(const DenseColumn &v)
Fractional SquaredNorm(const SparseColumn &v)
bool IsAllZero(const Container &input)
AccurateSum< Fractional > KahanSum
Fractional ScalarProduct(const DenseRowOrColumn1 &u, const DenseRowOrColumn2 &v)
void ComputeNonZeros(const StrictITIVector< IndexType, Fractional > &input, std::vector< IndexType > *non_zeros)
Fractional PreciseScalarProduct(const DenseRowOrColumn &u, const DenseRowOrColumn2 &v)
StrictITIVector< ColIndex, Fractional > DenseRow
Definition: lp_types.h:299
void RemoveNearZeroEntries(Fractional threshold, DenseRow *row)
SumWithOneMissing< false > SumWithNegativeInfiniteAndOneMissing
bool IsAllFalse(const BoolVector &v)
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 IsFinite(Fractional value)
Definition: lp_types.h:90
bool IsDominated(const ColumnView &column, const DenseColumn &radius)
void ClearAndResizeVectorWithNonZeros(IndexType size, ScatteredRowOrCol *v)
const DenseRow & Transpose(const DenseColumn &col)
RowIndex ColToRowIndex(ColIndex col)
Definition: lp_types.h:51
void ChangeSign(StrictITIVector< IndexType, Fractional > *data)
static Fractional Fractionality(Fractional f)
Fractional PartialScalarProduct(const DenseRowOrColumn &u, const SparseColumn &v, int max_index)
Fractional RestrictedInfinityNorm(const ColumnView &column, const DenseBooleanColumn &rows_to_consider, RowIndex *row_index)
StrictITIVector< RowIndex, Fractional > DenseColumn
Definition: lp_types.h:328
SumWithOneMissing< true > SumWithPositiveInfiniteAndOneMissing
StrictITIVector< RowIndex, bool > DenseBooleanColumn
Definition: lp_types.h:331
const double kInfinity
Definition: lp_types.h:83
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