Changeset 886

Timestamp:

03/29/10 23:01:49 (15 years ago)

Author:

smidl

Message:

epdf and emix now have _base classes

Location:

library

Files:

• bdm/estim/mixtures.cpp (modified) (1 diff)
• bdm/estim/mixtures.h (modified) (6 diffs)
• bdm/stat/emix.cpp (modified) (11 diffs)
• bdm/stat/emix.h (modified) (6 diffs)
• tests/stresssuite/merger_2d_stress.cpp (modified) (1 diff)
• tests/stresssuite/mixtures_stress.cpp (modified) (1 diff)
• tests/testsuite/epdf_test.cpp (modified) (2 diffs)

library/bdm/estim/mixtures.cpp

@@ -35 +35 @@
         }
 
-        //est already exists - must be deleted before build_est() can be used
-        delete est;
-        build_est();
 }
 

library/bdm/estim/mixtures.h

@@ -48 +48 @@
         //! Statistics for weights
         multiBM weights;
+        //aux
+        friend class eprod_mix;
         //!Posterior on component parameters
-        eprod* est;
+        class eprod_mix: public eprod_base {
+                protected:
+                const MixEF &mix; // pointer to parent.n
+                public:
+                eprod_mix(const MixEF &m):mix(m){}
+                const epdf* factor(int i) const {return (i==(mix.n-1)) ? &mix.weights.posterior() : &mix.Coms(i)->posterior();}
+                const int no_factors()const {return mix.n+1;}
+        } est;
         ////!Indeces of component rvc in common rvc
 
@@ -55 +64 @@
         MixEF_METHOD method;
 
-        //! Auxiliary function for use in constructors
-        void build_est() {
-                est = new eprod;
-                if ( n > 0 ) {
-                        Array<const epdf*> epdfs ( n + 1 );
-                        for ( int i = 0; i < Coms.length(); i++ ) {
-                                epdfs ( i ) = & ( Coms ( i )->posterior() );
-                        }
-                        // last in the product is the weight
-                        epdfs ( n ) = & ( weights.posterior() );
-                        est->set_parameters ( epdfs, false );
-                }
-        }
-
 public:
         //! Full constructor
         MixEF ( const Array<BMEF*> &Coms0, const vec &alpha0 ) :
                         BMEF ( ), n ( Coms0.length() ), Coms ( n ),
-                        weights (), method ( QB ) {
+                        weights (), est(*this), method ( QB ) {
                 for ( int i = 0; i < n; i++ ) {
                         Coms ( i ) = ( BMEF* ) Coms0 ( i )->_copy();
@@ -79 +74 @@
                 weights.set_parameters(alpha0);
                 weights.validate();
-                build_est();
         }
 
@@ -85 +79 @@
         MixEF () :
                         BMEF ( ), n ( 0 ), Coms ( 0 ),
-                        weights (), method ( QB ) {
-                build_est();
+                        weights (), est(*this), method ( QB ) {
         }
         //! Copy constructor
         MixEF ( const MixEF &M2 ) : BMEF ( ), n ( M2.n ), Coms ( n ),
-                        weights ( M2.weights ), method ( M2.method ) {
+                        weights ( M2.weights ), est(*this), method ( M2.method ) {
                 for ( int i = 0; i < n; i++ ) {
                         Coms ( i ) = (BMEF*) M2.Coms ( i )->_copy();
                 }
-                build_est();
         }
 
@@ -103 +95 @@
         void init ( BMEF* Com0, const mat &Data, const int c = 5 );
         //Destructor
-        ~MixEF() {
-                delete est;
-                for ( int i = 0; i < n; i++ ) {
-                        delete Coms ( i );
-                }
-        }
         //! Recursive EM-like algorithm (QB-variant), see Karny et. al, 2006
         void bayes ( const vec &yt, const vec &cond );
@@ -117 +103 @@
         double logpred ( const vec &yt ) const;
         //! return correctly typed posterior (covariant return)
-        const eprod& posterior() const {
-                return *est;
+        const eprod_mix& posterior() const {
+                return est;
         }
 

library/bdm/stat/emix.cpp

@@ -3 +3 @@
 namespace bdm {
 
-void emix::validate (){
-        bdm_assert ( Coms.length() > 0, "There has to be at least one component." );
-
-        bdm_assert ( Coms.length() == w.length(), "It is obligatory to define weights of all the components." );
+void emix_base::validate (){
+        bdm_assert ( no_coms() > 0, "There has to be at least one component." );
+
+        bdm_assert ( no_coms() == w.length(), "It is obligatory to define weights of all the components." );
 
         double sum_w = sum ( w );
@@ -12 +12 @@
         w = w / sum_w;
 
-        dim = Coms ( 0 )->dimension();
-        RV rv_tmp = Coms ( 0 )->_rv() ;
-        bool isnamed = Coms ( 0 )->isnamed();
-        for ( int i = 1; i < Coms.length(); i++ ) {
-                bdm_assert ( dim == ( Coms ( i )->dimension() ), "Component sizes do not match!" );
-                isnamed &= Coms(i)->isnamed() & Coms(i)->_rv().equal(rv_tmp);
+        dim = component ( 0 )->dimension();
+        RV rv_tmp = component ( 0 )->_rv() ;
+        bool isnamed = component( 0 )->isnamed();
+        for ( int i = 1; i < no_coms(); i++ ) {
+                bdm_assert ( dim == ( component ( i )->dimension() ), "Component sizes do not match!" );
+                isnamed &= component(i)->isnamed() & component(i)->_rv().equal(rv_tmp);
         }
         if (isnamed)
@@ -23 +23 @@
 }
 
-void emix::from_setting ( const Setting &set ) {
-        UI::get ( Coms, set, "pdfs", UI::compulsory );
-
-        if ( !UI::get ( w, set, "weights", UI::optional ) ) {
-                int len = Coms.length();
-                w.set_length ( len );
-                w = 1.0 / len;
-        }
-}
-
-
-vec emix::sample() const {
+
+
+vec emix_base::sample() const {
         //Sample which component
         vec cumDist = cumsum ( w );
@@ -46 +37 @@
         }
 
-        return Coms ( i )->sample();
-}
-
-vec emix::mean() const {
+        return component ( i )->sample();
+}
+
+vec emix_base::mean() const {
         int i;
         vec mu = zeros ( dim );
         for ( i = 0; i < w.length(); i++ ) {
-                mu += w ( i ) * Coms ( i )->mean();
+                mu += w ( i ) * component ( i )->mean();
         }
         return mu;
 }
 
-vec emix::variance() const {
+vec emix_base::variance() const {
         //non-central moment
         vec mom2 = zeros ( dim );
         for ( int i = 0; i < w.length(); i++ ) {
-                mom2 += w ( i ) * ( Coms ( i )->variance() + pow ( Coms ( i )->mean(), 2 ) );
+                mom2 += w ( i ) * ( component( i )->variance() + pow ( component ( i )->mean(), 2 ) );
         }
         //central moment
@@ -68 +59 @@
 }
 
-double emix::evallog ( const vec &val ) const {
+double emix_base::evallog ( const vec &val ) const {
         int i;
         double sum = 0.0;
         for ( i = 0; i < w.length(); i++ ) {
-                sum += w ( i ) * exp ( Coms ( i )->evallog ( val ) );
+                sum += w ( i ) * exp ( component ( i )->evallog ( val ) );
         }
         if ( sum == 0.0 ) {
@@ -82 +73 @@
 }
 
-vec emix::evallog_mat ( const mat &Val ) const {
+vec emix_base::evallog_mat ( const mat &Val ) const {
         vec x = zeros ( Val.cols() );
         for ( int i = 0; i < w.length(); i++ ) {
-                x += w ( i ) * exp ( Coms ( i )->evallog_mat ( Val ) );
+                x += w ( i ) * exp ( component( i )->evallog_mat ( Val ) );
         }
         return log ( x );
 };
 
-mat emix::evallog_coms ( const mat &Val ) const {
+mat emix_base::evallog_coms ( const mat &Val ) const {
         mat X ( w.length(), Val.cols() );
         for ( int i = 0; i < w.length(); i++ ) {
-                X.set_row ( i, w ( i ) *exp ( Coms ( i )->evallog_mat ( Val ) ) );
+                X.set_row ( i, w ( i ) *exp ( component( i )->evallog_mat ( Val ) ) );
         }
         return X;
 }
 
-shared_ptr<epdf> emix::marginal ( const RV &rv ) const {
+shared_ptr<epdf> emix_base::marginal ( const RV &rv ) const {
         emix *tmp = new emix();
         shared_ptr<epdf> narrow ( tmp );
@@ -105 +96 @@
 }
 
-void emix::marginal ( const RV &rv, emix &target ) const {
+void emix_base::marginal ( const RV &rv, emix &target ) const {
         bdm_assert ( isnamed(), "rvs are not assigned" );
 
-        Array<shared_ptr<epdf> > Cn ( Coms.length() );
-        for ( int i = 0; i < Coms.length(); i++ ) {
-                Cn ( i ) = Coms ( i )->marginal ( rv );
+        Array<shared_ptr<epdf> > Cn ( no_coms() );
+        for ( int i = 0; i < no_coms(); i++ ) {
+                Cn ( i ) = component ( i )->marginal ( rv );
         }
 
@@ -118 +109 @@
 }
 
-shared_ptr<pdf> emix::condition ( const RV &rv ) const {
+shared_ptr<pdf> emix_base::condition ( const RV &rv ) const {
         bdm_assert ( isnamed(), "rvs are not assigned" );
         mratio *tmp = new mratio ( this, rv );
@@ -124 +115 @@
 }
 
-void egiwmix::set_parameters ( const vec &w0, const Array<egiw*> &Coms0, bool copy ) {
-        w = w0 / sum ( w0 );
-        int i;
-        for ( i = 0; i < w.length(); i++ ) {
-                bdm_assert_debug ( dim == ( Coms0 ( i )->dimension() ), "Component sizes do not match!" );
-        }
-        if ( copy ) {
-                Coms.set_length ( Coms0.length() );
-                for ( i = 0; i < w.length(); i++ ) {
-                        bdm_error ( "Not implemented" );
-                        // *Coms ( i ) = *Coms0 ( i );
-                }
-                destroyComs = true;
-        } else {
-                Coms = Coms0;
-                destroyComs = false;
-        }
-}
-
-void    egiwmix::validate (){
-   dim = Coms ( 0 )->dimension();
-}
-
-vec egiwmix::sample() const {
-        //Sample which component
-        vec cumDist = cumsum ( w );
-        double u0;
-#pragma omp critical
-        u0 = UniRNG.sample();
-
-        int i = 0;
-        while ( ( cumDist ( i ) < u0 ) && ( i < ( w.length() - 1 ) ) ) {
-                i++;
-        }
-
-        return Coms ( i )->sample();
-}
-
-vec egiwmix::mean() const {
-        int i;
-        vec mu = zeros ( dim );
-        for ( i = 0; i < w.length(); i++ ) {
-                mu += w ( i ) * Coms ( i )->mean();
-        }
-        return mu;
-}
-
-vec egiwmix::variance() const {
-        // non-central moment
-        vec mom2 = zeros ( dim );
-        for ( int i = 0; i < w.length(); i++ ) {
-                // pow is overloaded, we have to use another approach
-                mom2 += w ( i ) * ( Coms ( i )->variance() + elem_mult ( Coms ( i )->mean(), Coms ( i )->mean() ) );
-        }
-        // central moment
-        // pow is overloaded, we have to use another approach
-        return mom2 - elem_mult ( mean(), mean() );
-}
-
-shared_ptr<epdf> egiwmix::marginal ( const RV &rv ) const {
-        emix *tmp = new emix();
-        shared_ptr<epdf> narrow ( tmp );
-        marginal ( rv, *tmp );
-        return narrow;
-}
-
-void egiwmix::marginal ( const RV &rv, emix &target ) const {
-        bdm_assert_debug ( isnamed(), "rvs are not assigned" );
-
-        Array<shared_ptr<epdf> > Cn ( Coms.length() );
-        for ( int i = 0; i < Coms.length(); i++ ) {
-                Cn ( i ) = Coms ( i )->marginal ( rv );
-        }
-
-        target._w() = w;
-        target._Coms() = Cn;
-        target.validate();
-}
-
-egiw*   egiwmix::approx() {
-        // NB: dimx == 1 !!!
-        // The following code might look a bit spaghetti-like,
-        // consult Dedecius, K. et al.: Partial forgetting in AR models.
-
-        double sumVecCommon;                            // common part for many terms in eq.
-        int len = w.length();                           // no. of mix components
-        int dimLS = Coms ( 1 )->_V()._D().length() - 1;         // dim of LS
-        vec vecNu ( len );                                      // vector of dfms of components
-        vec vecD ( len );                                       // vector of LS reminders of comps.
-        vec vecCommon ( len );                          // vector of common parts
-        mat matVecsTheta;                               // matrix which rows are theta vects.
-
-        // fill in the vectors vecNu, vecD and matVecsTheta
-        for ( int i = 0; i < len; i++ ) {
-                vecNu.shift_left ( Coms ( i )->_nu() );
-                vecD.shift_left ( Coms ( i )->_V()._D() ( 0 ) );
-                matVecsTheta.append_row ( Coms ( i )->est_theta() );
-        }
-
-        // calculate the common parts and their sum
-        vecCommon = elem_mult ( w, elem_div ( vecNu, vecD ) );
-        sumVecCommon = sum ( vecCommon );
-
-        // LS estimator of theta
-        vec aprEstTheta ( dimLS );
-        aprEstTheta.zeros();
-        for ( int i = 0; i < len; i++ ) {
-                aprEstTheta +=  matVecsTheta.get_row ( i ) * vecCommon ( i );
-        }
-        aprEstTheta /= sumVecCommon;
-
-
-        // LS estimator of dfm
-        double aprNu;
-        double A = log ( sumVecCommon );                // Term 'A' in equation
-
-        for ( int i = 0; i < len; i++ ) {
-                A += w ( i ) * ( log ( vecD ( i ) ) - psi ( 0.5 * vecNu ( i ) ) );
-        }
-
-        aprNu = ( 1 + sqrt ( 1 + 2 * ( A - LOG2 ) / 3 ) ) / ( 2 * ( A - LOG2 ) );
-
-
-        // LS reminder (term D(0,0) in C-syntax)
-        double aprD = aprNu / sumVecCommon;
-
-        // Aproximation of cov
-        // the following code is very numerically sensitive, thus
-        // we have to eliminate decompositions etc. as much as possible
-        mat aprC = zeros ( dimLS, dimLS );
-        for ( int i = 0; i < len; i++ ) {
-                aprC += Coms ( i )->est_theta_cov().to_mat() * w ( i );
-                vec tmp = ( matVecsTheta.get_row ( i ) - aprEstTheta );
-                aprC += vecCommon ( i ) * outer_product ( tmp, tmp );
-        }
-
-        // Construct GiW pdf :: BEGIN
-        ldmat aprCinv ( inv ( aprC ) );
-        vec D = concat ( aprD, aprCinv._D() );
-        mat L = eye ( dimLS + 1 );
-        L.set_submatrix ( 1, 0, aprCinv._L() * aprEstTheta );
-        L.set_submatrix ( 1, 1, aprCinv._L() );
-        ldmat aprLD ( L, D );
-
-        egiw* aprgiw = new egiw ( 1, aprLD, aprNu );
-        return aprgiw;
-};
+void emix::from_setting ( const Setting &set ) {
+        UI::get ( Coms, set, "pdfs", UI::compulsory );
+        UI::get ( w, set, "weights", UI::compulsory );
+}
+
+void    emix::validate (){
+        emix_base::validate();
+        dim = Coms ( 0 )->dimension();
+}
+
 
 double mprod::evallogcond ( const vec &val, const vec &cond ) {
@@ -378 +232 @@
 }
 
-vec eprod::mean() const {
+vec eprod_base::mean() const {
         vec tmp ( dim );
-        for ( int i = 0; i < epdfs.length(); i++ ) {
-                vec pom = epdfs ( i )->mean();
+        for ( int i = 0; i < no_factors(); i++ ) {
+                vec pom = factor( i )->mean();
                 dls ( i )->pushup ( tmp, pom );
         }
@@ -387 +241 @@
 }
 
-vec eprod::variance() const {
+vec eprod_base::variance() const {
         vec tmp ( dim ); //second moment
-        for ( int i = 0; i < epdfs.length(); i++ ) {
-                vec pom = epdfs ( i )->mean();
+        for ( int i = 0; i < no_factors(); i++ ) {
+                vec pom = factor ( i )->mean();
                 dls ( i )->pushup ( tmp, pow ( pom, 2 ) );
         }
         return tmp - pow ( mean(), 2 );
 }
-vec eprod::sample() const {
+vec eprod_base::sample() const {
         vec tmp ( dim );
-        for ( int i = 0; i < epdfs.length(); i++ ) {
-                vec pom = epdfs ( i )->sample();
+        for ( int i = 0; i < no_factors(); i++ ) {
+                vec pom = factor ( i )->sample();
                 dls ( i )->pushup ( tmp, pom );
         }
         return tmp;
 }
-double eprod::evallog ( const vec &val ) const {
+double eprod_base::evallog ( const vec &val ) const {
         double tmp = 0;
-        for ( int i = 0; i < epdfs.length(); i++ ) {
-                tmp += epdfs ( i )->evallog ( dls ( i )->pushdown ( val ) );
-        }
-        bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
-        return tmp;
-}
-
-}
-// mprod::mprod ( Array<pdf*> mFacs, bool overlap) : pdf ( RV(), RV() ), n ( mFacs.length() ), epdfs ( n ), pdfs ( mFacs ), rvinds ( n ), rvcinrv ( n ), irvcs_rvc ( n ) {
-//              int i;
-//              bool rvaddok;
-//              // Create rv
-//              for ( i = 0;i < n;i++ ) {
-//                      rvaddok=rv.add ( pdfs ( i )->_rv() ); //add rv to common rvs.
-//                      // If rvaddok==false, pdfs overlap => assert error.
-//                      epdfs ( i ) = & ( pdfs ( i )->posterior() ); // add pointer to epdf
-//              };
-//              // Create rvc
-//              for ( i = 0;i < n;i++ ) {
-//                      rvc.add ( pdfs ( i )->_rvc().subt ( rv ) ); //add rv to common rvs.
-//              };
-//
-// //           independent = true;
-//              //test rvc of pdfs and fill rvinds
-//              for ( i = 0;i < n;i++ ) {
-//                      // find ith rv in common rv
-//                      rvsinrv ( i ) = pdfs ( i )->_rv().dataind ( rv );
-//                      // find ith rvc in common rv
-//                      rvcinrv ( i ) = pdfs ( i )->_rvc().dataind ( rv );
-//                      // find ith rvc in common rv
-//                      irvcs_rvc ( i ) = pdfs ( i )->_rvc().dataind ( rvc );
-//                      //
-// /*                   if ( rvcinrv ( i ).length() >0 ) {independent = false;}
-//                      if ( irvcs_rvc ( i ).length() >0 ) {independent = false;}*/
-//              }
-//      };
-
+        for ( int i = 0; i < no_factors(); i++ ) {
+                tmp += factor ( i )->evallog ( dls ( i )->pushdown ( val ) );
+        }
+        //bdm_assert_debug ( std::isfinite ( tmp ), "Infinite" );
+        return tmp;
+}
+
+}
+

library/bdm/stat/emix.h

@@ -97 +97 @@
 
 
-
-
 private:
         // not implemented
         mratio ( const mratio & );
         mratio &operator= ( const mratio & );
+};
+
+class emix; //forward
+
+/*! Base class (Interface) for mixtures 
+*/
+class emix_base : public epdf {
+        protected:
+        //! reference to vector of weights
+        vec &w;
+        //! function returning ith component
+        virtual const epdf * component(const int &i) const=0;
+        
+        virtual int no_coms() const=0;
+        
+        public:
+                
+                emix_base(vec &w0): w(w0){}
+                
+                void validate ();
+                
+                vec sample() const;
+                
+                vec mean() const;
+                
+                vec variance() const;
+                
+                double evallog ( const vec &val ) const;
+                
+                vec evallog_mat ( const mat &Val ) const;
+                
+                //! Auxiliary function that returns pdflog for each component
+                mat evallog_coms ( const mat &Val ) const;
+                
+                shared_ptr<epdf> marginal ( const RV &rv ) const;
+                //! Update already existing marginal density  \c target
+                void marginal ( const RV &rv, emix &target ) const;
+                shared_ptr<pdf> condition ( const RV &rv ) const;
+                
+                //Access methods        
+                //! returns a reference to the internal weights. Use with Care!
+                vec& _w() {
+                        return w;
+                }
+                
+                const vec& _w() const {
+                        return w;
+                }
+                //!access
+                const epdf* _com(int i) const {return component(i);}
+                
 };
 
@@ -115 +164 @@
 
 */
-class emix : public epdf {
+class emix : public emix_base {
 protected:
         //! weights of the components
-        vec w;
+        vec weights;
 
         //! Component (epdfs)
@@ -125 +174 @@
 public:
         //! Default constructor
-        emix ( ) : epdf ( ) { }
-
-        virtual void validate ();
-
-        vec sample() const;
-
-        vec mean() const;
-
-        vec variance() const;
-
-        double evallog ( const vec &val ) const;
-
-        vec evallog_mat ( const mat &Val ) const;
-
-        //! Auxiliary function that returns pdflog for each component
-        mat evallog_coms ( const mat &Val ) const;
-
-        shared_ptr<epdf> marginal ( const RV &rv ) const;
-        //! Update already existing marginal density  \c target
-        void marginal ( const RV &rv, emix &target ) const;
-        shared_ptr<pdf> condition ( const RV &rv ) const;
-
-        //Access methods        
-        //! returns a reference to the internal weights. Use with Care!
-        vec& _w() {
-                return w;
-        }
-
-        /*!
-          \brief returns a reference to the internal array of components. Use with Care! Set components \c Coms
-
-          Shared pointers in Coms are kept inside this instance and
-          shouldn't be modified after being passed to this method.
-        */
-        Array<shared_ptr<epdf> >& _Coms ( ) {
-                return Coms;
-        }
-
-        //! returns a reference to the internal components specified by index i. Use with Care!
-        shared_ptr<epdf> _Coms ( int i ) {
-                return Coms ( i );
-        }
-
-        void set_rv ( const RV &rv ) {
-                epdf::set_rv ( rv );
-                for ( int i = 0; i < Coms.length(); i++ ) {
-                        Coms ( i )->set_rv ( rv );
-                }
-        }
+        emix ( ) : emix_base ( weights) { }
+        
+        const epdf* component(const int &i) const {return Coms(i).get();}
+        void validate();
+        
+
+        int no_coms() const {return Coms.length(); }
 
         //! Load from structure with elements:
@@ -184 +191 @@
         //!@}
         void from_setting ( const Setting &set );
+
+        void set_rv ( const RV &rv ) {
+                epdf::set_rv ( rv );
+                for ( int i = 0; i < no_coms(); i++ ) {
+                        Coms( i )->set_rv ( rv );
+                }
+        }
+        
+        Array<shared_ptr<epdf> >& _Coms ( ) {
+                        return Coms;
+                }
 };
 SHAREDPTR ( emix );
 UIREGISTER ( emix );
 
-/*!
-* \brief Mixture of egiws
-
-*/
-class egiwmix : public egiw {
-protected:
-        //! weights of the components
-        vec w;
-        //! Component (epdfs)
-        Array<egiw*> Coms;
-        //!Flag if owning Coms
-        bool destroyComs;
-public:
-        //!Default constructor
-        egiwmix ( ) : egiw ( ) {};
-
-        //! Set weights \c w and components \c Coms
-        //!By default Coms are copied inside. Parameter \c copy can be set to false if Coms live externally. Use method ownComs() if Coms should be destroyed by the destructor.
-        void set_parameters ( const vec &w, const Array<egiw*> &Coms, bool copy = false );
-        
-        //!return expected value
-        void validate();
-        
-        vec mean() const;
-
-        //!return a sample from the density
-        vec sample() const;
-
-        //!return the expected variance
-        vec variance() const;
-
-        // TODO!!! Defined to follow ANSI and/or for future development
-        void mean_mat ( mat &M, mat&R ) const {};
-        double evallog_nn ( const vec &val ) const {
-                return 0;
-        };
-        double lognc () const {
-                return 0;
-        }
-
-        shared_ptr<epdf> marginal ( const RV &rv ) const;
-        //! marginal density update
-        void marginal ( const RV &rv, emix &target ) const;
-
-//Access methods
-        //! returns a pointer to the internal mean value. Use with Care!
-        vec& _w() {
-                return w;
-        }
-        virtual ~egiwmix() {
-                if ( destroyComs ) {
-                        for ( int i = 0; i < Coms.length(); i++ ) {
-                                delete Coms ( i );
-                        }
-                }
-        }
-        //! Auxiliary function for taking ownership of the Coms()
-        void ownComs() {
-                destroyComs = true;
-        }
-
-        //!access function
-        egiw* _Coms ( int i ) {
-                return Coms ( i );
-        }
-
-        void set_rv ( const RV &rv ) {
-                egiw::set_rv ( rv );
-                for ( int i = 0; i < Coms.length(); i++ ) {
-                        Coms ( i )->set_rv ( rv );
-                }
-        }
-
-        //! Approximation of a GiW mix by a single GiW pdf
-        egiw* approx();
-};
 
 /*! \brief Chain rule decomposition of epdf
@@ -329 +271 @@
 SHAREDPTR ( mprod );
 
+
 //! Product of independent epdfs. For dependent pdfs, use mprod.
-class eprod: public epdf {
+class eprod_base: public epdf {
 protected:
-        //! Components (epdfs)
-        Array<const epdf*> epdfs;
         //! Array of indeces
         Array<datalink*> dls;
+        //! interface for a factor
+        virtual const epdf* factor(int i) const NOT_IMPLEMENTED(NULL); 
+        //!number of factors
+        virtual const int no_factors() const =0;
 public:
         //! Default constructor
-        eprod () : epdfs ( 0 ), dls ( 0 ) {};
+        eprod_base () :  dls (0) {};
         //! Set internal
-        void set_parameters ( const Array<const epdf*> &epdfs0, bool named = true ) {
-                epdfs = epdfs0;//.set_length ( epdfs0.length() );
-                dls.set_length ( epdfs.length() );
+        vec mean() const;
+
+        vec variance() const;
+
+        vec sample() const;
+
+        double evallog ( const vec &val ) const;
+
+        //!Destructor
+        ~eprod_base() {
+                for ( int i = 0; i < dls.length(); i++ ) {
+                        delete dls ( i );
+                }
+        }
+        void validate() {
+                dls.set_length ( no_factors() );
+                
                 bool independent = true;
-                if ( named ) {
-                        for ( int i = 0; i < epdfs.length(); i++ ) {
-                                independent = rv.add ( epdfs ( i )->_rv() );
-                                bdm_assert_debug ( independent, "eprod:: given components are not independent." );
-                        }
-                        dim = rv._dsize();
-                        
-                } else {
-                        dim = 0;
-                        for ( int i = 0; i < epdfs.length(); i++ ) {
-                                dim += epdfs ( i )->dimension();
-                        }
-                }
+                dim = 0;
+                for ( int i = 0; i < no_factors(); i++ ) {
+                        independent = rv.add ( factor ( i )->_rv() );
+                        dim += factor ( i )->dimension();
+                        bdm_assert_debug ( independent, "eprod:: given components are not independent." );
+                };
+                
                 //
                 int cumdim = 0;
                 int dimi = 0;
                 int i;
-                for ( i = 0; i < epdfs.length(); i++ ) {
+                for ( i = 0; i < no_factors(); i++ ) {
                         dls ( i ) = new datalink;
-                        if ( named ) {
-                                dls ( i )->set_connection ( epdfs ( i )->_rv() , rv );
-                        } else {
-                                dimi = epdfs ( i )->dimension();
+                        if ( isnamed() ) { // rvs are complete
+                                dls ( i )->set_connection ( factor ( i )->_rv() , rv );
+                        } else { //rvs are not reliable
+                                dimi = factor ( i )->dimension();
                                 dls ( i )->set_connection ( dimi, dim, linspace ( cumdim, cumdim + dimi - 1 ) );
                                 cumdim += dimi;
@@ -373 +326 @@
         
         }
-
-        
-        vec mean() const;
-
-        vec variance() const;
-
-        vec sample() const;
-
-        double evallog ( const vec &val ) const;
-
-        //!access function
-        const epdf* operator () ( int i ) const {
-                bdm_assert_debug ( i < epdfs.length(), "wrong index" );
-                return epdfs ( i );
-        }
-
-        //!Destructor
-        ~eprod() {
-                for ( int i = 0; i < epdfs.length(); i++ ) {
-                        delete dls ( i );
-                }
-        }
-};
-
+};
+
+class eprod: public eprod_base{
+        protected:
+                Array<shared_ptr<epdf> > factors;
+        public:
+                const epdf* factor(int i) const {return factors(i).get();}
+                const int no_factors() const {return factors.length();}
+        void set_parameters ( const Array<shared_ptr<epdf> > &epdfs0) {
+                factors = epdfs0;
+        }
+};
+
+//! similar to eprod but used only internally -- factors are external pointers
+class eprod_internal: public eprod_base{
+        protected:
+                Array<epdf* > factors;
+                const epdf* factor(int i) const {return factors(i);}
+                const int no_factors() const {return factors.length();}
+        public:
+                void set_parameters ( const Array<epdf *> &epdfs0) {
+                        factors = epdfs0;
+                }
+};
 
 /*! \brief Mixture of pdfs with constant weights, all pdfs are of equal RV and RVC

library/tests/stresssuite/merger_2d_stress.cpp

@@ -67 +67 @@
         it << Name ( "S1" ) << f1->evallog_mat ( Grid );
         it << Name ( "S2" ) << f2->evallog_mat ( Grid );
-        cout << ( ( enorm<ldmat>* ) ( MP->_Coms ( 0 ).get() ) )->_R().to_mat() << endl;
+//      cout << ( ( enorm<ldmat>* ) ( MP->_Coms ( 0 ).get() ) )->_R().to_mat() << endl;
 }

library/tests/stresssuite/mixtures_stress.cpp

@@ -15 +15 @@
         mat Var ( 2, 2 ), M ( 1, 2 );
         for ( i = 0; i < ep._rv().length() - 1; i++ ) { // -1 => last one is the weight
-                ( ( egiw* ) ep ( i ) )->mean_mat ( M, Var );
+                ( ( egiw* ) ep.factor ( i ) )->mean_mat ( M, Var );
                 cout << "Mean: " << M << endl << "Variance: " << endl << Var << endl;
         }
-        cout << "Weights: " << ep ( i )->mean() << endl;
+        cout << "Weights: " << ep.factor ( i )->mean() << endl;
 }
 

library/tests/testsuite/epdf_test.cpp

@@ -36 +36 @@
         g1->set_rv ( b );
         
-        Array<const epdf*> coms ( 2 );
-        coms ( 0 ) = g0.get();
-        coms ( 1 ) = g1.get();
+        Array<shared_ptr<epdf> > coms ( 2 );
+        coms ( 0 ) = g0;
+        coms ( 1 ) = g1;
         
         eprod p;
@@ -44 +44 @@
         // pointers in the array to outlast the eprod instance...
         p.set_parameters ( coms );
+        p.validate();
         
         CHECK_EQUAL ( vec ( "1 10 10" ), p.mean() );