root/applications/dual/src/iterativemc.cpp @ 650

/*!
\file
\brief Test file for dual control of ARX models

Todo, more info...
 */

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

#include "arx1_ctrl.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<1 ) mexErrMsgTxt ( "Usage:\n"
             "  result=iterativemc(config_str)\n" 
                 "with config_str with fields .b, .y0, .ndat, controller"
                 );
        
        RV::clear_all();
        //zde se zpracuji vstupni data
        UImxArray Cfg(input[0]);
        //DBG
        Cfg.writeFile ( "iterativemc.cfg" );
        
#else
int main ( int argc, char* argv[] ) {
        const char *fname;
        if ( argc>1 ) {
                fname = argv[1];
        } else {
                fname="iterativemc.cfg";
        }
        UIFile Cfg ( fname );
#endif

        // parametry vyskytujici se v algoritmu
        double b;
        double sigma;
        double ytm;
        double yr;
        int Ndat;
        
                // cteni ze souboru
                UI::get(b,Cfg,"b",UI::compulsory);
                UI::get(sigma,Cfg,"sigma",UI::compulsory);
                if(!UI::get(ytm,Cfg,"yt0",UI::optional)){
                        ytm=0.0; // start at zero
                }
                if(!UI::get(yr,Cfg,"yr",UI::optional)){
                        yr=1.0; // requested value
                }
                if(!UI::get(Ndat,Cfg,"ndat",UI::optional)){
                        Ndat=30;
                }
                shared_ptr<Controller> C=UI::build<Controller>(Cfg,"controller",UI::compulsory);
                
                int seed;
                if(UI::get(seed,Cfg,"seed",UI::optional)){
                        RNG_reset(seed);
                } else {
                        RNG_randomize();
                }
        // zpusob ukladani vysledku, pokud se zada jiny, lze zadat jako polocku logger v konfiguraci
        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 stdlog();
                #endif
        }
        
        // create ARX system
        mlnorm<fsqmat> Sys;
        mat theta=ones(1,2);
        theta(1) = b;
        Sys.set_parameters(theta,vec_1(0.0),sigma*sigma*eye(1));
        
        // decide what will be stored
        RV y("y",1); // name random variable
        RV u("u",1); // name random variable
        int L_yt=L->add(y);
        int L_ut=L->add(u);
        C->log_add(*L);
        L->init();
        
        vec psi(2);       // regressor
        double ut;        // input
        vec yt;
        
        for ( int tK=0;tK<Ndat;tK++ ) {
                //yt has now meaning of yt-1
                //design controller
                C->adapt(vec_1(ytm));
                ut = C->ctrlaction(vec_1(ytm))(0);   // the best controller there is
                
                //prepare regressor
                psi(0) = ytm; // now it is t-1
                psi(1) = ut;
                // simulate system
                yt =    Sys.samplecond(psi); 
                ytm = yt(0);
                
                //save results
                L->logit(L_yt, yt);
                L->logit(L_ut, vec_1(ut));
                C->logit(*L);
                L->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
}