LCOV - code coverage report
Current view: top level - src/include/linalg - PermutationMatrix.hpp (source / functions) Hit Total Coverage
Test: 98c60408be2d766eef282c11e0ebd9afde851352 Lines: 122 122 100.0 %
Date: 2021-02-20 15:44:35 Functions: 36 36 100.0 %
Legend: Lines: hit not hit

          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 : }

Generated by: LCOV version 1.15