Line data Source code
1 : #pragma once
2 :
3 : #include "Matrix.hpp"
4 :
5 : /// @addtogroup MatVec
6 : /// @{
7 :
8 : /// Class that represents matrices that permute the rows or columns of other
9 : /// matrices.
10 : /// Stored in an efficient manner with O(n) memory requirements.
11 : class PermutationMatrix {
12 :
13 : /// Container to store the elements of the permutation matrix internally.
14 : using storage_t = util::storage_t<size_t>;
15 :
16 : public:
17 : enum Type {
18 : Unspecified = 0, ///< Can be used for permuting rows or columns.
19 : RowPermutation = 1, ///< Can be used for permuting rows only.
20 : ColumnPermutation = 2, ///< Can be used for permuting columns only.
21 : };
22 :
23 : public:
24 : /// @name Constructors and assignment
25 : /// @{
26 :
27 : /// Default constructor.
28 : PermutationMatrix() = default;
29 :
30 : /// Create an empty permutation matrix with the given type.
31 14 : PermutationMatrix(Type type) : type(type) {}
32 :
33 : /// Create a permutation matrix without permutations.
34 3 : PermutationMatrix(size_t rows, Type type = Unspecified)
35 3 : : storage(rows), type(type) {
36 3 : fill_identity();
37 3 : }
38 :
39 : /// Create a permutation matrix with the given permutation.
40 : PermutationMatrix(std::initializer_list<size_t> init,
41 : Type type = Unspecified);
42 : /// Assign the given permutation to the matrix.
43 : PermutationMatrix &operator=(std::initializer_list<size_t> init);
44 :
45 : /// Default copy constructor.
46 5 : PermutationMatrix(const PermutationMatrix &) = default;
47 : /// Move constructor.
48 : PermutationMatrix(PermutationMatrix &&);
49 :
50 : /// Default copy assignment.
51 : PermutationMatrix &operator=(const PermutationMatrix &) = default;
52 : /// Move assignment.
53 : PermutationMatrix &operator=(PermutationMatrix &&);
54 :
55 : /// @}
56 :
57 : public:
58 : /// @name Matrix size
59 : /// @{
60 :
61 : /// Get the size of the permutation matrix.
62 1309 : size_t size() const { return storage.size(); }
63 : /// Get the number of rows of the permutation matrix.
64 : size_t rows() const { return size(); }
65 : /// Get the number of columns of the permutation matrix.
66 : size_t cols() const { return size(); }
67 : /// Get the number of elements in the matrix:
68 : size_t num_elems() const { return size(); }
69 : /// Resize the permutation matrix.
70 32 : void resize(size_t size) { storage.resize(size); }
71 :
72 : /// @}
73 :
74 : public:
75 : /// @name Element access
76 : /// @{
77 :
78 : /// Get the element at the given position in the swap sequence.
79 : /// If the k-th element is i, that is `P(k) == i`, this means that the k-th
80 : /// step of the swapping algorithm will swap `i` and `k`.
81 2270 : size_t &operator()(size_t index) { return storage[index]; }
82 : /// @copydoc operator()(size_t)
83 4466 : const size_t &operator()(size_t index) const { return storage[index]; }
84 :
85 : /// @}
86 :
87 : public:
88 : /// @name Transposition
89 : /// @{
90 :
91 : /// Reverse the order of the permutations.
92 5 : void reverse() { reverse_ = !reverse_; }
93 :
94 : /// Transpose or invert the permutation matrix.
95 5 : void transpose_inplace() { reverse(); }
96 :
97 : /// Check if the permutation should be applied in reverse.
98 26 : bool is_reversed() const { return reverse_; }
99 :
100 : /// Get the type of permutation matrix (whether it permutes rows or columns,
101 : /// or unspecified).
102 23 : Type get_type() const { return type; }
103 : /// Set the type of permutation matrix (whether it permutes rows or columns,
104 : /// or unspecified).
105 : void set_type(Type type) { this->type = type; }
106 :
107 : /// @}
108 :
109 : public:
110 : /// @name Conversion to a full matrix or a permutation
111 : /// @{
112 :
113 : /// Convert a permutation matrix into a full matrix.
114 : SquareMatrix to_matrix(Type type = Unspecified) const;
115 :
116 : /// Type that represents a permutation (in the mathematical sense, a
117 : /// permutation of the integers 0 through n-1).
118 : using Permutation = storage_t;
119 :
120 : /// Convert a permutation matrix into a mathematical permutation
121 : Permutation to_permutation() const;
122 :
123 : /// @}
124 :
125 : public:
126 : /// @name Applying the permutation to matrices
127 : /// @{
128 :
129 : /// Apply the permutation to the columns of matrix A.
130 : void permute_columns(Matrix &A) const;
131 : /// Apply the permutation to the rows of matrix A.
132 : void permute_rows(Matrix &A) const;
133 :
134 : /// @}
135 :
136 : public:
137 : /// @name Memory management
138 : /// @{
139 :
140 : /// Set the size to zero, and deallocate the storage.
141 : void clear_and_deallocate() {
142 : storage_t().swap(this->storage); // replace storage with empty storage
143 : // temporary storage goes out of scope and deallocates original storage
144 : }
145 :
146 : /// @}
147 :
148 : public:
149 : /// @name Generating permutations
150 : /// @{
151 :
152 : #ifndef NO_RANDOM_SUPPORT
153 : /// Return a random permutation of the integers 0 through length-1.
154 : static Permutation
155 : random_permutation(size_t length,
156 : std::default_random_engine::result_type seed =
157 : std::default_random_engine::default_seed);
158 : #endif
159 :
160 : /// Return the identity permutation (0, 1, 2, 3, ..., length-1).
161 : static Permutation identity_permutation(size_t length);
162 :
163 : /// @}
164 :
165 : public:
166 : /// @name Filling matrices
167 : /// @{
168 :
169 : /// Fill the matrix as an identity permutation.
170 17 : void fill_identity() { std::iota(begin(), end(), size_t(0)); }
171 :
172 : /// Create a permutation matrix from the given permutation.
173 : /// @note This isn't a very fast method, it's mainly used for tests.
174 : /// Internally, the permutation matrix is represented by a sequence
175 : /// of swap operations. Converting from this representation to a
176 : /// mathematical permutation is fast, but the other way around
177 : /// requires O(n²) operations (with the naive implementation used
178 : /// here, anyway).
179 : void fill_from_permutation(Permutation permutation);
180 :
181 : #ifndef NO_RANDOM_SUPPORT
182 : /// Fill the matrix with a random permutation.
183 : /// @note This isn't a very fast method, it's mainly used for tests.
184 1 : void fill_random(std::default_random_engine::result_type seed =
185 : std::default_random_engine::default_seed) {
186 1 : fill_from_permutation(random_permutation(size(), seed));
187 1 : }
188 : #endif
189 :
190 : /// @}
191 :
192 : public:
193 : /// @name Create special matrices
194 : /// @{
195 :
196 : /// Create an identity permutation matrix.
197 : static PermutationMatrix identity(size_t rows, Type type = Unspecified) {
198 : PermutationMatrix p(rows, type);
199 : return p;
200 : }
201 :
202 : /// @copydoc fill_from_permutation
203 2 : static PermutationMatrix from_permutation(Permutation permutation,
204 : Type type = Unspecified) {
205 2 : PermutationMatrix p(permutation.size(), type);
206 2 : p.fill_from_permutation(std::move(permutation));
207 2 : return p;
208 : }
209 :
210 : #ifndef NO_RANDOM_SUPPORT
211 : /// Create a random permutation matrix.
212 : /// @note This isn't a very fast method, it's mainly used for tests.
213 : static PermutationMatrix
214 1 : random(size_t rows, Type type = Unspecified,
215 : std::default_random_engine::result_type seed =
216 : std::default_random_engine::default_seed) {
217 1 : PermutationMatrix p(rows, type);
218 1 : p.fill_random(seed);
219 1 : return p;
220 : }
221 : #endif
222 :
223 : /// @}
224 :
225 : public:
226 : /// @name Iterators
227 : /// @{
228 :
229 : /// Get the iterator to the first element of the swapping sequence.
230 17 : storage_t::iterator begin() { return storage.begin(); }
231 : /// Get the iterator to the first element of the swapping sequence.
232 : storage_t::const_iterator begin() const { return storage.begin(); }
233 : /// Get the iterator to the first element of the swapping sequence.
234 : storage_t::const_iterator cbegin() const { return storage.begin(); }
235 :
236 : /// Get the iterator to the element past the end of the swapping sequence.
237 17 : storage_t::iterator end() { return storage.end(); }
238 : /// Get the iterator to the element past the end of the swapping sequence.
239 : storage_t::const_iterator end() const { return storage.end(); }
240 : /// Get the iterator to the element past the end of the swapping sequence.
241 : storage_t::const_iterator cend() const { return storage.end(); }
242 :
243 : /// @}
244 :
245 : public:
246 : /// @name Printing
247 : /// @{
248 :
249 : /// Print a permutation matrix.
250 : /// @param os
251 : /// The stream to print to.
252 : /// @param precision
253 : /// The number of significant figures to print.
254 : /// (0 = auto)
255 : /// @param width
256 : /// The width of each element (number of characters).
257 : /// (0 = auto)
258 : void print(std::ostream &os, uint8_t precision = 0,
259 : uint8_t width = 0) const;
260 :
261 : /// @}
262 :
263 : protected:
264 : storage_t storage;
265 : bool reverse_ = false;
266 : Type type = Unspecified;
267 : };
268 :
269 : /// @}
270 :
271 : /// Print a permutation matrix.
272 : /// @related PermutationMatrix
273 : std::ostream &operator<<(std::ostream &os, const PermutationMatrix &M);
274 :
275 : // :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: //
276 :
277 : /// @addtogroup MatMul
278 : /// @{
279 :
280 : /// Left application of permutation matrix (P permutes rows of A).
281 3 : inline Matrix operator*(const PermutationMatrix &P, const Matrix &A) {
282 3 : Matrix result = A;
283 3 : P.permute_rows(result);
284 3 : return result;
285 : }
286 : /// Left application of permutation matrix (P permutes rows of A).
287 1 : inline Matrix &&operator*(const PermutationMatrix &P, Matrix &&A) {
288 1 : P.permute_rows(A);
289 1 : return std::move(A);
290 : }
291 : /// Right application of permutation matrix (P permutes columns of A).
292 2 : inline Matrix operator*(const Matrix &A, const PermutationMatrix &P) {
293 2 : Matrix result = A;
294 2 : P.permute_columns(result);
295 2 : return result;
296 : }
297 : /// Right application of permutation matrix (P permutes columns of A).
298 1 : inline Matrix &&operator*(Matrix &&A, const PermutationMatrix &P) {
299 1 : P.permute_columns(A);
300 1 : return std::move(A);
301 : }
302 :
303 : /// Left application of permutation matrix (P permutes rows of A).
304 2 : inline SquareMatrix operator*(const PermutationMatrix &P,
305 : const SquareMatrix &A) {
306 2 : SquareMatrix result = A;
307 2 : P.permute_rows(result);
308 2 : return result;
309 : }
310 : /// Left application of permutation matrix (P permutes rows of A).
311 3 : inline SquareMatrix &&operator*(const PermutationMatrix &P, SquareMatrix &&A) {
312 3 : P.permute_rows(A);
313 3 : return std::move(A);
314 : }
315 : /// Right application of permutation matrix (P permutes columns of A).
316 1 : inline SquareMatrix operator*(const SquareMatrix &A,
317 : const PermutationMatrix &P) {
318 1 : SquareMatrix result = A;
319 1 : P.permute_columns(result);
320 1 : return result;
321 : }
322 : /// Right application of permutation matrix (P permutes columns of A).
323 1 : inline SquareMatrix &&operator*(SquareMatrix &&A, const PermutationMatrix &P) {
324 1 : P.permute_columns(A);
325 1 : return std::move(A);
326 : }
327 :
328 : /// Left application of permutation matrix (P permutes rows of v).
329 1 : inline Vector operator*(const PermutationMatrix &P, const Vector &v) {
330 1 : Vector result = v;
331 1 : P.permute_rows(result);
332 1 : return result;
333 : }
334 : /// Left application of permutation matrix (P permutes rows of v).
335 1 : inline Vector &&operator*(const PermutationMatrix &P, Vector &&v) {
336 1 : P.permute_rows(v);
337 1 : return std::move(v);
338 : }
339 :
340 : /// Right application of permutation matrix (P permutes columns of v).
341 1 : inline RowVector operator*(const RowVector &v, const PermutationMatrix &P) {
342 1 : RowVector result = v;
343 1 : P.permute_columns(result);
344 1 : return result;
345 : }
346 : /// Right application of permutation matrix (P permutes columns of v).
347 1 : inline RowVector &&operator*(RowVector &&v, const PermutationMatrix &P) {
348 1 : P.permute_columns(v);
349 1 : return std::move(v);
350 : }
351 :
352 : /// @}
353 :
354 : /// @addtogroup MatTrans
355 : /// @{
356 :
357 : /// Transpose a permutation matrix (inverse permutation).
358 3 : inline PermutationMatrix transpose(const PermutationMatrix &P) {
359 3 : PermutationMatrix result = P;
360 3 : result.transpose_inplace();
361 3 : return result;
362 : }
363 : /// Transpose a permutation matrix (inverse permutation).
364 2 : inline PermutationMatrix &&transpose(PermutationMatrix &&P) {
365 2 : P.transpose_inplace();
366 2 : return std::move(P);
367 : }
368 :
369 : /// @}
370 :
371 : // Implementations //
372 : // -------------------------------------------------------------------------- //
373 :
374 : inline PermutationMatrix::PermutationMatrix(PermutationMatrix &&other) {
375 : *this = std::move(other);
376 : }
377 :
378 : inline PermutationMatrix &
379 : PermutationMatrix::operator=(PermutationMatrix &&other) {
380 : // By explicitly defining move assignment, we can be sure that the object
381 : // that's being moved from has a consistent state.
382 : this->storage = std::move(other.storage);
383 : std::swap(this->type, other.type);
384 : std::swap(this->reverse_, other.reverse_);
385 : other.clear_and_deallocate();
386 : return *this;
387 : }
388 :
389 15 : inline PermutationMatrix::PermutationMatrix(std::initializer_list<size_t> init,
390 15 : Type type)
391 15 : : type(type) {
392 15 : *this = init;
393 15 : }
394 :
395 : inline SquareMatrix PermutationMatrix::to_matrix(Type type) const {
396 : // TODO: I'm sure this can be sped up
397 : Type actual_type = type == Unspecified ? this->type : type;
398 : assert(actual_type != Unspecified);
399 : if (actual_type == RowPermutation) {
400 : SquareMatrix P = SquareMatrix::identity(size());
401 : permute_rows(P);
402 : return P;
403 : } else if (actual_type == ColumnPermutation) {
404 : SquareMatrix P = SquareMatrix::identity(size());
405 : permute_columns(P);
406 : return P;
407 : }
408 : assert(false);
409 : return {};
410 : }
411 :
412 : inline PermutationMatrix::Permutation
413 3 : PermutationMatrix::to_permutation() const {
414 3 : Permutation p = identity_permutation(size());
415 3 : auto &This = *this;
416 3 : if (is_reversed()) {
417 : // Count down
418 1025 : for (size_t i = size(); i-- > 0;)
419 1024 : if (i != This(i))
420 1018 : std::swap(p[i], p[This(i)]);
421 : } else {
422 : // Count up
423 1154 : for (size_t i = 0; i < size(); ++i)
424 1152 : if (i != This(i))
425 1139 : std::swap(p[i], p[This(i)]);
426 : }
427 3 : return p;
428 : }
429 :
430 : inline PermutationMatrix &
431 15 : PermutationMatrix::operator=(std::initializer_list<size_t> init) {
432 30 : storage_t permutation(init.size());
433 15 : std::copy(init.begin(), init.end(), permutation.begin());
434 15 : fill_from_permutation(std::move(permutation));
435 :
436 30 : return *this;
437 : }
438 :
439 7 : inline void PermutationMatrix::permute_columns(Matrix &A) const {
440 7 : assert(A.cols() == size());
441 7 : assert(get_type() != RowPermutation);
442 7 : auto &This = *this;
443 7 : if (is_reversed()) {
444 : // Count down
445 5 : for (size_t i = size(); i-- > 0;)
446 4 : if (i != This(i))
447 2 : A.swap_columns(i, This(i));
448 : } else {
449 : // Count up
450 30 : for (size_t i = 0; i < size(); ++i)
451 24 : if (i != This(i))
452 12 : A.swap_columns(i, This(i));
453 : }
454 7 : }
455 :
456 16 : inline void PermutationMatrix::permute_rows(Matrix &A) const {
457 16 : assert(A.rows() == size());
458 16 : assert(get_type() != ColumnPermutation);
459 16 : auto &This = *this;
460 16 : if (is_reversed()) {
461 : // Count down
462 16 : for (size_t i = size(); i-- > 0;)
463 13 : if (i != This(i))
464 7 : A.swap_rows(i, This(i));
465 : } else {
466 : // Count up
467 61 : for (size_t i = 0; i < size(); ++i)
468 48 : if (i != This(i))
469 23 : A.swap_rows(i, This(i));
470 : }
471 16 : }
|