Changeset 1268 for applications/robust
- Timestamp:
- 01/05/11 18:01:55 (13 years ago)
- Location:
- applications/robust
- Files:
-
- 3 modified
-
main.cpp (modified) (4 diffs)
-
robustlib.cpp (modified) (2 diffs)
-
robustlib.h (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
applications/robust/main.cpp
r1267 r1268 13 13 using namespace bdm; 14 14 15 const int emlig_size = 3; 16 15 17 int main ( int argc, char* argv[] ) { 16 18 17 emlig* emlig1 = new emlig(4);18 19 20 emlig* emlig1 = new emlig(emlig_size); 21 19 22 /* 20 23 emlig1->set_correction_factors(4); … … 24 27 for(set<my_ivec>::iterator vec_ref = emlig1->correction_factors[j].begin();vec_ref!=emlig1->correction_factors[j].end();vec_ref++) 25 28 { 29 cout << j << " "; 30 26 31 for(int i=0;i<(*vec_ref).size();i++) 27 32 { … … 31 36 cout << endl; 32 37 } 33 } 34 */38 }*/ 39 35 40 36 41 //emlig1->step_me(0); … … 40 45 cout << "Step:" << i << endl; 41 46 42 double condition[ 4];47 double condition[emlig_size+1]; 43 48 44 condition[0] = rand()/1000.0; 45 condition[1] = rand()/1000.0; 46 condition[2] = rand()/1000.0; 47 condition[3] = rand()/1000.0; 48 /*condition[4] = rand()/1000.0; 49 condition[5] = rand()/1000.0; 50 condition[6] = rand()/1000.0; 51 condition[7] = rand()/1000.0; 52 condition[8] = rand()/1000.0;*/ 53 49 for(int k = 0;k<=emlig_size;k++) 50 { 51 condition[k] = rand()/1000.0; 52 } 53 54 54 55 vec* condition_vec = new vec(condition, 4);55 vec* condition_vec = new vec(condition,emlig_size+1); 56 56 emlig1->add_condition(*condition_vec); 57 57 58 58 59 emlig1->step_me(i); -
applications/robust/robustlib.cpp
r1262 r1268 20 20 if(should_integrate) 21 21 { 22 emlig* current_emlig; 23 24 if(this->my_emlig!=NULL) 25 { 26 current_emlig = this->my_emlig; 27 } 28 else 29 { 30 throw exception("The statistic of the polyhedron you are trying to integrate over doesn't belong to any emlig!"); 31 } 32 22 33 toprow* as_toprow = (toprow*)this; 23 34 … … 53 64 54 65 double int_value = 0; 66 double det_jacobian = det(jacobian); 55 67 for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++) 56 68 { 57 int fac_order = ((toprow*)this)->condition_order-as.size()-1;69 int condition_order = ((toprow*)this)->condition_order-1; 58 70 59 double fac_value = 1;71 double multiplier = det_jacobian; 60 72 if(a_0!=(*as_ref).first) 61 { 62 //fac_value = ((double)tgamma(fac_order)*det(jacobian))/(*as_ref).first/pow((*as_ref).first-a_0,fac_order); 63 64 fac_value /= (*as_ref).first; 73 { 74 int as1_order = (*as_ref).second; 65 75 76 current_emlig->set_correction_factors(as1_order); 77 78 vector<double> factors; 79 int number_of_factors = 0; 66 80 for(map<double,int>::iterator as2_ref = as.begin();as2_ref!=as.end();as2_ref++) 67 81 { 82 68 83 if(as2_ref!=as_ref) 84 { 85 for(int k = 0;k<(*as2_ref).second;k++) 86 { 87 factors.push_back((*as_ref).first-(*as2_ref).first); 88 } 89 90 multiplier /= pow((*as_ref).first-(*as2_ref).first,(*as2_ref).second); 91 number_of_factors += (*as2_ref).second; 92 } 93 else 69 94 { 70 double current_factor = (*as_ref).first-(*as2_ref).first;95 factors.push_back((*as_ref).first); 71 96 72 fac_value /= current_factor; 97 multiplier /= (*as_ref).first; 98 number_of_factors += 1; 99 } 100 } 73 101 74 75 } 76 } 102 double bracket = fact(condition_order-1-number_of_factors)/fact(as1_order-1)/*/pow((*as_ref).first-a_0,condition_order-1-number_of_factors)*/; 103 for(int k = 0;k < as1_order-1;k++) 104 { 105 double bracket_factor = pow((double)-1,k+1)*fact(condition_order-2-number_of_factors-k)/fact(as1_order-2-k)/pow((*as_ref).first-a_0,condition_order-2-number_of_factors-k); 106 107 for(set<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[k].begin();combi_ref!=this->my_emlig->correction_factors[k].end();combi_ref++) 108 { 109 double bracket_combination = 1; 110 for(int j = 0;j<(*combi_ref).size();j++) 111 { 112 bracket_combination /= factors[(*combi_ref)[j]]; 113 } 114 115 bracket+=bracket_factor*bracket_combination; 116 } 117 } 118 119 int_value += multiplier*bracket; 120 77 121 } 78 122 else 79 123 { 80 124 throw new exception("Should this happen? a_i=a_0 in the formula for integrating over a simplex! I think it can happen with 0 probability(when phi*2*v_0=phi*v_i)"); 81 } 82 83 84 int_value += fac_value; 125 } 126 85 127 } 86 128 -
applications/robust/robustlib.h
r1267 r1268 12 12 #include <vector> 13 13 #include <list> 14 #include <map>15 14 #include <set> 16 15 #include <algorithm> 17 #include <string>18 16 19 17 using namespace bdm; … … 21 19 using namespace itpp; 22 20 23 const double max_range = 1000.0;//numeric_limits<double>::max()/10e-10;21 const double max_range = 99999999.0;//numeric_limits<double>::max()/10e-10; 24 22 25 23 enum actions {MERGE, SPLIT}; … … 27 25 class polyhedron; 28 26 class vertex; 29 class toprow;30 27 31 28 /* … … 43 40 } 44 41 };*/ 42 43 class emlig; 44 45 45 46 46 /// A class describing a single polyhedron of the split complex. From a collection of such classes a Hasse diagram … … 60 60 61 61 public: 62 emlig* my_emlig; 63 62 64 /// A list of polyhedrons parents within the Hasse diagram. 63 65 list<polyhedron*> parents; … … 173 175 } 174 176 177 175 178 void triangulate(bool should_integrate); 176 177 179 178 180 … … 245 247 toprow(vec condition, int condition_order) 246 248 { 247 this->condition 249 this->condition = condition; 248 250 this->condition_order = condition_order; 249 } 251 } 250 252 251 253 }; 252 253 254 254 255 class condition … … 268 269 269 270 270 271 271 class c_statistic 272 272 { … … 274 274 polyhedron* start_poly; 275 275 276 public: 276 public: 277 277 278 vector<polyhedron*> rows; 278 279 … … 414 415 } 415 416 }; 416 417 417 418 418 class my_ivec : public ivec … … 557 557 558 558 559 560 559 561 //! Conditional(e) Multicriteria-Laplace-Inverse-Gamma distribution density 560 562 class emlig // : eEF 561 563 { 562 563 564 564 565 /// A statistic in a form of a Hasse diagram representing a complex of convex polyhedrons obtained as a result … … 573 574 574 575 double normalization_factor; 575 576 577 576 578 577 void alter_toprow_conditions(vec condition, bool should_be_added) … … 727 726 else 728 727 { 729 ((toprow*)current_parent)->condition_order++;730 731 728 if(current_parent->negativechildren.size()>0) 732 729 { 733 730 current_parent->set_state(-1, SPLIT); 734 731 735 ((toprow*)current_parent)->condition-=toadd; 736 732 ((toprow*)current_parent)->condition-=toadd; 737 733 } 738 734 else if(current_parent->positivechildren.size()>0) … … 740 736 current_parent->set_state(1, SPLIT); 741 737 742 ((toprow*)current_parent)->condition+=toadd; 738 ((toprow*)current_parent)->condition+=toadd; 743 739 } 744 740 else … … 790 786 emlig(c_statistic statistic) 791 787 { 792 this->statistic = statistic; 788 this->statistic = statistic; 793 789 } 794 790 … … 966 962 vertex* neutral_vertex = new vertex(new_coordinates); 967 963 964 neutral_vertex->my_emlig = this; 965 968 966 new_totally_neutral_child = neutral_vertex; 969 967 } … … 971 969 { 972 970 toprow* neutral_toprow = new toprow(); 971 972 neutral_toprow->my_emlig = this; 973 973 974 974 new_totally_neutral_child = neutral_toprow; … … 989 989 toprow* negative_poly = new toprow(((toprow*)current_polyhedron)->condition-toadd, ((toprow*)current_polyhedron)->condition_order+1); 990 990 991 positive_poly->my_emlig = this; 992 negative_poly->my_emlig = this; 993 991 994 for(list<polyhedron*>::iterator parent_ref = current_polyhedron->parents.begin();parent_ref!=current_polyhedron->parents.end();parent_ref++) 992 995 { … … 1041 1044 negative_poly->vertices.insert(current_polyhedron->negativeneutralvertices.begin(),current_polyhedron->negativeneutralvertices.end()); 1042 1045 negative_poly->vertices.insert(new_totally_neutral_child->vertices.begin(),new_totally_neutral_child->vertices.end()); 1043 1046 1044 1047 new_totally_neutral_child->triangulate(false); 1045 1048 1046 1049 positive_poly->triangulate(k==for_splitting.size()-1); 1047 1050 negative_poly->triangulate(k==for_splitting.size()-1); 1048 1049 statistic.append_polyhedron(k-1, new_totally_neutral_child); 1051 1052 statistic.append_polyhedron(k-1, new_totally_neutral_child); 1050 1053 1051 1054 statistic.insert_polyhedron(k, positive_poly, current_polyhedron); … … 1061 1064 } 1062 1065 1063 /*1066 1064 1067 vector<int> sizevector; 1065 1068 for(int s = 0;s<statistic.size();s++) 1066 1069 { 1067 1070 sizevector.push_back(statistic.row_size(s)); 1068 } */1071 } 1069 1072 1070 1073 … … 1136 1139 // We create an origin - this point will have all the coordinates zero, but now it has an empty vector of coords. 1137 1140 vertex *origin = new vertex(origin_coord); 1141 1142 origin->my_emlig = this; 1138 1143 1139 1144 /* … … 1149 1154 */ 1150 1155 1151 statistic = *(new c_statistic()); 1152 1156 statistic = *(new c_statistic()); 1157 1153 1158 statistic.append_polyhedron(0, origin); 1154 1159 … … 1169 1174 // Now we create the points 1170 1175 vertex* new_point1 = new vertex(origin_coord1); 1171 vertex* new_point2 = new vertex(origin_coord2); 1176 vertex* new_point2 = new vertex(origin_coord2); 1177 1178 new_point1->my_emlig = this; 1179 new_point2->my_emlig = this; 1172 1180 1173 1181 //********************************************************************************************************* … … 1251 1259 // We create a new toprow with the previously specified condition. 1252 1260 toprow* current_copy1 = new toprow(vec1, 0); 1253 toprow* current_copy2 = new toprow(vec2, 0); 1261 toprow* current_copy2 = new toprow(vec2, 0); 1262 1263 current_copy1->my_emlig = this; 1264 current_copy2->my_emlig = this; 1254 1265 1255 1266 // The vertices of the copies will be inherited, because there will be a parent/child relation … … 1371 1382 statistic.append_polyhedron(j+1,new_statistic1->rows[j],new_statistic1->row_ends[j]); 1372 1383 statistic.append_polyhedron(j+1,new_statistic2->rows[j],new_statistic2->row_ends[j]); 1373 } 1374 1375 1384 } 1376 1385 } 1377 1386 … … 1433 1442 1434 1443 1435 1436 1437 1438 1439 1440 1441 1444 #endif //TRAGE_H
