matrix.h

Go to the documentation of this file.

#ifndef MU_MATRIX_H_
#define MU_MATRIX_H_

#include <algorithm>
#include <array>
#include <cassert>
#include <type_traits>
#include <vector>

#include "mu/typetraits.h"
#include "mu/utility.h"
#include "vector.h"

namespace mu {

template <std::size_t N, std::size_t M, typename T>
class Matrix {
  static_assert(N != 0, "first matrix dimension (rows) cannot be zero");
  static_assert(M != 0, "second matrix dimension (columns) cannot be zero");
  static_assert(std::is_arithmetic<T>::value,
                "Matrix type T must be an arithmetic type");

 public:
  /* value and size type from the underlying container */
  using value_type = typename std::array<Vector<M, T>, N>::value_type;
  using size_type = typename std::array<Vector<M, T>, N>::size_type;
  /* use iterators from the underlying container */
  using iterator = typename std::array<Vector<M, T>, N>::iterator;
  using const_iterator = typename std::array<Vector<M, T>, N>::const_iterator;

  constexpr Matrix() = default;

  template <typename... TArgs,
            std::enable_if_t<
                sizeof...(TArgs) == N &&
                    (mu::conjunction<std::is_same<
                         T, std::remove_reference_t<TArgs>>::value...>::value),
                int> = 0>
  // NOLINTNEXTLINE(runtime/explicit) implicit to make copy-init. work
  Matrix(TArgs const(&&... rows)[M]) : data_{mu::to_array(rows)...} {}

  template <std::size_t Nn, std::size_t Mm, class U>
  // NOLINTNEXTLINE(runtime/explicit) implicit conversion is intentional
  Matrix(const Matrix<Nn, Mm, U> &m) {
    static_assert(N == Nn, "Matrix dimension mismatch (rows)");
    static_assert(M == Mm, "Matrix dimension mismatch (columns)");
    std::transform(m.begin(), m.end(), begin(),
                   [](Vector<M, U> data) { return data; });
  }

  // NOLINTNEXTLINE(runtime/explicit) implicit to make copy-init. work
  Matrix(const std::array<Vector<M, T>, N> &a) : data_{a} {}

  template <typename U = T>
  // NOLINTNEXTLINE(runtime/explicit) implicit to make copy-init. work
  Matrix(const std::array<Vector<M, U>, N> &a) {
    std::transform(a.begin(), a.end(), begin(),
                   [](Vector<M, U> data) { return data; });
  }

  template <typename U = T,
            std::enable_if_t<std::is_arithmetic<U>::value, int> = 0>
  // NOLINTNEXTLINE(runtime/explicit) implicit to make copy-init. work
  Matrix(const std::array<std::array<U, M>, N> &a) {
    std::transform(a.begin(), a.end(), begin(),
                   [](Vector<M, U> data) { return data; });
  }

  template <typename U = T,
            std::enable_if_t<std::is_arithmetic<U>::value, int> = 0>
  // NOLINTNEXTLINE(runtime/explicit) implicit to make copy-init. work
  Matrix(const U &value) {
    data_.fill(value);
  }

  ~Matrix() = default;

  Matrix(const Matrix &other) = default;

  Matrix(Matrix &&other) noexcept = default;

  Matrix &operator=(const Matrix &other) = default;

  Matrix &operator=(Matrix &&other) noexcept = default;

  Vector<M, T> &operator[](size_type idx) noexcept { return data_[idx]; }

  const Vector<M, T> &operator[](size_type idx) const noexcept {
    return data_[idx];
  }

  Vector<M, T> &at(size_type idx) { return data_.at(idx); }

  const Vector<M, T> &at(size_type idx) const { return data_.at(idx); }

  constexpr std::array<size_type, 2> size() const noexcept {
    return std::array<size_type, 2>{N, M};
  }

  constexpr size_type n_rows() const noexcept { return N; }

  constexpr size_type n_cols() const noexcept { return M; }

  iterator begin() noexcept { return data_.begin(); }

  const_iterator begin() const noexcept { return data_.begin(); }

  iterator end() noexcept { return data_.end(); }

  const_iterator end() const noexcept { return data_.end(); }

  Vector<M, T> row(const size_type &idx) const {
    /* runtime check for out-of-range. only in debug mode */
    assert(idx >= 0 && idx < N);
    return data_[idx];
  }

  Vector<N, T> col(const size_type &idx) const {
    /* runtime check for out-of-range. only in debug mode */
    assert(idx >= 0 && idx < M);
    Vector<N, T> ret;
    for (std::size_t i = 0; i < N; i++) {
      ret[i] = data_[i][idx];
    }
    return ret;
  }

  T min() const {
    T ret(data_[0][0]);
    for (const auto &row : data_) {
      ret = mu::min(ret, mu::min(row));
    }
    return ret;
  }

  T max() const {
    T ret(data_[0][0]);
    for (const auto &row : data_) {
      ret = mu::max(ret, mu::max(row));
    }
    return ret;
  }

  T sum() const {
    T ret{};
    for (const auto &row : data_) {
      ret += mu::sum(row);
    }
    return ret;
  }

  template <typename U = T>
  U mean() const {
    return U(sum()) / (N * M);
  }

  auto diag() const {
    constexpr std::size_t s = N < M ? N : M;
    Vector<s, T> ret;
    for (std::size_t i = 0; i < s; i++) {
      ret[i] = data_[i][i];
    }
    return ret;
  }

  T det() const {
    static_assert(N == M,
                  "Matrix dimensions must match to calculate the determinant");
    std::vector<std::vector<T>> vv;
    for (std::size_t i = 0; i < N; i++) {
      vv.push_back(std::vector<T>());
      for (const auto &item : data_[i]) {
        vv[i].push_back(item);
      }
    }
    return calc_det(vv);
  }

  void transpose() {
    static_assert(
        N == M,
        "Matrix dimensions must match to transpose this object. For Matrices "
        "with unequal dimensions, i.e. N != M, the \"transposed()\" method can "
        "be used instead");
    Matrix<N, M, T> copy = data_;
    for (std::size_t i = 0; i < N; i++) {
      for (std::size_t j = 0; j < M; j++) {
        data_[j][i] = copy[i][j];
      }
    }
  }

  Matrix<M, N, T> transposed() const {
    Matrix<M, N, T> ret;
    for (std::size_t i = 0; i < N; i++) {
      for (std::size_t j = 0; j < M; j++) {
        ret[j][i] = data_[i][j];
      }
    }
    return ret;
  }

  template <typename U = void, std::size_t N2, std::size_t M2, typename T2>
  std::conditional_t<std::is_same<U, void>::value, Matrix<N, M2, T>,
                     Matrix<N, M2, U>>
  dot(const Matrix<N2, M2, T2> &rhs) const {
    static_assert(
        M == N2,
        "Matrix dimension mismatch. Second dimension of first matrix must be "
        "equal to the first dimension of the second matrix");
    using U_ = std::conditional_t<!std::is_same<T, T2>::value, U, T>;
    static_assert(!std::is_same<U_, void>::value,
                  "Matrix types are different. please specify the return "
                  "type. e.g. \"mat1.dot<float>(mat2);\"");
    Matrix<N, M2, U_> ret;
    for (std::size_t i = 0; i < N; i++) {
      for (std::size_t j = 0; j < M2; j++) {
        U_ sum{0};
        for (std::size_t k = 0; k < M; k++) {
          sum += (data_[i][k] * rhs[k][j]);
        }
        ret[i][j] = sum;
      }
    }
    return ret;
  }

  template <typename U = void, std::size_t N2, typename T2>
  std::conditional_t<std::is_same<U, void>::value, Vector<N, T>, Vector<N, U>>
  dot(const Vector<N2, T2> &rhs) const {
    static_assert(
        M == N2,
        "Matrix-Vector dimension mismatch. Second dimension of the matrix "
        "must be equal to the vector size");
    using U_ = std::conditional_t<!std::is_same<T, T2>::value, U, T>;
    static_assert(
        !std::is_same<U_, void>::value,
        "Matrix and Vector types are different. please specify the return "
        "type. e.g. \"mat.dot<float>(vec);\"");
    Vector<N, U_> ret;
    for (std::size_t i = 0; i < N; i++) {
      U_ sum{0};
      for (std::size_t k = 0; k < M; k++) {
        sum += (data_[i][k] * rhs[k]);
      }
      ret[i] = sum;
    }
    return ret;
  }

  template <class U = T>
  U std() const {
    U sum{0};
    U m = mean<U>();
    for (const auto &row : data_) {
      for (const auto &item : row) {
        sum += mu::pow(item - m, 2);
      }
    }
    return U(mu::sqrt(sum / (N * M)));
  }

  /********************************* I/O ***********************************/

  template <std::size_t Nn, std::size_t Mm, class U>
  friend std::ostream &operator<<(std::ostream &os, const Matrix<Nn, Mm, U> &m);

  /*************************** matrix <> matrix ****************************/

  template <typename U = T>
  bool operator==(const Matrix<N, M, U> &rhs) const {
    for (std::size_t i = 0; i < N; i++) {
      /* forward comparison to Vector class */
      if (data_[i] != rhs[i]) {
        return false;
      }
    }
    return true;
  }

  template <typename U = T>
  bool operator!=(const Matrix<N, M, U> &rhs) const {
    return !operator==(rhs);
  }

  template <typename U = T>
  Matrix<N, M, T> &operator+=(const Matrix<N, M, U> &rhs) {
    for (std::size_t i = 0; i < N; i++) {
      data_[i] += rhs[i];
    }
    return *this;
  }

  template <typename U = T>
  Matrix<N, M, T> &operator-=(const Matrix<N, M, U> &rhs) {
    for (std::size_t i = 0; i < N; i++) {
      data_[i] -= rhs[i];
    }
    return *this;
  }

  template <typename U = T>
  Matrix<N, M, T> &operator*=(const Matrix<N, M, U> &rhs) {
    for (std::size_t i = 0; i < N; i++) {
      data_[i] *= rhs[i];
    }
    return *this;
  }

  template <typename U = T>
  Matrix<N, M, T> &operator/=(const Matrix<N, M, U> &rhs) {
    for (std::size_t i = 0; i < N; i++) {
      data_[i] /= rhs[i];
    }
    return *this;
  }

  /*************************** matrix <> scalar ****************************/

  /*
   * see Vector for more information
   */

  template <class TScalar>
  typename std::enable_if_t<std::is_arithmetic<TScalar>::value,
                            Matrix<N, M, T> &>
  operator+=(const TScalar &scalar) {
    for (auto &row : data_) {
      row += scalar;
    }
    return *this;
  }

  template <class TScalar>
  typename std::enable_if_t<std::is_arithmetic<TScalar>::value,
                            Matrix<N, M, T> &>
  operator-=(const TScalar &scalar) {
    for (auto &row : data_) {
      row -= scalar;
    }
    return *this;
  }

  template <class TScalar>
  typename std::enable_if_t<std::is_arithmetic<TScalar>::value,
                            Matrix<N, M, T> &>
  operator*=(const TScalar &scalar) {
    for (auto &row : data_) {
      row *= scalar;
    }
    return *this;
  }

  template <class TScalar>
  typename std::enable_if_t<std::is_arithmetic<TScalar>::value,
                            Matrix<N, M, T> &>
  operator/=(const TScalar &scalar) {
    /* a division by zero is forwarded to and handled by the Vector class */
    for (auto &row : data_) {
      row /= scalar;
    }
    return *this;
  }

  /*************************************************************************/

 protected:
  std::array<Vector<M, T>, N> data_;
};

/********************************** I/O ************************************/

template <std::size_t Nn, std::size_t Mm, class U>
std::ostream &operator<<(std::ostream &os, const Matrix<Nn, Mm, U> &m) {
  os << "[ ";
  for (std::size_t i = 0; i < Nn; i++) {
    os << m.data_[i];
    if (i < (Nn - 1)) {
      os << ",\n  ";
    } else {
      os << " ";
    }
  }
  os << "]";
  return os;
}

/**************************** matrix <> matrix *****************************/

template <std::size_t N, std::size_t M, class T, class U = T>
inline Matrix<N, M, T> operator+(const Matrix<N, M, T> &lhs,
                                 const Matrix<N, M, U> &rhs) {
  return Matrix<N, M, T>(lhs) += rhs;
}

template <std::size_t N, std::size_t M, class T, class U = T>
inline Matrix<N, M, T> operator-(const Matrix<N, M, T> &lhs,
                                 const Matrix<N, M, U> &rhs) {
  return Matrix<N, M, T>(lhs) -= rhs;
}

template <std::size_t N, std::size_t M, class T, class U = T>
inline Matrix<N, M, T> operator*(const Matrix<N, M, T> &lhs,
                                 const Matrix<N, M, U> &rhs) {
  return Matrix<N, M, T>(lhs) *= rhs;
}

template <std::size_t N, std::size_t M, class T, class U = T>
inline Matrix<N, M, T> operator/(const Matrix<N, M, T> &lhs,
                                 const Matrix<N, M, U> &rhs) {
  return Matrix<N, M, T>(lhs) /= rhs;
}

/**************************** matrix <> scalar *****************************/

template <std::size_t N, std::size_t M, class T, class TScalar>
typename std::enable_if_t<std::is_arithmetic<TScalar>::value,
                          Matrix<N, M, T>> inline
operator+(const Matrix<N, M, T> &lhs, const TScalar &rhs) {
  return Matrix<N, M, T>(lhs) += rhs;
}

template <std::size_t N, std::size_t M, class T, class TScalar>
typename std::enable_if_t<std::is_arithmetic<TScalar>::value,
                          Matrix<N, M, T>> inline
operator+(const TScalar &lhs, const Matrix<N, M, T> &rhs) {
  return Matrix<N, M, T>(rhs) += lhs;
}

template <std::size_t N, std::size_t M, class T, class TScalar>
typename std::enable_if_t<std::is_arithmetic<TScalar>::value,
                          Matrix<N, M, T>> inline
operator-(const Matrix<N, M, T> &lhs, const TScalar &rhs) {
  return Matrix<N, M, T>(lhs) -= rhs;
}

template <std::size_t N, std::size_t M, class T, class TScalar>
typename std::enable_if_t<std::is_arithmetic<TScalar>::value,
                          Matrix<N, M, T>> inline
operator*(const Matrix<N, M, T> &lhs, const TScalar &rhs) {
  return Matrix<N, M, T>(lhs) *= rhs;
}

template <std::size_t N, std::size_t M, class T, class TScalar>
typename std::enable_if_t<std::is_arithmetic<TScalar>::value,
                          Matrix<N, M, T>> inline
operator*(const TScalar &lhs, const Matrix<N, M, T> &rhs) {
  return Matrix<N, M, T>(rhs) *= lhs;
}

template <std::size_t N, std::size_t M, class T, class TScalar>
typename std::enable_if_t<std::is_arithmetic<TScalar>::value,
                          Matrix<N, M, T>> inline
operator/(const Matrix<N, M, T> &lhs, const TScalar &rhs) {
  return Matrix<N, M, T>(lhs) /= rhs;
}

/************************* convenience functions ***************************/

template <std::size_t N, std::size_t M, class T>
inline T min(const Matrix<N, M, T> &m) {
  return m.min();
}

template <std::size_t N, std::size_t M, class T>
inline T max(const Matrix<N, M, T> &m) {
  return m.max();
}

template <std::size_t N, std::size_t M, class T>
inline T sum(const Matrix<N, M, T> &m) {
  return m.sum();
}

template <class U = void, std::size_t N, std::size_t M, typename T>
inline std::conditional_t<std::is_same<U, void>::value, T, U> mean(
    const Matrix<N, M, T> &m) {
  return m
      .template mean<std::conditional_t<std::is_same<U, void>::value, T, U>>();
}

template <std::size_t N, std::size_t M, typename T>
    std::conditional_t <
    N<M, Vector<N, T>, Vector<M, T>> diag(const Matrix<N, M, T> &m) {
  return m.diag();
}

template <std::size_t N, typename T>
inline Matrix<N, N, T> diag(const Vector<N, T> &v) {
  Matrix<N, N, T> ret{};
  for (std::size_t i = 0; i < N; i++) {
    ret[i][i] = v[i];
  }
  return ret;
}

template <std::size_t N, std::size_t M, typename T>
T det(const Matrix<N, M, T> &m) {
  return m.det();
}

template <std::size_t N, std::size_t M, typename T>
// NOLINTNEXTLINE(runtime/references) intentional non-const reference
inline void transpose(Matrix<N, M, T> &m) {
  m.transpose();
}

template <std::size_t N, std::size_t M, typename T>
Matrix<M, N, T> transposed(const Matrix<N, M, T> &m) {
  return m.transposed();
}

template <typename U = void, std::size_t N1, std::size_t M1, typename T1,
          std::size_t N2, std::size_t M2, typename T2>
std::conditional_t<std::is_same<U, void>::value, Matrix<N1, M2, T1>,
                   Matrix<N1, M2, U>>
dot(const Matrix<N1, M1, T1> &lhs, const Matrix<N2, M2, T2> &rhs) {
  return lhs.template dot<U>(rhs);
}

template <typename U = void, std::size_t N, std::size_t M, typename T,
          std::size_t N2, typename T2>
std::conditional_t<std::is_same<U, void>::value, Vector<N, T>, Vector<N, U>>
dot(const Matrix<N, M, T> &lhs, const Vector<N2, T2> &rhs) {
  return lhs.template dot<U>(rhs);
}

template <std::size_t S, typename T = int>
inline Matrix<S, S, T> eye() {
  Matrix<S, S, T> ret{};
  for (std::size_t i = 0; i < S; i++) {
    ret[i][i] = T{1};
  }
  return ret;
}

template <std::size_t N, std::size_t M, typename T = int>
inline Matrix<N, M, T> ones() {
  return Matrix<N, M, T>{T{1}};
}

template <std::size_t N, std::size_t M, typename T = int>
inline Matrix<N, M, T> zeros() {
  return Matrix<N, M, T>{T{0}};
}

}  // namespace mu

namespace mu {

/* alias templates for Matrix */

template <std::size_t N, typename T>
using MatrixNxN = Matrix<N, N, T>;

template <typename T>
using Matrix2x2 = Matrix<2, 2, T>;

template <typename T>
using Matrix3x3 = Matrix<3, 3, T>;

}  // namespace mu
#endif  // MU_MATRIX_H_