root/applications/bdmtoolbox/mex/estimator.cpp @ 700

Revision 700, 6.4 kB (checked in by smidl, 15 years ago)

Making tutorial/userguide example work again (changes of mpdf and bayes)

RevLine 
[412]1/*!
2\file
3\brief Application Estimator
[293]4
[412]5The general task of estimation is defined on the following scheme:
6\dot
7digraph estimation{
8        node [shape=box];
9        {rank="same"; "Data Source"; "Bayesian Model"}
10        "Data Source" -> "Bayesian Model" [label="data"];
[622]11        "Bayesian Model" -> "Result Logger" [label="estimated\n statistics"];
[412]12        "Data Source" -> "Result Logger" [label="Simulated\n data"];
13}
14\enddot
15
[622]16Here,
[412]17\li Data Source is an object (class DS) providing sequential data, \f$ [d_1, d_2, \ldots d_t] \f$.
18\li Bayesian Model is an object (class BM) performing Bayesian filtering,
19\li Result Logger is an object (class logger) dedicated to storing important data from the experiment.
20
21\section  cmd Command-line usage
22Execute command:
23\code
24$> estimator config_file.cfg
25\endcode
26
27Full description of the experiment is in the file config_file.cfg which is expected to have the following structure:
28\code
29system = {type = "DS_offspring", ...};      // definition of a data source
30estimator = {type = "BM_offspring", ...};   // definition of an estimator
31logger = {type = "logger_type",...};        // definition of a logger
32experiment = {ndat = 11000; };              // definition of number of data records
33\endcode
34
35The above description must be specialized to specific classes. See, \subpage arx_ui how to do it for estimation of an ARX model.
36
37\section ex Matlab usage
38Execute command:
39\code
40>> estimator('config_file.cfg');
41\endcode
42when using loggers storing results on hard drives, and
43\code
44>> Res=estimator('config_file.cfg');
45\endcode
46when using logger of the type \c "mex_logger". The results will be stored in structure \c M.
47
48 */
49
[700]50#include <estim/arx.h>
51#include <stat/emix.h>
52#include <base/datasources.h>
53#include <base/loggers.h>
54#include <estim/particles.h>
55#include <estim/kalman.h>
[293]56
[384]57//#include "mex_datasource.h"
[373]58
[293]59using namespace bdm;
60
[412]61#ifdef MEX
62#include <itpp/itmex.h>
[700]63#include <mex/mex_BM.h>
64#include <mex/mex_logger.h>
65#include <mex/mex_datasource.h>
[412]66
[293]67void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) {
68        // Check the number of inputs and output arguments
[622]69        if ( n_input<2 ) mexErrMsgTxt ( "Usage:\n"
70                                                "result=estimator(system, estimators, experiment, logger)\n"
71                                                "  system     = struct('class','datasource',...);  % Estimated system\n"
72                                                "  estimators = {struct('class','estimator',...),  % Estimators\n"
73                                                "                struct('class','estimator',...),...} \n"
74                                                "  === optional ==="
75                                                "  experiment = struct('ndat',10);                 % number of data in experiment, full length of finite datasources, 10 otherwise \n"
76                                                "  logger     = struct('class','mexlogger');       % How to store results, default=mexlog, i.e. matlab structure\n\n"
77                                                "see documentation of classes datasource, BM, and mexlogger and their offsprings in BDM." );
78
[490]79        RV::clear_all();
[622]80        //CONFIG
[412]81        UImxArray Cfg;
[622]82        try {
83                Cfg.addGroup ( input[0],"system" );
84                Cfg.addList ( input[1],"estimators" );
85                if ( n_input>2 ) {
86                        Cfg.addGroup ( input[2],"experiment" );
87                }
88                if ( n_input>3 ) {
89                        Cfg.addGroup ( input[3],"logger" );
90                }/*else{
[412]91                // define logger as mexlog
92                Setting &S=Cfg.getRoot();
93                S.add("logger",Setting::TypeGroup);
94                S["logger"].add("class",Setting::TypeString);
95                S["logger"]["class"]="mexlog";
96                S["logger"].add("maxlen",Setting::TypeInt);
97                int maxlen;
98                S["experiment"].lookupValue("ndat",maxlen);
99                S["logger"]["maxlen"]=maxlen;
[609]100        }*/
[622]101        } catch ( SettingException e ) {
102                it_error ( "error: "+string ( e.getPath() ) );
103        }
[412]104
105        //DBG
[622]106        Cfg.writeFile ( "estimator.cfg" );
[412]107
108#else
109int main ( int argc, char* argv[] ) {
110        const char *fname;
111        if ( argc>1 ) {
112                fname = argv[1];
[293]113        } else {
[412]114                fname="estimator.cfg";
[293]115        }
[412]116        UIFile Cfg ( fname );
117#endif
[622]118
119        shared_ptr<DS> Ds = UI::build<DS> ( Cfg, "system" );
120        Array<shared_ptr<BM> > Es;
121        UI::get ( Es,Cfg, "estimators" );
[700]122        int Ndat=10;
[622]123        if ( Cfg.exists ( "experiment" ) ) {
[700]124                if ( Cfg.getRoot().lookupValue ( "experiment.ndat",Ndat ) ) {
[622]125                        bdm_assert ( Ndat<=Ds->max_length(), "Data source has less data then required" );
[609]126                };
127        } else {
[626]128                if ( Ds->max_length() < std::numeric_limits< int >::max() ) {
[609]129                        Ndat=Ds->max_length();
[622]130                }
131                ;// else Ndat=10;
[609]132        }
[622]133        shared_ptr<logger> L = UI::build<logger> ( Cfg, "logger",UI::optional );
134        if ( !L ) {
135#ifdef MEX
[609]136                //mex logger has only from_setting constructor - we  have to fill it...
[622]137                L=new mexlog ( Ndat );
138#else
[611]139                L=new stdlog();
[622]140#endif
[609]141        }
[622]142
[676]143        Ds->log_register ( *L, "DS" );
[412]144        string Ename;
145        Setting &S=Cfg;
[622]146        for ( int i=0; i<Es.length(); i++ ) {
[685]147                if (!UI::get ( Ename, S["estimators"][i], "name" ,UI::optional)){
[622]148                        Ename="Est"+num2str ( i );
[412]149                }
[622]150
[676]151                Es ( i )->log_register ( *L,Ename ); // estimate
[412]152        }
[622]153        L->init();
[293]154
[593]155        vec dt=zeros ( Ds->_drv()._dsize() );   //data variable
[622]156        Array<datalink*> Dls ( Es.length() );
[685]157        Array<datalink*> Dlsc ( Es.length() );
[626]158        Array<datalink_buffered*> Dls_buf (0);
[622]159        for ( int i=0; i<Es.length(); i++ ) {
[685]160                //connect actual data
161                Dls ( i ) = new datalink ( Es ( i )->_yrv(),Ds->_drv() ); //datalink between a datasource and estimator
162                //connect data in condition
163                if (Es ( i )->_rvc().mint()<0){ 
[631]164                        //delayed values are required
165                       
166                        //create delayed dl
[626]167                        int ith_buf=Dls_buf.size();
[631]168                        Dls_buf.set_size( ith_buf + 1, true);
[626]169                        Dls_buf(ith_buf) = new datalink_buffered(); 
[631]170                        //add dl to list of buffered DS
[685]171                        Dlsc(i) = Dls_buf(ith_buf);
172                        Dlsc(i)->set_connection ( Es ( i )->_rvc(),Ds->_drv() ); //datalink between a datasource and estimator
[626]173                       
174                } else {
[685]175                        Dlsc ( i ) = new datalink ( Es ( i )->_rvc(),Ds->_drv() ); //datalink between a datasource and estimator
[626]176                }
[412]177        }
[622]178
[609]179        for ( int tK=0;tK<Ndat;tK++ ) {
[593]180                Ds->getdata ( dt );                                     // read data
[676]181                Ds->log_write ( );
[622]182
183                for ( int i=0; i<Es.length(); i++ ) {
[685]184                        Es ( i )->bayes ( Dls ( i )->pushdown ( dt ), Dlsc(i) ->pushdown(dt) );         // update estimates
[676]185                        Es ( i )->log_write ();
[412]186                }
[626]187               
[293]188                L->step();
[609]189                Ds->step();                                                     // simulator step
[626]190                //update buffered fdat links
191                for (int i=0; i<Dls_buf.length(); i++){
192                        Dls_buf(i)->step(dt);
193                }
194                       
[293]195        }
196
197        L->finalize();
198        // ------------------ End of routine -----------------------------
199
[412]200#ifdef MEX
[622]201        mexlog* mL=dynamic_cast<mexlog*> ( L.get() );
[412]202
[622]203        if ( mL ) { // user wants output!!
[293]204                if ( n_output<1 ) mexErrMsgTxt ( "Wrong number of output variables!" );
205                output[0] = mL->toCell();
206        }
[412]207#endif
[685]208        for (int i=0;i<Dls.length(); i++){delete Dls(i); delete Dlsc(i);}
[293]209}
Note: See TracBrowser for help on using the browser.