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

/*!
\file
\brief Application Estimator

The general task of estimation is defined on the following scheme:
\dot
digraph estimation{
        node [shape=box];
        {rank="same"; "Data Source"; "Bayesian Model"}
        "Data Source" -> "Bayesian Model" [label="data"];
        "Bayesian Model" -> "Result Logger" [label="estimated\n statistics"]; 
        "Data Source" -> "Result Logger" [label="Simulated\n data"];
}
\enddot

Here, 
\li Data Source is an object (class DS) providing sequential data, \f$ [d_1, d_2, \ldots d_t] \f$.
\li Bayesian Model is an object (class BM) performing Bayesian filtering,
\li Result Logger is an object (class logger) dedicated to storing important data from the experiment.

\section  cmd Command-line usage
Execute command:
\code
$> estimator config_file.cfg
\endcode

Full description of the experiment is in the file config_file.cfg which is expected to have the following structure:
\code
system = {type = "DS_offspring", ...};      // definition of a data source
estimator = {type = "BM_offspring", ...};   // definition of an estimator
logger = {type = "logger_type",...};        // definition of a logger
experiment = {ndat = 11000; };              // definition of number of data records
\endcode

The above description must be specialized to specific classes. See, \subpage arx_ui how to do it for estimation of an ARX model.

\section ex Matlab usage
Execute command:
\code
>> estimator('config_file.cfg');
\endcode
when using loggers storing results on hard drives, and
\code
>> Res=estimator('config_file.cfg');
\endcode
when using logger of the type \c "mex_logger". The results will be stored in structure \c M.

 */

#include "estim/arx.h"
#include "stat/emix.h"
#include "base/datasources.h"
#include "base/loggers.h"

//#include "mex_datasource.h"

using namespace bdm;

#ifdef MEX
#include <itpp/itmex.h>
#include "mex/mex_BM.h"
#include "mex/mex_logger.h"
#include "mex/mex_datasource.h"

void mexFunction ( int n_output, mxArray *output[], int n_input, const mxArray *input[] ) {
        // Check the number of inputs and output arguments
        if(n_input<2) mexErrMsgTxt("Usage:\n"  
                "result=estimator(system, estimators, experiment, logger)\n"
                "  system     = struct('class','datasource',...);  % Estimated system\n"
                "  estimators = {struct('class','estimator',...),  % Estimators\n"
                "                struct('class','estimator',...),...} \n"
                "  === optional ==="
                "  experiment = struct('ndat',10);                 % number of data in experiment, full length of finite datasources, 10 otherwise \n"
                "  logger     = struct('class','mexlogger');       % How to store results, default=mexlog, i.e. matlab structure\n\n"
                "see documentation of classes datasource, BM, and mexlogger and their offsprings in BDM.");
        
        RV::clear_all();
        //CONFIG 
        UImxArray Cfg;
        try{
        Cfg.addGroup(input[0],"system");
        Cfg.addList(input[1],"estimators");
        if (n_input>2){
                Cfg.addGroup(input[2],"experiment");
        }
        if (n_input>3){
                Cfg.addGroup(input[3],"logger");
        }/*else{
                // define logger as mexlog
                Setting &S=Cfg.getRoot();
                S.add("logger",Setting::TypeGroup);
                S["logger"].add("class",Setting::TypeString);
                S["logger"]["class"]="mexlog";
                S["logger"].add("maxlen",Setting::TypeInt);
                int maxlen;
                S["experiment"].lookupValue("ndat",maxlen);
                S["logger"]["maxlen"]=maxlen;
        }*/
        } catch(SettingException e){it_error("error: "+string(e.getPath()));}

        //DBG
        Cfg.writeFile("estimator.cfg");

#else
int main ( int argc, char* argv[] ) {
        const char *fname;
        if ( argc>1 ) {
                fname = argv[1];
        } else {
                fname="estimator.cfg";
        }
        UIFile Cfg ( fname );
#endif
        
        shared_ptr<DS> Ds = UI::build<DS>( Cfg, "system" );
        Array<shared_ptr<BM> > Es;      UI::get(Es,Cfg, "estimators" );
        long Ndat=10;
        if (Cfg.exists("experiment")) {
                if (Cfg.lookupValue ( "experiment.ndat",Ndat )) {
                        bdm_assert(Ndat<=Ds->max_length(), "Data source has less data then required");
                };
        } else {
                if (Ds->max_length() < inf) {
                        Ndat=Ds->max_length();
                };// else Ndat=10;
        }
        shared_ptr<logger> L = UI::build<logger>( Cfg, "logger",UI::optional);
        if (!L) {
                #ifdef MEX
                //mex logger has only from_setting constructor - we  have to fill it...
                L=new mexlog(Ndat);
                #else
                L=new memlog(Ndat, "estimator_logger.it");
                #endif
        }
        
        Ds->log_add ( *L );
        string Ename;
        Setting &S=Cfg;
        for (int i=0; i<Es.length(); i++){
                try{
                        UI::get(Ename, S["estimators"][i], "name");
                } catch (UIException e){
                        Ename="Est"+num2str(i);
                }
                
                Es(i)->log_add(*L,Ename); // estimate
        }
                L->init();

        vec dt=zeros ( Ds->_drv()._dsize() );   //data variable
        Array<datalink*> Dls(Es.length()); 
        for (int i=0; i<Es.length(); i++){
                Dls(i)=new datalink( Es(i)->_drv(),Ds->_drv() ); //datalink between a datasource and estimator
        }
        
        for ( int tK=0;tK<Ndat;tK++ ) {
                Ds->getdata ( dt );                                     // read data
                Ds->logit ( *L );
                
                for (int i=0; i<Es.length(); i++){
                        Es(i)->bayes ( Dls(i)->pushdown ( dt ) );               // update estimates
                        Es(i)->logit (*L);
                }
                L->step();
                Ds->step();                                                     // simulator step
        }

        L->finalize();
        // ------------------ End of routine -----------------------------

#ifdef MEX
        mexlog* mL=dynamic_cast<mexlog*>(L.get());

        if (mL) { // user wants output!!
                if ( n_output<1 ) mexErrMsgTxt ( "Wrong number of output variables!" );
                output[0] = mL->toCell();
        }
#endif
}