/*!
 * \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"
#include "../bdmerror.h"

namespace bdm
{

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 );

	/*!
	  Some templates require this constructor to compile, but
	  it shouldn't actually be called.
	*/
	fsqmat ( const fsqmat &M, const ivec &perm ) {
		bdm_error ( "not implemented" );
	}

	//! 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;
	};
	
	//! cast to normal mat
	operator mat&() {return M;};
	
//		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