/*!
 * \file
 * \brief Matrices in decomposed forms (LDL', LU, UDU', etc).
 * \author Vaclav Smidl.
 *
 * -----------------------------------
 * BDM++ - C++ library for Bayesian Decision Making under Uncertainty
 *
 * Using IT++ for numerical operations
 * -----------------------------------
 */

#ifndef DC_H
#define DC_H


#include "../itpp_ext.h"

/*!
\defgroup math Auxiliary math functions 
@{
*/

using namespace itpp;

//! Auxiliary function dydr; dyadic reduction
void dydr( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx );

//! Auxiliary function ltuinv; inversion of a triangular matrix;
//TODO can be done via: dtrtri.f from lapack
mat ltuinv( const mat &L );

/*! \brief Abstract class for representation of double symmetric matrices in square-root form.

All operations defined on this class should be optimized for the chosen decomposition.
*/
class sqmat
{
	public:
		/*!
		 * Perfroms a rank-1 update by outer product of vectors: \f$V = V + w v v'\f$.
		 * @param  v Vector forming the outer product to be added
		 * @param  w weight of updating; can be negative

		 BLAS-2b operation.
		 */
		virtual void opupdt ( const vec &v, double w ) =0;

		/*! \brief Conversion to full matrix.
		*/

		virtual mat to_mat() const =0;

		/*! \brief Inplace symmetric multiplication by a SQUARE matrix \f$C\f$, i.e. \f$V = C*V*C'\f$
		@param C multiplying matrix,
		*/
		virtual void mult_sym ( const mat &C ) =0;
		
		/*! \brief Inplace symmetric multiplication by a SQUARE transpose of matrix \f$C\f$, i.e. \f$V = C'*V*C\f$
		@param C multiplying matrix,
		*/
		virtual void mult_sym_t ( const mat &C ) =0;


		/*!
		\brief Logarithm of a determinant.

		*/
		virtual double logdet() const =0;

		/*!
		\brief Multiplies square root of \f$V\f$ by vector \f$x\f$.

		Used e.g. in generating normal samples.
		*/
		virtual vec sqrt_mult (const vec &v ) const =0;

		/*!
		\brief Evaluates quadratic form \f$x= v'*V*v\f$;

		*/
		virtual double qform (const vec &v ) const =0;

		/*!
		\brief Evaluates quadratic form \f$x= v'*inv(V)*v\f$;

		*/
		virtual double invqform (const vec &v ) const =0;

//	//! easy version of the
//	sqmat inv();

		//! Clearing matrix so that it corresponds to zeros.
		virtual void clear() =0;

		//! Reimplementing common functions of mat: cols().
		int cols() const {return dim;};

		//! Reimplementing common functions of mat: rows().
		int rows() const {return dim;};

		//! Destructor for future use;
		virtual ~sqmat(){};
		//! Default constructor
		sqmat(const int dim0): dim(dim0){};
		//! Default constructor
		sqmat(): dim(0){};
	protected:
		//! dimension of the square matrix
		int dim;
};


/*! \brief Fake sqmat. This class maps sqmat operations to operations on full matrix.

This class can be used to compare performance of algorithms using decomposed matrices with perormance of the same algorithms using full matrices;
*/
class fsqmat: public sqmat
{
	protected:
		//! Full matrix on which the operations are performed
		mat M;
	public:
		void opupdt ( const vec &v, double w );
		mat to_mat() const;
		void mult_sym ( const mat &C);
		void mult_sym_t ( const mat &C);
		//! store result of \c mult_sym in external matrix \f$U\f$
		void mult_sym ( const mat &C, fsqmat &U) const;
		//! store result of \c mult_sym_t in external matrix \f$U\f$
		void mult_sym_t ( const mat &C, fsqmat &U) const;
		void clear();

		//! Default initialization
		fsqmat(){}; // mat will be initialized OK
		//! Default initialization with proper size
		fsqmat(const int dim0); // mat will be initialized OK
		//! Constructor
		fsqmat ( const mat &M );
		//! Constructor
		fsqmat ( const fsqmat &M, const ivec &perm ):sqmat(M.rows()){it_error("not implemneted");};
		//! Constructor
		fsqmat ( const vec &d ):sqmat(d.length()){M=diag(d);};

		//! Destructor for future use;
		virtual ~fsqmat(){};


		/*! \brief Matrix inversion preserving the chosen form.

		@param Inv a space where the inverse is stored.

		*/
		void inv ( fsqmat &Inv ) const;

		double logdet() const {return log ( det ( M ) );};
		double qform (const  vec &v ) const {return ( v* ( M*v ) );};
		double invqform (const  vec &v ) const {return ( v* ( itpp::inv(M)*v ) );};
		vec sqrt_mult (const vec &v ) const {mat Ch=chol(M); return Ch*v;};

		//! Add another matrix in fsq form with weight w
		void add ( const fsqmat &fsq2, double w=1.0 ){M+=fsq2.M;};

		//! Access functions
		void setD (const vec &nD){M=diag(nD);}
		//! Access functions
		vec getD (){return diag(M);}
		//! Access functions
		void setD (const vec &nD, int i){for(int j=i;j<nD.length();j++){M(j,j)=nD(j-i);}} //Fixme can be more general


		//! add another fsqmat matrix
		fsqmat& operator += ( const fsqmat &A ) {M+=A.M;return *this;};
		//! subtrack another fsqmat matrix
		fsqmat& operator -= ( const fsqmat &A ) {M-=A.M;return *this;};
		//! multiply by a scalar
		fsqmat& operator *= ( double x ) {M*=x;return *this;};
//		fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;};
		//! print full matrix
		friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq );

};

/*! \brief Matrix stored in LD form, (commonly known as UD) 

Matrix is decomposed as follows: \f[M = L'DL\f] where only \f$L\f$ and \f$D\f$ matrices are stored.
All inplace operations modifies only these and the need to compose and decompose the matrix is avoided.
*/
class ldmat: public sqmat
{
	public:
		//! Construct by copy of L and D.
		ldmat ( const mat &L, const vec &D );
		//! Construct by decomposition of full matrix V.
		ldmat (const mat &V );
		//! Construct by restructuring of V0 accordint to permutation vector perm.
		ldmat (const ldmat &V0, const ivec &perm):sqmat(V0.rows()){	ldform(V0.L.get_cols(perm), V0.D);};
		//! Construct diagonal matrix with diagonal D0
		ldmat ( vec D0 );
		//!Default constructor
		ldmat ();
		//! Default initialization with proper size
		ldmat(const int dim0);
		
		//! Destructor for future use;
		virtual ~ldmat(){};

		// Reimplementation of compulsory operatios

		void opupdt ( const vec &v, double w );
		mat to_mat() const;
		void mult_sym ( const mat &C);
		void mult_sym_t ( const mat &C);
		//! Add another matrix in LD form with weight w
		void add ( const ldmat &ld2, double w=1.0 );
		double logdet() const;
		double qform (const vec &v ) const;
		double invqform (const vec &v ) const;
		void clear();
		vec sqrt_mult ( const vec &v ) const;

		
		/*! \brief Matrix inversion preserving the chosen form.
		@param Inv a space where the inverse is stored.
		*/
		void inv ( ldmat &Inv ) const;

		/*! \brief Symmetric multiplication of \f$U\f$ by a general matrix \f$C\f$, result of which is stored in the current class.
		@param C matrix to multiply with
		@param U a space where the inverse is stored.
		*/
		void mult_sym ( const mat &C, ldmat &U) const;

		/*! \brief Symmetric multiplication of \f$U\f$ by a transpose of a general matrix \f$C\f$, result of which is stored in the current class.
		@param C matrix to multiply with
		@param U a space where the inverse is stored.
		*/
		void mult_sym_t ( const mat &C, ldmat &U) const;


		/*! \brief Transforms general \f$A'D0 A\f$ into pure \f$L'DL\f$

		The new decomposition fullfills: \f$A'*diag(D)*A = self.L'*diag(self.D)*self.L\f$
		@param A general matrix
		@param D0 general vector
		*/
		void ldform (const mat &A,const vec &D0 );

		//! Access functions
		void setD (const vec &nD){D=nD;}
		//! Access functions
		void setD (const vec &nD, int i){D.replace_mid(i,nD);} //Fixme can be more general
		//! Access functions
		void setL (const vec &nL){L=nL;}

		//! Access functions
		const vec& _D() const {return D;}
		//! Access functions
		const mat& _L() const {return L;}

		//! add another ldmat matrix
		ldmat& operator += ( const ldmat &ldA );
		//! subtract another ldmat matrix
		ldmat& operator -= ( const ldmat &ldA );
		//! multiply by a scalar
		ldmat& operator *= ( double x );

		//! print both \c L and \c D
		friend std::ostream &operator<< ( std::ostream &os, const ldmat &sq );

	protected:
		//! Positive vector \f$D\f$
		vec D;
		//! Lower-triangular matrix \f$L\f$
		mat L;

};

//////// Operations:
//!mapping of add operation to operators
inline ldmat& ldmat::operator += ( const ldmat &ldA )  {this->add ( ldA );return *this;}
//!mapping of negative add operation to operators
inline ldmat& ldmat::operator -= ( const ldmat &ldA )  {this->add ( ldA,-1.0 );return *this;}

/*! @} */

#endif // DC_H