libEF.h

Go to the documentation of this file.

#ifndef EF_H
#define EF_H


#include "libBM.h"
#include "../math/chmat.h"
//#include "../user_info.h"

namespace bdm
{


        extern Uniform_RNG UniRNG;
        extern Normal_RNG NorRNG;
        extern Gamma_RNG GamRNG;

        class eEF : public epdf
        {
                public:
//      eEF() :epdf() {};
                        eEF ( ) :epdf ( ) {};
                        virtual double lognc() const =0;
                        virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
                        virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
                        virtual double evallog ( const vec &val ) const {double tmp;tmp= evallog_nn ( val )-lognc();it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); return tmp;}
                        virtual vec evallog ( const mat &Val ) const
                        {
                                vec x ( Val.cols() );
                                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog_nn ( Val.get_col ( i ) ) ;}
                                return x-lognc();
                        }
                        virtual void pow ( double p ) {it_error ( "Not implemented" );};
        };

        class mEF : public mpdf
        {

                public:
                        mEF ( ) :mpdf ( ) {};
        };

        class BMEF : public BM
        {
                protected:
                        double frg;
                        double last_lognc;
                public:
                        BMEF ( double frg0=1.0 ) :BM ( ), frg ( frg0 ) {}
                        BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
                        virtual void set_statistics ( const BMEF* BM0 ) {it_error ( "Not implemented" );};
                        virtual void bayes ( const vec &data, const double w ) {};
                        //original Bayes
                        void bayes ( const vec &dt );
                        virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
//      virtual void flatten ( double nu0 ) {it_error ( "Not implemented" );}

                        BMEF* _copy_ () const {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
        };

        template<class sq_T>
        class mlnorm;

        template<class sq_T>
        class enorm : public eEF
        {
                protected:
                        vec mu;
                        sq_T R;
                public:

                        enorm ( ) :eEF ( ), mu ( ),R ( ) {};
                        enorm ( const vec &mu,const sq_T &R ) {set_parameters ( mu,R );}
                        void set_parameters ( const vec &mu,const sq_T &R );


                        void dupdate ( mat &v,double nu=1.0 );

                        vec sample() const;
                        mat sample ( int N ) const;
                        double evallog_nn ( const vec &val ) const;
                        double lognc () const;
                        vec mean() const {return mu;}
                        vec variance() const {return diag ( R.to_mat() );}
//      mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
                        mpdf* condition ( const RV &rvn ) const ;
        enorm<sq_T>* marginal ( const RV &rv ) const;
//                      epdf* marginal ( const RV &rv ) const;


                        vec& _mu() {return mu;}
                        void set_mu ( const vec mu0 ) { mu=mu0;}
                        sq_T& _R() {return R;}
                        const sq_T& _R() const {return R;}

        };

        class egiw : public eEF
        {
                protected:
                        ldmat V;
                        double nu;
                        int dimx;
                        int nPsi;
                public:
                        egiw() :eEF() {};
                        egiw ( int dimx0, ldmat V0, double nu0=-1.0 ) :eEF() {set_parameters ( dimx0,V0, nu0 );};

                        void set_parameters ( int dimx0, ldmat V0, double nu0=-1.0 )
                        {
                                dimx=dimx0;
                                nPsi = V0.rows()-dimx;
                                dim = dimx* ( dimx+nPsi ); // size(R) + size(Theta)

                                V=V0;
                                if ( nu0<0 )
                                {
                                        nu = 0.1 +nPsi +2*dimx +2; // +2 assures finite expected value of R
                                        // terms before that are sufficient for finite normalization
                                }
                                else
                                {
                                        nu=nu0;
                                }
                        }

                        vec sample() const;
                        vec mean() const;
                        vec variance() const;

                        vec est_theta() const;

                        ldmat est_theta_cov() const;

                        void mean_mat ( mat &M, mat&R ) const;
                        double evallog_nn ( const vec &val ) const;
                        double lognc () const;
                        void pow ( double p ) {V*=p;nu*=p;};


                        ldmat& _V() {return V;}
                        const ldmat& _V() const {return V;}
                        double& _nu()  {return nu;}
                        const double& _nu() const {return nu;}
        };

        class eDirich: public eEF
        {
                protected:
                        vec beta;
                public:

                        eDirich () : eEF ( ) {};
                        eDirich ( const eDirich &D0 ) : eEF () {set_parameters ( D0.beta );};
                        eDirich ( const vec &beta0 ) {set_parameters ( beta0 );};
                        void set_parameters ( const vec &beta0 )
                        {
                                beta= beta0;
                                dim = beta.length();
                        }

                        vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
                        vec mean() const {return beta/sum(beta);};
                        vec variance() const {double gamma =sum(beta); return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
                        double evallog_nn ( const vec &val ) const
                        {
                                double tmp; tmp= ( beta-1 ) *log ( val );               it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
                                return tmp;
                        };
                        double lognc () const
                        {
                                double tmp;
                                double gam=sum ( beta );
                                double lgb=0.0;
                                for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
                                tmp= lgb-lgamma ( gam );
                                it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
                                return tmp;
                        };
                        vec& _beta()  {return beta;}
        };

        class multiBM : public BMEF
        {
                protected:
                        eDirich est;
                        vec &beta;
                public:
                        multiBM ( ) : BMEF ( ),est ( ),beta ( est._beta() )
                        {
                                if ( beta.length() >0 ) {last_lognc=est.lognc();}
                                else{last_lognc=0.0;}
                        }
                        multiBM ( const multiBM &B ) : BMEF ( B ),est ( B.est ),beta ( est._beta() ) {}
                        void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
                        void bayes ( const vec &dt )
                        {
                                if ( frg<1.0 ) {beta*=frg;last_lognc=est.lognc();}
                                beta+=dt;
                                if ( evalll ) {ll=est.lognc()-last_lognc;}
                        }
                        double logpred ( const vec &dt ) const
                        {
                                eDirich pred ( est );
                                vec &beta = pred._beta();

                                double lll;
                                if ( frg<1.0 )
                                        {beta*=frg;lll=pred.lognc();}
                                else
                                        if ( evalll ) {lll=last_lognc;}
                                        else{lll=pred.lognc();}

                                beta+=dt;
                                return pred.lognc()-lll;
                        }
                        void flatten ( const BMEF* B )
                        {
                                const multiBM* E=dynamic_cast<const multiBM*> ( B );
                                // sum(beta) should be equal to sum(B.beta)
                                const vec &Eb=E->beta;//const_cast<multiBM*> ( E )->_beta();
                                beta*= ( sum ( Eb ) /sum ( beta ) );
                                if ( evalll ) {last_lognc=est.lognc();}
                        }
                        const epdf& posterior() const {return est;};
                        const eDirich* _e() const {return &est;};
                        void set_parameters ( const vec &beta0 )
                        {
                                est.set_parameters ( beta0 );
                                if ( evalll ) {last_lognc=est.lognc();}
                        }
        };

        class egamma : public eEF
        {
                protected:
                        vec alpha;
                        vec beta;
                public :
                        egamma ( ) :eEF ( ), alpha ( 0 ), beta ( 0 ) {};
                        egamma ( const vec &a, const vec &b ) {set_parameters ( a, b );};
                        void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;dim = alpha.length();};

                        vec sample() const;
//      mat sample ( int N ) const;
                        double evallog ( const vec &val ) const;
                        double lognc () const;
                        vec& _alpha() {return alpha;}
                        vec& _beta() {return beta;}
                        vec mean() const {return elem_div ( alpha,beta );}
                        vec variance() const {return elem_div ( alpha,elem_mult ( beta,beta ) ); }
        };

        class eigamma : public egamma
        {
                protected:
                public :

                        vec sample() const {return 1.0/egamma::sample();};
                        vec mean() const {return elem_div ( beta,alpha-1 );}
                        vec variance() const {vec mea=mean(); return elem_div ( elem_mult ( mea,mea ),alpha-2 );}
        };
        /*
        class emix : public epdf {
        protected:
                int n;
                vec &w;
                Array<epdf*> Coms;
        public:
                emix ( const RV &rv, vec &w0): epdf(rv), n(w0.length()), w(w0), Coms(n) {};
                void set_parameters( int &i, double wi, epdf* ep){w(i)=wi;Coms(i)=ep;}
                vec mean(){vec pom; for(int i=0;i<n;i++){pom+=Coms(i)->mean()*w(i);} return pom;};
                vec sample() {it_error ( "Not implemented" );return 0;}
        };
        */


        class euni: public epdf
        {
                protected:
                        vec low;
                        vec high;
                        vec distance;
                        double nk;
                        double lnk;
                public:
                        euni ( ) :epdf ( ) {}
                        euni ( const vec &low0, const vec &high0 ) {set_parameters ( low0,high0 );}
                        void set_parameters ( const vec &low0, const vec &high0 )
                        {
                                distance = high0-low0;
                                it_assert_debug ( min ( distance ) >0.0,"bad support" );
                                low = low0;
                                high = high0;
                                nk = prod ( 1.0/distance );
                                lnk = log ( nk );
                                dim = low.length();
                        }

                        double eval ( const vec &val ) const  {return nk;}
                        double evallog ( const vec &val ) const  {return lnk;}
                        vec sample() const
                        {
                                vec smp ( dim );
#pragma omp critical
                                UniRNG.sample_vector ( dim ,smp );
                                return low+elem_mult ( distance,smp );
                        }
                        vec mean() const {return ( high-low ) /2.0;}
                        vec variance() const {return ( pow ( high,2 ) +pow ( low,2 ) +elem_mult ( high,low ) ) /3.0;}
        };


        template<class sq_T>
        class mlnorm : public mEF
        {
                protected:
                        enorm<sq_T> epdf;
                        mat A;
                        vec mu_const;
                        vec& _mu; //cached epdf.mu;
                public:
                        mlnorm ( ) :mEF (),epdf ( ),A ( ),_mu ( epdf._mu() ) {ep =&epdf; };
                        mlnorm ( const  mat &A, const vec &mu0, const sq_T &R ) :epdf ( ),_mu ( epdf._mu() )
                        {
                                ep =&epdf; set_parameters ( A,mu0,R );
                        };
                        void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
                        void condition ( const vec &cond );

                        vec& _mu_const() {return mu_const;}
                        mat& _A() {return A;}
                        mat _R() {return epdf._R().to_mat();}

                        template<class sq_M>
                        friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M> &ml );
        };

        template<class sq_T>
        class mgnorm : public mEF
        {
                protected:
                        enorm<sq_T> epdf;
                        vec &mu;
                        fnc* g;
                public:
                        mgnorm() :mu ( epdf._mu() ) {ep=&epdf;}
                        void set_parameters ( fnc* g0, const sq_T &R0 ) {g=g0; epdf.set_parameters ( zeros ( g->dimension() ), R0 );}
                        void condition ( const vec &cond ) {mu=g->eval ( cond );};


                        void from_setting( const Setting &root ) 
                        {       
                                fnc* g = UI::build<fnc>( root, "g" );

                                mat R;
                                if ( root.exists( "dR" ) )
                                {
                                        vec dR;
                                        UI::get( dR, root, "dR" );
                                        R=diag(dR);
                                }
                                else 
                                        UI::get( R, root, "R");                                 
                
                                set_parameters(g,R);
                        }

                        /*void mgnorm::to_setting( Setting &root ) const
                        {       
                                Transport::to_setting( root );

                                Setting &kilometers_setting = root.add("kilometers", Setting::TypeInt );
                                kilometers_setting = kilometers;

                                UI::save( passengers, root, "passengers" );
                        }*/

        };

        UIREGISTER(mgnorm<chmat>);


        class mlstudent : public mlnorm<ldmat>
        {
                protected:
                        ldmat Lambda;
                        ldmat &_R;
                        ldmat Re;
                public:
                        mlstudent ( ) :mlnorm<ldmat> (),
                                        Lambda (),      _R ( epdf._R() ) {}
                        void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 )
                        {
                                it_assert_debug ( A0.rows() ==mu0.length(),"" );
                                it_assert_debug ( R0.rows() ==A0.rows(),"" );

                                epdf.set_parameters ( mu0,Lambda ); //
                                A = A0;
                                mu_const = mu0;
                                Re=R0;
                                Lambda = Lambda0;
                        }
                        void condition ( const vec &cond )
                        {
                                _mu = A*cond + mu_const;
                                double zeta;
                                //ugly hack!
                                if ( ( cond.length() +1 ) ==Lambda.rows() )
                                {
                                        zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
                                }
                                else
                                {
                                        zeta = Lambda.invqform ( cond );
                                }
                                _R = Re;
                                _R*= ( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
                        };

        };
        class mgamma : public mEF
        {
                protected:
                        egamma epdf;
                        double k;
                        vec &_beta;

                public:
                        mgamma ( ) : mEF ( ), epdf (), _beta ( epdf._beta() ) {ep=&epdf;};
                        void set_parameters ( double k, const vec &beta0 );
                        void condition ( const vec &val ) {_beta=k/val;};
        };

        class migamma : public mEF
        {
                protected:
                        eigamma epdf;
                        double k;
                        vec &_alpha;
                        vec &_beta;

                public:
                        migamma ( ) : mEF (), epdf ( ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
                        migamma ( const migamma &m ) : mEF (), epdf ( m.epdf ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};

                        void set_parameters ( int len, double k0 )
                        {
                                k=k0;
                                epdf.set_parameters ( ( 1.0/ ( k*k ) +2.0 ) *ones ( len ) /*alpha*/, ones ( len ) /*beta*/ );
                                dimc = dimension();
                        };
                        void condition ( const vec &val )
                        {
                                _beta=elem_mult ( val, ( _alpha-1.0 ) );
                        };
        };


        class mgamma_fix : public mgamma
        {
                protected:
                        double l;
                        vec refl;
                public:
                        mgamma_fix ( ) : mgamma ( ),refl () {};
                        void set_parameters ( double k0 , vec ref0, double l0 )
                        {
                                mgamma::set_parameters ( k0, ref0 );
                                refl=pow ( ref0,1.0-l0 );l=l0;
                                dimc=dimension();
                        };

                        void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); _beta=k/mean;};
        };


        class migamma_ref : public migamma
        {
                protected:
                        double l;
                        vec refl;
                public:
                        migamma_ref ( ) : migamma (),refl ( ) {};
                        void set_parameters ( double k0 , vec ref0, double l0 )
                        {
                                migamma::set_parameters ( ref0.length(), k0 );
                                refl=pow ( ref0,1.0-l0 );
                                l=l0;
                                dimc = dimension();
                        };

                        void condition ( const vec &val )
                        {
                                vec mean=elem_mult ( refl,pow ( val,l ) );
                                migamma::condition ( mean );
                        };

                        void from_setting( const Setting &root );

                        // TODO dodelat void to_setting( Setting &root ) const;
        };


        UIREGISTER(migamma_ref);

        class elognorm: public enorm<ldmat>
        {
                public:
                        vec sample() const {return exp ( enorm<ldmat>::sample() );};
                        vec mean() const {vec var=enorm<ldmat>::variance();return exp ( mu - 0.5*var );};

        };

        class mlognorm : public mpdf
        {
                protected:
                        elognorm eno;
                        double sig2;
                        vec &mu;
                public:
                        mlognorm ( ) : eno (), mu ( eno._mu() ) {ep=&eno;};
                        void set_parameters ( int size, double k )
                        {
                                sig2 = 0.5*log ( k*k+1 );
                                eno.set_parameters ( zeros ( size ),2*sig2*eye ( size ) );

                                dimc = size;
                        };

                        void condition ( const vec &val )
                        {
                                mu=log ( val )-sig2;//elem_mult ( refl,pow ( val,l ) );
                        };

                        void from_setting( const Setting &root );

                        // TODO dodelat void to_setting( Setting &root ) const;

        };
        
        UIREGISTER(mlognorm);

        class eWishartCh : public epdf
        {
                protected:
                        chmat Y;
                        int p;
                        double delta;
                public:
                        void set_parameters ( const mat &Y0, const double delta0 ) {Y=chmat ( Y0 );delta=delta0; p=Y.rows(); dim = p*p; }
                        mat sample_mat() const
                        {
                                mat X=zeros ( p,p );

                                //sample diagonal
                                for ( int i=0;i<p;i++ )
                                {
                                        GamRNG.setup ( 0.5* ( delta-i ) , 0.5 ); // no +1 !! index if from 0
#pragma omp critical
                                        X ( i,i ) =sqrt ( GamRNG() );
                                }
                                //do the rest
                                for ( int i=0;i<p;i++ )
                                {
                                        for ( int j=i+1;j<p;j++ )
                                        {
#pragma omp critical
                                                X ( i,j ) =NorRNG.sample();
                                        }
                                }
                                return X*Y._Ch();// return upper triangular part of the decomposition
                        }
                        vec sample () const
                        {
                                return vec ( sample_mat()._data(),p*p );
                        }
                        void setY ( const mat &Ch0 ) {copy_vector ( dim,Ch0._data(), Y._Ch()._data() );}
                        void _setY ( const vec &ch0 ) {copy_vector ( dim, ch0._data(), Y._Ch()._data() ); }
                        const chmat& getY()const {return Y;}
        };

        class eiWishartCh: public epdf
        {
                protected:
                        eWishartCh W;
                        int p;
                        double delta;
                public:
                        void set_parameters ( const mat &Y0, const double delta0) {
                                delta = delta0;
                                W.set_parameters ( inv ( Y0 ),delta0 ); 
                                dim = W.dimension(); p=Y0.rows();
                        }
                        vec sample() const {mat iCh; iCh=inv ( W.sample_mat() ); return vec ( iCh._data(),dim );}
                        void _setY ( const vec &y0 )
                        {
                                mat Ch ( p,p );
                                mat iCh ( p,p );
                                copy_vector ( dim, y0._data(), Ch._data() );
                                
                                iCh=inv ( Ch );
                                W.setY ( iCh );
                        }
                        virtual double evallog ( const vec &val ) const {
                                chmat X(p);
                                const chmat& Y=W.getY();
                                 
                                copy_vector(p*p,val._data(),X._Ch()._data());
                                chmat iX(p);X.inv(iX);
                                // compute  
//                              \frac{ |\Psi|^{m/2}|X|^{-(m+p+1)/2}e^{-tr(\Psi X^{-1})/2} }{ 2^{mp/2}\Gamma_p(m/2)},
                                mat M=Y.to_mat()*iX.to_mat();
                                
                                double log1 = 0.5*p*(2*Y.logdet())-0.5*(delta+p+1)*(2*X.logdet())-0.5*trace(M); 
                                //Fixme! Multivariate gamma omitted!! it is ok for sampling, but not otherwise!!
                                
/*                              if (0) {
                                        mat XX=X.to_mat();
                                        mat YY=Y.to_mat();
                                        
                                        double log2 = 0.5*p*log(det(YY))-0.5*(delta+p+1)*log(det(XX))-0.5*trace(YY*inv(XX)); 
                                        cout << log1 << "," << log2 << endl;
                                }*/
                                return log1;                            
                        };
                        
        };

        class rwiWishartCh : public mpdf
        {
                protected:
                        eiWishartCh eiW;
                        double sqd;
                        //reference point for diagonal
                        vec refl;
                        double l;
                        int p;
                public:
                        void set_parameters ( int p0, double k, vec ref0, double l0  )
                        {
                                p=p0;
                                double delta = 2/(k*k)+p+3;
                                sqd=sqrt ( delta-p-1 );
                                l=l0;
                                refl=pow(ref0,1-l);
                                
                                eiW.set_parameters ( eye ( p ),delta );
                                ep=&eiW;
                                dimc=eiW.dimension();
                        }
                        void condition ( const vec &c ) {
                                vec z=c;
                                int ri=0;
                                for(int i=0;i<p*p;i+=(p+1)){//trace diagonal element
                                        z(i) = pow(z(i),l)*refl(ri);
                                        ri++;
                                }

                                eiW._setY ( sqd*z );
                        }
        };

        enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
        class eEmp: public epdf
        {
                protected :
                        int n;
                        vec w;
                        Array<vec> samples;
                public:
                        eEmp ( ) :epdf ( ),w ( ),samples ( ) {};
                        eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};

                        void set_statistics ( const vec &w0, const epdf* pdf0 );
                        void set_statistics ( const epdf* pdf0 , int n ) {set_statistics ( ones ( n ) /n,pdf0 );};
                        void set_samples ( const epdf* pdf0 );
                        void set_parameters ( int n0, bool copy=true ) {n=n0; w.set_size ( n0,copy );samples.set_size ( n0,copy );};
                        vec& _w()  {return w;};
                        const vec& _w() const {return w;};
                        Array<vec>& _samples() {return samples;};
                        const Array<vec>& _samples() const {return samples;};
                        ivec resample ( RESAMPLING_METHOD method=SYSTEMATIC );
                        vec sample() const {it_error ( "Not implemented" );return 0;}
                        double evallog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
                        vec mean() const
                        {
                                vec pom=zeros ( dim );
                                for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
                                return pom;
                        }
                        vec variance() const
                        {
                                vec pom=zeros ( dim );
                                for ( int i=0;i<n;i++ ) {pom+=pow ( samples ( i ),2 ) *w ( i );}
                                return pom-pow ( mean(),2 );
                        }
                        void qbounds ( vec &lb, vec &ub, double perc=0.95 ) const
                        {
                                // lb in inf so than it will be pushed below;
                                lb.set_size ( dim );
                                ub.set_size ( dim );
                                lb = std::numeric_limits<double>::infinity();
                                ub = -std::numeric_limits<double>::infinity();
                                int j;
                                for ( int i=0;i<n;i++ )
                                {
                                        for ( j=0;j<dim; j++ )
                                        {
                                                if ( samples ( i ) ( j ) <lb ( j ) ) {lb ( j ) =samples ( i ) ( j );}
                                                if ( samples ( i ) ( j ) >ub ( j ) ) {ub ( j ) =samples ( i ) ( j );}
                                        }
                                }
                        }
        };



        template<class sq_T>
        void enorm<sq_T>::set_parameters ( const vec &mu0, const sq_T &R0 )
        {
//Fixme test dimensions of mu0 and R0;
                mu = mu0;
                R = R0;
                dim = mu0.length();
        };

        template<class sq_T>
        void enorm<sq_T>::dupdate ( mat &v, double nu )
        {
                //
        };

// template<class sq_T>
// void enorm<sq_T>::tupdate ( double phi, mat &vbar, double nubar ) {
//      //
// };

        template<class sq_T>
        vec enorm<sq_T>::sample() const
        {
                vec x ( dim );
#pragma omp critical
                NorRNG.sample_vector ( dim,x );
                vec smp = R.sqrt_mult ( x );

                smp += mu;
                return smp;
        };

        template<class sq_T>
        mat enorm<sq_T>::sample ( int N ) const
        {
                mat X ( dim,N );
                vec x ( dim );
                vec pom;
                int i;

                for ( i=0;i<N;i++ )
                {
#pragma omp critical
                        NorRNG.sample_vector ( dim,x );
                        pom = R.sqrt_mult ( x );
                        pom +=mu;
                        X.set_col ( i, pom );
                }

                return X;
        };

// template<class sq_T>
// double enorm<sq_T>::eval ( const vec &val ) const {
//      double pdfl,e;
//      pdfl = evallog ( val );
//      e = exp ( pdfl );
//      return e;
// };

        template<class sq_T>
        double enorm<sq_T>::evallog_nn ( const vec &val ) const
        {
                // 1.83787706640935 = log(2pi)
                double tmp=-0.5* ( R.invqform ( mu-val ) );// - lognc();
                return  tmp;
        };

        template<class sq_T>
        inline double enorm<sq_T>::lognc () const
        {
                // 1.83787706640935 = log(2pi)
                double tmp=0.5* ( R.cols() * 1.83787706640935 +R.logdet() );
                return tmp;
        };

        template<class sq_T>
        void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 )
        {
                it_assert_debug ( A0.rows() ==mu0.length(),"" );
                it_assert_debug ( A0.rows() ==R0.rows(),"" );

                epdf.set_parameters ( zeros ( A0.rows() ),R0 );
                A = A0;
                mu_const = mu0;
                dimc=A0.cols();
        }

// template<class sq_T>
// vec mlnorm<sq_T>::samplecond (const  vec &cond, double &lik ) {
//      this->condition ( cond );
//      vec smp = epdf.sample();
//      lik = epdf.eval ( smp );
//      return smp;
// }

// template<class sq_T>
// mat mlnorm<sq_T>::samplecond (const vec &cond, vec &lik, int n ) {
//      int i;
//      int dim = rv.count();
//      mat Smp ( dim,n );
//      vec smp ( dim );
//      this->condition ( cond );
//
//      for ( i=0; i<n; i++ ) {
//              smp = epdf.sample();
//              lik ( i ) = epdf.eval ( smp );
//              Smp.set_col ( i ,smp );
//      }
//
//      return Smp;
// }

        template<class sq_T>
        void mlnorm<sq_T>::condition ( const vec &cond )
        {
                _mu = A*cond + mu_const;
//R is already assigned;
        }

        template<class sq_T>
        enorm<sq_T>* enorm<sq_T>::marginal ( const RV &rvn ) const
        {
                it_assert_debug ( isnamed(), "rv description is not assigned" );
                ivec irvn = rvn.dataind ( rv );

                sq_T Rn ( R,irvn ); //select rows and columns of R

                enorm<sq_T>* tmp = new enorm<sq_T>;
                tmp->set_rv ( rvn );
                tmp->set_parameters ( mu ( irvn ), Rn );
                return tmp;
        }

        template<class sq_T>
        mpdf* enorm<sq_T>::condition ( const RV &rvn ) const
        {

                it_assert_debug ( isnamed(),"rvs are not assigned" );

                RV rvc = rv.subt ( rvn );
                it_assert_debug ( ( rvc._dsize() +rvn._dsize() ==rv._dsize() ),"wrong rvn" );
                //Permutation vector of the new R
                ivec irvn = rvn.dataind ( rv );
                ivec irvc = rvc.dataind ( rv );
                ivec perm=concat ( irvn , irvc );
                sq_T Rn ( R,perm );

                //fixme - could this be done in general for all sq_T?
                mat S=Rn.to_mat();
                //fixme
                int n=rvn._dsize()-1;
                int end=R.rows()-1;
                mat S11 = S.get ( 0,n, 0, n );
                mat S12 = S.get ( 0, n , rvn._dsize(), end );
                mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );

                vec mu1 = mu ( irvn );
                vec mu2 = mu ( irvc );
                mat A=S12*inv ( S22 );
                sq_T R_n ( S11 - A *S12.T() );

                mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( );
                tmp->set_rv ( rvn ); tmp->set_rvc ( rvc );
                tmp->set_parameters ( A,mu1-A*mu2,R_n );
                return tmp;
        }


        template<class sq_T>
        std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_T> &ml )
        {
                os << "A:"<< ml.A<<endl;
                os << "mu:"<< ml.mu_const<<endl;
                os << "R:" << ml.epdf._R().to_mat() <<endl;
                return os;
        };

}
#endif //EF_H

---