/***************************************************************************
 *   Copyright (C) 2005 by Joachim Friedrich                               *
 *   joachim.friedrich@biozentrum.uni-wuerzburg.de                         *
 *                                                                         *
 *   This file is part of profdist and cbcanalyzer                         *
 *                                                                         *
 *   Both profdist and cbcanalyzer are free software; you can redistribute * 
 *   it and/or modify it under the terms of the GNU General Public License * 
 *   as published by the Free Software Foundation; either version 2 of the * 
 *   License, or (at your option) any later version.                       *
 *                                                                         *
 *   Profdist and cbcanalyzer are distributed in the hope that it will be  *
 *   useful, but WITHOUT ANY WARRANTY; without even the implied warranty   *
 *   of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

//---------------------------------------------------------------------------
// tMatrix.h
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
// INCLUDES
//---------------------------------------------------------------------------

#include <complex>
#include <string>
#include <vector>
#include <cstring>
#include <map>
#include <iostream>
#include <fstream>
#include <math.h>
#include <iomanip>

#include "functions.h"
#include "types.h"

//---------------------------------------------------------------------------
// DEFINES, MAKROS
//---------------------------------------------------------------------------
#if !defined(tVector_H)
  #include "tVector.h"
#endif

#define USE_EXCEPTIONS

#if defined(USE_EXCEPTIONS)
  #include "exceptions.h"
#endif

#if !defined(tMatrix_H)
  #define tMatrix_H

using namespace std;

typedef std::vector < std::string > vecTypeNames;
typedef std::map < std::vector < int >, bool > mapTypeVecII;
typedef std::map < int, bool > mapTypeIB;

  #define VERY_SMALL	1.0e-12
  #define VERY_LARGE	1.0E+12

// mine

  #define zeroc complex<double>(0,0)

  #define ZERO VERY_SMALL
  #define SHIFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
  #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))

//---------------------------------------------------------------------------
// CLASS DEFINITIONS
//---------------------------------------------------------------------------



template < class T >
class tMatrix
{
public:
  typedef typename std::vector<T>::const_iterator const_iterator;
  typedef typename std::vector<T>::iterator iterator;
  typedef typename std::vector<T>::value_type value_type;

  const_iterator begin() const { return mMat.begin(); }
  iterator begin(){ return mMat.begin(); }
  const_iterator end() const{ return mMat.end(); }
  iterator end(){ return mMat.end(); }
  //-------------------------------------------------------------------------
  // Constructors, Destructors
  //-------------------------------------------------------------------------
  tMatrix() 
	: mRows(0), mCols(0)
  {
  }

  tMatrix( size_t row, size_t col ): mRows(row), mCols(col)
  {
    mMat.resize( row * col );
  }

  tMatrix( size_t row, size_t col, const std::vector< T > & sv ) 
	  : mRows( row ), mCols( col ), mMat( sv )
  {
  }

  tMatrix( size_t row, size_t col, T const& x ) 
	  : mRows( row ), mCols( col ), mMat( row*col, x )
  {
  }

  tMatrix( size_t row, size_t col, T const* x ) 
	  : mRows( row ), mCols( col ), mMat( row*col, *x )
  {
  }


  //-------------------------------------------------------------------------
  // Functions
  //-------------------------------------------------------------------------
  tMatrix < T > & readPhylip( std::istream & infile );
  tMatrix < T > & readPhylip( std::istream & infile, vecTypeNames & vecSeqNames );
  tMatrix < T > & setDiag( const T value );
  void writeFile( std::ostream & outfile ) const;
  void writeFile4( std::ostream & outfile ) const;
  void writeFile( std::ostream & outfile, vecTypeNames const& vecSeqNames ) const;
  void swap( size_t row1, size_t col1, size_t row2, size_t col2 );

  /**
   * @param mapIdentSeq
   * @param iThreshold 
   */
  void getIdentSeq( mapTypeVecII & mapIdentSeq, int iThreshold );
  std::string output() const;
  std::string output( vecTypeNames const& vecSeqNames ) const;
  T getMax() const;
  size_t getMaxLen() const;

  // Implementation
  inline size_t nRows() const {  return mRows;  }

  inline size_t nCols() const {  return mCols; }

  void resize( size_t rows, size_t cols )
  {
    mRows = rows;
    mCols = cols;
    mMat.resize( rows * cols );
  }

  void resize( size_t rows, size_t cols, T val )
  {
    mRows = rows;
    mCols = cols;
    mMat.resize( rows * cols );
    std::fill( mMat.begin(), mMat.end(), val );
  }

  void show()
  {
    long i, j;
    cout << mRows << "  " << mCols << endl;
    for ( i = 0; i < mRows; ++i )
    {
      for ( j = 0; j < mCols; ++j )
      {
        cout << mMat( i * mCols + j ) << "   ";
      }
      cout << endl;
    }
    cout << endl;

  }

  //-------------------------------------------------------------------------
  // Operators
  //-------------------------------------------------------------------------
  tMatrix & operator = ( T const& value ) // assignment
  {
    std::fill( mMat.begin(), mMat.end(), value);
    return * this;
  }

  tMatrix & operator *= ( T const& value ) // assignment
  {
    std::transform( mMat.begin(), mMat.end(), std::bind2nd( std::multiplies<T>(), value )  );
    return * this;
  }


  // indexing operators
  T & operator() ( size_t i ) // 1d indexing
  {
    return mMat[i];
  }

  T const& operator() ( size_t i ) const // 1d indexing
  {
    return mMat[i];
  }


  inline T & operator() ( size_t i, size_t j ) {  return mMat[ i * mCols + j ]; }
  inline T const& operator() ( size_t i, size_t j ) const{ return mMat[ i * mCols + j ]; }

protected:
  // Attributes
  std::size_t mRows; ///< number of rows
  std::size_t mCols; ///< number of columns
  std::vector< T > mMat; ///< base vector container
};

template<typename T> 
std::ostream & operator << ( std::ostream & s, tMatrix < T > const& m )
{
  s << m.nRows() << "  " << m.nCols() << endl;
  for ( int i = 0; i < m.nRows(); ++i )
  {
    for ( int j = 0; j < m.nCols(); ++j )
    {
      s << m( i, j ) << "  ";
    }
    s << endl;
  }
  s << endl;

  return s;
}


/**
 * @brief normalize the matrix, so that the sum of the rows is 1
 * @param n the matrix to normalize
 */
template<typename T>
void normalize_log_det( tMatrix<T> & n );

/**
 * @brief calculates M to the power of t
 * @param M base matrix
 * @param t scalar exponent
 * @param expm aproximation flag
 * You can control the calculation method using the third parameter
 * possible Values are Taylor and Pade, where Taylor is the default.
 */
template<typename T>
tMatrix<T> exp( tMatrix<T> const& M, T const& t, profdist::ExpMType expm = profdist::Taylor);

/**
 * @brief calculates M to the power of t using pade approximation
 * @param M base matrix
 * @param t scalar exponent
 */
template<typename T>
tMatrix<T> exp_pade( tMatrix<T> const& M, T const& t );

/**
 * @brief calculates M to the power of t using taylor approximation
 * @param M base matrix
 * @param t scalar exponent
 */
template<typename T>
tMatrix<T> exp_taylor( tMatrix<T> const& M, T const& t );

/**
 * @brief create transposed copy of parameter
 * @param M matrix 
 */
template<class T>
tMatrix<T> transpose( tMatrix<T> const& M );


/**
 * @brief calculate the determinant of a matrix
 */
template<class T>
T det( tMatrix<T> const& M );






#include "tMatrix.inl"

#endif // tMatrix_H
