Changeset 12

Timestamp:

01/25/08 10:53:47 (17 years ago)

Author:

smidl

Message:

opravy v libDC (stale nefunkcni)

Files:

• Makefile (modified) (2 diffs)
• libBM.cpp (modified) (4 diffs)
• libBM.h (modified) (3 diffs)
• libDC.cpp (modified) (7 diffs)
• libDC.h (modified) (7 diffs)
• libEF.h (modified) (10 diffs)

Makefile

@@ -2 +2 @@
 CPPFLAGS=-g
 
-all: testPF
+all: testSmp
 
 test: test0
@@ -13 +13 @@
 testKF: testKF.o libKF.o libBM.o libDC.o itpp_ext.o
 testPF: testPF.o libPF.o libBM.o libDC.o itpp_ext.o
-
-itpp_ext.o: itpp_ext.cpp itpp_ext.h
+testSmp: testSmp.o libEF.o libBM.o libDC.o itpp_ext.o

libBM.cpp

@@ -9 +9 @@
 void RV::init( ivec in_ids, Array<std::string> in_names, ivec in_sizes, ivec in_times, ivec in_obs ) {
         //
+        int i;
         len = in_ids.length();
         //PRUDENT_MODE
@@ -25 +26 @@
         times = in_times;
         obs = in_obs;
+        size = 0;
+        for(i=0;i<len;i++){size+=sizes(i);}
 };
 
@@ -32 +35 @@
 
 RV::RV () {};
-
-inline int RV::length () {return len;};
 
 RV::RV ( ivec in_ids ) {
@@ -48 +49 @@
 }
 
-RV RV::rvsubselect(ivec ind){
+RV RV::subselect(ivec ind){
         return RV(ids(ind), names(to_Arr(ind)), sizes(ind), times(ind), obs(ind));
 }

libBM.h

@@ -25 +25 @@
 */
 class RV {
-        int len;
+        int size,len;
         ivec ids;
         ivec sizes;
@@ -46 +46 @@
 
         //! Return length (number of scalars) of the RV.
-        int length();
+        int count() {return size;}
+        //TODO why not inline and later??
+        
         //! Find indexes of another rv in self
-        ivec rvfind(RV rv2);
+        ivec find(RV rv2);
         //! Add (concat) another variable to the current one
-        RV rvadd(RV rv2);
+        RV add(RV rv2);
         //! Subtract  another variable from the current one
-        RV rvsubt(RV rv2);
+        RV subt(RV rv2);
         //! Select only variables at indeces ind
-        RV rvsubselect(ivec ind);
+        RV subselect(ivec ind);
         //! Select only variables at indeces ind
         RV operator()(ivec ind);
@@ -110 +112 @@
 };
 
+
+
 #endif // BM_H

libDC.cpp

@@ -18 +18 @@
         D = exD;
         L = exL;
+        dim = exD.length();
 }
 
@@ -23 +24 @@
         vec D ;
         mat L;
+        dim = 0;
 }
 
@@ -28 +30 @@
 //TODO check if correct!! Based on heuristic observation of lu()
 
-        int dim = V.cols();
+        dim = V.cols();
         it_assert_debug( dim == V.rows(),"ldmat::ldmat matrix V is not square!" );
 
-        mat U( dim,dim );
-
-        L = V; //Allocate space for L
-        ivec p = ivec( dim ); //not clear why?
-
-        lu( V,L,U,p );
-
-//Now, if V is symmetric, L is what we seek and D is on diagonal of U
-        D = diag( U );
-
-//check if V was symmetric
-//TODO How? norm of L-U'?
-//it_assert_debug();
+        
 }
 
@@ -62 +52 @@
 }
 
-std::ostream &operator<< ( std::ostream &os,  sqmat &sq ) {
-        os << sq.to_mat() << endl;
+std::ostream &operator<< ( std::ostream &os,  ldmat &ld ) {
+        os << "L:" << ld.L << endl;
+        os << "D:" << ld.D << endl;
 }
 
@@ -85 +76 @@
                 }
         }
-        return V;
+        mat V2 = L.transpose()*diag( D )*L;
+        return V2;
 }
 
@@ -177 +169 @@
 }
 
-ldmat& ldmat::operator *= (double x){
-int i;
-for(i=0;i<D.length();i++){D(i)*=x;};
-}
-
+ldmat& ldmat::operator *= ( double x ) {
+        int i;
+        for ( i=0;i<D.length();i++ ){D( i )*=x;};
+}
+
+vec ldmat::sqrt_mult( vec &x ) {
+        int i,j;
+        vec res( dim );
+        double sum;
+        for ( i=0;i<dim;i++ ) {//for each element of result
+                res( i ) = 0.0;
+                for ( j=i;j<dim;j++ ) {//sum D(j)*L(:,i).*x
+                        res( i ) += sqrt( D( j ) )*L( j,i )*x( j );
+                }
+        }
+        vec res2 = L.transpose()*diag( sqrt( D ) )*x;
+        return res2;
+}
+
+void ldmat::ldform( mat &A,vec &D0 ) {
+        int m = A.rows();
+        int n = A.cols();
+        int mn = (m<n) ? m :n ;
+
+        it_assert_debug( A.cols()==dim,"ldmat::ldform A is not compatible" );
+        it_assert_debug( D.length()==A.rows(),"ldmat::ldform Vector D must have the length as row count of A" );
+
+        L=concat_vertical( zeros( n,n ), diag( sqrt( D0 ) )*A );
+        D=zeros( n+m );
+        
+        //unnecessary big L and D will be made smaller at the end of file
+        
+        vec w=zeros( n );
+        vec v;
+        double sum, beta, pom;
+
+        int cc=0;
+        int i=n; // +1 in .m
+        int ii,j,jj;
+        while (( i> ( n-mn+1-cc ) )&&( i>1 ) ) {
+                i--;
+                sum = 0.0;
+                v.set_size( m+i-( n-cc ) ); //prepare v
+                for ( ii=n-cc;ii<m+i;i++ ) {
+                        sum+= L( ii,i )*L( ii,i );
+                        v( ii )=L( ii,i ); //assign v
+                }
+                if ( L( m+i,i )==0 ) {
+                        beta = sqrt( sum );
+                } else {
+                        beta = L( m+i,i )+sign( L( m+i,i ) )*sqrt( sum ) ;
+                }
+                if ( std::fabs( beta )<eps ) {
+                        cc++;
+                        L.set_row( n+1-cc, L.get_row( m+i ) );
+                        L.set_row( m+i,0 );
+                        D( m+i )=0; L( m+i,i )=1;
+                        L.set_submatrix( n+1-cc,m+i,i,i,0 );
+                        continue;
+                }
+
+                sum-=v( v.length()-1 )*v( v.length()-1 ); //
+                sum/=beta*beta;
+                sum++;
+
+                v/=beta;
+                v( v.length()-1 )=1;
+
+                pom=-2/sum;
+                for ( j=i;i>=0;i-- ) {
+                        w( j )=0.0;
+                        for ( ii=n-cc;ii<m+i;ii++ ) {
+                                w( j )+= v( ii )*L( ii,j );
+                        }
+                        w( j )*=pom;
+                }
+
+                for ( ii=n-cc;ii<m+i;ii++ ) {
+                        for ( jj=0;jj<i-1;jj++ ) {
+                                L( ii,jj )+= v( ii )*w( jj );
+                        }
+                }
+                for ( ii=n-cc;ii<m+i;ii++ ) {
+                        L( ii,i )= 0;
+                }
+                L( m+i,i )+=w( i );
+
+                D( m+i )=L( m+i,i )*L( m+i,i );
+                for ( ii=0;ii<i;ii++ ) {
+                        L( m+i,ii )/=L( m+i,i );
+                }
+        }
+        if ( i>0 ) {
+                for ( ii=0;ii<i-1;ii++ ) {
+                        jj = D.length()-1-n+ii;
+                        L.set_row(jj,0);
+                        L(jj,jj)=1;
+                }
+        }
+
+        //cut-out L and D;
+        L.del_rows(0,m-1);
+        D.del(0,m-1);
+}
 
 //////// Auxiliary Functions
@@ -255 +346 @@
                 *kr = 0.0;
                 if ( *Df < -threshold ) {
-                it_warning( "Problem in dydr: subraction of dyad results in negative definitness. Likely mistake in calling function." );}
+                        it_warning( "Problem in dydr: subraction of dyad results in negative definitness. Likely mistake in calling function." );
+                }
                 *Df = 0.0;
         }

libDC.h

@@ -52 +52 @@
 
         /*!
+        \brief Multiplies square root of $V$ by vector $x$.
+        
+        Used e.g. in generating normal samples.
+        */
+        virtual vec sqrt_mult(vec &v) =0;
+        
+        /*!
         \brief Evaluates quadratic form $x= v'*V*v$;
         
@@ -59 +66 @@
 //      //! easy version of the 
 //      sqmat inv();
-
-        friend std::ostream &operator<< ( std::ostream &os, sqmat &sq );
 
         //! Clearing matrix so that it corresponds to zeros.
@@ -71 +76 @@
         virtual int rows() =0;
 
+protected:
+        int dim;
 };
 
@@ -86 +93 @@
         void clear();
         
+        
         //! Constructor
         fsqmat(const mat &M);
@@ -95 +103 @@
         */
         virtual void inv(fsqmat* Inv);
+        
+
 };
 
@@ -118 +128 @@
         int cols();
         int rows();
+        vec sqrt_mult(vec &v);
 
         /*! \brief Matrix inversion preserving the chosen form.
@@ -133 +144 @@
    void mult_sym( const mat &C, ldmat &U, bool trans=false );
 
+        /*! \brief Transforms general $A'D0A$ into pure $L'DL$
+        
+        The new decomposition fullfills: $A'*diag(D)*A = self.L'*diag(self.D)*self.L$
+
+        @param A general matrix
+        @param D0 general vector
+        
+        */
+   void ldform( mat &A, vec &D0 );
+
         ldmat& operator += (const ldmat &ldA);
         ldmat& operator -= (const ldmat &ldA);
         ldmat& operator *= (double x);
         
+        friend std::ostream &operator<< ( std::ostream &os, ldmat &sq );
+
 protected:
         vec D;

libEF.h

@@ -15 +15 @@
 
 #include <itpp/itbase.h>
+#include "libDC.h"
+#include "libBM.h"
 //#include <std>
 
@@ -24 +26 @@
 * More?...
 */
-class eEF :public epdf {
+class eEF : public epdf {
 
 public:
-        virtual void tupdate( double phi, mat &vbar, double nubar ){};
-        virtual void dupdate( mat &v,double nu=1.0 ){};
+        virtual void tupdate( double phi, mat &vbar, double nubar ) {};
+        virtual void dupdate( mat &v,double nu=1.0 ) {};
 };
 
-class mEF :public mpdf {
+class mEF : public mpdf {
 
 public:
@@ -44 +46 @@
 template<class sq_T>
 class enorm : public eEF {
+        int dim;
         vec mu;
         sq_T R;
@@ -52 +55 @@
         void dupdate( mat &v,double nu=1.0 );
         vec sample();
+        mat sample(int N);
         double eval( const vec &val );
-
+        Normal_RNG RNG;
 };
 
@@ -76 +80 @@
 template<class sq_T>
 enorm<sq_T>::enorm( RV &rv, vec &mu0, sq_T &R0 ) {
+        dim = rv.count();
         mu = mu0;
         R = R0;
@@ -92 +97 @@
 template<class sq_T>
 vec enorm<sq_T>::sample() {
-        //
+        vec x( dim );
+        RNG.sample_vector( dim,x );
+        vec smp = R.sqrt_mult( x );
+
+        smp += mu;
+        return smp;
+};
+
+template<class sq_T>
+mat enorm<sq_T>::sample( int N ) {
+        mat X( dim,N );
+        vec x( dim );
+        vec pom;
+        int i;
+        for ( i=0;i<N;i++ ) {
+                RNG.sample_vector( dim,x );
+                pom = R.sqrt_mult( x );
+                pom +=mu;
+                X.set_col( i, pom);
+        }
+        return X;
 };
 
@@ -106 +131 @@
 template<class sq_T>
 mlnorm<sq_T>::mlnorm( RV &rv,RV &rvc, mat &A, sq_T &R ) {
-        int dim = rv.length();
+        int dim = rv.count();
         vec mu( dim );
 
@@ -116 +141 @@
         this->condition( cond );
         vec smp = epdf.sample();
-        lik = epdf.eval(smp);
+        lik = epdf.eval( smp );
         return smp;
 }
@@ -123 +148 @@
 mat mlnorm<sq_T>::samplecond( vec &cond, vec &lik, int n ) {
         int i;
-        int dim = rv.length();
+        int dim = rv.count();
         mat Smp( dim,n );
-        vec smp( dim);
+        vec smp( dim );
         this->condition( cond );
         for ( i=0; i<dim; i++ ) {
                 smp = epdf.sample();
-                lik(i) = epdf.eval(smp);
-                Smp.set_col( i ,smp);
+                lik( i ) = epdf.eval( smp );
+                Smp.set_col( i ,smp );
         }
         return Smp;
@@ -136 +161 @@
 
 template<class sq_T>
-void mlnorm<sq_T>::condition( vec &cond) {
+void mlnorm<sq_T>::condition( vec &cond ) {
+        epdf.mu = A*cond;
+//R is already assigned;
 }