Matrix.cpp

  Calculate the derminant of the matrix
  @return Determinant of matrix.
*/
{
  if (nx != ny)
    throw Kernel::Exception::MisMatch<size_t>(
        nx, ny, "Determinant error :: Matrix is not NxN");

  Matrix<T> Mt(*this); // temp copy
  T D = Mt.factor();
  return D;
}

template <typename T>
T Matrix<T>::factor()
/**
   Gauss jordan diagonal factorisation
   The diagonal is left as the values,
   the lower part is zero.
   @return the factored matrix
*/
{
  if (nx != ny || nx < 1)
    throw std::runtime_error("Matrix::factor Matrix is not NxN");

  double deter = 1.0;
  for (int i = 0; i < static_cast<int>(nx) - 1; i++) // loop over each row
  {
    int jmax = i;
    double Pmax = fabs(V[i][i]);
    for (int j = i + 1; j < static_cast<int>(nx); j++) // find max in Row i
    {
      if (fabs(V[i][j]) > Pmax) {
        Pmax = fabs(V[i][j]);
        jmax = j;
      }
    }
    if (Pmax < 1e-8) // maxtrix signular
    {
      //          std::cerr<<"Matrix Singular"<<std::endl;
      return 0;
    }
    // Swap Columns
    if (i != jmax) {
      swapCols(i, jmax);
      deter *= -1; // change sign.
    }
    // zero all rows below diagonal
    Pmax = V[i][i];
    deter *= Pmax;
    for (int k = i + 1; k < static_cast<int>(nx); k++) // row index
    {
      const double scale = V[k][i] / Pmax;
      V[k][i] = static_cast<T>(0);
      for (int q = i + 1; q < static_cast<int>(nx); q++) // column index
        V[k][q] -= static_cast<T>(scale * V[i][q]);
    }
  }
  deter *= V[nx - 1][nx - 1];
  return static_cast<T>(deter);
}

template <typename T>
void Matrix<T>::normVert()
/**
  Normalise EigenVectors
  Assumes that they have already been calculated
*/
{
  for (size_t i = 0; i < nx; i++) {
    T sum = 0;
    for (size_t j = 0; j < ny; j++) {
      sum += V[i][j] * V[i][j];
    }
    sum = static_cast<T>(std::sqrt(static_cast<double>(sum)));
    for (size_t j = 0; j < ny; j++) {
      V[i][j] /= sum;
    }
  }
  return;
}

template <typename T>
T Matrix<T>::compSum() const
/**
  Add up each component sums for the matrix
  @return \f$ \sum_i \sum_j V_{ij}^2 \f$
 */
{
  T sum(0);
  for (size_t i = 0; i < nx; i++) {
    for (size_t j = 0; j < ny; j++) {
      sum += V[i][j] * V[i][j];
    }
  }
  return sum;
}

template <typename T>
void Matrix<T>::lubcmp(int *rowperm, int &interchange)
/**
  Find biggest pivot and move to top row. Then
  divide by pivot.
  @param interchange :: odd/even nterchange (+/-1)
  @param rowperm :: row permutations [nx values]
*/
{
  double sum, dum, big, temp;

  if (nx != ny || nx < 2) {
    std::cerr << "Error with lubcmp" << std::endl;
    return;
  }
  double *vv = new double[nx];
  interchange = 1;
  for (int i = 0; i < static_cast<int>(nx); i++) {
    big = 0.0;
    for (int j = 0; j < static_cast<int>(nx); j++)
      if ((temp = fabs(V[i][j])) > big)
        big = temp;

    if (big == 0.0) {
      delete[] vv;
      for (int j = 0; j < static_cast<int>(nx); j++) {
        rowperm[j] = j;
      }
      return;
    }
    vv[i] = 1.0 / big;
  }

  for (int j = 0; j < static_cast<int>(nx); j++) {
    for (int i = 0; i < j; i++) {
      sum = V[i][j];
      for (int k = 0; k < i; k++)
        sum -= V[i][k] * V[k][j];
      V[i][j] = static_cast<T>(sum);
    }
    big = 0.0;
    int imax = j;
    for (int i = j; i < static_cast<int>(nx); i++) {
      sum = V[i][j];
      for (int k = 0; k < j; k++)
        sum -= V[i][k] * V[k][j];
      V[i][j] = static_cast<T>(sum);
      if ((dum = vv[i] * fabs(sum)) >= big) {
        big = dum;
        imax = i;
      }
    }

    if (j != imax) {
      for (int k = 0; k < static_cast<int>(nx); k++) { // Interchange rows
        dum = V[imax][k];
        V[imax][k] = V[j][k];
        V[j][k] = static_cast<T>(dum);
      }
      interchange *= -1;
      vv[imax] = static_cast<T>(vv[j]);
    }
    rowperm[j] = imax;

    if (V[j][j] == 0.0)
      V[j][j] = static_cast<T>(1e-14);
    if (j != static_cast<int>(nx) - 1) {
      dum = 1.0 / (V[j][j]);
      for (int i = j + 1; i < static_cast<int>(nx); i++)
        V[i][j] *= static_cast<T>(dum);
    }
  }
  delete[] vv;
  return;
}

template <typename T>
void Matrix<T>::lubksb(const int *rowperm, double *b)
/**
  Impliments a separation of the Matrix
  into a triangluar matrix
*/
{
  int ii = -1;

  for (int i = 0; i < static_cast<int>(nx); i++) {
    int ip = rowperm[i];
    double sum = b[ip];
    b[ip] = b[i];
    if (ii != -1)
      for (int j = ii; j < i; j++)
        sum -= V[i][j] * b[j];
    else if (sum != 0.)
      ii = i;
    b[i] = sum;
  }

  for (int i = static_cast<int>(nx) - 1; i >= 0; i--) {
    double sum = static_cast<T>(b[i]);
    for (int j = i + 1; j < static_cast<int>(nx); j++)
      sum -= V[i][j] * b[j];
    b[i] = sum / V[i][i];
  }
  return;
}

template <typename T>
void Matrix<T>::averSymmetric()
/**
  Simple function to create an average symmetric matrix
  out of the Matrix
*/
{
  const size_t minSize = (nx > ny) ? ny : nx;
  for (size_t i = 0; i < minSize; i++) {
    for (size_t j = i + 1; j < minSize; j++) {
      V[i][j] = (V[i][j] + V[j][i]) / 2;
      V[j][i] = V[i][j];
    }
  }
}

template <typename T>
std::vector<T> Matrix<T>::Diagonal() const
/**
  Returns the diagonal form as a vector
  @return Diagonal elements
*/
{
  const size_t Msize = (ny > nx) ? nx : ny;
  std::vector<T> Diag(Msize);
  for (size_t i = 0; i < Msize; i++) {
    Diag[i] = V[i][i];
  }
  return Diag;
}

template <typename T>
T Matrix<T>::Trace() const
/**
  Calculates the trace of the matrix
  @return Trace of matrix
*/
{
  const size_t Msize = (ny > nx) ? nx : ny;
  T Trx = 0;
  for (size_t i = 0; i < Msize; i++) {
    Trx += V[i][i];
  }
  return Trx;
}

template <typename T>
void Matrix<T>::sortEigen(Matrix<T> &DiagMatrix)
/**
  Sorts the eigenvalues into increasing
  size. Moves the EigenVectors correspondingly
  @param DiagMatrix :: matrix of the EigenValues
*/
{
  if (ny != nx || nx != DiagMatrix.nx || nx != DiagMatrix.ny) {
    std::cerr << "Matrix not Eigen Form" << std::endl;
    throw(std::invalid_argument(" Matrix is not in an eigenvalue format"));
  }
  std::vector<int> index;
  std::vector<T> X = DiagMatrix.Diagonal();
  indexSort(X, index);
  Matrix<T> EigenVec(*this);
  for (size_t Icol = 0; Icol < nx; Icol++) {
    for (size_t j = 0; j < nx; j++) {
      V[j][Icol] = EigenVec[j][index[Icol]];
    }
    DiagMatrix[Icol][Icol] = X[index[Icol]];
  }

  return;
}

template <typename T>
int Matrix<T>::Diagonalise(Matrix<T> &EigenVec, Matrix<T> &DiagMatrix) const
/**
  Attempt to diagonalise the matrix IF symmetric
  @param EigenVec :: (output) the Eigenvectors matrix
  @param DiagMatrix :: the diagonal matrix of eigenvalues
  @return :: 1  on success 0 on failure
*/
{
  if (nx != ny || nx < 1) {
    std::cerr << "Matrix not square" << std::endl;
    return 0;
  }
  for (size_t i = 0; i < nx; i++)
    for (size_t j = i + 1; j < nx; j++)
      if (fabs(V[i][j] - V[j][i]) > 1e-6) {
        std::cerr << "Matrix not symmetric" << std::endl;
        std::cerr << (*this);
        return 0;
      }

  Matrix<T> A(*this);
  // Make V an identity matrix
  EigenVec.setMem(nx, nx);
  EigenVec.identityMatrix();
  DiagMatrix.setMem(nx, nx);
  DiagMatrix.zeroMatrix();

  std::vector<double> Diag(nx);
  std::vector<double> B(nx);
  std::vector<double> ZeroComp(nx);
  // set b and d to the diagonal elements o A
  for (size_t i = 0; i < nx; i++) {
    Diag[i] = B[i] = A.V[i][i];
    ZeroComp[i] = 0;
  }

  int iteration = 0;
  for (int i = 0; i < 100; i++) // max 50 iterations
  {
    double sm = 0.0; // sum of off-diagonal terms
    for (size_t ip = 0; ip < nx - 1; ip++)
      for (size_t iq = ip + 1; iq < nx; iq++)
        sm += fabs(A.V[ip][iq]);

    if (sm == 0.0) // Nothing to do return...
    {
      // Make OUTPUT -- D + A
      // sort Output::
      for (size_t ix = 0; ix < nx; ix++)
        DiagMatrix.V[ix][ix] = static_cast<T>(Diag[ix]);
      return 1;
    }

    // Threshold large for first 5 sweeps
    double tresh = (i < 6) ? 0.2 * sm / static_cast<int>(nx * nx) : 0.0;

    for (int ip = 0; ip < static_cast<int>(nx) - 1; ip++) {
      for (int iq = ip + 1; iq < static_cast<int>(nx); iq++) {
        double g = 100.0 * fabs(A.V[ip][iq]);
        // After 4 sweeps skip if off diagonal small
        if (i > 6 &&
            static_cast<float>(fabs(Diag[ip] + g)) ==
                static_cast<float>(fabs(Diag[ip])) &&
            static_cast<float>(fabs(Diag[iq] + g)) ==
                static_cast<float>(fabs(Diag[iq])))
          A.V[ip][iq] = 0;

        else if (fabs(A.V[ip][iq]) > tresh) {
          double tanAngle, cosAngle, sinAngle;
          double h = Diag[iq] - Diag[ip];
          if (static_cast<float>((fabs(h) + g)) == static_cast<float>(fabs(h)))
            tanAngle = A.V[ip][iq] / h; // tanAngle=1/(2theta)
          else {
            double theta = 0.5 * h / A.V[ip][iq];
            tanAngle = 1.0 / (fabs(theta) + sqrt(1.0 + theta * theta));
            if (theta < 0.0)
              tanAngle = -tanAngle;
          }
          cosAngle = 1.0 / sqrt(1 + tanAngle * tanAngle);
          sinAngle = tanAngle * cosAngle;
          double tau = sinAngle / (1.0 + cosAngle);
          h = tanAngle * A.V[ip][iq];
          ZeroComp[ip] -= h;
          ZeroComp[iq] += h;
          Diag[ip] -= h;
          Diag[iq] += h;
          A.V[ip][iq] = 0;
          // Rotations 0<j<p
          for (int j = 0; j < ip; j++)
            A.rotate(tau, sinAngle, j, ip, j, iq);
          for (int j = ip + 1; j < iq; j++)
            A.rotate(tau, sinAngle, ip, j, j, iq);
          for (int j = iq + 1; j < static_cast<int>(nx); j++)
            A.rotate(tau, sinAngle, ip, j, iq, j);
          for (int j = 0; j < static_cast<int>(nx); j++)
            EigenVec.rotate(tau, sinAngle, j, ip, j, iq);
          iteration++;
        }
      }
    }
    for (size_t j = 0; j < nx; j++) {
      B[j] += ZeroComp[j];
      Diag[j] = B[j];
      ZeroComp[j] = 0.0;
    }
  }
  std::cerr << "Error :: Iterations are a problem" << std::endl;
  return 0;
}

template <typename T>
bool Matrix<T>::isRotation() const
/** Check if a matrix represents a proper rotation
@ return :: true/false
*/
{
  if (this->nx != this->ny)
    throw(std::invalid_argument("matrix is not square"));
  //  std::cout << "Matrix determinant-1 is " << (this->determinant()-1) <<
  //  std::endl;
  if (fabs(this->determinant() - 1) > 1e-5) {
    return false;
  } else {
    Matrix<T> prod(nx, ny), ident(nx, ny, true);
    prod = this->operator*(this->Tprime());
    //    std::cout << "Matrix * Matrix' = " << std::endl << prod << std::endl;
    return prod.equals(ident, 1e-5);
  }
}

template <typename T>
bool Matrix<T>::isOrthogonal() const
/** Check if a matrix is orthogonal. Same as isRotation, but allows determinant
to be -1
@ return :: true/false
*/
{
  if (this->nx != this->ny)
    throw(std::invalid_argument("matrix is not square"));
  if (fabs(fabs(this->determinant()) - 1.) > 1e-5) {
    return false;
  } else {
    Matrix<T> prod(nx, ny), ident(nx, ny, true);
    prod = this->operator*(this->Tprime());
    return prod.equals(ident, 1e-7);
  }
}

template <typename T>
std::vector<T> Matrix<T>::toRotation()
/**
  Transform the matrix to a rotation matrix, by normalizing each column to 1
  @return :: a vector of scaling factors
  @throw :: std::invalid_argument if the absolute value of the determinant is
  less then 1e-10 or not square matrix
*/
{
  if (this->nx != this->ny)
    throw(std::invalid_argument("matrix is not square"));
  if (fabs(this->determinant()) < 1e-10)
    throw(std::invalid_argument("Determinant is too small"));
  // step 1: orthogonalize the matrix
  for (size_t i = 0; i < this->ny; ++i) {
    double spself = 0.;
    for (size_t j = 0; j < this->nx; ++j)
      spself += (V[j][i] * V[j][i]);
    for (size_t k = i + 1; k < this->ny; ++k) {
      double spother = 0;
      for (size_t j = 0; j < this->nx; ++j)
        spother += (V[j][i] * V[j][k]);
      for (size_t j = 0; j < this->nx; ++j)
        V[j][k] -= static_cast<T>(V[j][i] * spother / spself);
    }
  }
  // step 2: get scales and rescsale the matrix
  std::vector<T> scale(this->nx);
  T currentScale;
  for (size_t i = 0; i < this->ny; ++i) {
    currentScale = T(0.);
    for (size_t j = 0; j < this->nx; ++j)
      currentScale += (V[j][i] * V[j][i]);
    currentScale = static_cast<T>(sqrt(static_cast<double>(currentScale)));
    if (currentScale < 1e-10)
      throw(std::invalid_argument("Scale is too small"));
    scale[i] = currentScale;
  }
  Matrix<T> scalingMatrix(nx, ny), change(nx, ny, true);
  for (size_t i = 0; i < this->ny; ++i)
    scalingMatrix[i][i] = static_cast<T>(1.0 / scale[i]);
  *this = this->operator*(scalingMatrix);
  if (this->determinant() < 0.) {
    scale[0] = -scale[0];
    change[0][0] = (T)(-1);
    *this = this->operator*(change);
  }
  return scale;
}

template <typename T>
void Matrix<T>::print() const
/**
  Simple print out routine
 */
{
  write(std::cout, 10);
  return;
}

/** set matrix elements ito random values  in the range from  rMin to rMax*/
template <typename T>
void Matrix<T>::setRandom(size_t seed, double rMin, double rMax) {
  MersenneTwister rng(seed, rMin, rMax);

  for (size_t i = 0; i < nx; i++) {
    for (size_t j = 0; j < ny; j++) {
      V[i][j] = static_cast<T>(rng.nextValue());
    }
  }
}

template <typename T>
void Matrix<T>::write(std::ostream &Fh, const int blockCnt) const
/**
  Write out function for blocks of 10 Columns
  @param Fh :: file stream for output
  @param blockCnt :: number of columns per line (0 == full)
*/
{
  std::ios::fmtflags oldFlags = Fh.flags();
  Fh.setf(std::ios::floatfield, std::ios::scientific);
  const size_t blockNumber((blockCnt > 0) ? blockCnt : ny);
  size_t BCnt(0);
  do {
    const size_t ACnt = BCnt;
    BCnt += blockNumber;
    if (BCnt > ny) {
      BCnt = ny;
    }

    if (ACnt) {
      Fh << " ----- " << ACnt << " " << BCnt << " ------ " << std::endl;
    }
    for (size_t i = 0; i < nx; i++) {
      for (size_t j = ACnt; j < BCnt; j++) {
        Fh << std::setw(10) << V[i][j] << "  ";
      }
      Fh << std::endl;
    }
  } while (BCnt < ny);

  Fh.flags(oldFlags);
  return;
}

template <typename T>
std::string Matrix<T>::str() const
/**
  Convert the matrix into a simple linear string expression
  @return String value of output
*/
{
  std::ostringstream cx;
  for (size_t i = 0; i < nx; i++) {
    for (size_t j = 0; j < ny; j++) {
      cx << std::setprecision(6) << V[i][j] << " ";
    }
  }
  return cx.str();
}

/**
 * Write an object to a stream. Format will be
 * Matrix(nrows,ncols)x_00,x_01...,x_10,x_11
 * @param os :: output stream
 * @param matrix :: Matrix to write out
 * @return The output stream (of)
*/
template <typename T>
std::ostream &operator<<(std::ostream &os, const Matrix<T> &matrix) {
  dumpToStream(os, matrix, ',');
  return os;
}

/**
 * Write a Matrix to a stream. Format will be
 * Matrix(nrowsSEPncols)x_00SEPx_01...SEPx_10SEPx_11
 * @param os :: output stream
 * @param matrix :: Matrix to write out
 * @param delimiter :: A character to use as delimiter for the string
*/
template <typename T>
void dumpToStream(std::ostream &os, const Kernel::Matrix<T> &matrix,
                  const char delimiter) {
  size_t nrows(matrix.numRows()), ncols(matrix.numCols());
  os << "Matrix(" << nrows << delimiter << ncols << ")";
  for (size_t i = 0; i < nrows; ++i) {
    for (size_t j = 0; j < ncols; ++j) {
      os << matrix[i][j];
      if (i < nrows - 1 || j < ncols - 1)
        os << delimiter;
    }
  }
}

/**
* Fill an object from a stream. Format should be
* Matrix(nrows,ncols)x_00,x_01...,x_10,x_11
* @param is :: A stream object
* @param in :: An object to fill
* @returns A reference to the stream
*/
template <typename T>
std::istream &operator>>(std::istream &is, Kernel::Matrix<T> &in) {
  fillFromStream(is, in, ',');
  return is;
}

/**
* Fill a Matrix from a stream using the given separator. Format should be
* Matrix(nrowsSEPncols)x_00SEPx_01...SEPx_10SEPx_11
* where SEP is replaced by the given separator
* @param is :: A stream object
* @param in :: An Matrix object to fill
* @param delimiter :: A single character separator that delimits the entries
*/
template <typename T>
void fillFromStream(std::istream &is, Kernel::Matrix<T> &in,
                    const char delimiter) {
  // Stream should start with Matrix(
  char dump;
  std::string start(7, ' ');
  for (int i = 0; i < 7; ++i) {
    is >> dump;
    start[i] = dump;
    if (!is)
      throw std::invalid_argument(
          "Unexpected character when reading Matrix from stream.");
  }
  if (start != "Matrix(")
    throw std::invalid_argument("Incorrect input format for Matrix stream.");
  // Now read a nrows,ncols and )
  size_t nrows(0), ncols(0);
  is >> nrows;
  if (!is)
    throw std::invalid_argument("Expected number of rows when reading Matrix "
                                "from stream, found something else.");
  is >> dump;
  is >> ncols;
  if (!is)
    throw std::invalid_argument("Expected number of columns when reading "
                                "Matrix from stream, found something else.");
  is >> dump;
  if (dump != ')')
    throw std::invalid_argument("Expected closing parenthesis after ncols when "
                                "reading Matrix from stream, found something "
                                "else.");

  // Resize the matrix
  in.setMem(nrows, ncols);

  // Use getline with the delimiter set to "," to read
  std::string value_str;
  size_t row(0), col(0);
  while (!is.eof() && std::getline(is, value_str, delimiter)) {
    try {
      T value = boost::lexical_cast<T>(value_str);
      in.V[row][col] = value;
    } catch (boost::bad_lexical_cast &) {
      throw std::invalid_argument(
          "Unexpected type found while reading Matrix from stream: \"" +
          value_str + "\"");
    }
    ++col;
    if (col == ncols) // New row
    {
      col = 0;
      ++row;
    }
  }
}

///\cond TEMPLATE

// Symbol definitions for common types
template class MANTID_KERNEL_DLL Matrix<double>;
template class MANTID_KERNEL_DLL Matrix<int>;
template class MANTID_KERNEL_DLL Matrix<float>;

template MANTID_KERNEL_DLL std::ostream &operator<<(std::ostream &,
                                                    const DblMatrix &);
template MANTID_KERNEL_DLL void dumpToStream(std::ostream &, const DblMatrix &,
                                             const char);
template MANTID_KERNEL_DLL std::istream &operator>>(std::istream &,
                                                    DblMatrix &);
template MANTID_KERNEL_DLL void fillFromStream(std::istream &, DblMatrix &,
                                               const char);

template MANTID_KERNEL_DLL std::ostream &operator<<(std::ostream &,
                                                    const Matrix<float> &);
template MANTID_KERNEL_DLL void dumpToStream(std::ostream &,
                                             const Matrix<float> &, const char);
template MANTID_KERNEL_DLL std::istream &operator>>(std::istream &,
                                                    Matrix<float> &);
template MANTID_KERNEL_DLL void fillFromStream(std::istream &, Matrix<float> &,
                                               const char);

template MANTID_KERNEL_DLL std::ostream &operator<<(std::ostream &,
                                                    const IntMatrix &);
template MANTID_KERNEL_DLL void dumpToStream(std::ostream &, const IntMatrix &,
                                             const char);
template MANTID_KERNEL_DLL std::istream &operator>>(std::istream &,
                                                    IntMatrix &);
template MANTID_KERNEL_DLL void fillFromStream(std::istream &, IntMatrix &,
                                               const char);
///\endcond TEMPLATE

} // namespace Kernel
} // namespace