root/applications/robust/main.cpp @ 1375

/*!
\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 <time.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 int data_count      = 5;

static vec avg_vec;

HDDEDATA CALLBACK DdeCallback(
UINT uType,     // Transaction type.
UINT uFmt,      // Clipboard data format.
HCONV hconv,    // Handle to the conversation.
HSZ hsz1,       // Handle to a string.
HSZ hsz2,       // Handle to a string.
HDDEDATA hdata, // Handle to a global memory object.
DWORD dwData1,  // Transaction-specific data.
DWORD dwData2)  // Transaction-specific data.
{
        char  * pData;
        DWORD   data_len;
        
        switch (uType)
        {     
        case XTYP_ADVDATA: 

                if ( hdata == NULL )
                        cout << "No data in the message." << endl;
                else
                {
                        pData = (char*)DdeAccessData( hdata, &data_len );
                                                
                        avg_vec.ins(0,atof(pData));
                        
                        cout << pData << endl;

                        DdeUnaccessData( hdata );
                }
        }               
        return (HDDEDATA) NULL; 

}



void DDERequest(DWORD idInst, HCONV hConv, char* szItem)
{
        HSZ hszItem = DdeCreateStringHandle(idInst, szItem, 0);
        HDDEDATA hData = DdeClientTransaction(NULL,0,hConv,hszItem,CF_TEXT, 
                                                                        XTYP_ADVSTART,TIMEOUT_ASYNC , NULL); //TIMEOUT_ASYNC
        if (hData==NULL)
        {
                printf("Request failed: %s\n", szItem);
        }

        if (hData==0)
        {
                printf("Request failed: %s\n", szItem);
        }
}

DWORD WINAPI ThrdFunc( LPVOID n )
{    
    return 0;
}


class model
{
public:
        list<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;
        ARXwin* my_arx;

        bool has_constant;
        int  order;
        int  window_size;
        int  predicted_channel;
        int  number_of_updates;
        double prior_lognc;
        vec  lognc_history;
        mat* data_matrix;
        
        model(list<pair<int,int>> ar_components, 
                  bool robust, 
                  bool has_constant, 
                  int window_size, 
                  int predicted_channel,                  
                  mat* data_matrix)
        {
                this->ar_components.insert(this->ar_components.begin(),ar_components.begin(),ar_components.end());
                this->has_constant        = has_constant;
                this->window_size         = window_size;
                this->predicted_channel   = predicted_channel;          
                this->data_matrix         = data_matrix;
                number_of_updates = 0;

                order = 0;
                for(list<pair<int,int>>::iterator ar_ref = ar_components.begin();ar_ref!=ar_components.end();ar_ref++)
                {
                        if((*ar_ref).second > order)
                        {
                                order = (*ar_ref).second;
                        }
                }

                if(robust)
                {
                        if(has_constant)
                        {
                                my_rarx = new RARX(ar_components.size()+1,window_size,true);
                                my_arx  = NULL;
                        }
                        else
                        {
                                my_rarx = new RARX(ar_components.size(),window_size,false);
                                my_arx  = NULL;
                        }

                        prior_lognc = my_rarx->posterior->log_nc;
                }
                else
                {
                        my_rarx = NULL;
                        my_arx  = new ARXwin();
                        mat V0;

                        if(has_constant)
                        {                               
                                V0  = 0.001 * eye(ar_components.size()+2);
                                V0(0,0) = 1;
                                my_arx->set_constant(true);     
                                
                        }
                        else
                        {
                                
                                V0  = 0.001 * eye(ar_components.size()+1);
                                V0(0,0) = 1;
                                my_arx->set_constant(false);                            
                                
                        }

                        my_arx->set_statistics(1, V0); 
                        my_arx->set_parameters(window_size);
                        my_arx->validate();

                        prior_lognc = my_arx->posterior().lognc();
                }
        }

        void data_update(int time)
        {
                if(time > order)
                {
                        vec data_vector;
                        for(list<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));
                        }

                        // cout << "Update cond: " << data_vector << endl;

                        double cur_lognc;
                        if(my_rarx!=NULL)
                        {
                                data_vector.ins(0,(*data_matrix).get(predicted_channel,time));
                                my_rarx->bayes(data_vector);
                                cur_lognc = my_rarx->posterior->log_nc;
                        }
                        else
                        {
                                vec pred_vec;
                                pred_vec.ins(0,(*data_matrix).get(predicted_channel,time));
                                my_arx->bayes(pred_vec,data_vector);
                                cur_lognc = my_arx->posterior().lognc();
                        }

                        number_of_updates++;

                        if(number_of_updates>window_size)
                        {
                                lognc_history.ins(lognc_history.size(),cur_lognc-prior_lognc);
                        }
                }
                else
                {
                        lognc_history.ins(lognc_history.size(),0);
                }
        }

        pair<vec,vec> predict(int sample_size, int time, itpp::Laplace_RNG* LapRNG)
        {
                if(time > order)
                {
                        vec condition_vector;
                        for(list<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));
                        }

                        // cout << "Prediction cond: " << condition_vector << endl;

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

                                //cout << "Point estimate: " << (imp_samples.second*imp_samples.first)/(imp_samples.first*ones(imp_samples.first.size())) << 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(0,(*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);

                                //cout << "Point estimate: " << (samples*ones(samples.cols()))/samples.cols() << endl;

                                // cout << samples.get_col(1) << endl;
                                
                                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);
                        }
                }
                else
                {
                        
                        return pair<vec,vec>(zeros(1),ones(1));
                }
        
        }


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

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

                        return created_model_types;
                }
                else
                {
                        created_model_types = possible_models_recurse(max_order-1,number_of_channels);
                        list<list<pair<int,int>>> returned_types;
                        
                        for(list<list<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++)
                                        {
                                                list<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(returned_type.begin(),(*model_ref).begin(),(*model_ref).end());
                                                        returned_type.push_back(new_pair);
                                                        returned_types.push_back(returned_type);                                                        
                                                }
                                        }
                                }
                        }

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

                        return created_model_types;
                }               
        }
};







int main ( int argc, char* argv[] ) {   
        
        vec data_vec;
        mat data_matrix;
        DWORD Id;

        char  szApp[]   = "MT4";
        char  szTopic[] = "BID";        
        char  szItem[]  = "EURUSD";
        char* itRef     = &szItem[0];

        char* file_string = "c:\\rtdata";

        ofstream myfile;
        char fstring[80];                                       
        strcpy(fstring,file_string);
        strcat(fstring,itRef);
        strcat(fstring,"lognc.txt");

        list<list<pair<int,int>>> model_types = model::possible_models_recurse(max_model_order,1);
        
        int max_window = 50;
        int min_window = 1;

        list<model*> models;
        for(list<list<pair<int,int>>>::iterator model_type = model_types.begin();model_type!=model_types.end();model_type++)
        {
                for(int window_size = min_window;window_size < max_window;window_size++)
                {                       
                        //models.push_back(new model((*model_type),true,true,window_size,0,&data_matrix));
                        models.push_back(new model((*model_type),false,true,window_size,0,&data_matrix));                       
                        //models.push_back(new model((*model_type),true,false,window_size,0,&data_matrix));
                        //models.push_back(new model((*model_type),false,false,window_size,0,&data_matrix));                    
                }
        }
        
        mat result_lognc;
        mat result_preds;
        
        /*
        HANDLE hThrd = CreateThread( NULL, 0, (LPTHREAD_START_ROUTINE)ThrdFunc, (LPVOID)1, 0, &Id);
        if ( !hThrd )
    {
        cout<<"Error Creating Threads,,,,.exiting"<<endl;
        return -1;
    }
    Sleep ( 100 );
        */
        
                

        //DDE Initialization
        DWORD idInst = 0;
        UINT iReturn;
        iReturn = DdeInitialize(&idInst, (PFNCALLBACK)DdeCallback, 
                                                        APPCLASS_STANDARD | APPCMD_CLIENTONLY, 0 );
        if (iReturn!=DMLERR_NO_ERROR)
        {
                printf("DDE Initialization Failed: 0x%04x\n", iReturn);
                Sleep(1500);
                return 0;
        }

        //DDE Connect to Server using given AppName and topic.
        HSZ hszApp, hszTopic;
        HCONV hConv;
        hszApp = DdeCreateStringHandle(idInst, szApp, 0);
        hszTopic = DdeCreateStringHandle(idInst, szTopic, 0);
        hConv = DdeConnect(idInst, hszApp, hszTopic, NULL);
        //DdeFreeStringHandle(idInst, hszApp);
        //DdeFreeStringHandle(idInst, hszTopic);
        if (hConv == NULL)
        {
                printf("DDE Connection Failed.\n");
                Sleep(1500); DdeUninitialize(idInst);
                return 0;
        }

        //Execute commands/requests specific to the DDE Server.
        
        DDERequest(idInst, hConv, szItem);              
        
        itpp::Laplace_RNG LapRNG = Laplace_RNG();

        // time_t current_time = time(0);       

        while(1)
        {
                while(avg_vec.size() < data_count)
                {
                        MSG    msg;
                        BOOL   MsgReturn = GetMessage ( &msg , NULL , 0 , 0 );
                    
                        if(MsgReturn)
                        {
                                TranslateMessage(&msg);
                                DispatchMessage(&msg);                  
                        }
                        Sleep(500);
                }


                data_vec.ins(data_vec.size(),avg_vec*ones(avg_vec.size())/avg_vec.size());
                data_matrix = mat(data_vec).T();
                avg_vec.del(0,avg_vec.size()-1);                                
                                        
                vector<double> test_lognc;
                vec preds;

                for(list<model*>::iterator model_ref = models.begin();model_ref!=models.end();model_ref++)
                {
                        (*model_ref)->data_update(data_vec.size()-1);                   

                        if((*model_ref)->lognc_history.size()==0)
                        {
                                test_lognc.push_back(0);
                        }
                        else
                        {
                                test_lognc.push_back((ones((*model_ref)->lognc_history.size())*(*model_ref)->lognc_history)/(*model_ref)->lognc_history.size()*(max_window/(*model_ref)->window_size)); 
                        }

                        pair<vec,vec> predictions = (*model_ref)->predict(500,data_vec.size()-1,&LapRNG);

                        preds.ins(preds.size(),(predictions.first*predictions.second)/(predictions.first*ones(predictions.first.size())));
                        
                }

                //preds.ins(0,data_matrix.get(0,data_vec.size()));      
                
                result_preds.ins_col(result_preds.cols(),preds);                                                        

                myfile.open(fstring,ios::app);
                
                // myfile << my_rarx->posterior->minimal_vertex->get_coordinates()[0];
                
                for(int i = 0;i<test_lognc.size();i++)
                {
                        myfile << test_lognc[i] << ',';
                }

                myfile << endl;                         
                
                myfile.close();
                
        }
        
        //DDE Disconnect and Uninitialize.
        DdeDisconnect(hConv);
        DdeUninitialize(idInst);

        return 0;
}

                /*
                // EXPERIMENT: A moving window estimation and prediction of RARX is tested on data generated from 
                // y_t=0.95*y_(t-1)+0.05*y_(t-2)+0.2*e_t, where e_t is normally, student(4) and cauchy distributed. It
                // can be compared to the classical setup.
                
                vector<vector<string>> strings;

                char* file_string = "c:\\dataGCClosePercDiff"; 

                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;
                }
                */
                
                
                // cout << "Updated." << endl;
        
                /*
                // EXPERIMENT: 100 AR model generated time series of length of 30 from y_t=0.95*y_(t-1)+0.05*y_(t-2)+0.2*e_t, 
                // where e_t is normally, student(4) and cauchy distributed are tested using robust AR model, to obtain the 
                // variance of location parameter estimators and compare it to the classical setup.
                vector<vector<vector<string>>> string_lists;
                string_lists.push_back(vector<vector<string>>());
                string_lists.push_back(vector<vector<string>>());
                string_lists.push_back(vector<vector<string>>());

                char* file_strings[3] = {"c:\\ar_normal.txt", "c:\\ar_student.txt", "c:\\ar_cauchy.txt"};
                

                for(int i = 0;i<3;i++)
                {       
                        ifstream myfile(file_strings[i]);
                        if (myfile.is_open())
                        {
                                while ( myfile.good() )
                                {
                                        string line;
                                        getline(myfile,line);
                                        
                                        vector<string> parsed_line;
                                        while(line.find(',') != string::npos)
                                        {
                                                int loc = line.find(',');
                                                parsed_line.push_back(line.substr(0,loc));
                                                line.erase(0,loc+1);                                    
                                        }                               

                                        string_lists[i].push_back(parsed_line);
                                }
                                myfile.close();
                        }
                }

                for(int j = 0;j<string_lists.size();j++)
                { 
                        
                        for(int i = 0;i<string_lists[j].size()-1;i++)
                        {
                                vector<vec> conditions;
                                //emlig* emliga = new emlig(2);
                                RARX* my_rarx = new RARX(2,30);

                                for(int k = 1;k<string_lists[j][i].size();k++)
                                {
                                        vec condition;
                                        //condition.ins(0,1);                           
                                        condition.ins(0,string_lists[j][i][k]);                         
                                        conditions.push_back(condition);

                                        //cout << "orig:" << condition << endl;

                                        if(conditions.size()>1)
                                        {               
                                                conditions[k-2].ins(0,string_lists[j][i][k]);
                                                
                                        }

                                        if(conditions.size()>2)
                                        {
                                                conditions[k-3].ins(0,string_lists[j][i][k]);

                                                //cout << "modi:" << conditions[k-3] << endl;

                                                my_rarx->bayes(conditions[k-3]);

                                                
                                                //if(k>5)
                                                //{
                                                //      cout << "MaxLik coords:" << emliga->minimal_vertex->get_coordinates() << endl;
                                                //}
                                                
                                        }                               
                                        
                                }

                                //emliga->step_me(0);
                                /*
                                ofstream myfile;
                                myfile.open("c:\\robust_ar1.txt",ios::app);
                                myfile << my_rarx->minimal_vertex->get_coordinates()[0] << ";";
                                myfile.close();

                                myfile.open("c:\\robust_ar2.txt",ios::app);
                                myfile << emliga->minimal_vertex->get_coordinates()[1] << ";";
                                myfile.close();
                                

                                cout << "MaxLik coords:" << emliga->minimal_vertex->get_coordinates() << endl;
                                cout << "Step: " << i << endl;
                        }

                        cout << "One experiment finished." << endl;

                        ofstream myfile;
                        myfile.open("c:\\robust_ar1.txt",ios::app);
                        myfile << endl;
                        myfile.close();

                        myfile.open("c:\\robust_ar2.txt",ios::app);
                        myfile << endl;
                        myfile.close();
                }*/    

                /*
                vector<vec> conditions;
                //emlig* emliga = new emlig(2);
                RARX* my_rarx = new RARX(2,10,false);
                
                
                mat V0 = 0.0001 * eye ( 3 );
                ARX* my_arx = new ARX(0.85);
                my_arx->set_statistics ( 1, V0 ); //nu is default (set to have finite moments)
                my_arx->set_constant ( false );
                my_arx->validate();
                

                for(int k = 1;k<strings[j].size();k++)
                {
                        vec condition;
                        //condition.ins(0,1);                           
                        condition.ins(0,strings[j][k]);                         
                        conditions.push_back(condition);

                        //cout << "orig:" << condition << endl;

                        if(conditions.size()>1)
                        {               
                                conditions[k-2].ins(0,strings[j][k]);
                                        
                        }

                        if(conditions.size()>2)
                        {
                                conditions[k-3].ins(0,strings[j][k]);

                                // cout << "Condition:" << conditions[k-3] << endl;

                                my_rarx->bayes(conditions[k-3]);
                                //my_rarx->posterior->step_me(1);
                                
                                vec cond_vec;
                                cond_vec.ins(0,conditions[k-3][0]);
                                
                                my_arx->bayes(cond_vec,conditions[k-3].right(2));
                                        
                                /*
                                if(k>8)
                                {
                                        //my_rarx->posterior->step_me(0);

                                        //mat samples = my_rarx->posterior->sample_mat(10);

                                        pair<vec,mat> imp_samples = my_rarx->posterior->importance_sample(1000);

                                        //cout << imp_samples.first << endl;
                                        
                                        vec sample_prediction;
                                        vec averaged_params = zeros(imp_samples.second.rows());
                                        for(int t = 0;t<imp_samples.first.size();t++)
                                        {
                                                vec lap_sample = conditions[k-3].left(2);
                                                //lap_sample.ins(lap_sample.size(),1.0);
                                                
                                                lap_sample.ins(0,LapRNG());

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

                                                averaged_params += imp_samples.first[t]*imp_samples.second.get_col(t);
                                        }

                                        averaged_params = averaged_params*(1/(imp_samples.first*ones(imp_samples.first.size())));

                                        // cout << "Averaged estimated parameters: " << averaged_params << endl;
                                        
                                        vec sample_pow = sample_prediction;                                     
                                        
                                        // cout << sample_prediction << endl;
                                        vec poly_coefs;
                                        double prediction;
                                        bool stop_iteration = false;
                                        int en = 1;
                                        do
                                        {
                                                double poly_coef = imp_samples.first*sample_pow/(imp_samples.first*ones(imp_samples.first.size()));

                                                if(en==1)
                                                {
                                                        prediction = poly_coef;
                                                }

                                                poly_coef = poly_coef*en*fact(utility_constant-2+en)/fact(utility_constant-2);

                                                if(abs(poly_coef)>numeric_limits<double>::epsilon())
                                                {
                                                        sample_pow = elem_mult(sample_pow,sample_prediction);
                                                        poly_coefs.ins(0,pow(-1.0,en+1)*poly_coef);
                                                }
                                                else
                                                {
                                                        stop_iteration = true;
                                                }
                                                
                                                en++;

                                                if(en>20)
                                                {
                                                        stop_iteration = true;
                                                }
                                        }
                                        while(!stop_iteration);

                                        /*
                                        ofstream myfile_coef;                                           

                                        myfile_coef.open("c:\\coefs.txt",ios::app);
                                        
                                        for(int t = 0;t<poly_coefs.size();t++)
                                        {
                                                myfile_coef << poly_coefs[t] << ",";                                    
                                        }

                                        myfile_coef << endl;
                                        myfile_coef.close();
                                        */

                                        //cout << "Coefficients: " << poly_coefs << endl;
                                                                                
                                        /*
                                        vec bas_coef = vec("1.0 2.0 -8.0");
                                        cout << "Coefs: " << bas_coef << endl;
                                        cvec actions2 = roots(bas_coef);
                                        cout << "Roots: " << actions2 << endl;
                                        */
                                        
                                    /*

                                        cvec actions = roots(poly_coefs);
                                        

                                        bool is_max = false;
                                        for(int t = 0;t<actions.size();t++)
                                        {
                                                if(actions[t].imag() == 0)
                                                {
                                                        double second_derivative = 0;
                                                        for(int q = 1;q<poly_coefs.size();q++)
                                                        {
                                                                second_derivative+=poly_coefs[q]*pow(actions[t].real(),q-1)*q;
                                                        }

                                                        if(second_derivative<0)
                                                        {
                                                                cout << "Action:" << actions[t].real() << endl;

                                                                is_max = true;
                                                        }
                                                }
                                        }

                                        if(!is_max)
                                        {
                                                cout << "No maximum." << endl;
                                        }

                                        // cout << "MaxLik coords:" << my_rarx->posterior->minimal_vertex->get_coordinates() << endl;

                                        /*
                                        double prediction = 0;
                                        for(int s = 1;s<samples.rows();s++)
                                        {
                                                
                                                double avg_parameter = imp_samples.get_row(s)*ones(samples.cols())/samples.cols();

                                                prediction += avg_parameter*conditions[k-3][s-1];

                                                
                                                
                                                /*
                                                ofstream myfile;
                                                char fstring[80];
                                                strcpy(fstring,file_strings[j]);

                                                char es[5];
                                                strcat(fstring,itoa(s,es,10));

                                                strcat(fstring,"_res.txt");
                                                

                                                myfile.open(fstring,ios::app);
                                                
                                                //myfile << my_rarx->posterior->minimal_vertex->get_coordinates()[0];
                                                myfile << avg_parameter;
                                                
                                                if(k!=strings[j].size()-1)
                                                {
                                                        myfile << ",";
                                                }
                                                else
                                                {
                                                        myfile << endl;
                                                }
                                                myfile.close();
                                                */

                                        
                                        //}

                                        // cout << "Prediction: "<< prediction << endl;
                                        /*
                                        enorm<ldmat>* pred_mat = my_arx->epredictor(conditions[k-3].left(2));
                                        double prediction2 = pred_mat->mean()[0];
                                        */

                                        
                                        
                                        /*
                                        myfile.open(f2string,ios::app);
                                        myfile << prediction2;
                                        
                                        if(k!=strings[j].size()-1)
                                        {
                                                myfile << ",";
                                        }
                                        else
                                        {
                                                myfile << endl;
                                        }
                                        myfile.close();
                                        //*//*

                                }                                       
                        }       */
                        
                        //emliga->step_me(0);
                        /*
                        ofstream myfile;
                        myfile.open("c:\\robust_ar1.txt",ios::app);
                        myfile << my_rarx->minimal_vertex->get_coordinates()[0] << ";";
                        myfile.close();

                        myfile.open("c:\\robust_ar2.txt",ios::app);
                        myfile << emliga->minimal_vertex->get_coordinates()[1] << ";";
                        myfile.close();
                        

                        cout << "MaxLik coords:" << emliga->minimal_vertex->get_coordinates() << endl;
                        cout << "Step: " << i << endl;*/
                //}


        //}


        // EXPERIMENT: One step ahead price prediction. Comparison of classical and robust model using optimal trading
    //             with maximization of logarithm of one-step ahead wealth.

        
                
                /*
                cout << "One experiment finished." << endl;

                ofstream myfile;
                myfile.open("c:\\robust_ar1.txt",ios::app);
                myfile << endl;
                myfile.close();

                myfile.open("c:\\robust_ar2.txt",ios::app);
                myfile << endl;
                myfile.close();*/
        

        //emlig* emlig1 = new emlig(emlig_size);

        //emlig1->step_me(0);
        //emlig* emlig2 = new emlig(emlig_size);
        
        /*
        emlig1->set_correction_factors(4);

        for(int j = 0;j<emlig1->correction_factors.size();j++)
        {
                for(set<my_ivec>::iterator vec_ref = emlig1->correction_factors[j].begin();vec_ref!=emlig1->correction_factors[j].end();vec_ref++)
                {
                        cout << j << "    ";
                        
                        for(int i=0;i<(*vec_ref).size();i++)
                        {
                                cout << (*vec_ref)[i];
                        }

                        cout << endl;
                }
        }*/
        
        /*
    vec condition5 = "1.0 1.0 1.01";//"-0.3 1.7 1.5";

        emlig1->add_condition(condition5);
        //emlig1->step_me(0);


        vec condition1a = "-1.0 1.02 0.5";
        //vec condition1b = "1.0 1.0 1.01";
        emlig1->add_condition(condition1a);
        //emlig2->add_condition(condition1b);

        vec condition2a = "-0.3 1.7 1.5";
        //vec condition2b = "-1.0 1.0 1.0";
        emlig1->add_condition(condition2a);
        //emlig2->add_condition(condition2b);

        vec condition3a = "0.5 -1.01 1.0";
        //vec condition3b = "0.5 -1.01 1.0";

        emlig1->add_condition(condition3a);
        //emlig2->add_condition(condition3b);   

        vec condition4a = "-0.5 -1.0 1.0";
        //vec condition4b = "-0.5 -1.0 1.0";    

        emlig1->add_condition(condition4a);
        //cout << "************************************************" << endl;
        //emlig2->add_condition(condition4b);
        //cout << "************************************************" << endl;
        
        //cout << emlig1->minimal_vertex->get_coordinates();
        
        //emlig1->remove_condition(condition3a);
        //emlig1->step_me(0);
        //emlig1->remove_condition(condition2a);
        //emlig1->remove_condition(condition1a);
        //emlig1->remove_condition(condition5);
        

        //emlig1->step_me(0);
        //emlig2->step_me(0);
        

        // DA SE POUZIT PRO VYPIS DO SOUBORU
        // emlig1->step_me(0);

        //emlig1->remove_condition(condition1);
        
        

        
        
        /*
        for(int i = 0;i<100;i++)
        {
                cout << endl << "Step:" << i << endl;           

                double condition[emlig_size+1];         

                for(int k = 0;k<=emlig_size;k++)
                {
                        condition[k] = (rand()-RAND_MAX/2)/1000.0;              
                }
                        

                vec* condition_vec = new vec(condition,emlig_size+1);
                emlig1->add_condition(*condition_vec);

                /*
                for(polyhedron* toprow_ref = emlig1->statistic.rows[emlig_size]; toprow_ref != emlig1->statistic.end_poly; toprow_ref = toprow_ref->next_poly)
                {
                        cout << ((toprow*)toprow_ref)->probability << endl;
                }
                */
                /*
                cout << emlig1->statistic_rowsize(emlig_size) << endl << endl;
        
                /*
                if(i-emlig1->number_of_parameters >= 0)
                {
                        pause(30);
                }
                */

                // emlig1->step_me(i);
                
                /*
                vector<int> sizevector;
                for(int s = 0;s<=emlig1->number_of_parameters;s++)
                {
                        sizevector.push_back(emlig1->statistic_rowsize(s));
                }
                */
        //}
    


        
        /*
        emlig1->step_me(1);

        vec condition = "2.0 0.0 1.0";  

        emlig1->add_condition(condition);

        vector<int> sizevector;
        for(int s = 0;s<=emlig1->number_of_parameters;s++)
        {
                sizevector.push_back(emlig1->statistic_rowsize(s));
        }

        emlig1->step_me(2);

        condition = "2.0 1.0 0.0";

        emlig1->add_condition(condition);

        sizevector.clear();
        for(int s = 0;s<=emlig1->number_of_parameters;s++)
        {
                sizevector.push_back(emlig1->statistic_rowsize(s));
        }
        */