/***************************************************************************
 *   Copyright (C) 2005 by Andreas Pokorny                                 *
 *   andreas.pokorny@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.             *
 ***************************************************************************/

#ifndef PROFDIST_LIB_FIXED_MATRIX_H_INCLUDED
#define PROFDIST_LIB_FIXED_MATRIX_H_INCLUDED

#include <algorithm>
#include <stdexcept>
#include <iostream>
namespace profdist {

template<typename L,typename R> struct super_type;
template<typename B>  struct super_type<B,B>{ typedef B type; };
template<>  struct super_type<double,double>{ typedef double type; };
template<typename R>  struct super_type<double,R>{ typedef double type; };
template<typename L>  struct super_type<L,double>{ typedef double type; };
#define ST(A,B,C) template<>  struct super_type<A,B>{ typedef C type; }; \
template <> struct super_type<B,A>{ typedef C type; }
// float rules
ST(float,long long,float );
ST(float,long,float);
ST(float,int,float);
ST(float,short,float);
ST(float,char,float);
ST(float,unsigned long long,float );
ST(float,unsigned long,float);
ST(float,unsigned int,float);
ST(float,unsigned short,float);
ST(float,unsigned char,float);

// integer type rules
ST(long long,long,long long);
ST(long long,int,long long);
ST(long long,short,long long);
ST(long long,char,long long);
ST(long,int,long);
ST(long,short,long);
ST(long,char,long);
ST(int,short,int);
ST(int,char,int);
ST(short,char,short);

// unsigned integer type rules
ST(unsigned long long,unsigned long,unsigned long long);
ST(unsigned long long,unsigned int,unsigned long long);
ST(unsigned long long,unsigned short,unsigned long long);
ST(unsigned long long,unsigned char,unsigned long long);
ST(unsigned long,unsigned int,unsigned long);
ST(unsigned long,unsigned short,unsigned long);
ST(unsigned long,unsigned char,unsigned long);
ST(unsigned int,unsigned short,unsigned int);
ST(unsigned int,unsigned char,unsigned int);
ST(unsigned short,unsigned char,unsigned short);

template<typename T> struct create_error_msg {};
template<typename L,typename R> struct super_type { typedef typename create_error_msg<super_type<L,R> >::operation_not_allowed type;};


template<typename ScalarT> struct pick_next_math_type;
#define PICK_NEXT_MATH_TYPE(TYPEIN,TYPEOUT) template<> struct pick_next_math_type<TYPEIN>{ typedef TYPEOUT type; };

PICK_NEXT_MATH_TYPE(int,float)
PICK_NEXT_MATH_TYPE(long,float)
PICK_NEXT_MATH_TYPE(char,float)
PICK_NEXT_MATH_TYPE(short,float)
PICK_NEXT_MATH_TYPE(unsigned int,float)
PICK_NEXT_MATH_TYPE(unsigned long,float)
PICK_NEXT_MATH_TYPE(unsigned char,float)
PICK_NEXT_MATH_TYPE(unsigned short,float)
PICK_NEXT_MATH_TYPE(float,float)
PICK_NEXT_MATH_TYPE(double,double)
PICK_NEXT_MATH_TYPE(long double,long double)




template<typename LeftDerived, typename RightDerived, typename ResultValueT, size_t Rows, size_t Columns> struct add; 
template<typename LeftDerived, typename RightDerived, typename ResultValueT, size_t Rows, size_t Columns>  struct sub;
template<typename LeftDerived, typename RightDerived, typename ResultValueT, size_t Rows, size_t Columns> struct mul ; 
template<typename LeftDerived, typename ValueT, size_t Rows, size_t Columns> struct scale;
template<typename Derived, typename ValueT, size_t Rows, size_t Columns> struct neg;
template<typename Derived, typename ValueT, size_t Rows, size_t Columns> struct transp;
/**
 * Basetype of expression template system.
 */
template<typename Derived, typename ValueT, size_t Rows, size_t Columns> 
struct expr_type {
  static const size_t num_rows = Rows, num_columns = Columns;
  typedef ValueT value_type;
  inline Derived const& derived() const { return *static_cast<Derived const*>(this); }
  inline Derived & derived() { return *static_cast<Derived *>(this); }
};

template<typename LeftDerived,typename LeftValueT,size_t Rows,size_t Columns,typename RightDerived,typename RightValueT>
inline add<LeftDerived,RightDerived,typename super_type<LeftValueT,RightValueT>::type, Rows, Columns> 
//operator +( expr_type<LeftDerived, LeftValueT, Rows, Columns> const& left, expr_type<RightDerived, RightValueT, Rows, Columns> const& right )
operator +( expr_type<LeftDerived, LeftValueT, Rows, Columns> const& left, expr_type<RightDerived, RightValueT, Rows, Columns> const& right )
{ return add<LeftDerived,RightDerived,typename super_type<LeftValueT,RightValueT>::type, Rows, Columns >(left.derived(), right.derived()); }

template<typename LeftDerived,typename LeftValueT,size_t Rows, size_t Columns,typename RightDerived,typename RightValueT>
inline sub<LeftDerived,RightDerived,typename super_type<LeftValueT,RightValueT>::type, Rows, Columns> 
operator-(expr_type<LeftDerived, LeftValueT, Rows, Columns> const& left, expr_type<RightDerived, RightValueT, Rows, Columns> const& right )
{ return sub<LeftDerived,RightDerived,typename super_type<LeftValueT,RightValueT>::type, Rows, Columns >(left.derived(), right.derived()); }

template<typename LeftDerived,typename LeftValueT,size_t LeftRows, size_t LeftColumns,typename RightDerived,typename RightValueT, size_t RightColumns>
inline mul<LeftDerived,RightDerived,typename super_type<LeftValueT,RightValueT>::type, LeftRows, RightColumns> 
operator*(expr_type<LeftDerived, LeftValueT, LeftRows, LeftColumns> const& left, expr_type<RightDerived, RightValueT, LeftRows, RightColumns> const& right )
{ return mul<LeftDerived,RightDerived,typename super_type<LeftValueT,RightValueT>::type, LeftRows, RightColumns>(left.derived(), right.derived()); }

template<typename LeftDerived,typename LeftValueT,size_t LeftRows, size_t LeftColumns>
inline scale<LeftDerived,LeftValueT, LeftRows, LeftColumns> 
operator*(expr_type<LeftDerived, LeftValueT, LeftRows, LeftColumns> const& left, LeftValueT const& right )
{ return scale<LeftDerived,LeftValueT, LeftRows, LeftColumns> (left.derived(), right); }

template<typename Derived,typename ValueT,size_t Rows, size_t Columns>
inline neg<Derived,ValueT,Rows,Columns> 
operator-(expr_type<Derived,ValueT, Rows,Columns> const& left )
{ return neg<Derived,ValueT,Rows,Columns>(left.derived()); }

template<typename Derived,typename ValueT,size_t Rows, size_t Columns>
inline transp<Derived,ValueT,Rows,Columns> 
transpose(expr_type<Derived,ValueT, Rows,Columns> const& left )
{ return transp<Derived,ValueT,Rows,Columns>(left.derived()); }



  
/**
 * \brief fixed matrix is an optimized matrix, which uses a static stack array 
 * to store the fixed number of scalars.
 * This matrix should be used instead of tMatrix everywhere the    
 */
template<typename ValueT, size_t Rows, size_t Columns>
struct fixed_matrix : public expr_type<fixed_matrix<ValueT,Rows,Columns>, ValueT, Rows, Columns> {
  typedef ValueT value_type;
  value_type data[Rows][Columns];

  static const size_t num_rows = Rows, num_columns = Columns;

  explicit fixed_matrix( value_type const& initializer ) { std::fill(&(data[0][0]), &(data[0][0]) + Rows*Columns, initializer );  }

  template<typename OtherValueT>
  fixed_matrix( profdist::fixed_matrix<OtherValueT, Rows, Columns> const& cp ) {
    std::copy( &(cp.data[0][0]), &(cp.data[0][0]) + Rows*Columns, &(data[0][0]) );
  }

  
  template<typename Derived, typename OtherValueT>
  fixed_matrix( profdist::expr_type<Derived, OtherValueT, Rows, Columns> const& cp ) {
    std::size_t i=Rows;
    do { --i; std::size_t j = Columns;
      do { --j; 
        data[i][j] = cp.derived()(i,j);
      } while (j);
    } while(i);
  }
  
  inline value_type const& operator()(std::size_t const row, std::size_t const column ) const { return data[row][column]; }
  inline value_type &operator()(std::size_t const row, std::size_t const column ) { return data[row][column]; }

  inline value_type const *operator[](std::size_t const row ) const { return data[row]; }
  inline value_type *operator[](std::size_t const row ) { return data[row]; }

  inline fixed_matrix<ValueT, Rows, Columns>& operator =(fixed_matrix<ValueT, Rows, Columns> const& other ){
    std::copy( &(other.data[0][0]), &(other.data[0][0]) + Rows*Columns, &(data[0][0]) );
    return *this;
  }

  template<typename Derived, typename OtherValueT>
  inline  fixed_matrix<ValueT, Rows, Columns>& operator=( profdist::expr_type<Derived, OtherValueT, Rows, Columns> const& other ) {
    fixed_matrix<OtherValueT,Rows,Columns> temp( other.derived() );
    std::copy( &(temp.data[0][0]), &(temp.data[0][0]) + Rows*Columns, &(data[0][0]) );
    return *this;
  }

  template<typename Derived, typename OtherValueT>
  inline fixed_matrix<ValueT, Rows, Columns>& operator+=( profdist::expr_type<Derived, OtherValueT, Rows, Columns> const& other ) {
    fixed_matrix<OtherValueT,Rows,Columns> temp( *this + other.derived() );
    std::copy( &(temp.data[0][0]), &(temp.data[0][0]) + Rows*Columns, &(data[0][0]) );
    return *this;
  }

  template<typename Derived, typename OtherValueT>
  inline  fixed_matrix<ValueT, Rows, Columns>& operator-=( profdist::expr_type<Derived, OtherValueT, Rows, Columns> const& other ) {
    fixed_matrix<OtherValueT,Rows,Columns> temp( *this - other.derived() );
    std::copy( &(temp.data[0][0]), &(temp.data[0][0]) + Rows*Columns, &(data[0][0]) );
    return *this;
  }

  template<typename Derived, typename OtherValueT>
  inline fixed_matrix<ValueT, Rows, Columns>& operator*=( profdist::expr_type<Derived, OtherValueT, Rows, Columns> const& other) {
    fixed_matrix<ValueT,Rows,Columns> copy(*this);
    *this = copy*other;
    return *this;
  }


  inline  fixed_matrix<ValueT, Rows, Columns>& operator*=( ValueT const& other ) {
    std::size_t i=Rows;
    do {
      --i;
      std::size_t j = Columns;
      do {
        --j;
        data[i][j] *= other;
      } while (j);
    } while(i);
    return *this;
  }
};


template<typename ValueT>
struct fixed_matrix<ValueT,4,4> : public expr_type<fixed_matrix<ValueT,4,4>, ValueT, 4,4> {
  typedef ValueT value_type;
  value_type data[4][4];

  static const size_t num_rows = 4, num_columns = 4;

  explicit fixed_matrix( value_type const& initializer ) { std::fill(&(data[0][0]), &(data[0][0]) + 16, initializer );  }

   
#define UNROLL_44_EXPR_F(op,param) data[0][0] op param(0,0);\
  data[0][1] op param(0,1);\
  data[0][2] op param(0,2);\
  data[0][3] op param(0,3);\
  data[1][0] op param(1,0);\
  data[1][1] op param(1,1);\
  data[1][2] op param(1,2);\
  data[1][3] op param(1,3);\
  data[2][0] op param(2,0);\
  data[2][1] op param(2,1);\
  data[2][2] op param(2,2);\
  data[2][3] op param(2,3);\
  data[3][0] op param(3,0);\
  data[3][1] op param(3,1);\
  data[3][2] op param(3,2);\
  data[3][3] op param(3,3)

  template<typename OtherValueT>
  fixed_matrix( profdist::fixed_matrix<OtherValueT, 4, 4> const& cp ) {
    UNROLL_44_EXPR_F(=,cp);
  }

  template<typename Derived, typename OtherValueT>
  fixed_matrix( profdist::expr_type<Derived, OtherValueT, 4, 4> const& cp ) {
    UNROLL_44_EXPR_F(=,(cp.derived()));
  }
  
  inline value_type const& operator()(std::size_t const row, std::size_t const column ) const { return data[row][column]; }
  inline value_type &operator()(std::size_t const row, std::size_t const column ) { return data[row][column]; }

  inline value_type const *operator[](std::size_t const row ) const { return data[row]; }
  inline value_type *operator[](std::size_t const row ) { return data[row]; }

  inline fixed_matrix<ValueT, 4, 4>& operator =(fixed_matrix<ValueT, 4, 4> const& other ){
    UNROLL_44_EXPR_F(=,other);
    return *this;
  }

  template<typename Derived, typename OtherValueT>
  inline  fixed_matrix<ValueT, 4, 4>& operator=( profdist::expr_type<Derived, OtherValueT, 4, 4> const& other ) {
    fixed_matrix<ValueT,4,4> temp( other.derived() );
    std::copy( &(temp.data[0][0]), &(temp.data[0][0]) + 16 , &(data[0][0]) );
    return *this;
  }

  template<typename Derived, typename OtherValueT>
  inline fixed_matrix<ValueT, 4, 4>& operator+=( profdist::expr_type<Derived, OtherValueT, 4, 4> const& other ) {
    fixed_matrix<ValueT,4,4> temp( *this + other.derived() );
    std::copy( &(temp.data[0][0]), &(temp.data[0][0]) + 16 , &(data[0][0]) );
    return *this;
  }

  template<typename Derived, typename OtherValueT>
  inline  fixed_matrix<ValueT, 4, 4>& operator-=( profdist::expr_type<Derived, OtherValueT, 4, 4> const& other ) {
    fixed_matrix<ValueT,4,4> temp( *this - other.derived() );
    std::copy( &(temp.data[0][0]), &(temp.data[0][0]) + 16 , &(data[0][0]) );
    return *this;
  }

  template<typename Derived, typename OtherValueT>
  inline fixed_matrix<ValueT, 4, 4>& operator*=( profdist::expr_type<Derived, OtherValueT, 4, 4> const& other) {
    fixed_matrix<ValueT,4,4> copy(*this);
    *this = copy*other;
    return *this;
  }


  inline  fixed_matrix<ValueT, 4, 4>& operator*=( ValueT const& param ) {
    data[0][0] *= param;
    data[0][1] *= param;
    data[0][2] *= param;
    data[0][3] *= param;
    data[1][0] *= param;
    data[1][1] *= param;
    data[1][2] *= param;
    data[1][3] *= param;
    data[2][0] *= param;
    data[2][1] *= param;
    data[2][2] *= param;
    data[2][3] *= param;
    data[3][0] *= param;
    data[3][1] *= param;
    data[3][2] *= param;
    data[3][3] *= param;
    return *this;
  }
#undef UNROLL_44_EXPR_F
};


template<typename LeftDerived, typename RightDerived, typename ResultValueT, size_t Rows, size_t Columns> 
struct add : public expr_type<add<LeftDerived,RightDerived,ResultValueT,Rows,Columns>, ResultValueT, Rows, Columns> {
  LeftDerived const& left;
  RightDerived const& right ; 
  add( LeftDerived const & l, RightDerived const &r ) :left(l), right(r){}
  typedef ResultValueT value_type;
  inline value_type operator()(std::size_t const row, std::size_t const column ) const 
  { return value_type( left(row,column) ) + value_type( right(row, column) ); }
};

template<typename LeftDerived, typename RightDerived, typename ResultValueT, size_t Rows, size_t Columns> 
struct sub : public expr_type<sub<LeftDerived,RightDerived,ResultValueT,Rows,Columns>, ResultValueT, Rows, Columns > {
  LeftDerived const& left;
  RightDerived const& right ; 
  sub( LeftDerived const & l, RightDerived const &r ) :left(l), right(r){}
  typedef ResultValueT value_type;
  inline value_type operator()(std::size_t const row, std::size_t const column ) const 
  { return value_type( left(row,column) ) - value_type( right(row, column) ); }
};

template<typename LeftDerived, typename RightDerived, typename ResultValueT, size_t Rows, size_t Columns> 
struct mul : public expr_type<mul<LeftDerived,RightDerived,ResultValueT,Rows,Columns>, ResultValueT, Rows, Columns > {
  LeftDerived const& left;
  RightDerived const& right ; 
  mul( LeftDerived const & l, RightDerived const &r ) :left(l), right(r){}
  typedef ResultValueT value_type;
  inline value_type operator()(std::size_t const row, std::size_t const column ) const 
  { value_type ret = 0;
    std::size_t i = LeftDerived::num_columns;
    do{  --i;
      ret += value_type( left(row,i) ) * value_type( right(i, column) ); 
    }
    while(i);
    return ret;
  }
};

template<typename LeftDerived, typename ValueT, size_t Rows, size_t Columns> 
struct scale : public expr_type<scale<LeftDerived,ValueT,Rows,Columns>, ValueT, Rows, Columns > {
  LeftDerived const& left;
  ValueT const& right; 
  scale( LeftDerived const & l, ValueT const &r ) :left(l), right(r){}
  typedef ValueT value_type;
  inline value_type operator()(std::size_t const row, std::size_t const column ) const { return left(row,column) * right;  }
};

template<typename Derived, typename ValueT, size_t Rows, size_t Columns>
struct neg: public expr_type<neg<Derived,ValueT,Rows,Columns>, ValueT, Rows, Columns > {
  Derived const& left;
  neg( Derived const & l) :left(l){}
  typedef ValueT value_type;
  inline value_type operator()(std::size_t const row, std::size_t const column ) const { return -left(row,column);  }
};

template<typename Derived, typename ValueT, size_t Rows, size_t Columns>
struct transp: public expr_type<transp<Derived,ValueT,Rows,Columns>, ValueT, Rows, Columns > {
  Derived const& left;
  transp( Derived const & l) :left(l){}
  typedef ValueT value_type;
  inline value_type operator()(std::size_t const row, std::size_t const column ) const { return left(column,row);  }
};

template<typename Vt, size_t Rows, size_t Cols>
std::ostream & operator << ( std::ostream & obj, profdist::fixed_matrix<Vt,Rows,Cols> const& mat ) {
  for( std::size_t i = 0; i < Rows; ++i ) {
    for( std::size_t j = 0; j < Cols; ++j ) 
      std::cout << mat(i,j) << ", ";
    std::cout << std::endl;
  }
}


}


#endif

