root/applications/robust/main.cpp @ 1401

/*!
\file
\brief Robust
\author Vasek Smidl 

 */

#include "estim/arx.h"
#include "robustlib.h"
#include <vector>
#include <iostream>
#include <fstream>
//#include <itpp/itsignal.h>
#include "windows.h"
#include "ddeml.h"
#include "stdio.h"

//#include "DDEClient.h"
//#include <conio.h>


using namespace itpp;
using namespace bdm;

//const int emlig_size = 2;
//const int utility_constant = 5;

const int max_model_order = 2;
const double apriorno     = 0.01;
const int max_window_size = 30;

class model
{
public:
        set<pair<int,int>> ar_components;

        // Best thing would be to inherit the two models from a single souce, this is planned, but now structurally
        // problematic.
        RARX* my_rarx; //vzmenovane parametre pre triedu model
        ARXwin* my_arx;

        bool has_constant;
        int  window_size;  //musi byt vacsia ako pocet krokov ak to nema ovplyvnit
        int  predicted_channel;
        mat* data_matrix;
        vec  predictions;
        char name[80];
        
        model(set<pair<int,int>> ar_components, //funkcie treidz model-konstruktor
                  bool robust, 
                  bool has_constant, 
                  int window_size, 
                  int predicted_channel,
                  mat* data_matrix)
        {
                this->ar_components.insert(ar_components.begin(),ar_components.end());

                strcpy(name,"M");

                for(set<pair<int,int>>::iterator ar_ref = ar_components.begin();ar_ref!=ar_components.end();ar_ref++)
                {
                        char buffer1[2];
                        char buffer2[2];
                        itoa((*ar_ref).first,buffer1,10);
                        itoa((*ar_ref).second,buffer2,10);

                        strcat(name,buffer1);
                        strcat(name,buffer2);
                        strcat(name,"_");
                }

                this->has_constant      = has_constant;

                if(has_constant)
                {
                        strcat(name,"C");
                }

                this->window_size       = window_size;
                this->predicted_channel = predicted_channel;
                this->data_matrix       = data_matrix;

                if(robust)
                {
                        strcat(name,"R");

                        if(has_constant)
                        {
                                my_rarx = new RARX(ar_components.size()+1,window_size,true,sqrt(2*apriorno),sqrt(2*apriorno),ar_components.size()+4);
                                my_arx  = NULL;
                        }
                else
                        {
                                my_rarx = new RARX(ar_components.size(),window_size,false,sqrt(2*apriorno),sqrt(2*apriorno),ar_components.size()+3);
                                my_arx  = NULL;
                        }
                }
                else
                {
                        my_rarx = NULL;
                        my_arx  = new ARXwin();
                        mat V0;

                        if(has_constant)
                        {                               
                                V0  = apriorno * eye(ar_components.size()+2); //aj tu konst
                                //V0(0,0) = 0;
                                my_arx->set_constant(true);                             
                        }
                        else
                        {
                                
                                V0  = apriorno * eye(ar_components.size()+1);//menit konstantu
                                //V0(0,1) = -0.01;
                                //V0(1,0) = -0.01;
                                my_arx->set_constant(false);                            
                                
                        }

                        my_arx->set_statistics(1, V0, V0.rows()+2);                     
                        my_arx->set_parameters(window_size);
                        my_arx->validate();

                        vec mean = my_arx->posterior().mean();
                        cout << mean << endl;
                }
        }

        void data_update(int time) //vlozime cas a ono vlozi do data_vector podmineky(conditions) a predikce, ktore pouzije do bayes
        {
                vec data_vector;
                for(set<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++)
                {  
                        data_vector.ins(data_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second));
                }
                if(my_rarx!=NULL)
                {       
                        data_vector.ins(0,(*data_matrix).get(predicted_channel,time));
                        my_rarx->bayes(data_vector);
                }
                else
                {
                        vec pred_vec;
                        pred_vec.ins(0,(*data_matrix).get(predicted_channel,time));
                        my_arx->bayes(pred_vec,data_vector);
                }
        }

        pair<vec,vec> predict(int sample_size, int time, itpp::Laplace_RNG* LapRNG)  //nerozumiem, ale vraj to netreba, nepouziva to
        {
                vec condition_vector;
                for(set<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++)
                {
                        condition_vector.ins(condition_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second+1));
                }

                if(my_rarx!=NULL)
                {
                        pair<vec,mat> imp_samples = my_rarx->posterior->sample(sample_size,false);

                        //cout << imp_samples.first << endl;                    
                        
                        vec sample_prediction;                  
                        for(int t = 0;t<sample_size;t++)
                        {
                                vec lap_sample = condition_vector;
                                
                                if(has_constant)
                                {
                                        lap_sample.ins(lap_sample.size(),1.0);
                                }
                                
                                lap_sample.ins(lap_sample.size(),(*LapRNG)());

                                sample_prediction.ins(0,lap_sample*imp_samples.second.get_col(t));                              
                        }

                        return pair<vec,vec>(imp_samples.first,sample_prediction);
                }       
                else
                {
                        mat samples = my_arx->posterior().sample_mat(sample_size);                      
                        
                        vec sample_prediction;
                        for(int t = 0;t<sample_size;t++)
                        {
                                vec gau_sample = condition_vector;
                                
                                if(has_constant)
                                {
                                        gau_sample.ins(gau_sample.size(),1.0);
                                }
                                
                                gau_sample.ins(gau_sample.size(),randn());

                                sample_prediction.ins(0,gau_sample*samples.get_col(t));                         
                        }

                        return pair<vec,vec>(ones(sample_prediction.size()),sample_prediction);
                }
        
        }


        static set<set<pair<int,int>>> possible_models_recurse(int max_order,int number_of_channels)
        {
                set<set<pair<int,int>>> created_model_types;            

                if(max_order == 1)
                {                       
                        for(int channel = 0;channel<number_of_channels;channel++)
                        {
                                set<pair<int,int>> returned_type;
                                returned_type.insert(pair<int,int>(channel,1));  
                                created_model_types.insert(returned_type);
                        }

                        return created_model_types;
                }
                else
                {
                        created_model_types = possible_models_recurse(max_order-1,number_of_channels);
                        set<set<pair<int,int>>> returned_types;
                        
                        for(set<set<pair<int,int>>>::iterator model_ref = created_model_types.begin();model_ref!=created_model_types.end();model_ref++)
                        {                               
                                
                                for(int order = 1; order<=max_order; order++)
                                {
                                        for(int channel = 0;channel<number_of_channels;channel++)
                                        {
                                                set<pair<int,int>> returned_type;
                                                pair<int,int> new_pair = pair<int,int>(channel,order);
                                                if(find((*model_ref).begin(),(*model_ref).end(),new_pair)==(*model_ref).end()) 
                                                {
                                                        returned_type.insert((*model_ref).begin(),(*model_ref).end()); 
                                                        returned_type.insert(new_pair);
                                                        

                                                        returned_types.insert(returned_type);                                                   
                                                }
                                        }
                                }
                        }

                        created_model_types.insert(returned_types.begin(),returned_types.end());

                        return created_model_types;
                }               
        }
};




int main ( int argc, char* argv[] ) 
{
        vector<vector<string>> strings;

        char* file_string =  "C:\\results\\normalM"; // "C:\\dataADClosePercDiff"; //  

        char dfstring[80];
        strcpy(dfstring,file_string);
        strcat(dfstring,".txt");
        
        
        mat data_matrix;
        ifstream myfile(dfstring);
        if (myfile.is_open())
        {               
                string line;
                while(getline(myfile,line))
                {                       
                        vec data_vector;
                        while(line.find(',') != string::npos) 
                        {                               
                                int loc2 = line.find('\n');
                                int loc  = line.find(',');
                                data_vector.ins(data_vector.size(),atof(line.substr(0,loc).c_str()));                           
                                line.erase(0,loc+1);                                    
                        }

                        data_matrix.ins_row(data_matrix.rows(),data_vector);
                }               

                myfile.close(); 
        }
        else
        {
                cout << "Can't open data file!" << endl;
        }

        set<pair<int,int>> model_type;
        model_type.insert(pair<int,int>(0,1));
        model_type.insert(pair<int,int>(0,2));

        vector<model*> models;

        ofstream myfilew;
                

        while(data_matrix.rows()!=0)
        {
                for(int i=0;i<models.size();i++)
                {
                        delete models[i];
                }
                
                models.clear();
                models.push_back(new model(model_type,true,false,max_window_size,0,&data_matrix));
                models.push_back(new model(model_type,false,false,max_window_size,0,&data_matrix));

                for(int time = max_model_order;time<max_window_size;time++) //time<data_matrix.cols() 
                {               
                        vec cur_res_lognc;
                
                        vector<string> nazvy;
                        for(vector<model*>::iterator model_ref = models.begin();model_ref!=models.end();model_ref++)
                        {
                                (*model_ref)->data_update(time);

                                cout << "Updated:" << time << endl;     

                                if(time == max_window_size-1)
                                {
                                        char fstring[80];                                       
                                        strcpy(fstring,file_string);
                                        strcat(fstring,"ml");                                   
                                        strcat(fstring,(*model_ref)->name);
                                        strcat(fstring,".txt");

                                        vec coords;
                                        if((*model_ref)->my_arx!=NULL)
                                        {
                                                coords = (*model_ref)->my_arx->posterior().est_theta();
                                        }
                                        else
                                        {
                                                coords = (*model_ref)->my_rarx->posterior->minimal_vertex->get_coordinates();                                           
                                        }

                                        myfilew.open(fstring,ios::app);

                                        for(int i=0;i<coords.size();i++)
                                        {
                                                myfilew << coords.get(i) << ",";
                                        }
                                        myfilew << endl;

                                        myfilew.close();                                
                                }
                        }
                }

                data_matrix.del_row(0);
        }
        


        

        return 0;
}