//
// Some basic functions to make blitz::TinyMatrix and blitz::TinyVector more comfortable
//
//  2010 Olaf Ronneberger

#include <blitz/array.h>
#include <blitz/tinyvec-et.h>		// uncomment when using blitz++ version newer than 0.9

/*======================================================================*/
/*! 
 *   transpose a matrix of arbitrary dimension
 *
 *   \param a  matrix
 *
 *   \exception none
 *
 *   \return transposed matrix
 */
/*======================================================================*/
template< typename T, int I, int J>
inline
blitz::TinyMatrix<T,J,I> mytranspose(
    const blitz::TinyMatrix<T,I,J>& a)
{
  blitz::TinyMatrix<T,J,I> t;
  for( int i = 0; i < I; ++i)
  {
    for( int j = 0; j < J; ++j)
    {
      t(j,i) = a(i,j);
    }
  }
  return t;
}


/*======================================================================*/
/*! 
 *   product of two matrices with arbitrary dimensions
 *
 *   \param a  first matrix
 *   \param b  second matrix
 *
 *   \exception none
 *
 *   \return product of a and b
 */
/*======================================================================*/
template< typename T, int I, int J, int K>
inline
blitz::TinyMatrix<T,I,J> myproduct(
    const blitz::TinyMatrix<T,I,K>& a,
    const blitz::TinyMatrix<T,K,J>& b)
{
  blitz::TinyMatrix<T,I,J> c;
    
  for( int i = 0; i < I; ++i)
  {
    for( int j = 0; j < J; ++j)
    {
      T val = 0;
      for( int k = 0; k < K; ++k)
      {
        val += a(i,k) * b(k,j);
      }
      c(i,j) = val;
    }
  }
  return c;
}


/*======================================================================*/
/*! 
 *   product of  matrices with a vector of arbitrary dimensions
 *
 *   \param a  matrix
 *   \param b  vector
 *
 *   \exception none
 *
 *   \return product of a and b
 */
/*======================================================================*/
template< typename T, int I, int J>
inline
blitz::TinyVector<T,I> myproduct(
    const blitz::TinyMatrix<T,I,J>& a,
    const blitz::TinyVector<T,J>& b)
{
  blitz::TinyVector<T,I> c;
    
  for( int i = 0; i < I; ++i)
  {
    T val = 0;
    for( int j = 0; j < J; ++j)
    {
      val += a(i,j) * b(j);
    }
    c(i) = val;
    
  }
  return c;
}

/*======================================================================*/
/*! 
 *   Specialization of _bz_OneZeroTraits for
 *   blitz::TinyVector<float,3> to allow to use blitz::mean() for an
 *   Array of TinyVectors
 */
/*======================================================================*/
namespace blitz
{
  
template<>
struct _bz_OneZeroTraits<blitz::TinyVector<float,3> > {
  static inline blitz::TinyVector<float,3> zero() { return blitz::TinyVector<float,3>(0,0,0); }
  static inline blitz::TinyVector<float,3> one()  { return blitz::TinyVector<float,3>(1,1,1); }
};

template<>
struct _bz_OneZeroTraits<blitz::TinyVector<float,4> > {
  static inline blitz::TinyVector<float,4> zero() { return blitz::TinyVector<float,4>(0,0,0,0); }
  static inline blitz::TinyVector<float,4> one()  { return blitz::TinyVector<float,4>(1,1,1,1); }
};

}