Changeset 1376 for applications

Show
Ignore:
Timestamp:
06/21/11 19:09:56 (14 years ago)
Author:
sindj
Message:

Uprava integrace v souladu s clankem. JS

Location:
applications/robust
Files:
3 modified

Legend:

Unmodified
Added
Removed

  • applications/robust/main.cpp

    r1370 r1376  
    1212#include <iostream> 
    1313#include <fstream> 
    14 #include <itpp/itsignal.h> 
     14//#include <itpp/itsignal.h> 
    1515#include "windows.h" 
    1616#include "ddeml.h" 
    1717#include "stdio.h" 
    18 #include <time.h> 
    1918 
    2019//#include "DDEClient.h" 
    … …  
    2928 
    3029const int max_model_order = 2; 
    31 const int data_count      = 5; 
    32  
    33 static vec avg_vec; 
    34  
     30const double apriorno=0.02; 
     31 
     32/* 
    3533HDDEDATA CALLBACK DdeCallback( 
    36 UINT uType,     // Transaction type. 
    37 UINT uFmt,      // Clipboard data format. 
    38 HCONV hconv,    // Handle to the conversation. 
    39 HSZ hsz1,       // Handle to a string. 
    40 HSZ hsz2,       // Handle to a string. 
    41 HDDEDATA hdata, // Handle to a global memory object. 
    42 DWORD dwData1,  // Transaction-specific data. 
    43 DWORD dwData2)  // Transaction-specific data. 
     34        UINT uType,     // Transaction type. 
     35        UINT uFmt,      // Clipboard data format. 
     36        HCONV hconv,    // Handle to the conversation. 
     37        HSZ hsz1,       // Handle to a string. 
     38        HSZ hsz2,       // Handle to a string. 
     39        HDDEDATA hdata, // Handle to a global memory object. 
     40        DWORD dwData1,  // Transaction-specific data. 
     41        DWORD dwData2)  // Transaction-specific data. 
    4442{ 
    45         char  * pData; 
    46         DWORD   data_len; 
    47          
    48         switch (uType) 
    49         {      
    50         case XTYP_ADVDATA:  
    51  
    52                 if ( hdata == NULL ) 
    53                         cout << "No data in the message." << endl; 
    54                 else 
    55                 { 
    56                         pData = (char*)DdeAccessData( hdata, &data_len ); 
    57                                                  
    58                         avg_vec.ins(0,atof(pData)); 
    59                          
    60                         cout << pData << endl; 
    61  
    62                         DdeUnaccessData( hdata ); 
    63                 } 
    64         }                
    65         return (HDDEDATA) NULL;  
    66  
     43        return 0; 
    6744} 
    68  
    69  
    7045 
    7146void DDERequest(DWORD idInst, HCONV hConv, char* szItem) 
    … …  
    8964    return 0; 
    9065} 
    91  
     66*/ 
    9267 
    9368class model 
    9469{ 
    9570public: 
    96         list<pair<int,int>> ar_components; 
     71        set<pair<int,int>> ar_components; 
    9772 
    9873        // Best thing would be to inherit the two models from a single souce, this is planned, but now structurally 
    9974        // problematic. 
    100         RARX* my_rarx; 
    101         ARXwin* my_arx; 
     75        RARX* my_rarx; //vzmenovane parametre pre triedu model 
     76        ARX* my_arx; 
    10277 
    10378        bool has_constant; 
    104         int  order; 
    105         int  window_size; 
     79        int  window_size;  //musi byt vacsia ako pocet krokov ak to nema ovplyvnit 
    10680        int  predicted_channel; 
    107         int  number_of_updates; 
    108         double prior_lognc; 
    109         vec  lognc_history; 
    11081        mat* data_matrix; 
    11182         
    112         model(list<pair<int,int>> ar_components,  
     83        model(set<pair<int,int>> ar_components, //funkcie treidz model-konstruktor 
    11384                  bool robust,  
    11485                  bool has_constant,  
    11586                  int window_size,  
    116                   int predicted_channel,                   
     87                  int predicted_channel, 
    11788                  mat* data_matrix) 
    11889        { 
    119                 this->ar_components.insert(this->ar_components.begin(),ar_components.begin(),ar_components.end()); 
    120                 this->has_constant        = has_constant; 
    121                 this->window_size         = window_size; 
    122                 this->predicted_channel   = predicted_channel;           
    123                 this->data_matrix         = data_matrix; 
    124                 number_of_updates = 0; 
    125  
    126                 order = 0; 
    127                 for(list<pair<int,int>>::iterator ar_ref = ar_components.begin();ar_ref!=ar_components.end();ar_ref++) 
    128                 { 
    129                         if((*ar_ref).second > order) 
    130                         { 
    131                                 order = (*ar_ref).second; 
    132                         } 
    133                 } 
     90                this->ar_components.insert(ar_components.begin(),ar_components.end()); 
     91                this->has_constant      = has_constant; 
     92                this->window_size       = window_size; 
     93                this->predicted_channel = predicted_channel; 
     94                this->data_matrix       = data_matrix; 
    13495 
    13596                if(robust) 
    … …  
    140101                                my_arx  = NULL; 
    141102                        } 
    142                         else 
     103                else 
    143104                        { 
    144105                                my_rarx = new RARX(ar_components.size(),window_size,false); 
    145106                                my_arx  = NULL; 
    146107                        } 
    147  
    148                         prior_lognc = my_rarx->posterior->log_nc; 
    149108                } 
    150109                else 
    151110                { 
    152111                        my_rarx = NULL; 
    153                         my_arx  = new ARXwin(); 
     112                        my_arx  = new ARX(); 
    154113                        mat V0; 
    155114 
    156115                        if(has_constant) 
    157116                        {                                
    158                                 V0  = 0.001 * eye(ar_components.size()+2); 
     117                                V0  = apriorno * eye(ar_components.size()+2); //aj tu konst 
    159118                                V0(0,0) = 1; 
    160119                                my_arx->set_constant(true);      
    … …  
    164123                        { 
    165124                                 
    166                                 V0  = 0.001 * eye(ar_components.size()+1); 
     125                                V0  = apriorno * eye(ar_components.size()+1);//menit konstantu 
    167126                                V0(0,0) = 1; 
    168127                                my_arx->set_constant(false);                             
    … …  
    171130 
    172131                        my_arx->set_statistics(1, V0);  
    173                         my_arx->set_parameters(window_size); 
     132                        //my_arx->set_parameters(window_size); 
    174133                        my_arx->validate(); 
    175  
    176                         prior_lognc = my_arx->posterior().lognc(); 
    177                 } 
    178         } 
    179  
    180         void data_update(int time) 
    181         { 
    182                 if(time > order) 
    183                 { 
    184                         vec data_vector; 
    185                         for(list<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++) 
    186                         { 
    187                                 data_vector.ins(data_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second)); 
    188                         } 
    189  
    190                         // cout << "Update cond: " << data_vector << endl; 
    191  
    192                         double cur_lognc; 
    193                         if(my_rarx!=NULL) 
    194                         { 
    195                                 data_vector.ins(0,(*data_matrix).get(predicted_channel,time)); 
    196                                 my_rarx->bayes(data_vector); 
    197                                 cur_lognc = my_rarx->posterior->log_nc; 
    198                         } 
    199                         else 
    200                         { 
    201                                 vec pred_vec; 
    202                                 pred_vec.ins(0,(*data_matrix).get(predicted_channel,time)); 
    203                                 my_arx->bayes(pred_vec,data_vector); 
    204                                 cur_lognc = my_arx->posterior().lognc(); 
    205                         } 
    206  
    207                         number_of_updates++; 
    208  
    209                         if(number_of_updates>window_size) 
    210                         { 
    211                                 lognc_history.ins(lognc_history.size(),cur_lognc-prior_lognc); 
    212                         } 
     134                } 
     135        } 
     136 
     137        void data_update(int time) //vlozime cas a ono vlozi do data_vector podmineky(conditions) a predikce, ktore pouzije do bayes 
     138        { 
     139                vec data_vector; 
     140                for(set<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++) 
     141                {  //ar?iterator ide len od 1 pod 2, alebo niekedy len 1 
     142                        data_vector.ins(data_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second)); 
     143// do data vector vlozi pre dany typ regresoru prislusne cisla z data_matrix. Ale ako? preco time-ar_iterator->second? 
     144                } 
     145                if(my_rarx!=NULL) 
     146                {       //pre robusr priradi az tu do data_vector aj rpedikciu 
     147                        data_vector.ins(0,(*data_matrix).get(predicted_channel,time)); 
     148                        my_rarx->bayes(data_vector); 
    213149                } 
    214150                else 
    215151                { 
    216                         lognc_history.ins(lognc_history.size(),0); 
    217                 } 
    218         } 
    219  
    220         pair<vec,vec> predict(int sample_size, int time, itpp::Laplace_RNG* LapRNG) 
    221         { 
    222                 if(time > order) 
    223                 { 
    224                         vec condition_vector; 
    225                         for(list<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++) 
    226                         { 
    227                                 condition_vector.ins(condition_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second+1)); 
    228                         } 
    229  
    230                         // cout << "Prediction cond: " << condition_vector << endl; 
    231  
    232                         if(my_rarx!=NULL) 
    233                         { 
    234                                 pair<vec,mat> imp_samples = my_rarx->posterior->importance_sample(sample_size); 
    235  
    236                                 //cout << "Point estimate: " << (imp_samples.second*imp_samples.first)/(imp_samples.first*ones(imp_samples.first.size())) << endl; 
    237                                  
    238                                 vec sample_prediction;                   
    239                                 for(int t = 0;t<sample_size;t++) 
     152                        vec pred_vec;//tu sa predikcia zadava zvlast 
     153                        pred_vec.ins(0,(*data_matrix).get(predicted_channel,time)); 
     154                        my_arx->bayes(pred_vec,data_vector); 
     155                } 
     156        } 
     157 
     158        pair<vec,vec> predict(int sample_size, int time, itpp::Laplace_RNG* LapRNG)  //nerozumiem, ale vraj to netreba, nepouziva to 
     159        { 
     160                vec condition_vector; 
     161                for(set<pair<int,int>>::iterator ar_iterator = ar_components.begin();ar_iterator!=ar_components.end();ar_iterator++) 
     162                { 
     163                        condition_vector.ins(condition_vector.size(),(*data_matrix).get(ar_iterator->first,time-ar_iterator->second+1)); 
     164                } 
     165 
     166                if(my_rarx!=NULL) 
     167                { 
     168                        pair<vec,mat> imp_samples = my_rarx->posterior->importance_sample(sample_size); 
     169 
     170                        //cout << imp_samples.first << endl; 
     171                         
     172                        vec sample_prediction;                   
     173                        for(int t = 0;t<sample_size;t++) 
     174                        { 
     175                                vec lap_sample = condition_vector; 
     176                                 
     177                                if(has_constant) 
    240178                                { 
    241                                         vec lap_sample = condition_vector; 
    242                                          
    243                                         if(has_constant) 
    244                                         { 
    245                                                 lap_sample.ins(lap_sample.size(),1.0); 
    246                                         } 
    247                                          
    248                                         lap_sample.ins(0,(*LapRNG)()); 
    249  
    250                                         sample_prediction.ins(0,lap_sample*imp_samples.second.get_col(t));                               
     179                                        lap_sample.ins(lap_sample.size(),1.0); 
    251180                                } 
    252  
    253                                 return pair<vec,vec>(imp_samples.first,sample_prediction); 
    254                         } 
    255                         else 
    256                         { 
    257                                 mat samples = my_arx->posterior().sample_mat(sample_size); 
    258  
    259                                 //cout << "Point estimate: " << (samples*ones(samples.cols()))/samples.cols() << endl; 
    260  
    261                                 // cout << samples.get_col(1) << endl; 
    262                                  
    263                                 vec sample_prediction; 
    264                                 for(int t = 0;t<sample_size;t++) 
     181                                 
     182                                lap_sample.ins(0,(*LapRNG)()); 
     183 
     184                                sample_prediction.ins(0,lap_sample*imp_samples.second.get_col(t));                               
     185                        } 
     186 
     187                        return pair<vec,vec>(imp_samples.first,sample_prediction); 
     188                } 
     189                else 
     190                { 
     191                        mat samples = my_arx->posterior().sample_mat(sample_size); 
     192                         
     193                        vec sample_prediction; 
     194                        for(int t = 0;t<sample_size;t++) 
     195                        { 
     196                                vec gau_sample = condition_vector; 
     197                                 
     198                                if(has_constant) 
    265199                                { 
    266                                         vec gau_sample = condition_vector; 
    267                                          
    268                                         if(has_constant) 
    269                                         { 
    270                                                 gau_sample.ins(gau_sample.size(),1.0); 
    271                                         } 
    272                                          
    273                                         gau_sample.ins(gau_sample.size(),randn()); 
    274  
    275                                         sample_prediction.ins(0,gau_sample*samples.get_col(t));                          
     200                                        gau_sample.ins(gau_sample.size(),1.0); 
    276201                                } 
    277  
    278                                 return pair<vec,vec>(ones(sample_prediction.size()),sample_prediction); 
    279                         } 
     202                                 
     203                                gau_sample.ins(0,randn()); 
     204 
     205                                sample_prediction.ins(0,gau_sample*samples.get_col(t));                          
     206                        } 
     207 
     208                        return pair<vec,vec>(ones(sample_prediction.size()),sample_prediction); 
     209                } 
     210         
     211        } 
     212 
     213 
     214        static set<set<pair<int,int>>> possible_models_recurse(int max_order,int number_of_channels) 
     215        { 
     216                set<set<pair<int,int>>> created_model_types;             
     217 
     218                if(max_order == 1)//ukoncovacia vetva 
     219                {                        
     220                        for(int channel = 0;channel<number_of_channels;channel++)//pre AR 1 model vytvori kombinace kanalov v prvom kroku poyadu 
     221                        { 
     222                                set<pair<int,int>> returned_type; 
     223                                returned_type.insert(pair<int,int>(channel,1));  //?? 
     224                                created_model_types.insert(returned_type); 
     225                        } 
     226 
     227                        return created_model_types; 
    280228                } 
    281229                else 
    282230                { 
     231                        created_model_types = possible_models_recurse(max_order-1,number_of_channels);//tu uz mame ulozene kombinace o jeden krok dozadu  //rekuryivne volanie 
     232                        set<set<pair<int,int>>> returned_types; 
    283233                         
    284                         return pair<vec,vec>(zeros(1),ones(1)); 
    285                 } 
    286          
    287         } 
    288  
    289  
    290         static list<list<pair<int,int>>> possible_models_recurse(int max_order,int number_of_channels) 
    291         { 
    292                 list<list<pair<int,int>>> created_model_types;           
    293  
    294                 if(max_order == 1) 
    295                 {                        
    296                         for(int channel = 0;channel<number_of_channels;channel++) 
    297                         { 
    298                                 list<pair<int,int>> returned_type; 
    299                                 returned_type.push_back(pair<int,int>(channel,1)); 
    300                                 created_model_types.push_back(returned_type); 
    301                         } 
    302  
    303                         return created_model_types; 
    304                 } 
    305                 else 
    306                 { 
    307                         created_model_types = possible_models_recurse(max_order-1,number_of_channels); 
    308                         list<list<pair<int,int>>> returned_types; 
    309                          
    310                         for(list<list<pair<int,int>>>::iterator model_ref = created_model_types.begin();model_ref!=created_model_types.end();model_ref++) 
     234                        for(set<set<pair<int,int>>>::iterator model_ref = created_model_types.begin();model_ref!=created_model_types.end();model_ref++) 
    311235                        {                                
    312236                                 
    … …  
    315239                                        for(int channel = 0;channel<number_of_channels;channel++) 
    316240                                        { 
    317                                                 list<pair<int,int>> returned_type; 
    318                                                 pair<int,int> new_pair = pair<int,int>(channel,order); 
    319                                                 if(find((*model_ref).begin(),(*model_ref).end(),new_pair)==(*model_ref).end()) 
     241                                                set<pair<int,int>> returned_type; 
     242                                                pair<int,int> new_pair = pair<int,int>(channel,order);//?? 
     243                                                if(find((*model_ref).begin(),(*model_ref).end(),new_pair)==(*model_ref).end()) //?? 
    320244                                                { 
    321                                                         returned_type.insert(returned_type.begin(),(*model_ref).begin(),(*model_ref).end()); 
    322                                                         returned_type.push_back(new_pair); 
    323                                                         returned_types.push_back(returned_type);                                                         
     245                                                        returned_type.insert((*model_ref).begin(),(*model_ref).end()); //co vlozi na zaciatok retuned_type? 
     246                                                        returned_type.insert(new_pair); 
     247                                                         
     248 
     249                                                        returned_types.insert(returned_type);                                                    
    324250                                                } 
    325251                                        } 
    … …  
    327253                        } 
    328254 
    329                         created_model_types.insert(created_model_types.end(),returned_types.begin(),returned_types.end()); 
     255                        created_model_types.insert(returned_types.begin(),returned_types.end()); 
    330256 
    331257                        return created_model_types; 
    … …  
    337263 
    338264 
    339  
    340  
    341  
    342265int main ( int argc, char* argv[] ) {    
    343266         
    344         vec data_vec; 
    345         mat data_matrix; 
     267        /* 
    346268        DWORD Id; 
    347  
    348         char  szApp[]   = "MT4"; 
    349         char  szTopic[] = "BID";         
    350         char  szItem[]  = "EURUSD"; 
    351         char* itRef     = &szItem[0]; 
    352  
    353         char* file_string = "c:\\rtdata"; 
    354  
    355         ofstream myfile; 
    356         char fstring[80];                                        
    357         strcpy(fstring,file_string); 
    358         strcat(fstring,itRef); 
    359         strcat(fstring,"lognc.txt"); 
    360  
    361         list<list<pair<int,int>>> model_types = model::possible_models_recurse(max_model_order,1); 
    362          
    363         int max_window = 50; 
    364         int min_window = 1; 
    365  
    366         list<model*> models; 
    367         for(list<list<pair<int,int>>>::iterator model_type = model_types.begin();model_type!=model_types.end();model_type++) 
    368         { 
    369                 for(int window_size = min_window;window_size < max_window;window_size++) 
    370                 {                        
    371                         //models.push_back(new model((*model_type),true,true,window_size,0,&data_matrix)); 
    372                         models.push_back(new model((*model_type),false,true,window_size,0,&data_matrix));                        
    373                         //models.push_back(new model((*model_type),true,false,window_size,0,&data_matrix)); 
    374                         //models.push_back(new model((*model_type),false,false,window_size,0,&data_matrix));                     
    375                 } 
    376         } 
    377          
    378         mat result_lognc; 
    379         mat result_preds; 
    380          
    381         /* 
    382269        HANDLE hThrd = CreateThread( NULL, 0, (LPTHREAD_START_ROUTINE)ThrdFunc, (LPVOID)1, 0, &Id); 
    383         if ( !hThrd ) 
     270   
     271        if ( !hThrd ) 
    384272    { 
    385273        cout<<"Error Creating Threads,,,,.exiting"<<endl; 
    … …  
    387275    } 
    388276    Sleep ( 100 ); 
    389         */ 
    390          
    391                  
     277 
     278         
     279        char szApp[] = "MT4"; 
     280        char szTopic[] = "QUOTE";        
     281        char szItem1[] = "EURUSD";       
    392282 
    393283        //DDE Initialization 
    394         DWORD idInst = 0; 
     284        DWORD idInst=0; 
    395285        UINT iReturn; 
    396286        iReturn = DdeInitialize(&idInst, (PFNCALLBACK)DdeCallback,  
    … …  
    420310        //Execute commands/requests specific to the DDE Server. 
    421311         
    422         DDERequest(idInst, hConv, szItem);               
    423          
    424         itpp::Laplace_RNG LapRNG = Laplace_RNG(); 
    425  
    426         // time_t current_time = time(0);        
    427  
     312        DDERequest(idInst, hConv, szItem1);              
     313         
    428314        while(1) 
    429315        { 
    430                 while(avg_vec.size() < data_count) 
    431                 { 
    432                         MSG    msg; 
    433                         BOOL   MsgReturn = GetMessage ( &msg , NULL , 0 , 0 ); 
    434                      
    435                         if(MsgReturn) 
    436                         { 
    437                                 TranslateMessage(&msg); 
    438                                 DispatchMessage(&msg);                   
    439                         } 
    440                         Sleep(500); 
    441                 } 
    442  
    443  
    444                 data_vec.ins(data_vec.size(),avg_vec*ones(avg_vec.size())/avg_vec.size()); 
    445                 data_matrix = mat(data_vec).T(); 
    446                 avg_vec.del(0,avg_vec.size()-1);                                 
    447                                          
    448                 vector<double> test_lognc; 
    449                 vec preds; 
    450  
    451                 for(list<model*>::iterator model_ref = models.begin();model_ref!=models.end();model_ref++) 
    452                 { 
    453                         (*model_ref)->data_update(data_vec.size()-1);                    
    454  
    455                         if((*model_ref)->lognc_history.size()==0) 
    456                         { 
    457                                 test_lognc.push_back(0); 
    458                         } 
    459                         else 
    460                         { 
    461                                 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));  
    462                         } 
    463  
    464                         pair<vec,vec> predictions = (*model_ref)->predict(500,data_vec.size()-1,&LapRNG); 
    465  
    466                         preds.ins(preds.size(),(predictions.first*predictions.second)/(predictions.first*ones(predictions.first.size()))); 
    467                          
    468                 } 
    469  
    470                 //preds.ins(0,data_matrix.get(0,data_vec.size()));       
    471                  
    472                 result_preds.ins_col(result_preds.cols(),preds);                                                         
    473  
    474                 myfile.open(fstring,ios::app); 
    475                  
    476                 // myfile << my_rarx->posterior->minimal_vertex->get_coordinates()[0]; 
    477                  
    478                 for(int i = 0;i<test_lognc.size();i++) 
    479                 { 
    480                         myfile << test_lognc[i] << ','; 
    481                 } 
    482  
    483                 myfile << endl;                          
    484                  
    485                 myfile.close(); 
    486                  
    487         } 
    488          
     316                MSG    msg; 
     317                BOOL   MsgReturn = GetMessage ( &msg , NULL , 0 , 0 ); 
     318             
     319                if(MsgReturn) 
     320                { 
     321                        TranslateMessage(&msg); 
     322                        DispatchMessage(&msg);                   
     323                } 
     324        } 
     325 
    489326        //DDE Disconnect and Uninitialize. 
    490327        DdeDisconnect(hConv); 
    491328        DdeUninitialize(idInst); 
    492  
    493         return 0; 
    494 } 
    495  
    496                 /* 
    497                 // EXPERIMENT: A moving window estimation and prediction of RARX is tested on data generated from  
    498                 // 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 
    499                 // can be compared to the classical setup. 
    500                  
    501                 vector<vector<string>> strings; 
    502  
    503                 char* file_string = "c:\\dataGCClosePercDiff";  
    504  
    505                 char dfstring[80]; 
    506                 strcpy(dfstring,file_string); 
    507                 strcat(dfstring,".txt"); 
    508                          
    509                  
    510                 mat data_matrix; 
    511                 ifstream myfile(dfstring); 
     329        */ 
     330 
     331         
     332 
     333        /* 
     334        // 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,  
     335        // where e_t is normally, student(4) and cauchy distributed are tested using robust AR model, to obtain the  
     336        // variance of location parameter estimators and compare it to the classical setup. 
     337        vector<vector<vector<string>>> string_lists; 
     338        string_lists.push_back(vector<vector<string>>()); 
     339        string_lists.push_back(vector<vector<string>>()); 
     340        string_lists.push_back(vector<vector<string>>()); 
     341 
     342        char* file_strings[3] = {"c:\\ar_normal.txt", "c:\\ar_student.txt", "c:\\ar_cauchy.txt"}; 
     343         
     344 
     345        for(int i = 0;i<3;i++) 
     346        {        
     347                ifstream myfile(file_strings[i]); 
    512348                if (myfile.is_open()) 
    513                 {                
    514                         string line; 
    515                         while(getline(myfile,line)) 
    516                         {                        
    517                                 vec data_vector; 
     349                { 
     350                        while ( myfile.good() ) 
     351                        { 
     352                                string line; 
     353                                getline(myfile,line); 
     354                                 
     355                                vector<string> parsed_line; 
    518356                                while(line.find(',') != string::npos) 
    519357                                { 
    520                                         int loc2 = line.find('\n'); 
    521                                         int loc  = line.find(','); 
    522                                         data_vector.ins(data_vector.size(),atof(line.substr(0,loc).c_str()));                            
     358                                        int loc = line.find(','); 
     359                                        parsed_line.push_back(line.substr(0,loc)); 
    523360                                        line.erase(0,loc+1);                                     
     361                                }                                
     362 
     363                                string_lists[i].push_back(parsed_line); 
     364                        } 
     365                        myfile.close(); 
     366                } 
     367        } 
     368 
     369        for(int j = 0;j<string_lists.size();j++) 
     370        {  
     371                 
     372                for(int i = 0;i<string_lists[j].size()-1;i++) 
     373                { 
     374                        vector<vec> conditions; 
     375                        //emlig* emliga = new emlig(2); 
     376                        RARX* my_rarx = new RARX(2,30); 
     377 
     378                        for(int k = 1;k<string_lists[j][i].size();k++) 
     379                        { 
     380                                vec condition; 
     381                                //condition.ins(0,1);                            
     382                                condition.ins(0,string_lists[j][i][k]);                          
     383                                conditions.push_back(condition); 
     384 
     385                                //cout << "orig:" << condition << endl; 
     386 
     387                                if(conditions.size()>1) 
     388                                {                
     389                                        conditions[k-2].ins(0,string_lists[j][i][k]); 
     390                                         
    524391                                } 
    525392 
    526                                 data_matrix.ins_row(data_matrix.rows(),data_vector); 
    527                         }                
    528  
    529                         myfile.close();  
    530                 } 
    531                 else 
    532                 { 
    533                         cout << "Can't open data file!" << endl; 
    534                 } 
    535                 */ 
    536                  
    537                  
    538                 // cout << "Updated." << endl; 
    539          
    540                 /* 
    541                 // 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,  
    542                 // where e_t is normally, student(4) and cauchy distributed are tested using robust AR model, to obtain the  
    543                 // variance of location parameter estimators and compare it to the classical setup. 
    544                 vector<vector<vector<string>>> string_lists; 
    545                 string_lists.push_back(vector<vector<string>>()); 
    546                 string_lists.push_back(vector<vector<string>>()); 
    547                 string_lists.push_back(vector<vector<string>>()); 
    548  
    549                 char* file_strings[3] = {"c:\\ar_normal.txt", "c:\\ar_student.txt", "c:\\ar_cauchy.txt"}; 
    550                  
    551  
    552                 for(int i = 0;i<3;i++) 
    553                 {        
    554                         ifstream myfile(file_strings[i]); 
    555                         if (myfile.is_open()) 
    556                         { 
    557                                 while ( myfile.good() ) 
     393                                if(conditions.size()>2) 
    558394                                { 
    559                                         string line; 
    560                                         getline(myfile,line); 
    561                                          
    562                                         vector<string> parsed_line; 
    563                                         while(line.find(',') != string::npos) 
    564                                         { 
    565                                                 int loc = line.find(','); 
    566                                                 parsed_line.push_back(line.substr(0,loc)); 
    567                                                 line.erase(0,loc+1);                                     
    568                                         }                                
    569  
    570                                         string_lists[i].push_back(parsed_line); 
    571                                 } 
    572                                 myfile.close(); 
    573                         } 
    574                 } 
    575  
    576                 for(int j = 0;j<string_lists.size();j++) 
    577                 {  
    578                          
    579                         for(int i = 0;i<string_lists[j].size()-1;i++) 
    580                         { 
    581                                 vector<vec> conditions; 
    582                                 //emlig* emliga = new emlig(2); 
    583                                 RARX* my_rarx = new RARX(2,30); 
    584  
    585                                 for(int k = 1;k<string_lists[j][i].size();k++) 
    586                                 { 
    587                                         vec condition; 
    588                                         //condition.ins(0,1);                            
    589                                         condition.ins(0,string_lists[j][i][k]);                          
    590                                         conditions.push_back(condition); 
    591  
    592                                         //cout << "orig:" << condition << endl; 
    593  
    594                                         if(conditions.size()>1) 
    595                                         {                
    596                                                 conditions[k-2].ins(0,string_lists[j][i][k]); 
    597                                                  
    598                                         } 
    599  
    600                                         if(conditions.size()>2) 
    601                                         { 
    602                                                 conditions[k-3].ins(0,string_lists[j][i][k]); 
    603  
    604                                                 //cout << "modi:" << conditions[k-3] << endl; 
    605  
    606                                                 my_rarx->bayes(conditions[k-3]); 
    607  
    608                                                  
    609                                                 //if(k>5) 
    610                                                 //{ 
    611                                                 //      cout << "MaxLik coords:" << emliga->minimal_vertex->get_coordinates() << endl; 
    612                                                 //} 
    613                                                  
    614                                         }                                
    615                                          
    616                                 } 
    617  
    618                                 //emliga->step_me(0); 
    619                                 /* 
    620                                 ofstream myfile; 
    621                                 myfile.open("c:\\robust_ar1.txt",ios::app); 
    622                                 myfile << my_rarx->minimal_vertex->get_coordinates()[0] << ";"; 
    623                                 myfile.close(); 
    624  
    625                                 myfile.open("c:\\robust_ar2.txt",ios::app); 
    626                                 myfile << emliga->minimal_vertex->get_coordinates()[1] << ";"; 
    627                                 myfile.close(); 
    628                                  
    629  
    630                                 cout << "MaxLik coords:" << emliga->minimal_vertex->get_coordinates() << endl; 
    631                                 cout << "Step: " << i << endl; 
    632                         } 
    633  
    634                         cout << "One experiment finished." << endl; 
    635  
     395                                        conditions[k-3].ins(0,string_lists[j][i][k]); 
     396 
     397                                        //cout << "modi:" << conditions[k-3] << endl; 
     398 
     399                                        my_rarx->bayes(conditions[k-3]); 
     400 
     401                                         
     402                                        //if(k>5) 
     403                                        //{ 
     404                                        //      cout << "MaxLik coords:" << emliga->minimal_vertex->get_coordinates() << endl; 
     405                                        //} 
     406                                         
     407                                }                                
     408                                 
     409                        } 
     410 
     411                        //emliga->step_me(0); 
     412                        /* 
    636413                        ofstream myfile; 
    637414                        myfile.open("c:\\robust_ar1.txt",ios::app); 
    638                         myfile << endl; 
     415                        myfile << my_rarx->minimal_vertex->get_coordinates()[0] << ";"; 
    639416                        myfile.close(); 
    640417 
    641418                        myfile.open("c:\\robust_ar2.txt",ios::app); 
    642                         myfile << endl; 
     419                        myfile << emliga->minimal_vertex->get_coordinates()[1] << ";"; 
    643420                        myfile.close(); 
    644                 }*/     
    645  
     421                         
     422 
     423                        cout << "MaxLik coords:" << emliga->minimal_vertex->get_coordinates() << endl; 
     424                        cout << "Step: " << i << endl; 
     425                } 
     426 
     427                cout << "One experiment finished." << endl; 
     428 
     429                ofstream myfile; 
     430                myfile.open("c:\\robust_ar1.txt",ios::app); 
     431                myfile << endl; 
     432                myfile.close(); 
     433 
     434                myfile.open("c:\\robust_ar2.txt",ios::app); 
     435                myfile << endl; 
     436                myfile.close(); 
     437        }*/     
     438         
     439        // EXPERIMENT: A moving window estimation and prediction of RARX is tested on data generated from  
     440    // 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 
     441    // can be compared to the classical setup. 
     442         
     443        itpp::Laplace_RNG LapRNG = Laplace_RNG();        
     444 
     445        vector<vector<string>> strings; 
     446 
     447        char* file_string = "c:\\ar_cauchy_single";  
     448 
     449        char dfstring[80]; 
     450        strcpy(dfstring,file_string); 
     451        strcat(dfstring,".txt"); 
     452         
     453         
     454        mat data_matrix; 
     455        ifstream myfile(dfstring); 
     456        if (myfile.is_open()) 
     457        {                
     458                string line; 
     459                while(getline(myfile,line)) 
     460                {                        
     461                        vec data_vector; 
     462                        while(line.find(',') != string::npos) //zmenil som ciarku za medzeru 
     463                        { 
     464                                line.erase(0,1); // toto som sem pridal 
     465                                int loc2 = line.find('\n'); 
     466                                int loc  = line.find(','); 
     467                                data_vector.ins(data_vector.size(),atof(line.substr(0,loc).c_str()));                            
     468                                line.erase(0,loc+1);                                     
     469                        } 
     470 
     471                        data_matrix.ins_row(data_matrix.rows(),data_vector); 
     472                }                
     473 
     474                myfile.close();  
     475        } 
     476        else 
     477        { 
     478                cout << "Can't open data file!" << endl; 
     479        } 
     480         
     481        //konec nacitavania dat 
     482        set<set<pair<int,int>>> model_types = model::possible_models_recurse(max_model_order,data_matrix.rows()); //volanie funkce kde robi kombinace modelov 
     483                                                                                                //to priradime do model_types, data_matrix.row urcuje pocet kanalov dat 
     484        vector<model*> models; 
     485        for(set<set<pair<int,int>>>::iterator model_type = model_types.begin();model_type!=model_types.end();model_type++) 
     486        {// prechadza rozne typy kanalov, a poctu regresorov 
     487                for(int window_size = 100;window_size < 101;window_size++) 
     488                { 
     489                        models.push_back(new model((*model_type),true,true,window_size,0,&data_matrix));   // to su len konstruktory, len inicializujeme rozne typy 
     490                        models.push_back(new model((*model_type),false,true,window_size,0,&data_matrix)); 
     491                        models.push_back(new model((*model_type),true,false,window_size,0,&data_matrix)); 
     492                        models.push_back(new model((*model_type),false,false,window_size,0,&data_matrix));               
     493                } 
     494 
     495                set<pair<int,int>> empty_list; 
     496                models.push_back(new model(empty_list,false,true,100,0,&data_matrix)); 
     497        } 
     498         
     499        mat result_lognc; 
     500        // mat result_preds; 
     501 
     502        for(int time = max_model_order;time<data_matrix.cols();time++) //time<data_matrix.cols()  
     503        {       //pocet stlpcov data_matrix je pocet casovych krokov 
     504                vec cur_res_lognc; 
     505                // vec preds; 
     506                vector<string> nazvy; 
     507                for(vector<model*>::iterator model_ref = models.begin();model_ref!=models.end();model_ref++) 
     508                {//posuvam s apo models, co je pole modelov urobene o cyklus vyssie. Teda som v case time a robim to tam pre vsetky typy modelov, kombinace regresorov 
     509                        (*model_ref)->data_update(time); //pozret sa preco je toto tu nutne 
     510                        //if (time = max_model_order) nazvy.push_back(models.model_ref]);// ako by som mohol dostat nazov modelu? 
     511 
     512                        if((*model_ref)->my_rarx!=NULL) //vklada normalizacnz faktor do cur_res_lognc 
     513                        { 
     514                                cur_res_lognc.ins(cur_res_lognc.size(),(*model_ref)->my_rarx->posterior->log_nc); 
     515                        } 
     516                        else 
     517                        { 
     518                                cur_res_lognc.ins(cur_res_lognc.size(),(*model_ref)->my_arx->posterior().lognc()); 
     519                        } 
     520 
     521                        // pair<vec,vec> predictions = (*model_ref)->predict(200,time,&LapRNG); 
     522 
     523                        // preds.ins(preds.size(),(predictions.first*predictions.second)/(predictions.first*ones(predictions.first.size()))); 
     524                        // preds.ins(0,data_matrix.get(0,time+1)); 
     525 
     526                } 
     527 
     528                result_lognc.ins_col(result_lognc.cols(),cur_res_lognc); 
     529                // result_preds.ins_col(result_preds.cols(),preds); 
     530 
     531                // cout << "Updated." << endl; 
     532         
    646533                /* 
    647534                vector<vec> conditions; 
    … …  
    855742 
    856743                                         
    857                                          
     744                                        ofstream myfile; 
     745                                        char fstring[80];                                        
     746                                        strcpy(fstring,file_string); 
     747                                         
     748                                        strcat(fstring,"lognc.txt");                                     
     749 
     750                                        myfile.open(fstring,ios::app); 
     751                                         
     752                                        // myfile << my_rarx->posterior->minimal_vertex->get_coordinates()[0]; 
     753                                         
     754                                        if(time == max_model_order) 
     755                                        {  
     756                                                for(int i = 0;i<cur_res_lognc.size();i++) 
     757                                                { 
     758                                                        for(set<pair<int,int>>::iterator ar_ref = models[i]->ar_components.begin();ar_ref != models[i]->ar_components.end();ar_ref++) 
     759                                                        { 
     760                                                                myfile << (*ar_ref).second << (*ar_ref).first;                                                   
     761                                                        } 
     762 
     763                                                        myfile << "."; 
     764 
     765                                                        if(models[i]->my_arx == NULL) 
     766                                                        { 
     767                                                                myfile << "1"; 
     768                                                        } 
     769                                                        else 
     770                                                        { 
     771                                                                myfile << "0"; 
     772                                                        } 
     773                                                                 
     774                                                        if(models[i]->has_constant) 
     775                                                        { 
     776                                                                myfile << "1"; 
     777                                                        } 
     778                                                        else 
     779                                                        { 
     780                                                                myfile << "0"; 
     781                                                        } 
     782 
     783                                                        myfile << ","; 
     784                                                } 
     785 
     786                                                myfile << endl; 
     787                                        } 
     788                                         
     789                                        for(int i = 0;i<cur_res_lognc.size();i++) 
     790                                        { 
     791                                                myfile << cur_res_lognc[i] << ' ';//zmenil som ciarku ze medzeru 
     792                                        } 
     793 
     794                                        myfile << endl;                          
     795                                         
     796                                        myfile.close(); 
     797                        } 
    858798                                        /* 
    859799                                        myfile.open(f2string,ios::app); 
    … …  
    10601000        */ 
    10611001 
    1062          
    1063  
    1064  
     1002        return 0; 
     1003} 
     1004 
     1005 
     1006 

  • applications/robust/robustlib.cpp

    r1367 r1376  
    121121                        }                                                
    122122 
    123                         toprow* as_toprow = (toprow*)this; 
    124  
    125                         vec cur_condition = as_toprow->condition_sum.get(1,as_toprow->condition_sum.size()-1); 
    126  
    127                         // cout << as_toprow->condition << endl;                         
     123                        toprow* as_toprow = (toprow*)this;                                               
    128124 
    129125                        int dimension = simplex->vertices.size()-1; 
    … …  
    131127                        mat jacobian(dimension,dimension);                       
    132128 
    133                         double a_0; 
    134129                        map<double,int> as; 
    135                         vertex* base_vertex; 
    136                         set<vertex*>::iterator base_ref = simplex->vertices.begin(); 
    137                         bool order_correct; 
    138                                                  
    139  
    140                         do 
    141                         { 
    142                                 order_correct = true; 
    143                                 as.clear(); 
    144  
    145                                 base_vertex = (*base_ref); 
    146                          
    147                                 /* 
    148                                 cout << "Base coords: " << base_vertex->get_coordinates() << endl; 
    149                                 cout << "Condition: " << as_toprow->condition_sum << endl; 
    150                                 pause(0.1); 
    151                                 */ 
    152  
    153                                 a_0 = -base_vertex->get_coordinates()*cur_condition+as_toprow->condition_sum[0];                 
     130                        vertex* base_vertex = simplex->vertices.begin(); 
     131                                                         
    154132                                 
    155                                 int row_count = 0; 
    156                                 for(set<vertex*>::iterator vert_ref = simplex->vertices.begin(); vert_ref!=simplex->vertices.end();vert_ref++) 
    157                                 { 
    158                                         if(vert_ref != base_ref) 
    159                                         { 
    160                                                 vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates();                                             
    161  
    162                                                 jacobian.set_row(row_count,relative_coords); 
    163  
    164                                                 double new_a = -relative_coords*cur_condition; 
    165  
    166                                                 if(new_a + a_0 == 0 || new_a == 0) 
    167                                                 { 
    168                                                         base_ref++; 
    169  
    170                                                         if(base_ref == simplex->vertices.end()&&simplex->vertices.size()!=2) 
    171                                                         { 
    172                                                                 throw new exception("Equal local conditions are paired. If this ever occurs, the software has to be modified to include multiplied a_0!!"); 
    173                                                         } 
    174                                                         /* 
    175                                                         else if(implex->vertices.size()==2) 
    176                                                         { 
    177                                                                 return relative_coords[0]*exp 
    178                                                         } 
    179                                                         */ 
    180                                                          
    181                                                         order_correct = false;                                           
    182  
    183                                                         break; 
    184                                                 }                                                
    185                                                  
    186                                                  
    187  
    188                                                 if(a_0<current_emlig->min_ll) 
    189                                                 { 
    190                                                         current_emlig->minimal_vertex = base_vertex; 
    191                                                         current_emlig->min_ll = a_0; 
    192                                                 } 
    193  
    194                                                 double a_m = -(*vert_ref)->get_coordinates()*cur_condition+as_toprow->condition_sum[0]; 
    195                                                 if(a_m<current_emlig->min_ll) 
    196                                                 { 
    197                                                         current_emlig->minimal_vertex = (*vert_ref); 
    198                                                         current_emlig->min_ll = a_m;                                             
    199                                                 } 
    200  
    201                                                 //cout << "a0:" << a_0 << " a0 coords:" << base_vertex->get_coordinates() << " am:" << a_m << " am coords:" << (*vert_ref)->get_coordinates() << endl; 
    202  
    203                                                 // cout << "Absolute coords:(V"  << row_count << ")" << (*vert_ref)->get_coordinates() << endl; 
    204                                                 //cout << "Relative coords:(V"  << row_count << ")" << relative_coords << endl; 
    205  
    206                                                 pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(new_a,1)); 
    207                                                 if(returned.second == false) 
    208                                                 { 
    209                                                         (*returned.first).second++; 
    210                                                 } 
    211                                                 /* 
    212                                                 else 
    213                                                 { 
    214                                                         cout << "a[" << row_count << "] = " << new_a << endl; 
    215                                                 } 
    216                                                 */ 
    217                                                  
    218                                                 //as.ins(as.size(),new_a);                                                       
    219                                                  
    220                                                 row_count++; 
    221                                         } 
    222                                 } 
    223                         } 
    224                         while(!order_correct); 
     133                        for(set<vertex*>::iterator vert_ref = simplex->vertices.begin(); vert_ref!=simplex->vertices.end();vert_ref++) 
     134                        { 
     135                                if(vert_ref!=simplex->vertices.begin()) 
     136                                { 
     137                                        vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates(); 
     138                                        jacobian.set_row(row_count,relative_coords); 
     139                                } 
     140 
     141                                double a = -((*vert_ref)->get_coordinates().ins(0,-1)*cur_condition); 
     142                                if(a<current_emlig->min_ll) 
     143                                { 
     144                                        current_emlig->minimal_vertex = (*vert_ref); 
     145                                        current_emlig->min_ll = a;                                               
     146                                } 
     147 
     148                                //cout << "a0:" << a_0 << " a0 coords:" << base_vertex->get_coordinates() << " am:" << a_m << " am coords:" << (*vert_ref)->get_coordinates() << endl; 
     149 
     150                                // cout << "Absolute coords:(V"  << row_count << ")" << (*vert_ref)->get_coordinates() << endl; 
     151                                //cout << "Relative coords:(V"  << row_count << ")" << relative_coords << endl; 
     152 
     153                                pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(a,1)); 
     154                                if(returned.second == false) 
     155                                { 
     156                                        (*returned.first).second++; 
     157                                }                                                
     158                        } 
     159                         
    225160                         
    226161                        /* 
    … …  
    239174 
    240175                        double det_jacobian    = abs(det(jacobian)); 
    241                         double correction_term = det_jacobian;                   
    242                         for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++) 
    243                         { 
    244                                 double multiplier = det_jacobian;                                
     176                        double correction_term; 
     177                        for(map<double,int>::iterator a_ref = as.begin();as_ref!=as.end();as_ref++) 
     178                        { 
     179                                double multiplier = det_jacobian/fact(jacobian.rows());                          
    245180                                 
    246                                 int as1_order = (*as_ref).second; 
    247                                  
    248                                 correction_term /= pow(-(*as_ref).first,as1_order);                                      
    249                                  
    250                                 current_emlig->set_correction_factors(as1_order); 
     181                                int a_order = (*a_ref).second;                   
     182                                current_emlig->set_correction_factors(a_order); 
    251183 
    252184                                vector<double> factors; 
    253185                                int number_of_factors = 0;                                                               
    254                                 for(map<double,int>::iterator as2_ref = as.begin();as2_ref!=as.end();as2_ref++) 
    255                                 { 
    256                                          
    257                                         if(as2_ref!=as_ref) 
     186                                for(map<double,int>::iterator a2_ref = as.begin();a2_ref!=as.end();a2_ref++) 
     187                                {                                        
     188                                        if(a2_ref!=a_ref) 
    258189                                        {                                                                                
    259                                                 for(int k = 0;k<(*as2_ref).second;k++) 
     190                                                for(int k = 0;k<(*a2_ref).second;k++) 
    260191                                                { 
    261                                                         factors.push_back((*as_ref).first-(*as2_ref).first); 
     192                                                        factors.push_back((*a_ref).first-(*a2_ref).first); 
    262193                                                } 
    263194 
    264                                                 multiplier        /= pow((*as_ref).first-(*as2_ref).first,(*as2_ref).second); 
    265                                                 number_of_factors += (*as2_ref).second; 
    266                                         } 
    267                                         else 
    268                                         { 
    269                                                 factors.push_back((*as_ref).first); 
    270  
    271                                                 multiplier        /= (*as_ref).first; 
    272                                                 number_of_factors += 1; 
    273                                         }                                        
     195                                                multiplier        /= pow((*a_ref).first-(*a2_ref).first,(*a2_ref).second); 
     196                                                number_of_factors += (*a2_ref).second; 
     197                                        }                                                                                
    274198                                }        
    275199 
    276                                 double cur_as_ref = (*as_ref).first;                             
     200                                double cur_as_ref = (*a_ref).first;                              
    277201 
    278202                                double gamma_multiplier = -cur_as_ref-a_0; 

  • applications/robust/robustlib.h

    r1367 r1376  
    2222using namespace itpp; 
    2323 
    24 const double max_range = 10;//numeric_limits<double>::max()/10e-10; 
     24const double max_range = 5;//numeric_limits<double>::max()/10e-10; 
    2525 
    2626/// An enumeration of possible actions performed on the polyhedrons. We can merge them or split them. 
    … …  
    753753                        { 
    754754                                vertex_ref++; 
    755                         } 
    756                         while((*vertex_ref)->parentconditions.find(condition)==(*vertex_ref)->parentconditions.end()); 
    757  
    758                         double product = (*vertex_ref)->get_coordinates()*condition->value; 
     755 
     756                                if(vertex_ref==horiz_ref->vertices.end()) 
     757                                { 
     758                                        return; 
     759                                } 
     760                        } 
     761                        while((*vertex_ref)->parentconditions.find(condition)!=(*vertex_ref)->parentconditions.end()); 
     762 
     763                         
     764                         
     765                        vec appended_coords = (*vertex_ref)->get_coordinates(); 
     766                        appended_coords.ins(0,-1.0); 
     767                         
     768                        double product = appended_coords*condition->value; 
    759769 
    760770                        if(should_be_added) 
    … …  
    16771687                        } 
    16781688 
     1689                        for(int i = 1;i<for_merging.size();i++) 
     1690                        { 
     1691                                for(list<polyhedron*>::iterator merge_ref = for_merging[i].begin();merge_ref!=for_merging[i].end();merge_ref++) 
     1692                                { 
     1693                                        delete (*merge_ref); 
     1694                                } 
     1695                        } 
     1696                         
    16791697                        for(set<vertex*>::iterator vert_ref = vertices_to_be_reduced.begin();vert_ref!=vertices_to_be_reduced.end();vert_ref++) 
    16801698                        {