Changeset 270 for bdm/stat

Timestamp:

02/16/09 10:02:08 (15 years ago)

Author:

smidl

Message:

Changes in the very root classes!
* rv and rvc are no longer compulsory,
* samplecond does not return ll
* BM has drv

Location:

bdm/stat

Files:

• emix.cpp (modified) (3 diffs)
• emix.h (modified) (12 diffs)
• libBM.cpp (modified) (11 diffs)
• libBM.h (modified) (3 diffs)
• libDS.cpp (modified) (3 diffs)
• libDS.h (modified) (6 diffs)
• libEF.cpp (modified) (12 diffs)
• libEF.h (modified) (38 diffs)
• libFN.cpp (modified) (1 diff)
• libFN.h (modified) (5 diffs)
• loggers.cpp (modified) (1 diff)
• loggers.h (modified) (1 diff)

bdm/stat/emix.cpp

@@ -5 +5 @@
 void emix::set_parameters ( const vec &w0, const Array<epdf*> &Coms0, bool copy ) {
         w = w0/sum ( w0 );
+        dim = Coms0(0)->dimension();
         int i;
         for ( i=0;i<w.length();i++ ) {
-                it_assert_debug ( rv.equal ( Coms0 ( i )->_rv() ),"RVs do not match!" );
+                it_assert_debug ( dim== ( Coms0 ( i )->dimension() ),"Component sizes do not match!" );
         }
         if ( copy ) {
@@ -35 +36 @@
 
 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(rv);
+        emix* tmp = new emix();
         tmp->set_parameters(w,Cn,false);
         tmp->ownComs();
@@ -44 +47 @@
 
 mratio* emix::condition(const RV &rv) const{
+        it_assert_debug(isnamed(), "rvs are not assigned");
         return new mratio(this,rv);
 };

bdm/stat/emix.h

@@ -47 +47 @@
         //!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 ( rv,nom0->_rv().subt ( rv ) ), dl ( rv,rvc,nom0->_rv() ) {
+        mratio ( const epdf* nom0, const RV &rv, bool copy=false ) :mpdf ( ), dl ( rv,rvc,nom0->_rv() ) {
+                ep->set_rv(_rv());
+                set_rvc(nom0->_rv().subt ( rv ) );
                 if ( copy ) {
 //                      nom = nom0->_copy_();
@@ -62 +64 @@
         double evallogcond ( const vec &val, const vec &cond ) {
                 double tmp;
-                vec nom_val ( rv.count() +rvc.count() );
-                dl.fill_val_cond ( nom_val,val,cond );
+                vec nom_val ( ep->dimension() + dimc );
+                dl.pushup_cond ( nom_val,val,cond );
                 tmp = exp ( nom->evallog ( nom_val ) - den->evallog ( cond ) );
                 it_assert_debug(std::isfinite(tmp),"Infinite value");
@@ -94 +96 @@
 public:
         //!Default constructor
-        emix ( const RV &rv ) : epdf ( rv ) {};
+        emix (  ) : epdf ( ) {};
         //! 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=true );
+        void set_parameters ( const vec &w, const Array<epdf*> &Coms, bool copy=false );
 
         vec sample() const;
         vec mean() const {
-                int i; vec mu = zeros ( rv.count() );
+                int i; vec mu = zeros ( dim);
                 for ( i = 0;i < w.length();i++ ) {mu += w ( i ) * Coms ( i )->mean(); }
                 return mu;
@@ -107 +109 @@
         vec variance() const {
                 //non-central moment
-                vec mom2 = zeros(rv.count());
+                vec mom2 = zeros(dim);
                 for ( int i = 0;i < w.length();i++ ) {mom2 += w ( i ) * pow(Coms ( i )->mean(),2); }
                 //central moment
@@ -165 +167 @@
         //! Data link for each mpdfs
         Array<datalink_m2m*> dls;
+        //! dummy ep
+        epdf dummy;
 public:
         /*!\brief Constructor from list of mFacs,
         */
-        mprod ( Array<mpdf*> mFacs ) : compositepdf ( mFacs ), mpdf ( getrv ( true ),RV() ), epdfs ( n ), dls ( n ) {
-                setrvc ( rv,rvc );
+        mprod ( Array<mpdf*> mFacs ) : compositepdf ( mFacs ), mpdf (), epdfs ( n ), dls ( n ) {
+                ep=&dummy;
+                RV rv=getrv(true);
+                set_rv(rv);dummy.set_parameters(rv._dsize());
+                setrvc ( ep->_rv(),rvc );
                 // rv and rvc established = > we can link them with mpdfs
                 for ( int i = 0;i < n;i++ ) {
-                        dls ( i ) = new datalink_m2m ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), rv, rvc );
+                        dls ( i ) = new datalink_m2m ( mpdfs ( i )->_rv(), mpdfs ( i )->_rvc(), _rv(), _rvc() );
                 }
 
@@ -182 +189 @@
         double evallogcond ( const vec &val, const vec &cond ) {
                 int i;
-                double res = 1.0;
+                double res = 0.0;
                 for ( i = n - 1;i >= 0;i-- ) {
                         /*                      if ( mpdfs(i)->_rvc().count() >0) {
@@ -188 +195 @@
                                                 }
                                                 // add logarithms
-                                                res += epdfs ( i )->evallog ( dls ( i )->get_val ( val ) );*/
-                        res *= mpdfs ( i )->evallogcond (
-                                   dls ( i )->get_val ( val ),
+                                                res += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );*/
+                        res += mpdfs ( i )->evallogcond (
+                                   dls ( i )->pushdown ( val ),
                                    dls ( i )->get_cond ( val, cond )
                                );
@@ -196 +203 @@
                 return res;
         }
-        vec samplecond ( const vec &cond, double &ll ) {
+        //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 ( rv.count() );
+                vec smp= std::numeric_limits<double>::infinity() * ones ( ep->dimension() );
                 vec smpi;
-                ll = 0;
                 // Hard assumption here!!! We are going backwards, to assure that samples that are needed from smp are already generated!
                 for ( int i = ( n - 1 );i >= 0;i-- ) {
-                        if ( mpdfs ( i )->_rvc().count() ) {
+                        if ( mpdfs ( i )->dimensionc() ) {
                                 mpdfs ( i )->condition ( dls ( i )->get_cond ( smp ,cond ) ); // smp is val here!!
                         }
                         smpi = epdfs ( i )->sample();
                         // copy contribution of this pdf into smp
-                        dls ( i )->fill_val ( smp, smpi );
+                        dls ( i )->pushup ( smp, smpi );
                         // add ith likelihood contribution
-                        ll+=epdfs ( i )->evallog ( smpi );
                 }
                 return smp;
         }
-        mat samplecond ( const vec &cond, vec &ll, int N ) {
-                mat Smp ( rv.count(),N );
-                for ( int i=0;i<N;i++ ) {Smp.set_col ( i,samplecond ( cond,ll ( i ) ) );}
+        mat samplecond ( const vec &cond,  int N ) {
+                mat Smp ( dimension(),N );
+                for ( int i=0;i<N;i++ ) {Smp.set_col ( i,samplecond ( cond ) );}
                 return Smp;
         }
@@ -229 +235 @@
         Array<const epdf*> epdfs;
         //! Array of indeces
-        Array<datalink_e2e*> dls;
-public:
-        eprod ( const Array<const epdf*> epdfs0 ) : epdf ( RV() ),epdfs ( epdfs0 ),dls ( epdfs.length() ) {
+        Array<datalink*> dls;
+public:
+        eprod ( const Array<const epdf*> epdfs0 ) : epdf (  ),epdfs ( epdfs0 ),dls ( epdfs.length() ) {
                 bool independent=true;
                 for ( int i=0;i<epdfs.length();i++ ) {
@@ -238 +244 @@
                 }
                 for ( int i=0;i<epdfs.length();i++ ) {
-                        dls ( i ) = new datalink_e2e ( epdfs ( i )->_rv() , rv );
+                        dls ( i ) = new datalink ( epdfs ( i )->_rv() , rv );
                 }
         }
 
         vec mean() const {
-                vec tmp ( rv.count() );
+                vec tmp ( dim );
                 for ( int i=0;i<epdfs.length();i++ ) {
                         vec pom = epdfs ( i )->mean();
-                        dls ( i )->fill_val ( tmp, pom );
+                        dls ( i )->pushup ( tmp, pom );
                 }
                 return tmp;
         }
         vec variance() const {
-                vec tmp ( rv.count() ); //second moment
+                vec tmp ( dim ); //second moment
                 for ( int i=0;i<epdfs.length();i++ ) {
                         vec pom = epdfs ( i )->mean();
-                        dls ( i )->fill_val ( tmp, pow(pom,2) );
+                        dls ( i )->pushup ( tmp, pow(pom,2) );
                 }
                 return tmp-pow(mean(),2);
         }
         vec sample() const {
-                vec tmp ( rv.count() );
+                vec tmp ( dim );
                 for ( int i=0;i<epdfs.length();i++ ) {
                         vec pom = epdfs ( i )->sample();
-                        dls ( i )->fill_val ( tmp, pom );
+                        dls ( i )->pushup ( tmp, pom );
                 }
                 return tmp;
@@ -269 +275 @@
                 double tmp=0;
                 for ( int i=0;i<epdfs.length();i++ ) {
-                        tmp+=epdfs ( i )->evallog ( dls ( i )->get_val ( val ) );
+                        tmp+=epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
                 }
                 it_assert_debug(std::isfinite(tmp),"Infinite");
@@ -293 +299 @@
 public:
         //!Default constructor
-        mmix ( RV &rv, RV &rvc ) : mpdf ( rv, rvc ), Epdf ( rv ) {ep = &Epdf;};
+        mmix ( ) : mpdf ( ), Epdf () {ep = &Epdf;};
         //! Set weights \c w and components \c R
         void set_parameters ( const vec &w, const Array<mpdf*> &Coms ) {

bdm/stat/libBM.cpp

@@ -5 +5 @@
 //! Space of basic BDM structures
 namespace bdm {
-
-using std::cout;
-
-static int RVcounter=0;
+        
+const int RV_BUFFER_STEP=1;
+RVmap RV_MAP;
+Array<string> RV_NAMES(RV_BUFFER_STEP);
+ivec RV_SIZES(RV_BUFFER_STEP);
 
 RV RV0=RV();
 
-void RV::init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times ) {
+int RV::init ( const string &name, int size ) {
         //Refer
-        int i;
-        len = in_ids.length();
-        //PRUDENT_MODE
-        // All vetors should be of same length
-        if ( ( len != in_names.length() ) || \
-                ( len != in_sizes.length() ) || \
-                ( len != in_times.length() ) ) {
-                it_error ( "RV::RV inconsistent length of input vectors." );
-        }
-
-        ids = in_ids;
-        names = in_names;
-        sizes = in_sizes;
+        int id;
+        RVmap::const_iterator iter = RV_MAP.find ( name );
+        if ( iter == RV_MAP.end() ) {
+                id=RV_MAP.size()+1;
+                RV_MAP.insert ( make_pair ( name,id ) ); //add new rv
+                if (id>=RV_NAMES.length()){
+                        RV_NAMES.set_length(id+RV_BUFFER_STEP,true);
+                        RV_SIZES.set_length(id+RV_BUFFER_STEP,true);
+                }
+                RV_NAMES(id)=name;
+                RV_SIZES(id)=size;
+        }
+        else {
+                id = iter->second;
+                it_warning ("RV of name " + name + "already exists" );
+        }
+        return id;
+};
+
+int RV::countsize() const{ int tmp=0; for(int i=0;i<len;i++){tmp+=RV_SIZES(ids(i));} return tmp;}
+
+void RV::init ( Array<std::string> in_names, ivec in_sizes,ivec in_times ) {
+        len = in_names.length();
+        it_assert_debug ( in_names.length() ==in_times.length(), "check \"times\" " );
+        it_assert_debug ( in_names.length() ==in_sizes.length(), "check \"sizes\" " );
+        
+        times.set_length ( len );
+        ids.set_length ( len );
+        int id;
+        for ( int i=0; i<len; i++ ) {
+                id = init ( in_names ( i ), in_sizes ( i ) );
+                ids ( i ) = id;
+        }
         times = in_times;
-        tsize = 0;
-        for ( i = 0;i < len;i++ ) {tsize += sizes ( i );}
-};
-
-RV::RV ( Array<std::string> in_names, ivec in_sizes, ivec in_times ) {
-        int len = in_names.length();
-        init ( linspace ( RVcounter+1, RVcounter+len ), in_names, in_sizes, in_times );
-        RVcounter+=len;
-}
-
-RV::RV ( Array<std::string> in_names, ivec in_sizes ) {
-        int len = in_names.length();
-        init ( linspace ( RVcounter+1, RVcounter+len ), in_names, in_sizes, zeros_i ( in_names.length() ) );
-        RVcounter+=len;
-}
-
-RV::RV ( Array<std::string> in_names ) {
-        int len = in_names.length();
-        init ( linspace ( RVcounter+1, RVcounter+len ), in_names, ones_i ( in_names.length() ), zeros_i ( in_names.length() ) );
-        RVcounter+=len;
-}
-
-RV::RV() : tsize ( 0 ), len ( 0 ), ids ( 0 ),  sizes ( 0 ), times ( 0 ),names ( 0 ) {};
-
-RV::RV(string name, int id, int sz, int tm){
-        if (id>RVcounter) {RVcounter=id;};
-        Array<string> A(1); A(0)=name;
-        init(vec_1(id),A,vec_1(sz),vec_1(tm));
+        dsize = countsize();
+}
+
+RV::RV ( string name, int sz, int tm ) {
+        Array<string> A ( 1 ); A ( 0 ) =name;
+        init (A,vec_1 ( sz ),vec_1 ( tm ) );
 }
 
@@ -64 +63 @@
                 ivec index = itpp::find ( ind==-1 );
 
-
                 if ( index.length() < rv2.len ) { //conflict
                         ids = concat ( ids, rv2.ids ( index ) );
-                        sizes = concat ( sizes, rv2.sizes ( index ) );
                         times = concat ( times, rv2.times ( index ) );
-                        names = concat ( names, rv2.names ( to_Arr ( index ) ) );
                 }
                 else {
                         ids = concat ( ids, rv2.ids );
-                        sizes = concat ( sizes, rv2.sizes );
                         times = concat ( times, rv2.times );
-                        names = concat ( names, rv2.names );
-                }
-                tsize = sum ( sizes );
+                }
                 len = ids.length();
-                return ( index.length() ==rv2.len );
+                dsize = countsize();
+                return ( index.length() ==rv2.len ); //conflict or not
         }
         else { //rv2 is empty
-                return true;
-        }
-//      return *this;
+                return true; // no conflict
+        }
 };
 
-// RV::RV ( ivec in_ids ) {
-//
-//      len = in_ids.length();
-//      Array<std::string> A ( len );
-//      std::string rvstr = "rv";
-//
-//      for ( int i = 0; i < len;i++ ) {
-//              A ( i ) = rvstr + to_str ( i );
-//      }
-//
-//      init ( in_ids, A, ones_i ( len ), zeros_i ( len ) );
-// }
-
-RV RV::subselect (const ivec &ind ) const {
+RV RV::subselect ( const ivec &ind ) const {
         RV ret;
-        ret.init ( ids ( ind ), names ( to_Arr ( ind ) ), sizes ( ind ), times ( ind ) );
+        ret.ids = ids(ind);
+        ret.times= times(ind);
+        ret.len = ind.length();
+        ret.dsize=ret.countsize();
         return ret;
 }
@@ -108 +91 @@
 void RV::t ( int delta ) { times += delta;}
 
-RV RV::operator() (const ivec &ind ) const {
-        RV ret;
-        if ( ind.length() >0 ) {
-                ret.init ( ids ( ind ), names ( to_Arr ( ind ) ), sizes ( ind ), times ( ind ) );
-        }
-        return ret;
-}
-
 bool RV::equal ( const RV &rv2 ) const {
-        return ( ids == rv2.ids ) && ( times == rv2.times ) && ( sizes == rv2.sizes );
+        return ( ids == rv2.ids ) && ( times == rv2.times );
 }
 
 mat epdf::sample_m ( int N ) const {
-        mat X = zeros ( rv.count(), N );
+        mat X = zeros ( dim, N );
         for ( int i = 0;i < N;i++ ) X.set_col ( i, this->sample() );
         return X;
@@ -128 +103 @@
 
 std::ostream &operator<< ( std::ostream &os, const RV &rv ) {
-
+        int id;
         for ( int i = 0; i < rv.len ;i++ ) {
-                os << rv.ids ( i ) << "(" << rv.sizes ( i ) << ")" <<  // id(size)=
-                "=" << rv.names ( i )  << "_{"  << rv.times ( i ) << "}; "; //name_{time}
+                id=rv.ids ( i );
+                os << id << "(" << RV_SIZES ( id ) << ")" <<  // id(size)=
+                "=" << RV_NAMES ( id )  << "_{"  << rv.times ( i ) << "}; "; //name_{time}
         }
         return os;
@@ -137 +113 @@
 
 str RV::tostr() const {
-        ivec idlist ( tsize );
-        ivec tmlist ( tsize );
+        ivec idlist ( dsize );
+        ivec tmlist ( dsize );
         int i;
         int pos = 0;
         for ( i = 0;i < len;i++ ) {
-                idlist.set_subvector ( pos, pos + sizes ( i ) - 1, ids ( i ) );
-                tmlist.set_subvector ( pos, pos + sizes ( i ) - 1, times ( i ) );
-                pos += sizes ( i );
+                idlist.set_subvector ( pos, pos + size ( ids(i) ) - 1, ids ( i ) );
+                tmlist.set_subvector ( pos, pos + size ( ids(i) ) - 1, times ( i ) );
+                pos += size ( ids(i) );
         }
         return str ( idlist, tmlist );
 }
 
-ivec RV::dataind (const RV &rv2 ) const {
+ivec RV::dataind ( const RV &rv2 ) const {
         ivec res ( 0 );
-        if ( rv2.count()>0 ) {
+        if ( rv2._dsize() >0 ) {
                 str str2 = rv2.tostr();
                 ivec part;
@@ -160 +136 @@
                 }
         }
-        it_assert_debug(res.length()==tsize,"this rv is not fully present in crv!");
+        it_assert_debug ( res.length() ==dsize,"this rv is not fully present in crv!" );
         return res;
 
 }
 
-void RV::dataind (const RV &rv2, ivec &selfi, ivec &rv2i) const {
+void RV::dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const {
         //clean results
-        selfi.set_size(0);
-        rv2i.set_size(0);
-        
+        selfi.set_size ( 0 );
+        rv2i.set_size ( 0 );
+
         // just in case any rv is empty
-        if ((len==0)||(rv2.length()==0)){return;} 
-        
+        if ( ( len==0 ) || ( rv2.length() ==0 ) ) {return;}
+
         //find comon rv
-        ivec cids=itpp::find(this->findself(rv2)>=0);
-        
-        // index of 
-        if ( cids.length()>0 ) {
+        ivec cids=itpp::find ( this->findself ( rv2 ) >=0 );
+
+        // index of
+        if ( cids.length() >0 ) {
                 str str1 = tostr();
-                str str2 = rv2.tostr(); 
-                
+                str str2 = rv2.tostr();
+
                 ivec part1;
                 ivec part2;
@@ -186 +162 @@
                 // find common rv in strs
                 for ( j=0; j < cids.length();j++ ) {
-                        i = cids(j);
+                        i = cids ( j );
                         part1 = itpp::find ( ( str1.ids == ids ( i ) ) & ( str1.times == times ( i ) ) );
                         part2 = itpp::find ( ( str2.ids == ids ( i ) ) & ( str2.times == times ( i ) ) );
@@ -193 +169 @@
                 }
         }
-        it_assert_debug(selfi.length() == rv2i.length(),"this should not happen!");
+        it_assert_debug ( selfi.length() == rv2i.length(),"this should not happen!" );
 }
 
@@ -199 +175 @@
         ivec res = this->findself ( rv2 ); // nonzeros
         ivec valid;
-        if (tsize>0) {valid= itpp::find ( res == -1 );} //-1 => value not found => it remains
+        if ( dsize>0 ) {valid= itpp::find ( res == -1 );} //-1 => value not found => it remains
         return ( *this ) ( valid ); //keep those that were not found in rv2
 }
@@ -223 +199 @@
 }
 
-void RV::newids(){ids=linspace ( RVcounter+1, RVcounter+len ),RVcounter+=len;}
-
-RV compositepdf::getrv(bool checkoverlap){
+RV compositepdf::getrv ( bool checkoverlap ) {
         RV rv; //empty rv
         bool rvaddok;
-        for (int i = 0;i < n;i++ ) {
+        for ( int i = 0;i < n;i++ ) {
                 rvaddok=rv.add ( mpdfs ( i )->_rv() ); //add rv to common rvs.
-                        // If rvaddok==false, mpdfs overlap => assert error.
-                it_assert_debug(rvaddok||(!checkoverlap),"mprod::mprod() input mpdfs overlap in rv!");
+                // If rvaddok==false, mpdfs overlap => assert error.
+                it_assert_debug ( rvaddok|| ( !checkoverlap ),"mprod::mprod() input mpdfs overlap in rv!" );
         };
         return rv;
 }
 
-void compositepdf::setrvc(const RV &rv, RV &rvc){
-        for (int i = 0;i < n;i++ ) {
+void compositepdf::setrvc ( const RV &rv, RV &rvc ) {
+        for ( int i = 0;i < n;i++ ) {
                 RV rvx = mpdfs ( i )->_rvc().subt ( rv );
                 rvc.add ( rvx ); //add rv to common rvc
@@ -243 +217 @@
 }
 
-void BM::bayesB(const mat &Data){
-        for(int t=0;t<Data.cols();t++){bayes(Data.get_col(t));}
-}
-}
+void BM::bayesB ( const mat &Data ) {
+        for ( int t=0;t<Data.cols();t++ ) {bayes ( Data.get_col ( t ) );}
+}
+}

bdm/stat/libBM.h

@@ -14 +14 @@
 \defgroup core Core BDM classes
 @{
+\defgroup const Constructors
+\defgroup math Mathematical operations
 */
 #ifndef BM_H
@@ -20 +22 @@
 
 #include "../itpp_ext.h"
-//#include <std>
+#include <map>
 
 namespace bdm {
-        using namespace itpp;
-        using std::string;
+using namespace itpp;
+using namespace std;
 
 //! Root class of BDM objects
-        class bdmroot {
-                virtual void print() {}
+class bdmroot {
+public:
+        //! make sure this is a virtual object
+        virtual ~bdmroot() {}
+};
+
+typedef std::map<string, int> RVmap;
+extern ivec RV_SIZES;
+extern Array<string> RV_NAMES;
+
+//! Structure of RV (used internally), i.e. expanded RVs
+class str {
+public:
+        //! vector id ids (non-unique!)
+        ivec ids;
+        //! vector of times
+        ivec times;
+        //!Default constructor
+        str ( ivec ids0, ivec times0 ) :ids ( ids0 ),times ( times0 ) {
+                it_assert_debug ( times0.length() ==ids0.length(),"Incompatible input" );
         };
-
-//! Structure of RV (used internally), i.e. expanded RVs
-        class str {
-        public:
-                //! vector id ids (non-unique!)
-                ivec ids;
-                //! vector of times
-                ivec times;
-                //!Default constructor
-                str ( ivec ids0, ivec times0 ) :ids ( ids0 ),times ( times0 ) {
-                        it_assert_debug ( times0.length() ==ids0.length(),"Incompatible input" );
-                };
+};
+
+/*!
+* \brief Class representing variables, most often random variables
+
+The purpose of this class is to decribe a vector of data. Such description is used for connecting various vectors between each other, see class datalink.
+
+The class is implemented using global variables to assure uniqueness of description:
+
+ In is a vector
+\dot
+digraph datalink {
+rankdir=LR;
+subgraph cluster0 {
+node [shape=record];
+label = "RV_MAP \n std::map<string,int>";
+map [label="{{\"a\"| \"b\" | \"c\"} | {<3> 3 |<1> 1|<2> 2}}"];
+color = "white"
+}
+subgraph cluster1{
+node [shape=record];
+label = "RV_NAMES";
+names [label="{<1> \"b\" | <2> \"c\" | <3>\"a\" }"];
+color = "white"
+}
+subgraph cluster2{
+node [shape=record];
+label = "RV_SIZES";
+labelloc = b;
+sizes [label="{<1>1 |<2> 4 |<3> 1}"];
+color = "white"
+}
+map:1 -> names:1;
+map:1 -> sizes:1;
+map:3 -> names:3;
+map:3 -> sizes:3;
+}
+\enddot
+*/
+
+class RV :public bdmroot {
+protected:
+        //! size of the data vector
+        int dsize;
+        //! number of individual rvs
+        int len;
+        //! Vector of unique IDs
+        ivec ids;
+        //! Vector of shifts from current time
+        ivec times;
+
+private:
+        //! auxiliary function used in constructor
+        void init ( Array<std::string> in_names, ivec in_sizes, ivec in_times );
+        int init ( const  string &name, int size );
+public:
+        //! \name Constructors 
+        //!@{
+        
+        //! Full constructor 
+        RV ( Array<std::string> in_names, ivec in_sizes, ivec in_times ) {init ( in_names,in_sizes,in_times );};
+        //! Constructor with times=0 
+        RV ( Array<std::string> in_names, ivec in_sizes ) {init ( in_names,in_sizes,zeros_i ( in_names.length() ) );};
+        //! Constructor with sizes=1, times=0 
+        RV ( Array<std::string> in_names ) {init ( in_names,ones_i ( in_names.length() ),zeros_i ( in_names.length() ) );}
+        //! Constructor of empty RV 
+        RV () :dsize ( 0 ),len ( 0 ),ids ( 0 ),times ( 0 ) {};
+        //! Constructor of a single RV with given id
+        RV ( string name, int sz, int tm=0 );
+        //!@}
+        
+        //! \name Access functions
+        //!@{
+        
+        //! Printing output e.g. for debugging.
+        friend std::ostream &operator<< ( std::ostream &os, const RV &rv );
+        int _dsize() const {return dsize;} ;
+        //! Recount size of the corresponding data vector
+        int countsize() const;
+        int length() const {return len;} ;
+        int id ( int at ) const{return ids ( at );};
+        int size ( int at ) const {return RV_SIZES ( at );};
+        int time ( int at ) const{return times ( at );};
+        std::string name ( int at ) const {return RV_NAMES ( at );};
+        void set_time ( int at, int time0 ) {times ( at ) =time0;};
+        //!@}
+        
+        //TODO why not inline and later??
+
+        //! \name Algebra on Random Variables
+        //!@{
+        
+        //! Find indices of self in another rv, \return ivec of the same size as self.
+        ivec findself ( const RV &rv2 ) const;
+        //! Compare if \c rv2 is identical to this \c RV
+        bool equal ( const RV &rv2 ) const;
+        //! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict
+        bool add ( const RV &rv2 );
+        //! Subtract  another variable from the current one
+        RV subt ( const RV &rv2 ) const;
+        //! Select only variables at indeces ind
+        RV subselect ( const ivec &ind ) const;
+        //! Select only variables at indeces ind
+        RV operator() ( const ivec &ind ) const {return subselect ( ind );};
+        //! Shift \c time shifted by delta.
+        void t ( int delta );
+        //!@}
+        
+        //!\name Relation to vectors 
+        //!@{
+        
+        //! generate \c str from rv, by expanding sizes
+        str tostr() const;
+        //! when this rv is a part of bigger rv, this function returns indeces of self in the data vector of the bigger crv.
+        //! Then, data can be copied via: data_of_this = cdata(ind);
+        ivec dataind ( const RV &crv ) const;
+        //! generate mutual indeces when copying data betwenn self and crv.
+        //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i)
+        void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const;
+        //! Minimum time-offset
+        int mint () const {return min ( times );};
+        //!@}
+        
+};
+
+
+//! Concat two random variables
+RV concat ( const RV &rv1, const RV &rv2 );
+
+//!Default empty RV that can be used as default argument
+extern RV RV0;
+
+//! Class representing function \f$f(x)\f$ of variable \f$x\f$ represented by \c rv
+
+class fnc :public bdmroot {
+protected:
+        //! Length of the output vector
+        int dimy;
+public:
+        //!default constructor
+        fnc ( ) {};
+        //! function evaluates numerical value of \f$f(x)\f$ at \f$x=\f$ \c cond
+        virtual vec eval ( const vec &cond ) {
+                return vec ( 0 );
         };
 
-        /*!
-        * \brief Class representing variables, most often random variables
-
-        * More?...
-        */
-
-        class RV :public bdmroot {
-        protected:
-                //! size = sum of sizes
-                int tsize;
-                //! len = number of individual rvs
-                int len;
-                //! Vector of unique IDs
-                ivec ids;
-                //! Vector of sizes
-                ivec sizes;
-                //! Vector of shifts from current time
-                ivec times;
-                //! Array of names
-                Array<std::string> names;
-
-        private:
-                //! auxiliary function used in constructor
-                void init ( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times );
-        public:
-                //! Full constructor
-                RV ( Array<std::string> in_names, ivec in_sizes, ivec in_times );
-                //! Constructor with times=0
-                RV ( Array<std::string> in_names, ivec in_sizes );
-                //! Constructor with sizes=1, times=0
-                RV ( Array<std::string> in_names );
-                //! Constructor of empty RV
-                RV ();
-                //! Constructor of a single RV with given id
-                RV (string name, int id, int sz=1, int tm=0);
-
-                //! Printing output e.g. for debugging.
-                friend std::ostream &operator<< ( std::ostream &os, const RV &rv );
-
-                //! Return number of scalars in the RV.
-                int count() const {return tsize;} ;
-                //! Return length (number of entries) of the RV.
-                int length() const {return len;} ;
-
-                //TODO why not inline and later??
-
-                //! Find indexes of self in another rv, \return ivec of the same size as self.
-                ivec findself ( const RV &rv2 ) const;
-                //! Compare if \c rv2 is identical to this \c RV
-                bool equal ( const RV &rv2 ) const;
-                //! Add (concat) another variable to the current one, \return true if all rv2 were added, false if rv2 is in conflict
-                bool add ( const RV &rv2 );
-                //! Subtract  another variable from the current one
-                RV subt ( const RV &rv2 ) const;
-                //! Select only variables at indeces ind
-                RV subselect ( const ivec &ind ) const;
-                //! Select only variables at indeces ind
-                RV operator() ( const ivec &ind ) const;
-                //! Shift \c time shifted by delta.
-                void t ( int delta );
-                //! generate \c str from rv, by expanding sizes
-                str tostr() const;
-                //! when this rv is a part of bigger rv, this function returns indeces of self in the data vector of the bigger crv.
-                //! Then, data can be copied via: data_of_this = cdata(ind);
-                ivec dataind ( const RV &crv ) const;
-                //! generate mutual indeces when copying data betwenn self and crv.
-                //! Data are copied via: data_of_this(selfi) = data_of_rv2(rv2i)
-                void dataind ( const RV &rv2, ivec &selfi, ivec &rv2i ) const;
-                //! Minimum time-offset
-                int mint () const {return min ( times );};
-
-                //!access function
-                Array<std::string>& _names() {return names;};
-
-                //!access function
-                int id ( int at ) {return ids ( at );};
-                //!access function
-                int size ( int at ) {return sizes ( at );};
-                //!access function
-                int time ( int at ) {return times ( at );};
-                //!access function
-                std::string name ( int at ) {return names ( at );};
-
-                //!access function
-                void set_id ( int at, int id0 ) {ids ( at ) =id0;};
-                //!access function
-                void set_size ( int at, int size0 ) {sizes ( at ) =size0; tsize=sum ( sizes );};
-                //!access function
-                void set_time ( int at, int time0 ) {times ( at ) =time0;};
-
-                //!Assign unused ids to this rv
-                void newids();
-        };
-
-//! Concat two random variables
-        RV concat ( const RV &rv1, const RV &rv2 );
-
-//!Default empty RV that can be used as default argument
-        extern RV RV0;
-
-//! Class representing function \f$f(x)\f$ of variable \f$x\f$ represented by \c rv
-
-        class fnc :public bdmroot {
-        protected:
-                //! Length of the output vector
-                int dimy;
-        public:
-                //!default constructor
-                fnc ( int dy ) :dimy ( dy ) {};
-                //! function evaluates numerical value of \f$f(x)\f$ at \f$x=\f$ \c cond
-                virtual vec eval ( const vec &cond ) {
-                        return vec ( 0 );
-                };
-
-                //! function substitutes given value into an appropriate position
-                virtual void condition ( const vec &val ) {};
-
-                //! access function
-                int _dimy() const{return dimy;}
-
-                //! Destructor for future use;
-                virtual ~fnc() {};
-        };
-
-        class mpdf;
+        //! function substitutes given value into an appropriate position
+        virtual void condition ( const vec &val ) {};
+
+        //! access function
+        int _dimy() const{return dimy;}
+};
+
+class mpdf;
 
 //! Probability density function with numerical statistics, e.g. posterior density.
 
-        class epdf :public bdmroot {
-        protected:
-                //! Identified of the random variable
-                RV rv;
-        public:
-                //!default constructor
-                epdf() :rv ( ) {};
-
-                //!default constructor
-                epdf ( const RV &rv0 ) :rv ( rv0 ) {};
-
-//      //! Returns the required moment of the epdf
-//      virtual vec moment ( const int order = 1 );
-
-                //! Returns a sample, \f$x\f$ from density \f$epdf(rv)\f$
-                virtual vec sample () const =0;
-                //! Returns N samples from density \f$epdf(rv)\f$
-                virtual mat sample_m ( int N ) const;
-
-                //! Compute log-probability of argument \c val
-                virtual double evallog ( const vec &val ) const =0;
-
-                //! Compute log-probability of multiple values argument \c val
-                virtual vec evallog_m ( const mat &Val ) const {
-                        vec x ( Val.cols() );
-                        for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog ( Val.get_col ( i ) ) ;}
-                        return x;
-                }
-                //! Return conditional density on the given RV, the remaining rvs will be in conditioning
-                virtual mpdf* condition ( const RV &rv ) const  {it_warning ( "Not implemented" ); return NULL;}
-                //! Return marginal density on the given RV, the remainig rvs are intergrated out
-                virtual epdf* marginal ( const RV &rv ) const {it_warning ( "Not implemented" ); return NULL;}
-
-                //! return expected value
-                virtual vec mean() const =0;
-
-                //! return expected variance (not covariance!)
-                virtual vec variance() const = 0;
-
-                //! Destructor for future use;
-                virtual ~epdf() {};
-                //! access function, possibly dangerous!
-                const RV& _rv() const {return rv;}
-                //! modifier function - useful when copying epdfs
-                void _renewrv ( const RV &in_rv ) {rv=in_rv;}
-                //!
-        };
+class epdf :public bdmroot {
+protected:
+        //! dimension of the random variable
+        int dim;
+        //! Description of the random variable
+        RV rv;
+
+public:
+        /*! \name Constructors
+         Construction of each epdf should support two types of constructors: 
+        \li empty constructor, 
+        \li copy constructor,
+        
+        The following constructors should be supported for convenience:
+        \li constructor followed by calling \c set_parameters() 
+        \li constructor accepting random variables calling \c set_rv()
+        
+         All internal data structures are constructed as empty. Their values (including sizes) will be set by method \c set_parameters(). This way references can be initialized in constructors.
+        @{*/
+        epdf() :dim(0),rv ( ) {};
+        epdf(const epdf &e) :dim(e.dim),rv (e.rv) {};
+        epdf(const RV &rv0) {set_rv(rv0);};
+        void set_parameters(int dim0){dim=dim0;}
+        //!@}
+        
+        //! \name Matematical Operations
+        //!@{
+        
+        //! Returns a sample, \f$ x \f$ from density \f$ f_x()\f$
+        virtual vec sample () const {it_error("not implemneted");return vec(0);};
+        //! Returns N samples, \f$ [x_1 , x_2 , \ldots \ \f$  from density \f$ f_x(rv)\f$
+        virtual mat sample_m ( int N ) const;
+        //! Compute log-probability of argument \c val
+        virtual double evallog ( const vec &val ) const {it_error("not implemneted");return 0.0;};
+        //! Compute log-probability of multiple values argument \c val
+        virtual vec evallog_m ( const mat &Val ) const {
+                vec x ( Val.cols() );
+                for ( int i=0;i<Val.cols();i++ ) {x ( i ) =evallog ( Val.get_col ( i ) ) ;}
+                return x;
+        }
+        //! Return conditional density on the given RV, the remaining rvs will be in conditioning
+        virtual mpdf* condition ( const RV &rv ) const  {it_warning ( "Not implemented" ); return NULL;}
+        //! Return marginal density on the given RV, the remainig rvs are intergrated out
+        virtual epdf* marginal ( const RV &rv ) const {it_warning ( "Not implemented" ); return NULL;}
+        //! return expected value
+        virtual vec mean() const {it_error("not implemneted");return vec(0);};
+        //! return expected variance (not covariance!)
+        virtual vec variance() const {it_error("not implemneted");return vec(0);};
+        //!@}
+        
+        //! \name Connection to other classes
+        //! Description of the random quantity via attribute \c rv is optional. 
+        //! For operations such as sampling \c rv does not need to be set. However, for \c marginalization 
+        //! and \c conditioning \c rv has to be set. NB: 
+        //! @{
+        
+        //!Name its rv
+        void set_rv ( const RV &rv0 ) {rv = rv0; }//it_assert_debug(isnamed(),""); };
+        //! True if rv is assigned 
+        bool isnamed() const {bool b=( dim==rv._dsize() );return b;}
+        //! Return name (fails when isnamed is false)
+        const RV& _rv() const {it_assert_debug ( isnamed(),"" ); return rv;}
+        //!@}
+        
+        //! \name Access to attributes
+        //! @{
+        
+        //! Size of the random variable
+        int dimension() const {return dim;}
+        //!@}
+        
+};
 
 
@@ -224 +284 @@
 //TODO Samplecond can be generalized
 
-        class mpdf : public bdmroot {
-        protected:
-                //! modeled random variable
-                RV rv;
-                //! random variable in condition
-                RV rvc;
-                //! pointer to internal epdf
-                epdf* ep;
-        public:
-
-                //! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv 
-                virtual vec samplecond ( const vec &cond) {
-                        this->condition ( cond );
-                        vec temp= ep->sample();
-                        return temp;
-                };
-                //! Returns \param N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv \param ll is a return value of log-likelihood of the sample.
-                virtual mat samplecond_m ( const vec &cond, vec &ll, int N ) {
-                        this->condition ( cond );
-                        mat temp ( rv.count(),N ); vec smp ( rv.count() );
-                        for ( int i=0;i<N;i++ ) {smp=ep->sample() ;temp.set_col ( i, smp );ll ( i ) =ep->evallog ( smp );}
-                        return temp;
-                };
-                //! Update \c ep so that it represents this mpdf conditioned on \c rvc = cond
-                virtual void condition ( const vec &cond ) {it_error ( "Not implemented" );};
-
-                //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
-                virtual double evallogcond ( const vec &dt, const vec &cond ) {
-                        double tmp; this->condition ( cond );tmp = ep->evallog ( dt );          it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); return tmp;
-                };
-
-                //! Matrix version of evallogcond
-                virtual vec evallogcond_m ( const mat &Dt, const vec &cond ) {this->condition ( cond );return ep->evallog_m ( Dt );};
-
-                //! Destructor for future use;
-                virtual ~mpdf() {};
-
-                //! Default constructor
-                mpdf ( const RV &rv0, const RV &rvc0 ) :rv ( rv0 ),rvc ( rvc0 ) {};
-                //! access function
-                RV _rvc() const {return rvc;}
-                //! access function
-                RV _rv() const {return rv;}
-                //!access function
-                epdf& _epdf() {return *ep;}
-                //!access function
-                epdf* _e() {return ep;}
+class mpdf : public bdmroot {
+protected:
+        //!dimension of the condition
+        int dimc;
+        //! random variable in condition
+        RV rvc;
+        //! pointer to internal epdf
+        epdf* ep;
+public:
+        //! \name Constructors
+        //! @{
+        
+        mpdf ( ) :dimc(0),rvc ( ) {};
+        //! copy constructor does not set pointer \c ep - has to be done in offsprings!
+        mpdf (const mpdf &m ) :dimc(m.dimc),rvc (m.rvc ) {};
+        //!@}
+
+        //! \name Matematical operations
+        //!@{
+        
+        //! Returns a sample from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv
+        virtual vec samplecond ( const vec &cond ) {
+                this->condition ( cond );
+                vec temp= ep->sample();
+                return temp;
         };
-
-        /*! \brief DataLink is a connection between two data vectors Up and Down
-
-        Up can be longer than Down. Down must be fully present in Up (TODO optional)
-        See chart:
-        \dot
-        digraph datalink {
-                node [shape=record];
-                subgraph cluster0 {
-                        label = "Up";
-                up [label="<1>|<2>|<3>|<4>|<5>"];
-                        color = "white"
-        }
-                subgraph cluster1{
-                        label = "Down";
-                        labelloc = b;
-                down [label="<1>|<2>|<3>"];
-                        color = "white"
-        }
-            up:1 -> down:1;
-            up:3 -> down:2;
-            up:5 -> down:3;
-        }
-        \enddot
-
-        */
-        class datalink_e2e {
-        protected:
-                //! Remember how long val should be
-                int valsize;
-                //! Remember how long val of "Up" should be
-                int valupsize;
-                //! val-to-val link, indeces of the upper val
-                ivec v2v_up;
-        public:
-                //! Constructor
-                datalink_e2e ( const RV &rv, const RV &rv_up ) :
-                                valsize ( rv.count() ), valupsize ( rv_up.count() ), v2v_up ( rv.dataind ( rv_up ) )  {
-                        it_assert_debug ( v2v_up.length() ==valsize,"rv is not fully in rv_up" );
-                }
-                //! Get val for myself from val of "Up"
-                vec get_val ( const vec &val_up ) {it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" ); return get_vec ( val_up,v2v_up );}
-                //! Fill val of "Up" by my pieces
-                void fill_val ( vec &val_up, const vec &val ) {
-                        it_assert_debug ( valsize==val.length(),"Wrong val" );
-                        it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" );
-                        set_subvector ( val_up, v2v_up, val );
-                }
+        //! Returns \param N samples from the density conditioned on \c cond, \f$x \sim epdf(rv|cond)\f$. \param cond is numeric value of \c rv \param ll is a return value of log-likelihood of the sample.
+        virtual mat samplecond_m ( const vec &cond, int N ) {
+                this->condition ( cond );
+                mat temp ( ep->dimension(),N ); vec smp ( ep->dimension() );
+                for ( int i=0;i<N;i++ ) {smp=ep->sample() ;temp.set_col ( i, smp );}
+                return temp;
         };
+        //! Update \c ep so that it represents this mpdf conditioned on \c rvc = cond
+        virtual void condition ( const vec &cond ) {it_error ( "Not implemented" );};
+
+        //! Shortcut for conditioning and evaluation of the internal epdf. In some cases,  this operation can be implemented efficiently.
+        virtual double evallogcond ( const vec &dt, const vec &cond ) {
+                double tmp; this->condition ( cond );tmp = ep->evallog ( dt );          it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); return tmp;
+        };
+
+        //! Matrix version of evallogcond
+        virtual vec evallogcond_m ( const mat &Dt, const vec &cond ) {this->condition ( cond );return ep->evallog_m ( Dt );};
+
+        //! \name Access to attributes
+        //! @{
+        
+        RV _rv() {return ep->_rv();}
+        RV _rvc() {it_assert_debug ( isnamed(),"" ); return rvc;}
+        int dimension() {return ep->dimension();}
+        int dimensionc() {return dimc;}
+        epdf& _epdf() {return *ep;}
+        epdf* _e() {return ep;}
+        //!@}
+        
+        //! \name Connection to other objects
+        //!@{
+        void set_rvc ( const RV &rvc0 ) {rvc=rvc0;}
+        void set_rv ( const RV &rv0 ) {ep->set_rv(rv0);}
+        bool isnamed() {return (ep->isnamed())&&(dimc==rvc._dsize());}
+        //!@}
+};
+
+/*! \brief DataLink is a connection between two data vectors Up and Down
+
+Up can be longer than Down. Down must be fully present in Up (TODO optional)
+See chart:
+\dot
+digraph datalink {
+        node [shape=record];
+        subgraph cluster0 {
+                label = "Up";
+        up [label="<1>|<2>|<3>|<4>|<5>"];
+                color = "white"
+}
+        subgraph cluster1{
+                label = "Down";
+                labelloc = b;
+        down [label="<1>|<2>|<3>"];
+                color = "white"
+}
+    up:1 -> down:1;
+    up:3 -> down:2;
+    up:5 -> down:3;
+}
+\enddot
+
+*/
+class datalink {
+protected:
+        //! Remember how long val should be
+        int downsize;
+        //! Remember how long val of "Up" should be
+        int upsize;
+        //! val-to-val link, indeces of the upper val
+        ivec v2v_up;
+public:
+        //! Constructor
+        datalink ( const RV &rv, const RV &rv_up ) :
+                        downsize ( rv._dsize() ), upsize ( rv_up._dsize() ), v2v_up ( rv.dataind ( rv_up ) )  {
+                it_assert_debug ( v2v_up.length() ==downsize,"rv is not fully in rv_up" );
+        }
+        //! Get val for myself from val of "Up"
+        vec pushdown ( const vec &val_up ) {
+                it_assert_debug ( upsize==val_up.length(),"Wrong val_up" );
+                return get_vec ( val_up,v2v_up );
+        }
+        //! Fill val of "Up" by my pieces
+        void pushup ( vec &val_up, const vec &val ) {
+                it_assert_debug ( downsize==val.length(),"Wrong val" );
+                it_assert_debug ( upsize==val_up.length(),"Wrong val_up" );
+                set_subvector ( val_up, v2v_up, val );
+        }
+};
 
 //! data link between
-        class datalink_m2e: public datalink_e2e {
-        protected:
-                //! Remember how long cond should be
-                int condsize;
-                //!upper_val-to-local_cond link, indeces of the upper val
-                ivec v2c_up;
-                //!upper_val-to-local_cond link, ideces of the local cond
-                ivec v2c_lo;
-
-        public:
-                //! Constructor
-                datalink_m2e ( const RV &rv,  const RV &rvc, const RV &rv_up ) :
-                                datalink_e2e ( rv,rv_up ), condsize ( rvc.count() ) {
-                        //establish v2c connection
-                        rvc.dataind ( rv_up, v2c_lo, v2c_up );
-                }
-                //!Construct condition
-                vec get_cond ( const vec &val_up ) {
-                        vec tmp ( condsize );
-                        set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) );
-                        return tmp;
-                }
-                void fill_val_cond ( vec &val_up, const vec &val, const vec &cond ) {
-                        it_assert_debug ( valsize==val.length(),"Wrong val" );
-                        it_assert_debug ( valupsize==val_up.length(),"Wrong val_up" );
-                        set_subvector ( val_up, v2v_up, val );
-                        set_subvector ( val_up, v2c_up, cond );
-                }
-        };
+class datalink_m2e: public datalink {
+protected:
+        //! Remember how long cond should be
+        int condsize;
+        //!upper_val-to-local_cond link, indeces of the upper val
+        ivec v2c_up;
+        //!upper_val-to-local_cond link, ideces of the local cond
+        ivec v2c_lo;
+
+public:
+        //! Constructor
+        datalink_m2e ( const RV &rv,  const RV &rvc, const RV &rv_up ) :
+                        datalink ( rv,rv_up ), condsize ( rvc._dsize() ) {
+                //establish v2c connection
+                rvc.dataind ( rv_up, v2c_lo, v2c_up );
+        }
+        //!Construct condition
+        vec get_cond ( const vec &val_up ) {
+                vec tmp ( condsize );
+                set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) );
+                return tmp;
+        }
+        void pushup_cond ( vec &val_up, const vec &val, const vec &cond ) {
+                it_assert_debug ( downsize==val.length(),"Wrong val" );
+                it_assert_debug ( upsize==val_up.length(),"Wrong val_up" );
+                set_subvector ( val_up, v2v_up, val );
+                set_subvector ( val_up, v2c_up, cond );
+        }
+};
 //!DataLink is a connection between mpdf and its superordinate (Up)
 //! This class links
-        class datalink_m2m: public datalink_m2e {
-        protected:
-                //!cond-to-cond link, indeces of the upper cond
-                ivec c2c_up;
-                //!cond-to-cond link, indeces of the local cond
-                ivec c2c_lo;
-        public:
-                //! Constructor
-                datalink_m2m ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) :
-                                datalink_m2e ( rv, rvc, rv_up ) {
-                        //establish c2c connection
-                        rvc.dataind ( rvc_up, c2c_lo, c2c_up );
-                        it_assert_debug ( c2c_lo.length() +v2c_lo.length() ==condsize, "cond is not fully given" );
-                }
-                //! Get cond for myself from val and cond of "Up"
-                vec get_cond ( const vec &val_up, const vec &cond_up ) {
-                        vec tmp ( condsize );
-                        set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) );
-                        set_subvector ( tmp,c2c_lo,cond_up ( c2c_up ) );
-                        return tmp;
-                }
-                //! Fill
-
-        };
-
-        /*!
-        @brief Class for storing results (and semi-results) of an experiment
-
-        This class abstracts logging of results from implementation. This class replaces direct logging of results (e.g. to files or to global variables) by calling methods of a logger. Specializations of this abstract class for specific storage method are designed.
-         */
-        class logger : public bdmroot {
-                protected:
-                //! RVs of all logged variables.
-                        Array<RV> entries;
-                //! Names of logged quantities, e.g. names of algorithm variants
-                        Array<string> names;
-                public:
-                //!Default constructor
-                        logger ( ) : entries ( 0 ),names ( 0 ) {}
-
-                //! returns an identifier which will be later needed for calling the log() function
-                        virtual int add ( const RV &rv, string name="" ) {
-                                int id=entries.length();
-                                names=concat ( names, name ); // diff
-                                entries.set_length ( id+1,true );
-                                entries ( id ) = rv;
-                                return id; // identifier of the last entry
-                        }
-
-                //! log this vector
-                        virtual void logit ( int id, const vec &v ) =0;
-
-                //! Shifts storage position for another time step.
-                        virtual void step() =0;
-
-                //! Finalize storing information
-                        virtual void finalize() {};
-
-                //! Initialize the storage
-                        virtual void init() {};
-
-                //! for future use
-                        virtual ~logger() {};
-        };
-
-        /*! \brief Unconditional mpdf, allows using epdf in the role of mpdf.
-
-        */
-        class mepdf : public mpdf {
-        public:
-                //!Default constructor
-                mepdf ( const epdf* em ) :mpdf ( em->_rv(),RV() ) {ep=const_cast<epdf*> ( em );};
-                void condition ( const vec &cond ) {}
-        };
+class datalink_m2m: public datalink_m2e {
+protected:
+        //!cond-to-cond link, indeces of the upper cond
+        ivec c2c_up;
+        //!cond-to-cond link, indeces of the local cond
+        ivec c2c_lo;
+public:
+        //! Constructor
+        datalink_m2m ( const RV &rv, const RV &rvc, const RV &rv_up, const RV &rvc_up ) :
+                        datalink_m2e ( rv, rvc, rv_up ) {
+                //establish c2c connection
+                rvc.dataind ( rvc_up, c2c_lo, c2c_up );
+                it_assert_debug ( c2c_lo.length() +v2c_lo.length() ==condsize, "cond is not fully given" );
+        }
+        //! Get cond for myself from val and cond of "Up"
+        vec get_cond ( const vec &val_up, const vec &cond_up ) {
+                vec tmp ( condsize );
+                set_subvector ( tmp,v2c_lo,val_up ( v2c_up ) );
+                set_subvector ( tmp,c2c_lo,cond_up ( c2c_up ) );
+                return tmp;
+        }
+        //! Fill
+
+};
+
+/*!
+@brief Class for storing results (and semi-results) of an experiment
+
+This class abstracts logging of results from implementation. This class replaces direct logging of results (e.g. to files or to global variables) by calling methods of a logger. Specializations of this abstract class for specific storage method are designed.
+ */
+class logger : public bdmroot {
+protected:
+        //! RVs of all logged variables.
+        Array<RV> entries;
+        //! Names of logged quantities, e.g. names of algorithm variants
+        Array<string> names;
+public:
+        //!Default constructor
+        logger ( ) : entries ( 0 ),names ( 0 ) {}
+
+        //! returns an identifier which will be later needed for calling the log() function
+        virtual int add ( const RV &rv, string name="" ) {
+                int id=entries.length();
+                names=concat ( names, name ); // diff
+                entries.set_length ( id+1,true );
+                entries ( id ) = rv;
+                return id; // identifier of the last entry
+        }
+
+        //! log this vector
+        virtual void logit ( int id, const vec &v ) =0;
+
+        //! Shifts storage position for another time step.
+        virtual void step() =0;
+
+        //! Finalize storing information
+        virtual void finalize() {};
+
+        //! Initialize the storage
+        virtual void init() {};
+
+};
+
+/*! \brief Unconditional mpdf, allows using epdf in the role of mpdf.
+
+*/
+class mepdf : public mpdf {
+public:
+        //!Default constructor
+        mepdf ( const epdf* em ) :mpdf ( ) {ep=const_cast<epdf*> ( em );};
+        void condition ( const vec &cond ) {}
+};
 
 //!\brief Abstract composition of pdfs, will be used for specific classes
 //!this abstract class is common to epdf and mpdf
-        class compositepdf {
-        protected:
-                //!Number of mpdfs in the composite
-                int n;
-                //! Elements of composition
-                Array<mpdf*> mpdfs;
-        public:
-                compositepdf ( Array<mpdf*> A0 ) : n ( A0.length() ), mpdfs ( A0 ) {};
-                //! find common rv, flag \param checkoverlap modifies whether overlaps are acceptable
-                RV getrv ( bool checkoverlap=false );
-                //! common rvc of all mpdfs is written to rvc
-                void setrvc ( const RV &rv, RV &rvc );
-        };
-
-        /*! \brief Abstract class for discrete-time sources of data.
-
-        The class abstracts operations of: (i) data aquisition, (ii) data-preprocessing, (iii) scaling of data, and (iv) data resampling from the task of estimation and control.
-        Moreover, for controlled systems, it is able to receive the desired control action and perform it in the next step. (Or as soon as possible).
-
+class compositepdf {
+protected:
+        //!Number of mpdfs in the composite
+        int n;
+        //! Elements of composition
+        Array<mpdf*> mpdfs;
+public:
+        compositepdf ( Array<mpdf*> A0 ) : n ( A0.length() ), mpdfs ( A0 ) {};
+        //! find common rv, flag \param checkoverlap modifies whether overlaps are acceptable
+        RV getrv ( bool checkoverlap=false );
+        //! common rvc of all mpdfs is written to rvc
+        void setrvc ( const RV &rv, RV &rvc );
+};
+
+/*! \brief Abstract class for discrete-time sources of data.
+
+The class abstracts operations of: (i) data aquisition, (ii) data-preprocessing, (iii) scaling of data, and (iv) data resampling from the task of estimation and control.
+Moreover, for controlled systems, it is able to receive the desired control action and perform it in the next step. (Or as soon as possible).
+
+*/
+
+class DS : public bdmroot {
+protected:
+        int dtsize;
+        int utsize;
+        //!Description of data returned by \c getdata().
+        RV Drv;
+        //!Description of data witten by by \c write().
+        RV Urv; //
+        //! Remember its own index in Logger L
+        int L_dt, L_ut;
+public:
+        //! default constructors
+        DS() :Drv (  ),Urv ( ) {};
+        //! Returns full vector of observed data=[output, input]
+        virtual void getdata ( vec &dt ) {it_error ( "abstract class" );};
+        //! Returns data records at indeces.
+        virtual void getdata ( vec &dt, const ivec &indeces ) {it_error ( "abstract class" );};
+        //! Accepts action variable and schedule it for application.
+        virtual void write ( vec &ut ) {it_error ( "abstract class" );};
+        //! Accepts action variables at specific indeces
+        virtual void write ( vec &ut, const ivec &indeces ) {it_error ( "abstract class" );};
+
+        //! Moves from \f$ t \f$ to \f$ t+1 \f$, i.e. perfroms the actions and reads response of the system.
+        virtual void step() =0;
+
+        //! Register DS for logging into logger L
+        virtual void log_add ( logger &L ) {
+                it_assert_debug ( dtsize==Drv._dsize(),"" );
+                it_assert_debug ( utsize==Urv._dsize(),"" );
+
+                L_dt=L.add ( Drv,"" );
+                L_ut=L.add ( Urv,"" );
+        }
+        //! Register DS for logging into logger L
+        virtual void logit ( logger &L ) {
+                vec tmp ( Drv._dsize() +Urv._dsize() );
+                getdata ( tmp );
+                // d is first in getdata
+                L.logit ( L_dt,tmp.left ( Drv._dsize() ) );
+                // u follows after d in getdata
+                L.logit ( L_ut,tmp.mid ( Drv._dsize(), Urv._dsize() ) );
+        }
+        //!access function
+        virtual RV _drv() const {return concat ( Drv,Urv );}
+        //!access function
+        const RV& _urv() const {return Urv;}
+};
+
+/*! \brief Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities.
+
+*/
+
+class BM :public bdmroot {
+protected:
+        //! Random variable of the data (optional)
+        RV drv;
+        //!Logarithm of marginalized data likelihood.
+        double ll;
+        //!  If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save computational time.
+        bool evalll;
+public:
+        //! \name Constructors
+        //! @{
+        
+        BM () :ll (0),evalll ( false) {};
+        BM ( const BM &B ) :  drv(B.drv), ll ( B.ll ), evalll ( B.evalll ) {}
+        //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
+        //! Prototype: \code BM* _copy_(){return new BM(*this);} \endcode
+        virtual BM* _copy_ () {return NULL;};
+        //!@}
+
+        //! \name Mathematical operations
+        //!@{
+        
+        /*! \brief Incremental Bayes rule
+        @param dt vector of input data
         */
-
-        class DS : public bdmroot {
-        protected:
-                //!Observed variables, returned by \c getdata().
-                RV Drv;
-                //!Action variables, accepted by \c write().
-                RV Urv; //
-                //! Remember its own index in Logger L
-                int L_dt, L_ut;
-        public:
-                DS() :Drv ( RV0 ),Urv ( RV0 ) {};
-                DS ( const RV &Drv0, const RV &Urv0 ) :Drv ( Drv0 ),Urv ( Urv0 ) {};
-                //! Returns full vector of observed data=[output, input]
-                virtual void getdata ( vec &dt ) {it_error ( "abstract class" );};
-                //! Returns data records at indeces.
-                virtual void getdata ( vec &dt, const ivec &indeces ) {it_error ( "abstract class" );};
-                //! Accepts action variable and schedule it for application.
-                virtual void write ( vec &ut ) {it_error ( "abstract class" );};
-                //! Accepts action variables at specific indeces
-                virtual void write ( vec &ut, const ivec &indeces ) {it_error ( "abstract class" );};
-
-                //! Moves from \f$t\f$ to \f$t+1\f$, i.e. perfroms the actions and reads response of the system.
-                virtual void step() =0;
-
-                //! Register DS for logging into logger L
-                virtual void log_add ( logger &L ) {
-                        L_dt=L.add ( Drv,"" );
-                        L_ut=L.add ( Urv,"" );
-                }
-                //! Register DS for logging into logger L
-                virtual void logit ( logger &L ) {
-                        vec tmp(Drv.count()+Urv.count());
-                        getdata(tmp);
-                        // d is first in getdata
-                        L.logit ( L_dt,tmp.left ( Drv.count() ) );
-                        // u follows after d in getdata
-                        L.logit ( L_ut,tmp.mid ( Drv.count(), Urv.count() ) );
-                }
-                //!access function 
-                virtual RV _drv() const {return concat(Drv,Urv);}
-                //!access function 
-                const RV& _urv() const {return Urv;}
-        };
-
-        /*! \brief Bayesian Model of a system, i.e. all uncertainty is modeled by probabilities.
-
-        */
-
-        class BM :public bdmroot {
-        protected:
-                //!Random variable of the posterior
-                RV rv;
-                //! Random variable of the data (optional)
-                RV drv;
-                //!Logarithm of marginalized data likelihood.
-                double ll;
-                //!  If true, the filter will compute likelihood of the data record and store it in \c ll . Set to false if you want to save computational time.
-                bool evalll;
-        public:
-
-                //!Default constructor
-                BM ( const RV &rv0, double ll0=0,bool evalll0=true ) :rv ( rv0 ), ll ( ll0 ),evalll ( evalll0 ) {//Fixme: test rv
-                };
-                //!Copy constructor
-                BM ( const BM &B ) : rv ( B.rv ), ll ( B.ll ), evalll ( B.evalll ) {}
-
-                /*! \brief Incremental Bayes rule
-                @param dt vector of input data
-                */
-                virtual void bayes ( const vec &dt ) = 0;
-                //! Batch Bayes rule (columns of Dt are observations)
-                virtual void bayesB ( const mat &Dt );
-                //! Returns a reference to the epdf representing posterior density on parameters.
-                virtual const epdf& _epdf() const =0;
-
-                //! Returns a pointer to the epdf representing posterior density on parameters. Use with care!
-                virtual const epdf* _e() const =0;
-
-                //! Evaluates predictive log-likelihood of the given data record
-                //! I.e. marginal likelihood of the data with the posterior integrated out.
-                virtual double logpred ( const vec &dt ) const{it_error ( "Not implemented" );return 0.0;}
-                //! Matrix version of logpred
-                vec logpred_m ( const mat &dt ) const{vec tmp ( dt.cols() );for ( int i=0;i<dt.cols();i++ ) {tmp ( i ) =logpred ( dt.get_col ( i ) );}return tmp;}
-
-                //!Constructs a predictive density (marginal density on data)
-                virtual epdf* predictor ( const RV &rv ) const {it_error ( "Not implemented" );return NULL;};
-
-                //! Destructor for future use;
-                virtual ~BM() {};
-                //!access function
-                const RV& _rv() const {return rv;}
-                //!access function 
-                const RV& _drv() const {return drv;}
-                //!set drv
-                void set_drv(const RV &rv){drv=rv;}
-                //!access function
-                double _ll() const {return ll;}
-                //!access function
-                void set_evalll ( bool evl0 ) {evalll=evl0;}
-
-                //! Copy function required in vectors, Arrays of BM etc. Have to be DELETED manually!
-                //! Prototype: BM* _copy_(){BM Tmp*=new Tmp(this*);  return Tmp; }
-                virtual BM* _copy_ ( bool changerv=false ) {it_error ( "function _copy_ not implemented for this BM" ); return NULL;};
-        };
-
-        /*!
-        \brief Conditional Bayesian Filter
-
-        Evaluates conditional filtering density \f$f(rv|rvc,data)\f$ for a given \c rvc which is specified in each step by calling function \c condition.
-
-        This is an interface class used to assure that certain BM has operation \c condition .
-
-        */
-
-        class BMcond :public bdmroot {
-        protected:
-                //! Identificator of the conditioning variable
-                RV rvc;
-        public:
-                //! Substitute \c val for \c rvc.
-                virtual void condition ( const vec &val ) =0;
-                //! Default constructor
-                BMcond ( RV &rv0 ) :rvc ( rv0 ) {};
-                //! Destructor for future use
-                virtual ~BMcond() {};
-                //! access function
-                const RV& _rvc() const {return rvc;}
-        };
+        virtual void bayes ( const vec &dt ) = 0;
+        //! Batch Bayes rule (columns of Dt are observations)
+        virtual void bayesB ( const mat &Dt );
+        //! Evaluates predictive log-likelihood of the given data record
+        //! I.e. marginal likelihood of the data with the posterior integrated out.
+        virtual double logpred ( const vec &dt ) const{it_error ( "Not implemented" );return 0.0;}
+        //! Matrix version of logpred
+        vec logpred_m ( const mat &dt ) const{vec tmp ( dt.cols() );for ( int i=0;i<dt.cols();i++ ) {tmp ( i ) =logpred ( dt.get_col ( i ) );}return tmp;}
+
+        //!Constructs a predictive density \f$ f(d_{t+1} |d_{t}, \ldots d_{0}) \f$
+        virtual epdf* epredictor (  ) const {it_error ( "Not implemented" );return NULL;};
+        //!Constructs a conditional density 1-step ahead predictor \f$ f(d_{t+1} |d_{t+h-1}, \ldots d_{t})
+        virtual mpdf* predictor ( ) const {it_error ( "Not implemented" );return NULL;};
+        //!@}
+        
+        //! \name Access to attributes
+        //!@{
+        
+        const RV& _drv() const {return drv;}
+        void set_drv ( const RV &rv ) {drv=rv;}
+        double _ll() const {return ll;}
+        void set_evalll ( bool evl0 ) {evalll=evl0;}
+        virtual const epdf& _epdf() const =0;
+        virtual const epdf* _e() const =0;
+        //!@}
+
+};
+
+/*!
+\brief Conditional Bayesian Filter
+
+Evaluates conditional filtering density \f$f(rv|rvc,data)\f$ for a given \c rvc which is specified in each step by calling function \c condition.
+
+This is an interface class used to assure that certain BM has operation \c condition .
+
+*/
+
+class BMcond :public bdmroot {
+protected:
+        //! Identificator of the conditioning variable
+        RV rvc;
+public:
+        //! Substitute \c val for \c rvc.
+        virtual void condition ( const vec &val ) =0;
+        //! Default constructor
+        BMcond ( ) :rvc ( ) {};
+        //! Destructor for future use
+        virtual ~BMcond() {};
+        //! access function
+        const RV& _rvc() const {return rvc;}
+};
 
 }; //namespace

bdm/stat/libDS.cpp

@@ -26 +26 @@
 }
 
-void MemDS::linkrvs ( RV &drv, RV &urv ) {
-        it_assert_debug ( drv.count() ==rowid.length(),"MemDS::linkrvs incompatible drv" );
-        it_assert_debug ( urv.count() ==0,"MemDS does not support urv." );
+void MemDS::set_rvs ( RV &drv, RV &urv ) {
+        it_assert_debug ( drv._dsize() ==rowid.length(),"MemDS::set_rvs incompatible drv" );
+        it_assert_debug ( urv._dsize() ==0,"MemDS does not support urv." );
 
         Drv = drv;
-        Urv = urv;
 }
 
@@ -44 +43 @@
 void ArxDS::step() {
         //shift history
-        H.shift_right ( 0, Drv.count() +Urv.count() );
+        H.shift_right ( 0, Drv._dsize() +Urv._dsize() );
 
-        H.set_subvector ( Drv.count(), U ); // write U after Drv
+        H.set_subvector ( Drv._dsize(), U ); // write U after Drv
 
         //get regressor
-        rgr = rgrlnk.get_val ( H );
+        rgr = rgrlnk->pushdown ( H );
         // Eval Y
         H.set_subvector ( 0, model.samplecond ( rgr ) );
@@ -66 +65 @@
         return T;
 }
-
-ArxDS::ArxDS ( RV &drv, RV &urv, RV &rrv ) :
-                DS ( drv,urv ), Rrv ( rrv ), Hrv ( concat(drv,fullrgr ( drv,urv,rrv )) ), //RVs
-                H ( Hrv.count() ) ,U ( urv.count() ),rgr ( rrv.count() ),  //tmp variables
-                rgrlnk ( rrv, Hrv ) ,model ( drv,rrv ) {
-}

bdm/stat/libDS.h

@@ -23 +23 @@
         * \brief Memory storage of off-line data column-wise
 
-        The data are stored in an internal matrix \c Data . Each column of Data corresponds to one discrete time observation \f$t\f$. Access to this matrix is via indexes \c rowid and \c delays.
+        The data are stored in an internal matrix \c Data . Each column of Data corresponds to one discrete time observation \f$t\f$. Access to this matrix is via indices \c rowid and \c delays.
 
         The data can be loaded from a file.
@@ -40 +40 @@
                 void getdata ( vec &dt );
                 void getdata ( vec &dt, const ivec &indeces );
-                void linkrvs ( RV &drv, RV &urv );
+                void set_rvs ( RV &drv, RV &urv );
                 void write ( vec &ut ) {it_error ( "MemDS::write is not supported" );}
-                void write ( vec &ut,ivec &indexes ) {it_error ( "MemDS::write is not supported" );}
+                void write ( vec &ut,ivec &indices ) {it_error ( "MemDS::write is not supported" );}
                 void step();
                 //!Default constructor
@@ -65 +65 @@
                 vec rgr;
                 //! data link: H -> rgr
-                datalink_e2e rgrlnk;
+                datalink *rgrlnk;
                 //! model of Y - linear Gaussian
                 mlnorm<chmat> model;
@@ -75 +75 @@
         public:
                 void getdata ( vec &dt ) {
-                        it_assert_debug ( dt.length() ==Drv.count(),"ArxDS" );
-                        dt=H.left ( Urv.count() +Drv.count() );
+                        //it_assert_debug ( dt.length() ==Drv.count(),"ArxDS" );
+                        dt=H.left ( Urv._dsize() +Drv._dsize() );
                 };
-                void getdata ( vec &dt, const ivec &indexes ) {
-                        it_assert_debug ( dt.length() ==indeces.length(),"ArxDS" );
-                        dt=H ( indexes );
+                void getdata ( vec &dt, const ivec &indices ) {
+                        it_assert_debug ( dt.length() ==indices.length(),"ArxDS" );
+                        dt=H ( indices );
                 };
                 void write ( vec &ut ) {
-                        it_assert_debug ( ut.length() ==Urv.count(),"ArxDS" );
+                        //it_assert_debug ( ut.length() ==Urv.count(),"ArxDS" );
                         U=ut;
                 };
-                void write ( vec &ut, const ivec &indexes ) {
-                        it_assert_debug ( ut.length() ==indeces.length(),"ArxDS" );
-                        set_subvector ( U, indexes,ut );
+                void write ( vec &ut, const ivec &indices ) {
+                        it_assert_debug ( ut.length() ==indices.length(),"ArxDS" );
+                        set_subvector ( U, indices,ut );
                 };
                 void step();
                 //!Default constructor
-                ArxDS ( RV &drv, RV &urv, RV &rrv );
+                ArxDS ( ){};
                 //! Set parameters of the internal model
                 void set_parameters ( const mat &Th0, const vec mu0, const chmat &sqR0 )
@@ -119 +119 @@
         class ARXDS : public ArxDS {
         public:
-                ARXDS ( RV &drv, RV &urv, RV &rrv ) : ArxDS ( drv,urv,rrv ) {}
+                ARXDS ( ) : ArxDS ( ) {}
 
                 void getdata ( vec &dt ) {dt=H;}
@@ -145 +145 @@
                 void getdata ( vec &dt0, const ivec &indeces ) {dt0=dt ( indeces );}
 
-                stateDS ( mpdf* IM0, mpdf* OM0, RV &Urv0 ) :DS ( OM0->_rv(),Urv0 ),IM ( IM0 ),OM ( OM0 ),
-                                dt ( OM0->_rv().count() ), xt ( IM0->_rv().count() ), ut ( Urv0.count() ) {}
+                stateDS ( mpdf* IM0, mpdf* OM0, int usize ) :DS ( ),IM ( IM0 ),OM ( OM0 ),
+                                dt ( OM0->dimension() ), xt ( IM0->dimension() ), ut ( usize) {}
                 ~stateDS() {delete IM; delete OM;}
                 virtual void step() {
-                        double tmp;
                         xt=IM->samplecond(concat ( xt,ut ));
                         dt=OM->samplecond(concat ( xt,ut ));

bdm/stat/libEF.cpp

@@ -22 +22 @@
         int vend = val.length()-1;
 
-        if ( xdim==1 ) { //same as the following, just quicker.
+        if ( dimx==1 ) { //same as the following, just quicker.
                 double r = val ( vend ); //last entry!
-                vec Psi ( nPsi+xdim );
+                vec Psi ( nPsi+dimx );
                 Psi ( 0 ) = -1.0;
                 Psi.set_subvector ( 1,val ( 0,vend-1 ) ); // fill the rest
@@ -32 +32 @@
         }
         else {
-                mat Th= reshape ( val ( 0,nPsi*xdim-1 ),nPsi,xdim );
-                fsqmat R ( reshape ( val ( nPsi*xdim,vend ),xdim,xdim ) );
-                mat Tmp=concat_vertical ( -eye ( xdim ),Th );
-                fsqmat iR ( xdim );
+                mat Th= reshape ( val ( 0,nPsi*dimx-1 ),nPsi,dimx );
+                fsqmat R ( reshape ( val ( nPsi*dimx,vend ),dimx,dimx ) );
+                mat Tmp=concat_vertical ( -eye ( dimx ),Th );
+                fsqmat iR ( dimx );
                 R.inv ( iR );
 
@@ -45 +45 @@
         const vec& D = V._D();
 
-        double m = nu - nPsi -xdim-1;
+        double m = nu - nPsi -dimx-1;
 #define log2  0.693147180559945286226763983
 #define logpi 1.144729885849400163877476189
@@ -51 +51 @@
 #define Inf std::numeric_limits<double>::infinity()
 
-        double nkG = 0.5* xdim* ( -nPsi *log2pi + sum ( log ( D ( xdim,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<xdim;i++ ) {lg+=lgamma ( 0.5* ( m-i ) );}
-
-        double nkW = 0.5* ( m*sum ( log ( D ( 0,xdim-1 ) ) ) ) \
-                     - 0.5*xdim* ( m*log2 + 0.5* ( xdim-1 ) *log2pi )  - lg;
+        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" );
@@ -65 +65 @@
 vec egiw::mean() const {
 
-        if ( xdim==1 ) {
+        if ( dimx==1 ) {
                 const mat &L= V._L();
                 const vec &D= V._D();
                 int end = L.rows()-1;
 
-                vec m ( rv.count() );
-                mat iLsub = ltuinv ( L ( xdim,end,xdim,end ) );
+                vec m ( dim );
+                mat iLsub = ltuinv ( L ( dimx,end,dimx,end ) );
 
                 vec L0 = L.get_col ( 0 );
 
                 m.set_subvector ( 0,iLsub*L0 ( 1,end ) );
-                m ( end ) = D ( 0 ) / ( nu -nPsi -2*xdim -2 );
+                m ( end ) = D ( 0 ) / ( nu -nPsi -2*dimx -2 );
                 return m;
         }
@@ -90 +90 @@
 vec egiw::variance() const {
 
-        if ( xdim==1 ) {
+        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*xdim -2 );
+                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*xdim -2 );
+                var(l-1)=cove*cove/( nu -nPsi -2*dimx -2 );
                 return var;
         }
@@ -110 +110 @@
         int end = L.rows()-1;
 
-        ldmat ldR ( L ( 0,xdim-1,0,xdim-1 ), D ( 0,xdim-1 ) / ( nu -nPsi -2*xdim -2 ) ); //exp val of R
-        mat iLsub=ltuinv ( L ( xdim,end,xdim,end ) );
+        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 ( xdim,end,0,xdim-1 );
+        mat Lpsi = L ( dimx,end,0,dimx-1 );
         M= iLsub*Lpsi;
         R= ldR.to_mat()  ;
@@ -120 +120 @@
 
 vec egamma::sample() const {
-        vec smp ( rv.count() );
+        vec smp ( dim);
         int i;
 
-        for ( i=0; i<rv.count(); i++ ) {
+        for ( i=0; i<dim; i++ ) {
                 if ( beta ( i ) >std::numeric_limits<double>::epsilon() ) {
                         GamRNG.setup ( alpha ( i ),beta ( i ) );
@@ -156 +156 @@
         int i;
 
-        for ( i=0; i<rv.count(); i++ ) {
+        for ( i=0; i<dim; i++ ) {
                 res += ( alpha ( i ) - 1 ) *std::log ( val ( i ) ) - beta ( i ) *val ( i );
         }
@@ -168 +168 @@
         int i;
 
-        for ( i=0; i<rv.count(); i++ ) {
+        for ( i=0; i<dim; i++ ) {
                 res += lgamma ( alpha ( i ) ) - alpha ( i ) *std::log ( beta ( i ) ) ;
         }
@@ -176 +176 @@
 
 //MGamma
-void mgamma::set_parameters ( double k0 ) {
+void mgamma::set_parameters ( double k0, const vec &beta0 ) {
         k=k0;
         ep = &epdf;
-        epdf.set_parameters ( k*ones ( rv.count() ),*_beta );
+        epdf.set_parameters ( k*ones ( beta0.length()),beta0 );
+        dimc = ep->dimension();
 };
 
@@ -273 +274 @@
         n=w.length();
         samples.set_size ( n );
+        dim = epdf0->dimension();
 
         for ( int i=0;i<n;i++ ) {samples ( i ) =epdf0->sample();}

bdm/stat/libEF.h

@@ -19 +19 @@
 //#include <std>
 
-namespace bdm{
+namespace bdm {
 
 
@@ -39 +39 @@
 //      eEF() :epdf() {};
         //! default constructor
-        eEF ( const RV &rv ) :epdf ( rv ) {};
+        eEF ( ) :epdf ( ) {};
         //! logarithm of the normalizing constant, \f$\mathcal{I}\f$
         virtual double lognc() const =0;
@@ -47 +47 @@
         virtual double evallog_nn ( const vec &val ) const{it_error ( "Not implemented" );return 0.0;};
         //!Evaluate normalized log-probability
-        virtual double evallog ( const vec &val ) const {double tmp;tmp= evallog_nn ( val )-lognc();it_assert_debug(std::isfinite(tmp),"Infinite value"); return tmp;}
+        virtual double evallog ( const vec &val ) const {double tmp;tmp= evallog_nn ( val )-lognc();it_assert_debug ( std::isfinite ( tmp ),"Infinite value" ); return tmp;}
         //!Evaluate normalized log-probability for many samples
         virtual vec evallog ( const mat &Val ) const {
@@ -68 +68 @@
 public:
         //! Default constructor
-        mEF ( const RV &rv0, const RV &rvc0 ) :mpdf ( rv0,rvc0 ) {};
+        mEF ( ) :mpdf ( ) {};
 };
 
@@ -79 +79 @@
         double last_lognc;
 public:
-        //! Default constructor
-        BMEF ( const RV &rv, double frg0=1.0 ) :BM ( rv ), frg ( frg0 ) {}
+        //! Default constructor (=empty constructor)
+        BMEF ( double frg0=1.0 ) :BM ( ), frg ( frg0 ) {}
         //! Copy constructor
         BMEF ( const BMEF &B ) :BM ( B ), frg ( B.frg ), last_lognc ( B.last_lognc ) {}
@@ -112 +112 @@
         //! Covariance matrix in decomposed form
         sq_T R;
-        //! dimension (redundant from rv.count() for easier coding )
-        int dim;
-public:
-        //!Default constructor
-        enorm ( const RV &rv );
-        //! Set mean value \c mu and covariance \c R
+public:
+        //!\name Constructors
+        //!@{
+
+        enorm ( ):eEF ( ), mu ( ),R ( ) {};
+        enorm ( const vec &mu,const sq_T &R ) {set_parameters ( mu,R );}
         void set_parameters ( const vec &mu,const sq_T &R );
-        ////! tupdate in exponential form (not really handy)
-        //void tupdate ( double phi, mat &vbar, double nubar );
+        //!@}
+
+        //! \name Mathematical operations
+        //!@{
+
         //! dupdate in exponential form (not really handy)
         void dupdate ( mat &v,double nu=1.0 );
 
         vec sample() const;
-        //! TODO is it used?
         mat sample ( int N ) const;
-//      double eval ( const vec &val ) const ;
         double evallog_nn ( const vec &val ) const;
         double lognc () const;
         vec mean() const {return mu;}
-        vec variance() const {return diag(R.to_mat());}
+        vec variance() const {return diag ( R.to_mat() );}
 //      mlnorm<sq_T>* condition ( const RV &rvn ) const ;
         mpdf* condition ( const RV &rvn ) const ;
 //      enorm<sq_T>* marginal ( const RV &rv ) const;
         epdf* marginal ( const RV &rv ) const;
-//Access methods
-        //! returns a pointer to the internal mean value. Use with Care!
+        //!@}
+
+        //! \name Access to attributes
+        //!@{
+
         vec& _mu() {return mu;}
-
-        //! access function
         void set_mu ( const vec mu0 ) { mu=mu0;}
-
-        //! returns pointers to the internal variance and its inverse. Use with Care!
         sq_T& _R() {return R;}
         const sq_T& _R() const {return R;}
-
-        //! access method
-//      mat getR () {return R.to_mat();}
+        //!@}
+
 };
 
@@ -164 +163 @@
         double nu;
         //! Dimension of the output
-        int xdim;
+        int dimx;
         //! Dimension of the regressor
         int nPsi;
 public:
-        //!Default constructor, if nu0<0 a minimal nu0 will be computed
-        egiw ( RV rv, mat V0, double nu0=-1.0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
-                xdim = rv.count() /V.rows();
-                it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
-                nPsi = V.rows()-xdim;
-                //set mu to have proper normalization and 
-                if (nu0<0){
-                        nu = 0.1 +nPsi +2*xdim +2; // +2 assures finite expected value of R
+        //!\name Constructors
+        //!@{
+        egiw() :eEF() {};
+        egiw ( int dimx0, ldmat V0, double nu0=-1.0 ) :eEF() {set_parameters ( dimx0,V0, nu0 );};
+
+        void set_parameters ( int dimx0, ldmat V0, double nu0=-1.0 ) {
+                dimx=dimx0;
+                nPsi = V0.rows()-dimx;
+                dim = dimx* ( dimx+nPsi ); // size(R) + size(Theta)
+
+                V=V0;
+                if ( nu0<0 ) {
+                        nu = 0.1 +nPsi +2*dimx +2; // +2 assures finite expected value of R
                         // terms before that are sufficient for finite normalization
                 }
-        }
-        //!Full constructor for V in ldmat form
-        egiw ( RV rv, ldmat V0, double nu0=-1.0 ) : eEF ( rv ), V ( V0 ), nu ( nu0 ) {
-                xdim = rv.count() /V.rows();
-                it_assert_debug ( rv.count() ==xdim*V.rows(),"Incompatible V0." );
-                nPsi = V.rows()-xdim;
-                if (nu0<0){
-                        nu = 0.1 +nPsi +2*xdim +2; // +2 assures finite expected value of R
-                        // terms before that are sufficient for finite normalization
+                else {
+                        nu=nu0;
                 }
         }
+        //!@}
 
         vec sample() const;
@@ -197 +195 @@
         double evallog_nn ( const vec &val ) const;
         double lognc () const;
-
-        //Access
-        //! returns a pointer to the internal statistics. Use with Care!
+        void pow ( double p ) {V*=p;nu*=p;};
+
+        //! \name Access attributes
+        //!@{
+
         ldmat& _V() {return V;}
-        //! returns a pointer to the internal statistics. Use with Care!
         const ldmat& _V() const {return V;}
-        //! returns a pointer to the internal statistics. Use with Care!
         double& _nu()  {return nu;}
         const double& _nu() const {return nu;}
-        void pow ( double p ) {V*=p;nu*=p;};
+        //!@}
 };
 
@@ -224 +222 @@
         double gamma;
 public:
-        //!Default constructor
-        eDirich ( const RV &rv, const vec &beta0 ) : eEF ( rv ),beta ( beta0 ) {it_assert_debug ( rv.count() ==beta.length(),"Incompatible statistics" ); };
-        //! Copy constructor
-        eDirich ( const eDirich &D0 ) : eEF ( D0.rv ),beta ( D0.beta ) {};
+        //!\name Constructors
+        //!@{
+
+        eDirich () : eEF ( ) {};
+        eDirich ( const eDirich &D0 ) : eEF () {set_parameters ( D0.beta );};
+        eDirich ( const vec &beta0 ) {set_parameters ( beta0 );};
+        void set_parameters ( const vec &beta0 ) {
+                beta= beta0;
+                dim = beta.length();
+                gamma = sum ( beta );
+        }
+        //!@}
+
         vec sample() const {it_error ( "Not implemented" );return vec_1 ( 0.0 );};
         vec mean() const {return beta/gamma;};
-        vec variance() const {return elem_mult(beta,(beta+1))/ (gamma*(gamma+1));}
+        vec variance() const {return elem_mult ( beta, ( beta+1 ) ) / ( gamma* ( gamma+1 ) );}
         //! In this instance, val is ...
-        double evallog_nn ( const vec &val ) const {double tmp; tmp=( beta-1 ) *log ( val );            it_assert_debug(std::isfinite(tmp),"Infinite value");
-        return tmp;};
+        double evallog_nn ( const vec &val ) const {
+                double tmp; tmp= ( beta-1 ) *log ( val );               it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
+                return tmp;
+        };
         double lognc () const {
                 double tmp;
@@ -240 +249 @@
                 for ( int i=0;i<beta.length();i++ ) {lgb+=lgamma ( beta ( i ) );}
                 tmp= lgb-lgamma ( gam );
-                it_assert_debug(std::isfinite(tmp),"Infinite value");
+                it_assert_debug ( std::isfinite ( tmp ),"Infinite value" );
                 return tmp;
         };
@@ -246 +255 @@
         vec& _beta()  {return beta;}
         //!Set internal parameters
-        void set_parameters ( const vec &beta0 ) {
-                if ( beta0.length() !=beta.length() ) {
-                        it_assert_debug ( rv.length() ==1,"Undefined" );
-                        rv.set_size ( 0,beta0.length() );
-                }
-                beta= beta0;
-                gamma = sum(beta);
-        }
 };
 
@@ -265 +266 @@
 public:
         //!Default constructor
-        multiBM ( const RV &rv, const vec beta0 ) : BMEF ( rv ),est ( rv,beta0 ),beta ( est._beta() ) {if(beta.length()>0){last_lognc=est.lognc();}else{last_lognc=0.0;}}
+        multiBM ( ) : BMEF ( ),est ( ),beta ( est._beta() ) {if ( beta.length() >0 ) {last_lognc=est.lognc();}else{last_lognc=0.0;}}
         //!Copy constructor
-        multiBM ( const multiBM &B ) : BMEF ( B ),est ( rv,B.beta ),beta ( est._beta() ) {}
+        multiBM ( const multiBM &B ) : BMEF ( B ),est ( B.est ),beta ( est._beta() ) {}
         //! Sets sufficient statistics to match that of givefrom mB0
         void set_statistics ( const BM* mB0 ) {const multiBM* mB=dynamic_cast<const multiBM*> ( mB0 ); beta=mB->beta;}
@@ -300 +301 @@
         void set_parameters ( const vec &beta0 ) {
                 est.set_parameters ( beta0 );
-                rv = est._rv();
                 if ( evalll ) {last_lognc=est.lognc();}
         }
@@ -321 +321 @@
         vec beta;
 public :
-        //! Default constructor
-        egamma ( const RV &rv ) :eEF ( rv ), alpha(rv.count()), beta(rv.count()) {};
-        //! Sets parameters
-        void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;};
+        //! \name Constructors
+        //!@{
+        egamma ( ) :eEF ( ), alpha(), beta() {};
+        egamma ( const vec &a, const vec &b ) {set_parameters ( a, b );};
+        void set_parameters ( const vec &a, const vec &b ) {alpha=a,beta=b;dim = alpha.length();};
+        //!@}
+
         vec sample() const;
         //! TODO: is it used anywhere?
@@ -331 +334 @@
         double lognc () const;
         //! Returns poiter to alpha and beta. Potentially dengerous: use with care!
-        void _param ( vec* &a, vec* &b ) {a=&alpha;b=&beta;};
-        vec mean() const {return elem_div(alpha,beta);}
-        vec variance() const {return elem_div(alpha,elem_mult(beta,beta)); }
+        vec& _alpha() {return alpha;}
+        vec& _beta() {return beta;}
+        vec mean() const {return elem_div ( alpha,beta );}
+        vec variance() const {return elem_div ( alpha,elem_mult ( beta,beta ) ); }
 };
 
@@ -351 +355 @@
 
 class eigamma : public eEF {
-        protected:
+protected:
+        //!internal egamma
+        egamma eg;
         //! Vector \f$\alpha\f$
-                vec* alpha;
+        vec &alpha;
         //! Vector \f$\beta\f$ (in fact it is 1/beta as used in definition of iG)
-                vec* beta;
-                //!internal egamma
-                egamma eg;
-        public :
-        //! Default constructor
-                eigamma ( const RV &rv ) :eEF ( rv ), eg(rv) {eg._param(alpha,beta);};
-        //! Sets parameters
-                void set_parameters ( const vec &a, const vec &b ) {*alpha=a,*beta=b;};
-                vec sample() const {return 1.0/eg.sample();};
+        vec &beta;
+public :
+        //! \name Constructors
+        //!@{
+        eigamma ( ) :eEF ( ), eg(),alpha ( eg._alpha() ), beta ( eg._beta() ) {};
+        eigamma ( const vec &a, const vec &b ):eEF ( ), eg(),alpha ( eg._alpha() ), beta ( eg._beta() ) {eg.set_parameters ( a,b );};
+        void set_parameters ( const vec &a, const vec &b ) {eg.set_parameters ( a,b );};
+        //!@}
+
+        vec sample() const {return 1.0/eg.sample();};
         //! TODO: is it used anywhere?
 //      mat sample ( int N ) const;
-                double evallog ( const vec &val ) const {return eg.evallog(val);};
-                double lognc () const {return eg.lognc();};
+        double evallog ( const vec &val ) const {return eg.evallog ( val );};
+        double lognc () const {return eg.lognc();};
         //! Returns poiter to alpha and beta. Potentially dangerous: use with care!
-                void _param ( vec* &a, vec* &b ) {a=alpha;b=beta;};
-                vec mean() const {return elem_div(*beta,*alpha-1);}
-                vec variance() const {vec mea=mean(); return elem_div(elem_mult(mea,mea),*alpha-2);}
+        vec mean() const {return elem_div ( beta,alpha-1 );}
+        vec variance() const {vec mea=mean(); return elem_div ( elem_mult ( mea,mea ),alpha-2 );}
+        vec& _alpha() {return alpha;}
+        vec& _beta() {return beta;}
 };
 /*
@@ -404 +412 @@
         double lnk;
 public:
-        //! Defualt constructor
-        euni ( const RV rv ) :epdf ( rv ) {}
-        double eval ( const vec &val ) const  {return nk;}
-        double evallog ( const vec &val ) const  {return lnk;}
-        vec sample() const {
-                vec smp ( rv.count() );
-#pragma omp critical
-                UniRNG.sample_vector ( rv.count(),smp );
-                return low+elem_mult ( distance,smp );
-        }
-        //! set values of \c low and \c high
+        //! \name Constructors
+        //!@{
+        euni ( ) :epdf ( ) {}
+        euni ( const vec &low0, const vec &high0 ) {set_parameters ( low0,high0 );}
         void set_parameters ( const vec &low0, const vec &high0 ) {
                 distance = high0-low0;
@@ -422 +423 @@
                 nk = prod ( 1.0/distance );
                 lnk = log ( nk );
-        }
-        vec mean() const {return (high-low)/2.0;}
-        vec variance() const {return (pow(high,2)+pow(low,2)+elem_mult(high,low))/3.0;}
+                dim = low.length();
+        }
+        //!@}
+
+        double eval ( const vec &val ) const  {return nk;}
+        double evallog ( const vec &val ) const  {return lnk;}
+        vec sample() const {
+                vec smp ( dim );
+#pragma omp critical
+                UniRNG.sample_vector ( dim ,smp );
+                return low+elem_mult ( distance,smp );
+        }
+        //! set values of \c low and \c high
+        vec mean() const {return ( high-low ) /2.0;}
+        vec variance() const {return ( pow ( high,2 ) +pow ( low,2 ) +elem_mult ( high,low ) ) /3.0;}
 };
 
@@ -442 +455 @@
         vec& _mu; //cached epdf.mu;
 public:
-        //! Constructor
-        mlnorm ( const RV &rv, const RV &rvc );
+        //! \name Constructors
+        //!@{
+        mlnorm ( ) :mEF (),epdf ( ),A ( ),_mu ( epdf._mu() ) {ep =&epdf; };
+        mlnorm ( const  mat &A, const vec &mu0, const sq_T &R ) :epdf ( ),_mu ( epdf._mu() ) {
+                ep =&epdf; set_parameters ( A,mu0,R );
+        };
         //! Set \c A and \c R
         void set_parameters ( const  mat &A, const vec &mu0, const sq_T &R );
-//      //!Generate one sample of the posterior
-//      vec samplecond (const vec &cond, double &lik );
-//      //!Generate matrix of samples of the posterior
-//      mat samplecond (const vec &cond, vec &lik, int n );
+        //!@}
         //! Set value of \c rvc . Result of this operation is stored in \c epdf use function \c _ep to access it.
         void condition ( const vec &cond );
@@ -466 +480 @@
 /*! (Approximate) Student t density with linear function of mean value
 
-The internal epdf of this class is of the type of a Gaussian (enorm). 
+The internal epdf of this class is of the type of a Gaussian (enorm).
 However, each conditioning is trying to assure the best possible approximation by taking into account the zeta function. See [] for reference.
 
-Perhaps a moment-matching technique? 
+Perhaps a moment-matching technique?
 */
 class mlstudent : public mlnorm<ldmat> {
@@ -477 +491 @@
         ldmat Re;
 public:
-        mlstudent ( const RV &rv0, const RV &rvc0 ) :mlnorm<ldmat> ( rv0,rvc0 ),
-                        Lambda ( rv0.count() ),
-                        _R ( epdf._R() ) {}
-        void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0) {
-                epdf.set_parameters ( zeros ( rv.count() ),Lambda );
+        mlstudent ( ) :mlnorm<ldmat> (),
+                        Lambda (),      _R ( epdf._R() ) {}
+        void set_parameters ( const mat &A0, const vec &mu0, const ldmat &R0, const ldmat& Lambda0 ) {
+                it_assert_debug ( A0.rows() ==mu0.length(),"" );
+                it_assert_debug ( R0.rows() ==A0.rows(),"" );
+
+                epdf.set_parameters ( mu0,Lambda ); //
                 A = A0;
                 mu_const = mu0;
@@ -491 +507 @@
                 double zeta;
                 //ugly hack!
-                if ((cond.length()+1)==Lambda.rows()){
-                        zeta = Lambda.invqform ( concat(cond, vec_1(1.0)) );
-                } else {
+                if ( ( cond.length() +1 ) ==Lambda.rows() ) {
+                        zeta = Lambda.invqform ( concat ( cond, vec_1 ( 1.0 ) ) );
+                }
+                else {
                         zeta = Lambda.invqform ( cond );
                 }
                 _R = Re;
-                _R*=( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
+                _R*= ( 1+zeta );// / ( nu ); << nu is in Re!!!!!!
         };
 
@@ -517 +534 @@
         double k;
         //! cache of epdf.beta
-        vec* _beta;
+        vec &_beta;
 
 public:
         //! Constructor
-        mgamma ( const RV &rv,const RV &rvc ): mEF ( rv,rvc ), epdf ( rv ) {vec* tmp; epdf._param ( tmp,_beta );ep=&epdf;};
+        mgamma ( ) : mEF ( ), epdf (), _beta ( epdf._beta() ) {ep=&epdf;};
         //! Set value of \c k
-        void set_parameters ( double k );
-        void condition ( const vec &val ) {*_beta=k/val;};
+        void set_parameters ( double k, const vec &beta0 );
+        void condition ( const vec &val ) {_beta=k/val;};
 };
 
@@ -530 +547 @@
 \brief  Inverse-Gamma random walk
 
-Mean value, \f$\mu\f$, of this density is given by \c rvc .
-Standard deviation of the random walk is proportional to one \f$k\f$-th the mean.
-This is achieved by setting \f$\alpha=\mu/k+2\f$ and \f$\beta=\mu(\alpha-1)\f$.
-
-The standard deviation of the walk is then: \f$\mu/\sqrt(k)\f$.
+Mean value, \f$ \mu \f$, of this density is given by \c rvc .
+Standard deviation of the random walk is proportional to one \f$ k \f$-th the mean.
+This is achieved by setting \f$ \alpha=\mu/k^2+2 \f$ and \f$ \beta=\mu(\alpha-1)\f$.
+
+The standard deviation of the walk is then: \f$ \mu/\sqrt(k)\f$.
  */
 class migamma : public mEF {
-        protected:
+protected:
         //! Internal epdf that arise by conditioning on \c rvc
-                eigamma epdf;
+        eigamma epdf;
         //! Constant \f$k\f$
-                double k;
+        double k;
+        //! cache of epdf.alpha
+        vec &_alpha;
         //! cache of epdf.beta
-                vec* _beta;
-                //! chaceh of epdf.alpha
-                vec* _alpha;
-
-        public:
-        //! Constructor
-                migamma ( const RV &rv,const RV &rvc ): mEF ( rv,rvc ), epdf ( rv ) {epdf._param ( _alpha,_beta );ep=&epdf;};
+        vec &_beta;
+
+public:
+        //! \name Constructors
+        //!@{
+        migamma ( ) : mEF (), epdf ( ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
+        migamma ( const migamma &m ) : mEF (), epdf ( m.epdf ), _alpha ( epdf._alpha() ), _beta ( epdf._beta() ) {ep=&epdf;};
+        //!@}
+
         //! Set value of \c k
-                void set_parameters ( double k0 ){k=k0;*_alpha=1.0/(k*k)+2;};
-                void condition ( const vec &val ) {
-                        *_beta=elem_mult(val,(*_alpha-1));
-                };
+        void set_parameters ( int len, double k0 ) {
+                k=k0;
+                epdf.set_parameters ( ( 1.0/ ( k*k ) +2.0 ) *ones ( len ) /*alpha*/, ones ( len ) /*beta*/ );
+                dimc = dimension();
+        };
+        void condition ( const vec &val ) {
+                _beta=elem_mult ( val, ( _alpha-1.0 ) );
+        };
 };
 
@@ -576 +601 @@
 public:
         //! Constructor
-        mgamma_fix ( const RV &rv,const RV &rvc ) : mgamma ( rv,rvc ),refl ( rv.count() ) {};
+        mgamma_fix ( ) : mgamma ( ),refl () {};
         //! Set value of \c k
         void set_parameters ( double k0 , vec ref0, double l0 ) {
-                mgamma::set_parameters ( k0 );
+                mgamma::set_parameters ( k0, ref0 );
                 refl=pow ( ref0,1.0-l0 );l=l0;
+                dimc=dimension();
         };
 
-        void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); *_beta=k/mean;};
+        void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); _beta=k/mean;};
 };
 
@@ -600 +626 @@
  */
 class migamma_fix : public migamma {
-        protected:
+protected:
         //! parameter l
-                double l;
+        double l;
         //! reference vector
-                vec refl;
-        public:
+        vec refl;
+public:
         //! Constructor
-                migamma_fix ( const RV &rv,const RV &rvc ) : migamma ( rv,rvc ),refl ( rv.count() ) {};
+        migamma_fix ( ) : migamma (),refl ( ) {};
         //! Set value of \c k
-                void set_parameters ( double k0 , vec ref0, double l0 ) {
-                        migamma::set_parameters ( k0 );
-                        refl=pow ( ref0,1.0-l0 );l=l0;
-                };
-
-                void condition ( const vec &val ) {vec mean=elem_mult ( refl,pow ( val,l ) ); migamma::condition(mean);};
+        void set_parameters ( double k0 , vec ref0, double l0 ) {
+                migamma::set_parameters ( ref0.length(), k0 );
+                refl=pow ( ref0,1.0-l0 );
+                l=l0;
+                dimc = dimension();
+        };
+
+        void condition ( const vec &val ) {
+                vec mean=elem_mult ( refl,pow ( val,l ) );
+                migamma::condition ( mean );
+        };
 };
 //! Switch between various resampling methods.
@@ -632 +663 @@
         Array<vec> samples;
 public:
-        //! Default constructor
-        eEmp ( const RV &rv0 ,int n0 ) :epdf ( rv0 ),n ( n0 ),w ( n ),samples ( n ) {};
+        //! \name Constructors 
+        //!@{
+        eEmp ( ) :epdf ( ),w (  ),samples ( ) {};
+        eEmp (const eEmp &e ): epdf(e), w(e.w), samples(e.samples) {};
+        //!@}
+        
         //! Set samples and weights
         void set_parameters ( const vec &w0, const epdf* pdf0 );
@@ -639 +674 @@
         void set_samples ( const epdf* pdf0 );
         //! Set sample
-        void set_n ( int n0, bool copy=true ){w.set_size(n0,copy);samples.set_size(n0,copy);};
+        void set_n ( int n0, bool copy=true ) {w.set_size ( n0,copy );samples.set_size ( n0,copy );};
         //! Potentially dangerous, use with care.
         vec& _w()  {return w;};
@@ -655 +690 @@
         double evallog ( const vec &val ) const {it_error ( "Not implemented" );return 0.0;}
         vec mean() const {
-                vec pom=zeros ( rv.count() );
+                vec pom=zeros ( dim );
                 for ( int i=0;i<n;i++ ) {pom+=samples ( i ) *w ( i );}
                 return pom;
         }
         vec variance() const {
-                vec pom=zeros ( rv.count() );
-                for ( int i=0;i<n;i++ ) {pom+=pow(samples ( i ),2) *w ( i );}
-                return pom-pow(mean(),2);
+                vec pom=zeros ( dim );
+                for ( int i=0;i<n;i++ ) {pom+=pow ( samples ( i ),2 ) *w ( i );}
+                return pom-pow ( mean(),2 );
         }
 };
@@ -668 +703 @@
 
 ////////////////////////
-
-template<class sq_T>
-enorm<sq_T>::enorm ( const RV &rv ) :eEF ( rv ), mu ( rv.count() ),R ( rv.count() ),dim ( rv.count() ) {};
 
 template<class sq_T>
@@ -677 +709 @@
         mu = mu0;
         R = R0;
+        dim = mu0.length();
 };
 
@@ -692 +725 @@
 vec enorm<sq_T>::sample() const {
         vec x ( dim );
-        #pragma omp critical 
+#pragma omp critical
         NorRNG.sample_vector ( dim,x );
         vec smp = R.sqrt_mult ( x );
@@ -708 +741 @@
 
         for ( i=0;i<N;i++ ) {
-        #pragma omp critical 
+#pragma omp critical
                 NorRNG.sample_vector ( dim,x );
                 pom = R.sqrt_mult ( x );
@@ -741 +774 @@
 
 template<class sq_T>
-mlnorm<sq_T>::mlnorm ( const RV &rv0, const RV &rvc0 ) :mEF ( rv0,rvc0 ),epdf ( rv0 ),A ( rv0.count(),rv0.count() ),_mu ( epdf._mu() ) {
-        ep =&epdf;
-}
-
-template<class sq_T>
 void mlnorm<sq_T>::set_parameters ( const mat &A0, const vec &mu0, const sq_T &R0 ) {
-        epdf.set_parameters ( zeros ( rv.count() ),R0 );
+        it_assert_debug ( A0.rows() ==mu0.length(),"" );
+        it_assert_debug ( A0.rows() ==R0.rows(),"" );
+
+        epdf.set_parameters ( zeros ( A0.rows() ),R0 );
         A = A0;
         mu_const = mu0;
+        dimc=A0.cols();
 }
 
@@ -785 +817 @@
 template<class sq_T>
 epdf* enorm<sq_T>::marginal ( const RV &rvn ) const {
+        it_assert_debug(isnamed(), "rv description is not assigned");
         ivec irvn = rvn.dataind ( rv );
 
-        sq_T Rn ( R,irvn );
-        enorm<sq_T>* tmp = new enorm<sq_T> ( rvn );
+        sq_T Rn ( R,irvn ); //select rows and columns of R
+        
+        enorm<sq_T>* tmp = new enorm<sq_T>;
+        tmp->set_rv ( rvn );
         tmp->set_parameters ( mu ( irvn ), Rn );
         return tmp;
@@ -796 +831 @@
 mpdf* enorm<sq_T>::condition ( const RV &rvn ) const {
 
+        it_assert_debug ( isnamed(),"rvs are not assigned" );
+
         RV rvc = rv.subt ( rvn );
-        it_assert_debug ( ( rvc.count() +rvn.count() ==rv.count() ),"wrong rvn" );
+        it_assert_debug ( ( rvc._dsize() +rvn._dsize() ==rv._dsize() ),"wrong rvn" );
         //Permutation vector of the new R
         ivec irvn = rvn.dataind ( rv );
@@ -807 +844 @@
         mat S=Rn.to_mat();
         //fixme
-        int n=rvn.count()-1;
+        int n=rvn._dsize()-1;
         int end=R.rows()-1;
         mat S11 = S.get ( 0,n, 0, n );
-        mat S12 = S.get ( 0, n , rvn.count(), end );
-        mat S22 = S.get ( rvn.count(), end, rvn.count(), end );
+        mat S12 = S.get ( 0, n , rvn._dsize(), end );
+        mat S22 = S.get ( rvn._dsize(), end, rvn._dsize(), end );
 
         vec mu1 = mu ( irvn );
@@ -818 +855 @@
         sq_T R_n ( S11 - A *S12.T() );
 
-        mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( rvn,rvc );
-
+        mlnorm<sq_T>* tmp=new mlnorm<sq_T> ( );
+        tmp->set_rv(rvn); tmp->set_rvc(rvc);
         tmp->set_parameters ( A,mu1-A*mu2,R_n );
         return tmp;

bdm/stat/libFN.cpp

@@ -5 +5 @@
 using namespace bdm;
 
-bilinfn::bilinfn ( const RV &rvx0, const RV &rvu0, const mat &A0, const mat &B0 ) : diffbifn (A0.rows(), rvx0,rvu0 )
-{
-        //check input
-        it_assert_debug ( ( A0.cols() ==dimx ) & ( A0.rows() ==B0.rows() ), "linfn:: wrong A" );
-        it_assert_debug ( ( B0.cols() ==dimu ), "linfn:: wrong B" );
-
-        // set dimensions
-        dimy = A0.rows();
-
-        //set internals
-        A = A0;
-        B = B0;
-};
 
 inline vec bilinfn::eval ( const  vec &x0, const vec &u0 )

bdm/stat/libFN.h

@@ -28 +28 @@
                 vec eval ( const vec &cond ) {return val;};
                 //!Default constructor
-                constfn ( const vec &val0 ) :fnc(val0.length()), val ( val0 ) {};
+                constfn ( const vec &val0 ) :fnc(), val ( val0 ) {dimy=val.length();};
 };
 
@@ -38 +38 @@
                 //! Matrix A
                 mat A;
-                //! Matrix B
+                //! vector B
                 vec B;
         public :
-                vec eval (const vec &cond ) {it_assert_debug ( cond.length() ==rv.count(), "linfn::eval Wrong cond." );return A*cond+B;};
+                vec eval (const vec &cond ) {it_assert_debug ( cond.length() ==A.cols(), "linfn::eval Wrong cond." );return A*cond+B;};
 
 //              linfn evalsome ( ivec &rvind );
                 //!default constructor
-                linfn ( const RV &rv0 ) : fnc(rv0.count()), rv ( rv0 ),A ( eye ( rv0.count() ) ),B ( zeros ( rv0.count() ) ) { };
+                linfn ( ) : fnc(), A ( ),B () { };
                 //! Set values of \c A and \c B
-                void set_parameters ( const mat &A0 , const vec &B0 ) {A=A0; B=B0;};
+                void set_parameters ( const mat &A0 , const vec &B0 ) {A=A0; B=B0; dimy=A.rows();};
 };
 
@@ -86 +86 @@
                 virtual void dfdu_cond ( const vec &x0, const vec &u0, mat &A, bool full=true ) {};
                 //!Default constructor (dimy is not set!)
-                diffbifn (int dimy, const RV rvx0, const RV rvu0 ) : fnc(dimy), rvx ( rvx0 ),rvu ( rvu0 ) {dimx=rvx.count();dimu=rvu.count();};
+                diffbifn () : fnc() {};
                 //! access function
                 int _dimx() const{return dimx;}
@@ -100 +100 @@
                 mat B;
         public :
+                //!\name Constructors
+                //!@{
+                
+                bilinfn () : diffbifn () ,A() ,B()      {};
+                bilinfn (const mat A0, const mat B0) {set_parameters(A0,B0);};
+                //! Alternative constructor
+                void set_parameters(const mat A0, const mat B0){
+                        it_assert_debug(A0.rows()==B0.rows(),"");
+                        A=A0;B=B0;
+                        dimy=A.rows();
+                        dimx=A.cols();
+                        dimu=B.cols();
+                }
+                //!@}
+                
+                //!\name Mathematical operations
+                //!@{
                 vec eval ( const  vec &x0, const vec &u0 );
-
-                //! Default constructor
-                bilinfn ( const RV &rvx0, const RV &rvu0 ) : diffbifn (dimx, rvx0,rvu0 ) ,A ( eye ( dimx ) ),B ( zeros ( dimx,dimu ) )  {};
-                //! Alternative constructor
-                bilinfn ( const RV &rvx0, const RV &rvu0, const mat &A0, const mat &B0 );
-                //!
                 void dfdx_cond ( const vec &x0, const vec &u0, mat &F, bool full )
                 {
@@ -118 +129 @@
                         if ( full ) F=B;        //else : nothing has changed no need to regenerate
                 }
+                //!@}
 };
 

bdm/stat/loggers.cpp

@@ -27 +27 @@
         int i,j,k;
         int nsc=0;
-        for ( i=0;i<entries.length();i++ ) {nsc+=entries ( i ).count();}
+        for ( i=0;i<entries.length();i++ ) {nsc+=entries ( i )._dsize();}
         ; //all entries!!
 

bdm/stat/loggers.h

@@ -42 +42 @@
                 int i; int n =entries.length();
                 vectors.set_size ( n ); 
-                for ( i=0;i<n;i++ ) {vectors(i).set_size (maxlen,entries(i).count() );}
+                for ( i=0;i<n;i++ ) {vectors(i).set_size (maxlen,entries(i)._dsize() );}
                 ;
         }