PermutationMatrix.hpp

Go to the documentation of this file.

#pragma once
 
#include "Matrix.hpp"
 
 
class PermutationMatrix {
 
    using storage_t = util::storage_t<size_t>;
 
  public:
    enum Type {
        Unspecified = 0,       
        RowPermutation = 1,    
        ColumnPermutation = 2, 
    };
 
  public:
 
    PermutationMatrix() = default;
 
    PermutationMatrix(Type type) : type(type) {}
 
    PermutationMatrix(size_t rows, Type type = Unspecified)
        : storage(rows), type(type) {
        fill_identity();
    }
 
    PermutationMatrix(std::initializer_list<size_t> init,
                      Type type = Unspecified);
    PermutationMatrix &operator=(std::initializer_list<size_t> init);
 
    PermutationMatrix(const PermutationMatrix &) = default;
    PermutationMatrix(PermutationMatrix &&);
 
    PermutationMatrix &operator=(const PermutationMatrix &) = default;
    PermutationMatrix &operator=(PermutationMatrix &&);
 
 
  public:
 
    size_t size() const { return storage.size(); }
    size_t rows() const { return size(); }
    size_t cols() const { return size(); }
    size_t num_elems() const { return size(); }
    void resize(size_t size) { storage.resize(size); }
 
 
  public:
 
    size_t &operator()(size_t index) { return storage[index]; }
    const size_t &operator()(size_t index) const { return storage[index]; }
 
 
  public:
 
    void reverse() { reverse_ = !reverse_; }
 
    void transpose_inplace() { reverse(); }
 
    bool is_reversed() const { return reverse_; }
 
    Type get_type() const { return type; }
    void set_type(Type type) { this->type = type; }
 
 
  public:
 
    SquareMatrix to_matrix(Type type = Unspecified) const;
 
    using Permutation = storage_t;
 
    Permutation to_permutation() const;
 
 
  public:
 
    void permute_columns(Matrix &A) const;
    void permute_rows(Matrix &A) const;
 
 
  public:
 
    void clear_and_deallocate() {
        storage_t().swap(this->storage); // replace storage with empty storage
        // temporary storage goes out of scope and deallocates original storage
    }
 
 
  public:
 
#ifndef NO_RANDOM_SUPPORT
    static Permutation
    random_permutation(size_t length,
                       std::default_random_engine::result_type seed =
                           std::default_random_engine::default_seed);
#endif
 
    static Permutation identity_permutation(size_t length);
 
 
  public:
 
    void fill_identity() { std::iota(begin(), end(), size_t(0)); }
 
    void fill_from_permutation(Permutation permutation);
 
#ifndef NO_RANDOM_SUPPORT
    void fill_random(std::default_random_engine::result_type seed =
                         std::default_random_engine::default_seed) {
        fill_from_permutation(random_permutation(size(), seed));
    }
#endif
 
 
  public:
 
    static PermutationMatrix identity(size_t rows, Type type = Unspecified) {
        PermutationMatrix p(rows, type);
        return p;
    }
 
    static PermutationMatrix from_permutation(Permutation permutation,
                                              Type type = Unspecified) {
        PermutationMatrix p(permutation.size(), type);
        p.fill_from_permutation(std::move(permutation));
        return p;
    }
 
#ifndef NO_RANDOM_SUPPORT
    static PermutationMatrix
    random(size_t rows, Type type = Unspecified,
           std::default_random_engine::result_type seed =
               std::default_random_engine::default_seed) {
        PermutationMatrix p(rows, type);
        p.fill_random(seed);
        return p;
    }
#endif
 
 
  public:
 
    storage_t::iterator begin() { return storage.begin(); }
    storage_t::const_iterator begin() const { return storage.begin(); }
    storage_t::const_iterator cbegin() const { return storage.begin(); }
 
    storage_t::iterator end() { return storage.end(); }
    storage_t::const_iterator end() const { return storage.end(); }
    storage_t::const_iterator cend() const { return storage.end(); }
 
 
  public:
 
    void print(std::ostream &os, uint8_t precision = 0,
               uint8_t width = 0) const;
 
 
  protected:
    storage_t storage;
    bool reverse_ = false;
    Type type = Unspecified;
};
 
 
std::ostream &operator<<(std::ostream &os, const PermutationMatrix &M);
 
// :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: //
 
 
inline Matrix operator*(const PermutationMatrix &P, const Matrix &A) {
    Matrix result = A;
    P.permute_rows(result);
    return result;
}
inline Matrix &&operator*(const PermutationMatrix &P, Matrix &&A) {
    P.permute_rows(A);
    return std::move(A);
}
inline Matrix operator*(const Matrix &A, const PermutationMatrix &P) {
    Matrix result = A;
    P.permute_columns(result);
    return result;
}
inline Matrix &&operator*(Matrix &&A, const PermutationMatrix &P) {
    P.permute_columns(A);
    return std::move(A);
}
 
inline SquareMatrix operator*(const PermutationMatrix &P,
                              const SquareMatrix &A) {
    SquareMatrix result = A;
    P.permute_rows(result);
    return result;
}
inline SquareMatrix &&operator*(const PermutationMatrix &P, SquareMatrix &&A) {
    P.permute_rows(A);
    return std::move(A);
}
inline SquareMatrix operator*(const SquareMatrix &A,
                              const PermutationMatrix &P) {
    SquareMatrix result = A;
    P.permute_columns(result);
    return result;
}
inline SquareMatrix &&operator*(SquareMatrix &&A, const PermutationMatrix &P) {
    P.permute_columns(A);
    return std::move(A);
}
 
inline Vector operator*(const PermutationMatrix &P, const Vector &v) {
    Vector result = v;
    P.permute_rows(result);
    return result;
}
inline Vector &&operator*(const PermutationMatrix &P, Vector &&v) {
    P.permute_rows(v);
    return std::move(v);
}
 
inline RowVector operator*(const RowVector &v, const PermutationMatrix &P) {
    RowVector result = v;
    P.permute_columns(result);
    return result;
}
inline RowVector &&operator*(RowVector &&v, const PermutationMatrix &P) {
    P.permute_columns(v);
    return std::move(v);
}
 
 
 
inline PermutationMatrix transpose(const PermutationMatrix &P) {
    PermutationMatrix result = P;
    result.transpose_inplace();
    return result;
}
inline PermutationMatrix &&transpose(PermutationMatrix &&P) {
    P.transpose_inplace();
    return std::move(P);
}
 
 
//                              Implementations                               //
// -------------------------------------------------------------------------- //
 
inline PermutationMatrix::PermutationMatrix(PermutationMatrix &&other) {
    *this = std::move(other);
}
 
inline PermutationMatrix &
PermutationMatrix::operator=(PermutationMatrix &&other) {
    // By explicitly defining move assignment, we can be sure that the object
    // that's being moved from has a consistent state.
    this->storage = std::move(other.storage);
    std::swap(this->type, other.type);
    std::swap(this->reverse_, other.reverse_);
    other.clear_and_deallocate();
    return *this;
}
 
inline PermutationMatrix::PermutationMatrix(std::initializer_list<size_t> init,
                                            Type type)
    : type(type) {
    *this = init;
}
 
inline SquareMatrix PermutationMatrix::to_matrix(Type type) const {
    // TODO: I'm sure this can be sped up
    Type actual_type = type == Unspecified ? this->type : type;
    assert(actual_type != Unspecified);
    if (actual_type == RowPermutation) {
        SquareMatrix P = SquareMatrix::identity(size());
        permute_rows(P);
        return P;
    } else if (actual_type == ColumnPermutation) {
        SquareMatrix P = SquareMatrix::identity(size());
        permute_columns(P);
        return P;
    }
    assert(false);
    return {};
}
 
inline PermutationMatrix::Permutation
PermutationMatrix::to_permutation() const {
    Permutation p = identity_permutation(size());
    auto &This = *this;
    if (is_reversed()) {
        // Count down
        for (size_t i = size(); i-- > 0;)
            if (i != This(i))
                std::swap(p[i], p[This(i)]);
    } else {
        // Count up
        for (size_t i = 0; i < size(); ++i)
            if (i != This(i))
                std::swap(p[i], p[This(i)]);
    }
    return p;
}
 
inline PermutationMatrix &
PermutationMatrix::operator=(std::initializer_list<size_t> init) {
    storage_t permutation(init.size());
    std::copy(init.begin(), init.end(), permutation.begin());
    fill_from_permutation(std::move(permutation));
 
    return *this;
}
 
inline void PermutationMatrix::permute_columns(Matrix &A) const {
    assert(A.cols() == size());
    assert(get_type() != RowPermutation);
    auto &This = *this;
    if (is_reversed()) {
        // Count down
        for (size_t i = size(); i-- > 0;)
            if (i != This(i))
                A.swap_columns(i, This(i));
    } else {
        // Count up
        for (size_t i = 0; i < size(); ++i)
            if (i != This(i))
                A.swap_columns(i, This(i));
    }
}
 
inline void PermutationMatrix::permute_rows(Matrix &A) const {
    assert(A.rows() == size());
    assert(get_type() != ColumnPermutation);
    auto &This = *this;
    if (is_reversed()) {
        // Count down
        for (size_t i = size(); i-- > 0;)
            if (i != This(i))
                A.swap_rows(i, This(i));
    } else {
        // Count up
        for (size_t i = 0; i < size(); ++i)
            if (i != This(i))
                A.swap_rows(i, This(i));
    }
}