EigenMatrix.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidCurveFitting/EigenMatrix.h"
#include <iostream>
 
namespace Mantid::CurveFitting {
 
EigenMatrix::EigenMatrix(const size_t nx, const size_t ny)
    : m_data(nx * ny), m_view(EigenMatrix_View(m_data.data(), nx, ny)) {}
 
EigenMatrix::EigenMatrix(std::initializer_list<std::initializer_list<double>> ilist)
    : EigenMatrix(ilist.size(), ilist.begin()->size()) {
  for (auto row = ilist.begin(); row != ilist.end(); ++row) { // loop through row by row
    if (row->size() != size2()) {                             // check number of cells in row equals number of cols
      throw std::runtime_error("All rows in initializer list must have the same size.");
    }
    auto i = static_cast<size_t>(std::distance(ilist.begin(), row));
    for (auto cell = row->begin(); cell != row->end(); ++cell) {
      auto j = static_cast<size_t>(std::distance(row->begin(), cell));
      set(i, j, *cell);
    }
  }
}
 
EigenMatrix::EigenMatrix(const EigenMatrix &M)
    : m_data(M.m_data), m_view(EigenMatrix_View(m_data.data(), M.size1(), M.size2())) {}
 
EigenMatrix::EigenMatrix(EigenMatrix &M, size_t row, size_t col, size_t nRows, size_t nCols) {
  if (row + nRows > M.size1() || col + nCols > M.size2()) {
    throw std::runtime_error("Submatrix exceeds matrix size.");
  }
  size_t newSize = nRows * nCols;
  m_data.resize(newSize);
  m_view = EigenMatrix_View(M.mutator(), nRows, nCols, row, col);
  std::copy(M.mutator().data(), M.mutator().data() + newSize, m_data.data());
}
 
EigenMatrix::EigenMatrix(const Kernel::Matrix<double> &M)
    : m_data(M.numRows() * M.numCols()), m_view(EigenMatrix_View(m_data.data(), M.numRows(), M.numCols())) {
 
  // Mantid::Kernel::Matrix stores row by row, so have to transpose initial matrix.
  // To do this we copy into an eign matrix, then tranpose.
  auto tempVec = M.getVector();
  Eigen::MatrixXd tempMatrTr = EigenMatrix_View(tempVec.data(), M.numCols(), M.numRows()).matrix_mutator();
  Eigen::MatrixXd tempMatr = tempMatrTr.transpose();
  std::copy(tempMatr.data(), tempMatr.data() + tempMatr.size(), m_data.data());
}
 
EigenMatrix::EigenMatrix(const Kernel::Matrix<double> &M, size_t row, size_t col, size_t nRows, size_t nCols) {
  if (row + nRows > M.numRows() || col + nCols > M.numCols()) {
    throw std::runtime_error("Submatrix exceeds matrix size.");
  }
  resize(nRows, nCols);
 
  // Mantid::Kernel::Matrix stores row by row, so have to transpose initial matrix before taking view.
  // To do this we copy into an eign matrix, tranpose, then take a sub matrix.
  auto tempVec = M.getVector();
  Eigen::MatrixXd tempMatrTr = EigenMatrix_View(tempVec.data(), M.numCols(), M.numRows()).matrix_mutator();
  Eigen::MatrixXd tempMatr = tempMatrTr.transpose();
  auto tempMatrSubView = EigenMatrix_View(tempMatr.data(), tempMatr.rows(), tempMatr.cols(), nRows, nCols, row, col);
  Eigen::MatrixXd tempMatrSub = tempMatrSubView.matrix_inspector();
  std::copy(tempMatrSub.data(), tempMatrSub.data() + tempMatrSub.size(), m_data.data());
}
 
EigenMatrix::EigenMatrix(std::vector<double> &&data, size_t nx, size_t ny)
    : m_data(std::move(data)), m_view(EigenMatrix_View(m_data.data(), nx, ny)) {}
 
EigenMatrix &EigenMatrix::operator=(const EigenMatrix &M) {
  m_data = M.m_data;
  m_view = EigenMatrix_View(m_data.data(), M.size1(), M.size2());
  return *this;
}
 
EigenMatrix &EigenMatrix::operator=(const Eigen::MatrixXd m) {
  m_data.resize(m.rows() * m.cols());
  for (size_t i = 0; i < (size_t)m.size(); i++) {
    m_data[i] = *(m.data() + i);
  }
  m_view = EigenMatrix_View(m_data.data(), m.rows(), m.cols());
  return *this;
}
 
bool EigenMatrix::isEmpty() const { return m_data.empty(); }
 
void EigenMatrix::resize(const size_t nx, const size_t ny) {
  if (nx == 0 && ny == 0) {
    // Matrix as minimum is 1x1, retained from gsl for consistency.
    m_data.resize(2);
    m_view = EigenMatrix_View(m_data.data(), 1, 1);
 
  } else {
    m_data.resize(nx * ny);
    m_view = EigenMatrix_View(m_data.data(), nx, ny);
  }
}
 
size_t EigenMatrix::size1() const { return m_view.rows(); }
 
size_t EigenMatrix::size2() const { return m_view.cols(); }
 
void EigenMatrix::set(size_t i, size_t j, double value) {
  if (isEmpty()) {
    throw std::out_of_range("Matrix is empty.");
  }
  if (i < m_view.rows() && j < m_view.cols())
    m_view.matrix_mutator()(i, j) = value;
  else {
    throw std::out_of_range("EigenMatrix indices are out of range.");
  }
}
 
double EigenMatrix::get(size_t i, size_t j) const {
  if (isEmpty()) {
    throw std::out_of_range("Matrix is empty.");
  }
  if (i < m_view.rows() && j < m_view.cols())
    return m_view.matrix_inspector()(i, j);
  else {
    throw std::out_of_range("EigenMatrix indices are out of range.");
  }
}
 
void EigenMatrix::identity() { m_view.matrix_mutator().setIdentity(); }
 
void EigenMatrix::zero() { m_view.matrix_mutator().setZero(); }
 
void EigenMatrix::diag(const EigenVector &d) {
  const auto n = d.size();
  resize(n, n);
  zero();
  for (size_t i = 0; i < n; ++i) {
    set(i, i, d.get(i));
  }
}
 
EigenMatrix &EigenMatrix::operator+=(const EigenMatrix &M) {
  m_view.matrix_mutator() += M.inspector();
  return *this;
}
EigenMatrix &EigenMatrix::operator+=(const double &d) {
  m_view.matrix_mutator().array() += d;
  return *this;
}
EigenMatrix &EigenMatrix::operator-=(const EigenMatrix &M) {
  m_view.matrix_mutator() -= M.inspector();
  return *this;
}
EigenMatrix &EigenMatrix::operator-=(const double &d) {
  m_view.matrix_mutator().array() -= d;
  return *this;
}
EigenMatrix &EigenMatrix::operator*=(const double &d) {
  m_view.matrix_mutator() *= d;
  return *this;
}
 
EigenVector EigenMatrix::operator*(const EigenVector &v) const {
  if (v.size() != size2()) {
    throw std::invalid_argument("Matrix by vector multiplication: wrong size of vector.");
  }
 
  EigenVector res;
  res = inspector() * v.inspector();
  return res;
}
 
EigenMatrix EigenMatrix::operator*(const EigenMatrix &m) const {
  if (m.size1() != size2()) {
    throw std::invalid_argument("Matrix by matrix multiplication: matricies are of incompatible sizes.");
  }
 
  EigenMatrix res;
  res = inspector() * m.inspector();
  return res;
}
 
void EigenMatrix::solve(const EigenVector &rhs, EigenVector &x) {
  if (size1() != size2()) {
    throw std::invalid_argument("System of linear equations: the matrix must be square.");
  }
  size_t n = size1();
  if (rhs.size() != n) {
    throw std::invalid_argument("System of linear equations: right-hand side vector has wrong size.");
  }
  if (det() == 0) {
    throw std::invalid_argument("Matrix A is singular.");
  }
 
  Eigen::MatrixXd b = rhs.inspector();
  Eigen::ColPivHouseholderQR<Eigen::MatrixXd> dec(inspector());
  auto res = dec.solve(b);
  x = res;
 
  // Commented out soltution error, as gsl implementation didn't provide this.
  // if (!rhs.inspector().isApprox(inspector() * x.inspector())) {
  //  //throw std::runtime_error("Matrix Solution Error: solution does not exist.");
  //}
}
 
void EigenMatrix::invert() {
  if (size1() != size2()) {
    throw std::runtime_error("Matrix inverse: the matrix must be square.");
  }
  mutator() = inspector().inverse();
}
 
double EigenMatrix::det() const {
  if (size1() != size2()) {
    throw std::runtime_error("Matrix inverse: the matrix must be square.");
  }
 
  return inspector().determinant();
}
 
void EigenMatrix::eigenSystem(EigenVector &eigenValues, EigenMatrix &eigenVectors) {
  size_t n = size1();
  if (n != size2()) {
    throw std::runtime_error("Matrix eigenSystem: the matrix must be square.");
  }
 
  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> solver;
  solver.compute(inspector());
 
  // previously gsl used "gsl_eigen_symmv" to calculate the eigenSystem. This function only returned
  // real eigenvalues/vectors. The eigen function used can handle and return complex eigen vectors.
  // this causes errors when returning to a non-complex EigenVector and EigenMatrix classes.
  // this check is here to alert to this error.
  if (!solver.eigenvectors().imag().isZero()) {
    throw std::invalid_argument("eigensystem has complex eigenvalues or eigenvectors");
  } else {
    eigenValues = solver.eigenvalues().real();
    eigenVectors = solver.eigenvectors().real();
  }
}
 
EigenVector EigenMatrix::copyRow(size_t i) const {
  if (i >= size1()) {
    throw std::out_of_range("EigenMatrix row index is out of range.");
  }
 
  EigenMatrix_View rowView = EigenMatrix_View(copy_view().data(), size1(), size2(), 1, size1(), i, 0);
  Eigen::VectorXd rowVec = Eigen::Map<Eigen::VectorXd, 0, Eigen::InnerStride<>>(
      rowView.matrix_mutator().data(), rowView.cols(), Eigen::InnerStride<>(rowView.outerStride()));
  return EigenVector(&rowVec);
}
 
EigenVector EigenMatrix::copyColumn(size_t i) const {
  if (i >= size2()) {
    throw std::out_of_range("EigenMatrix column index is out of range.");
  }
 
  EigenMatrix_View colView = EigenMatrix_View(copy_view().data(), size1(), size2(), size1(), 1, 0, i);
  Eigen::VectorXd colVec = Eigen::Map<Eigen::VectorXd, 0, Eigen::InnerStride<>>(
      colView.matrix_mutator().data(), colView.rows(), Eigen::InnerStride<>(colView.innerStride()));
  return EigenVector(&colVec);
}
 
EigenMatrix EigenMatrix::move() { return EigenMatrix(std::move(m_data), size1(), size2()); }
 
EigenVector EigenMatrix::multiplyByVector(const EigenVector &v) const {
  if (v.size() != size2()) {
    throw std::invalid_argument("Matrix by vector multiplication: wrong size of vector.");
  }
 
  EigenVector res;
  res = inspector() * v.inspector();
  return res;
}
 
EigenMatrix EigenMatrix::tr() const {
  EigenMatrix res;
  res = inspector().transpose();
  return res;
}
 
} // namespace Mantid::CurveFitting