Changeset 787

Timestamp:

01/14/10 21:37:57 (14 years ago)

Author:

smidl

Message:

arena experiment + numerical fixes

Files:

• applications/bdmtoolbox/tutorial/mpdm/dist_ctrl_acc.m (modified) (2 diffs)
• applications/bdmtoolbox/tutorial/mpdm/dist_ctrl_example.m (modified) (4 diffs)
• library/bdm/itpp_ext.cpp (modified) (1 diff)
• library/bdm/mpdm/arx_agent.h (modified) (3 diffs)
• library/bdm/stat/merger.cpp (modified) (2 diffs)

applications/bdmtoolbox/tutorial/mpdm/dist_ctrl_acc.m

@@ -50 +50 @@
 %P1.merger.dbg_file = 'mp.it';
 P1.merger.ncoms = 20;
+P1.merger.stop_niter= 5;
 P1.neighbours = {};%{'P2'};
 
@@ -88 +89 @@
     Y = [M.DS_y1 M.DS_y2 M.DS_y3];
     Yreq = ones(size(M.DS_y1))*[0 1 0];
-    loss_non_coop(i) = trace((Y-Yreq)'*0.01*(Y-Yreq)) + M.DS_u1'*C1.Qu*M.DS_u1 + M.DS_u2'*C1.Qu*M.DS_u2;
+    loss_coop(i) = trace((Y-Yreq)'*0.01*(Y-Yreq)) + M.DS_u1'*C1.Qu*M.DS_u1 + M.DS_u2'*C1.Qu*M.DS_u2;
     if loss_coop(i)>100
         %keyboard

applications/bdmtoolbox/tutorial/mpdm/dist_ctrl_example.m

@@ -31 +31 @@
 C1.Qy = 1*eye(1);
 C1.yreq = 1;
-C1.horizon = 5;
+C1.horizon = 1;
 
 C2=C1;
@@ -44 +44 @@
 %P1.merger.dbg_file = 'mp.it';
 P1.merger.ncoms = 20;
+P1.merger.stop_niter= 5;
 P1.neighbours = {};%{'P2'};
 
@@ -60 +61 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MONTE CARLO %%%%%%%%%%%%%%%%%%%
 
-Ntrials = 100;
+Ntrials = 3;
 loss_non_coop = zeros(1,Ntrials);
 for i=1:Ntrials
     M= arena(DS,{P1,P2},exper);
 
-    loss_non_coop(i) = M.DS_y'*C1.Qy*M.DS_y + M.DS_u1'*C1.Qu*M.DS_u1 + M.DS_u2'*C1.Qu*M.DS_u2;
+    loss_non_coop(i) = (M.DS_y-C1.yreq)'*C1.Qy*(M.DS_y-C1.yreq) + M.DS_u1'*C1.Qu*M.DS_u1 + M.DS_u2'*C1.Qu*M.DS_u2;
     if loss_non_coop(i)>100
         %keyboard
@@ -78 +79 @@
     M= arena(DS,{P1,P2},exper);
 
-    loss_coop(i) = M.DS_y'*C1.Qy*M.DS_y + M.DS_u1'*C1.Qu*M.DS_u1 + M.DS_u2'*C1.Qu*M.DS_u2;
+    loss_coop(i) = (M.DS_y-C1.yreq)'*C1.Qy*(M.DS_y-C1.yreq) + M.DS_u1'*C1.Qu*M.DS_u1 + M.DS_u2'*C1.Qu*M.DS_u2;
     if loss_coop(i)>100
         %keyboard

library/bdm/itpp_ext.cpp

@@ -390 +390 @@
 void triu ( mat &A ) {
         for ( int i = 1; i < A.rows(); i++ ) { // row cycle
-                for ( int j = 0; j < i; j++ ) {
+                for ( int j = 0; (j < i) && (j<A.cols()); j++ ) {
                         A ( i, j ) = 0;
                 }

library/bdm/mpdm/arx_agent.h

@@ -65 +65 @@
                 }
                 void adapt(const vec &glob_dt){                 
-                        dlDt.filldown(glob_dt, dt);
+                        dlDt.filldown(glob_dt, dt); 
+                        dlDt.store_data(glob_dt); // do not ever use dlDt again! 
+                        
                         lqg_arx->adapt(dt);
                         // remove old predictors
@@ -80 +82 @@
                                 merger->merge();
                                 enorm<chmat> joint_pred; 
-                                joint_pred.set_parameters(merger->mean(), merger->covariance());
+                                mat Cov=merger->covariance();
+                                if (sumsum(Cov)==0.0){
+                                        bdm_error("merging failed");
+                                }
+                                joint_pred.set_parameters(merger->mean(), Cov);
                                 joint_pred.set_rv(merger->_rv());
                                 
@@ -96 +102 @@
                         ut=lqg_arx->ctrlaction(dt);
                         dlU.filldown(ut,glob_ut);
-                        // save dat for next time
-                        dlDt.store_data(dt);
                 }
                 virtual void ds_register(const DS &ds){

library/bdm/stat/merger.cpp

@@ -163 +163 @@
         mat Tmp = zeros ( dim, dim );
         for ( int i = 0; i < Npoints; i++ ) {
-                Tmp += w ( i ) * outer_product ( S ( i ), S ( i ) );
-        }
-        return Tmp - outer_product ( mea, mea );
+                vec tmp=S ( i )-mea; //inefficient but numerically stable
+                Tmp += w ( i ) * outer_product (tmp , tmp );
+        }
+        return Tmp;
 }
 
@@ -174 +175 @@
 void merger_mix::merge ( ) {
         if(Npoints<1){
-                set_support(enorm<fsqmat>(zeros(dim), eye(dim)), 100);
+                set_support(enorm<fsqmat>(zeros(dim), eye(dim)), 1000);
         }