Changeset 11

Timestamp:

01/25/08 10:31:08 (17 years ago)

Author:

smidl

Message:

resampling in PF + initial port of functions

Files:

• Makefile (modified) (2 diffs)
• itpp_ext.cpp (modified) (1 diff)
• libFN.cpp (added)
• libFN.h (added)
• libPF.cpp (modified) (3 diffs)
• libPF.h (modified) (2 diffs)
• mixpp.kdevelop (modified) (3 diffs)
• mixpp.kdevelop.filelist (modified) (1 diff)
• testKF.cpp (modified) (1 diff)
• testPF.cpp (modified) (1 diff)

Makefile

@@ -1 @@
-LDFLAGS=-litpp
+LDFLAGS=-Ldebug/usr/lib/ -litpp 
 CPPFLAGS=-g
 
-all: testKF
+all: testPF
 
 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

itpp_ext.cpp

@@ -23 +23 @@
   }
 
+
 }

libPF.cpp

@@ -5 +5 @@
 using std::endl;
 
+PF::PF(vec w0){w=w0/sum(w0); n=w.length();}
+
 ivec PF::resample( RESAMPLING_METHOD method ) {
-        ivec ind( n );
+        ivec ind=zeros_i( n );
+        ivec N_babies = zeros_i( n );
+        vec cumDist = cumsum( w );
+        vec u( n );
+        int i,j,parent;
+        double u0;
+
+        switch ( method ) {
+        case MULTINOMIAL:
+                u( n - 1 ) = pow( URNG.sample(), 1.0 / n );
+                for ( i = n - 2;i >= 0;i-- ) {
+                        u( i ) = u( i + 1 ) * pow( URNG.sample(), 1.0 / ( i + 1 ) );
+                }
+                break;
+        case STRATIFIED:
+                for ( i = 0;i < n;i++ ) {
+                        u( i ) = ( i + URNG.sample() ) / n;
+                }
+                break;
+        case SYSTEMATIC:
+                u0 = URNG.sample();
+                for ( i = 0;i < n;i++ ) {
+                        u( i ) = ( i + u0 ) / n;
+                }
+                break;
+        default:
+                it_error( "PF::resample(): Unknown resampling method" );
+        }
+        
+//      // <<<<<<<<<<<
+//      using std::cout;
+//      cout << "resampling:" << endl;
+//      cout << w <<endl;
+//      cout << cumDist <<endl;
+//      cout << u <<endl;
+
+        // U is now full
+        j = 0;
+        for ( i = 0;i < n;i++ ) {
+                while ( u( i ) > cumDist( j ) ) j++;
+                N_babies( j )++;
+        }
+
+//      // <<<<<<<<<<<
+//      cout << N_babies <<endl;
+
+        // We have assigned new babies for each Particle
+        // Now, we fill the resulting index such that:
+        // * particles with at least one baby should not move *
+        // This assures that reassignment can be done inplace;
+
+        // find the first parent;
+        parent=0; while(N_babies(parent)==0) parent++;
+        // Build index
+        for ( i = 0;i < n;i++ ) {
+                if ( N_babies( i ) > 0 ){
+                        ind( i ) = i;
+                        N_babies( i ) --; //this index was now replicated;
+                        }
+                else {
+                        // test if the parent has been fully replicated 
+                        // if yes, find the next one
+                        while((N_babies(parent)==0) || (N_babies(parent)==1 && parent>i) ) parent++;
+                        // Replicate parent
+                        ind (i) = parent;
+                        N_babies( parent ) --; //this index was now replicated;
+                }
+/*                      // <<<<<<<<<<<
+                        cout <<ind << endl;*/
+        }
         return ind;
 }
 
-TrivialPF::TrivialPF(mpdf &par0, mpdf &obs0, mpdf &prop0, int n0){
+TrivialPF::TrivialPF( mpdf &par0, mpdf &obs0, mpdf &prop0, int n0 ) : PF(1/n*ones(n)) {
         is_proposal = true;
         prop = &prop0;
@@ -16 +87 @@
         obs = &obs0;
 }
-TrivialPF::TrivialPF(mpdf &par0, mpdf &obs0,  int n0){
+TrivialPF::TrivialPF( mpdf &par0, mpdf &obs0,  int n0 ) : PF(1/n*ones(n)) {
         is_proposal = false;
         par = &par0;
@@ -22 +93 @@
 }
 
-void TrivialPF::bayes( const vec &dt , bool evalll) {
+void TrivialPF::bayes( const vec &dt , bool evalll ) {
         int i;
         vec oldp;
-        double ll, gl, sum=0.0;
+        double ll, gl, sum = 0.0;
         Sort<double> S;
         ivec ind, iw;

libPF.h

@@ -20 +20 @@
 using namespace itpp;
 
-enum RESAMPLING_METHOD { MULTINOMIAL = 0, DETERMINISTIC = 1, RESIDUAL = 2, SYSTEMATIC = 3 };
+enum RESAMPLING_METHOD { MULTINOMIAL = 0, STRATIFIED = 1, SYSTEMATIC = 3 };
 
 /*!
@@ -31 +31 @@
         int n; //number of particles
         vec w; //particle weights
+        Uniform_RNG URNG; //used for resampling
         
 public:
-        //! Returns indexes of particles that should be resampled
+        //! Returns indexes of particles that should be resampled. The ordering MUST guarantee inplace replacement. (Important for MPF.)
         ivec resample(RESAMPLING_METHOD method = SYSTEMATIC);
-        //TODO get them on the web
+        PF (vec w);
+        //TODO remove or implement bayes()!
+        void bayes(const vec &dt, bool evell){};
 };
 

mixpp.kdevelop

@@ -18 +18 @@
     <run>
       <directoryradio>executable</directoryradio>
-      <mainprogram>./testKF</mainprogram>
+      <mainprogram>/home/smidl/work/mixpp/testPF</mainprogram>
       <programargs></programargs>
       <globaldebugarguments></globaldebugarguments>
@@ -200 +200 @@
       <orientation>Vertical</orientation>
     </splitheadersource>
+    <references/>
   </kdevcppsupport>
   <cppsupportpart>
@@ -215 +216 @@
   </kdevdocumentation>
   <kdevfileview>
-    <groups/>
+    <groups>
+      <hidenonprojectfiles>false</hidenonprojectfiles>
+      <hidenonlocation>false</hidenonlocation>
+    </groups>
+    <tree>
+      <showvcsfields>false</showvcsfields>
+      <hidenonprojectfiles>true</hidenonprojectfiles>
+      <hidepatterns>*.o,*.lo,CVS</hidepatterns>
+    </tree>
   </kdevfileview>
   <ctagspart>
-    <customArguments></customArguments>
+    <customArguments/>
     <customTagfilePath>/home/smidl/work/mixpp/tags</customTagfilePath>
     <activeTagsFiles/>

mixpp.kdevelop.filelist

@@ -1 +1 @@
 # KDevelop Custom Project File List
 Makefile
-arx.cpp
 doc
 doc/latex

testKF.cpp

@@ -17 +17 @@
         // input from Matlab
         it_file fin( "testKF.it" );
+
         mat Dt, Xt,Xt2;
         int Ndat;

testPF.cpp

@@ -19 +19 @@
         ldmat R("1");
         
+        
+        vec ptc = randn(10);
+        vec w = exp(-0.5*(pow(ptc,2))/0.2);
+        
+        cout << "p:" << ptc << endl;
+        cout << "w:" << w << endl;
+        
+        PF pf(w);
+        ivec ind = pf.resample();
+        
+        cout << ind << endl;
+        /*
         mlnorm<ldmat> obs(x,xm,A,R);
         mlnorm<ldmat> par(y,x,A,R);
         
         TrivialPF TPF(obs,par,10);
+        */
         //Exit program:
         return 0;