Matrix.cpp

#include "MantidKernel/Matrix.h"
#include "MantidKernel/V3D.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/MersenneTwister.h"
#include <iomanip>
#include <iostream>

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>());
      return;
    }
    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>&);
  }



template<typename T>
std::vector<T> 
Matrix<T>::getVector()const
{
  std::vector<T> rez(nx*ny);
  size_t ic(0);
  for(size_t i=0;i<nx;i++)
  {
    for(size_t j=0;j<ny;j++)
    {
      rez[ic]=V[i][j];
      ic++;
    }
  }
  return rez;
}
//
template<typename T>
Matrix<T>::Matrix(const size_t nrow, const size_t ncol, const bool makeIdentity)
  : nx(0),ny(0),V(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:: nx,ny zeroed so setMem always works
  setMem(nrow,ncol);
  zeroMatrix();
  if(makeIdentity==true) identityMatrix();
}

template<typename T>
Matrix<T>::Matrix(const std::vector<T>& A,const std::vector<T>& B)
  : nx(0),ny(0),V(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:: nx,ny zeroed so setMem always works
  setMem(A.size(),B.size());
  for(size_t i=0;i<nx;i++)
  {
    for(size_t j=0;j<ny;j++)
    {
      V[i][j]=A[i]*B[j];
    }
  }  
}
//
template<typename T>
Matrix<T>::Matrix(const std::vector<T>&data) : nx(0),ny(0),V(0)
{
        size_t numel = data.size();
        size_t nxt =(size_t)sqrt(double(numel));
        size_t test = nxt*nxt;
        if(test != numel){
                throw(std::invalid_argument("number of elements in input vector have to be square of some value"));
        }

        setMem(nxt,nxt);

        size_t ic(0);
        for(size_t i=0;i<nx;i++)
        {
                for(size_t j=0;j<ny;j++)
                {
                        V[i][j]=data[ic];
                        ic++;
                }
        }  
}

template<typename T>
Matrix<T>::Matrix(const Matrix<T>& A,const size_t nrow,const size_t ncol)
  : nx(A.nx-1),ny(A.ny-1),V(0)
  /**
    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
  */
{
  // Note:: nx,ny zeroed so setMem always works
  if (nrow>nx)
    throw Kernel::Exception::IndexError(nrow,A.nx, "Matrix::Constructor without col");
  if (ncol>ny)
    throw Kernel::Exception::IndexError(ncol,A.ny, "Matrix::Constructor without col");
  setMem(nx,ny);
  size_t iR(0);
  for(size_t i=0;i<=nx;i++)
  {
    if (i!=nrow)
    {
      size_t jR(0);
      for(size_t j=0;j<=ny;j++)
      {
        if (j!=ncol)
        {
          V[iR][jR]=A.V[i][j];
          jR++;
        }
      }
          iR++;
    }
  }
}


template<typename T>
Matrix<T>::Matrix(const Matrix<T>& A)
  : nx(0),ny(0),V(0)