Changeset 477 for library/bdm/stat

Timestamp:

08/05/09 14:40:03 (15 years ago)

Author:

mido

Message:

panove, vite, jak jsem peclivej na upravu kodu.. snad se vam bude libit:) konfigurace je v souboru /system/astylerc

Location:

library/bdm/stat

Files:

• emix.cpp (modified) (7 diffs)
• emix.h (modified) (19 diffs)
• exp_family.cpp (modified) (15 diffs)
• merger.cpp (modified) (1 diff)
• merger.h (modified) (2 diffs)

library/bdm/stat/emix.cpp

@@ -1 +1 @@
 #include "emix.h"
 
-namespace bdm{
+namespace bdm {
 
 void emix::set_parameters ( const vec &w0, const Array<epdf*> &Coms0, bool copy ) {
-        w = w0/sum ( w0 );
-        dim = Coms0(0)->dimension();
-        bool isnamed = Coms0(0)->isnamed();
+        w = w0 / sum ( w0 );
+        dim = Coms0 ( 0 )->dimension();
+        bool isnamed = Coms0 ( 0 )->isnamed();
         int i;
         RV tmp_rv;
-        if (isnamed) tmp_rv=Coms0(0)->_rv();
-        
-        for ( i=0;i<w.length();i++ ) {
-                it_assert_debug ( dim== ( Coms0 ( i )->dimension() ),"Component sizes do not match!" ); 
-                it_assert_debug ( !isnamed || tmp_rv.equal( Coms0 ( i )->_rv() ),"Component RVs do not match!" );
+        if ( isnamed ) tmp_rv = Coms0 ( 0 )->_rv();
+
+        for ( i = 0; i < w.length(); i++ ) {
+                it_assert_debug ( dim == ( Coms0 ( i )->dimension() ), "Component sizes do not match!" );
+                it_assert_debug ( !isnamed || tmp_rv.equal ( Coms0 ( i )->_rv() ), "Component RVs do not match!" );
         }
         if ( copy ) {
-                Coms.set_length(Coms0.length());
-                for ( i=0;i<w.length();i++ ) {it_error("Not imp...");
-                        *Coms ( i ) =*Coms0 ( i );}
-                destroyComs=true;
-        }
-        else {
+                Coms.set_length ( Coms0.length() );
+                for ( i = 0; i < w.length(); i++ ) {
+                        it_error ( "Not imp..." );
+                        *Coms ( i ) = *Coms0 ( i );
+                }
+                destroyComs = true;
+        } else {
                 Coms = Coms0;
-                destroyComs=false;
-        }
-        if (isnamed) epdf::set_rv(tmp_rv); //coms aer already OK, no need for set_rv
+                destroyComs = false;
+        }
+        if ( isnamed ) epdf::set_rv ( tmp_rv ); //coms aer already OK, no need for set_rv
 }
 
@@ -32 +33 @@
         vec cumDist = cumsum ( w );
         double u0;
-        #pragma omp critical
+#pragma omp critical
         u0 = UniRNG.sample();
 
-        int i=0;
-        while ( ( cumDist ( i ) <u0 ) && ( i< ( w.length()-1 ) ) ) {i++;}
+        int i = 0;
+        while ( ( cumDist ( i ) < u0 ) && ( i < ( w.length() - 1 ) ) ) {
+                i++;
+        }
 
         return Coms ( i )->sample();
 }
 
-emix* emix::marginal(const RV &rv) const{
-        it_assert_debug(isnamed(), "rvs are not assigned");
-                        
-        Array<epdf*> Cn(Coms.length());
-        for(int i=0;i<Coms.length();i++){Cn(i)=Coms(i)->marginal(rv);}
+emix* emix::marginal ( const RV &rv ) const {
+        it_assert_debug ( isnamed(), "rvs are not assigned" );
+
+        Array<epdf*> Cn ( Coms.length() );
+        for ( int i = 0; i < Coms.length(); i++ ) {
+                Cn ( i ) = Coms ( i )->marginal ( rv );
+        }
         emix* tmp = new emix();
-        tmp->set_parameters(w,Cn,false);
+        tmp->set_parameters ( w, Cn, false );
         tmp->ownComs();
         return tmp;
 }
 
-mratio* emix::condition(const RV &rv) const{
-        it_assert_debug(isnamed(), "rvs are not assigned");
-        return new mratio(this,rv);
+mratio* emix::condition ( const RV &rv ) const {
+        it_assert_debug ( isnamed(), "rvs are not assigned" );
+        return new mratio ( this, rv );
 };
 
 void egiwmix::set_parameters ( const vec &w0, const Array<egiw*> &Coms0, bool copy ) {
-        w = w0/sum ( w0 );
-        dim = Coms0(0)->dimension();
+        w = w0 / sum ( w0 );
+        dim = Coms0 ( 0 )->dimension();
         int i;
-        for ( i=0;i<w.length();i++ ) {
-                it_assert_debug ( dim== ( Coms0 ( i )->dimension() ),"Component sizes do not match!" );
+        for ( i = 0; i < w.length(); i++ ) {
+                it_assert_debug ( dim == ( Coms0 ( i )->dimension() ), "Component sizes do not match!" );
         }
         if ( copy ) {
-                Coms.set_length(Coms0.length());
-                for ( i=0;i<w.length();i++ ) {it_error("Not imp...");
-                        *Coms ( i ) =*Coms0 ( i );}
-                destroyComs=true;
-        }
-        else {
+                Coms.set_length ( Coms0.length() );
+                for ( i = 0; i < w.length(); i++ ) {
+                        it_error ( "Not imp..." );
+                        *Coms ( i ) = *Coms0 ( i );
+                }
+                destroyComs = true;
+        } else {
                 Coms = Coms0;
-                destroyComs=false;
+                destroyComs = false;
         }
 }
@@ -80 +86 @@
         vec cumDist = cumsum ( w );
         double u0;
-        #pragma omp critical
+#pragma omp critical
         u0 = UniRNG.sample();
 
-        int i=0;
-        while ( ( cumDist ( i ) <u0 ) && ( i< ( w.length()-1 ) ) ) {i++;}
+        int i = 0;
+        while ( ( cumDist ( i ) < u0 ) && ( i < ( w.length() - 1 ) ) ) {
+                i++;
+        }
 
         return Coms ( i )->sample();
@@ -90 +98 @@
 
 vec egiwmix::mean() const {
-        int i; vec mu = zeros ( dim );
-        for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); }
+        int i;
+        vec mu = zeros ( dim );
+        for ( i = 0; i < w.length(); i++ ) {
+                mu += w ( i ) * Coms ( i )->mean();
+        }
         return mu;
 }
@@ -98 +109 @@
         // non-central moment
         vec mom2 = zeros ( dim );
-        for ( int i = 0;i < w.length();i++ ) {
+        for ( int i = 0; i < w.length(); i++ ) {
                 // pow is overloaded, we have to use another approach
-                mom2 += w ( i ) * (Coms(i)->variance() + elem_mult ( Coms(i)->mean(), Coms(i)->mean() )); 
+                mom2 += w ( i ) * ( Coms ( i )->variance() + elem_mult ( Coms ( i )->mean(), Coms ( i )->mean() ) );
         }
         // central moment
         // pow is overloaded, we have to use another approach
-        return mom2 - elem_mult( mean(), mean() );
-}
-
-emix* egiwmix::marginal(const RV &rv) const{
-        it_assert_debug(isnamed(), "rvs are not assigned");
-                        
-        Array<epdf*> Cn(Coms.length());
-        for(int i=0;i<Coms.length();i++){Cn(i)=Coms(i)->marginal(rv);}
+        return mom2 - elem_mult ( mean(), mean() );
+}
+
+emix* egiwmix::marginal ( const RV &rv ) const {
+        it_assert_debug ( isnamed(), "rvs are not assigned" );
+
+        Array<epdf*> Cn ( Coms.length() );
+        for ( int i = 0; i < Coms.length(); i++ ) {
+                Cn ( i ) = Coms ( i )->marginal ( rv );
+        }
         emix* tmp = new emix();
-        tmp->set_parameters(w,Cn,false);
+        tmp->set_parameters ( w, Cn, false );
         tmp->ownComs();
         return tmp;
@@ -124 +137 @@
 
         double sumVecCommon;                            // common part for many terms in eq.
-        int len = w.length();                           // no. of mix components        
-        int dimLS = Coms(1)->_V()._D().length() - 1;    // dim of LS
-        vec vecNu(len);                                 // vector of dfms of components
-        vec vecD(len);                                  // vector of LS reminders of comps.
-        vec vecCommon(len);                             // vector of common parts
+        int len = w.length();                           // no. of mix components
+        int dimLS = Coms ( 1 )->_V()._D().length() - 1;         // dim of LS
+        vec vecNu ( len );                                      // vector of dfms of components
+        vec vecD ( len );                                       // vector of LS reminders of comps.
+        vec vecCommon ( len );                          // vector of common parts
         mat matVecsTheta;                               // matrix which rows are theta vects.
 
         // fill in the vectors vecNu, vecD and matVecsTheta
-        for ( int i=0; i<len; i++ ) {
-                vecNu.shift_left( Coms(i)->_nu() );
-                vecD.shift_left( Coms(i)->_V()._D()(0) );
-                matVecsTheta.append_row( Coms(i)->est_theta()  );
+        for ( int i = 0; i < len; i++ ) {
+                vecNu.shift_left ( Coms ( i )->_nu() );
+                vecD.shift_left ( Coms ( i )->_V()._D() ( 0 ) );
+                matVecsTheta.append_row ( Coms ( i )->est_theta() );
         }
 
         // calculate the common parts and their sum
-        vecCommon = elem_mult ( w, elem_div(vecNu, vecD) );
-        sumVecCommon = sum(vecCommon);
+        vecCommon = elem_mult ( w, elem_div ( vecNu, vecD ) );
+        sumVecCommon = sum ( vecCommon );
 
         // LS estimator of theta
-        vec aprEstTheta(dimLS);  aprEstTheta.zeros();
-        for ( int i=0; i<len; i++ ) {
-                aprEstTheta +=  matVecsTheta.get_row( i ) * vecCommon ( i );
+        vec aprEstTheta ( dimLS );
+        aprEstTheta.zeros();
+        for ( int i = 0; i < len; i++ ) {
+                aprEstTheta +=  matVecsTheta.get_row ( i ) * vecCommon ( i );
         }
         aprEstTheta /= sumVecCommon;
-        
-        
+
+
         // LS estimator of dfm
         double aprNu;
-        double A = log( sumVecCommon );         // Term 'A' in equation
-
-        for ( int i=0; i<len; i++ ) {
-                A += w(i) * ( log( vecD(i) ) - psi( 0.5 * vecNu(i) ) );
-        }
-
-        aprNu = ( 1 + sqrt(1 + 2 * (A - LOG2)/3 ) ) / ( 2 * (A - LOG2) );
+        double A = log ( sumVecCommon );                // Term 'A' in equation
+
+        for ( int i = 0; i < len; i++ ) {
+                A += w ( i ) * ( log ( vecD ( i ) ) - psi ( 0.5 * vecNu ( i ) ) );
+        }
+
+        aprNu = ( 1 + sqrt ( 1 + 2 * ( A - LOG2 ) / 3 ) ) / ( 2 * ( A - LOG2 ) );
 
 
@@ -167 +181 @@
         // the following code is very numerically sensitive, thus
         // we have to eliminate decompositions etc. as much as possible
-        mat aprC = zeros(dimLS, dimLS);
-        for ( int i=0; i<len; i++ ) {
-                aprC += Coms(i)->est_theta_cov().to_mat() * w(i); 
-                vec tmp = ( matVecsTheta.get_row(i) - aprEstTheta );
-                aprC += vecCommon(i) * outer_product( tmp, tmp);
+        mat aprC = zeros ( dimLS, dimLS );
+        for ( int i = 0; i < len; i++ ) {
+                aprC += Coms ( i )->est_theta_cov().to_mat() * w ( i );
+                vec tmp = ( matVecsTheta.get_row ( i ) - aprEstTheta );
+                aprC += vecCommon ( i ) * outer_product ( tmp, tmp );
         }
 
         // Construct GiW pdf :: BEGIN
-        ldmat aprCinv ( inv(aprC) );
-        vec D = concat( aprD, aprCinv._D() );
-        mat L = eye(dimLS+1);
-        L.set_submatrix(1,0, aprCinv._L() * aprEstTheta);
-        L.set_submatrix(1,1, aprCinv._L());
-        ldmat aprLD (L, D);
-
-        egiw* aprgiw = new egiw(1, aprLD, aprNu);
+        ldmat aprCinv ( inv ( aprC ) );
+        vec D = concat ( aprD, aprCinv._D() );
+        mat L = eye ( dimLS + 1 );
+        L.set_submatrix ( 1, 0, aprCinv._L() * aprEstTheta );
+        L.set_submatrix ( 1, 1, aprCinv._L() );
+        ldmat aprLD ( L, D );
+
+        egiw* aprgiw = new egiw ( 1, aprLD, aprNu );
         return aprgiw;
 };

library/bdm/stat/emix.h

@@ -14 +14 @@
 #define EMIX_H
 
-#define LOG2  0.69314718055995  
+#define LOG2  0.69314718055995
 
 #include "../shared_ptr.h"
@@ -48 +48 @@
         //!Default constructor. By default, the given epdf is not copied!
         //! It is assumed that this function will be used only temporarily.
-        mratio ( const epdf* nom0, const RV &rv, bool copy=false ) :mpdf ( ), dl ( ) {
+        mratio ( const epdf* nom0, const RV &rv, bool copy = false ) : mpdf ( ), dl ( ) {
                 // adjust rv and rvc
                 rvc = nom0->_rv().subt ( rv );
                 dimc = rvc._dsize();
-                set_ep(shared_ptr<epdf>(new epdf));
-                e()->set_parameters(rv._dsize());
-                e()->set_rv(rv);
-                
+                set_ep ( shared_ptr<epdf> ( new epdf ) );
+                e()->set_parameters ( rv._dsize() );
+                e()->set_rv ( rv );
+
                 //prepare data structures
-                if ( copy ) {it_error ( "todo" ); destroynom=true; }
-                else { nom = nom0; destroynom = false; }
-                it_assert_debug ( rvc.length() >0,"Makes no sense to use this object!" );               
-                
+                if ( copy ) {
+                        it_error ( "todo" );
+                        destroynom = true;
+                } else {
+                        nom = nom0;
+                        destroynom = false;
+                }
+                it_assert_debug ( rvc.length() > 0, "Makes no sense to use this object!" );
+
                 // build denominator
                 den = nom->marginal ( rvc );
-                dl.set_connection(rv,rvc,nom0->_rv());
+                dl.set_connection ( rv, rvc, nom0->_rv() );
         };
         double evallogcond ( const vec &val, const vec &cond ) {
                 double tmp;
                 vec nom_val ( e()->dimension() + dimc );
-                dl.pushup_cond ( nom_val,val,cond );
+                dl.pushup_cond ( nom_val, val, cond );
                 tmp = exp ( nom->evallog ( nom_val ) - den->evallog ( cond ) );
-                it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+                it_assert_debug ( std::isfinite ( tmp ), "Infinite value" );
                 return tmp;
         }
         //! Object takes ownership of nom and will destroy it
-        void ownnom() {destroynom=true;}
+        void ownnom() {
+                destroynom = true;
+        }
         //! Default destructor
-        ~mratio() {delete den; if ( destroynom ) {delete nom;}}
+        ~mratio() {
+                delete den;
+                if ( destroynom ) {
+                        delete nom;
+                }
+        }
 };
 
@@ -102 +114 @@
         //! Set weights \c w and components \c Coms
         //!By default Coms are copied inside. Parameter \c copy can be set to false if Coms live externally. Use method ownComs() if Coms should be destroyed by the destructor.
-        void set_parameters ( const vec &w, const Array<epdf*> &Coms, bool copy=false );
+        void set_parameters ( const vec &w, const Array<epdf*> &Coms, bool copy = false );
 
         vec sample() const;
         vec mean() const {
-                int i; vec mu = zeros ( dim );
-                for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); }
+                int i;
+                vec mu = zeros ( dim );
+                for ( i = 0; i < w.length(); i++ ) {
+                        mu += w ( i ) * Coms ( i )->mean();
+                }
                 return mu;
         }
@@ -113 +128 @@
                 //non-central moment
                 vec mom2 = zeros ( dim );
-                for ( int i = 0;i < w.length();i++ ) {mom2 += w ( i ) * (Coms(i)->variance() + pow ( Coms ( i )->mean(),2 )); }
+                for ( int i = 0; i < w.length(); i++ ) {
+                        mom2 += w ( i ) * ( Coms ( i )->variance() + pow ( Coms ( i )->mean(), 2 ) );
+                }
                 //central moment
-                return mom2-pow ( mean(),2 );
+                return mom2 - pow ( mean(), 2 );
         }
         double evallog ( const vec &val ) const {
                 int i;
                 double sum = 0.0;
-                for ( i = 0;i < w.length();i++ ) {sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );}
-                if ( sum==0.0 ) {sum=std::numeric_limits<double>::epsilon();}
-                double tmp=log ( sum );
-                it_assert_debug ( std::isfinite ( tmp ),"Infinite" );
+                for ( i = 0; i < w.length(); i++ ) {
+                        sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );
+                }
+                if ( sum == 0.0 ) {
+                        sum = std::numeric_limits<double>::epsilon();
+                }
+                double tmp = log ( sum );
+                it_assert_debug ( std::isfinite ( tmp ), "Infinite" );
                 return tmp;
         };
         vec evallog_m ( const mat &Val ) const {
-                vec x=zeros ( Val.cols() );
+                vec x = zeros ( Val.cols() );
                 for ( int i = 0; i < w.length(); i++ ) {
-                        x+= w ( i ) *exp ( Coms ( i )->evallog_m ( Val ) );
+                        x += w ( i ) * exp ( Coms ( i )->evallog_m ( Val ) );
                 }
                 return log ( x );
@@ -147 +168 @@
 //Access methods
         //! returns a pointer to the internal mean value. Use with Care!
-        vec& _w() {return w;}
-        virtual ~emix() {if ( destroyComs ) {for ( int i=0;i<Coms.length();i++ ) {delete Coms ( i );}}}
+        vec& _w() {
+                return w;
+        }
+        virtual ~emix() {
+                if ( destroyComs ) {
+                        for ( int i = 0; i < Coms.length(); i++ ) {
+                                delete Coms ( i );
+                        }
+                }
+        }
         //! Auxiliary function for taking ownership of the Coms()
-        void ownComs() {destroyComs=true;}
+        void ownComs() {
+                destroyComs = true;
+        }
 
         //!access function
-        epdf* _Coms ( int i ) {return Coms ( i );}
-        void set_rv(const RV &rv){
-                epdf::set_rv(rv);
-                for(int i=0;i<Coms.length();i++){Coms(i)->set_rv(rv);}
+        epdf* _Coms ( int i ) {
+                return Coms ( i );
+        }
+        void set_rv ( const RV &rv ) {
+                epdf::set_rv ( rv );
+                for ( int i = 0; i < Coms.length(); i++ ) {
+                        Coms ( i )->set_rv ( rv );
+                }
         }
 };
@@ -179 +214 @@
         //! Set weights \c w and components \c Coms
         //!By default Coms are copied inside. Parameter \c copy can be set to false if Coms live externally. Use method ownComs() if Coms should be destroyed by the destructor.
-        void set_parameters ( const vec &w, const Array<egiw*> &Coms, bool copy=false );
+        void set_parameters ( const vec &w, const Array<egiw*> &Coms, bool copy = false );
 
         //!return expected value
@@ -188 +223 @@
 
         //!return the expected variance
-        vec variance() const;
+        vec variance() const;
 
         // TODO!!! Defined to follow ANSI and/or for future development
         void mean_mat ( mat &M, mat&R ) const {};
-        double evallog_nn ( const vec &val ) const {return 0;};
-        double lognc () const {return 0;};
+        double evallog_nn ( const vec &val ) const {
+                return 0;
+        };
+        double lognc () const {
+                return 0;
+        };
         emix* marginal ( const RV &rv ) const;
 
 //Access methods
         //! returns a pointer to the internal mean value. Use with Care!
-        vec& _w() {return w;}
-        virtual ~egiwmix() {if ( destroyComs ) {for ( int i=0;i<Coms.length();i++ ) {delete Coms ( i );}}}
+        vec& _w() {
+                return w;
+        }
+        virtual ~egiwmix() {
+                if ( destroyComs ) {
+                        for ( int i = 0; i < Coms.length(); i++ ) {
+                                delete Coms ( i );
+                        }
+                }
+        }
         //! Auxiliary function for taking ownership of the Coms()
-        void ownComs() {destroyComs=true;}
+        void ownComs() {
+                destroyComs = true;
+        }
 
         //!access function
-        egiw* _Coms ( int i ) {return Coms ( i );}
-
-        void set_rv(const RV &rv){
-                egiw::set_rv(rv);
-                for(int i=0;i<Coms.length();i++){Coms(i)->set_rv(rv);}
+        egiw* _Coms ( int i ) {
+                return Coms ( i );
+        }
+
+        void set_rv ( const RV &rv ) {
+                egiw::set_rv ( rv );
+                for ( int i = 0; i < Coms.length(); i++ ) {
+                        Coms ( i )->set_rv ( rv );
+                }
         }
 
@@ -236 +289 @@
         /*!\brief Constructor from list of mFacs,
         */
-        mprod():dummy(new epdf()) { }
-        mprod (Array<mpdf*> mFacs):
-            dummy(new epdf()) {
-            set_elements(mFacs);
-        }
-
-        void set_elements(Array<mpdf*> mFacs , bool own=false) {
-                
-                compositepdf::set_elements(mFacs,own);
-                dls.set_size(mFacs.length());
-                epdfs.set_size(mFacs.length());
-                                
-                set_ep(dummy);
-                RV rv=getrv ( true );
+        mprod() : dummy ( new epdf() ) { }
+        mprod ( Array<mpdf*> mFacs ) :
+                        dummy ( new epdf() ) {
+                set_elements ( mFacs );
+        }
+
+        void set_elements ( Array<mpdf*> mFacs , bool own = false ) {
+
+                compositepdf::set_elements ( mFacs, own );
+                dls.set_size ( mFacs.length() );
+                epdfs.set_size ( mFacs.length() );
+
+                set_ep ( dummy );
+                RV rv = getrv ( true );
                 set_rv ( rv );
-                dummy->set_parameters(rv._dsize());
-                setrvc(e()->_rv(), rvc);
+                dummy->set_parameters ( rv._dsize() );
+                setrvc ( e()->_rv(), rvc );
                 // rv and rvc established = > we can link them with mpdfs
-                for ( int i = 0;i < mpdfs.length(); i++ ) {
+                for ( int i = 0; i < mpdfs.length(); i++ ) {
                         dls ( i ) = new datalink_m2m;
-                        dls(i)->set_connection( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), _rv(), _rvc() );
-                }
-
-                for ( int i=0; i<mpdfs.length(); i++ ) {
-                        epdfs(i) = mpdfs(i)->e();
+                        dls ( i )->set_connection ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), _rv(), _rvc() );
+                }
+
+                for ( int i = 0; i < mpdfs.length(); i++ ) {
+                        epdfs ( i ) = mpdfs ( i )->e();
                 }
         };
@@ -267 +320 @@
                 int i;
                 double res = 0.0;
-                for ( i = mpdfs.length() - 1;i >= 0;i-- ) {
+                for ( i = mpdfs.length() - 1; i >= 0; i-- ) {
                         /*                      if ( mpdfs(i)->_rvc().count() >0) {
                                                         mpdfs ( i )->condition ( dls ( i )->get_cond ( val,cond ) );
@@ -280 +333 @@
                 return res;
         }
-        vec evallogcond_m(const mat &Dt, const vec &cond) {
-                vec tmp(Dt.cols());
-                for(int i=0;i<Dt.cols(); i++){
-                        tmp(i) = evallogcond(Dt.get_col(i),cond);
-                }
-                return tmp;
-        };
-        vec evallogcond_m(const Array<vec> &Dt, const vec &cond) {
-                vec tmp(Dt.length());
-                for(int i=0;i<Dt.length(); i++){
-                        tmp(i) = evallogcond(Dt(i),cond);
-                }
-                return tmp;             
+        vec evallogcond_m ( const mat &Dt, const vec &cond ) {
+                vec tmp ( Dt.cols() );
+                for ( int i = 0; i < Dt.cols(); i++ ) {
+                        tmp ( i ) = evallogcond ( Dt.get_col ( i ), cond );
+                }
+                return tmp;
+        };
+        vec evallogcond_m ( const Array<vec> &Dt, const vec &cond ) {
+                vec tmp ( Dt.length() );
+                for ( int i = 0; i < Dt.length(); i++ ) {
+                        tmp ( i ) = evallogcond ( Dt ( i ), cond );
+                }
+                return tmp;
         };
 
@@ -298 +351 @@
         //TODO smarter...
         vec samplecond ( const vec &cond ) {
-            //! Ugly hack to help to discover if mpfs are not in proper order. Correct solution = check that explicitely.
-            vec smp= std::numeric_limits<double>::infinity() * ones ( e()->dimension() );
+                //! Ugly hack to help to discover if mpfs are not in proper order. Correct solution = check that explicitely.
+                vec smp = std::numeric_limits<double>::infinity() * ones ( e()->dimension() );
                 vec smpi;
                 // Hard assumption here!!! We are going backwards, to assure that samples that are needed from smp are already generated!
-                for ( int i = ( mpdfs.length() - 1 );i >= 0;i-- ) {
+                for ( int i = ( mpdfs.length() - 1 ); i >= 0; i-- ) {
                         if ( mpdfs ( i )->dimensionc() ) {
-                                mpdfs ( i )->condition ( dls ( i )->get_cond ( smp ,cond ) ); // smp is val here!!
+                                mpdfs ( i )->condition ( dls ( i )->get_cond ( smp , cond ) ); // smp is val here!!
                         }
                         smpi = epdfs ( i )->sample();
@@ -314 +367 @@
         }
         mat samplecond ( const vec &cond,  int N ) {
-                mat Smp ( dimension(),N );
-                for ( int i=0;i<N;i++ ) {Smp.set_col ( i,samplecond ( cond ) );}
+                mat Smp ( dimension(), N );
+                for ( int i = 0; i < N; i++ ) {
+                        Smp.set_col ( i, samplecond ( cond ) );
+                }
                 return Smp;
         }
@@ -327 +382 @@
         //! \endcode
         //!@}
-        void from_setting(const Setting &set){
+        void from_setting ( const Setting &set ) {
                 Array<mpdf*> Atmp; //temporary Array
-                UI::get(Atmp,set, "mpdfs", UI::compulsory);
-                set_elements(Atmp,true);
-        }
-        
+                UI::get ( Atmp, set, "mpdfs", UI::compulsory );
+                set_elements ( Atmp, true );
+        }
+
 };
-UIREGISTER(mprod);
+UIREGISTER ( mprod );
 
 //! Product of independent epdfs. For dependent pdfs, use mprod.
@@ -344 +399 @@
         Array<datalink*> dls;
 public:
-        eprod () : epdfs ( 0 ),dls ( 0 ) {};
-        void set_parameters ( const Array<const epdf*> &epdfs0, bool named=true ) {
-                epdfs=epdfs0;//.set_length ( epdfs0.length() );
+        eprod () : epdfs ( 0 ), dls ( 0 ) {};
+        void set_parameters ( const Array<const epdf*> &epdfs0, bool named = true ) {
+                epdfs = epdfs0;//.set_length ( epdfs0.length() );
                 dls.set_length ( epdfs.length() );
 
-                bool independent=true;
+                bool independent = true;
                 if ( named ) {
-                        for ( int i=0;i<epdfs.length();i++ ) {
-                                independent=rv.add ( epdfs ( i )->_rv() );
-                                it_assert_debug ( independent==true, "eprod:: given components are not independent." );
+                        for ( int i = 0; i < epdfs.length(); i++ ) {
+                                independent = rv.add ( epdfs ( i )->_rv() );
+                                it_assert_debug ( independent == true, "eprod:: given components are not independent." );
                         }
-                        dim=rv._dsize();
-                }
-                else {
-                        dim =0; for ( int i=0;i<epdfs.length();i++ ) {
-                                dim+=epdfs ( i )->dimension();
+                        dim = rv._dsize();
+                } else {
+                        dim = 0;
+                        for ( int i = 0; i < epdfs.length(); i++ ) {
+                                dim += epdfs ( i )->dimension();
                         }
                 }
                 //
-                int cumdim=0;
-                int dimi=0;
+                int cumdim = 0;
+                int dimi = 0;
                 int i;
-                for ( i=0;i<epdfs.length();i++ ) {
+                for ( i = 0; i < epdfs.length(); i++ ) {
                         dls ( i ) = new datalink;
-                        if ( named ) {dls ( i )->set_connection ( epdfs ( i )->_rv() , rv );}
-                        else {
+                        if ( named ) {
+                                dls ( i )->set_connection ( epdfs ( i )->_rv() , rv );
+                        } else {
                                 dimi = epdfs ( i )->dimension();
-                                dls ( i )->set_connection ( dimi, dim, linspace ( cumdim,cumdim+dimi-1 ) );
-                                cumdim+=dimi;
+                                dls ( i )->set_connection ( dimi, dim, linspace ( cumdim, cumdim + dimi - 1 ) );
+                                cumdim += dimi;
                         }
                 }
@@ -379 +435 @@
         vec mean() const {
                 vec tmp ( dim );
-                for ( int i=0;i<epdfs.length();i++ ) {
+                for ( int i = 0; i < epdfs.length(); i++ ) {
                         vec pom = epdfs ( i )->mean();
                         dls ( i )->pushup ( tmp, pom );
@@ -387 +443 @@
         vec variance() const {
                 vec tmp ( dim ); //second moment
-                for ( int i=0;i<epdfs.length();i++ ) {
+                for ( int i = 0; i < epdfs.length(); i++ ) {
                         vec pom = epdfs ( i )->mean();
-                        dls ( i )->pushup ( tmp, pow ( pom,2 ) );
-                }
-                return tmp-pow ( mean(),2 );
+                        dls ( i )->pushup ( tmp, pow ( pom, 2 ) );
+                }
+                return tmp - pow ( mean(), 2 );
         }
         vec sample() const {
                 vec tmp ( dim );
-                for ( int i=0;i<epdfs.length();i++ ) {
+                for ( int i = 0; i < epdfs.length(); i++ ) {
                         vec pom = epdfs ( i )->sample();
                         dls ( i )->pushup ( tmp, pom );
@@ -402 +458 @@
         }
         double evallog ( const vec &val ) const {
-                double tmp=0;
-                for ( int i=0;i<epdfs.length();i++ ) {
-                        tmp+=epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
-                }
-                it_assert_debug ( std::isfinite ( tmp ),"Infinite" );
+                double tmp = 0;
+                for ( int i = 0; i < epdfs.length(); i++ ) {
+                        tmp += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
+                }
+                it_assert_debug ( std::isfinite ( tmp ), "Infinite" );
                 return tmp;
         }
         //!access function
-        const epdf* operator () ( int i ) const {it_assert_debug ( i<epdfs.length(),"wrong index" );return epdfs ( i );}
+        const epdf* operator () ( int i ) const {
+                it_assert_debug ( i < epdfs.length(), "wrong index" );
+                return epdfs ( i );
+        }
 
         //!Destructor
-        ~eprod() {for ( int i=0;i<epdfs.length();i++ ) {delete dls ( i );}}
+        ~eprod() {
+                for ( int i = 0; i < epdfs.length(); i++ ) {
+                        delete dls ( i );
+                }
+        }
 };
 
@@ -430 +493 @@
 public:
         //!Default constructor
-        mmix():iepdf(new emix()) {
-            set_ep(iepdf);
+        mmix() : iepdf ( new emix() ) {
+                set_ep ( iepdf );
         }
 
@@ -438 +501 @@
                 Array<epdf*> Eps ( Coms.length() );
 
-                for ( int i = 0;i < Coms.length();i++ ) {
-                        Eps(i) = Coms(i)->e();
-                }
-
-                iepdf->set_parameters(w, Eps);
+                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 );}
+                for ( int i = 0; i < Coms.length(); i++ ) {
+                        Coms ( i )->condition ( cond );
+                }
         };
 };

library/bdm/stat/exp_family.cpp

@@ -4 +4 @@
 #include "exp_family.h"
 
-namespace bdm{
+namespace bdm {
 
 Uniform_RNG UniRNG;
@@ -12 +12 @@
 using std::cout;
 
-void BMEF::bayes ( const vec &dt ) {this->bayes ( dt,1.0 );};
+void BMEF::bayes ( const vec &dt ) {
+        this->bayes ( dt, 1.0 );
+};
 
 vec egiw::sample() const {
@@ -20 +22 @@
 
 double egiw::evallog_nn ( const vec &val ) const {
-        int vend = val.length()-1;
-
-        if ( dimx==1 ) { //same as the following, just quicker.
+        int vend = val.length() - 1;
+
+        if ( dimx == 1 ) { //same as the following, just quicker.
                 double r = val ( vend ); //last entry!
-                if (r<0) return -inf;
-                vec Psi ( nPsi+dimx );
+                if ( r < 0 ) return -inf;
+                vec Psi ( nPsi + dimx );
                 Psi ( 0 ) = -1.0;
-                Psi.set_subvector ( 1,val ( 0,vend-1 ) ); // fill the rest
-
-                double Vq=V.qform ( Psi );
-                return -0.5* ( nu*log ( r ) + Vq /r );
-        }
-        else {
-                mat Th= reshape ( val ( 0,nPsi*dimx-1 ),nPsi,dimx );
-                fsqmat R ( reshape ( val ( nPsi*dimx,vend ),dimx,dimx ) );
-                double ldetR=R.logdet();
-                if (ldetR) return -inf;
-                mat Tmp=concat_vertical ( -eye ( dimx ),Th );
+                Psi.set_subvector ( 1, val ( 0, vend - 1 ) ); // fill the rest
+
+                double Vq = V.qform ( Psi );
+                return -0.5* ( nu*log ( r ) + Vq / r );
+        } else {
+                mat Th = reshape ( val ( 0, nPsi * dimx - 1 ), nPsi, dimx );
+                fsqmat R ( reshape ( val ( nPsi*dimx, vend ), dimx, dimx ) );
+                double ldetR = R.logdet();
+                if ( ldetR ) return -inf;
+                mat Tmp = concat_vertical ( -eye ( dimx ), Th );
                 fsqmat iR ( dimx );
                 R.inv ( iR );
@@ -48 +49 @@
         const vec& D = V._D();
 
-        double m = nu - nPsi -dimx-1;
+        double m = nu - nPsi - dimx - 1;
 #define log2  0.693147180559945286226763983
 #define logpi 1.144729885849400163877476189
@@ -54 +55 @@
 #define Inf std::numeric_limits<double>::infinity()
 
-        double nkG = 0.5* dimx* ( -nPsi *log2pi + sum ( log ( D ( dimx,D.length()-1 ) ) ) );
+        double nkG = 0.5 * dimx * ( -nPsi * log2pi + sum ( log ( D ( dimx, D.length() - 1 ) ) ) );
         // temporary for lgamma in Wishart
-        double lg=0;
-        for ( int i =0; i<dimx;i++ ) {lg+=lgamma ( 0.5* ( m-i ) );}
-
-        double nkW = 0.5* ( m*sum ( log ( D ( 0,dimx-1 ) ) ) ) \
-                     - 0.5*dimx* ( m*log2 + 0.5* ( dimx-1 ) *log2pi )  - lg;
-
-        it_assert_debug ( ( ( -nkG-nkW ) >-Inf ) && ( ( -nkG-nkW ) <Inf ), "ARX improper" );
-        return -nkG-nkW;
+        double lg = 0;
+        for ( int i = 0; i < dimx; i++ ) {
+                lg += lgamma ( 0.5 * ( m - i ) );
+        }
+
+        double nkW = 0.5 * ( m * sum ( log ( D ( 0, dimx - 1 ) ) ) ) \
+                     - 0.5 * dimx * ( m * log2 + 0.5 * ( dimx - 1 ) * log2pi )  - lg;
+
+        it_assert_debug ( ( ( -nkG - nkW ) > -Inf ) && ( ( -nkG - nkW ) < Inf ), "ARX improper" );
+        return -nkG - nkW;
 }
 
 vec egiw::est_theta() const {
-        if ( dimx==1 ) {
+        if ( dimx == 1 ) {
                 const mat &L = V._L();
                 int end = L.rows() - 1;
 
-                mat iLsub = ltuinv ( L ( dimx,end,dimx,end ) );
-
-                vec L0 = L.get_col(0);
-
-                return iLsub * L0(1,end);
-        }
-        else {
-                it_error("ERROR: est_theta() not implemented for dimx>1");
+                mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) );
+
+                vec L0 = L.get_col ( 0 );
+
+                return iLsub * L0 ( 1, end );
+        } else {
+                it_error ( "ERROR: est_theta() not implemented for dimx>1" );
                 return 0;
         }
@@ -89 +91 @@
                 int end = D.length() - 1;
 
-                mat Lsub = L ( 1, end, 1, end);
-                mat Dsub = diag(D(1, end));
-
-                return inv( transpose(Lsub) * Dsub * Lsub );
-
-        }
-        else {
-                it_error("ERROR: est_theta_cov() not implemented for dimx>1");
+                mat Lsub = L ( 1, end, 1, end );
+                mat Dsub = diag ( D ( 1, end ) );
+
+                return inv ( transpose ( Lsub ) * Dsub * Lsub );
+
+        } else {
+                it_error ( "ERROR: est_theta_cov() not implemented for dimx>1" );
                 return 0;
         }
@@ -104 +105 @@
 vec egiw::mean() const {
 
-        if ( dimx==1 ) {
-                const vec &D= V._D();
-                int end = D.length()-1;
+        if ( dimx == 1 ) {
+                const vec &D = V._D();
+                int end = D.length() - 1;
 
                 vec m ( dim );
                 m.set_subvector ( 0, est_theta() );
-                m ( end ) = D ( 0 ) / ( nu -nPsi -2*dimx -2 );
+                m ( end ) = D ( 0 ) / ( nu - nPsi - 2 * dimx - 2 );
                 return m;
-        }
-        else {
+        } else {
                 mat M;
                 mat R;
-                mean_mat ( M,R );
-                return cvectorize ( concat_vertical ( M,R ) );
+                mean_mat ( M, R );
+                return cvectorize ( concat_vertical ( M, R ) );
         }
 
@@ -124 +124 @@
 vec egiw::variance() const {
 
-        if ( dimx==1 ) {
-                int l=V.rows();
-                const ldmat tmp(V,linspace(1,l-1));
-                ldmat itmp(l);
-                tmp.inv(itmp);
-                double cove = V._D() ( 0 ) / ( nu -nPsi -2*dimx -2 );
-                
-                vec var(l);
-                var.set_subvector(0,diag(itmp.to_mat())*cove);
-                var(l-1)=cove*cove/( nu -nPsi -2*dimx -2 );
+        if ( dimx == 1 ) {
+                int l = V.rows();
+                const ldmat tmp ( V, linspace ( 1, l - 1 ) );
+                ldmat itmp ( l );
+                tmp.inv ( itmp );
+                double cove = V._D() ( 0 ) / ( nu - nPsi - 2 * dimx - 2 );
+
+                vec var ( l );
+                var.set_subvector ( 0, diag ( itmp.to_mat() ) *cove );
+                var ( l - 1 ) = cove * cove / ( nu - nPsi - 2 * dimx - 2 );
                 return var;
-        }
-        else {it_error("not implemented"); return vec(0);}
+        } else {
+                it_error ( "not implemented" );
+                return vec ( 0 );
+        }
 }
 
 void egiw::mean_mat ( mat &M, mat&R ) const {
-        const mat &L= V._L();
-        const vec &D= V._D();
-        int end = L.rows()-1;
-
-        ldmat ldR ( L ( 0,dimx-1,0,dimx-1 ), D ( 0,dimx-1 ) / ( nu -nPsi -2*dimx -2 ) ); //exp val of R
-        mat iLsub=ltuinv ( L ( dimx,end,dimx,end ) );
+        const mat &L = V._L();
+        const vec &D = V._D();
+        int end = L.rows() - 1;
+
+        ldmat ldR ( L ( 0, dimx - 1, 0, dimx - 1 ), D ( 0, dimx - 1 ) / ( nu - nPsi - 2*dimx - 2 ) ); //exp val of R
+        mat iLsub = ltuinv ( L ( dimx, end, dimx, end ) );
 
         // set mean value
-        mat Lpsi = L ( dimx,end,0,dimx-1 );
-        M= iLsub*Lpsi;
-        R= ldR.to_mat()  ;
+        mat Lpsi = L ( dimx, end, 0, dimx - 1 );
+        M = iLsub * Lpsi;
+        R = ldR.to_mat()  ;
 }
 
 vec egamma::sample() const {
-        vec smp ( dim);
+        vec smp ( dim );
         int i;
 
-        for ( i=0; i<dim; i++ ) {
-                if ( beta ( i ) >std::numeric_limits<double>::epsilon() ) {
-                        GamRNG.setup ( alpha ( i ),beta ( i ) );
-                }
-                else {
-                        GamRNG.setup ( alpha ( i ),std::numeric_limits<double>::epsilon() );
+        for ( i = 0; i < dim; i++ ) {
+                if ( beta ( i ) > std::numeric_limits<double>::epsilon() ) {
+                        GamRNG.setup ( alpha ( i ), beta ( i ) );
+                } else {
+                        GamRNG.setup ( alpha ( i ), std::numeric_limits<double>::epsilon() );
                 }
 #pragma omp critical
@@ -190 +191 @@
         int i;
 
-        if (any(val<=0.)) return -inf;
-        if (any(beta<=0.)) return -inf;
-        for ( i=0; i<dim; i++ ) {
-                res += ( alpha ( i ) - 1 ) *std::log ( val ( i ) ) - beta ( i ) *val ( i );
-        }
-        double tmp=res-lognc();;
-        it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+        if ( any ( val <= 0. ) ) return -inf;
+        if ( any ( beta <= 0. ) ) return -inf;
+        for ( i = 0; i < dim; i++ ) {
+                res += ( alpha ( i ) - 1 ) * std::log ( val ( i ) ) - beta ( i ) * val ( i );
+        }
+        double tmp = res - lognc();;
+        it_assert_debug ( std::isfinite ( tmp ), "Infinite value" );
         return tmp;
 }
@@ -204 +205 @@
         int i;
 
-        for ( i=0; i<dim; i++ ) {
-                res += lgamma ( alpha ( i ) ) - alpha ( i ) *std::log ( beta ( i ) ) ;
+        for ( i = 0; i < dim; i++ ) {
+                res += lgamma ( alpha ( i ) ) - alpha ( i ) * std::log ( beta ( i ) ) ;
         }
 
@@ -211 +212 @@
 }
 
-void mgamma::set_parameters(double k0, const vec &beta0) {
+void mgamma::set_parameters ( double k0, const vec &beta0 ) {
         k = k0;
-        iepdf->set_parameters(k * ones(beta0.length()), beta0);
+        iepdf->set_parameters ( k * ones ( beta0.length() ), beta0 );
         dimc = e()->dimension();
 }
 
-ivec eEmp::resample (RESAMPLING_METHOD method) {
-        ivec ind=zeros_i ( n );
+ivec eEmp::resample ( RESAMPLING_METHOD method ) {
+        ivec ind = zeros_i ( n );
         ivec N_babies = zeros_i ( n );
         vec cumDist = cumsum ( w );
         vec u ( n );
-        int i,j,parent;
+        int i, j, parent;
         double u0;
 
         switch ( method ) {
-                case MULTINOMIAL:
-                        u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );
-
-                        for ( i = n - 2;i >= 0;i-- ) {
-                                u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
-                        }
-
-                        break;
-
-                case STRATIFIED:
-
-                        for ( i = 0;i < n;i++ ) {
-                                u ( i ) = ( i + UniRNG.sample() ) / n;
-                        }
-
-                        break;
-
-                case SYSTEMATIC:
-                        u0 = UniRNG.sample();
-
-                        for ( i = 0;i < n;i++ ) {
-                                u ( i ) = ( i + u0 ) / n;
-                        }
-
-                        break;
-
-                default:
-                        it_error ( "PF::resample(): Unknown resampling method" );
+        case MULTINOMIAL:
+                u ( n - 1 ) = pow ( UniRNG.sample(), 1.0 / n );
+
+                for ( i = n - 2; i >= 0; i-- ) {
+                        u ( i ) = u ( i + 1 ) * pow ( UniRNG.sample(), 1.0 / ( i + 1 ) );
+                }
+
+                break;
+
+        case STRATIFIED:
+
+                for ( i = 0; i < n; i++ ) {
+                        u ( i ) = ( i + UniRNG.sample() ) / n;
+                }
+
+                break;
+
+        case SYSTEMATIC:
+                u0 = UniRNG.sample();
+
+                for ( i = 0; i < n; i++ ) {
+                        u ( i ) = ( i + u0 ) / n;
+                }
+
+                break;
+
+        default:
+                it_error ( "PF::resample(): Unknown resampling method" );
         }
 
@@ -259 +260 @@
         j = 0;
 
-        for ( i = 0;i < n;i++ ) {
+        for ( i = 0; i < n; i++ ) {
                 while ( u ( i ) > cumDist ( j ) ) j++;
 
@@ -270 +271 @@
 
         // find the first parent;
-        parent=0; while ( N_babies ( parent ) ==0 ) parent++;
+        parent = 0;
+        while ( N_babies ( parent ) == 0 ) parent++;
 
         // Build index
-        for ( i = 0;i < n;i++ ) {
+        for ( i = 0; i < n; i++ ) {
                 if ( N_babies ( i ) > 0 ) {
                         ind ( i ) = i;
                         N_babies ( i ) --; //this index was now replicated;
-                }
-                else {
+                } else {
                         // test if the parent has been fully replicated
                         // if yes, find the next one
-                        while ( ( N_babies ( parent ) ==0 ) || ( N_babies ( parent ) ==1 && parent>i ) ) parent++;
+                        while ( ( N_babies ( parent ) == 0 ) || ( N_babies ( parent ) == 1 && parent > i ) ) parent++;
 
                         // Replicate parent
@@ -292 +293 @@
 
         // copy the internals according to ind
-        for ( i=0;i<n;i++ ) {
-                if ( ind ( i ) !=i ) {
-                        samples ( i ) =samples ( ind ( i ) );
-                }
-                w ( i ) = 1.0/n;
+        for ( i = 0; i < n; i++ ) {
+                if ( ind ( i ) != i ) {
+                        samples ( i ) = samples ( ind ( i ) );
+                }
+                w ( i ) = 1.0 / n;
         }
 
@@ -305 +306 @@
         //it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
         dim = epdf0->dimension();
-        w=w0;
-        w/=sum ( w0 );//renormalize
-        n=w.length();
+        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();}
+        for ( int i = 0; i < n; i++ ) {
+                samples ( i ) = epdf0->sample();
+        }
 }
 
 void eEmp::set_samples ( const epdf* epdf0 ) {
         //it_assert_debug(rv==epdf0->rv(),"Wrong epdf0");
-        w=1;
-        w/=sum ( w );//renormalize
-
-        for ( int i=0;i<n;i++ ) {samples ( i ) =epdf0->sample();}
-}
-
-void migamma_ref::from_setting( const Setting &set ) 
-{                               
+        w = 1;
+        w /= sum ( w );//renormalize
+
+        for ( int i = 0; i < n; i++ ) {
+                samples ( i ) = epdf0->sample();
+        }
+}
+
+void migamma_ref::from_setting ( const Setting &set ) {
         vec ref;
-        UI::get( ref, set, "ref" , UI::compulsory);
-        set_parameters(set["k"],ref,set["l"]);
-}
-
-void mlognorm::from_setting( const Setting &set ) 
-{       
-        vec mu0;        
-        UI::get( mu0, set, "mu0", UI::compulsory);
-        set_parameters(mu0.length(),set["k"]);
-        condition(mu0);
+        UI::get ( ref, set, "ref" , UI::compulsory );
+        set_parameters ( set["k"], ref, set["l"] );
+}
+
+void mlognorm::from_setting ( const Setting &set ) {
+        vec mu0;
+        UI::get ( mu0, set, "mu0", UI::compulsory );
+        set_parameters ( mu0.length(), set["k"] );
+        condition ( mu0 );
 }
 

library/bdm/stat/merger.cpp

@@ -3 +3 @@
 #include "../estim/arx.h"
 
-namespace bdm
-{
-
-        merger_base::merger_base(const Array<mpdf*> &S, bool own) {
-            DBG = false;
-            dbg_file = NULL;
-            set_sources(S, own);
-        }
-
-        vec merger_base::merge_points ( mat &lW ) {
-                int nu=lW.rows();
-                vec result;
-                ivec indW;
-                bool infexist;
-                switch ( METHOD ) {
-                        case ARITHMETIC:
-                                result= log ( sum ( exp ( lW ) ) ); //ugly!
-                                break;
-                        case GEOMETRIC:
-                                result= sum ( lW ) /nu;
-                                break;
-                        case LOGNORMAL:
-                                vec sumlW=sum ( lW ) ;
-                                indW=find((sumlW<inf) & (sumlW>-inf));
-                                infexist=(indW.size()<lW.cols());
-                                vec mu;
-                                vec lam;
-                                if (infexist){
-                                        mu = sumlW(indW) /nu; //mean of logs
-                                        // 
-                                        mat validlW=lW.get_cols(indW);
-                                        lam = sum ( pow ( validlW-outer_product ( ones ( validlW.rows() ),mu ),2 ) );
-                                }
-                                else {
-                                        mu = sum ( lW ) /nu; //mean of logs
-                                        lam = sum ( pow ( lW-outer_product ( ones ( lW.rows() ),mu ),2 ) );
-                                }
-                                //
-                                double coef=0.0;
-                                vec sq2bl=sqrt ( 2*beta*lam ); //this term is everywhere
-                                switch ( nu ) {
-                                        case 2:
-                                                coef= ( 1-0.5*sqrt ( ( 4.0*beta-3.0 ) /beta ) );
-                                                result =coef*sq2bl + mu ;
-                                                break;
-                                        // case 4: == can be done similar to case 2 - is it worth it???
-                                        default: // see accompanying document merge_lognorm_derivation.lyx
-                                                coef = sqrt(1-(nu+1)/(2*beta*nu));
-                                                result =log(besselk((nu-3)/2, sq2bl*coef))-log(besselk((nu-3)/2, sq2bl)) + mu;
-                                                break;
-                                }
-                                break;
-                }
-                if (infexist){
-                        vec tmp =-inf*ones(lW.cols());
-                        set_subvector(tmp, indW, result);
-                        return tmp;
-                }
-                else {return result;}
-        }
-
-        void merger_mix::merge ( )
-        {
-                Array<vec> &Smp = eSmp._samples(); //aux
-                vec &w = eSmp._w(); //aux
-
-                mat Smp_ex =ones ( dim +1,Npoints ); // Extended samples for the ARX model - the last row is ones
-                for ( int i=0;i<Npoints;i++ ) { set_col_part ( Smp_ex,i,Smp ( i ) );}
-
-                if ( DBG )      *dbg_file << Name ( "Smp_0" ) << Smp_ex;
-
-                // Stuff for merging
-                vec lw_src ( Npoints );         // weights of the ith source
-                vec lw_mix ( Npoints );         // weights of the approximating mixture
-                vec lw ( Npoints );                     // tmp
-                mat lW=zeros ( Nsources,Npoints ); // array of weights of all sources
-                vec vec0 ( 0 );
-
-                //initialize importance weights
-                lw_mix = 1.0; // assuming uniform grid density -- otherwise 
-
-                // Initial component in the mixture model
-                mat V0=1e-8*eye ( dim +1 );
-                ARX A0; A0.set_statistics ( dim, V0 ); //initial guess of Mix:
-
-                Mix.init ( &A0, Smp_ex, Ncoms );
-                //Preserve initial mixture for repetitive estimation via flattening
-                MixEF Mix_init ( Mix );
-
-                // ============= MAIN LOOP ==================
-                bool converged=false;
-                int niter = 0;
-                char dbg_str[100];
-
-                emix* Mpred=Mix.epredictor ( );
-                vec Mix_pdf ( Npoints );
-                while ( !converged )
-                {
-                        //Re-estimate Mix
-                        //Re-Initialize Mixture model
-                        Mix.flatten ( &Mix_init );
-                        Mix.bayesB ( Smp_ex, w*Npoints );
-                        delete Mpred;
-                        Mpred = Mix.epredictor ( ); // Allocation => must be deleted at the end!!
-                        Mpred->set_rv ( rv ); //the predictor predicts rv of this merger
-
-                        // This will be active only later in iterations!!!
-                        if (  1./sum_sqr ( w ) <effss_coef*Npoints ) 
-                        {
-                                // Generate new samples
-                                eSmp.set_samples ( Mpred );
-                                for ( int i=0;i<Npoints;i++ )
-                                {
-                                        //////////// !!!!!!!!!!!!!
-                                        //if ( Smp ( i ) ( 2 ) <0 ) {Smp ( i ) ( 2 ) = 0.01; }
-                                        set_col_part ( Smp_ex,i,Smp ( i ) );
-                                        //Importance of the mixture
-                                        //lw_mix ( i ) =Mix.logpred (Smp_ex.get_col(i) );
-                                        lw_mix ( i ) = Mpred->evallog ( Smp ( i ) );
-                                }
-                                if ( DBG )
-                                {
-                                        cout<<"Resampling =" << 1./sum_sqr ( w ) << endl;
-                                        cout << Mix._e()->mean() <<endl;
-                                        cout << sum ( Smp_ex,2 ) /Npoints <<endl;
-                                        cout << Smp_ex*Smp_ex.T() /Npoints << endl;
-                                }
-                        }
-                        if ( DBG )
-                        {
-                                sprintf ( dbg_str,"Mpred_mean%d",niter );
-                                *dbg_file << Name ( dbg_str ) << Mpred->mean();
-                                sprintf ( dbg_str,"Mpred_var%d",niter );
-                                *dbg_file << Name ( dbg_str ) << Mpred->variance();
-
-
-                                sprintf ( dbg_str,"Mpdf%d",niter );
-                                for ( int i=0;i<Npoints;i++ ) {Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );}
-                                *dbg_file << Name ( dbg_str ) << Mix_pdf;
-
-                                sprintf ( dbg_str,"Smp%d",niter );
-                                *dbg_file << Name ( dbg_str ) << Smp_ex;
-
-                        }
-                        //Importace weighting
-                        for ( int i=0;i<mpdfs.length();i++ )
-                        {
-                                lw_src=0.0;
-                                //======== Same RVs ===========
-                                //Split according to dependency in rvs
-                                if ( mpdfs ( i )->dimension() ==dim )
-                                {
-                                        // no need for conditioning or marginalization
-                                        lw_src = mpdfs ( i )->e()->evallog_m ( Smp  );
-                                }
-                                else
-                                {
-                                        // compute likelihood of marginal on the conditional variable
-                                        if ( mpdfs ( i )->dimensionc() >0 )
-                                        {
-                                                // Make marginal on rvc_i
-                                                epdf* tmp_marg = Mpred->marginal ( mpdfs ( i )->_rvc() );
-                                                //compute vector of lw_src
-                                                for ( int k=0;k<Npoints;k++ )
-                                                {
-                                                        // Here val of tmp_marg = cond of mpdfs(i) ==> calling dls->get_cond
-                                                        lw_src ( k ) += tmp_marg->evallog ( dls ( i )->get_cond ( Smp ( k ) ) );
-                                                }
-                                                delete tmp_marg;
+namespace bdm {
+
+merger_base::merger_base ( const Array<mpdf*> &S, bool own ) {
+        DBG = false;
+        dbg_file = NULL;
+        set_sources ( S, own );
+}
+
+vec merger_base::merge_points ( mat &lW ) {
+        int nu = lW.rows();
+        vec result;
+        ivec indW;
+        bool infexist;
+        switch ( METHOD ) {
+        case ARITHMETIC:
+                result = log ( sum ( exp ( lW ) ) ); //ugly!
+                break;
+        case GEOMETRIC:
+                result = sum ( lW ) / nu;
+                break;
+        case LOGNORMAL:
+                vec sumlW = sum ( lW ) ;
+                indW = find ( ( sumlW < inf ) & ( sumlW > -inf ) );
+                infexist = ( indW.size() < lW.cols() );
+                vec mu;
+                vec lam;
+                if ( infexist ) {
+                        mu = sumlW ( indW ) / nu; //mean of logs
+                        //
+                        mat validlW = lW.get_cols ( indW );
+                        lam = sum ( pow ( validlW - outer_product ( ones ( validlW.rows() ), mu ), 2 ) );
+                } else {
+                        mu = sum ( lW ) / nu; //mean of logs
+                        lam = sum ( pow ( lW - outer_product ( ones ( lW.rows() ), mu ), 2 ) );
+                }
+                //
+                double coef = 0.0;
+                vec sq2bl = sqrt ( 2 * beta * lam ); //this term is everywhere
+                switch ( nu ) {
+                case 2:
+                        coef = ( 1 - 0.5 * sqrt ( ( 4.0 * beta - 3.0 ) / beta ) );
+                        result = coef * sq2bl + mu ;
+                        break;
+                        // case 4: == can be done similar to case 2 - is it worth it???
+                default: // see accompanying document merge_lognorm_derivation.lyx
+                        coef = sqrt ( 1 - ( nu + 1 ) / ( 2 * beta * nu ) );
+                        result = log ( besselk ( ( nu - 3 ) / 2, sq2bl * coef ) ) - log ( besselk ( ( nu - 3 ) / 2, sq2bl ) ) + mu;
+                        break;
+                }
+                break;
+        }
+        if ( infexist ) {
+                vec tmp = -inf * ones ( lW.cols() );
+                set_subvector ( tmp, indW, result );
+                return tmp;
+        } else {
+                return result;
+        }
+}
+
+void merger_mix::merge ( ) {
+        Array<vec> &Smp = eSmp._samples(); //aux
+        vec &w = eSmp._w(); //aux
+
+        mat Smp_ex = ones ( dim + 1, Npoints ); // Extended samples for the ARX model - the last row is ones
+        for ( int i = 0; i < Npoints; i++ ) {
+                set_col_part ( Smp_ex, i, Smp ( i ) );
+        }
+
+        if ( DBG )      *dbg_file << Name ( "Smp_0" ) << Smp_ex;
+
+        // Stuff for merging
+        vec lw_src ( Npoints );         // weights of the ith source
+        vec lw_mix ( Npoints );         // weights of the approximating mixture
+        vec lw ( Npoints );                     // tmp
+        mat lW = zeros ( Nsources, Npoints ); // array of weights of all sources
+        vec vec0 ( 0 );
+
+        //initialize importance weights
+        lw_mix = 1.0; // assuming uniform grid density -- otherwise
+
+        // Initial component in the mixture model
+        mat V0 = 1e-8 * eye ( dim + 1 );
+        ARX A0;
+        A0.set_statistics ( dim, V0 ); //initial guess of Mix:
+
+        Mix.init ( &A0, Smp_ex, Ncoms );
+        //Preserve initial mixture for repetitive estimation via flattening
+        MixEF Mix_init ( Mix );
+
+        // ============= MAIN LOOP ==================
+        bool converged = false;
+        int niter = 0;
+        char dbg_str[100];
+
+        emix* Mpred = Mix.epredictor ( );
+        vec Mix_pdf ( Npoints );
+        while ( !converged ) {
+                //Re-estimate Mix
+                //Re-Initialize Mixture model
+                Mix.flatten ( &Mix_init );
+                Mix.bayesB ( Smp_ex, w*Npoints );
+                delete Mpred;
+                Mpred = Mix.epredictor ( ); // Allocation => must be deleted at the end!!
+                Mpred->set_rv ( rv ); //the predictor predicts rv of this merger
+
+                // This will be active only later in iterations!!!
+                if ( 1. / sum_sqr ( w ) < effss_coef*Npoints ) {
+                        // Generate new samples
+                        eSmp.set_samples ( Mpred );
+                        for ( int i = 0; i < Npoints; i++ ) {
+                                //////////// !!!!!!!!!!!!!
+                                //if ( Smp ( i ) ( 2 ) <0 ) {Smp ( i ) ( 2 ) = 0.01; }
+                                set_col_part ( Smp_ex, i, Smp ( i ) );
+                                //Importance of the mixture
+                                //lw_mix ( i ) =Mix.logpred (Smp_ex.get_col(i) );
+                                lw_mix ( i ) = Mpred->evallog ( Smp ( i ) );
+                        }
+                        if ( DBG ) {
+                                cout << "Resampling =" << 1. / sum_sqr ( w ) << endl;
+                                cout << Mix._e()->mean() << endl;
+                                cout << sum ( Smp_ex, 2 ) / Npoints << endl;
+                                cout << Smp_ex*Smp_ex.T() / Npoints << endl;
+                        }
+                }
+                if ( DBG ) {
+                        sprintf ( dbg_str, "Mpred_mean%d", niter );
+                        *dbg_file << Name ( dbg_str ) << Mpred->mean();
+                        sprintf ( dbg_str, "Mpred_var%d", niter );
+                        *dbg_file << Name ( dbg_str ) << Mpred->variance();
+
+
+                        sprintf ( dbg_str, "Mpdf%d", niter );
+                        for ( int i = 0; i < Npoints; i++ ) {
+                                Mix_pdf ( i ) = Mix.logpred ( Smp_ex.get_col ( i ) );
+                        }
+                        *dbg_file << Name ( dbg_str ) << Mix_pdf;
+
+                        sprintf ( dbg_str, "Smp%d", niter );
+                        *dbg_file << Name ( dbg_str ) << Smp_ex;
+
+                }
+                //Importace weighting
+                for ( int i = 0; i < mpdfs.length(); i++ ) {
+                        lw_src = 0.0;
+                        //======== Same RVs ===========
+                        //Split according to dependency in rvs
+                        if ( mpdfs ( i )->dimension() == dim ) {
+                                // no need for conditioning or marginalization
+                                lw_src = mpdfs ( i )->e()->evallog_m ( Smp );
+                        } else {
+                                // compute likelihood of marginal on the conditional variable
+                                if ( mpdfs ( i )->dimensionc() > 0 ) {
+                                        // Make marginal on rvc_i
+                                        epdf* tmp_marg = Mpred->marginal ( mpdfs ( i )->_rvc() );
+                                        //compute vector of lw_src
+                                        for ( int k = 0; k < Npoints; k++ ) {
+                                                // Here val of tmp_marg = cond of mpdfs(i) ==> calling dls->get_cond
+                                                lw_src ( k ) += tmp_marg->evallog ( dls ( i )->get_cond ( Smp ( k ) ) );
+                                        }
+                                        delete tmp_marg;
 
 //                                      sprintf ( str,"marg%d",niter );
 //                                      *dbg << Name ( str ) << lw_src;
 
+                                }
+                                // Compute likelihood of the missing variable
+                                if ( dim > ( mpdfs ( i )->dimension() + mpdfs ( i )->dimensionc() ) ) {
+                                        ///////////////
+                                        // There are variales unknown to mpdfs(i) : rvzs
+                                        mpdf* tmp_cond = Mpred->condition ( rvzs ( i ) );
+                                        // Compute likelihood
+                                        vec lw_dbg = lw_src;
+                                        for ( int k = 0; k < Npoints; k++ ) {
+                                                lw_src ( k ) += log (
+                                                                    tmp_cond->evallogcond (
+                                                                        zdls ( i )->pushdown ( Smp ( k ) ),
+                                                                        zdls ( i )->get_cond ( Smp ( k ) ) ) );
+                                                if ( !std::isfinite ( lw_src ( k ) ) ) {
+                                                        lw_src ( k ) = -1e16;
+                                                        cout << "!";
+                                                }
                                         }
-                                        // Compute likelihood of the missing variable
-                                        if ( dim > ( mpdfs ( i )->dimension() + mpdfs ( i )->dimensionc() ) )
-                                        {
-                                                ///////////////
-                                                // There are variales unknown to mpdfs(i) : rvzs
-                                                mpdf* tmp_cond = Mpred->condition ( rvzs ( i ) );
-                                                // Compute likelihood
-                                                vec lw_dbg=lw_src;
-                                                for ( int k= 0; k<Npoints; k++ )
-                                                {
-                                                        lw_src ( k ) += log (
-                                                                            tmp_cond->evallogcond (
-                                                                                zdls ( i )->pushdown ( Smp ( k ) ),
-                                                                                zdls ( i )->get_cond ( Smp ( k ) ) ) );
-                                                        if ( !std::isfinite ( lw_src ( k ) ) )
-                                                        {
-                                                                lw_src ( k ) = -1e16; cout << "!";
-                                                        }
-                                                }
-                                                delete tmp_cond;
-                                        }
-                                        // Compute likelihood of the partial source
-                                        for ( int k= 0; k<Npoints; k++ )
-                                        {
-                                                mpdfs ( i )->condition ( dls ( i )->get_cond ( Smp ( k ) ) );
-                                                lw_src ( k ) += mpdfs ( i )->e()->evallog ( dls ( i )->pushdown ( Smp ( k ) ) );
-                                        }
-
+                                        delete tmp_cond;
                                 }
-                //                      it_assert_debug(std::isfinite(sum(lw_src)),"bad");
-                                lW.set_row ( i, lw_src ); // do not divide by mix
-                        }
-                        lw = merger_base::merge_points ( lW ); //merge
-
-                        //Importance weighting
-                        lw -=  lw_mix; // hoping that it is not numerically sensitive...
-                        w = exp ( lw-max ( lw ) );
-
-                        //renormalize
-                        double sumw =sum ( w );
-                        if ( std::isfinite ( sumw ) )
-                        {
-                                w = w/sumw;
-                        }
-                        else
-                        {
-                                it_file itf ( "merg_err.it" );
-                                itf << Name ( "w" ) << w;
-                        }
-
-                        if ( DBG )
-                        {
-                                sprintf ( dbg_str,"lW%d",niter );
-                                *dbg_file << Name ( dbg_str ) << lW;
-                                sprintf ( dbg_str,"w%d",niter );
-                                *dbg_file << Name ( dbg_str ) << w;
-                                sprintf ( dbg_str,"lw_m%d",niter );
-                                *dbg_file << Name ( dbg_str ) << lw_mix;
-                        }
-                        // ==== stopping rule ===
-                        niter++;
-                        converged = ( niter>stop_niter );
-                }
-                delete Mpred;
+                                // Compute likelihood of the partial source
+                                for ( int k = 0; k < Npoints; k++ ) {
+                                        mpdfs ( i )->condition ( dls ( i )->get_cond ( Smp ( k ) ) );
+                                        lw_src ( k ) += mpdfs ( i )->e()->evallog ( dls ( i )->pushdown ( Smp ( k ) ) );
+                                }
+
+                        }
+                        //                      it_assert_debug(std::isfinite(sum(lw_src)),"bad");
+                        lW.set_row ( i, lw_src ); // do not divide by mix
+                }
+                lw = merger_base::merge_points ( lW ); //merge
+
+                //Importance weighting
+                lw -=  lw_mix; // hoping that it is not numerically sensitive...
+                w = exp ( lw - max ( lw ) );
+
+                //renormalize
+                double sumw = sum ( w );
+                if ( std::isfinite ( sumw ) ) {
+                        w = w / sumw;
+                } else {
+                        it_file itf ( "merg_err.it" );
+                        itf << Name ( "w" ) << w;
+                }
+
+                if ( DBG ) {
+                        sprintf ( dbg_str, "lW%d", niter );
+                        *dbg_file << Name ( dbg_str ) << lW;
+                        sprintf ( dbg_str, "w%d", niter );
+                        *dbg_file << Name ( dbg_str ) << w;
+                        sprintf ( dbg_str, "lw_m%d", niter );
+                        *dbg_file << Name ( dbg_str ) << lw_mix;
+                }
+                // ==== stopping rule ===
+                niter++;
+                converged = ( niter > stop_niter );
+        }
+        delete Mpred;
 //              cout << endl;
 
-        }
-
-        // DEFAULTS FOR MERGER_BASE 
-        const MERGER_METHOD merger_base::DFLT_METHOD=LOGNORMAL;
-        const double merger_base::DFLT_beta=1.2;
-        // DEFAULTS FOR MERGER_MIX
-        const int merger_mix::DFLT_Ncoms=10;
-        const double merger_mix::DFLT_effss_coef=0.5;
-
-}
+}
+
+// DEFAULTS FOR MERGER_BASE
+const MERGER_METHOD merger_base::DFLT_METHOD = LOGNORMAL;
+const double merger_base::DFLT_beta = 1.2;
+// DEFAULTS FOR MERGER_MIX
+const int merger_mix::DFLT_Ncoms = 10;
+const double merger_mix::DFLT_effss_coef = 0.5;
+
+}

library/bdm/stat/merger.h

@@ -17 +17 @@
 #include "../estim/mixtures.h"
 
-namespace bdm
-{
+namespace bdm {
 using std::string;
 
@@ -43 +42 @@
 */
 
-class merger_base : public compositepdf, public epdf
-{
-        protected:
-                //! Data link for each mpdf in mpdfs
-                Array<datalink_m2e*> dls;
-                //! Array of rvs that are not modelled by mpdfs at all, \f$ z_i \f$
-                Array<RV> rvzs;
-                //! Data Links for extension \f$ f(z_i|x_i,y_i) \f$
-                Array<datalink_m2e*> zdls;
-                //! number of support points
-                int Npoints;
-                //! number of sources
-                int Nsources;
-
-                //! switch of the methoh used for merging
-                MERGER_METHOD METHOD;
-                //! Default for METHOD
-                static const MERGER_METHOD DFLT_METHOD;
-                
-                //!Prior on the log-normal merging model
-                double beta;
-                //! default for beta
-                static const double DFLT_beta;
-                
-                //! Projection to empirical density (could also be piece-wise linear)
-                eEmp eSmp;
-
-                //! debug or not debug
-                bool DBG;
-
-                //! debugging file
-                it_file* dbg_file;
-        public:
-                //! \name Constructors
-                //! @{
-
-                //!Empty constructor
-                merger_base () : compositepdf() {DBG = false;dbg_file = NULL;}
-
-                //!Constructor from sources
-                merger_base (const Array<mpdf*> &S, bool own=false);
-
-                //! Function setting the main internal structures
-                void set_sources (const Array<mpdf*> &Sources, bool own) {
-                        compositepdf::set_elements (Sources,own);
-                        Nsources=mpdfs.length();
-                        //set sizes
-                        dls.set_size (Sources.length());
-                        rvzs.set_size (Sources.length());
-                        zdls.set_size (Sources.length());
-
-                        rv = getrv (/* checkoverlap = */ false);
-                        RV rvc; setrvc (rv, rvc);  // Extend rv by rvc!
-                        // join rv and rvc - see descriprion
-                        rv.add (rvc);
-                        // get dimension
-                        dim = rv._dsize();
-
-                        // create links between sources and common rv
-                        RV xytmp;
-                        for (int i = 0;i < mpdfs.length();i++) {
-                                //Establich connection between mpdfs and merger
-                                dls (i) = new datalink_m2e;
-                                dls (i)->set_connection (mpdfs (i)->_rv(), mpdfs (i)->_rvc(), rv);
-
-                                // find out what is missing in each mpdf
-                                xytmp = mpdfs (i)->_rv();
-                                xytmp.add (mpdfs (i)->_rvc());
-                                // z_i = common_rv-xy
-                                rvzs (i) = rv.subt (xytmp);
-                                //establish connection between extension (z_i|x,y)s and common rv
-                                zdls (i) = new datalink_m2e; zdls (i)->set_connection (rvzs (i), xytmp, rv) ;
-                        };
-                }
-                //! Rectangular support  each vector of XYZ specifies (begining-end) interval for each dimension. Same number of points, \c dimsize, in each dimension. 
-                void set_support (const Array<vec> &XYZ, const int dimsize) {
-                        set_support(XYZ,dimsize*ones_i(XYZ.length()));
-                }
-                //! Rectangular support  each vector of XYZ specifies (begining-end) interval for each dimension. \c gridsize specifies number of points is each dimension.
-                void set_support (const Array<vec> &XYZ, const ivec &gridsize) {
-                        int dim = XYZ.length();  //check with internal dim!!
-                        Npoints = prod (gridsize); 
-                        eSmp.set_parameters (Npoints, false);
-                        Array<vec> &samples = eSmp._samples();
-                        eSmp._w() = ones (Npoints) / Npoints; //unifrom size of bins
-                        //set samples
-                        ivec ind = zeros_i (dim);      //indeces of dimensions in for cycle;
-                        vec smpi (dim);            // ith sample
-                        vec steps =zeros(dim);            // ith sample
-                        // first corner
-                        for (int j = 0; j < dim; j++) {
-                                smpi (j) = XYZ (j) (0); /* beginning of the interval*/ 
-                                it_assert(gridsize(j)!=0.0,"Zeros in gridsize!");
-                                steps (j) = ( XYZ(j)(1)-smpi(j) )/gridsize(j);
-                        }
-                        // fill samples
-                        for (int i = 0; i < Npoints; i++) {
-                                // copy 
-                                samples(i) = smpi; 
-                                // go through all dimensions
-                                for (int j = 0;j < dim;j++) {
-                                        if (ind (j) == gridsize (j) - 1) { //j-th index is full
-                                                ind (j) = 0; //shift back
-                                                smpi(j) = XYZ(j)(0);
-                                                
-                                                if (i<Npoints-1) {
-                                                        ind (j + 1) ++; //increase the next dimension;
-                                                        smpi(j+1) += steps(j+1);
-                                                        break;
-                                                }
-                                                
-                                        } else {
-                                                ind (j) ++; 
-                                                smpi(j) +=steps(j);
+class merger_base : public compositepdf, public epdf {
+protected:
+        //! Data link for each mpdf in mpdfs
+        Array<datalink_m2e*> dls;
+        //! Array of rvs that are not modelled by mpdfs at all, \f$ z_i \f$
+        Array<RV> rvzs;
+        //! Data Links for extension \f$ f(z_i|x_i,y_i) \f$
+        Array<datalink_m2e*> zdls;
+        //! number of support points
+        int Npoints;
+        //! number of sources
+        int Nsources;
+
+        //! switch of the methoh used for merging
+        MERGER_METHOD METHOD;
+        //! Default for METHOD
+        static const MERGER_METHOD DFLT_METHOD;
+
+        //!Prior on the log-normal merging model
+        double beta;
+        //! default for beta
+        static const double DFLT_beta;
+
+        //! Projection to empirical density (could also be piece-wise linear)
+        eEmp eSmp;
+
+        //! debug or not debug
+        bool DBG;
+
+        //! debugging file
+        it_file* dbg_file;
+public:
+        //! \name Constructors
+        //! @{
+
+        //!Empty constructor
+        merger_base () : compositepdf() {
+                DBG = false;
+                dbg_file = NULL;
+        }
+
+        //!Constructor from sources
+        merger_base ( const Array<mpdf*> &S, bool own = false );
+
+        //! Function setting the main internal structures
+        void set_sources ( const Array<mpdf*> &Sources, bool own ) {
+                compositepdf::set_elements ( Sources, own );
+                Nsources = mpdfs.length();
+                //set sizes
+                dls.set_size ( Sources.length() );
+                rvzs.set_size ( Sources.length() );
+                zdls.set_size ( Sources.length() );
+
+                rv = getrv ( /* checkoverlap = */ false );
+                RV rvc;
+                setrvc ( rv, rvc );  // Extend rv by rvc!
+                // join rv and rvc - see descriprion
+                rv.add ( rvc );
+                // get dimension
+                dim = rv._dsize();
+
+                // create links between sources and common rv
+                RV xytmp;
+                for ( int i = 0; i < mpdfs.length(); i++ ) {
+                        //Establich connection between mpdfs and merger
+                        dls ( i ) = new datalink_m2e;
+                        dls ( i )->set_connection ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv );
+
+                        // find out what is missing in each mpdf
+                        xytmp = mpdfs ( i )->_rv();
+                        xytmp.add ( mpdfs ( i )->_rvc() );
+                        // z_i = common_rv-xy
+                        rvzs ( i ) = rv.subt ( xytmp );
+                        //establish connection between extension (z_i|x,y)s and common rv
+                        zdls ( i ) = new datalink_m2e;
+                        zdls ( i )->set_connection ( rvzs ( i ), xytmp, rv ) ;
+                };
+        }
+        //! Rectangular support  each vector of XYZ specifies (begining-end) interval for each dimension. Same number of points, \c dimsize, in each dimension.
+        void set_support ( const Array<vec> &XYZ, const int dimsize ) {
+                set_support ( XYZ, dimsize*ones_i ( XYZ.length() ) );
+        }
+        //! Rectangular support  each vector of XYZ specifies (begining-end) interval for each dimension. \c gridsize specifies number of points is each dimension.
+        void set_support ( const Array<vec> &XYZ, const ivec &gridsize ) {
+                int dim = XYZ.length();  //check with internal dim!!
+                Npoints = prod ( gridsize );
+                eSmp.set_parameters ( Npoints, false );
+                Array<vec> &samples = eSmp._samples();
+                eSmp._w() = ones ( Npoints ) / Npoints; //unifrom size of bins
+                //set samples
+                ivec ind = zeros_i ( dim );    //indeces of dimensions in for cycle;
+                vec smpi ( dim );          // ith sample
+                vec steps = zeros ( dim );        // ith sample
+                // first corner
+                for ( int j = 0; j < dim; j++ ) {
+                        smpi ( j ) = XYZ ( j ) ( 0 ); /* beginning of the interval*/
+                        it_assert ( gridsize ( j ) != 0.0, "Zeros in gridsize!" );
+                        steps ( j ) = ( XYZ ( j ) ( 1 ) - smpi ( j ) ) / gridsize ( j );
+                }
+                // fill samples
+                for ( int i = 0; i < Npoints; i++ ) {
+                        // copy
+                        samples ( i ) = smpi;
+                        // go through all dimensions
+                        for ( int j = 0; j < dim; j++ ) {
+                                if ( ind ( j ) == gridsize ( j ) - 1 ) { //j-th index is full
+                                        ind ( j ) = 0; //shift back
+                                        smpi ( j ) = XYZ ( j ) ( 0 );
+
+                                        if ( i < Npoints - 1 ) {
+                                                ind ( j + 1 ) ++; //increase the next dimension;
+                                                smpi ( j + 1 ) += steps ( j + 1 );
                                                 break;
                                         }
+
+                                } else {
+                                        ind ( j ) ++;
+                                        smpi ( j ) += steps ( j );
+                                        break;
                                 }
                         }
                 }
-                //! set debug file
-                void set_debug_file (const string fname) {
-                        if (DBG) delete dbg_file;
-                        dbg_file = new it_file (fname);
-                        if (dbg_file) DBG = true;
-                }
-                //! Set internal parameters used in approximation
-                void set_method (MERGER_METHOD MTH=DFLT_METHOD, double beta0 = DFLT_beta) {
-                        METHOD = MTH;
-                        beta = beta0;
-                }
-                //! Set support points from a pdf by drawing N samples
-                void set_support (const epdf &overall, int N) {
-                        eSmp.set_statistics (&overall, N);
-                        Npoints = N;
-                }
-
-                //! Destructor
-                virtual ~merger_base() {
-                        for (int i = 0; i < Nsources; i++) {
-                                delete dls (i);
-                                delete zdls (i);
-                        }
-                        if (DBG) delete dbg_file;
-                };
-                //!@}
-
-                //! \name Mathematical operations
-                //!@{
-
-                //!Merge given sources in given points
-                virtual void merge () {
-                        validate();
-
-                        //check if sources overlap:
-                        bool OK = true;
-                        for (int i = 0;i < mpdfs.length(); i++) {
-                                OK &= (rvzs (i)._dsize() == 0); // z_i is empty
-                                OK &= (mpdfs (i)->_rvc()._dsize() == 0); // y_i is empty
-                        }
-
-                        if (OK) {
-                                mat lW = zeros (mpdfs.length(), eSmp._w().length());
-
-                                vec emptyvec (0);
-                                for (int i = 0; i < mpdfs.length(); i++) {
-                                        for (int j = 0; j < eSmp._w().length(); j++) {
-                                                lW (i, j) = mpdfs (i)->evallogcond (eSmp._samples() (j), emptyvec);
-                                        }
-                                }
-
-                                vec w_nn=merge_points (lW);
-                                vec wtmp = exp (w_nn-max(w_nn));
-                                //renormalize
-                                eSmp._w() = wtmp / sum (wtmp);
-                        } else {
-                                it_error ("Sources are not compatible - use merger_mix");
-                        }
-                };
-
-
-                //! Merge log-likelihood values in points using method specified by parameter METHOD
-                vec merge_points (mat &lW);
-
-
-                //! sample from merged density
-//! weight w is a
-                vec mean() const {
-                        const Vec<double> &w = eSmp._w();
-                        const Array<vec> &S = eSmp._samples();
-                        vec tmp = zeros (dim);
-                        for (int i = 0; i < Npoints; i++) {
-                                tmp += w (i) * S (i);
-                        }
-                        return tmp;
-                }
-                mat covariance() const {
-                        const vec &w = eSmp._w();
-                        const Array<vec> &S = eSmp._samples();
-
-                        vec mea = mean();
-
-//                      cout << sum (w) << "," << w*w << endl;
-
-                        mat Tmp = zeros (dim, dim);
-                        for (int i = 0; i < Npoints; i++) {
-                                Tmp += w (i) * outer_product (S (i), S (i));
-                        }
-                        return Tmp -outer_product (mea, mea);
-                }
-                vec variance() const {
-                        const vec &w = eSmp._w();
-                        const Array<vec> &S = eSmp._samples();
-
-                        vec tmp = zeros (dim);
-                        for (int i = 0; i < Nsources; i++) {
-                                tmp += w (i) * pow (S (i), 2);
-                        }
-                        return tmp -pow (mean(), 2);
-                }
-                //!@}
-
-                //! \name Access to attributes
-                //! @{
-
-                //! Access function
-                eEmp& _Smp() {return eSmp;}
-
-                //! load from setting
-                void from_setting (const Setting& set) {
-                        // get support
-                        // find which method to use
-                        string meth_str;
-                        UI::get<string> (meth_str, set, "method", UI::compulsory);
-                        if (!strcmp (meth_str.c_str(), "arithmetic"))
-                                set_method (ARITHMETIC);
-                        else {
-                                if (!strcmp (meth_str.c_str(), "geometric"))
-                                        set_method (GEOMETRIC);
-                                else if (!strcmp (meth_str.c_str(), "lognormal")) {
-                                        set_method (LOGNORMAL);
-                                        set.lookupValue( "beta",beta);
+        }
+        //! set debug file
+        void set_debug_file ( const string fname ) {
+                if ( DBG ) delete dbg_file;
+                dbg_file = new it_file ( fname );
+                if ( dbg_file ) DBG = true;
+        }
+        //! Set internal parameters used in approximation
+        void set_method ( MERGER_METHOD MTH = DFLT_METHOD, double beta0 = DFLT_beta ) {
+                METHOD = MTH;
+                beta = beta0;
+        }
+        //! Set support points from a pdf by drawing N samples
+        void set_support ( const epdf &overall, int N ) {
+                eSmp.set_statistics ( &overall, N );
+                Npoints = N;
+        }
+
+        //! Destructor
+        virtual ~merger_base() {
+                for ( int i = 0; i < Nsources; i++ ) {
+                        delete dls ( i );
+                        delete zdls ( i );
+                }
+                if ( DBG ) delete dbg_file;
+        };
+        //!@}
+
+        //! \name Mathematical operations
+        //!@{
+
+        //!Merge given sources in given points
+        virtual void merge () {
+                validate();
+
+                //check if sources overlap:
+                bool OK = true;
+                for ( int i = 0; i < mpdfs.length(); i++ ) {
+                        OK &= ( rvzs ( i )._dsize() == 0 ); // z_i is empty
+                        OK &= ( mpdfs ( i )->_rvc()._dsize() == 0 ); // y_i is empty
+                }
+
+                if ( OK ) {
+                        mat lW = zeros ( mpdfs.length(), eSmp._w().length() );
+
+                        vec emptyvec ( 0 );
+                        for ( int i = 0; i < mpdfs.length(); i++ ) {
+                                for ( int j = 0; j < eSmp._w().length(); j++ ) {
+                                        lW ( i, j ) = mpdfs ( i )->evallogcond ( eSmp._samples() ( j ), emptyvec );
                                 }
                         }
-                        string dbg_file;
-                        if (UI::get(dbg_file, set, "dbg_file"))
-                                set_debug_file(dbg_file);
-                        //validate() - not used
-                }
-
-                void validate() {
-                        it_assert (eSmp._w().length() > 0, "Empty support, use set_support().");
-                        it_assert (dim == eSmp._samples() (0).length(), "Support points and rv are not compatible!");
-                        it_assert (isnamed(),"mergers must be named");
-                }
-                //!@}
+
+                        vec w_nn = merge_points ( lW );
+                        vec wtmp = exp ( w_nn - max ( w_nn ) );
+                        //renormalize
+                        eSmp._w() = wtmp / sum ( wtmp );
+                } else {
+                        it_error ( "Sources are not compatible - use merger_mix" );
+                }
+        };
+
+
+        //! Merge log-likelihood values in points using method specified by parameter METHOD
+        vec merge_points ( mat &lW );
+
+
+        //! sample from merged density
+//! weight w is a
+        vec mean() const {
+                const Vec<double> &w = eSmp._w();
+                const Array<vec> &S = eSmp._samples();
+                vec tmp = zeros ( dim );
+                for ( int i = 0; i < Npoints; i++ ) {
+                        tmp += w ( i ) * S ( i );
+                }
+                return tmp;
+        }
+        mat covariance() const {
+                const vec &w = eSmp._w();
+                const Array<vec> &S = eSmp._samples();
+
+                vec mea = mean();
+
+//                      cout << sum (w) << "," << w*w << endl;
+
+                mat Tmp = zeros ( dim, dim );
+                for ( int i = 0; i < Npoints; i++ ) {
+                        Tmp += w ( i ) * outer_product ( S ( i ), S ( i ) );
+                }
+                return Tmp - outer_product ( mea, mea );
+        }
+        vec variance() const {
+                const vec &w = eSmp._w();
+                const Array<vec> &S = eSmp._samples();
+
+                vec tmp = zeros ( dim );
+                for ( int i = 0; i < Nsources; i++ ) {
+                        tmp += w ( i ) * pow ( S ( i ), 2 );
+                }
+                return tmp - pow ( mean(), 2 );
+        }
+        //!@}
+
+        //! \name Access to attributes
+        //! @{
+
+        //! Access function
+        eEmp& _Smp() {
+                return eSmp;
+        }
+
+        //! load from setting
+        void from_setting ( const Setting& set ) {
+                // get support
+                // find which method to use
+                string meth_str;
+                UI::get<string> ( meth_str, set, "method", UI::compulsory );
+                if ( !strcmp ( meth_str.c_str(), "arithmetic" ) )
+                        set_method ( ARITHMETIC );
+                else {
+                        if ( !strcmp ( meth_str.c_str(), "geometric" ) )
+                                set_method ( GEOMETRIC );
+                        else if ( !strcmp ( meth_str.c_str(), "lognormal" ) ) {
+                                set_method ( LOGNORMAL );
+                                set.lookupValue ( "beta", beta );
+                        }
+                }
+                string dbg_file;
+                if ( UI::get ( dbg_file, set, "dbg_file" ) )
+                        set_debug_file ( dbg_file );
+                //validate() - not used
+        }
+
+        void validate() {
+                it_assert ( eSmp._w().length() > 0, "Empty support, use set_support()." );
+                it_assert ( dim == eSmp._samples() ( 0 ).length(), "Support points and rv are not compatible!" );
+                it_assert ( isnamed(), "mergers must be named" );
+        }
+        //!@}
 };
-UIREGISTER(merger_base);
-
-class merger_mix : public merger_base
-{
-        protected:
-                //!Internal mixture of EF models
-                MixEF Mix;
-                //!Number of components in a mixture
-                int Ncoms;
-                //! coefficient of resampling [0,1]
-                double effss_coef;
-                //! stop after niter iterations
-                int stop_niter;
-                
-                //! default value for Ncoms
-                static const int DFLT_Ncoms;
-                //! default value for efss_coef;
-                static const double DFLT_effss_coef;
-
-        public:
-                //!\name Constructors
-                //!@{
-                merger_mix () {};
-                merger_mix (const Array<mpdf*> &S,bool own=false) {set_sources (S,own);};
-                //! Set sources and prepare all internal structures
-                void set_sources (const Array<mpdf*> &S, bool own) {
-                        merger_base::set_sources (S,own);
-                        Nsources = S.length();
-                }
-                //! Set internal parameters used in approximation
-                void set_parameters (int Ncoms0 = DFLT_Ncoms, double effss_coef0 = DFLT_effss_coef) {
-                        Ncoms = Ncoms0;
-                        effss_coef = effss_coef0;
-                }
-                //!@}
-
-                //! \name Mathematical operations
-                //!@{
-
-                //!Merge values using mixture approximation
-                void merge ();
-
-                //! sample from the approximating mixture
-                vec sample () const { return Mix.posterior().sample();}
-                //! loglikelihood computed on mixture models
-                double evallog (const vec &dt) const {
-                        vec dtf = ones (dt.length() + 1);
-                        dtf.set_subvector (0, dt);
-                        return Mix.logpred (dtf);
-                }
-                //!@}
-
-                //!\name Access functions
-                //!@{
+UIREGISTER ( merger_base );
+
+class merger_mix : public merger_base {
+protected:
+        //!Internal mixture of EF models
+        MixEF Mix;
+        //!Number of components in a mixture
+        int Ncoms;
+        //! coefficient of resampling [0,1]
+        double effss_coef;
+        //! stop after niter iterations
+        int stop_niter;
+
+        //! default value for Ncoms
+        static const int DFLT_Ncoms;
+        //! default value for efss_coef;
+        static const double DFLT_effss_coef;
+
+public:
+        //!\name Constructors
+        //!@{
+        merger_mix () {};
+        merger_mix ( const Array<mpdf*> &S, bool own = false ) {
+                set_sources ( S, own );
+        };
+        //! Set sources and prepare all internal structures
+        void set_sources ( const Array<mpdf*> &S, bool own ) {
+                merger_base::set_sources ( S, own );
+                Nsources = S.length();
+        }
+        //! Set internal parameters used in approximation
+        void set_parameters ( int Ncoms0 = DFLT_Ncoms, double effss_coef0 = DFLT_effss_coef ) {
+                Ncoms = Ncoms0;
+                effss_coef = effss_coef0;
+        }
+        //!@}
+
+        //! \name Mathematical operations
+        //!@{
+
+        //!Merge values using mixture approximation
+        void merge ();
+
+        //! sample from the approximating mixture
+        vec sample () const {
+                return Mix.posterior().sample();
+        }
+        //! loglikelihood computed on mixture models
+        double evallog ( const vec &dt ) const {
+                vec dtf = ones ( dt.length() + 1 );
+                dtf.set_subvector ( 0, dt );
+                return Mix.logpred ( dtf );
+        }
+        //!@}
+
+        //!\name Access functions
+        //!@{
 //! Access function
-                MixEF& _Mix() {return Mix;}
-                //! Access function
-                emix* proposal() {emix* tmp = Mix.epredictor(); tmp->set_rv (rv); return tmp;}
-                //! from_settings
-                void from_setting(const Setting& set){
-                        merger_base::from_setting(set);
-                        set.lookupValue("ncoms",Ncoms);
-                        set.lookupValue("effss_coef",effss_coef);
-                        set.lookupValue("stop_niter",stop_niter);
-                }
-                
-                //! @}
+        MixEF& _Mix() {
+                return Mix;
+        }
+        //! Access function
+        emix* proposal() {
+                emix* tmp = Mix.epredictor();
+                tmp->set_rv ( rv );
+                return tmp;
+        }
+        //! from_settings
+        void from_setting ( const Setting& set ) {
+                merger_base::from_setting ( set );
+                set.lookupValue ( "ncoms", Ncoms );
+                set.lookupValue ( "effss_coef", effss_coef );
+                set.lookupValue ( "stop_niter", stop_niter );
+        }
+
+        //! @}
 
 };
-UIREGISTER(merger_mix);
+UIREGISTER ( merger_mix );
 
 }