#ifndef __MATRIX_HXX__
#define __MATRIX_HXX__

#include <iostream>
#include <utility> //for pair<T,U>()
#include <cmath>   //for sqrt()
#include <exception>
using namespace std;

#ifdef USE_ALTIVEC
#include <vecLib/vDSP.h>
#endif

class BadFFT : public exception
{
  public:
    virtual ~BadFFT() throw() {};
    virtual const char* what(void) const throw()
            {return "Bad FFT : input and output matrices should have same size,\n"
	            "which should be a power of 2";}
};

//class Matrix
template<typename T> class Matrix;
template<typename T> ostream& operator<<(ostream&,const Matrix<T>&);
template<typename T> istream& operator>>(istream&,Matrix<T>&);

template<typename T = float>
class Matrix
{
  public:
    //Constructors
    Matrix(size_t aHeight=0, size_t aWidth=0, const T& fill=T(),bool aBinaryMode = true)
	 :height(aHeight),width(aWidth),data(0),binaryMode(aBinaryMode)
	 {data = allocation(height,width,fill);}

    Matrix(const Matrix<T>& other)
         :height(other.height),width(other.width),data(0),binaryMode(other.binaryMode)
         {data=allocation(height,width);copy(other);}
    template<typename U> Matrix(const Matrix<U>& other)
                              :height(other.getHeight()),width(other.getWidth()),data(0),
                               binaryMode(other.getBinaryMode())
			 {data=allocation(height,width);copy(other);}
    
    template <typename U>
    Matrix(const Matrix<U>& other, size_t fromRow, size_t fromCol)
         :height(other.getHeight()-fromRow),width(other.getWidth()-fromCol),data(0),
          binaryMode(other.getBinaryMode())
	 {data=allocation(height,width);
          copy(other, fromRow, fromCol, height, width);}
    
    template <typename U>
    Matrix(const Matrix<U>& other, size_t fromRow, size_t fromCol,
                                   size_t aHeight, size_t aWidth,
				   const T& fill = T(), bool loop = false)
         :height(aHeight),width(aWidth),data(0),binaryMode(other.getBinaryMode())
	 {data=allocation(height,width,fill);
	  copy(other, fromRow, fromCol, height, width, loop);}
    
    //operator=
    Matrix<T>& operator=(const Matrix<T>& other)
	       {if (((void*)&other) != ((void*)this))
		   {clean();height=other.height;width=other.width;
		    data=allocation(height,width);copy(other);} return *this;}
    template<typename U> Matrix<T>& operator=(const Matrix<U>& other)
	       {if (((void*)&other) != ((void*)this))
		   {clean();height=other.getHeight();width=other.getWidth();
		    data=allocation(height,width);copy(other);} return *this;}
    
    //destructor
    virtual ~Matrix() {clean();}

  public:
    template<typename U> operator Matrix<U>() const;

  public:
    T*        operator[](size_t index)       {return data[index];}
    const T*  operator[](size_t index) const {return data[index];}
    
  public:
//    friend class <> Matrix;
    template<typename U> Matrix<T>& operator+=(const Matrix<U>& other);
    template<typename U> Matrix<T>& operator-=(const Matrix<U>& other);
    template<typename U> Matrix<T>& operator*=(const Matrix<U>& other);
    template<typename U> Matrix<T>& operator/=(const Matrix<U>& other);
    template<typename U> Matrix<T>& operator+=(const U& x);
    template<typename U> Matrix<T>& operator-=(const U& x);
    template<typename U> Matrix<T>& operator*=(const U& x);
    template<typename U> Matrix<T>& operator/=(const U& x);
    template<typename U> Matrix<T>& operator&=(const Matrix<U>& other);
    template<typename U> Matrix<T>& operator|=(const Matrix<U>& other);
    template<typename U> Matrix<T>& operator^=(const Matrix<U>& other);
    template<typename U> Matrix<T>& operator&=(const U& x);
    template<typename U> Matrix<T>& operator|=(const U& x);
    template<typename U> Matrix<T>& operator^=(const U& x);
    Matrix<T>  operator~(void);
    Matrix<T>& operator++(void) {return (*this)+=static_cast<T>(1);}
    Matrix<T>  operator++(int) {Matrix<T> save(*this); ++(*this); return save;}
    Matrix<T>& operator--(void) {return (*this)-=static_cast<T>(1);}
    Matrix<T>  operator--(int) {Matrix<T> save(*this); --(*this); return save;}
    Matrix<T> operator-() const {Matrix<T> m(*this); return m.oppose();}
    T sum(void) const;
    T product(void) const;

    Matrix<T>& transpose(void);

  public:
    //resize
    Matrix<T>& resize(size_t height, size_t width, const T& fill=T())
               {return resize(height,width,0,0,fill);}
    Matrix<T>& resize(size_t height, size_t width,
		      size_t topMargin, size_t leftMargin, const T& fill=T());

    //subMatrix
    Matrix<T>  subMatrix(size_t fromRow, size_t fromCol,
                         bool loop = false) const
	       {if (loop) {if (height) fromRow%=height; if (width) fromCol%=width;}
				  return subMatrix(fromRow,fromCol,
			                     height>fromRow ? height-fromRow : 0,
			                     width >fromCol ? width-fromCol  : 0,
													 T(),loop);}
    Matrix<T>  subMatrix(size_t fromRow, size_t fromCol,
		         size_t height,  size_t width,
			 const T& fill = T(), bool loop = false) const
	       {return Matrix(*this,fromRow,fromCol,height,width,fill,loop);}
    
    //fillWith
    Matrix<T>& reset(const T& fill = T()) {return fillWith(fill);}
    Matrix<T>& fillWith(const T& fill) {if (data) memorySet(data[0],height*width,fill);return *this;}
    
    template<typename U> Matrix<T>& fillWith(const Matrix<U>& other, bool loop = false)
            {return fillWith(other,0,0,loop);}
    template<typename U> Matrix<T>& fillWith(const Matrix<U>& other,
		                             size_t fromRow, size_t fromCol,
					     bool loop = false)
            {return fillWith(0,0,other,fromRow,fromCol,height,width,loop);}
    
    template<typename U> Matrix<T>& fillWith(const Matrix<U>& other,
		                             size_t fromRow, size_t fromCol,
					     size_t height,  size_t width,
					     bool loop = false)
            {return fillWith(0,0,other,fromRow,fromCol,height,width,loop);}
    
    Matrix<T>& fillWith(size_t targetRow, size_t targetCol,
		        const Matrix<T>& other, bool loop = false)
	       {return fillWith(targetRow,targetCol,other,0,0,loop);}
    
    template <typename U>
    Matrix<T>& fillWith(size_t targetRow, size_t targetCol,
		        const Matrix<U>& other, bool loop = false)
	       {return fillWith(targetRow,targetCol,other,0,0,loop);}

    Matrix<T>& fillWith(size_t targetRow, size_t targetCol,
		        const Matrix<T>& other,
		        size_t fromRow=0, size_t fromCol=0,bool loop = false)
	       {return fillWith(targetRow,targetCol,other,fromRow,fromCol,
			        loop ? (height) :
				       ((other.getHeight()>fromRow) ? other.getHeight()-fromRow : 0),
			        loop ? (width) :
				       ((other.getWidth()>fromCol) ? other.getWidth()-fromCol : 0),
				loop);}
    
    template<typename U> Matrix<T>& fillWith(size_t targetRow, size_t targetCol,
		                             const Matrix<U>& other,
		                             size_t fromRow=0, size_t fromCol=0,
					     bool loop = false)
	       {return fillWith(targetRow,targetCol,other,fromRow,fromCol,
			        loop ? (height) :
				       ((other.getHeight()>fromRow) ? other.getHeight()-fromRow : 0),
			        loop ? (width) :
				       ((other.getWidth()>fromCol) ? other.getWidth()-fromCol : 0),
				loop);}

    Matrix<T>& fillWith(size_t targetRow, size_t targetCol,
		        const Matrix<T>& other,
		        size_t fromRow, size_t fromCol,
			size_t height,  size_t width,
			bool loop = false);
    
    template<typename U> Matrix<T>& fillWith(size_t targetRow, size_t targetCol,
		                             const Matrix<U>& other,
		                             size_t fromRow, size_t fromCol,
					     size_t height,  size_t width,
					     bool loop = false);

    Matrix<T>& shift(int deltaHeight, int deltaWidth, const T& fill=T());
    Matrix<T>&  loop(int deltaHeight, int deltaWidth);

    Matrix<T>& flip(bool bottomUp = true, bool leftRight = true);
    
    template <typename U,typename V> Matrix<T>& operator()( U (*f)(const V&));
    template <typename U,typename V> Matrix<T>& operator()( U (*f)(V) );
    template <typename U> Matrix<T>& operator()(const U&);
    pair<T,T> meanAndStandardDeviation(void) const;
    T mean(void) const;
    T standardDeviation(void) const {return meanAndStandardDeviation().second;}
    T max(void) const;
    T min(void) const;
    T norm1(void) const;
    T norm2(void) const;
    T normInf(void) const;
    Matrix<T>& oppose(void);
    Matrix<T>& inverse(void);
    Matrix<T>& sqr(void);
    Matrix<T>& sqrt(void) {return (*this)( (double (*)(double)) &std::sqrt);}
    Matrix<T>& abs(void);
		template<typename U,typename V>
    Matrix<T>& mulAdd(const Matrix<U>& m1, const Matrix<V>& m2);
		Matrix<T>& pow(long n);
		Matrix<T>& pow(unsigned long n);
		Matrix<T>& pow(int n)            {return pow(static_cast<long>(n));}
		Matrix<T>& pow(unsigned int  n)  {return pow(static_cast<unsigned long>(n));}
		Matrix<T>& pow(short n)          {return pow(static_cast<long>(n));}
		Matrix<T>& pow(unsigned short n) {return pow(static_cast<unsigned long>(n));}
		Matrix<T>& pow(char n)           {return pow(static_cast<long>(n));}
		Matrix<T>& pow(unsigned char n)  {return pow(static_cast<unsigned long>(n));}
		Matrix<T>& pow(double y);
    bool any(void) const;
    bool none(void) const {return !any();}
    bool every(void) const;
    template<typename U> T dotProduct(const Matrix<U>& other) const;
    template<typename U> Matrix<T> mulByVec(const Matrix<U>& lineMatrix) const;

  public:
    Matrix<T>& writeToStream(ostream&) const;
    Matrix<T>& readFromStream(istream&);
    friend ostream& operator<< <> (ostream&,const Matrix<T>&);
    friend istream& operator>> <> (istream&,Matrix<T>&);
    
  public:
    size_t getHeight()     const {return height;}
    size_t getWidth()      const {return width;}
    T**          getData()       const {return data;}
    size_t getDataLength() const {return height*width;}
    bool  getBinaryMode(void)    const {return binaryMode;}
    void  setBinaryMode(bool mode)     {binaryMode = mode;}

  public:
    static T DotProduct(const Matrix<T>& m1, const Matrix<T>& m2)
		       {return m1.dotProduct(m2);}

    static void FFT2D(Matrix<T>* Output_real, Matrix<T>* Output_imag,
		      const Matrix<T>& Input_real, const Matrix<T>& Input_imag);
    static void IFFT2D(Matrix<T>* Output_real, Matrix<T>* Output_imag,
		       const Matrix<T>& Input_real, const Matrix<T>& Input_imag);

    #ifdef USE_ALTIVEC
    static void FFT2D(Matrix<T>* Output_real, Matrix<T>* Output_imag,
		      const Matrix<T>& Input_real, const Matrix<T>& Input_imag,
		      const FFTSetup& setup);
    static void IFFT2D(Matrix<T>* Output_real, Matrix<T>* Output_imag,
		       const Matrix<T>& Input_real, const Matrix<T>& Input_imag,
		       const FFTSetup& setup);
    #endif

  private:
    size_t height;
    size_t width;
    T** data;
    bool binaryMode;

  private:
    T**  allocation(size_t height, size_t width, const T& fill= T());
    void clean(void);
    void memorySet(T* data,size_t size, const T& fill = T());
    void memoryCopy(T* dest, const T* src, size_t size);
    void memoryMove(T* dest, const T* src, size_t size);

    void copy(const Matrix<T>& other);
    template<typename U> void copy(const Matrix<U>& other);

    void copy(const Matrix<T>& other, size_t fromRow, size_t fromCol,
		                      size_t height,  size_t width,
				      bool loop = false);
    template <typename U>
    void copy(const Matrix<U>& other, size_t fromRow, size_t fromCol,
		                      size_t height,  size_t width,
				      bool loop = false);
  private:
    static T abs(const T&);
    static bool CheckFFTOk(const Matrix<T>* Output_real, const Matrix<T>* Output_imag,
                           const Matrix<T>& Input_real , const Matrix<T>& Input_imag );
    static void Four2(Matrix<T>* fftr, Matrix<T>* ffti,
		      Matrix<T>* datar,Matrix<T>* datai, int isign);
    static void Four1(Matrix<T>* data, int isign);
    template<typename U> static T DotProduct(const T* data1, const U* data2, size_t size);
    static size_t LogTwoOfCeilingPowerOfTwo(size_t x);
};

template <typename T>
inline Matrix<T> operator+(const Matrix<T>& m1, const Matrix<T>& m2) {Matrix<T> m(m1); return m+=m2;}
template <typename T>
inline Matrix<T> operator-(const Matrix<T>& m1, const Matrix<T>& m2) {Matrix<T> m(m1); return m-=m2;}
template <typename T>
inline Matrix<T> operator*(const Matrix<T>& m1, const Matrix<T>& m2) {Matrix<T> m(m1); return m*=m2;}
template <typename T>
inline Matrix<T> operator/(const Matrix<T>& m1, const Matrix<T>& m2) {Matrix<T> m(m1); return m/=m2;}

template <typename T>
inline Matrix<T> operator&(const Matrix<T>& m1, const Matrix<T>& m2) {Matrix<T> m(m1); return m&=m2;}
template <typename T>
inline Matrix<T> operator|(const Matrix<T>& m1, const Matrix<T>& m2) {Matrix<T> m(m1); return m|=m2;}
template <typename T>
inline Matrix<T> operator^(const Matrix<T>& m1, const Matrix<T>& m2) {Matrix<T> m(m1); return m^=m2;}

template <typename T,typename U>
inline Matrix<T> operator+(const Matrix<T>& m1, U x) {Matrix<T> m(m1); return m+=x;}
template <typename T,typename U>
inline Matrix<T> operator+(U x, const Matrix<T>& m1) {Matrix<T> m(m1); return m+=x;}

template <typename T,typename U>
inline Matrix<T> operator-(const Matrix<T>& m1, U x) {Matrix<T> m(m1); return m-=x;}
template <typename T,typename U>
inline Matrix<T> operator-(U x, const Matrix<T>& m1) {Matrix<T> m(m1); return m-=x;}

template <typename T,typename U>
inline Matrix<T> operator*(const Matrix<T>& m1, U x) {Matrix<T> m(m1); return m*=x;}
template <typename T,typename U>
inline Matrix<T> operator*(U x, const Matrix<T>& m1) {Matrix<T> m(m1); return m*=x;}

template <typename T,typename U>
inline Matrix<T> operator/(const Matrix<T>& m1, U x) {Matrix<T> m(m1); return m/=x;}
template <typename T,typename U>
inline Matrix<T> operator/(U x, const Matrix<T>& m1) {Matrix<T> m(m1); return m/=x;}

template <typename T,typename U>
inline Matrix<T> operator&(const Matrix<T>& m1, U x) {Matrix<T> m(m1); return m&=x;}
template <typename T,typename U>
inline Matrix<T> operator&(U x, const Matrix<T>& m1) {Matrix<T> m(m1); return m&=x;}
template <typename T,typename U>
inline Matrix<T> operator|(const Matrix<T>& m1, U x) {Matrix<T> m(m1); return m|=x;}
template <typename T,typename U>
inline Matrix<T> operator|(U x, const Matrix<T>& m1) {Matrix<T> m(m1); return m|=x;}
template <typename T,typename U>
inline Matrix<T> operator^(const Matrix<T>& m1, U x) {Matrix<T> m(m1); return m^=x;}
template <typename T,typename U>
inline Matrix<T> operator^(U x, const Matrix<T>& m1) {Matrix<T> m(m1); return m^=x;}

template <typename T, typename U>
inline bool operator==(const Matrix<T>& m1, const Matrix<U>& m2);
template <typename T, typename U>
inline bool operator!=(const Matrix<T>& m1, const Matrix<U>& m2) {return !(m1==m2);}
template <typename T, typename U>
inline bool operator<(const Matrix<T>& m1, const Matrix<U>& m2);
template <typename T, typename U>
inline bool operator<=(const Matrix<T>& m1, const Matrix<U>& m2);
template <typename T, typename U>
inline bool operator>(const Matrix<T>& m1, const Matrix<U>& m2);
template <typename T, typename U>
inline bool operator>=(const Matrix<T>& m1, const Matrix<U>& m2);

template<typename T> inline ostream& operator<<(ostream&,const Matrix<T>&);
template<typename T> inline istream& operator>>(istream& stream, Matrix<T>& m)
                     {m.readFromStream(stream);return stream;}

#endif