Matrix.cpp

#include "MantidKernel/Matrix.h"

#include "MantidKernel/Exception.h"
#include "MantidKernel/make_unique.h"
#include "MantidKernel/MersenneTwister.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/V3D.h"

#include <algorithm>
#include <memory>
#include <sstream>

using Mantid::Kernel::TimeSeriesProperty;

namespace Mantid {

namespace Kernel {
//
#define fabs(x) std::fabs((x)*1.0)

namespace {
//-------------------------------------------------------------------------
// Utility methods and function objects in anonymous namespace
//-------------------------------------------------------------------------
/**
\class PIndex
\author S. Ansell
\date Aug 2005
\version 1.0
\brief Class  to fill an index with a progressive count
*/
template <typename T> struct PIndex {
private:
  int count; ///< counter
public:
  /// Constructor
  PIndex() : count(0) {}
  /// functional
  std::pair<T, int> operator()(const T &A) {
    return std::pair<T, int>(A, count++);
  }
};

/**
\class PSep
\author S. Ansell
\date Aug 2005
\version 1.0
\brief Class to access the second object in index pair.
*/
template <typename T> struct PSep {
  /// Functional to the second object
  int operator()(const std::pair<T, int> &A) { return A.second; }
};

/**
* Function to take a vector and sort the vector
* so as to produce an index. Leaves the vector unchanged.
* @param pVec :: Input vector
* @param Index :: Output vector
*/
template <typename T>
void indexSort(const std::vector<T> &pVec, std::vector<int> &Index) {
  Index.resize(pVec.size());
  std::vector<typename std::pair<T, int>> PartList;
  PartList.resize(pVec.size());

  transform(pVec.begin(), pVec.end(), PartList.begin(), PIndex<T>());
  sort(PartList.begin(), PartList.end());
  transform(PartList.begin(), PartList.end(), Index.begin(), PSep<T>());
}
template void indexSort(const std::vector<double> &, std::vector<int> &);
template void indexSort(const std::vector<float> &, std::vector<int> &);
template void indexSort(const std::vector<int> &, std::vector<int> &);
} // End of Anonymous name space

template <typename T> std::vector<T> Matrix<T>::getVector() const {
  std::vector<T> rez(m_numRows * m_numColumns);
  size_t ic(0);
  for (size_t i = 0; i < m_numRows; i++) {
    for (size_t j = 0; j < m_numColumns; j++) {
      rez[ic] = m_rawData[i][j];
      ic++;
    }
  }
  return rez;
}
//
template <typename T>
Matrix<T>::Matrix(const size_t nrow, const size_t ncol, const bool makeIdentity)
    : m_numRows(0), m_numColumns(0)
/**
  Constructor with pre-set sizes. Matrix is zeroed
  @param nrow :: number of rows
  @param ncol :: number of columns
  @param makeIdentity :: flag for the constructor to return an identity matrix
*/
{
  // Note:: m_numRows, m_numColumns zeroed so setMem always works
  setMem(nrow, ncol);
  zeroMatrix();
  if (makeIdentity)
    identityMatrix();
}

template <typename T>
Matrix<T>::Matrix(const std::vector<T> &A, const std::vector<T> &B)
    : m_numRows(0), m_numColumns(0)
/**
  Constructor to take two vectors and multiply them to
  construct a matrix. (assuming that we have columns x row
  vector.
  @param A :: Column vector to multiply
  @param B :: Row vector to multiply
*/
{
  // Note:: m_numRows,m_numColumns zeroed so setMem always works
  setMem(A.size(), B.size());
  for (size_t i = 0; i < m_numRows; i++) {
    for (size_t j = 0; j < m_numColumns; j++) {
      m_rawData[i][j] = A[i] * B[j];
    }
  }
}
//
template <typename T>
Matrix<T>::Matrix(const std::vector<T> &data)
    : m_numRows(0), m_numColumns(0) {
  size_t numElements = data.size();
  size_t numRows = static_cast<size_t>(sqrt(double(numElements)));
  size_t numRowsSquare = numRows * numRows;
  if (numElements != numRowsSquare) {
    throw(std::invalid_argument(
        "number of elements in input vector have to be square of some value"));
  }

  setMem(numRows, numRows);

  size_t ic(0);
  for (size_t i = 0; i < m_numRows; i++) {
    for (size_t j = 0; j < m_numColumns; j++) {
      m_rawData[i][j] = data[ic];
      ic++;
    }
  }
}

template <typename T>
Matrix<T>::Matrix(const std::vector<T> &data, const size_t nrow,
                  const size_t ncol)
    : m_numRows(0), m_numColumns(0) {
  size_t numel = data.size();
  size_t test = nrow * ncol;
  if (test != numel) {
    throw(std::invalid_argument("number of elements in input vector have is "
                                "incompatible with the number of rows and "
                                "columns"));
  }

  setMem(nrow, ncol);

  size_t ic(0);
  for (size_t i = 0; i < m_numRows; i++) {
    for (size_t j = 0; j < m_numColumns; j++) {
      m_rawData[i][j] = data[ic];
      ic++;
    }
  }
}

template <typename T>
Matrix<T>::Matrix(const Matrix<T> &A, const size_t nrow, const size_t ncol)
    : m_numRows(A.m_numRows - 1), m_numColumns(A.m_numColumns - 1)
/**
  Constructor with for a missing row/column.
  @param A :: The input matrix
  @param nrow :: number of row to miss
  @param ncol :: number of column to miss
*/
{
  if (nrow > m_numRows)
    throw Kernel::Exception::IndexError(nrow, A.m_numRows,
                                        "Matrix::Constructor without row");
  if (ncol > m_numColumns)
    throw Kernel::Exception::IndexError(ncol, A.m_numColumns,
                                        "Matrix::Constructor without col");
  setMem(m_numRows, m_numColumns);
  size_t iR(0);
  for (size_t i = 0; i <= m_numRows; i++) {
    if (i != nrow) {
      size_t jR(0);
      for (size_t j = 0; j <= m_numColumns; j++) {
        if (j != ncol) {
          m_rawData[iR][jR] = A.m_rawData[i][j];
          jR++;
        }
      }
      iR++;
    }
  }