Changeset 1301 for applications/robust
- Timestamp:
- 03/21/11 09:01:57 (14 years ago)
- Location:
- applications/robust
- Files:
-
- 3 modified
Legend:
- Unmodified
- Added
- Removed
-
applications/robust/main.cpp
r1300 r1301 21 21 22 22 /* 23 // 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, 24 // where e_t is normally, student(4) and cauchy distributed are tested using robust AR model, to obtain the 25 // variance of location parameter estimators and compare it to the classical setup. 23 26 vector<vector<vector<string>>> string_lists; 24 27 string_lists.push_back(vector<vector<string>>()); … … 55 58 for(int j = 0;j<string_lists.size();j++) 56 59 { 57 /*60 58 61 for(int i = 0;i<string_lists[j].size()-1;i++) 59 62 { 60 63 vector<vec> conditions; 61 emlig* emliga = new emlig(2); 64 //emlig* emliga = new emlig(2); 65 RARX* my_rarx = new RARX(2,30); 66 62 67 for(int k = 1;k<string_lists[j][i].size();k++) 63 68 { … … 81 86 //cout << "modi:" << conditions[k-3] << endl; 82 87 83 emliga->add_condition(conditions[k-3]);88 my_rarx->bayes(conditions[k-3]); 84 89 85 90 … … 94 99 95 100 //emliga->step_me(0); 101 /* 96 102 ofstream myfile; 97 103 myfile.open("c:\\robust_ar1.txt",ios::app); 98 myfile << emliga->minimal_vertex->get_coordinates()[0] << ";";104 myfile << my_rarx->minimal_vertex->get_coordinates()[0] << ";"; 99 105 myfile.close(); 100 106 … … 102 108 myfile << emliga->minimal_vertex->get_coordinates()[1] << ";"; 103 109 myfile.close(); 110 104 111 105 112 cout << "MaxLik coords:" << emliga->minimal_vertex->get_coordinates() << endl; … … 117 124 myfile << endl; 118 125 myfile.close(); 119 } 120 */ 121 122 123 emlig* emlig1 = new emlig(emlig_size); 124 emlig* emlig2 = new emlig(emlig_size); 126 }*/ 127 128 129 // EXPERIMENT: A moving window estimation and prediction of RARX is tested on data generated from 130 // 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 131 // can be compared to the classical setup. 132 133 134 vector<vector<string>> strings; 135 136 char* file_strings[3] = {"c:\\ar_student_single.txt", "c:\\ar_cauchy_single.txt","c:\\ar_normal_single.txt"}; 137 138 for(int i = 0;i<3;i++) 139 { 140 ifstream myfile(file_strings[i]); 141 if (myfile.is_open()) 142 { 143 string line; 144 getline(myfile,line); 145 146 vector<string> parsed_line; 147 while(line.find(',') != string::npos) 148 { 149 int loc = line.find(','); 150 parsed_line.push_back(line.substr(0,loc)); 151 line.erase(0,loc+1); 152 } 153 154 strings.push_back(parsed_line); 155 156 myfile.close(); 157 } 158 } 159 160 for(int j = 0;j<strings.size();j++) 161 { 162 vector<vec> conditions; 163 //emlig* emliga = new emlig(2); 164 RARX* my_rarx = new RARX(2,0); 165 166 for(int k = 1;k<170;k++) 167 { 168 vec condition; 169 //condition.ins(0,1); 170 condition.ins(0,strings[j][k]); 171 conditions.push_back(condition); 172 173 //cout << "orig:" << condition << endl; 174 175 if(conditions.size()>1) 176 { 177 conditions[k-2].ins(0,strings[j][k]); 178 179 } 180 181 if(conditions.size()>2) 182 { 183 conditions[k-3].ins(0,strings[j][k]); 184 185 // cout << "modi:" << conditions[k-3] << endl; 186 187 my_rarx->bayes(conditions[k-3]); 188 189 190 191 if(k>5) 192 { 193 cout << "MaxLik coords:" << my_rarx->posterior->minimal_vertex->get_coordinates() << endl; 194 195 ofstream myfile; 196 char fstring[80]; 197 strcpy(fstring,file_strings[j]); 198 strcat(fstring,"_res.txt"); 199 200 myfile.open(fstring,ios::app); 201 myfile << my_rarx->posterior->minimal_vertex->get_coordinates()[0]; 202 if(k!=strings[j].size()-1) 203 { 204 myfile << ","; 205 } 206 else 207 { 208 myfile << endl; 209 } 210 myfile.close(); 211 } 212 213 } 214 215 //emliga->step_me(0); 216 /* 217 ofstream myfile; 218 myfile.open("c:\\robust_ar1.txt",ios::app); 219 myfile << my_rarx->minimal_vertex->get_coordinates()[0] << ";"; 220 myfile.close(); 221 222 myfile.open("c:\\robust_ar2.txt",ios::app); 223 myfile << emliga->minimal_vertex->get_coordinates()[1] << ";"; 224 myfile.close(); 225 226 227 cout << "MaxLik coords:" << emliga->minimal_vertex->get_coordinates() << endl; 228 cout << "Step: " << i << endl;*/ 229 } 230 } 231 232 /* 233 cout << "One experiment finished." << endl; 234 235 ofstream myfile; 236 myfile.open("c:\\robust_ar1.txt",ios::app); 237 myfile << endl; 238 myfile.close(); 239 240 myfile.open("c:\\robust_ar2.txt",ios::app); 241 myfile << endl; 242 myfile.close();*/ 243 244 245 //emlig* emlig1 = new emlig(emlig_size); 246 247 //emlig1->step_me(0); 248 //emlig* emlig2 = new emlig(emlig_size); 125 249 126 250 /* … … 142 266 }*/ 143 267 144 268 /* 145 269 vec condition5 = "1.0 1.0 1.01";//"-0.3 1.7 1.5"; 146 270 147 271 emlig1->add_condition(condition5); 148 149 //vec condition1a = "1.0 1.0 1.01"; 272 //emlig1->step_me(0); 273 274 275 vec condition1a = "-1.0 1.02 0.5"; 150 276 //vec condition1b = "1.0 1.0 1.01"; 151 //emlig1->add_condition(condition1a);277 emlig1->add_condition(condition1a); 152 278 //emlig2->add_condition(condition1b); 153 279 154 //vec condition2a = "-1.0 1.0 1.0";280 vec condition2a = "-0.3 1.7 1.5"; 155 281 //vec condition2b = "-1.0 1.0 1.0"; 156 //emlig1->add_condition(condition2a);282 emlig1->add_condition(condition2a); 157 283 //emlig2->add_condition(condition2b); 158 284 159 //vec condition3a = "0.5 -1.01 1.0";285 vec condition3a = "0.5 -1.01 1.0"; 160 286 //vec condition3b = "0.5 -1.01 1.0"; 161 287 162 //emlig1->add_condition(condition3a);288 emlig1->add_condition(condition3a); 163 289 //emlig2->add_condition(condition3b); 164 290 165 //vec condition4a = "-0.5 -1.0 1.0";291 vec condition4a = "-0.5 -1.0 1.0"; 166 292 //vec condition4b = "-0.5 -1.0 1.0"; 167 293 168 //emlig1->add_condition(condition4a);294 emlig1->add_condition(condition4a); 169 295 //cout << "************************************************" << endl; 170 296 //emlig2->add_condition(condition4b); … … 173 299 //cout << emlig1->minimal_vertex->get_coordinates(); 174 300 175 emlig1->remove_condition(condition5); 301 //emlig1->remove_condition(condition3a); 302 //emlig1->step_me(0); 303 //emlig1->remove_condition(condition2a); 304 //emlig1->remove_condition(condition1a); 305 //emlig1->remove_condition(condition5); 306 176 307 177 308 //emlig1->step_me(0); … … 183 314 184 315 //emlig1->remove_condition(condition1); 185 186 187 188 189 /*190 316 317 318 319 320 321 /* 191 322 for(int i = 0;i<100;i++) 192 323 { -
applications/robust/robustlib.cpp
r1300 r1301 4 4 void polyhedron::triangulate(bool should_integrate) 5 5 { 6 t his->triangulation.clear();6 triangulation.clear(); 7 7 8 8 if(should_integrate) … … 11 11 } 12 12 13 if(vertices.size()==1) 14 { 15 set<vertex*> vert_simplex; 16 vert_simplex.insert((vertex*)this); 17 18 triangulation.push_back(vert_simplex); 19 } 20 13 21 for(list<polyhedron*>::iterator child_ref = children.begin();child_ref!=children.end();child_ref++) 14 22 { … … 16 24 { 17 25 set<vertex*> new_simplex; 26 new_simplex.insert((*t_ref).begin(),(*t_ref).end()); 27 28 pair<set<vertex*>::iterator,bool> ret_val = new_simplex.insert(*vertices.begin()); 29 30 if(ret_val.second == true) 31 { 32 triangulation.push_back(new_simplex); 33 34 if(should_integrate) 35 { 36 ((toprow *)this)->probability += ((toprow *)this)->integrate_simplex(new_simplex, 'S'); 37 } 38 } 39 } 40 } 41 42 /* 43 for(list<polyhedron*>::iterator child_ref = children.begin();child_ref!=children.end();child_ref++) 44 { 45 for(list<set<vertex*>>::iterator t_ref = (*child_ref)->triangulation.begin();t_ref!=(*child_ref)->triangulation.end();t_ref++) 46 { 47 set<vertex*> new_simplex; 18 48 new_simplex.insert((*t_ref).begin(),(*t_ref).end()); 19 49 20 50 for(set<vertex*>::iterator suitable_vert_ref = vertices.begin();*suitable_vert_ref<*(*t_ref).begin();suitable_vert_ref++) 21 51 { 22 new_simplex.insert(*suitable_vert_ref); 23 24 triangulation.push_back(new_simplex); 52 set<vertex*> suitable_simplex; 53 suitable_simplex.insert(new_simplex.begin(),new_simplex.end()); 54 55 suitable_simplex.insert(*suitable_vert_ref); 56 57 triangulation.push_back(suitable_simplex); 25 58 26 59 if(should_integrate) 27 60 { 28 ((toprow *)this)->probability += ((toprow *)this)->integrate_simplex( new_simplex, 'S');61 ((toprow *)this)->probability += ((toprow *)this)->integrate_simplex(suitable_simplex, 'S'); 29 62 } 30 63 } 31 64 } 32 } 65 }*/ 33 66 } 34 67 … … 89 122 90 123 91 cout << endl << "Base coords:" << base_vertex->get_coordinates() << endl;124 // cout << endl << "Base coords:" << base_vertex->get_coordinates() << endl; 92 125 93 126 a_0 = base_vertex->get_coordinates()*cur_condition-as_toprow->condition_sum[0]; … … 136 169 //cout << "a0:" << a_0 << " a0 coords:" << base_vertex->get_coordinates() << " am:" << a_m << " am coords:" << (*vert_ref)->get_coordinates() << endl; 137 170 138 cout << "Absolute coords:(V" << row_count << ")" << (*vert_ref)->get_coordinates() << endl;171 // cout << "Absolute coords:(V" << row_count << ")" << (*vert_ref)->get_coordinates() << endl; 139 172 //cout << "Relative coords:(V" << row_count << ")" << relative_coords << endl; 140 173 -
applications/robust/robustlib.h
r1300 r1301 21 21 using namespace itpp; 22 22 23 const double max_range = 10;//numeric_limits<double>::max()/10e-10; 24 23 const double max_range = 100;//numeric_limits<double>::max()/10e-10; 24 25 /// An enumeration of possible actions performed on the polyhedrons. We can merge them or split them. 25 26 enum actions {MERGE, SPLIT}; 26 27 27 28 28 // Forward declaration of polyhedron, vertex and emlig 29 29 class polyhedron; 30 30 class vertex; 31 class emlig; 31 32 32 33 /* … … 45 46 };*/ 46 47 47 class emlig; 48 class polyhedron; 49 48 /// A class representing a single condition that can be added to the emlig. A condition represents data entries in a statistical model. 50 49 class condition 51 50 { 52 51 public: 52 /// Value of the condition representing the data 53 53 vec value; 54 54 55 /// Mulitplicity of the given condition may represent multiple occurences of same data entry. 55 56 int multiplicity; 56 57 58 /// Default constructor of condition class takes the value of data entry and creates a condition with multiplicity 1 (first occurence of the data). 57 59 condition(vec value) 58 60 { … … 72 74 int multiplicity; 73 75 76 /// A property representing the position of the polyhedron related to current condition with relation to which we 77 /// are splitting the parameter space (new data has arrived). This property is setup within a classification procedure and 78 /// is only valid while the new condition is being added. It has to be reset when new condition is added and new classification 79 /// has to be performed. 74 80 int split_state; 75 81 82 /// A property representing the position of the polyhedron related to current condition with relation to which we 83 /// are merging the parameter space (data is being deleted usually due to a moving window model which is more adaptive and 84 /// steps in for the forgetting in a classical Gaussian AR model). This property is setup within a classification procedure and 85 /// is only valid while the new condition is being removed. It has to be reset when new condition is removed and new classification 86 /// has to be performed. 76 87 int merge_state; 77 88 … … 79 90 80 91 public: 92 /// A pointer to the multi-Laplace inverse gamma distribution this polyhedron belongs to. 81 93 emlig* my_emlig; 82 94 … … 90 102 set<vertex*> vertices; 91 103 104 /// The conditions that gave birth to the polyhedron. If some of them is removed, the polyhedron ceases to exist. 92 105 set<condition*> parentconditions; 93 106 … … 101 114 list<polyhedron*> neutralchildren; 102 115 116 /// A set of grandchildren of the polyhedron that when new condition is added lie exactly on the condition hyperplane. These grandchildren 117 /// behave differently from other grandchildren, when the polyhedron is split. New grandchild is not necessarily created on the crossection of 118 /// the polyhedron and new condition. 103 119 set<polyhedron*> totallyneutralgrandchildren; 104 120 121 /// A set of children of the polyhedron that when new condition is added lie exactly on the condition hyperplane. These children 122 /// behave differently from other children, when the polyhedron is split. New child is not necessarily created on the crossection of 123 /// the polyhedron and new condition. 105 124 set<polyhedron*> totallyneutralchildren; 106 125 126 /// Reverse relation to the totallyneutralgrandchildren set is needed for merging of already existing polyhedrons to keep 127 /// totallyneutralgrandchildren list up to date. 107 128 set<polyhedron*> grandparents; 108 129 130 /// Vertices of the polyhedron classified as positive related to an added condition. When the polyhderon is split by the new condition, 131 /// these vertices will belong to the positive part of the splitted polyhedron. 109 132 set<vertex*> positiveneutralvertices; 110 133 134 /// Vertices of the polyhedron classified as negative related to an added condition. When the polyhderon is split by the new condition, 135 /// these vertices will belong to the negative part of the splitted polyhedron. 111 136 set<vertex*> negativeneutralvertices; 112 137 138 /// A bool specifying if the polyhedron lies exactly on the newly added condition or not. 113 139 bool totally_neutral; 114 140 141 /// When two polyhedrons are merged, there always exists a child lying on the former border of the polyhedrons. This child manages the merge 142 /// of the two polyhedrons. This property gives us the address of the mediator child. 115 143 polyhedron* mergechild; 116 144 145 /// If the polyhedron serves as a mergechild for two of its parents, we need to have the address of the parents to access them. This 146 /// is the pointer to the positive parent being merged. 117 147 polyhedron* positiveparent; 118 148 149 /// If the polyhedron serves as a mergechild for two of its parents, we need to have the address of the parents to access them. This 150 /// is the pointer to the negative parent being merged. 119 151 polyhedron* negativeparent; 120 152 153 /// Adressing withing the statistic. Next_poly is a pointer to the next polyhedron in the statistic on the same level (if this is a point, 154 /// next_poly will be a point etc.). 121 155 polyhedron* next_poly; 122 156 157 /// Adressing withing the statistic. Prev_poly is a pointer to the previous polyhedron in the statistic on the same level (if this is a point, 158 /// next_poly will be a point etc.). 123 159 polyhedron* prev_poly; 124 160 161 /// A property counting the number of messages obtained from children within a classification procedure of position of the polyhedron related 162 /// an added/removed condition. If the message counter reaches the number of children, we know the polyhedrons' position has been fully classified. 125 163 int message_counter; 126 164 … … 173 211 174 212 175 213 /// A setter of state of current polyhedron relative to the action specified in the argument. The three possible states of the 214 /// polyhedron are -1 - NEGATIVE, 0 - NEUTRAL, 1 - POSITIVE. Neutral state means that either the state has been reset or the polyhedron is 215 /// ready to be split/merged. 176 216 int set_state(double state_indicator, actions action) 177 217 { … … 187 227 } 188 228 229 /// A getter of state of current polyhedron relative to the action specified in the argument. The three possible states of the 230 /// polyhedron are -1 - NEGATIVE, 0 - NEUTRAL, 1 - POSITIVE. Neutral state means that either the state has been reset or the polyhedron is 231 /// ready to be split/merged. 189 232 int get_state(actions action) 190 233 { … … 200 243 } 201 244 245 /// Method for obtaining the number of children of given polyhedron. 202 246 int number_of_children() 203 247 { … … 205 249 } 206 250 207 251 /// A method for triangulation of given polyhedron. 208 252 void triangulate(bool should_integrate); 209 253 }; … … 218 262 219 263 public: 220 264 /// A property specifying the value of the density (ted nevim, jestli je to jakoby log nebo ne) above the vertex. 221 265 double function_value; 222 266 … … 255 299 256 300 257 /// A class representing a polyhedron in a top row of the complex. Such polyhedron has a condition that differ itiates301 /// A class representing a polyhedron in a top row of the complex. Such polyhedron has a condition that differen tiates 258 302 /// it from polyhedrons in other rows. 259 303 class toprow : public polyhedron … … 282 326 toprow(vec condition_sum, int condition_order) 283 327 { 284 this->condition_sum = condition_sum; 328 this->condition_sum = condition_sum; 329 this->condition_order = condition_order; 285 330 } 286 331 … … 625 670 double product = (*vertex_ref)->get_coordinates()*condition->value; 626 671 627 if((product>0 && should_be_added)||(product<0 && !should_be_added)) 628 { 629 ((toprow*) horiz_ref)->condition_sum += condition->value; 672 if(should_be_added) 673 { 674 ((toprow*) horiz_ref)->condition_order++; 675 676 if(product>0) 677 { 678 ((toprow*) horiz_ref)->condition_sum += condition->value; 679 } 680 else 681 { 682 ((toprow*) horiz_ref)->condition_sum -= condition->value; 683 } 630 684 } 631 685 else 632 { 633 ((toprow*) horiz_ref)->condition_sum -= condition->value; 686 { 687 ((toprow*) horiz_ref)->condition_order--; 688 689 if(product<0) 690 { 691 ((toprow*) horiz_ref)->condition_sum += condition->value; 692 } 693 else 694 { 695 ((toprow*) horiz_ref)->condition_sum -= condition->value; 696 } 634 697 } 635 698 } … … 641 704 { 642 705 643 bool shouldmerge = (toremove == NULL);644 bool shouldsplit = (toadd == NULL);706 bool shouldmerge = (toremove != NULL); 707 bool shouldsplit = (toadd != NULL); 645 708 646 709 if(shouldsplit||shouldmerge) … … 675 738 if(is_last) 676 739 { 677 if(current_parent->mergechild !=NULL)740 if(current_parent->mergechild != NULL) 678 741 { 679 742 if(current_parent->mergechild->get_multiplicity()==1) … … 688 751 current_parent->mergechild->negativeparent = current_parent; 689 752 } 753 } 754 } 755 else 756 { 757 if(parent_state == 1) 758 { 759 ((toprow*)current_parent)->condition_sum-=toremove->value; 760 ((toprow*)current_parent)->condition_order--; 690 761 } 691 } 692 762 763 if(parent_state == -1) 764 { 765 ((toprow*)current_parent)->condition_sum+=toremove->value; 766 ((toprow*)current_parent)->condition_order--; 767 } 768 769 //current_parent->set_state(0,MERGE); 770 } 693 771 694 772 if(parent_state == 0) 695 773 { 696 for_merging[level+1].push_back(current_parent); 774 for_merging[level+1].push_back(current_parent); 775 // current_parent->parentconditions.erase(toremove); 776 } 777 778 current_parent->mergechild = NULL; 779 current_parent->message_counter = 0; 780 781 if(level == number_of_parameters - 1) 782 { 783 toprow* cur_par_toprow = ((toprow*)current_parent); 784 cur_par_toprow->probability = 0.0; 785 786 for(list<set<vertex*>>::iterator t_ref = current_parent->triangulation.begin();t_ref!=current_parent->triangulation.end();t_ref++) 787 { 788 cur_par_toprow->probability += cur_par_toprow->integrate_simplex(*t_ref,'C'); 789 } 790 } 791 } 792 } 793 794 if(shouldsplit) 795 { 796 current_parent->totallyneutralgrandchildren.insert(sender->totallyneutralchildren.begin(),sender->totallyneutralchildren.end()); 797 798 for(set<polyhedron*>::iterator tot_child_ref = sender->totallyneutralchildren.begin();tot_child_ref!=sender->totallyneutralchildren.end();tot_child_ref++) 799 { 800 (*tot_child_ref)->grandparents.insert(current_parent); 801 } 802 803 switch(sender->get_state(SPLIT)) 804 { 805 case 1: 806 current_parent->positivechildren.push_back(sender); 807 current_parent->positiveneutralvertices.insert(sender->vertices.begin(),sender->vertices.end()); 808 break; 809 case 0: 810 current_parent->neutralchildren.push_back(sender); 811 current_parent->positiveneutralvertices.insert(sender->positiveneutralvertices.begin(),sender->positiveneutralvertices.end()); 812 current_parent->negativeneutralvertices.insert(sender->negativeneutralvertices.begin(),sender->negativeneutralvertices.end()); 813 814 if(current_parent->totally_neutral == NULL) 815 { 816 current_parent->totally_neutral = sender->totally_neutral; 697 817 } 698 818 else 699 819 { 700 current_parent->set_state(0,MERGE); 701 } 702 703 current_parent->mergechild = NULL; 704 } 705 } 706 707 if(shouldsplit) 708 { 709 current_parent->totallyneutralgrandchildren.insert(sender->totallyneutralchildren.begin(),sender->totallyneutralchildren.end()); 710 711 for(set<polyhedron*>::iterator tot_child_ref = sender->totallyneutralchildren.begin();tot_child_ref!=sender->totallyneutralchildren.end();tot_child_ref++) 712 { 713 (*tot_child_ref)->grandparents.insert(current_parent); 714 } 715 716 switch(sender->get_state(SPLIT)) 717 { 718 case 1: 719 current_parent->positivechildren.push_back(sender); 720 current_parent->positiveneutralvertices.insert(sender->vertices.begin(),sender->vertices.end()); 721 break; 722 case 0: 723 current_parent->neutralchildren.push_back(sender); 724 current_parent->positiveneutralvertices.insert(sender->positiveneutralvertices.begin(),sender->positiveneutralvertices.end()); 725 current_parent->negativeneutralvertices.insert(sender->negativeneutralvertices.begin(),sender->negativeneutralvertices.end()); 726 727 if(current_parent->totally_neutral == NULL) 820 current_parent->totally_neutral = current_parent->totally_neutral && sender->totally_neutral; 821 } 822 823 if(sender->totally_neutral) 824 { 825 current_parent->totallyneutralchildren.insert(sender); 826 } 827 828 break; 829 case -1: 830 current_parent->negativechildren.push_back(sender); 831 current_parent->negativeneutralvertices.insert(sender->vertices.begin(),sender->vertices.end()); 832 break; 833 } 834 835 if(is_last) 836 { 837 if((current_parent->negativechildren.size()>0&¤t_parent->positivechildren.size()>0)|| 838 (current_parent->neutralchildren.size()>0&¤t_parent->totally_neutral==false)) 839 { 840 841 for_splitting[level+1].push_back(current_parent); 842 843 current_parent->set_state(0, SPLIT); 844 } 845 else 846 { 847 848 849 if(current_parent->negativechildren.size()>0) 728 850 { 729 current_parent->totally_neutral = sender->totally_neutral; 851 current_parent->set_state(-1, SPLIT); 852 853 ((toprow*)current_parent)->condition_sum-=toadd->value; 854 855 856 } 857 else if(current_parent->positivechildren.size()>0) 858 { 859 current_parent->set_state(1, SPLIT); 860 861 ((toprow*)current_parent)->condition_sum+=toadd->value; 730 862 } 731 863 else 732 864 { 733 current_parent-> totally_neutral = current_parent->totally_neutral && sender->totally_neutral;865 current_parent->raise_multiplicity(); 734 866 } 735 867 736 if(sender->totally_neutral) 868 ((toprow*)current_parent)->condition_order++; 869 870 if(level == number_of_parameters - 1) 737 871 { 738 current_parent->totallyneutralchildren.insert(sender); 872 toprow* cur_par_toprow = ((toprow*)current_parent); 873 cur_par_toprow->probability = 0.0; 874 875 for(list<set<vertex*>>::iterator t_ref = current_parent->triangulation.begin();t_ref!=current_parent->triangulation.end();t_ref++) 876 { 877 cur_par_toprow->probability += cur_par_toprow->integrate_simplex(*t_ref,'C'); 878 } 739 879 } 740 741 break; 742 case -1: 743 current_parent->negativechildren.push_back(sender); 744 current_parent->negativeneutralvertices.insert(sender->vertices.begin(),sender->vertices.end()); 745 break; 746 } 747 748 if(is_last) 749 { 750 if((current_parent->negativechildren.size()>0&¤t_parent->positivechildren.size()>0)|| 751 (current_parent->neutralchildren.size()>0&¤t_parent->totally_neutral==false)) 752 { 753 754 for_splitting[level+1].push_back(current_parent); 755 756 current_parent->set_state(0, SPLIT); 757 } 758 else 759 { 760 761 762 if(current_parent->negativechildren.size()>0) 763 { 764 current_parent->set_state(-1, SPLIT); 765 766 ((toprow*)current_parent)->condition_sum-=toadd->value; 767 768 769 } 770 else if(current_parent->positivechildren.size()>0) 771 { 772 current_parent->set_state(1, SPLIT); 773 774 ((toprow*)current_parent)->condition_sum+=toadd->value; 775 } 776 else 777 { 778 current_parent->raise_multiplicity(); 779 } 780 781 ((toprow*)current_parent)->condition_order++; 782 783 if(level == number_of_parameters - 1) 784 { 785 toprow* cur_par_toprow = ((toprow*)current_parent); 786 cur_par_toprow->probability = 0.0; 787 788 for(list<set<vertex*>>::iterator t_ref = current_parent->triangulation.begin();t_ref!=current_parent->triangulation.end();t_ref++) 789 { 790 cur_par_toprow->probability += cur_par_toprow->integrate_simplex(*t_ref,'C'); 791 } 792 } 793 794 current_parent->positivechildren.clear(); 795 current_parent->negativechildren.clear(); 796 current_parent->neutralchildren.clear(); 797 current_parent->totallyneutralchildren.clear(); 798 current_parent->totallyneutralgrandchildren.clear(); 799 current_parent->positiveneutralvertices.clear(); 800 current_parent->negativeneutralvertices.clear(); 801 current_parent->totally_neutral = NULL; 802 current_parent->kids_rel_addresses.clear(); 803 current_parent->message_counter = 0; 804 } 805 } 806 } 807 808 if(is_last) 809 { 810 send_state_message(current_parent,toadd,toremove,level+1); 811 } 880 881 current_parent->positivechildren.clear(); 882 current_parent->negativechildren.clear(); 883 current_parent->neutralchildren.clear(); 884 current_parent->totallyneutralchildren.clear(); 885 current_parent->totallyneutralgrandchildren.clear(); 886 // current_parent->grandparents.clear(); 887 current_parent->positiveneutralvertices.clear(); 888 current_parent->negativeneutralvertices.clear(); 889 current_parent->totally_neutral = NULL; 890 current_parent->kids_rel_addresses.clear(); 891 current_parent->message_counter = 0; 892 } 893 } 894 } 895 896 if(is_last) 897 { 898 send_state_message(current_parent,toadd,toremove,level+1); 899 } 812 900 813 901 } … … 849 937 void step_me(int marker) 850 938 { 851 /*939 852 940 for(int i = 0;i<statistic.size();i++) 853 941 { 854 942 for(polyhedron* horiz_ref = statistic.rows[i];horiz_ref!=statistic.get_end();horiz_ref=horiz_ref->next_poly) 855 943 { 944 856 945 if(i==statistic.size()-1) 857 946 { 858 cout << ((toprow*)horiz_ref)->condition<< " " << ((toprow*)horiz_ref)->probability << endl;947 //cout << ((toprow*)horiz_ref)->condition_sum << " " << ((toprow*)horiz_ref)->probability << endl; 859 948 cout << "Order:" << ((toprow*)horiz_ref)->condition_order << endl; 860 949 } 950 if(i==0) 951 { 952 cout << ((vertex*)horiz_ref)->get_coordinates() << endl; 953 } 954 861 955 char* string = "Checkpoint"; 862 956 } 863 957 } 864 */958 865 959 866 960 /* … … 954 1048 955 1049 list<condition*>::iterator toremove_ref = conditions.end(); 956 bool condition_should_be_added = false;1050 bool condition_should_be_added = should_add; 957 1051 958 1052 for(list<condition*>::iterator ref = conditions.begin();ref!=conditions.end();ref++) … … 986 1080 987 1081 should_add = false; 988 } 989 else 990 { 991 condition_should_be_added = true; 992 } 993 } 994 } 1082 1083 condition_should_be_added = false; 1084 } 1085 } 1086 } 995 1087 996 1088 condition* condition_to_remove = NULL; … … 998 1090 if(toremove_ref!=conditions.end()) 999 1091 { 1000 conditions.erase(toremove_ref);1001 1092 condition_to_remove = *toremove_ref; 1093 conditions.erase(toremove_ref); 1002 1094 } 1003 1095 … … 1006 1098 if(condition_should_be_added) 1007 1099 { 1008 conditions.push_back(new condition(toadd)); 1009 } 1010 1011 1100 condition* new_condition = new condition(toadd); 1101 1102 conditions.push_back(new_condition); 1103 condition_to_add = new_condition; 1104 } 1012 1105 1013 1106 for(polyhedron* horizontal_position = statistic.rows[0];horizontal_position!=statistic.get_end();horizontal_position=horizontal_position->next_poly) … … 1117 1210 { 1118 1211 (*child_ref)->parents.remove(current_negative); 1119 (*child_ref)->parents.push_back(current_positive); 1212 (*child_ref)->parents.push_back(current_positive); 1120 1213 } 1121 1214 … … 1140 1233 1141 1234 current_positive->parents.insert(current_positive->parents.begin(),current_negative->parents.begin(),current_negative->parents.end()); 1142 unique(current_positive->parents.begin(),current_positive->parents.end());1235 // unique(current_positive->parents.begin(),current_positive->parents.end()); 1143 1236 1144 1237 for(list<polyhedron*>::iterator parent_ref = current_negative->parents.begin();parent_ref!=current_negative->parents.end();parent_ref++) … … 1184 1277 { 1185 1278 (*grand_ref)->totallyneutralgrandchildren.erase(current_positive); 1186 } 1187 1188 current_positive->grandparents.clear(); 1279 } 1189 1280 } 1190 1281 else … … 1194 1285 (*grand_ref)->totallyneutralgrandchildren.erase(current_negative); 1195 1286 (*grand_ref)->totallyneutralgrandchildren.insert(current_positive); 1196 } 1197 1198 current_positive->grandparents.insert(current_negative->grandparents.begin(),current_negative->grandparents.end()); 1287 } 1199 1288 } 1200 1289 } … … 1210 1299 } 1211 1300 1301 current_positive->grandparents.clear(); 1302 1212 1303 current_positive->totally_neutral = (current_positive->totally_neutral && current_negative->totally_neutral); 1213 1304 … … 1236 1327 } 1237 1328 1329 1238 1330 for(set<polyhedron*>::iterator grand_p_ref = (*merge_ref)->grandparents.begin();grand_p_ref!=(*merge_ref)->grandparents.end();grand_p_ref++) 1239 1331 { 1240 1332 (*grand_p_ref)->totallyneutralgrandchildren.erase(*merge_ref); 1241 1333 } 1334 1242 1335 1243 1336 statistic.delete_polyhedron(k-1,*merge_ref); … … 1315 1408 } 1316 1409 1410 new_totally_neutral_child->parentconditions.insert(current_polyhedron->parentconditions.begin(),current_polyhedron->parentconditions.end()); 1411 new_totally_neutral_child->parentconditions.insert(condition_to_add); 1412 1317 1413 new_totally_neutral_child->my_emlig = this; 1318 1414 … … 1334 1430 { 1335 1431 (*grand_ref)->parents.push_back(new_totally_neutral_child); 1336 (*grand_ref)->grandparents.insert(positive_poly); 1337 (*grand_ref)->grandparents.insert(negative_poly); 1432 1433 // tohle tu nebylo. ma to tu byt? 1434 //positive_poly->totallyneutralgrandchildren.insert(*grand_ref); 1435 //negative_poly->totallyneutralgrandchildren.insert(*grand_ref); 1436 1437 //(*grand_ref)->grandparents.insert(positive_poly); 1438 //(*grand_ref)->grandparents.insert(negative_poly); 1338 1439 1339 1440 new_totally_neutral_child->vertices.insert((*grand_ref)->vertices.begin(),(*grand_ref)->vertices.end()); … … 1347 1448 { 1348 1449 (*parent_ref)->totallyneutralgrandchildren.insert(new_totally_neutral_child); 1349 new_totally_neutral_child->grandparents.insert(*parent_ref);1450 // new_totally_neutral_child->grandparents.insert(*parent_ref); 1350 1451 1351 1452 (*parent_ref)->neutralchildren.remove(current_polyhedron); … … 1421 1522 { 1422 1523 sizevector.push_back(statistic.row_size(s)); 1423 } 1524 cout << statistic.row_size(s) << ", "; 1525 } 1526 1527 cout << endl; 1424 1528 1425 1529 /* … … 1489 1593 void create_statistic(int number_of_parameters) 1490 1594 { 1595 /* 1491 1596 for(int i = 0;i<number_of_parameters;i++) 1492 1597 { … … 1498 1603 conditions.push_back(new_condition); 1499 1604 } 1605 */ 1500 1606 1501 1607 // An empty vector of coordinates. … … 1791 1897 private: 1792 1898 1899 1900 1901 int window_size; 1902 1903 list<vec> conditions; 1904 1905 public: 1793 1906 emlig* posterior; 1794 1907 1795 int window_size;1796 1797 public:1798 1908 RARX(int number_of_parameters, const int window_size)//:BM() 1799 1909 { 1800 1910 posterior = new emlig(number_of_parameters); 1801 1911 1802 this->window_size = window_size; 1912 this->window_size = window_size; 1803 1913 }; 1804 1914 1805 1915 void bayes(const itpp::vec &yt, const itpp::vec &cond = "") 1806 1916 { 1917 conditions.push_back(yt); 1918 1919 //posterior->step_me(0); 1807 1920 1921 if(conditions.size()>window_size && window_size!=0) 1922 { 1923 posterior->add_and_remove_condition(yt,conditions.front()); 1924 conditions.pop_front(); 1925 1926 //posterior->step_me(1); 1927 } 1928 else 1929 { 1930 posterior->add_condition(yt); 1931 } 1932 1808 1933 } 1809 1934