root/library/bdm/stat/exp_family.h @ 956

/*!
  \file
  \brief Probability distributions for Exponential Family models.
  \author Vaclav Smidl.

  -----------------------------------
  BDM++ - C++ library for Bayesian Decision Making under Uncertainty

  Using IT++ for numerical operations
  -----------------------------------
*/

#ifndef EF_H
#define EF_H


#include "../shared_ptr.h"
#include "../base/bdmbase.h"
#include "../math/chmat.h"

namespace bdm {


//! Global Uniform_RNG
extern Uniform_RNG UniRNG;
//! Global Normal_RNG
extern Normal_RNG NorRNG;
//! Global Gamma_RNG
extern Gamma_RNG GamRNG;

/*!
* \brief General conjugate exponential family posterior density.

* More?...
*/

class eEF : public epdf {
public:
//      eEF() :epdf() {};
        //! default constructor
        eEF () : epdf () {};
        //! logarithm of the normalizing constant, \f$\mathcal{I}\f$
        virtual double lognc() const = 0;

        //!Evaluate normalized log-probability
        virtual double evallog_nn ( const vec &val ) const NOT_IMPLEMENTED(0);

        //!Evaluate normalized log-probability
        virtual double evallog ( const vec &val ) const {
                double tmp;
                tmp = evallog_nn ( val ) - lognc();
                return tmp;
        }
        //!Evaluate normalized log-probability for many samples
        virtual vec evallog_mat ( 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();
        }
        //!Evaluate normalized log-probability for many samples
        virtual vec evallog_mat ( const Array<vec> &Val ) const {
                vec x ( Val.length() );
                for ( int i = 0; i < Val.length(); i++ ) {
                        x ( i ) = evallog_nn ( Val ( i ) ) ;
                }
                return x - lognc();
        }

        //!Power of the density, used e.g. to flatten the density
        virtual void pow ( double p ) NOT_IMPLEMENTED_VOID;
};


//! Estimator for Exponential family
class BMEF : public BM {
protected:
        //! forgetting factor
        double frg;
        //! cached value of lognc() in the previous step (used in evaluation of \c ll )
        double last_lognc;
public:
        //! Default constructor (=empty constructor)
        BMEF ( double frg0 = 1.0 ) : BM (), frg ( frg0 ) {}
        //! Copy constructor
        BMEF ( const BMEF &B ) : BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
        //!get statistics from another model
        virtual void set_statistics ( const BMEF* BM0 ) NOT_IMPLEMENTED_VOID;

        //! Weighted update of sufficient statistics (Bayes rule)
        virtual void bayes_weighted ( const vec &data, const vec &cond = empty_vec, const double w = 1.0 ) {};
        //original Bayes
        void bayes ( const vec &yt, const vec &cond = empty_vec );

        //!Flatten the posterior according to the given BMEF (of the same type!)
        virtual void flatten ( const BMEF * B ) NOT_IMPLEMENTED_VOID;


        void to_setting ( Setting &set ) const
        {                       
                BM::to_setting( set );
                UI::save(frg, set, "frg");
        } 

        void from_setting( const Setting &set) {
                BM::from_setting(set);
                if ( !UI::get ( frg, set, "frg" ) )
                        frg = 1.0;
        }

        void validate() {
                BM::validate(); 
        }

};

/*! Dirac delta density with predefined transformation

Density of the type:\f[ f(x_t | y_t) = \delta (x_t - g(y_t)) \f]
where \f$ x_t \f$ is the \c rv, \f$ y_t \f$ is the \c rvc and g is a deterministic transformation of class fn.
*/
class mgdirac: public pdf{
        protected:
        shared_ptr<fnc> g;
        public:
                vec samplecond(const vec &cond) {
                        bdm_assert_debug(cond.length()==g->dimensionc(),"given cond in not compatible with g");
                        vec tmp = g->eval(cond);
                        return tmp;
                } 
                double evallogcond ( const vec &yt, const vec &cond ){
                        return std::numeric_limits< double >::max();
                }
                void from_setting(const Setting& set);
                void to_setting(Setting &set) const;
                void validate();
};
UIREGISTER(mgdirac);


template<class sq_T, template <typename> class TEpdf>
class mlnorm;

/*!
* \brief Gaussian density with positive definite (decomposed) covariance matrix.

* More?...
*/
template<class sq_T>
class enorm : public eEF {
protected:
        //! mean value
        vec mu;
        //! Covariance matrix in decomposed form
        sq_T R;
public:
        //!\name Constructors
        //!@{

        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 );
        /*! Create Normal density
        \f[ f(rv) = N(\mu, R) \f]
        from structure
        \code
        class = 'enorm<ldmat>', (OR) 'enorm<chmat>', (OR) 'enorm<fsqmat>';
        mu    = [];                  // mean value
        R     = [];                  // variance, square matrix of appropriate dimension
        \endcode
        */
        void from_setting ( const Setting &root );
        void to_setting ( Setting &root ) const ;
        
        void validate();
        //!@}

        //! \name Mathematical operations
        //!@{

        //! dupdate in exponential form (not really handy)
        void dupdate ( mat &v, double nu = 1.0 );

        //! evaluate bhattacharya distance
        double bhattacharyya(const enorm<sq_T> &e2){
                bdm_assert(dim == e2.dimension(), "enorms of differnt dimensions");
                sq_T P=R;
                P.add(e2._R());
                
                double tmp = 0.125*P.invqform(mu - e2._mu()) + 0.5*(P.logdet() - 0.5*(R.logdet() + e2._R().logdet()));
                return tmp;
        }
        
        vec sample() const;

        double evallog_nn ( const vec &val ) const;
        double lognc () const;
        vec mean() const {
                return mu;
        }
        vec variance() const {
                return diag ( R.to_mat() );
        }
        mat covariance() const {
                return R.to_mat();
        }
        //      mlnorm<sq_T>* condition ( const RV &rvn ) const ; <=========== fails to cmpile. Why?
        shared_ptr<pdf> condition ( const RV &rvn ) const;

        // target not typed to mlnorm<sq_T, enorm<sq_T> > &
        // because that doesn't compile (perhaps because we
        // haven't finished defining enorm yet), but the type
        // is required
        void condition ( const RV &rvn, pdf &target ) const;

        shared_ptr<epdf> marginal ( const RV &rvn ) const;
        void marginal ( const RV &rvn, enorm<sq_T> &target ) const;
        //!@}

        //! \name Access to attributes
        //!@{

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

};
UIREGISTER2 ( enorm, chmat );
SHAREDPTR2 ( enorm, chmat );
UIREGISTER2 ( enorm, ldmat );
SHAREDPTR2 ( enorm, ldmat );
UIREGISTER2 ( enorm, fsqmat );
SHAREDPTR2 ( enorm, fsqmat );

//! \class bdm::egauss 
//!\brief Gaussian (Normal) distribution. Same as enorm<fsqmat>.
typedef enorm<ldmat> egauss;
UIREGISTER(egauss);


//forward declaration
class mstudent;

/*! distribution of multivariate Student t density

Based on article by Genest and Zidek, 
*/
template<class sq_T>
class estudent : public eEF{
        protected:
                //! mena value
                vec mu;
                //! matrix H
                sq_T H;
                //! degrees of freedom
                double delta;
        public:
                double evallog_nn(const vec &val) const{
                        double tmp = -0.5*H.logdet() - 0.5*(delta + dim) * log(1+ H.invqform(val - mu)/delta);
                        return tmp;
                }
                double lognc() const {
                        //log(pi) = 1.14472988584940
                        double tmp = -lgamma(0.5*(delta+dim))+lgamma(0.5*delta) + 0.5*dim*(log(delta) + 1.14472988584940);
                        return tmp;
                }
                void marginal (const RV &rvm, estudent<sq_T> &marg) const {
                        ivec ind = rvm.findself_ids(rv); // indices of rvm in rv
                        marg._mu() = mu(ind);
                        marg._H() = sq_T(H,ind);
                        marg._delta() = delta;
                        marg.validate();
                }
                shared_ptr<epdf> marginal(const RV &rvm) const {
                        shared_ptr<estudent<sq_T> > tmp = new estudent<sq_T>;
                        marginal(rvm, *tmp);
                        return tmp;
                }
                vec sample() const NOT_IMPLEMENTED(vec(0))
                
                vec mean() const {return mu;}
                mat covariance() const {return delta/(delta-2)*H.to_mat();}
                vec variance() const {return diag(covariance());}
                        //! \name access
                //! @{
                        //! access function
                        vec& _mu() {return mu;}
                        //! access function
                        sq_T& _H() {return H;}
                        //! access function
                        double& _delta() {return delta;}
                        //!@}
                        //! todo
                        void from_setting(const Setting &set){
                                epdf::from_setting(set);
                                mat H0;
                                UI::get(H0,set, "H");
                                H= H0; // conversion!!
                                UI::get(delta,set,"delta");
                                UI::get(mu,set,"mu");
                        }
                        void to_setting(Setting &set) const{
                                epdf::to_setting(set);
                                UI::save(H.to_mat(), set, "H");
                                UI::save(delta, set, "delta");
                                UI::save(mu, set, "mu");
                        }
                        void validate() {
                                eEF::validate();
                                dim = H.rows();
                        }
};
UIREGISTER2(estudent,fsqmat);
UIREGISTER2(estudent,ldmat);
UIREGISTER2(estudent,chmat);

/*!
* \brief Gauss-inverse-Wishart density stored in LD form

* For \f$p\f$-variate densities, given rv.count() should be \f$p\times\f$ V.rows().
*
*/
class egiw : public eEF {
        //! \var log_level_enums logvariance
        //! TODO DOPLNIT

        //! \var log_level_enums logmean
        //! TODO DOPLNIT
        LOG_LEVEL(egiw,logmean, logvariance);

protected:
        //! Extended information matrix of sufficient statistics
        ldmat V;
        //! Number of data records (degrees of freedom) of sufficient statistics
        double nu;
        //! Dimension of the output
        int dimx;
        //! Dimension of the regressor
        int nPsi;
public:
        //!\name Constructors
        //!@{
        egiw() : eEF() {};
        egiw ( int dimx0, ldmat V0, double nu0 = -1.0 ) : eEF() {
                set_parameters ( dimx0, V0, nu0 );
                validate();
        };

        void set_parameters ( int dimx0, ldmat V0, double nu0 = -1.0 );
        //!@}

        vec sample() const;
        mat sample_mat ( int n ) const;
        vec mean() const;
        vec variance() const;
        void sample_mat ( mat &Mi, chmat &Ri ) const;

        void factorize ( mat &M, ldmat &Vz, ldmat &Lam ) const;
        //! LS estimate of \f$\theta\f$
        vec est_theta() const;

        //! Covariance of the LS estimate
        ldmat est_theta_cov() const;

        //! expected values of the linear coefficient and the covariance matrix are written to \c M and \c R , respectively
        void mean_mat ( mat &M, mat&R ) const;
        //! In this instance, val= [theta, r]. For multivariate instances, it is stored columnwise val = [theta_1 theta_2 ... r_1 r_2 ]
        double evallog_nn ( const vec &val ) const;
        double lognc () const;
        void pow ( double p ) {
                V *= p;
                nu *= p;
        };

        //! marginal density (only student for now)
        shared_ptr<epdf> marginal(const RV &rvm) const {
                bdm_assert(dimx==1, "Not supported");
                //TODO - this is too trivial!!!
                ivec ind = rvm.findself_ids(rv);
                if (min(ind)==0) { //assume it si 
                        shared_ptr<estudent<ldmat> > tmp = new estudent<ldmat>;
                        mat M;
                        ldmat Vz;
                        ldmat Lam;
                        factorize(M,Vz,Lam);
                        
                        tmp->_mu() = M.get_col(0);
                        ldmat H;
                        Vz.inv(H);
                        H *=Lam._D()(0)/nu;
                        tmp->_H() = H;
                        tmp->_delta() = nu;
                        tmp->validate();
                        return tmp;
                }
                return NULL;
        }
        //! \name Access attributes
        //!@{

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

        /*! Create Gauss-inverse-Wishart density
        \f[ f(rv) = GiW(V,\nu) \f]
        from structure
        \code
        class = 'egiw';
        V     = [];               // square matrix
        dV    = [];               // vector of diagonal of V (when V not given)
        nu    = [];               // scalar \nu ((almost) degrees of freedom)
                                                          // when missing, it will be computed to obtain proper pdf
        dimx  = [];               // dimension of the wishart part
        rv = RV({'name'})         // description of RV
        rvc = RV({'name'})        // description of RV in condition
        log_level = 'tri';                // set the level of logged details
        \endcode

        \sa log_level_enums
        */
        void from_setting ( const Setting &set );
        void to_setting ( Setting& set ) const;
        void validate();
        void log_register ( bdm::logger& L, const string& prefix );

        void log_write() const;
        //!@}
};
UIREGISTER ( egiw );
SHAREDPTR ( egiw );

/*! \brief Dirichlet posterior density

Continuous Dirichlet density of \f$n\f$-dimensional variable \f$x\f$
\f[
f(x|\beta) = \frac{\Gamma[\gamma]}{\prod_{i=1}^{n}\Gamma(\beta_i)} \prod_{i=1}^{n}x_i^{\beta_i-1}
\f]
where \f$\gamma=\sum_i \beta_i\f$.
*/
class eDirich: public eEF {
protected:
        //!sufficient statistics
        vec beta;
public:
        //!\name Constructors
        //!@{

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

        //! using sampling procedure from wikipedia
        vec sample() const {
                vec y ( beta.length() );
                for ( int i = 0; i < beta.length(); i++ ) {
                        GamRNG.setup ( beta ( i ), 1 );
#pragma omp critical
                        y ( i ) = GamRNG();
                }
                return y / sum ( y );
        }

        vec mean() const {
                return beta / sum ( beta );
        };
        vec variance() const {
                double gamma = sum ( beta );
                return elem_mult ( beta, ( gamma - beta ) ) / ( gamma*gamma* ( gamma + 1 ) );
        }
        //! In this instance, val is ...
        double evallog_nn ( const vec &val ) const {
                double tmp;
                tmp = ( beta - 1 ) * log ( val );
                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 );
                return tmp;
        }

        //!access function
        vec& _beta()  {
                return beta;
        }
        /*! configuration structure
        \code
        class = 'eDirich';
        beta  = [];           //parametr beta
        \endcode
        */
        void from_setting ( const Setting &set );
        void validate();
        void to_setting ( Setting &set ) const;
};
UIREGISTER ( eDirich );

/*! Random Walk on Dirichlet
Using simple assignment
 \f[ \beta = rvc / k + \beta_c \f]
 hence, mean value = rvc, variance = (k+1)*mean*mean;

 The greater k is, the greater is the variance of the random walk;

 \f$ \beta_c \f$ is used as regularizing element to avoid corner cases, i.e. when one element of rvc is zero.
 By default is it set to 0.1;
*/

class mDirich: public pdf_internal<eDirich> {
protected:
        //! constant \f$ k \f$ of the random walk
        double k;
        //! cache of beta_i
        vec &_beta;
        //! stabilizing coefficient \f$ \beta_c \f$
        vec betac;
public:
        mDirich() : pdf_internal<eDirich>(), _beta ( iepdf._beta() ) {};
        void condition ( const vec &val ) {
                _beta =  val / k + betac;
        };
        /*! Create Dirichlet random walk
        \f[ f(rv|rvc) = Di(rvc*k) \f]
        from structure
        \code
        class = 'mDirich';
        k = 1;                      // multiplicative constant k
        --- optional ---
        rv = RV({'name'},size)      // description of RV
        beta0 = [];                 // initial value of beta
        betac = [];                 // initial value of beta
        \endcode
        */
        void from_setting ( const Setting &set );
        void    to_setting  (Setting  &set) const;
        void validate();
};
UIREGISTER ( mDirich );

//! \brief Estimator for Multinomial density
class multiBM : public BMEF {
protected:
        //! Conjugate prior and posterior
        eDirich est;
        //! Pointer inside est to sufficient statistics
        vec &beta;
public:
        //!Default constructor
        multiBM () : BMEF (), est (), beta ( est._beta() ) {
                if ( beta.length() > 0 ) {
                        last_lognc = est.lognc();
                } else {
                        last_lognc = 0.0;
                }
        }
        //!Copy constructor
        multiBM ( const multiBM &B ) : BMEF ( B ), est ( B.est ), beta ( est._beta() ) {}
        //! Sets sufficient statistics to match that of givefrom mB0
        void set_statistics ( const BM* mB0 ) {
                const multiBM* mB = dynamic_cast<const multiBM*> ( mB0 );
                beta = mB->beta;
        }
        void bayes ( const vec &yt, const vec &cond = empty_vec );

        double logpred ( const vec &yt ) const;

        void flatten ( const BMEF* B );

        //! return correctly typed posterior (covariant return)
        const eDirich& posterior() const {
                return est;
        };
        //! constructor function
        void set_parameters ( const vec &beta0 ) {
                est.set_parameters ( beta0 );
                est.validate();
                if ( evalll ) {
                        last_lognc = est.lognc();
                }
        }

        void to_setting ( Setting &set ) const {
                BMEF::to_setting ( set );
                UI::save( &est, set, "prior" );
        }
};
UIREGISTER( multiBM );

/*!
 \brief Gamma posterior density

 Multivariate Gamma density as product of independent univariate densities.
 \f[
 f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
 \f]
*/

class egamma : public eEF {
protected:
        //! Vector \f$\alpha\f$
        vec alpha;
        //! Vector \f$\beta\f$
        vec beta;
public :
        //! \name Constructors
        //!@{
        egamma () : eEF (), alpha ( 0 ), beta ( 0 ) {};
        egamma ( const vec &a, const vec &b ) {
                set_parameters ( a, b );
                validate();
        };
        void set_parameters ( const vec &a, const vec &b ) {
                alpha = a, beta = b;
        };
        //!@}

        vec sample() const;
        double evallog ( const vec &val ) const;
        double lognc () const;
        //! Returns pointer to internal alpha. Potentially dengerous: use with care!
        vec& _alpha() {
                return alpha;
        }
        //! Returns pointer to internal beta. Potentially dengerous: use with care!
        vec& _beta() {
                return beta;
        }
        vec mean() const {
                return elem_div ( alpha, beta );
        }
        vec variance() const {
                return elem_div ( alpha, elem_mult ( beta, beta ) );
        }

        /*! Create Gamma density
        \f[ f(rv|rvc) = \Gamma(\alpha, \beta) \f]
        from structure
        \code
        class = 'egamma';
        alpha = [...];         // vector of alpha
        beta = [...];          // vector of beta
        rv = RV({'name'})      // description of RV
        \endcode
        */
        void from_setting ( const Setting &set );
        void to_setting ( Setting &set ) const;
        void validate(); 
};
UIREGISTER ( egamma );
SHAREDPTR ( egamma );

/*!
 \brief Inverse-Gamma posterior density

 Multivariate inverse-Gamma density as product of independent univariate densities.
 \f[
 f(x|\alpha,\beta) = \prod f(x_i|\alpha_i,\beta_i)
 \f]

Vector \f$\beta\f$ has different meaning (in fact it is 1/beta as used in definition of iG)

 Inverse Gamma can be converted to Gamma using \f[
 x\sim iG(a,b) => 1/x\sim G(a,1/b)
\f]
This relation is used in sampling.
 */

class eigamma : public egamma {
protected:
public :
        //! \name Constructors
        //! All constructors are inherited
        //!@{
        //!@}

        vec sample() const {
                return 1.0 / egamma::sample();
        };
        //! Returns poiter to alpha and beta. Potentially dangerous: use with care!
        vec mean() const {
                return elem_div ( beta, alpha - 1 );
        }
        vec variance() const {
                vec mea = mean();
                return elem_div ( elem_mult ( mea, mea ), alpha - 2 );
        }
};
/*
//! Weighted mixture of epdfs with external owned components.
class emix : public epdf {
protected:
        int n;
        vec &w;
        Array<epdf*> Coms;
public:
//! Default constructor
        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;};
};
*/

//!  Uniform distributed density on a rectangular support

class euni: public epdf {
protected:
//! lower bound on support
        vec low;
//! upper bound on support
        vec high;
//! internal
        vec distance;
//! normalizing coefficients
        double nk;
//! cache of log( \c nk )
        double lnk;
public:
        //! \name Constructors
        //!@{
        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;
                low = low0;
                high = high0;
                nk = prod ( 1.0 / distance );
                lnk = log ( nk );
        }
        //!@}

        double evallog ( const vec &val ) const  {
                if ( any ( val < low ) && any ( val > high ) ) {
                        return -inf;
                } else return lnk;
        }
        vec sample() const {
                vec smp ( dim );
#pragma omp critical
                UniRNG.sample_vector ( dim , smp );
                return low + elem_mult ( distance, smp );
        }
        //! set values of \c low and \c high
        vec mean() const {
                return ( high - low ) / 2.0;
        }
        vec variance() const {
                return ( pow ( high, 2 ) + pow ( low, 2 ) + elem_mult ( high, low ) ) / 3.0;
        }
        /*! Create Uniform density
        \f[ f(rv) = U(low,high) \f]
        from structure
         \code
         class = 'euni'
         high = [...];          // vector of upper bounds
         low = [...];           // vector of lower bounds
         rv = RV({'name'});     // description of RV
         \endcode
         */
        void from_setting ( const Setting &set );
        void    to_setting  (Setting  &set) const;
        void validate();
};
UIREGISTER ( euni );

//! Uniform density with conditional mean value
class mguni : public pdf_internal<euni> {
        //! function of the mean value
        shared_ptr<fnc> mean;
        //! distance from mean to both sides
        vec delta;
public:
        void condition ( const vec &cond ) {
                vec mea = mean->eval ( cond );
                iepdf.set_parameters ( mea - delta, mea + delta );
        }
        //! load from
        void from_setting ( const Setting &set ) {
                pdf::from_setting ( set ); //reads rv and rvc
                UI::get ( delta, set, "delta", UI::compulsory );
                mean = UI::build<fnc> ( set, "mean", UI::compulsory );
                iepdf.set_parameters ( -delta, delta );
        }
        void    to_setting  (Setting  &set) const {
                pdf::to_setting ( set ); 
                UI::save( iepdf.mean(), set, "delta");
                UI::save(mean, set, "mean");
        }
        void validate(){
                pdf_internal<euni>::validate();
                dimc = mean->dimensionc();
                
        }

};
UIREGISTER ( mguni );
/*!
 \brief Normal distributed linear function with linear function of mean value;

 Mean value \f$ \mu=A*\mbox{rvc}+\mu_0 \f$.
*/
template < class sq_T, template <typename> class TEpdf = enorm >
class mlnorm : public pdf_internal< TEpdf<sq_T> > {
protected:
        //! Internal epdf that arise by conditioning on \c rvc
        mat A;
        //! Constant additive term
        vec mu_const;
//                      vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
public:
        //! \name Constructors
        //!@{
        mlnorm() : pdf_internal< TEpdf<sq_T> >() {};
        mlnorm ( const mat &A, const vec &mu0, const sq_T &R ) : pdf_internal< TEpdf<sq_T> >() {
                set_parameters ( A, mu0, R );
                validate();
        }

        //! Set \c A and \c R
        void set_parameters ( const  mat &A0, const vec &mu0, const sq_T &R0 ) {
                this->iepdf.set_parameters ( zeros ( A0.rows() ), R0 );
                A = A0;
                mu_const = mu0;
        }

        //!@}
        //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
        void condition ( const vec &cond ) {
                this->iepdf._mu() = A * cond + mu_const;
//R is already assigned;
        }

        //!access function
        const vec& _mu_const() const {
                return mu_const;
        }
        //!access function
        const mat& _A() const {
                return A;
        }
        //!access function
        mat _R() const {
                return this->iepdf._R().to_mat();
        }
        //!access function
        sq_T __R() const {
                return this->iepdf._R();
        }

        //! Debug stream
        template<typename sq_M>
        friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M, enorm> &ml );

        /*! Create Normal density with linear function of mean value
         \f[ f(rv|rvc) = N(A*rvc+const, R) \f]
        from structure
        \code
        class = 'mlnorm<ldmat>', (OR) 'mlnorm<chmat>', (OR) 'mlnorm<fsqmat>';
        A     = [];                  // matrix or vector of appropriate dimension
        R     = [];                  // square matrix of appropriate dimension
        --- optional ---
        const = zeros(A.rows);       // vector of constant term
        \endcode
        */
        void from_setting ( const Setting &set ) {
                pdf::from_setting ( set );

                UI::get ( A, set, "A", UI::compulsory );
                UI::get ( mu_const, set, "const", UI::optional);
                mat R0;
                UI::get ( R0, set, "R", UI::compulsory );
                set_parameters ( A, mu_const, R0 );
        }

void to_setting (Setting &set) const {
                pdf::to_setting(set);
                UI::save ( A, set, "A");
                UI::save ( mu_const, set, "const");
                UI::save ( _R(), set, "R");
        }

void validate() {
                pdf_internal<TEpdf<sq_T> >::validate();
                if (mu_const.length()==0) { // default in from_setting
                        mu_const=zeros(A.rows());
                }
                bdm_assert ( A.rows() == mu_const.length(), "mlnorm: A vs. mu mismatch" );
                bdm_assert ( A.rows() == _R().rows(), "mlnorm: A vs. R mismatch" );
                this->dimc = A.cols();

        }
};
UIREGISTER2 ( mlnorm, ldmat );
SHAREDPTR2 ( mlnorm, ldmat );
UIREGISTER2 ( mlnorm, fsqmat );
SHAREDPTR2 ( mlnorm, fsqmat );
UIREGISTER2 ( mlnorm, chmat );
SHAREDPTR2 ( mlnorm, chmat );

//! \class mlgauss 
//!\brief Normal distribution with linear function of mean value. Same as mlnorm<fsqmat>.
typedef mlnorm<fsqmat> mlgauss;
UIREGISTER(mlgauss);

//! pdf with general function for mean value
template<class sq_T>
class mgnorm : public pdf_internal< enorm< sq_T > > {
private:
//                      vec &mu; WHY NOT?
        shared_ptr<fnc> g;

public:
        //!default constructor
        mgnorm() : pdf_internal<enorm<sq_T> >() { }
        //!set mean function
        inline void set_parameters ( const shared_ptr<fnc> &g0, const sq_T &R0 );
        inline void condition ( const vec &cond );


        /*! Create Normal density with given function of mean value
        \f[ f(rv|rvc) = N( g(rvc), R) \f]
        from structure
        \code
        class = 'mgnorm';
        g.class =  'fnc';      // function for mean value evolution
        g._fields_of_fnc = ...;

        R = [1, 0;             // covariance matrix
                0, 1];
                --OR --
        dR = [1, 1];           // diagonal of cavariance matrix

        rv = RV({'name'})      // description of RV
        rvc = RV({'name'})     // description of RV in condition
        \endcode
        */


void from_setting ( const Setting &set ) {
                pdf::from_setting ( set );
                shared_ptr<fnc> g = UI::build<fnc> ( set, "g", UI::compulsory );

                mat R;
                vec dR;
                if ( UI::get ( dR, set, "dR" ) )
                        R = diag ( dR );
                else
                        UI::get ( R, set, "R", UI::compulsory );

                set_parameters ( g, R );
                //validate();
        }
        

void to_setting  (Setting  &set) const {
                UI::save( g,set, "g");
                UI::save(iepdf._R().to_mat(),set, "R");
        
        }



void validate() {
                this->iepdf.validate();
                bdm_assert ( g->dimension() == this->iepdf.dimension(), "incompatible function" );
                this->dim = g->dimension();
                this->dimc = g->dimensionc();
                this->iepdf.validate();
        }

};

UIREGISTER2 ( mgnorm, chmat );
UIREGISTER2 ( mgnorm, ldmat );
SHAREDPTR2 ( mgnorm, chmat );


/*! (Approximate) Student t density with linear function of mean value

The internal epdf of this class is of the type of a Gaussian (enorm).
However, each conditioning is trying to assure the best possible approximation by taking into account the zeta function. See [] for reference.

Perhaps a moment-matching technique?
*/
class mlstudent : public mlnorm<ldmat, enorm> {
protected:
        //! Variable \f$ \Lambda \f$ from theory
        ldmat Lambda;
        //! Reference to variable \f$ R \f$
        ldmat &_R;
        //! Variable \f$ R_e \f$
        ldmat Re;
public:
        mlstudent () : mlnorm<ldmat, enorm> (),
                        Lambda (),      _R ( iepdf._R() ) {}
        //! constructor function
        void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 ) {
                iepdf.set_parameters ( mu0, R0 );// was Lambda, why?
                A = A0;
                mu_const = mu0;
                Re = R0;
                Lambda = Lambda0;
        }

        void condition ( const vec &cond ); 

        void validate() {
                mlnorm<ldmat, enorm>::validate();
                bdm_assert ( A.rows() == mu_const.length(), "mlstudent: A vs. mu mismatch" );
                bdm_assert ( _R.rows() == A.rows(), "mlstudent: A vs. R mismatch" );

        }
};

/*!
\brief  Gamma random walk

Mean value, \f$\mu\f$, of this density is given by \c rvc .
Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.

The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
*/
class mgamma : public pdf_internal<egamma> {
protected:

        //! Constant \f$k\f$
        double k;

        //! cache of iepdf.beta
        vec &_beta;

public:
        //! Constructor
        mgamma() : pdf_internal<egamma>(), k ( 0 ),
                        _beta ( iepdf._beta() ) {
        }

        //! Set value of \c k
        void set_parameters ( double k, const vec &beta0 );

        void condition ( const vec &val ) {
                _beta = k / val;
        };
        /*! Create Gamma density with conditional mean value
        \f[ f(rv|rvc) = \Gamma(k, k/rvc) \f]
        from structure
        \code
          class = 'mgamma';
          beta = [...];          // vector of initial alpha
          k = 1.1;               // multiplicative constant k
          rv = RV({'name'})      // description of RV
          rvc = RV({'name'})     // description of RV in condition
         \endcode
        */
        void from_setting ( const Setting &set );
        void    to_setting  (Setting  &set) const;
        void validate();
};
UIREGISTER ( mgamma );
SHAREDPTR ( mgamma );

/*!
\brief  Inverse-Gamma random walk

Mean value, \f$ \mu \f$, of this density is given by \c rvc .
Standard deviation of the random walk is proportional to one \f$ k \f$-th the mean.
This is achieved by setting \f$ \alpha=\mu/k^2+2 \f$ and \f$ \beta=\mu(\alpha-1)\f$.

The standard deviation of the walk is then: \f$ \mu/\sqrt(k)\f$.
 */
class migamma : public pdf_internal<eigamma> {
protected:
        //! Constant \f$k\f$
        double k;

        //! cache of iepdf.alpha
        vec &_alpha;

        //! cache of iepdf.beta
        vec &_beta;

public:
        //! \name Constructors
        //!@{
        migamma() : pdf_internal<eigamma>(),
                        k ( 0 ),
                        _alpha ( iepdf._alpha() ),
                        _beta ( iepdf._beta() ) {
        }

        migamma ( const migamma &m ) : pdf_internal<eigamma>(),
                        k ( 0 ),
                        _alpha ( iepdf._alpha() ),
                        _beta ( iepdf._beta() ) {
        }
        //!@}

        //! Set value of \c k
        void set_parameters ( int len, double k0 ) {
                k = k0;
                iepdf.set_parameters ( ( 1.0 / ( k*k ) + 2.0 ) *ones ( len ) /*alpha*/, ones ( len ) /*beta*/ );                
        };

        void    validate (){
                pdf_internal<eigamma>::validate();
                dimc = dimension();
};

        void condition ( const vec &val ) {
                _beta = elem_mult ( val, ( _alpha - 1.0 ) );
        };
};


/*!
\brief  Gamma random walk around a fixed point

Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]

Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.

The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
*/
class mgamma_fix : public mgamma {
protected:
        //! parameter l
        double l;
        //! reference vector
        vec refl;
public:
        //! Constructor
        mgamma_fix () : mgamma (), refl () {};
        //! Set value of \c k
        void set_parameters ( double k0 , vec ref0, double l0 ) {
                mgamma::set_parameters ( k0, ref0 );
                refl = pow ( ref0, 1.0 - l0 );
                l = l0; 
        };

        void    validate (){
                mgamma::validate();
                dimc = dimension();
        };

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


/*!
\brief  Inverse-Gamma random walk around a fixed point

Mean value, \f$\mu\f$, of this density is given by a geometric combination of \c rvc and given fixed point, \f$p\f$. \f$l\f$ is the coefficient of the geometric combimation
\f[ \mu = \mu_{t-1} ^{l} p^{1-l}\f]

==== Check == vv =
Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
This is achieved by setting \f$\alpha=k\f$ and \f$\beta=k/\mu\f$.

The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
 */
class migamma_ref : public migamma {
protected:
        //! parameter l
        double l;
        //! reference vector
        vec refl;
public:
        //! Constructor
        migamma_ref () : migamma (), refl () {};
        
        //! Set value of \c k
        void set_parameters ( double k0 , vec ref0, double l0 ) {
                migamma::set_parameters ( ref0.length(), k0 );
                refl = pow ( ref0, 1.0 - l0 );
                l = l0;
        };

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


        /*! Create inverse-Gamma density with conditional mean value
        \f[ f(rv|rvc) = i\Gamma(k, k/(rvc^l \circ ref^{(1-l)}) \f]
        from structure
        \code
        class = 'migamma_ref';
        ref = [1e-5; 1e-5; 1e-2 1e-3];            // reference vector
        l = 0.999;                                // constant l
        k = 0.1;                                  // constant k
        rv = RV({'name'})                         // description of RV
        rvc = RV({'name'})                        // description of RV in condition
        \endcode
         */
        void from_setting ( const Setting &set );
        void to_setting  (Setting  &set) const; 
        // TODO dodelat void to_setting( Setting &set ) const;
};


UIREGISTER ( migamma_ref );
SHAREDPTR ( migamma_ref );

/*! Log-Normal probability density
 only allow diagonal covariances!

Density of the form \f$ \log(x)\sim \mathcal{N}(\mu,\sigma^2) \f$ , i.e.
\f[
x \sim \frac{1}{x\sigma\sqrt{2\pi}}\exp{-\frac{1}{2\sigma^2}(\log(x)-\mu)}
\f]

Function from_setting loads mu and R in the same way as it does for enorm<>!
*/
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 );
        };

};

/*!
\brief  Log-Normal random walk

Mean value, \f$\mu\f$, is...

 */
class mlognorm : public pdf_internal<elognorm> {
protected:
        //! parameter 1/2*sigma^2
        double sig2;

        //! access
        vec &mu;
public:
        //! Constructor
        mlognorm() : pdf_internal<elognorm>(),
                        sig2 ( 0 ),
                        mu ( iepdf._mu() ) {
        }

        //! Set value of \c k
        void set_parameters ( int size, double k ) {
                sig2 = 0.5 * log ( k * k + 1 );
                iepdf.set_parameters ( zeros ( size ), 2*sig2*eye ( size ) );
        };
        
        void validate(){
                pdf_internal<elognorm>::validate();
                dimc = iepdf.dimension();
        }

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

        /*! Create logNormal random Walk
        \f[ f(rv|rvc) = log\mathcal{N}( \log(rvc)-0.5\log(k^2+1), k I) \f]
        from structure
        \code
        class = 'mlognorm';
        k   = 0.1;               // "variance" k
        mu0 = 0.1;               // Initial value of mean
        rv  = RV({'name'})       // description of RV
        rvc = RV({'name'})       // description of RV in condition
        \endcode
        */
        void from_setting ( const Setting &set );
        void to_setting  (Setting  &set) const; 

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

};

UIREGISTER ( mlognorm );
SHAREDPTR ( mlognorm );

/*! inverse Wishart density defined on Choleski decomposition

*/
class eWishartCh : public epdf {
protected:
        //! Upper-Triagle of Choleski decomposition of \f$ \Psi \f$
        chmat Y;
        //! dimension of matrix  \f$ \Psi \f$
        int p;
        //! degrees of freedom  \f$ \nu \f$
        double delta;
public:
        //! Set internal structures
        void set_parameters ( const mat &Y0, const double delta0 ) {
                Y = chmat ( Y0 );
                delta = delta0;
                p = Y.rows();   
        }
        //! Set internal structures
        void set_parameters ( const chmat &Y0, const double delta0 ) {
                Y = Y0;
                delta = delta0;
                p = Y.rows();
                }
        
        virtual void validate (){
                epdf::validate();
                dim = p * p;
        }

        //! Sample matrix argument
        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 );
        }

        virtual vec mean() const NOT_IMPLEMENTED(0);

        //! return expected variance (not covariance!)
        virtual vec variance() const NOT_IMPLEMENTED(0);

        virtual double evallog ( const vec &val ) const NOT_IMPLEMENTED(0);

        //! fast access function y0 will be copied into Y.Ch.
        void setY ( const mat &Ch0 ) {
                copy_vector ( dim, Ch0._data(), Y._Ch()._data() );
        }

        //! fast access function y0 will be copied into Y.Ch.
        void _setY ( const vec &ch0 ) {
                copy_vector ( dim, ch0._data(), Y._Ch()._data() );
        }

        //! access function
        const chmat& getY() const {
                return Y;
        }
};

//! Inverse Wishart on Choleski decomposition
/*! Being computed by conversion from `standard' Wishart
*/
class eiWishartCh: public epdf {
protected:
        //! Internal instance of Wishart density
        eWishartCh W;
        //! size of Ch
        int p;
        //! parameter delta
        double delta;
public:
        //! constructor function
        void set_parameters ( const mat &Y0, const double delta0 ) {
                delta = delta0;
                W.set_parameters ( inv ( Y0 ), delta0 );
                p = Y0.rows();
        }
        
        virtual void    validate (){
                epdf::validate();
                W.validate();
                dim = W.dimension();
        }
        
        
        vec sample() const {
                mat iCh;
                iCh = inv ( W.sample_mat() );
                return vec ( iCh._data(), dim );
        }
        //! access function
        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;
        };

        virtual vec mean() const NOT_IMPLEMENTED(0);

        //! return expected variance (not covariance!)
        virtual vec variance() const NOT_IMPLEMENTED(0);
};

//! Random Walk on inverse Wishart
class rwiWishartCh : public pdf_internal<eiWishartCh> {
protected:
        //!square root of \f$ \nu-p-1 \f$ - needed for computation of \f$ \Psi \f$ from conditions
        double sqd;
        //!reference point for diagonal
        vec refl;
        //! power of the reference
        double l;
        //! dimension
        int p;

public:
        rwiWishartCh() : sqd ( 0 ), l ( 0 ), p ( 0 ) {}
        //! constructor function
        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 );
                iepdf.set_parameters ( eye ( p ), delta );
        };
        
        void validate(){
                pdf_internal<eiWishartCh>::validate();
                dimc = iepdf.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++;
                }

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

//! Switch between various resampling methods.
enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };

//! Shortcut for multinomial.sample(int n). Various simplifications are allowed see RESAMPLING_METHOD 
void resample(const vec &w, ivec &ind, RESAMPLING_METHOD=SYSTEMATIC);

/*!
\brief Weighted empirical density

Used e.g. in particle filters.
*/
class eEmp: public epdf {
protected :
        //! Number of particles
        int n;
        //! Sample weights \f$w\f$
        vec w;
        //! Samples \f$x^{(i)}, i=1..n\f$
        Array<vec> samples;
public:
        //! \name Constructors
        //!@{
        eEmp () : epdf (), w (), samples () {};
        //! copy constructor
        eEmp ( const eEmp &e ) : epdf ( e ), w ( e.w ), samples ( e.samples ) {};
        //!@}

        //! Set samples and weights
        void set_statistics ( const vec &w0, const epdf &pdf0 );
        //! Set samples and weights
        void set_statistics ( const epdf &pdf0 , int n ) {
                set_statistics ( ones ( n ) / n, pdf0 );
        };
        //! Set sample
        void set_samples ( const epdf* pdf0 );
        //! Set sample
        void set_parameters ( int n0, bool copy = true ) {
                n = n0;
                w.set_size ( n0, copy );
                samples.set_size ( n0, copy );
        };
        //! Set samples
        void set_parameters ( const Array<vec> &Av ) {
                n = Av.size();
                w = 1 / n * ones ( n );
                samples = Av;
        };
        virtual void    validate ();
        //! Potentially dangerous, use with care.
        vec& _w()  {
                return w;
        };
        //! Potentially dangerous, use with care.
        const vec& _w() const {
                return w;
        };
        //! access function
        Array<vec>& _samples() {
                return samples;
        };
        //! access function
        const vec& _sample ( int i ) const {
                return samples ( i );
        };
        //! access function
        const Array<vec>& _samples() const {
                return samples;
        };
        //! Function performs resampling, i.e. removal of low-weight samples and duplication of high-weight samples such that the new samples represent the same density.
        void resample ( RESAMPLING_METHOD method = SYSTEMATIC );

        //! inherited operation : NOT implemented
        vec sample() const NOT_IMPLEMENTED(0);

        //! inherited operation : NOT implemented
        double evallog ( const vec &val ) const NOT_IMPLEMENTED(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 );
        }
        //! For this class, qbounds are minimum and maximum value of the population!
        void qbounds ( vec &lb, vec &ub, double perc = 0.95 ) const;

        void to_setting ( Setting &set ) const;
        void from_setting ( const Setting &set );

};
UIREGISTER(eEmp);


////////////////////////

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;
        validate();
};

template<class sq_T>
void enorm<sq_T>::from_setting ( const Setting &set ) {
        epdf::from_setting ( set ); //reads rv

        UI::get ( mu, set, "mu", UI::compulsory );
        mat Rtmp;// necessary for conversion
        UI::get ( Rtmp, set, "R", UI::compulsory );
        R = Rtmp; // conversion
}

template<class sq_T>
void enorm<sq_T>::validate() {
                eEF::validate();
                bdm_assert ( mu.length() == R.rows(), "mu and R parameters do not match" );
                dim = mu.length();
        }

template<class sq_T>
void enorm<sq_T>::to_setting ( Setting &set ) const {
        epdf::to_setting ( set ); //reads rv    
        UI::save ( mu, set, "mu");
        UI::save ( R.to_mat(), set, "R");
}



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>
// 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>
// 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>
shared_ptr<epdf> enorm<sq_T>::marginal ( const RV &rvn ) const {
        enorm<sq_T> *tmp = new enorm<sq_T> ();
        shared_ptr<epdf> narrow ( tmp );
        marginal ( rvn, *tmp );
        return narrow;
}

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

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

        target.set_rv ( rvn );
        target.set_parameters ( mu ( irvn ), Rn );
}

template<class sq_T>
shared_ptr<pdf> enorm<sq_T>::condition ( const RV &rvn ) const {
        mlnorm<sq_T> *tmp = new mlnorm<sq_T> ();
        shared_ptr<pdf> narrow ( tmp );
        condition ( rvn, *tmp );
        return narrow;
}

template<class sq_T>
void enorm<sq_T>::condition ( const RV &rvn, pdf &target ) const {
        typedef mlnorm<sq_T> TMlnorm;

        bdm_assert ( isnamed(), "rvs are not assigned" );
        TMlnorm &uptarget = dynamic_cast<TMlnorm &> ( target );

        RV rvc = rv.subt ( rvn );
        bdm_assert ( ( 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() );

        uptarget.set_rv ( rvn );
        uptarget.set_rvc ( rvc );
        uptarget.set_parameters ( A, mu1 - A*mu2, R_n );
        uptarget.validate();
}

/*! Dirac delta function distribution */
class dirac: public epdf{
        public: 
                vec point;
        public:
                double evallog (const vec &dt) const {return -inf;}
                vec mean () const {return point;}
                vec variance () const {return zeros(point.length());}
                void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const { lb = point; ub = point;}
                //! access
                const vec& _point() {return point;}
                void set_point(const vec& p){point =p; dim=p.length();}
                vec sample() const {return point;}
        };

////
///////
template<class sq_T>
void mgnorm<sq_T >::set_parameters ( const shared_ptr<fnc> &g0, const sq_T &R0 ) {
        g = g0;
        this->iepdf.set_parameters ( zeros ( g->dimension() ), R0 );
}

template<class sq_T>
void mgnorm<sq_T >::condition ( const vec &cond ) {
        this->iepdf._mu() = g->eval ( cond );
};

//! \todo unify this stuff with to_string()
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._R() << endl;
        return os;
};

}
#endif //EF_H