Changeset 488 for library/bdm

Timestamp:

08/08/09 13:43:13 (15 years ago)

Author:

smidl

Message:

changes in mpdf -> compile OK, broken tests!

Location:

library/bdm

Files:

• base/bdmbase.cpp (modified) (1 diff)
• base/bdmbase.h (modified) (11 diffs)
• estim/particles.h (modified) (2 diffs)
• stat/emix.cpp (modified) (1 diff)
• stat/emix.h (modified) (1 diff)
• stat/exp_family.cpp (modified) (1 diff)
• stat/exp_family.h (modified) (4 diffs)
• stat/merger.h (modified) (1 diff)

library/bdm/base/bdmbase.cpp

@@ -140 +140 @@
 }
 
-template<class EPDF>
-vec mpdf_internal<EPDF>::samplecond ( const vec &cond ) {
-        condition ( cond );
-        vec temp = iepdf.sample();
-        return temp;
-}
-
-template<class EPDF>
-mat mpdf_internal<EPDF>::samplecond_m ( const vec &cond, int N ) {
-        condition ( cond );
-        mat temp ( dimension(), N );
-        vec smp ( dimension() );
-        for ( int i = 0; i < N; i++ ) {
-                smp = iepdf.sample();
-                temp.set_col ( i, smp );
-        }
-
-        return temp;
-}
-
-template<class EPDF>
-double mpdf_internal<EPDF>::evallogcond ( const vec &dt, const vec &cond ) {
-        double tmp;
-        condition ( cond );
-        tmp = iepdf.evallog ( dt );
-        // it_assert_debug(std::isfinite(tmp), "Infinite value");
-        return tmp;
-}
-
-template<class EPDF>
-vec mpdf_internal<EPDF>::evallogcond_m ( const mat &Dt, const vec &cond ) {
-        condition ( cond );
-        return iepdf.evallog_m ( Dt );
-}
-
-template<class EPDF>
-                vec mpdf_internal<EPDF>::evallogcond_m ( const Array<vec> &Dt, const vec &cond ) {
-        condition ( cond );
-        return iepdf.evallog_m ( Dt );
-}
 
 void mpdf::from_setting ( const Setting &set ) {

library/bdm/base/bdmbase.h

@@ -25 +25 @@
 using namespace std;
 
-namespace bdm {
+namespace bdm
+{
 
 typedef std::map<string, int> RVmap;
@@ -32 +33 @@
 
 //! Structure of RV, i.e. RVs expanded into a flat list of IDs, used for debugging.
-class str {
-public:
-        //! vector id ids (non-unique!)
-        ivec ids;
-        //! vector of times
-        ivec times;
-        //!Default constructor
-        str ( ivec ids0, ivec times0 ) : ids ( ids0 ), times ( times0 ) {
-                it_assert_debug ( times0.length() == ids0.length(), "Incompatible input" );
-        };
+class str
+{
+        public:
+                //! vector id ids (non-unique!)
+                ivec ids;
+                //! vector of times
+                ivec times;
+                //!Default constructor
+                str (ivec ids0, ivec times0) : ids (ids0), times (times0) {
+                        it_assert_debug (times0.length() == ids0.length(), "Incompatible input");
+                };
 };
 
@@ -82 +84 @@
 */
 
-class RV : public root {
-protected:
-        //! size of the data vector
-        int dsize;
-        //! number of individual rvs
-        int len;
-        //! Vector of unique IDs
-        ivec ids;
-        //! Vector of shifts from current time
-        ivec times;
-
-private:
-        //! auxiliary function used in constructor
-        void init ( const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times );
-        int init ( const string &name, int size );
-public:
-        //! \name Constructors
-        //!@{
-
-        //! Full constructor
-        RV ( const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times ) {
-                init ( in_names, in_sizes, in_times );
-        }
-
-        //! Constructor with times=0
-        RV ( const Array<std::string> &in_names, const ivec &in_sizes ) {
-                init ( in_names, in_sizes, zeros_i ( in_names.length() ) );
-        }
-
-        //! Constructor with sizes=1, times=0
-        RV ( const Array<std::string> &in_names ) {
-                init ( in_names, ones_i ( in_names.length() ), zeros_i ( in_names.length() ) );
-        }
-
-        //! Constructor of empty RV
-        RV() : dsize ( 0 ), len ( 0 ), ids ( 0 ), times ( 0 ) {}
-        //! Constructor of a single RV with given id
-        RV ( string name, int sz, int tm = 0 );
-        //!@}
-
-        //! \name Access functions
-        //!@{
-
-        //! State output, e.g. for debugging.
-        friend std::ostream &operator<< ( std::ostream &os, const RV &rv );
-
-        int _dsize() const {
-                return dsize;
-        }
-
-        //! Recount size of the corresponding data vector
-        int countsize() const;
-        ivec cumsizes() const;
-        int length() const {
-                return len;
-        }
-        int id ( int at ) const {
-                return ids ( at );
-        }
-        int size ( int at ) const {
-                return RV_SIZES ( ids ( at ) );
-        }
-        int time ( int at ) const {
-                return times ( at );
-        }
-        std::string name ( int at ) const {
-                return RV_NAMES ( ids ( at ) );
-        }
-        void set_time ( int at, int time0 ) {
-                times ( at ) = time0;
-        }
-        //!@}
-
-        //TODO why not inline and later??
-
-        //! \name Algebra on Random Variables
-        //!@{
-
-        //! Find indices of self in another rv, \return ivec of the same size as self.
-        ivec findself ( const RV &rv2 ) const;
-        //! Compare if \c rv2 is identical to this \c RV
-        bool equal ( const RV &rv2 ) const;
-        //! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict
-        bool add ( const RV &rv2 );
-        //! Subtract  another variable from the current one
-        RV subt ( const RV &rv2 ) const;
-        //! Select only variables at indices ind
-        RV subselect ( const ivec &ind ) const;
-
-        //! Select only variables at indices ind
-        RV operator() ( const ivec &ind ) const {
-                return subselect ( ind );
-        }
-
-        //! Select from data vector starting at di1 to di2
-        RV operator() ( int di1, int di2 ) const;
-
-        //! Shift \c time by delta.
-        void t ( int delta );
-        //!@}
-
-        //!\name Relation to vectors
-        //!@{
-
-        //! generate \c str from rv, by expanding sizes TODO to_string..
-        str tostr() const;
-        //! when this rv is a part of bigger rv, this function returns indices of self in the data vector of the bigger crv.
-        //! Then, data can be copied via: data_of_this = cdata(ind);
-        ivec dataind ( const RV &crv ) const;
-        //! generate mutual indices when copying data between self and crv.
-        //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i)
-        void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const;
-        //! Minimum time-offset
-        int mint() const {
-                return min ( times );
-        }
-        //!@}
-
-        // TODO aktualizovat dle soucasneho UI
-        /*! \brief UI for class RV (description of data vectors)
-
-        \code
-        rv = {
-          type = "rv"; //identifier of the description
-          // UNIQUE IDENTIFIER same names = same variable
-          names = ["a", "b", "c", ...];   // which will be used e.g. in loggers
-
-          //optional arguments
-          sizes = [1, 2, 3, ...];         // (optional) default = ones()
-          times = [-1, -2, 0, ...];       // time shifts with respect to current time (optional) default = zeros()
-        }
-        \endcode
-        */
-        void from_setting ( const Setting &set );
-
-        // TODO dodelat void to_setting( Setting &set ) const;
-
-        //! Invalidate all named RVs. Use before initializing any RV instances, with care...
-        static void clear_all();
-};
-UIREGISTER ( RV );
+class RV : public root
+{
+        protected:
+                //! size of the data vector
+                int dsize;
+                //! number of individual rvs
+                int len;
+                //! Vector of unique IDs
+                ivec ids;
+                //! Vector of shifts from current time
+                ivec times;
+
+        private:
+                //! auxiliary function used in constructor
+                void init (const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times);
+                int init (const string &name, int size);
+        public:
+                //! \name Constructors
+                //!@{
+
+                //! Full constructor
+                RV (const Array<std::string> &in_names, const ivec &in_sizes, const ivec &in_times) {
+                        init (in_names, in_sizes, in_times);
+                }
+
+                //! Constructor with times=0
+                RV (const Array<std::string> &in_names, const ivec &in_sizes) {
+                        init (in_names, in_sizes, zeros_i (in_names.length()));
+                }
+
+                //! Constructor with sizes=1, times=0
+                RV (const Array<std::string> &in_names) {
+                        init (in_names, ones_i (in_names.length()), zeros_i (in_names.length()));
+                }
+
+                //! Constructor of empty RV
+                RV() : dsize (0), len (0), ids (0), times (0) {}
+                //! Constructor of a single RV with given id
+                RV (string name, int sz, int tm = 0);
+                //!@}
+
+                //! \name Access functions
+                //!@{
+
+                //! State output, e.g. for debugging.
+                friend std::ostream &operator<< (std::ostream &os, const RV &rv);
+
+                int _dsize() const {
+                        return dsize;
+                }
+
+                //! Recount size of the corresponding data vector
+                int countsize() const;
+                ivec cumsizes() const;
+                int length() const {
+                        return len;
+                }
+                int id (int at) const {
+                        return ids (at);
+                }
+                int size (int at) const {
+                        return RV_SIZES (ids (at));
+                }
+                int time (int at) const {
+                        return times (at);
+                }
+                std::string name (int at) const {
+                        return RV_NAMES (ids (at));
+                }
+                void set_time (int at, int time0) {
+                        times (at) = time0;
+                }
+                //!@}
+
+                //TODO why not inline and later??
+
+                //! \name Algebra on Random Variables
+                //!@{
+
+                //! Find indices of self in another rv, \return ivec of the same size as self.
+                ivec findself (const RV &rv2) const;
+                //! Compare if \c rv2 is identical to this \c RV
+                bool equal (const RV &rv2) const;
+                //! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict
+                bool add (const RV &rv2);
+                //! Subtract  another variable from the current one
+                RV subt (const RV &rv2) const;
+                //! Select only variables at indices ind
+                RV subselect (const ivec &ind) const;
+
+                //! Select only variables at indices ind
+                RV operator() (const ivec &ind) const {
+                        return subselect (ind);
+                }
+
+                //! Select from data vector starting at di1 to di2
+                RV operator() (int di1, int di2) const;
+
+                //! Shift \c time by delta.
+                void t (int delta);
+                //!@}
+
+                //!\name Relation to vectors
+                //!@{
+
+                //! generate \c str from rv, by expanding sizes TODO to_string..
+                str tostr() const;
+                //! when this rv is a part of bigger rv, this function returns indices of self in the data vector of the bigger crv.
+                //! Then, data can be copied via: data_of_this = cdata(ind);
+                ivec dataind (const RV &crv) const;
+                //! generate mutual indices when copying data between self and crv.
+                //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i)
+                void dataind (const RV &rv2, ivec &selfi, ivec &rv2i) const;
+                //! Minimum time-offset
+                int mint() const {
+                        return min (times);
+                }
+                //!@}
+
+                // TODO aktualizovat dle soucasneho UI
+                /*! \brief UI for class RV (description of data vectors)
+
+                \code
+                rv = {
+                  type = "rv"; //identifier of the description
+                  // UNIQUE IDENTIFIER same names = same variable
+                  names = ["a", "b", "c", ...];   // which will be used e.g. in loggers
+
+                  //optional arguments
+                  sizes = [1, 2, 3, ...];         // (optional) default = ones()
+                  times = [-1, -2, 0, ...];       // time shifts with respect to current time (optional) default = zeros()
+                }
+                \endcode
+                */
+                void from_setting (const Setting &set);
+
+                // TODO dodelat void to_setting( Setting &set ) const;
+
+                //! Invalidate all named RVs. Use before initializing any RV instances, with care...
+                static void clear_all();
+};
+UIREGISTER (RV);
 
 //! Concat two random variables
-RV concat ( const RV &rv1, const RV &rv2 );
+RV concat (const RV &rv1, const RV &rv2);
 
 //!Default empty RV that can be used as default argument
@@ -232 +235 @@
 //! Class representing function \f$f(x)\f$ of variable \f$x\f$ represented by \c rv
 
-class fnc : public root {
-protected:
-        //! Length of the output vector
-        int dimy;
-public:
-        //!default constructor
-        fnc() {};
-        //! function evaluates numerical value of \f$f(x)\f$ at \f$x=\f$ \c cond
-        virtual vec eval ( const vec &cond ) {
-                return vec ( 0 );
-        };
-
-        //! function substitutes given value into an appropriate position
-        virtual void condition ( const vec &val ) {};
-
-        //! access function
-        int dimension() const {
-                return dimy;
-        }
+class fnc : public root
+{
+        protected:
+                //! Length of the output vector
+                int dimy;
+        public:
+                //!default constructor
+                fnc() {};
+                //! function evaluates numerical value of \f$f(x)\f$ at \f$x=\f$ \c cond
+                virtual vec eval (const vec &cond) {
+                        return vec (0);
+                };
+
+                //! function substitutes given value into an appropriate position
+                virtual void condition (const vec &val) {};
+
+                //! access function
+                int dimension() const {
+                        return dimy;
+                }
 };
 
@@ -257 +261 @@
 //! Probability density function with numerical statistics, e.g. posterior density.
 
-class epdf : public root {
-protected:
-        //! dimension of the random variable
-        int dim;
-        //! Description of the random variable
-        RV rv;
-
-public:
-        /*! \name Constructors
-         Construction of each epdf should support two types of constructors:
-        \li empty constructor,
-        \li copy constructor,
-
-        The following constructors should be supported for convenience:
-        \li constructor followed by calling \c set_parameters()
-        \li constructor accepting random variables calling \c set_rv()
-
-         All internal data structures are constructed as empty. Their values (including sizes) will be set by method \c set_parameters(). This way references can be initialized in constructors.
-        @{*/
-        epdf() : dim ( 0 ), rv() {};
-        epdf ( const epdf &e ) : dim ( e.dim ), rv ( e.rv ) {};
-        epdf ( const RV &rv0 ) : dim ( rv0._dsize() ) {
-                set_rv ( rv0 );
-        };
-        void set_parameters ( int dim0 ) {
-                dim = dim0;
-        }
-        //!@}
-
-        //! \name Matematical Operations
-        //!@{
-
-        //! Returns a sample, \f$ x \f$ from density \f$ f_x()\f$
-        virtual vec sample() const {
-                it_error ( "not implemneted" );
-                return vec ( 0 );
-        };
-        //! Returns N samples, \f$ [x_1 , x_2 , \ldots \ \f$  from density \f$ f_x(rv)\f$
-        virtual mat sample_m ( int N ) const;
-        //! Compute log-probability of argument \c val
-        //! In case the argument is out of suport return -Infinity
-        virtual double evallog ( const vec &val ) const {
-                it_error ( "not implemneted" );
-                return 0.0;
-        };
-        //! Compute log-probability of multiple values argument \c val
-        virtual vec evallog_m ( const mat &Val ) const {
-                vec x ( Val.cols() );
-                for ( int i = 0; i < Val.cols(); i++ ) {
-                        x ( i ) = evallog ( Val.get_col ( i ) ) ;
-                }
-                return x;
-        }
-        //! Compute log-probability of multiple values argument \c val
-        virtual vec evallog_m ( const Array<vec> &Avec ) const {
-                vec x ( Avec.size() );
-                for ( int i = 0; i < Avec.size(); i++ ) {
-                        x ( i ) = evallog ( Avec ( i ) ) ;
-                }
-                return x;
-        }
-        //! Return conditional density on the given RV, the remaining rvs will be in conditioning
-        virtual mpdf* condition ( const RV &rv ) const  {
-                it_warning ( "Not implemented" );
-                return NULL;
-        }
-
-        //! Return marginal density on the given RV, the remainig rvs are intergrated out
-        virtual epdf* marginal ( const RV &rv ) const {
-                it_warning ( "Not implemented" );
-                return NULL;
-        }
-
-        //! return expected value
-        virtual vec mean() const {
-                it_error ( "not implemneted" );
-                return vec ( 0 );
-        };
-
-        //! return expected variance (not covariance!)
-        virtual vec variance() const {
-                it_error ( "not implemneted" );
-                return vec ( 0 );
-        };
-        //! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default
-        virtual void qbounds ( vec &lb, vec &ub, double percentage = 0.95 ) const {
-                vec mea = mean();
-                vec std = sqrt ( variance() );
-                lb = mea - 2 * std;
-                ub = mea + 2 * std;
-        };
-        //!@}
-
-        //! \name Connection to other classes
-        //! Description of the random quantity via attribute \c rv is optional.
-        //! For operations such as sampling \c rv does not need to be set. However, for \c marginalization
-        //! and \c conditioning \c rv has to be set. NB:
-        //! @{
-
-        //!Name its rv
-        void set_rv ( const RV &rv0 ) {
-                rv = rv0;
-        }   //it_assert_debug(isnamed(),""); };
-        //! True if rv is assigned
-        bool isnamed() const {
-                bool b = ( dim == rv._dsize() );
-                return b;
-        }
-        //! Return name (fails when isnamed is false)
-        const RV& _rv() const {
-                it_assert_debug ( isnamed(), "" );
-                return rv;
-        }
-        //!@}
-
-        //! \name Access to attributes
-        //! @{
-
-        //! Size of the random variable
-        int dimension() const {
-                return dim;
-        }
-        //! Load from structure with elements:
-        //!  \code
-        //! { rv = {class="RV", names=(...),}; // RV describing meaning of random variable
-        //!   // elements of offsprings
-        //! }
-        //! \endcode
-        //!@}
-        void from_setting ( const Setting &set ) {
-                RV* r = UI::build<RV> ( set, "rv" );
-                if ( r ) {
-                        set_rv ( *r );
-                        delete r;
-                }
-        }
+class epdf : public root
+{
+        protected:
+                //! dimension of the random variable
+                int dim;
+                //! Description of the random variable
+                RV rv;
+
+        public:
+                /*! \name Constructors
+                 Construction of each epdf should support two types of constructors:
+                \li empty constructor,
+                \li copy constructor,
+
+                The following constructors should be supported for convenience:
+                \li constructor followed by calling \c set_parameters()
+                \li constructor accepting random variables calling \c set_rv()
+
+                 All internal data structures are constructed as empty. Their values (including sizes) will be set by method \c set_parameters(). This way references can be initialized in constructors.
+                @{*/
+                epdf() : dim (0), rv() {};
+                epdf (const epdf &e) : dim (e.dim), rv (e.rv) {};
+                epdf (const RV &rv0) : dim (rv0._dsize()) {
+                        set_rv (rv0);
+                };
+                void set_parameters (int dim0) {
+                        dim = dim0;
+                }
+                //!@}
+
+                //! \name Matematical Operations
+                //!@{
+
+                //! Returns a sample, \f$ x \f$ from density \f$ f_x()\f$
+                virtual vec sample() const {
+                        it_error ("not implemneted");
+                        return vec (0);
+                };
+                //! Returns N samples, \f$ [x_1 , x_2 , \ldots \ \f$  from density \f$ f_x(rv)\f$
+                virtual mat sample_m (int N) const;
+                //! Compute log-probability of argument \c val
+                //! In case the argument is out of suport return -Infinity
+                virtual double evallog (const vec &val) const {
+                        it_error ("not implemneted");
+                        return 0.0;
+                };
+                //! Compute log-probability of multiple values argument \c val
+                virtual vec evallog_m (const mat &Val) const {
+                        vec x (Val.cols());
+                        for (int i = 0; i < Val.cols(); i++) {
+                                x (i) = evallog (Val.get_col (i)) ;
+                        }
+                        return x;
+                }
+                //! Compute log-probability of multiple values argument \c val
+                virtual vec evallog_m (const Array<vec> &Avec) const {
+                        vec x (Avec.size());
+                        for (int i = 0; i < Avec.size(); i++) {
+                                x (i) = evallog (Avec (i)) ;
+                        }
+                        return x;
+                }
+                //! Return conditional density on the given RV, the remaining rvs will be in conditioning
+                virtual mpdf* condition (const RV &rv) const  {
+                        it_warning ("Not implemented");
+                        return NULL;
+                }
+
+                //! Return marginal density on the given RV, the remainig rvs are intergrated out
+                virtual epdf* marginal (const RV &rv) const {
+                        it_warning ("Not implemented");
+                        return NULL;
+                }
+
+                //! return expected value
+                virtual vec mean() const {
+                        it_error ("not implemneted");
+                        return vec (0);
+                };
+
+                //! return expected variance (not covariance!)
+                virtual vec variance() const {
+                        it_error ("not implemneted");
+                        return vec (0);
+                };
+                //! Lower and upper bounds of \c percentage % quantile, returns mean-2*sigma as default
+                virtual void qbounds (vec &lb, vec &ub, double percentage = 0.95) const {
+                        vec mea = mean();
+                        vec std = sqrt (variance());
+                        lb = mea - 2 * std;
+                        ub = mea + 2 * std;
+                };
+                //!@}
+
+                //! \name Connection to other classes
+                //! Description of the random quantity via attribute \c rv is optional.
+                //! For operations such as sampling \c rv does not need to be set. However, for \c marginalization
+                //! and \c conditioning \c rv has to be set. NB:
+                //! @{
+
+                //!Name its rv
+                void set_rv (const RV &rv0) {
+                        rv = rv0;
+                }   //it_assert_debug(isnamed(),""); };
+                //! True if rv is assigned
+                bool isnamed() const {
+                        bool b = (dim == rv._dsize());
+                        return b;
+                }
+                //! Return name (fails when isnamed is false)
+                const RV& _rv() const {
+                        it_assert_debug (isnamed(), "");
+                        return rv;
+                }
+                //!@}
+
+                //! \name Access to attributes
+                //! @{
+
+                //! Size of the random variable
+                int dimension() const {
+                        return dim;
+                }
+                //! Load from structure with elements:
+                //!  \code
+                //! { rv = {class="RV", names=(...),}; // RV describing meaning of random variable
+                //!   // elements of offsprings
+                //! }
+                //! \endcode
+                //!@}
+                void from_setting (const Setting &set) {
+                        RV* r = UI::build<RV> (set, "rv");
+                        if (r) {
+                                set_rv (*r);
+                                delete r;
+                        }
+                }
 
 };
@@ -399 +404 @@
 //! Conditional probability density, e.g. modeling some dependencies.
 //TODO Samplecond can be generalized
-class mpdf : public root {
-protected:
-        //!dimension of the condition
-        int dimc;
-        //! random variable in condition
-        RV rvc;
-private:
-        //! internal epdf, used only as cache to avoid virtual calls of \c _rv() and \c _dimension()
-        epdf* ep;
-        
-protected:
-        void set_ep(epdf &iepdf) {
-                ep = &iepdf;
-        }
-
-public:
-        //! \name Constructors
-        //! @{
-
-        mpdf() : dimc ( 0 ), rvc(), ep(NULL) { }
-
-        mpdf ( const mpdf &m ) : dimc ( m.dimc ), rvc ( m.rvc ), ep ( m.ep ) { }
-        //!@}
-
-        //! \name Matematical operations
-        //!@{
-
-        //! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv
-        virtual vec samplecond ( const vec &cond ){it_error("Not implemented");return vec(0);};
-
-        //! Returns \param N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv
-        virtual mat samplecond_m ( const vec &cond, int N ){it_error("Not implemented");return mat();}
-
-
-        //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
-        virtual double evallogcond ( const vec &dt, const vec &cond ){it_error("Not implemented");return 0.0;}
-
-        //! Matrix version of evallogcond
-        virtual vec evallogcond_m ( const mat &Dt, const vec &cond ){it_error("Not implemented");return vec(0);}
-
-        //! Array<vec> version of evallogcond
-        virtual vec evallogcond_m ( const Array<vec> &Dt, const vec &cond ){it_error("Not implemented");return vec(0);}
-
-        //! \name Access to attributes
-        //! @{
-
-        RV _rv() {
-                return ep->_rv();
-        }
-        RV _rvc() {
-                it_assert_debug ( isnamed(), "" );
-                return rvc;
-        }
-        int dimension() {
-                return ep->dimension();
-        }
-        int dimensionc() {
-                return dimc;
-        }
-
-        //! Load from structure with elements:
-        //!  \code
-        //! { class = "mpdf_offspring",
-        //!   rv = {class="RV", names=(...),}; // RV describing meaning of random variable
-        //!   rvc= {class="RV", names=(...),}; // RV describing meaning of random variable in condition
-        //!   // elements of offsprings
-        //! }
-        //! \endcode
-        //!@}
-        void from_setting ( const Setting &set );
-        //!@}
-
-        //! \name Connection to other objects
-        //!@{
-        void set_rvc ( const RV &rvc0 ) {
-                rvc = rvc0;
-        }
-        void set_rv ( const RV &rv0 ) {
-                ep->set_rv ( rv0 );
-        }
-        bool isnamed() {
-                return ( ep->isnamed() ) && ( dimc == rvc._dsize() );
-        }
-        //!@}
+class mpdf : public root
+{
+        protected:
+                //!dimension of the condition
+                int dimc;
+                //! random variable in condition
+                RV rvc;
+        private:
+                //! internal epdf, used only as cache to avoid virtual calls of \c _rv() and \c _dimension()
+                epdf* ep;
+
+        protected:
+                void set_ep (epdf &iepdf) {
+                        ep = &iepdf;
+                }
+
+        public:
+                //! \name Constructors
+                //! @{
+
+                mpdf() : dimc (0), rvc(), ep (NULL) { }
+
+                mpdf (const mpdf &m) : dimc (m.dimc), rvc (m.rvc), ep (m.ep) { }
+                //!@}
+
+                //! \name Matematical operations
+                //!@{
+
+                //! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv
+                virtual vec samplecond (const vec &cond) {it_error ("Not implemented");return vec (0);};
+
+                //! Returns \param N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv
+                virtual mat samplecond_m (const vec &cond, int N) {it_error ("Not implemented");return mat();}
+
+
+                //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
+                virtual double evallogcond (const vec &dt, const vec &cond) {it_error ("Not implemented");return 0.0;}
+
+                //! Matrix version of evallogcond
+                virtual vec evallogcond_m (const mat &Dt, const vec &cond) {it_error ("Not implemented");return vec (0);}
+
+                //! Array<vec> version of evallogcond
+                virtual vec evallogcond_m (const Array<vec> &Dt, const vec &cond) {it_error ("Not implemented");return vec (0);}
+
+                //! \name Access to attributes
+                //! @{
+
+                RV _rv() {
+                        return ep->_rv();
+                }
+                RV _rvc() {
+                        it_assert_debug (isnamed(), "");
+                        return rvc;
+                }
+                int dimension() {
+                        return ep->dimension();
+                }
+                int dimensionc() {
+                        return dimc;
+                }
+
+                //! Load from structure with elements:
+                //!  \code
+                //! { class = "mpdf_offspring",
+                //!   rv = {class="RV", names=(...),}; // RV describing meaning of random variable
+                //!   rvc= {class="RV", names=(...),}; // RV describing meaning of random variable in condition
+                //!   // elements of offsprings
+                //! }
+                //! \endcode
+                //!@}
+                void from_setting (const Setting &set);
+                //!@}
+
+                //! \name Connection to other objects
+                //!@{
+                void set_rvc (const RV &rvc0) {
+                        rvc = rvc0;
+                }
+                void set_rv (const RV &rv0) {
+                        ep->set_rv (rv0);
+                }
+                bool isnamed() {
+                        return (ep->isnamed()) && (dimc == rvc._dsize());
+                }
+                //!@}
 };
 
 template <class EPDF>
-class mpdf_internal: public mpdf{
+class mpdf_internal: public mpdf
+{
         protected :
                 EPDF iepdf;
         public:
                 //! constructor
-                mpdf_internal(): mpdf(),iepdf(){set_ep(iepdf);}
+                mpdf_internal() : mpdf(), iepdf() {set_ep (iepdf);}
                 //! Update \c iepdf so that it represents this mpdf conditioned on \c rvc = cond
                 //! This function provides convenient reimplementation in offsprings
-                        virtual void condition ( const vec &cond ) {
-                                it_error ( "Not implemented" );
-                        };
+                virtual void condition (const vec &cond) {
+                        it_error ("Not implemented");
+                };
                 //!access function to iepdf
-                EPDF& e(){return iepdf;}
-                        
+                EPDF& e() {return iepdf;}
+
                 //! Reimplements samplecond using \c condition()
-                vec samplecond ( const vec &cond );
+                vec samplecond (const vec &cond);
                 //! Reimplements evallogcond using \c condition()
-                double evallogcond ( const vec &val, const vec &cond );
-                //! Efficient version of evallogcond for matrices 
-                virtual vec evallogcond_m ( const mat &Dt, const vec &cond );
+                double evallogcond (const vec &val, const vec &cond);
+                //! Efficient version of evallogcond for matrices
+                virtual vec evallogcond_m (const mat &Dt, const vec &cond);
                 //! Efficient version of evallogcond for Array<vec>
-                virtual vec evallogcond_m ( const Array<vec> &Dt, const vec &cond );
-                //! Efficient version of samplecond 
-                virtual mat samplecond_m ( const vec &cond, int N );
+                virtual vec evallogcond_m (const Array<vec> &Dt, const vec &cond);
+                //! Efficient version of samplecond
+                virtual mat samplecond_m (const vec &cond, int N);
 };
 
@@ -537 +544 @@
 
 */
-class datalink {
-protected:
-        //! Remember how long val should be
-        int downsize;
-
-        //! Remember how long val of "Up" should be
-        int upsize;
-
-        //! val-to-val link, indices of the upper val
-        ivec v2v_up;
-
-public:
-        //! Constructor
-        datalink() : downsize ( 0 ), upsize ( 0 ) { }
-        datalink ( const RV &rv, const RV &rv_up ) {
-                set_connection ( rv, rv_up );
-        }
-
-        //! set connection, rv must be fully present in rv_up
-        void set_connection ( const RV &rv, const RV &rv_up ) {
-                downsize = rv._dsize();
-                upsize = rv_up._dsize();
-                v2v_up = rv.dataind ( rv_up );
-
-                it_assert_debug ( v2v_up.length() == downsize, "rv is not fully in rv_up" );
-        }
-
-        //! set connection using indices
-        void set_connection ( int ds, int us, const ivec &upind ) {
-                downsize = ds;
-                upsize = us;
-                v2v_up = upind;
-
-                it_assert_debug ( v2v_up.length() == downsize, "rv is not fully in rv_up" );
-        }
-
-        //! Get val for myself from val of "Up"
-        vec pushdown ( const vec &val_up ) {
-                it_assert_debug ( upsize == val_up.length(), "Wrong val_up" );
-                return get_vec ( val_up, v2v_up );
-        }
-
-        //! Fill val of "Up" by my pieces
-        void pushup ( vec &val_up, const vec &val ) {
-                it_assert_debug ( downsize == val.length(), "Wrong val" );
-                it_assert_debug ( upsize == val_up.length(), "Wrong val_up" );
-                set_subvector ( val_up, v2v_up, val );
-        }
+class datalink
+{
+        protected:
+                //! Remember how long val should be
+                int downsize;
+
+                //! Remember how long val of "Up" should be
+                int upsize;
+
+                //! val-to-val link, indices of the upper val
+                ivec v2v_up;
+
+        public:
+                //! Constructor
+                datalink() : downsize (0), upsize (0) { }
+                datalink (const RV &rv, const RV &rv_up) {
+                        set_connection (rv, rv_up);
+                }
+
+                //! set connection, rv must be fully present in rv_up
+                void set_connection (const RV &rv, const RV &rv_up) {
+                        downsize = rv._dsize();
+                        upsize = rv_up._dsize();
+                        v2v_up = rv.dataind (rv_up);
+
+                        it_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up");
+                }
+
+                //! set connection using indices
+                void set_connection (int ds, int us, const ivec &upind) {
+                        downsize = ds;
+                        upsize = us;
+                        v2v_up = upind;
+
+                        it_assert_debug (v2v_up.length() == downsize, "rv is not fully in rv_up");
+                }
+
+                //! Get val for myself from val of "Up"
+                vec pushdown (const vec &val_up) {
+                        it_assert_debug (upsize == val_up.length(), "Wrong val_up");
+                        return get_vec (val_up, v2v_up);
+                }
+
+                //! Fill val of "Up" by my pieces
+                void pushup (vec &val_up, const vec &val) {
+                        it_assert_debug (downsize == val.length(), "Wrong val");
+                        it_assert_debug (upsize == val_up.length(), "Wrong val_up");
+                        set_subvector (val_up, v2v_up, val);
+                }
 };
 
 //! Data link with a condition.
-class datalink_m2e: public datalink {
-protected:
-        //! Remember how long cond should be
-        int condsize;
-
-        //!upper_val-to-local_cond link, indices of the upper val
-        ivec v2c_up;
-
-        //!upper_val-to-local_cond link, indices of the local cond
-        ivec v2c_lo;
-
-public:
-        //! Constructor
-        datalink_m2e() : condsize ( 0 ) { }
-
-        void set_connection ( const RV &rv, const RV &rvc, const RV &rv_up ) {
-                datalink::set_connection ( rv, rv_up );
-                condsize = rvc._dsize();
-                //establish v2c connection
-                rvc.dataind ( rv_up, v2c_lo, v2c_up );
-        }
-
-        //!Construct condition
-        vec get_cond ( const vec &val_up ) {
-                vec tmp ( condsize );
-                set_subvector ( tmp, v2c_lo, val_up ( v2c_up ) );
-                return tmp;
-        }
-
-        void pushup_cond ( vec &val_up, const vec &val, const vec &cond ) {
-                it_assert_debug ( downsize == val.length(), "Wrong val" );
-                it_assert_debug ( upsize == val_up.length(), "Wrong val_up" );
-                set_subvector ( val_up, v2v_up, val );
-                set_subvector ( val_up, v2c_up, cond );
-        }
+class datalink_m2e: public datalink
+{
+        protected:
+                //! Remember how long cond should be
+                int condsize;
+
+                //!upper_val-to-local_cond link, indices of the upper val
+                ivec v2c_up;
+
+                //!upper_val-to-local_cond link, indices of the local cond
+                ivec v2c_lo;
+
+        public:
+                //! Constructor
+                datalink_m2e() : condsize (0) { }
+
+                void set_connection (const RV &rv, const RV &rvc, const RV &rv_up) {
+                        datalink::set_connection (rv, rv_up);
+                        condsize = rvc._dsize();
+                        //establish v2c connection
+                        rvc.dataind (rv_up, v2c_lo, v2c_up);
+                }
+
+                //!Construct condition
+                vec get_cond (const vec &val_up) {
+                        vec tmp (condsize);
+                        set_subvector (tmp, v2c_lo, val_up (v2c_up));
+                        return tmp;
+                }
+
+                void pushup_cond (vec &val_up, const vec &val, const vec &cond) {
+                        it_assert_debug (downsize == val.length(), "Wrong val");
+                        it_assert_debug (upsize == val_up.length(), "Wrong val_up");
+                        set_subvector (val_up, v2v_up, val);
+                        set_subvector (val_up, v2c_up, cond);
+                }
 };
 
 //!DataLink is a connection between mpdf and its superordinate (Up)
 //! This class links
-class datalink_m2m: public datalink_m2e {
-protected:
-        //!cond-to-cond link, indices of the upper cond
-        ivec c2c_up;
-        //!cond-to-cond link, indices of the local cond
-        ivec c2c_lo;
-
-public:
-        //! Constructor
-        datalink_m2m() {};
-        void set_connection ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) {
-                datalink_m2e::set_connection ( rv, rvc, rv_up );
-                //establish c2c connection
-                rvc.dataind ( rvc_up, c2c_lo, c2c_up );
-                it_assert_debug ( c2c_lo.length() + v2c_lo.length() == condsize, "cond is not fully given" );
-        }
-
-        //! Get cond for myself from val and cond of "Up"
-        vec get_cond ( const vec &val_up, const vec &cond_up ) {
-                vec tmp ( condsize );
-                set_subvector ( tmp, v2c_lo, val_up ( v2c_up ) );
-                set_subvector ( tmp, c2c_lo, cond_up ( c2c_up ) );
-                return tmp;
-        }
-        //! Fill
+class datalink_m2m: public datalink_m2e
+{
+        protected:
+                //!cond-to-cond link, indices of the upper cond
+                ivec c2c_up;
+                //!cond-to-cond link, indices of the local cond
+                ivec c2c_lo;
+
+        public:
+                //! Constructor
+                datalink_m2m() {};
+                void set_connection (const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up) {
+                        datalink_m2e::set_connection (rv, rvc, rv_up);
+                        //establish c2c connection
+                        rvc.dataind (rvc_up, c2c_lo, c2c_up);
+                        it_assert_debug (c2c_lo.length() + v2c_lo.length() == condsize, "cond is not fully given");
+                }
+
+                //! Get cond for myself from val and cond of "Up"
+                vec get_cond (const vec &val_up, const vec &cond_up) {
+                        vec tmp (condsize);
+                        set_subvector (tmp, v2c_lo, val_up (v2c_up));
+                        set_subvector (tmp, c2c_lo, cond_up (c2c_up));
+                        return tmp;
+                }
+                //! Fill
 
 };
@@ -660 +670 @@
 This class abstracts logging of results from implementation. This class replaces direct logging of results (e.g. to files or to global variables) by calling methods of a logger. Specializations of this abstract class for specific storage method are designed.
  */
-class logger : public root {
-protected:
-        //! RVs of all logged variables.
-        Array<RV> entries;
-        //! Names of logged quantities, e.g. names of algorithm variants
-        Array<string> names;
-public:
-        //!Default constructor
-        logger() : entries ( 0 ), names ( 0 ) {}
-
-        //! returns an identifier which will be later needed for calling the \c logit() function
-        //! For empty RV it returns -1, this entry will be ignored by \c logit().
-        virtual int add ( const RV &rv, string prefix = "" ) {
-                int id;
-                if ( rv._dsize() > 0 ) {
-                        id = entries.length();
-                        names = concat ( names, prefix ); // diff
-                        entries.set_length ( id + 1, true );
-                        entries ( id ) = rv;
-                } else {
-                        id = -1;
-                }
-                return id; // identifier of the last entry
-        }
-
-        //! log this vector
-        virtual void logit ( int id, const vec &v ) = 0;
-        //! log this double
-        virtual void logit ( int id, const double &d ) = 0;
-
-        //! Shifts storage position for another time step.
-        virtual void step() = 0;
-
-        //! Finalize storing information
-        virtual void finalize() {};
-
-        //! Initialize the storage
-        virtual void init() {};
+class logger : public root
+{
+        protected:
+                //! RVs of all logged variables.
+                Array<RV> entries;
+                //! Names of logged quantities, e.g. names of algorithm variants
+                Array<string> names;
+        public:
+                //!Default constructor
+                logger() : entries (0), names (0) {}
+
+                //! returns an identifier which will be later needed for calling the \c logit() function
+                //! For empty RV it returns -1, this entry will be ignored by \c logit().
+                virtual int add (const RV &rv, string prefix = "") {
+                        int id;
+                        if (rv._dsize() > 0) {
+                                id = entries.length();
+                                names = concat (names, prefix);   // diff
+                                entries.set_length (id + 1, true);
+                                entries (id) = rv;
+                        } else {
+                                id = -1;
+                        }
+                        return id; // identifier of the last entry
+                }
+
+                //! log this vector
+                virtual void logit (int id, const vec &v) = 0;
+                //! log this double
+                virtual void logit (int id, const double &d) = 0;
+
+                //! Shifts storage position for another time step.
+                virtual void step() = 0;
+
+                //! Finalize storing information
+                virtual void finalize() {};
+
+                //! Initialize the storage
+                virtual void init() {};
 
 };
@@ -704 +715 @@
 
 */
-class mepdf : public mpdf {
-
-shared_ptr<epdf> ipdf;
-public:
-        //!Default constructor
-        mepdf() { }
-
-        mepdf ( shared_ptr<epdf> em ) {
-                ipdf = em;
-                set_ep(*ipdf.get());
-                dimc = 0;
-        }
-
-        //! empty
-        void condition ( const vec &cond );
-
-        //! Load from structure with elements:
-        //!  \code
-        //! { class = "mepdf",
-        //!   epdf = {class="epdf_offspring",...}
-        //! }
-        //! \endcode
-        //!@}
-        void from_setting ( const Setting &set );
-};
-UIREGISTER ( mepdf );
+class mepdf : public mpdf
+{
+
+                shared_ptr<epdf> ipdf;
+        public:
+                //!Default constructor
+                mepdf() { }
+
+                mepdf (shared_ptr<epdf> em) {
+                        ipdf = em;
+                        set_ep (*ipdf.get());
+                        dimc = 0;
+                }
+
+                //! empty
+                void condition (const vec &cond);
+
+                //! Load from structure with elements:
+                //!  \code
+                //! { class = "mepdf",
+                //!   epdf = {class="epdf_offspring",...}
+                //! }
+                //! \endcode
+                //!@}
+                void from_setting (const Setting &set);
+};
+UIREGISTER (mepdf);
 
 //!\brief Chain rule of pdfs - abstract part common for mprod and merger.
 //!this abstract class is common to epdf and mpdf
 //!\todo Think of better design - global functions rv=get_rv(Array<mpdf*> mpdfs); ??
-class compositepdf {
-protected:
-        //! Elements of composition
-        Array<mpdf*> mpdfs;
-        bool owning_mpdfs;
-public:
-        compositepdf() : mpdfs ( 0 ) {};
-        compositepdf ( Array<mpdf*> A0, bool own = false ) {
-                set_elements ( A0, own );
-        };
-        void set_elements ( Array<mpdf*> A0, bool own = false ) {
-                mpdfs = A0;
-                owning_mpdfs = own;
-        };
-        //! find common rv, flag \param checkoverlap modifies whether overlaps are acceptable
-        RV getrv ( bool checkoverlap = false );
-        //! common rvc of all mpdfs is written to rvc
-        void setrvc ( const RV &rv, RV &rvc );
-        ~compositepdf() {
-                if ( owning_mpdfs ) for ( int i = 0; i < mpdfs.length(); i++ ) {
-                                delete mpdfs ( i );
-                        }
-        };
+class compositepdf
+{
+        protected:
+                //! Elements of composition
+                Array<mpdf*> mpdfs;
+                bool owning_mpdfs;
+        public:
+                compositepdf() : mpdfs (0) {};
+                compositepdf (Array<mpdf*> A0, bool own = false) {
+                        set_elements (A0, own);
+                };
+                void set_elements (Array<mpdf*> A0, bool own = false) {
+                        mpdfs = A0;
+                        owning_mpdfs = own;
+                };
+                //! find common rv, flag \param checkoverlap modifies whether overlaps are acceptable
+                RV getrv (bool checkoverlap = false);
+                //! common rvc of all mpdfs is written to rvc
+                void setrvc (const RV &rv, RV &rvc);
+                ~compositepdf() {
+                        if (owning_mpdfs) for (int i = 0; i < mpdfs.length(); i++) {
+                                        delete mpdfs (i);
+                                }
+                };
 };
 
@@ -766 +779 @@
 */
 
-class DS : public root {
-protected:
-        int dtsize;
-        int utsize;
-        //!Description of data returned by \c getdata().
-        RV Drv;
-        //!Description of data witten by by \c write().
-        RV Urv; //
-        //! Remember its own index in Logger L
-        int L_dt, L_ut;
-public:
-        //! default constructors
-        DS() : Drv(), Urv() {};
-        //! Returns full vector of observed data=[output, input]
-        virtual void getdata ( vec &dt ) {
-                it_error ( "abstract class" );
-        };
-        //! Returns data records at indeces.
-        virtual void getdata ( vec &dt, const ivec &indeces ) {
-                it_error ( "abstract class" );
-        };
-        //! Accepts action variable and schedule it for application.
-        virtual void write ( vec &ut ) {
-                it_error ( "abstract class" );
-        };
-        //! Accepts action variables at specific indeces
-        virtual void write ( vec &ut, const ivec &indeces ) {
-                it_error ( "abstract class" );
-        };
-
-        //! Moves from \f$ t \f$ to \f$ t+1 \f$, i.e. perfroms the actions and reads response of the system.
-        virtual void step() = 0;
-
-        //! Register DS for logging into logger L
-        virtual void log_add ( logger &L ) {
-                it_assert_debug ( dtsize == Drv._dsize(), "" );
-                it_assert_debug ( utsize == Urv._dsize(), "" );
-
-                L_dt = L.add ( Drv, "" );
-                L_ut = L.add ( Urv, "" );
-        }
-        //! Register DS for logging into logger L
-        virtual void logit ( logger &L ) {
-                vec tmp ( Drv._dsize() + Urv._dsize() );
-                getdata ( tmp );
-                // d is first in getdata
-                L.logit ( L_dt, tmp.left ( Drv._dsize() ) );
-                // u follows after d in getdata
-                L.logit ( L_ut, tmp.mid ( Drv._dsize(), Urv._dsize() ) );
-        }
-        //!access function
-        virtual RV _drv() const {
-                return concat ( Drv, Urv );
-        }
-        //!access function
-        const RV& _urv() const {
-                return Urv;
-        }
-        //! set random rvariables
-        virtual void set_drv ( const  RV &drv, const RV &urv ) {
-                Drv = drv;
-                Urv = urv;
-        }
+class DS : public root
+{
+        protected:
+                int dtsize;
+                int utsize;
+                //!Description of data returned by \c getdata().
+                RV Drv;
+                //!Description of data witten by by \c write().
+                RV Urv; //
+                //! Remember its own index in Logger L
+                int L_dt, L_ut;
+        public:
+                //! default constructors
+                DS() : Drv(), Urv() {};
+                //! Returns full vector of observed data=[output, input]
+                virtual void getdata (vec &dt) {
+                        it_error ("abstract class");
+                };
+                //! Returns data records at indeces.
+                virtual void getdata (vec &dt, const ivec &indeces) {
+                        it_error ("abstract class");
+                };
+                //! Accepts action variable and schedule it for application.
+                virtual void write (vec &ut) {
+                        it_error ("abstract class");
+                };
+                //! Accepts action variables at specific indeces
+                virtual void write (vec &ut, const ivec &indeces) {
+                        it_error ("abstract class");
+                };
+
+                //! Moves from \f$ t \f$ to \f$ t+1 \f$, i.e. perfroms the actions and reads response of the system.
+                virtual void step() = 0;
+
+                //! Register DS for logging into logger L
+                virtual void log_add (logger &L) {
+                        it_assert_debug (dtsize == Drv._dsize(), "");
+                        it_assert_debug (utsize == Urv._dsize(), "");
+
+                        L_dt = L.add (Drv, "");
+                        L_ut = L.add (Urv, "");
+                }
+                //! Register DS for logging into logger L
+                virtual void logit (logger &L) {
+                        vec tmp (Drv._dsize() + Urv._dsize());
+                        getdata (tmp);
+                        // d is first in getdata
+                        L.logit (L_dt, tmp.left (Drv._dsize()));
+                        // u follows after d in getdata
+                        L.logit (L_ut, tmp.mid (Drv._dsize(), Urv._dsize()));
+                }
+                //!access function
+                virtual RV _drv() const {
+                        return concat (Drv, Urv);
+                }
+                //!access function
+                const RV& _urv() const {
+                        return Urv;
+                }
+                //! set random rvariables
+                virtual void set_drv (const  RV &drv, const RV &urv) {
+                        Drv = drv;
+                        Urv = urv;
+                }
 };
 
@@ -852 +866 @@
 */
 
-class BM : public root {
-protected:
-        //! Random variable of the data (optional)
-        RV drv;
-        //!Logarithm of marginalized data likelihood.
-        double ll;
-        //!  If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save computational time.
-        bool evalll;
-public:
-        //! \name Constructors
-        //! @{
-
-        BM() : ll ( 0 ), evalll ( true ), LIDs ( 4 ), LFlags ( 4 ) {
-                LIDs = -1;/*empty IDs*/
-                LFlags = 0;
-                LFlags ( 0 ) = 1;/*log only mean*/
-        };
-        BM ( const BM &B ) :  drv ( B.drv ), ll ( B.ll ), evalll ( B.evalll ) {}
-        //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
-        //! Prototype: \code BM* _copy_() const {return new BM(*this);} \endcode
-        virtual BM* _copy_() const {
-                return NULL;
-        };
-        //!@}
-
-        //! \name Mathematical operations
-        //!@{
-
-        /*! \brief Incremental Bayes rule
-        @param dt vector of input data
-        */
-        virtual void bayes ( const vec &dt ) = 0;
-        //! Batch Bayes rule (columns of Dt are observations)
-        virtual void bayesB ( const mat &Dt );
-        //! Evaluates predictive log-likelihood of the given data record
-        //! I.e. marginal likelihood of the data with the posterior integrated out.
-        virtual double logpred ( const vec &dt ) const {
-                it_error ( "Not implemented" );
-                return 0.0;
+class BM : public root
+{
+        protected:
+                //! Random variable of the data (optional)
+                RV drv;
+                //!Logarithm of marginalized data likelihood.
+                double ll;
+                //!  If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save computational time.
+                bool evalll;
+        public:
+                //! \name Constructors
+                //! @{
+
+                BM() : ll (0), evalll (true), LIDs (4), LFlags (4) {
+                        LIDs = -1;/*empty IDs*/
+                        LFlags = 0;
+                        LFlags (0) = 1;  /*log only mean*/
+                };
+                BM (const BM &B) :  drv (B.drv), ll (B.ll), evalll (B.evalll) {}
+                //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
+                //! Prototype: \code BM* _copy_() const {return new BM(*this);} \endcode
+                virtual BM* _copy_() const {
+                        return NULL;
+                };
+                //!@}
+
+                //! \name Mathematical operations
+                //!@{
+
+                /*! \brief Incremental Bayes rule
+                @param dt vector of input data
+                */
+                virtual void bayes (const vec &dt) = 0;
+                //! Batch Bayes rule (columns of Dt are observations)
+                virtual void bayesB (const mat &Dt);
+                //! Evaluates predictive log-likelihood of the given data record
+                //! I.e. marginal likelihood of the data with the posterior integrated out.
+                virtual double logpred (const vec &dt) const {
+                        it_error ("Not implemented");
+                        return 0.0;
+                }
+                //! Matrix version of logpred
+                vec logpred_m (const mat &dt) const {
+                        vec tmp (dt.cols());
+                        for (int i = 0; i < dt.cols(); i++) {
+                                tmp (i) = logpred (dt.get_col (i));
+                        }
+                        return tmp;
+                }
+
+                //!Constructs a predictive density \f$ f(d_{t+1} |d_{t}, \ldots d_{0}) \f$
+                virtual epdf* epredictor() const {
+                        it_error ("Not implemented");
+                        return NULL;
+                };
+                //!Constructs a conditional density 1-step ahead predictor \f$ f(d_{t+1} |d_{t+h-1}, \ldots d_{t})
+                virtual mpdf* predictor() const {
+                        it_error ("Not implemented");
+                        return NULL;
+                };
+                //!@}
+
+                //! \name Extension to conditional BM
+                //! This extension is useful e.g. in Marginalized Particle Filter (\ref bdm::MPF).
+                //! Alternatively, it can be used for automated connection to DS when the condition is observed
+                //!@{
+
+                //! Name of extension variable
+                RV rvc;
+                //! access function
+                const RV& _rvc() const {
+                        return rvc;
+                }
+
+                //! Substitute \c val for \c rvc.
+                virtual void condition (const vec &val) {
+                        it_error ("Not implemented!");
+                };
+
+                //!@}
+
+
+                //! \name Access to attributes
+                //!@{
+
+                const RV& _drv() const {
+                        return drv;
+                }
+                void set_drv (const RV &rv) {
+                        drv = rv;
+                }
+                void set_rv (const RV &rv) {
+                        const_cast<epdf&> (posterior()).set_rv (rv);
+                }
+                double _ll() const {
+                        return ll;
+                }
+                void set_evalll (bool evl0) {
+                        evalll = evl0;
+                }
+                virtual const epdf& posterior() const = 0;
+                virtual const epdf* _e() const = 0;
+                //!@}
+
+                //! \name Logging of results
+                //!@{
+
+                //! Set boolean options from a string, recognized are: "logbounds,logll"
+                virtual void set_options (const string &opt) {
+                        LFlags (0) = 1;
+                        if (opt.find ("logbounds") != string::npos) {
+                                LFlags (1) = 1;
+                                LFlags (2) = 1;
+                        }
+                        if (opt.find ("logll") != string::npos) {
+                                LFlags (3) = 1;
+                        }
+                }
+                //! IDs of storages in loggers 4:[1=mean,2=lb,3=ub,4=ll]
+                ivec LIDs;
+
+                //! Flags for logging - same size as LIDs, each entry correspond to the same in LIDs
+                ivec LFlags;
+                //! Add all logged variables to a logger
+                virtual void log_add (logger &L, const string &name = "") {
+                        // internal
+                        RV r;
+                        if (posterior().isnamed()) {
+                                r = posterior()._rv();
+                        } else {
+                                r = RV ("est", posterior().dimension());
+                        };
+
+                        // Add mean value
+                        if (LFlags (0)) LIDs (0) = L.add (r, name + "mean_");
+                        if (LFlags (1)) LIDs (1) = L.add (r, name + "lb_");
+                        if (LFlags (2)) LIDs (2) = L.add (r, name + "ub_");
+                        if (LFlags (3)) LIDs (3) = L.add (RV ("ll", 1), name);              //TODO: "local" RV
+                }
+                virtual void logit (logger &L) {
+                        L.logit (LIDs (0), posterior().mean());
+                        if (LFlags (1) || LFlags (2)) {        //if one of them is off, its LID==-1 and will not be stored
+                                vec ub, lb;
+                                posterior().qbounds (lb, ub);
+                                L.logit (LIDs (1), lb);
+                                L.logit (LIDs (2), ub);
+                        }
+                        if (LFlags (3)) L.logit (LIDs (3), ll);
+                }
+                //!@}
+};
+
+template<class EPDF>
+vec mpdf_internal<EPDF>::samplecond (const vec &cond)
+{
+        condition (cond);
+        vec temp = iepdf.sample();
+        return temp;
+}
+
+template<class EPDF>
+mat mpdf_internal<EPDF>::samplecond_m (const vec &cond, int N)
+{
+        condition (cond);
+        mat temp (dimension(), N);
+        vec smp (dimension());
+        for (int i = 0; i < N; i++) {
+                smp = iepdf.sample();
+                temp.set_col (i, smp);
         }
-        //! Matrix version of logpred
-        vec logpred_m ( const mat &dt ) const {
-                vec tmp ( dt.cols() );
-                for ( int i = 0; i < dt.cols(); i++ ) {
-                        tmp ( i ) = logpred ( dt.get_col ( i ) );
-                }
-                return tmp;
-        }
-
-        //!Constructs a predictive density \f$ f(d_{t+1} |d_{t}, \ldots d_{0}) \f$
-        virtual epdf* epredictor() const {
-                it_error ( "Not implemented" );
-                return NULL;
-        };
-        //!Constructs a conditional density 1-step ahead predictor \f$ f(d_{t+1} |d_{t+h-1}, \ldots d_{t})
-        virtual mpdf* predictor() const {
-                it_error ( "Not implemented" );
-                return NULL;
-        };
-        //!@}
-
-        //! \name Extension to conditional BM
-        //! This extension is useful e.g. in Marginalized Particle Filter (\ref bdm::MPF).
-        //! Alternatively, it can be used for automated connection to DS when the condition is observed
-        //!@{
-
-        //! Name of extension variable
-        RV rvc;
-        //! access function
-        const RV& _rvc() const {
-                return rvc;
-        }
-
-        //! Substitute \c val for \c rvc.
-        virtual void condition ( const vec &val ) {
-                it_error ( "Not implemented!" );
-        };
-
-        //!@}
-
-
-        //! \name Access to attributes
-        //!@{
-
-        const RV& _drv() const {
-                return drv;
-        }
-        void set_drv ( const RV &rv ) {
-                drv = rv;
-        }
-        void set_rv ( const RV &rv ) {
-                const_cast<epdf&> ( posterior() ).set_rv ( rv );
-        }
-        double _ll() const {
-                return ll;
-        }
-        void set_evalll ( bool evl0 ) {
-                evalll = evl0;
-        }
-        virtual const epdf& posterior() const = 0;
-        virtual const epdf* _e() const = 0;
-        //!@}
-
-        //! \name Logging of results
-        //!@{
-
-        //! Set boolean options from a string, recognized are: "logbounds,logll"
-        virtual void set_options ( const string &opt ) {
-                LFlags ( 0 ) = 1;
-                if ( opt.find ( "logbounds" ) != string::npos ) {
-                        LFlags ( 1 ) = 1;
-                        LFlags ( 2 ) = 1;
-                }
-                if ( opt.find ( "logll" ) != string::npos ) {
-                        LFlags ( 3 ) = 1;
-                }
-        }
-        //! IDs of storages in loggers 4:[1=mean,2=lb,3=ub,4=ll]
-        ivec LIDs;
-
-        //! Flags for logging - same size as LIDs, each entry correspond to the same in LIDs
-        ivec LFlags;
-        //! Add all logged variables to a logger
-        virtual void log_add ( logger &L, const string &name = "" ) {
-                // internal
-                RV r;
-                if ( posterior().isnamed() ) {
-                        r = posterior()._rv();
-                } else {
-                        r = RV ( "est", posterior().dimension() );
-                };
-
-                // Add mean value
-                if ( LFlags ( 0 ) ) LIDs ( 0 ) = L.add ( r, name + "mean_" );
-                if ( LFlags ( 1 ) ) LIDs ( 1 ) = L.add ( r, name + "lb_" );
-                if ( LFlags ( 2 ) ) LIDs ( 2 ) = L.add ( r, name + "ub_" );
-                if ( LFlags ( 3 ) ) LIDs ( 3 ) = L.add ( RV ( "ll", 1 ), name );    //TODO: "local" RV
-        }
-        virtual void logit ( logger &L ) {
-                L.logit ( LIDs ( 0 ), posterior().mean() );
-                if ( LFlags ( 1 ) || LFlags ( 2 ) ) {  //if one of them is off, its LID==-1 and will not be stored
-                        vec ub, lb;
-                        posterior().qbounds ( lb, ub );
-                        L.logit ( LIDs ( 1 ), lb );
-                        L.logit ( LIDs ( 2 ), ub );
-                }
-                if ( LFlags ( 3 ) ) L.logit ( LIDs ( 3 ), ll );
-        }
-        //!@}
-};
-
+
+        return temp;
+}
+
+template<class EPDF>
+double mpdf_internal<EPDF>::evallogcond (const vec &dt, const vec &cond)
+{
+        double tmp;
+        condition (cond);
+        tmp = iepdf.evallog (dt);
+        // it_assert_debug(std::isfinite(tmp), "Infinite value");
+        return tmp;
+}
+
+template<class EPDF>
+vec mpdf_internal<EPDF>::evallogcond_m (const mat &Dt, const vec &cond)
+{
+        condition (cond);
+        return iepdf.evallog_m (Dt);
+}
+
+template<class EPDF>
+vec mpdf_internal<EPDF>::evallogcond_m (const Array<vec> &Dt, const vec &cond)
+{
+        condition (cond);
+        return iepdf.evallog_m (Dt);
+}
 
 }; //namespace

library/bdm/estim/particles.h

@@ -67 +67 @@
                 resmethod = rm;
         };
-        void set_statistics ( const vec w0, epdf *epdf0 ) {
+        void set_statistics ( const vec w0, const epdf &epdf0 ) {
                 est.set_statistics ( w0, epdf0 );
         };
@@ -201 +201 @@
                 BMs.set_length ( n0 );
         }
-        void set_statistics ( epdf *epdf0, const BM_T* BMcond0 ) {
+        void set_statistics ( const epdf &epdf0, const BM_T* BMcond0 ) {
 
                 PF::set_statistics ( ones ( n ) / n, epdf0 );

library/bdm/stat/emix.cpp

@@ -199 +199 @@
         return aprgiw;
 };
+
+vec mmix::samplecond(const vec &cond) {
+        //Sample which component
+        vec cumDist = cumsum ( w );
+        double u0;
+#pragma omp critical
+        u0 = UniRNG.sample();
+
+        int i = 0;
+        while ( ( cumDist ( i ) < u0 ) && ( i < ( w.length() - 1 ) ) ) {
+                i++;
+        }
+
+        return Coms ( i )->samplecond(cond);
+}
 
 }

library/bdm/stat/emix.h

@@ -472 +472 @@
 
 
-/*! \brief Mixture of mpdfs with constant weights, all mpdfs are of equal type
+/*! \brief Mixture of mpdfs with constant weights, all mpdfs are of equal RV and RVC
 
 */
-// class mmix : public mpdf {
-// protected:
-//      //! Component (epdfs)
-//      Array<mpdf*> Coms;
-// 
-//      //!Internal epdf
-//      emix iepdf;
-// 
-// public:
-//      //!Default constructor
-//      mmix() : iepdf ( ) {
-//              set_ep ( iepdf );
-//      }
-// 
-//      //! Set weights \c w and components \c R
-//      void set_parameters ( const vec &w, const Array<mpdf*> &Coms ) {
-//              Array<epdf*> Eps ( Coms.length() );
-// 
-//              for ( int i = 0; i < Coms.length(); i++ ) {
-//                      Eps ( i ) = Coms ( i )->e();
-//              }
-// 
-//              iepdf->set_parameters ( w, Eps );
-//      }
-// 
-//      void condition ( const vec &cond ) {
-//              for ( int i = 0; i < Coms.length(); i++ ) {
-//                      Coms ( i )->condition ( cond );
-//              }
-//      };
-// };
+class mmix : public mpdf {
+protected:
+        //! Component (mpdfs)
+        Array<shared_ptr<mpdf> > Coms;
+        //!weights of the components
+        vec w;
+        //! dummy epdfs
+        epdf dummy_epdf;
+public:
+        //!Default constructor
+        mmix() : Coms(0), dummy_epdf() { set_ep(dummy_epdf);    }
+
+        //! Set weights \c w and components \c R
+        void set_parameters ( const vec &w0, const Array<shared_ptr<mpdf> > &Coms0 ) {
+                //!\TODO check if all components are OK
+                Coms = Coms0;
+                w=w0;   
+                
+                set_rv(Coms(0)->_rv());
+                set_rvc(Coms(0)->_rvc());
+        }
+        double evallogcond (const vec &dt, const vec &cond) {
+                double ll=0.0;
+                for (int i=0;i<Coms.length();i++){
+                        ll+=Coms(i)->evallogcond(dt,cond);
+                }
+                return ll;
+        }
+
+        vec samplecond (const vec &cond);
+
+};
 
 }

library/bdm/stat/exp_family.cpp

@@ -305 +305 @@
 }
 
-void eEmp::set_statistics ( const vec &w0, const epdf* epdf0 ) {
+void eEmp::set_statistics ( const vec &w0, const epdf &epdf0 ) {
         //it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
-        dim = epdf0->dimension();
+        dim = epdf0.dimension();
         w = w0;
         w /= sum ( w0 );//renormalize
         n = w.length();
         samples.set_size ( n );
-        dim = epdf0->dimension();
 
         for ( int i = 0; i < n; i++ ) {
-                samples ( i ) = epdf0->sample();
+                samples ( i ) = epdf0.sample();
         }
 }

library/bdm/stat/exp_family.h

@@ -24 +24 @@
 
 //! Global Uniform_RNG
-        extern Uniform_RNG UniRNG;
+extern Uniform_RNG UniRNG;
 //! Global Normal_RNG
-        extern Normal_RNG NorRNG;
+extern Normal_RNG NorRNG;
 //! Global Gamma_RNG
-        extern Gamma_RNG GamRNG;
-
-        /*!
-        * \brief General conjugate exponential family posterior density.
-
-        * More?...
-        */
-
-        class eEF : public epdf
-        {
-                public:
+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;
-                        //!TODO decide if it is really needed
-                        virtual void dupdate ( mat &v ) {it_error ( "Not implemented" );};
-                        //!Evaluate normalized log-probability
-                        virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
-                        //!Evaluate normalized log-probability
-                        virtual double evallog ( const vec &val ) const {
-                                double tmp;
-                                tmp= evallog_nn ( val )-lognc();
-//                              it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); 
-                                return tmp;}
-                        //!Evaluate normalized log-probability for many samples
-                        virtual vec evallog_m ( 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();
+                //! default constructor
+                eEF () : epdf () {};
+                //! logarithm of the normalizing constant, \f$\mathcal{I}\f$
+                virtual double lognc() const = 0;
+                //!TODO decide if it is really needed
+                virtual void dupdate (mat &v) {it_error ("Not implemented");};
+                //!Evaluate normalized log-probability
+                virtual double evallog_nn (const vec &val) const{it_error ("Not implemented");return 0.0;};
+                //!Evaluate normalized log-probability
+                virtual double evallog (const vec &val) const {
+                        double tmp;
+                        tmp = evallog_nn (val) - lognc();
+//                              it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+                        return tmp;
+                }
+                //!Evaluate normalized log-probability for many samples
+                virtual vec evallog_m (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_m (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) {it_error ("Not implemented");};
+};
+
+
+//! 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) {it_error ("Not implemented");};
+                //! Weighted update of sufficient statistics (Bayes rule)
+                virtual void bayes (const vec &data, const double w) {};
+                //original Bayes
+                void bayes (const vec &dt);
+                //!Flatten the posterior according to the given BMEF (of the same type!)
+                virtual void flatten (const BMEF * B) {it_error ("Not implemented");}
+                //!Flatten the posterior as if to keep nu0 data
+//      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, 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);
+                void from_setting (const Setting &root);
+                void validate() {
+                        it_assert (mu.length() == R.rows(), "parameters mismatch");
+                        dim = mu.length();
+                }
+                //!@}
+
+                //! \name Mathematical operations
+                //!@{
+
+                //! dupdate in exponential form (not really handy)
+                void dupdate (mat &v, double nu = 1.0);
+
+                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());}
+//      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;
+                //!@}
+
+                //! \name Access to attributes
+                //!@{
+
+                vec& _mu() {return mu;}
+                void set_mu (const vec mu0) { mu = mu0;}
+                sq_T& _R() {return R;}
+                const sq_T& _R() const {return R;}
+                //!@}
+
+};
+UIREGISTER (enorm<chmat>);
+UIREGISTER (enorm<ldmat>);
+UIREGISTER (enorm<fsqmat>);
+
+
+/*!
+* \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
+{
+        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);};
+
+                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;
                         }
-                        //!Power of the density, used e.g. to flatten the density
-                        virtual void pow ( double p ) {it_error ( "Not implemented" );};
-        };
-
-
-//! 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 ) {it_error ( "Not implemented" );};
-                        //! Weighted update of sufficient statistics (Bayes rule)
-                        virtual void bayes ( const vec &data, const double w ) {};
-                        //original Bayes
-                        void bayes ( const vec &dt );
-                        //!Flatten the posterior according to the given BMEF (of the same type!)
-                        virtual void flatten ( const BMEF * B ) {it_error ( "Not implemented" );}
-                        //!Flatten the posterior as if to keep nu0 data
-//      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, 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 );
-                        void from_setting(const Setting &root);
-                        void validate() {
-                                it_assert(mu.length()==R.rows(),"parameters mismatch");
-                                dim = mu.length();
+                }
+                //!@}
+
+                vec sample() const;
+                vec mean() const;
+                vec variance() const;
+
+                //! LS estimate of \f$\theta\f$
+                vec est_theta() const;
+
+                //! Covariance of the LS estimate
+                ldmat est_theta_cov() const;
+
+                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;};
+
+                //! \name Access attributes
+                //!@{
+
+                ldmat& _V() {return V;}
+                const ldmat& _V() const {return V;}
+                double& _nu()  {return nu;}
+                const double& _nu() const {return nu;}
+                void from_setting (const Setting &set) {
+                        UI::get (nu, set, "nu", UI::compulsory);
+                        UI::get (dimx, set, "dimx", UI::compulsory);
+                        mat V;
+                        UI::get (V, set, "V", UI::compulsory);
+                        set_parameters (dimx, V, nu);
+                        RV* rv = UI::build<RV> (set, "rv", UI::compulsory);
+                        set_rv (*rv);
+                        delete rv;
+                }
+                //!@}
+};
+UIREGISTER (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);};
+                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));}
+                //! In this instance, val is ...
+                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;
+                };
+                //!access function
+                vec& _beta()  {return beta;}
+                //!Set internal parameters
+};
+
+//! \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 &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();}
+                }
+};
+
+/*!
+ \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);};
+                void set_parameters (const vec &a, const vec &b) {alpha = a, beta = b;dim = alpha.length();};
+                //!@}
+
+                vec sample() const;
+                //! TODO: is it used anywhere?
+//      mat sample ( int N ) const;
+                double evallog (const vec &val) const;
+                double lognc () const;
+                //! Returns poiter to alpha and beta. Potentially dengerous: use with care!
+                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)); }
+
+                //! Load from structure with elements:
+                //!  \code
+                //! { alpha = [...];         // vector of alpha
+                //!   beta = [...];          // vector of beta
+                //!   rv = {class="RV",...}  // description
+                //! }
+                //! \endcode
+                //!@}
+                void from_setting (const Setting &set) {
+                        epdf::from_setting (set); // reads rv
+                        UI::get (alpha, set, "alpha", UI::compulsory);
+                        UI::get (beta, set, "beta", UI::compulsory);
+                        validate();
+                }
+                void validate() {
+                        it_assert (alpha.length() == beta.length(), "parameters do not match");
+                        dim = alpha.length();
+                }
+};
+UIREGISTER (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;};
+        vec sample() {it_error ( "Not implemented" );return 0;}
+};
+*/
+
+//!  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;
+                        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  {
+                        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;}
+                //! Load from structure with elements:
+                //!  \code
+                //! { high = [...];          // vector of upper bounds
+                //!   low = [...];           // vector of lower bounds
+                //!   rv = {class="RV",...}  // description of RV
+                //! }
+                //! \endcode
+                //!@}
+                void from_setting (const Setting &set) {
+                        epdf::from_setting (set); // reads rv and rvc
+
+                        UI::get (high, set, "high", UI::compulsory);
+                        UI::get (low, set, "low", UI::compulsory);
+                }
+};
+
+
+/*!
+ \brief Normal distributed linear function with linear function of mean value;
+
+ Mean value \f$mu=A*rvc+mu_0\f$.
+*/
+template < class sq_T, template <typename> class TEpdf = enorm >
+class mlnorm : public mpdf_internal< TEpdf<sq_T> >
+{
+        protected:
+                //! Internal epdf that arise by conditioning on \c rvc
+                mat A;
+                vec mu_const;
+//                      vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
+        public:
+                //! \name Constructors
+                //!@{
+                mlnorm() : mpdf_internal< TEpdf<sq_T> >() {};
+                mlnorm (const mat &A, const vec &mu0, const sq_T &R) : mpdf_internal< TEpdf<sq_T> >() {
+                        set_parameters (A, mu0, R);
+                }
+
+                //! Set \c A and \c R
+                void 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(), "");
+
+                        this->iepdf.set_parameters (zeros (A0.rows()), R0);
+                        A = A0;
+                        mu_const = mu0;
+                        this->dimc = A0.cols();
+                }
+                //!@}
+                //! 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
+                vec& _mu_const() {return mu_const;}
+                //!access function
+                mat& _A() {return A;}
+                //!access function
+                mat _R() { return this->iepdf._R().to_mat(); }
+
+                template<typename sq_M>
+                friend std::ostream &operator<< (std::ostream &os,  mlnorm<sq_M, enorm> &ml);
+
+                void from_setting (const Setting &set) {
+                        mpdf::from_setting (set);
+
+                        UI::get (A, set, "A", UI::compulsory);
+                        UI::get (mu_const, set, "const", UI::compulsory);
+                        mat R0;
+                        UI::get (R0, set, "R", UI::compulsory);
+                        set_parameters (A, mu_const, R0);
+                };
+};
+UIREGISTER (mlnorm<ldmat>);
+UIREGISTER (mlnorm<fsqmat>);
+UIREGISTER (mlnorm<chmat>);
+
+//! Mpdf with general function for mean value
+template<class sq_T>
+class mgnorm : public mpdf_internal< enorm< sq_T > >
+{
+        protected:
+//                      vec &mu; WHY NOT?
+                fnc* g;
+        public:
+                //!default constructor
+                mgnorm() : mpdf_internal<enorm<sq_T> >() { }
+                //!set mean function
+                inline void set_parameters (fnc* g0, const sq_T &R0);
+                inline void condition (const vec &cond);
+
+
+                /*! UI for mgnorm
+
+                The mgnorm is constructed from a structure with fields:
+                \code
+                system = {
+                        type = "mgnorm";
+                        // function for mean value evolution
+                        g = {type="fnc"; ... }
+
+                        // variance
+                        R = [1, 0,
+                                 0, 1];
+                        // --OR --
+                        dR = [1, 1];
+
+                        // == OPTIONAL ==
+
+                        // description of y variables
+                        y = {type="rv"; names=["y", "u"];};
+                        // description of u variable
+                        u = {type="rv"; names=[];}
+                };
+                \endcode
+
+                Result if
+                */
+
+                void from_setting (const Setting &set) {
+                        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);
+                }
+};
+
+UIREGISTER (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:
+                ldmat Lambda;
+                ldmat &_R;
+                ldmat Re;
+        public:
+                mlstudent () : mlnorm<ldmat, enorm> (),
+                                Lambda (),      _R (iepdf._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(), "");
+
+                        iepdf.set_parameters (mu0, Lambda);  //
+                        A = A0;
+                        mu_const = mu0;
+                        Re = R0;
+                        Lambda = Lambda0;
+                }
+                void condition (const vec &cond) {
+                        iepdf._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);
                         }
-                        //!@}
-
-                        //! \name Mathematical operations
-                        //!@{
-
-                        //! dupdate in exponential form (not really handy)
-                        void dupdate ( mat &v,double nu=1.0 );
-
-                        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() );}
-//      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;
-                        //!@}
-
-                        //! \name Access to attributes
-                        //!@{
-
-                        vec& _mu() {return mu;}
-                        void set_mu ( const vec mu0 ) { mu=mu0;}
-                        sq_T& _R() {return R;}
-                        const sq_T& _R() const {return R;}
-                        //!@}
-
-        };
-        UIREGISTER(enorm<chmat>);
-        UIREGISTER(enorm<ldmat>);
-        UIREGISTER(enorm<fsqmat>);
-
-
-        /*!
-        * \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
-        {
-                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 );};
-
-                        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;
+                        _R = Re;
+                        _R *= (1 + zeta);// / ( nu ); << nu is in Re!!!!!!
+                };
+
+};
+/*!
+\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 mpdf_internal<egamma>
+{
+        protected:
+
+                //! Constant \f$k\f$
+                double k;
+
+                //! cache of iepdf.beta
+                vec &_beta;
+
+        public:
+                //! Constructor
+                mgamma() : mpdf_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;};
+                //! Load from structure with elements:
+                //!  \code
+                //! { alpha = [...];         // vector of alpha
+                //!   k = 1.1;               // multiplicative constant k
+                //!   rv = {class="RV",...}  // description of RV
+                //!   rvc = {class="RV",...} // description of RV in condition
+                //! }
+                //! \endcode
+                //!@}
+                void from_setting (const Setting &set) {
+                        mpdf::from_setting (set); // reads rv and rvc
+                        vec betatmp; // ugly but necessary
+                        UI::get (betatmp, set, "beta", UI::compulsory);
+                        UI::get (k, set, "k", UI::compulsory);
+                        set_parameters (k, betatmp);
+                }
+};
+UIREGISTER (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 mpdf_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() : mpdf_internal<eigamma>(),
+                                k (0),
+                                _alpha (iepdf._alpha()),
+                                _beta (iepdf._beta()) {
+                }
+
+                migamma (const migamma &m) : mpdf_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*/);
+                        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;
+                        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;
+                        dimc = dimension();
+                };
+
+                void condition (const vec &val) {
+                        vec mean = elem_mult (refl, pow (val, l));
+                        migamma::condition (mean);
+                };
+
+                /*! UI for migamma_ref
+
+                The migamma_ref is constructed from a structure with fields:
+                \code
+                system = {
+                        type = "migamma_ref";
+                        ref = [1e-5; 1e-5; 1e-2 1e-3];            // reference vector
+                        l = 0.999;                                // constant l
+                        k = 0.1;                                  // constant k
+                        
+                        // == OPTIONAL ==
+                        // description of y variables
+                        y = {type="rv"; names=["y", "u"];};
+                        // description of u variable
+                        u = {type="rv"; names=[];}
+                };
+                \endcode
+
+                Result if
+                 */
+                void from_setting (const Setting &set);
+
+                // TODO dodelat void to_setting( Setting &set ) const;
+};
+
+
+UIREGISTER (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]
+
+*/
+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...
+
+==== 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 mlognorm : public mpdf_internal<elognorm>
+{
+        protected:
+                //! parameter 1/2*sigma^2
+                double sig2;
+
+                //! access
+                vec &mu;
+        public:
+                //! Constructor
+                mlognorm() : mpdf_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));
+
+                        dimc = size;
+                };
+
+                void condition (const vec &val) {
+                        mu = log (val) - sig2;//elem_mult ( refl,pow ( val,l ) );
+                };
+
+                /*! UI for mlognorm
+
+                The mlognorm is constructed from a structure with fields:
+                \code
+                system = {
+                        type = "mlognorm";
+                        k = 0.1;                                  // constant k
+                        mu0 = [1., 1.];
+                        
+                        // == OPTIONAL ==
+                        // description of y variables
+                        y = {type="rv"; names=["y", "u"];};
+                        // description of u variable
+                        u = {type="rv"; names=[];}
+                };
+                \endcode
+
+                 */
+                void from_setting (const Setting &set);
+
+                // TODO dodelat void to_setting( Setting &set ) const;
+
+};
+
+UIREGISTER (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:
+                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();
                                 }
                         }
-                        //!@}
-
-                        vec sample() const;
-                        vec mean() const;
-                        vec variance() const;
-
-                        //! LS estimate of \f$\theta\f$
-                        vec est_theta() const;
-
-                        //! Covariance of the LS estimate
-                        ldmat est_theta_cov() const;
-
-                        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;};
-
-                        //! \name Access attributes
-                        //!@{
-
-                        ldmat& _V() {return V;}
-                        const ldmat& _V() const {return V;}
-                        double& _nu()  {return nu;}
-                        const double& _nu() const {return nu;}
-                        void from_setting(const Setting &set){
-                                UI::get(nu, set, "nu", UI::compulsory);
-                                UI::get(dimx, set, "dimx", UI::compulsory);
-                                mat V;
-                                UI::get(V,set,"V", UI::compulsory);
-                                set_parameters(dimx, V, nu);
-                                RV* rv=UI::build<RV>(set,"rv", UI::compulsory);
-                                set_rv(*rv);
-                                delete rv;
+                        return X*Y._Ch();// return upper triangular part of the decomposition
+                }
+                vec sample () const {
+                        return vec (sample_mat()._data(), p*p);
+                }
+                //! 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;}
+};
+
+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;
+                //!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;
+                double l;
+                int p;
+
+        public:
+                rwiWishartCh() : eiW(),
+                                sqd (0), l (0), p (0) {
+                        set_ep (eiW);
+                }
+
+                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);
+                        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++;
                         }
-                        //!@}
-        };
-        UIREGISTER(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 );};
-                        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 ) );}
-                        //! In this instance, val is ...
-                        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;
-                        };
-                        //!access function
-                        vec& _beta()  {return beta;}
-                        //!Set internal parameters
-        };
-
-//! \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 &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();}
-                        }
-        };
-
-        /*!
-         \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 );};
-                        void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;dim = alpha.length();};
-                        //!@}
-
-                        vec sample() const;
-                        //! TODO: is it used anywhere?
-//      mat sample ( int N ) const;
-                        double evallog ( const vec &val ) const;
-                        double lognc () const;
-                        //! Returns poiter to alpha and beta. Potentially dengerous: use with care!
-                        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 ) ); }
-                        
-                        //! Load from structure with elements:
-                        //!  \code
-                        //! { alpha = [...];         // vector of alpha
-                        //!   beta = [...];          // vector of beta
-                        //!   rv = {class="RV",...}  // description
-                        //! }
-                        //! \endcode
-                        //!@}
-                        void from_setting(const Setting &set){
-                                epdf::from_setting(set); // reads rv
-                                UI::get(alpha,set,"alpha", UI::compulsory);
-                                UI::get(beta,set,"beta", UI::compulsory);
-                                validate();
-                        }
-                        void validate(){
-                                it_assert(alpha.length() ==beta.length(), "parameters do not match");
-                                dim =alpha.length();
-                        }
-        };
-UIREGISTER(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:
+
+                        eiW._setY (sqd*z);
+                }
+};
+
+//! Switch between various resampling methods.
+enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
+/*!
+\brief Weighted empirical density
+
+Used e.g. in particle filters.
+*/
+class eEmp: public epdf
+{
+        protected :
+                //! Number of particles
                 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;};
-                vec sample() {it_error ( "Not implemented" );return 0;}
-        };
-        */
-
-//!  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;
-                                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  {
-                                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;}
-                        //! Load from structure with elements:
-                        //!  \code
-                        //! { high = [...];          // vector of upper bounds
-                        //!   low = [...];           // vector of lower bounds
-                        //!   rv = {class="RV",...}  // description of RV
-                        //! }
-                        //! \endcode
-                        //!@}
-                        void from_setting(const Setting &set){
-                                epdf::from_setting(set); // reads rv and rvc
-
-                                UI::get(high,set,"high", UI::compulsory);
-                                UI::get(low,set,"low", UI::compulsory);
-                        }
-        };
-
-
-        /*!
-         \brief Normal distributed linear function with linear function of mean value;
-
-         Mean value \f$mu=A*rvc+mu_0\f$.
-        */
-        template<class sq_T, template <typename> class TEpdf =enorm>
-        class mlnorm : public mpdf_internal< TEpdf<sq_T> >
-        {
-                protected:
-                        //! Internal epdf that arise by conditioning on \c rvc
-                        mat A;
-                        vec mu_const;
-//                      vec& _mu; //cached epdf.mu; !!!!!! WHY NOT?
-                public:
-                        //! \name Constructors
-                        //!@{
-                        mlnorm():mpdf_internal< TEpdf<sq_T> >() {};
-                        mlnorm (const mat &A, const vec &mu0, const sq_T &R ) : mpdf_internal< TEpdf<sq_T> >()
-                        {
-                                set_parameters ( A,mu0,R );
-                        }
-
-                        //! Set \c A and \c R
-                        void 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(),"" );
-
-                                this->iepdf.set_parameters(zeros(A0.rows()), R0);
-                                A = A0;
-                                mu_const = mu0;
-                                this->dimc = A0.cols();
-                        }
-                        //!@}
-                        //! 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
-                        vec& _mu_const() {return mu_const;}
-                        //!access function
-                        mat& _A() {return A;}
-                        //!access function
-                        mat _R(){ return this->iepdf._R().to_mat(); }
-                
-                        template<typename sq_M>
-                        friend std::ostream &operator<< ( std::ostream &os,  mlnorm<sq_M,enorm> &ml );
-                        
-                        void from_setting(const Setting &set){
-                                mpdf::from_setting(set);        
-
-                                UI::get(A,set,"A", UI::compulsory);
-                                UI::get(mu_const,set,"const", UI::compulsory);
-                                mat R0;
-                                UI::get(R0,set,"R", UI::compulsory);
-                                set_parameters(A,mu_const,R0);
-                        };
-        };
-        UIREGISTER(mlnorm<ldmat>);
-        UIREGISTER(mlnorm<fsqmat>);
-        UIREGISTER(mlnorm<chmat>);
-
-//! Mpdf with general function for mean value
-        template<class sq_T>
-        class mgnorm : public mpdf_internal< enorm< sq_T > >
-        {
-                protected:
-//                      vec &mu; WHY NOT?
-                        fnc* g;
-                public:
-                        //!default constructor
-                        mgnorm():mpdf_internal<enorm<sq_T> >() { }
-                        //!set mean function
-                        inline void set_parameters ( fnc* g0, const sq_T &R0 );
-                        inline void condition ( const vec &cond );
-
-
-                        /*! UI for mgnorm
-
-                        The mgnorm is constructed from a structure with fields:
-                        \code
-                        system = {
-                                type = "mgnorm";
-                                // function for mean value evolution
-                                g = {type="fnc"; ... }
-
-                                // variance
-                                R = [1, 0,
-                                         0, 1];
-                                // --OR --
-                                dR = [1, 1];
-
-                                // == OPTIONAL ==
-
-                                // description of y variables
-                                y = {type="rv"; names=["y", "u"];};
-                                // description of u variable
-                                u = {type="rv"; names=[];}
-                        };
-                        \endcode
-
-                        Result if
-                        */
-
-                        void from_setting( const Setting &set ) 
-                        {       
-                                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);
-                        }
-        };
-
-        UIREGISTER(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:
-                        ldmat Lambda;
-                        ldmat &_R;
-                        ldmat Re;
-                public:
-                        mlstudent ( ) :mlnorm<ldmat,enorm> (),
-                                        Lambda (),      _R ( iepdf._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(),"" );
-
-                                iepdf.set_parameters ( mu0,Lambda ); //
-                                A = A0;
-                                mu_const = mu0;
-                                Re=R0;
-                                Lambda = Lambda0;
-                        }
-                        void condition ( const vec &cond )
-                        {
-                                iepdf._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!!!!!!
-                        };
-
-        };
-        /*!
-        \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 mpdf_internal<egamma>
-        {
-                protected:
-
-                        //! Constant \f$k\f$
-                        double k;
-
-                        //! cache of iepdf.beta
-                        vec &_beta;
-
-                public:
-                        //! Constructor
-                        mgamma():mpdf_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;};
-                        //! Load from structure with elements:
-                        //!  \code
-                        //! { alpha = [...];         // vector of alpha
-                        //!   k = 1.1;               // multiplicative constant k
-                        //!   rv = {class="RV",...}  // description of RV
-                        //!   rvc = {class="RV",...} // description of RV in condition
-                        //! }
-                        //! \endcode
-                        //!@}
-                        void from_setting(const Setting &set){
-                                mpdf::from_setting(set); // reads rv and rvc
-                                vec betatmp; // ugly but necessary
-                                UI::get(betatmp,set,"beta", UI::compulsory);
-                                UI::get(k,set,"k", UI::compulsory);
-                                set_parameters(k,betatmp);
-                        }
-        };
-        UIREGISTER(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 mpdf_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():mpdf_internal<eigamma>(),
-                            k(0),
-                            _alpha(iepdf._alpha()),
-                            _beta(iepdf._beta()) {
-                        }
-
-                        migamma(const migamma &m):mpdf_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*/ );
-                                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;
-                                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;
-                                dimc = dimension();
-                        };
-
-                        void condition ( const vec &val )
-                        {
-                                vec mean=elem_mult ( refl,pow ( val,l ) );
-                                migamma::condition ( mean );
-                        };
-
-                        /*! UI for migamma_ref
-
-                        The migamma_ref is constructed from a structure with fields:
-                        \code
-                        system = {
-                                type = "migamma_ref";
-                                ref = [1e-5; 1e-5; 1e-2 1e-3];            // reference vector
-                                l = 0.999;                                // constant l
-                                k = 0.1;                                  // constant k
-                                
-                                // == OPTIONAL ==
-                                // description of y variables
-                                y = {type="rv"; names=["y", "u"];};
-                                // description of u variable
-                                u = {type="rv"; names=[];}
-                        };
-                        \endcode
-
-                        Result if
-                         */
-                        void from_setting( const Setting &set );
-
-                        // TODO dodelat void to_setting( Setting &set ) const;
-        };
-
-
-        UIREGISTER(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]
-
-        */
-        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...
-
-        ==== 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 mlognorm : public mpdf_internal<elognorm>
-        {
-                protected:
-                        //! parameter 1/2*sigma^2
-                        double sig2;
-
-                        //! access
-                        vec &mu;
-                public:
-                        //! Constructor
-                        mlognorm():mpdf_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 ) );
-
-                                dimc = size;
-                        };
-
-                        void condition ( const vec &val )
-                        {
-                                mu=log ( val )-sig2;//elem_mult ( refl,pow ( val,l ) );
-                        };
-
-                        /*! UI for mlognorm
-
-                        The mlognorm is constructed from a structure with fields:
-                        \code
-                        system = {
-                                type = "mlognorm";
-                                k = 0.1;                                  // constant k
-                                mu0 = [1., 1.];
-                                
-                                // == OPTIONAL ==
-                                // description of y variables
-                                y = {type="rv"; names=["y", "u"];};
-                                // description of u variable
-                                u = {type="rv"; names=[];}
-                        };
-                        \endcode
-
-                         */
-                        void from_setting( const Setting &set );
-
-                        // TODO dodelat void to_setting( Setting &set ) const;
-
-        };
-        
-        UIREGISTER(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:
-                        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 );
-                        }
-                        //! 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;}
-        };
-
-        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;
-                        //!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;
-                        double l;
-                        int p;
-
-                public:
-                        rwiWishartCh():eiW(),
-                            sqd(0), l(0), p(0) {
-                            set_ep(eiW);
-                        }
-
-                        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 );
-                                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 );
-                        }
-        };
-
-//! Switch between various resampling methods.
-        enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
-        /*!
-        \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 );};
-                        //! 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 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.
-                        ivec resample ( RESAMPLING_METHOD method=SYSTEMATIC );
-                        //! inherited operation : NOT implemneted
-                        vec sample() const {it_error ( "Not implemented" );return 0;}
-                        //! inherited operation : NOT implemneted
-                        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 );
-                        }
-                        //! For this class, qbounds are minimum and maximum value of the population!
-                        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 );}
-                                        }
+                //! 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);};
+                //! 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 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.
+                ivec resample (RESAMPLING_METHOD method = SYSTEMATIC);
+                //! inherited operation : NOT implemneted
+                vec sample() const {it_error ("Not implemented");return 0;}
+                //! inherited operation : NOT implemneted
+                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);
+                }
+                //! For this class, qbounds are minimum and maximum value of the population!
+                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 )
-        {
+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
-                validate();
-        }
-
-        template<class sq_T>
-        void enorm<sq_T>::dupdate ( mat &v, double nu )
-        {
-                //
-        };
+        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
+        validate();
+}
+
+template<class sq_T>
+void enorm<sq_T>::dupdate (mat &v, double nu)
+{
+        //
+};
 
 // template<class sq_T>
@@ -1204 +1171 @@
 // };
 
-        template<class sq_T>
-        vec enorm<sq_T>::sample() const
-        {
-                vec x ( dim );
+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;
-        };
+        NorRNG.sample_vector (dim, x);
+        vec smp = R.sqrt_mult (x);
+
+        smp += mu;
+        return smp;
+};
 
 // template<class sq_T>
@@ -1224 +1191 @@
 // };
 
-        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>
+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;
+};
 
 
@@ -1267 +1234 @@
 
 
-        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>                            
-                        void mgnorm<sq_T >::set_parameters ( 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 );};
-
-        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.iepdf->_R().to_mat() <<endl;
-                return os;
-        };
+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>
+void mgnorm<sq_T >::set_parameters (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);};
+
+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;
+};
 
 }

library/bdm/stat/merger.h

@@ -178 +178 @@
         //! Set support points from a pdf by drawing N samples
         void set_support ( const epdf &overall, int N ) {
-                eSmp.set_statistics ( &overall, N );
+                eSmp.set_statistics ( overall, N );
                 Npoints = N;
         }