root/library/bdm/math/square_mat.h @ 495

/*!
 * \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"

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