OR-Tools  8.2
matrix_utils.cc
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
15
16#include <algorithm>
17
18#include "ortools/base/hash.h"
19
20namespace operations_research {
21namespace glop {
22
23namespace {
24
25// Returns true iff the two given sparse columns are proportional. The two
26// sparse columns must be ordered by row and must not contain any zero entry.
27//
28// See the header comment on FindProportionalColumns() for the exact definition
29// of two proportional columns with a given tolerance.
30bool AreColumnsProportional(const SparseColumn& a, const SparseColumn& b,
31 Fractional tolerance) {
32 DCHECK(a.IsCleanedUp());
33 DCHECK(b.IsCleanedUp());
34 if (a.num_entries() != b.num_entries()) return false;
35 Fractional multiple = 0.0;
36 bool a_is_larger = true;
37 for (const EntryIndex i : a.AllEntryIndices()) {
38 if (a.EntryRow(i) != b.EntryRow(i)) return false;
39 const Fractional coeff_a = a.EntryCoefficient(i);
40 const Fractional coeff_b = b.EntryCoefficient(i);
41 if (multiple == 0.0) {
42 a_is_larger = std::abs(coeff_a) > std::abs(coeff_b);
43 multiple = a_is_larger ? coeff_a / coeff_b : coeff_b / coeff_a;
44 } else {
45 if (a_is_larger) {
46 if (std::abs(coeff_a / coeff_b - multiple) > tolerance) return false;
47 } else {
48 if (std::abs(coeff_b / coeff_a - multiple) > tolerance) return false;
49 }
50 }
51 }
52 return true;
53}
54
55// A column index together with its fingerprint. See ComputeFingerprint().
56struct ColumnFingerprint {
57 ColumnFingerprint(ColIndex _col, int64 _hash, double _value)
58 : col(_col), hash(_hash), value(_value) {}
59 ColIndex col;
61 double value;
62
63 // This order has the property that if AreProportionalCandidates() is true for
64 // two given columns, then in a sorted list of columns
65 // AreProportionalCandidates() will be true for all the pairs of columns
66 // between the two given ones (included).
67 bool operator<(const ColumnFingerprint& other) const {
68 if (hash == other.hash) {
69 return value < other.value;
70 }
71 return hash < other.hash;
72 }
73};
74
75// Two columns can be proportional only if:
76// - Their non-zero pattern hashes are the same.
77// - Their double fingerprints are close to each other.
78bool AreProportionalCandidates(ColumnFingerprint a, ColumnFingerprint b,
79 Fractional tolerance) {
80 if (a.hash != b.hash) return false;
81 return std::abs(a.value - b.value) < tolerance;
82}
83
84// The fingerprint of a column has two parts:
85// - A hash value of the column non-zero pattern.
86// - A double value which should be the same for two proportional columns
87// modulo numerical errors.
88ColumnFingerprint ComputeFingerprint(ColIndex col, const SparseColumn& column) {
89 int64 non_zero_pattern_hash = 0;
91 Fractional max_abs = 0.0;
92 Fractional sum = 0.0;
93 for (const SparseColumn::Entry e : column) {
94 non_zero_pattern_hash =
95 util_hash::Hash(e.row().value(), non_zero_pattern_hash);
96 sum += e.coefficient();
97 min_abs = std::min(min_abs, std::abs(e.coefficient()));
98 max_abs = std::max(max_abs, std::abs(e.coefficient()));
99 }
100
101 // The two scaled values are in [0, 1].
102 // TODO(user): A better way to discriminate columns would be to take the
103 // scalar product with a constant but random vector scaled by max_abs.
104 DCHECK_NE(0.0, max_abs);
105 const double inverse_dynamic_range = min_abs / max_abs;
106 const double scaled_average =
107 std::abs(sum) /
108 (static_cast<double>(column.num_entries().value()) * max_abs);
109 return ColumnFingerprint(col, non_zero_pattern_hash,
110 inverse_dynamic_range + scaled_average);
111}
112
113} // namespace
114
116 Fractional tolerance) {
117 const ColIndex num_cols = matrix.num_cols();
118 ColMapping mapping(num_cols, kInvalidCol);
119
120 // Compute the fingerprint of each columns and sort them.
121 std::vector<ColumnFingerprint> fingerprints;
122 for (ColIndex col(0); col < num_cols; ++col) {
123 if (!matrix.column(col).IsEmpty()) {
124 fingerprints.push_back(ComputeFingerprint(col, matrix.column(col)));
125 }
126 }
127 std::sort(fingerprints.begin(), fingerprints.end());
128
129 // Find a representative of each proportional columns class. This only
130 // compares columns with a close-enough fingerprint.
131 for (int i = 0; i < fingerprints.size(); ++i) {
132 const ColIndex col_a = fingerprints[i].col;
133 if (mapping[col_a] != kInvalidCol) continue;
134 for (int j = i + 1; j < fingerprints.size(); ++j) {
135 const ColIndex col_b = fingerprints[j].col;
136 if (mapping[col_b] != kInvalidCol) continue;
137
138 // Note that we use the same tolerance for the fingerprints.
139 // TODO(user): Derive precise bounds on what this tolerance should be so
140 // that no proportional columns are missed.
141 if (!AreProportionalCandidates(fingerprints[i], fingerprints[j],
142 tolerance)) {
143 break;
144 }
145 if (AreColumnsProportional(matrix.column(col_a), matrix.column(col_b),
146 tolerance)) {
147 mapping[col_b] = col_a;
148 }
149 }
150 }
151
152 // Sort the mapping so that the representative of each class is the smallest
153 // column. To achieve this, the current representative is used as a pointer
154 // to the new one, a bit like in an union find algorithm.
155 for (ColIndex col(0); col < num_cols; ++col) {
156 if (mapping[col] == kInvalidCol) continue;
157 const ColIndex new_representative = mapping[mapping[col]];
158 if (new_representative != kInvalidCol) {
159 mapping[col] = new_representative;
160 } else {
161 if (mapping[col] > col) {
162 mapping[mapping[col]] = col;
163 mapping[col] = kInvalidCol;
164 }
165 }
166 }
167
168 return mapping;
169}
170
172 const SparseMatrix& matrix, Fractional tolerance) {
173 const ColIndex num_cols = matrix.num_cols();
174 ColMapping mapping(num_cols, kInvalidCol);
175 for (ColIndex col_a(0); col_a < num_cols; ++col_a) {
176 if (matrix.column(col_a).IsEmpty()) continue;
177 if (mapping[col_a] != kInvalidCol) continue;
178 for (ColIndex col_b(col_a + 1); col_b < num_cols; ++col_b) {
179 if (matrix.column(col_b).IsEmpty()) continue;
180 if (mapping[col_b] != kInvalidCol) continue;
181 if (AreColumnsProportional(matrix.column(col_a), matrix.column(col_b),
182 tolerance)) {
183 mapping[col_b] = col_a;
184 }
185 }
186 }
187 return mapping;
188}
189
190bool AreFirstColumnsAndRowsExactlyEquals(RowIndex num_rows, ColIndex num_cols,
191 const SparseMatrix& matrix_a,
192 const CompactSparseMatrix& matrix_b) {
193 // TODO(user): Also DCHECK() that matrix_b is ordered by rows.
194 DCHECK(matrix_a.IsCleanedUp());
195 if (num_rows > matrix_a.num_rows() || num_rows > matrix_b.num_rows() ||
196 num_cols > matrix_a.num_cols() || num_cols > matrix_b.num_cols()) {
197 return false;
198 }
199 for (ColIndex col(0); col < num_cols; ++col) {
200 const SparseColumn& col_a = matrix_a.column(col);
201 const ColumnView& col_b = matrix_b.column(col);
202 const EntryIndex end = std::min(col_a.num_entries(), col_b.num_entries());
203 if (end < col_a.num_entries() && col_a.EntryRow(end) < num_rows) {
204 return false;
205 }
206 if (end < col_b.num_entries() && col_b.EntryRow(end) < num_rows) {
207 return false;
208 }
209 for (EntryIndex i(0); i < end; ++i) {
210 if (col_a.EntryRow(i) != col_b.EntryRow(i)) {
211 if (col_a.EntryRow(i) < num_rows || col_b.EntryRow(i) < num_rows) {
212 return false;
213 } else {
214 break;
215 }
216 }
217 if (col_a.EntryCoefficient(i) != col_b.EntryCoefficient(i)) {
218 return false;
219 }
220 if (col_a.num_entries() > end && col_a.EntryRow(end) < num_rows) {
221 return false;
222 }
223 if (col_b.num_entries() > end && col_b.EntryRow(end) < num_rows) {
224 return false;
225 }
226 }
227 }
228 return true;
229}
230
232 DCHECK(matrix.IsCleanedUp());
233 if (matrix.num_rows().value() > matrix.num_cols().value()) return false;
234 const ColIndex first_identity_col =
235 matrix.num_cols() - RowToColIndex(matrix.num_rows());
236 for (ColIndex col = first_identity_col; col < matrix.num_cols(); ++col) {
237 const SparseColumn& column = matrix.column(col);
238 if (column.num_entries() != 1 ||
239 column.EntryCoefficient(EntryIndex(0)) != 1.0) {
240 return false;
241 }
242 }
243 return true;
244}
245
246} // namespace glop
247} // namespace operations_research
int64 min
Definition: alldiff_cst.cc:138
int64 max
Definition: alldiff_cst.cc:139
#define DCHECK_NE(val1, val2)
Definition: base/logging.h:886
#define DCHECK(condition)
Definition: base/logging.h:884
Fractional EntryCoefficient(EntryIndex i) const
Definition: sparse_column.h:83
RowIndex EntryRow(EntryIndex i) const
Definition: sparse_column.h:89
ColumnView column(ColIndex col) const
Definition: sparse.h:364
Fractional EntryCoefficient(EntryIndex i) const
Definition: sparse_column.h:52
RowIndex EntryRow(EntryIndex i) const
Definition: sparse_column.h:51
const SparseColumn & column(ColIndex col) const
Definition: sparse.h:180
int64_t int64
ColIndex col
Definition: matrix_utils.cc:59
int64 hash
Definition: matrix_utils.cc:60
double value
Definition: matrix_utils.cc:61
bool IsRightMostSquareMatrixIdentity(const SparseMatrix &matrix)
ColMapping FindProportionalColumnsUsingSimpleAlgorithm(const SparseMatrix &matrix, Fractional tolerance)
const ColIndex kInvalidCol(-1)
ColIndex RowToColIndex(RowIndex row)
Definition: lp_types.h:48
bool AreFirstColumnsAndRowsExactlyEquals(RowIndex num_rows, ColIndex num_cols, const SparseMatrix &matrix_a, const CompactSparseMatrix &matrix_b)
ColMapping FindProportionalColumns(const SparseMatrix &matrix, Fractional tolerance)
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
uint64 Hash(uint64 num, uint64 c)
Definition: hash.h:150