[976] | 1 | /*! |
---|
| 2 | \file |
---|
| 3 | \brief Robust Bayesian auto-regression model |
---|
| 4 | \author Jan Sindelar. |
---|
| 5 | */ |
---|
| 6 | |
---|
| 7 | #ifndef ROBUST_H |
---|
| 8 | #define ROBUST_H |
---|
| 9 | |
---|
| 10 | #include <stat/exp_family.h> |
---|
[1171] | 11 | #include <limits> |
---|
[1172] | 12 | #include <vector> |
---|
| 13 | #include <algorithm> |
---|
[976] | 14 | |
---|
| 15 | using namespace bdm; |
---|
| 16 | using namespace std; |
---|
[1204] | 17 | using namespace itpp; |
---|
[976] | 18 | |
---|
[1172] | 19 | const double max_range = numeric_limits<double>::max()/10e-5; |
---|
[1171] | 20 | |
---|
[1208] | 21 | enum actions {MERGE, SPLIT}; |
---|
| 22 | |
---|
[1171] | 23 | class polyhedron; |
---|
| 24 | class vertex; |
---|
| 25 | |
---|
[1172] | 26 | /// A class describing a single polyhedron of the split complex. From a collection of such classes a Hasse diagram |
---|
| 27 | /// of the structure in the exponent of a Laplace-Inverse-Gamma density will be created. |
---|
[1171] | 28 | class polyhedron |
---|
[976] | 29 | { |
---|
[1172] | 30 | /// A property having a value of 1 usually, with higher value only if the polyhedron arises as a coincidence of |
---|
| 31 | /// more than just the necessary number of conditions. For example if a newly created line passes through an already |
---|
| 32 | /// existing point, the points multiplicity will rise by 1. |
---|
| 33 | int multiplicity; |
---|
[976] | 34 | |
---|
[1204] | 35 | int split_state; |
---|
| 36 | |
---|
| 37 | int merge_state; |
---|
| 38 | |
---|
| 39 | |
---|
| 40 | |
---|
[1172] | 41 | public: |
---|
| 42 | /// A list of polyhedrons parents within the Hasse diagram. |
---|
| 43 | vector<polyhedron*> parents; |
---|
[1171] | 44 | |
---|
[1172] | 45 | /// A list of polyhedrons children withing the Hasse diagram. |
---|
| 46 | vector<polyhedron*> children; |
---|
[1171] | 47 | |
---|
[1172] | 48 | /// All the vertices of the given polyhedron |
---|
| 49 | vector<vertex*> vertices; |
---|
[1171] | 50 | |
---|
[1172] | 51 | /// A list used for storing children that lie in the positive region related to a certain condition |
---|
| 52 | vector<polyhedron*> positivechildren; |
---|
[1171] | 53 | |
---|
[1172] | 54 | /// A list used for storing children that lie in the negative region related to a certain condition |
---|
| 55 | vector<polyhedron*> negativechildren; |
---|
[1171] | 56 | |
---|
[1172] | 57 | /// Children intersecting the condition |
---|
| 58 | vector<polyhedron*> neutralchildren; |
---|
[1171] | 59 | |
---|
[1211] | 60 | vector<polyhedron*> totallyneutralgrandchildren; |
---|
| 61 | |
---|
| 62 | vector<polyhedron*> totallyneutralchildren; |
---|
| 63 | |
---|
[1208] | 64 | bool totally_neutral; |
---|
| 65 | |
---|
[1204] | 66 | vector<polyhedron*> mergechildren; |
---|
| 67 | |
---|
| 68 | polyhedron* positiveparent; |
---|
| 69 | |
---|
| 70 | polyhedron* negativeparent; |
---|
| 71 | |
---|
| 72 | int message_counter; |
---|
| 73 | |
---|
[1172] | 74 | /// List of triangulation polyhedrons of the polyhedron given by their relative vertices. |
---|
| 75 | vector<vector<vertex*>> triangulations; |
---|
[1171] | 76 | |
---|
[1172] | 77 | /// A list of relative addresses serving for Hasse diagram construction. |
---|
[1171] | 78 | vector<int> kids_rel_addresses; |
---|
| 79 | |
---|
[1172] | 80 | /// Default constructor |
---|
[1171] | 81 | polyhedron() |
---|
| 82 | { |
---|
[1204] | 83 | multiplicity = 1; |
---|
| 84 | |
---|
| 85 | message_counter = 0; |
---|
[1208] | 86 | |
---|
| 87 | totally_neutral = NULL; |
---|
[1171] | 88 | } |
---|
| 89 | |
---|
[1172] | 90 | /// Setter for raising multiplicity |
---|
[1204] | 91 | void raise_multiplicity() |
---|
[1171] | 92 | { |
---|
| 93 | multiplicity++; |
---|
| 94 | } |
---|
| 95 | |
---|
[1172] | 96 | /// Setter for lowering multiplicity |
---|
[1204] | 97 | void lower_multiplicity() |
---|
[1171] | 98 | { |
---|
| 99 | multiplicity--; |
---|
| 100 | } |
---|
| 101 | |
---|
[1172] | 102 | /// An obligatory operator, when the class is used within a C++ STL structure like a vector |
---|
| 103 | int operator==(polyhedron polyhedron2) |
---|
| 104 | { |
---|
| 105 | return true; |
---|
| 106 | } |
---|
| 107 | |
---|
| 108 | /// An obligatory operator, when the class is used within a C++ STL structure like a vector |
---|
| 109 | int operator<(polyhedron polyhedron2) |
---|
| 110 | { |
---|
| 111 | return false; |
---|
| 112 | } |
---|
[1204] | 113 | |
---|
[1208] | 114 | |
---|
| 115 | |
---|
[1204] | 116 | void set_state(double state_indicator, actions action) |
---|
| 117 | { |
---|
| 118 | switch(action) |
---|
| 119 | { |
---|
| 120 | case MERGE: |
---|
| 121 | merge_state = (int)sign(state_indicator); |
---|
| 122 | break; |
---|
| 123 | case SPLIT: |
---|
| 124 | split_state = (int)sign(state_indicator); |
---|
| 125 | break; |
---|
| 126 | } |
---|
| 127 | } |
---|
| 128 | |
---|
| 129 | int get_state(actions action) |
---|
| 130 | { |
---|
| 131 | switch(action) |
---|
| 132 | { |
---|
| 133 | case MERGE: |
---|
| 134 | return merge_state; |
---|
| 135 | break; |
---|
| 136 | case SPLIT: |
---|
| 137 | return split_state; |
---|
| 138 | break; |
---|
| 139 | } |
---|
| 140 | } |
---|
| 141 | |
---|
| 142 | int number_of_children() |
---|
| 143 | { |
---|
[1207] | 144 | return children.size(); |
---|
[1204] | 145 | } |
---|
| 146 | |
---|
[1207] | 147 | |
---|
[976] | 148 | }; |
---|
| 149 | |
---|
[1186] | 150 | /// A class for representing 0-dimensional polyhedron - a vertex. It will be located in the bottom row of the Hasse |
---|
| 151 | /// diagram representing a complex of polyhedrons. It has its coordinates in the parameter space. |
---|
[1172] | 152 | class vertex : public polyhedron |
---|
[976] | 153 | { |
---|
[1186] | 154 | /// A dynamic array representing coordinates of the vertex |
---|
[1204] | 155 | vec coordinates; |
---|
[976] | 156 | |
---|
[1208] | 157 | |
---|
[1204] | 158 | |
---|
[976] | 159 | public: |
---|
[1171] | 160 | |
---|
[1204] | 161 | |
---|
| 162 | |
---|
[1186] | 163 | /// Default constructor |
---|
[1171] | 164 | vertex(); |
---|
| 165 | |
---|
[1186] | 166 | /// Constructor of a vertex from a set of coordinates |
---|
[1204] | 167 | vertex(vec coordinates) |
---|
[1171] | 168 | { |
---|
[1172] | 169 | this->coordinates = coordinates; |
---|
[1171] | 170 | } |
---|
| 171 | |
---|
[1186] | 172 | /// A method that widens the set of coordinates of given vertex. It is used when a complex in a parameter |
---|
| 173 | /// space of certain dimension is established, but the dimension is not known when the vertex is created. |
---|
[1171] | 174 | void push_coordinate(double coordinate) |
---|
| 175 | { |
---|
[1204] | 176 | coordinates = concat(coordinates,coordinate); |
---|
[1171] | 177 | } |
---|
| 178 | |
---|
[1186] | 179 | /// A method obtaining the set of coordinates of a vertex. These coordinates are not obtained as a pointer |
---|
| 180 | /// (not given by reference), but a new copy is created (they are given by value). |
---|
[1204] | 181 | vec get_coordinates() |
---|
| 182 | { |
---|
| 183 | return coordinates; |
---|
| 184 | } |
---|
[1172] | 185 | |
---|
[1204] | 186 | |
---|
[1172] | 187 | }; |
---|
[976] | 188 | |
---|
[1186] | 189 | /// A class representing a polyhedron in a top row of the complex. Such polyhedron has a condition that differitiates |
---|
| 190 | /// it from polyhedrons in other rows. |
---|
[1172] | 191 | class toprow : public polyhedron |
---|
[1171] | 192 | { |
---|
[1204] | 193 | |
---|
| 194 | public: |
---|
[1186] | 195 | /// A condition used for determining the function of a Laplace-Inverse-Gamma density resulting from Bayesian estimation |
---|
[1204] | 196 | vec condition; |
---|
[976] | 197 | |
---|
[1186] | 198 | /// Default constructor |
---|
[1171] | 199 | toprow(); |
---|
[976] | 200 | |
---|
[1186] | 201 | /// Constructor creating a toprow from the condition |
---|
[1204] | 202 | toprow(vec condition) |
---|
[1171] | 203 | { |
---|
[1172] | 204 | this->condition = condition; |
---|
[1171] | 205 | } |
---|
| 206 | |
---|
[1172] | 207 | }; |
---|
[1171] | 208 | |
---|
[1204] | 209 | class condition |
---|
| 210 | { |
---|
| 211 | public: |
---|
| 212 | vec value; |
---|
[1171] | 213 | |
---|
[1204] | 214 | int multiplicity; |
---|
[1171] | 215 | |
---|
[1204] | 216 | condition(vec value) |
---|
| 217 | { |
---|
| 218 | this->value = value; |
---|
| 219 | multiplicity = 1; |
---|
| 220 | } |
---|
[1208] | 221 | }; |
---|
[1171] | 222 | |
---|
[1204] | 223 | |
---|
[976] | 224 | //! Conditional(e) Multicriteria-Laplace-Inverse-Gamma distribution density |
---|
[1186] | 225 | class emlig // : eEF |
---|
[1172] | 226 | { |
---|
[976] | 227 | |
---|
[1186] | 228 | /// A statistic in a form of a Hasse diagram representing a complex of convex polyhedrons obtained as a result |
---|
| 229 | /// of data update from Bayesian estimation or set by the user if this emlig is a prior density |
---|
[1172] | 230 | vector<vector<polyhedron*>> statistic; |
---|
[1204] | 231 | |
---|
[1207] | 232 | vector<vector<polyhedron*>> for_splitting; |
---|
| 233 | |
---|
| 234 | vector<vector<polyhedron*>> for_merging; |
---|
| 235 | |
---|
[1204] | 236 | vector<condition*> conditions; |
---|
| 237 | |
---|
| 238 | double normalization_factor; |
---|
| 239 | |
---|
| 240 | void alter_toprow_conditions(vec condition, bool should_be_added) |
---|
| 241 | { |
---|
| 242 | for(vector<polyhedron*>::iterator horiz_ref = statistic[statistic.size()-1].begin();horiz_ref<statistic[statistic.size()-1].end();horiz_ref++) |
---|
| 243 | { |
---|
| 244 | double product = 0; |
---|
| 245 | |
---|
| 246 | vector<vertex*>::iterator vertex_ref = (*horiz_ref)->vertices.begin(); |
---|
| 247 | |
---|
| 248 | do |
---|
| 249 | { |
---|
[1208] | 250 | product = (*vertex_ref)->get_coordinates()*condition; |
---|
[1204] | 251 | } |
---|
[1208] | 252 | while(product == 0); |
---|
[1204] | 253 | |
---|
| 254 | if((product>0 && should_be_added)||(product<0 && !should_be_added)) |
---|
| 255 | { |
---|
| 256 | ((toprow*) (*horiz_ref))->condition += condition; |
---|
| 257 | } |
---|
| 258 | else |
---|
| 259 | { |
---|
| 260 | ((toprow*) (*horiz_ref))->condition -= condition; |
---|
| 261 | } |
---|
| 262 | } |
---|
| 263 | } |
---|
[1207] | 264 | |
---|
| 265 | |
---|
| 266 | void send_state_message(polyhedron* sender, bool shouldsplit, bool shouldmerge, int level) |
---|
[1211] | 267 | { |
---|
| 268 | |
---|
[1207] | 269 | if(shouldsplit||shouldmerge) |
---|
| 270 | { |
---|
| 271 | for(vector<polyhedron*>::iterator parent_iterator = sender->parents.begin();parent_iterator<sender->parents.end();parent_iterator++) |
---|
| 272 | { |
---|
| 273 | polyhedron* current_parent = *parent_iterator; |
---|
| 274 | |
---|
| 275 | current_parent->message_counter++; |
---|
| 276 | |
---|
| 277 | bool is_last = (current_parent->message_counter == current_parent->number_of_children()); |
---|
| 278 | |
---|
| 279 | if(shouldmerge) |
---|
| 280 | { |
---|
| 281 | int child_state = sender->get_state(MERGE); |
---|
| 282 | int parent_state = current_parent->get_state(MERGE); |
---|
| 283 | |
---|
| 284 | if(parent_state == 0) |
---|
| 285 | { |
---|
| 286 | current_parent->set_state(child_state, MERGE); |
---|
| 287 | |
---|
| 288 | if(child_state == 0) |
---|
| 289 | { |
---|
| 290 | current_parent->mergechildren.push_back(sender); |
---|
| 291 | } |
---|
| 292 | } |
---|
| 293 | else |
---|
| 294 | { |
---|
| 295 | if(child_state == 0) |
---|
| 296 | { |
---|
| 297 | if(parent_state > 0) |
---|
| 298 | { |
---|
| 299 | sender->positiveparent = current_parent; |
---|
| 300 | } |
---|
| 301 | else |
---|
| 302 | { |
---|
| 303 | sender->negativeparent = current_parent; |
---|
| 304 | } |
---|
| 305 | } |
---|
| 306 | } |
---|
| 307 | |
---|
| 308 | if(is_last) |
---|
| 309 | { |
---|
| 310 | if(parent_state > 0) |
---|
| 311 | { |
---|
| 312 | for(vector<polyhedron*>::iterator merge_child = current_parent->mergechildren.begin(); merge_child < current_parent->mergechildren.end();merge_child++) |
---|
| 313 | { |
---|
| 314 | (*merge_child)->positiveparent = current_parent; |
---|
| 315 | } |
---|
| 316 | } |
---|
| 317 | |
---|
| 318 | if(parent_state < 0) |
---|
| 319 | { |
---|
| 320 | for(vector<polyhedron*>::iterator merge_child = current_parent->mergechildren.begin(); merge_child < current_parent->mergechildren.end();merge_child++) |
---|
| 321 | { |
---|
| 322 | (*merge_child)->negativeparent = current_parent; |
---|
| 323 | } |
---|
| 324 | } |
---|
| 325 | |
---|
| 326 | if(parent_state == 0) |
---|
| 327 | { |
---|
| 328 | for_merging[level+1].push_back(current_parent); |
---|
| 329 | } |
---|
| 330 | |
---|
| 331 | current_parent->mergechildren.clear(); |
---|
| 332 | } |
---|
| 333 | |
---|
[1208] | 334 | |
---|
| 335 | } |
---|
| 336 | |
---|
| 337 | if(shouldsplit) |
---|
[1207] | 338 | { |
---|
[1211] | 339 | current_parent->totallyneutralgrandchildren.insert(current_parent->totallyneutralgrandchildren.end(),sender->totallyneutralchildren.begin(),sender->totallyneutralchildren.end()); |
---|
| 340 | |
---|
[1207] | 341 | switch(sender->get_state(SPLIT)) |
---|
| 342 | { |
---|
| 343 | case 1: |
---|
| 344 | current_parent->positivechildren.push_back(sender); |
---|
| 345 | break; |
---|
| 346 | case 0: |
---|
| 347 | current_parent->neutralchildren.push_back(sender); |
---|
[1208] | 348 | |
---|
| 349 | if(current_parent->totally_neutral == NULL) |
---|
| 350 | { |
---|
| 351 | current_parent->totally_neutral = sender->totally_neutral; |
---|
| 352 | } |
---|
| 353 | else |
---|
| 354 | { |
---|
| 355 | current_parent->totally_neutral = current_parent->totally_neutral && sender->totally_neutral; |
---|
| 356 | } |
---|
[1211] | 357 | |
---|
| 358 | if(sender->totally_neutral) |
---|
| 359 | { |
---|
| 360 | current_parent->totallyneutralchildren.push_back(sender); |
---|
| 361 | } |
---|
[1208] | 362 | |
---|
[1207] | 363 | break; |
---|
| 364 | case -1: |
---|
| 365 | current_parent->negativechildren.push_back(sender); |
---|
| 366 | break; |
---|
| 367 | } |
---|
| 368 | |
---|
| 369 | if(is_last) |
---|
| 370 | { |
---|
[1211] | 371 | unique(current_parent->totallyneutralgrandchildren.begin(),current_parent->totallyneutralgrandchildren.end()); |
---|
| 372 | |
---|
[1208] | 373 | if((current_parent->negativechildren.size()>0&¤t_parent->positivechildren.size()>0)|| |
---|
| 374 | (current_parent->neutralchildren.size()>0&¤t_parent->totally_neutral==false)) |
---|
[1207] | 375 | { |
---|
[1208] | 376 | |
---|
[1207] | 377 | for_splitting[level+1].push_back(current_parent); |
---|
| 378 | |
---|
| 379 | current_parent->set_state(0, SPLIT); |
---|
[1208] | 380 | } |
---|
[1211] | 381 | else |
---|
[1208] | 382 | { |
---|
[1211] | 383 | if(current_parent->negativechildren.size()>0) |
---|
| 384 | { |
---|
| 385 | current_parent->set_state(-1, SPLIT); |
---|
| 386 | } |
---|
| 387 | else if(current_parent->positivechildren.size()>0) |
---|
| 388 | { |
---|
[1207] | 389 | current_parent->set_state(1, SPLIT); |
---|
[1211] | 390 | } |
---|
| 391 | else |
---|
| 392 | { |
---|
| 393 | current_parent->raise_multiplicity(); |
---|
| 394 | } |
---|
[1207] | 395 | |
---|
| 396 | current_parent->positivechildren.clear(); |
---|
[1211] | 397 | current_parent->negativechildren.clear(); |
---|
[1207] | 398 | current_parent->neutralchildren.clear(); |
---|
[1211] | 399 | current_parent->totallyneutralchildren.clear(); |
---|
| 400 | current_parent->totallyneutralgrandchildren.clear(); |
---|
| 401 | current_parent->totally_neutral = NULL; |
---|
[1207] | 402 | } |
---|
| 403 | } |
---|
| 404 | } |
---|
| 405 | |
---|
| 406 | if(is_last) |
---|
| 407 | { |
---|
| 408 | send_state_message(current_parent,shouldsplit,shouldmerge,level+1); |
---|
| 409 | } |
---|
| 410 | |
---|
| 411 | } |
---|
| 412 | |
---|
| 413 | } |
---|
| 414 | } |
---|
[1171] | 415 | |
---|
| 416 | public: |
---|
[976] | 417 | |
---|
[1186] | 418 | /// A default constructor creates an emlig with predefined statistic representing only the range of the given |
---|
| 419 | /// parametric space, where the number of parameters of the needed model is given as a parameter to the constructor. |
---|
[1171] | 420 | emlig(int number_of_parameters) |
---|
| 421 | { |
---|
[1172] | 422 | create_statistic(number_of_parameters); |
---|
[1207] | 423 | |
---|
| 424 | for(int i = 0;i<statistic.size();i++) |
---|
| 425 | { |
---|
| 426 | vector<polyhedron*> empty_split; |
---|
| 427 | vector<polyhedron*> empty_merge; |
---|
| 428 | |
---|
| 429 | for_splitting.push_back(empty_split); |
---|
| 430 | for_merging.push_back(empty_merge); |
---|
| 431 | } |
---|
[1171] | 432 | } |
---|
| 433 | |
---|
[1186] | 434 | /// A constructor for creating an emlig when the user wants to create the statistic by himself. The creation of a |
---|
| 435 | /// statistic is needed outside the constructor. Used for a user defined prior distribution on the parameters. |
---|
[1172] | 436 | emlig(vector<vector<polyhedron*>> statistic) |
---|
[1171] | 437 | { |
---|
[1172] | 438 | this->statistic = statistic; |
---|
[1171] | 439 | } |
---|
| 440 | |
---|
[1208] | 441 | void add_condition(vec toadd) |
---|
[1204] | 442 | { |
---|
[1208] | 443 | vec null_vector = ""; |
---|
| 444 | |
---|
| 445 | add_and_remove_condition(toadd, null_vector); |
---|
| 446 | } |
---|
| 447 | |
---|
| 448 | void remove_condition(vec toremove) |
---|
| 449 | { |
---|
| 450 | vec null_vector = ""; |
---|
| 451 | |
---|
| 452 | add_and_remove_condition(null_vector, toremove); |
---|
| 453 | |
---|
| 454 | } |
---|
| 455 | |
---|
| 456 | void add_and_remove_condition(vec toadd, vec toremove) |
---|
| 457 | { |
---|
| 458 | bool should_remove = (toremove.size() != 0); |
---|
| 459 | bool should_add = (toadd.size() != 0); |
---|
| 460 | |
---|
[1204] | 461 | vector<condition*>::iterator toremove_ref = conditions.end(); |
---|
| 462 | bool condition_should_be_added = false; |
---|
| 463 | |
---|
| 464 | for(vector<condition*>::iterator ref = conditions.begin();ref<conditions.end();ref++) |
---|
| 465 | { |
---|
[1208] | 466 | if(should_remove) |
---|
[1204] | 467 | { |
---|
| 468 | if((*ref)->value == toremove) |
---|
| 469 | { |
---|
[1208] | 470 | if((*ref)->multiplicity>1) |
---|
[1204] | 471 | { |
---|
[1208] | 472 | (*ref)->multiplicity--; |
---|
[1204] | 473 | |
---|
| 474 | alter_toprow_conditions(toremove,false); |
---|
| 475 | |
---|
[1208] | 476 | should_remove = false; |
---|
[1204] | 477 | } |
---|
| 478 | else |
---|
| 479 | { |
---|
| 480 | toremove_ref = ref; |
---|
| 481 | } |
---|
| 482 | } |
---|
| 483 | } |
---|
| 484 | |
---|
[1208] | 485 | if(should_add) |
---|
[1204] | 486 | { |
---|
[1208] | 487 | if((*ref)->value == toadd) |
---|
[1204] | 488 | { |
---|
[1208] | 489 | (*ref)->multiplicity++; |
---|
[1204] | 490 | |
---|
| 491 | alter_toprow_conditions(toadd,true); |
---|
| 492 | |
---|
[1208] | 493 | should_add = false; |
---|
[1204] | 494 | } |
---|
| 495 | else |
---|
| 496 | { |
---|
| 497 | condition_should_be_added = true; |
---|
| 498 | } |
---|
| 499 | } |
---|
| 500 | } |
---|
| 501 | |
---|
| 502 | if(toremove_ref!=conditions.end()) |
---|
| 503 | { |
---|
| 504 | conditions.erase(toremove_ref); |
---|
| 505 | } |
---|
| 506 | |
---|
| 507 | if(condition_should_be_added) |
---|
| 508 | { |
---|
| 509 | conditions.push_back(new condition(toadd)); |
---|
| 510 | } |
---|
| 511 | |
---|
[1207] | 512 | |
---|
[1204] | 513 | |
---|
| 514 | for(vector<polyhedron*>::iterator horizontal_position = statistic[0].begin();horizontal_position<statistic[0].end();horizontal_position++) |
---|
| 515 | { |
---|
[1208] | 516 | vertex* current_vertex = (vertex*)(*horizontal_position); |
---|
[1204] | 517 | |
---|
[1208] | 518 | if(should_add||should_remove) |
---|
[1204] | 519 | { |
---|
[1208] | 520 | vec appended_vec = current_vertex->get_coordinates(); |
---|
| 521 | appended_vec.ins(0,1.0); |
---|
| 522 | |
---|
| 523 | if(should_add) |
---|
| 524 | { |
---|
[1211] | 525 | double local_condition = toadd*appended_vec; |
---|
| 526 | |
---|
| 527 | current_vertex->set_state(local_condition,SPLIT); |
---|
| 528 | |
---|
| 529 | if(local_condition == 0) |
---|
| 530 | { |
---|
| 531 | current_vertex->totally_neutral = true; |
---|
| 532 | |
---|
| 533 | current_vertex->multiplicity++; |
---|
| 534 | } |
---|
[1208] | 535 | } |
---|
[1204] | 536 | |
---|
[1208] | 537 | if(should_remove) |
---|
| 538 | { |
---|
[1211] | 539 | double local_condition = toremove*appended_vec; |
---|
[1207] | 540 | |
---|
[1211] | 541 | current_vertex->set_state(local_condition,MERGE); |
---|
| 542 | |
---|
| 543 | if(local_condition == 0) |
---|
| 544 | { |
---|
| 545 | for_merging[0].push_back(current_vertex); |
---|
| 546 | } |
---|
| 547 | } |
---|
[1204] | 548 | } |
---|
| 549 | |
---|
[1208] | 550 | send_state_message(current_vertex, should_add, should_remove, 0); |
---|
| 551 | } |
---|
[1207] | 552 | |
---|
[1208] | 553 | for(vector<vector<polyhedron*>>::iterator vert_ref = for_splitting.begin();vert_ref<for_splitting.end();vert_ref++) |
---|
| 554 | { |
---|
| 555 | int original_size = (*vert_ref).size(); |
---|
| 556 | |
---|
| 557 | for(int split_counter = 0;split_counter<original_size;split_counter++) |
---|
| 558 | { |
---|
| 559 | polyhedron* current_polyhedron = (*vert_ref)[original_size-1-split_counter]; |
---|
| 560 | |
---|
[1211] | 561 | |
---|
[1208] | 562 | } |
---|
[1204] | 563 | } |
---|
[1208] | 564 | |
---|
| 565 | |
---|
[1204] | 566 | } |
---|
| 567 | |
---|
[1171] | 568 | protected: |
---|
| 569 | |
---|
[1186] | 570 | /// A method for creating plain default statistic representing only the range of the parameter space. |
---|
[1172] | 571 | void create_statistic(int number_of_parameters) |
---|
[1171] | 572 | { |
---|
[1186] | 573 | // An empty vector of coordinates. |
---|
[1204] | 574 | vec origin_coord; |
---|
[1171] | 575 | |
---|
[1186] | 576 | // We create an origin - this point will have all the coordinates zero, but now it has an empty vector of coords. |
---|
[1172] | 577 | vertex *origin = new vertex(origin_coord); |
---|
[1171] | 578 | |
---|
[1186] | 579 | // It has itself as a vertex. There will be a nice use for this when the vertices of its parents are searched in |
---|
| 580 | // the recursive creation procedure below. |
---|
[1172] | 581 | origin->vertices.push_back(origin); |
---|
[1171] | 582 | |
---|
[1186] | 583 | // As a statistic, we have to create a vector of vectors of polyhedron pointers. It will then represent the Hasse |
---|
| 584 | // diagram. First we create a vector of polyhedrons.. |
---|
[1172] | 585 | vector<polyhedron*> origin_vec; |
---|
[1171] | 586 | |
---|
[1186] | 587 | // ..we fill it with the origin.. |
---|
[1172] | 588 | origin_vec.push_back(origin); |
---|
| 589 | |
---|
[1186] | 590 | // ..and we fill the statistic with the created vector. |
---|
[1171] | 591 | statistic.push_back(origin_vec); |
---|
| 592 | |
---|
[1186] | 593 | // Now we have a statistic for a zero dimensional space. Regarding to how many dimensional space we need to |
---|
| 594 | // describe, we have to widen the descriptional default statistic. We use an iterative procedure as follows: |
---|
[1172] | 595 | for(int i=0;i<number_of_parameters;i++) |
---|
[1171] | 596 | { |
---|
[1186] | 597 | // We first will create two new vertices. These will be the borders of the parameter space in the dimension |
---|
| 598 | // of newly added parameter. Therefore they will have all coordinates except the last one zero. We get the |
---|
| 599 | // right amount of zero cooridnates by reading them from the origin |
---|
[1204] | 600 | vec origin_coord = origin->get_coordinates(); |
---|
[1171] | 601 | |
---|
[1186] | 602 | // And we incorporate the nonzero coordinates into the new cooordinate vectors |
---|
[1204] | 603 | vec origin_coord1 = concat(origin_coord,max_range); |
---|
| 604 | vec origin_coord2 = concat(origin_coord,-max_range); |
---|
[1172] | 605 | |
---|
[1186] | 606 | // Now we create the points |
---|
[1172] | 607 | vertex *new_point1 = new vertex(origin_coord1); |
---|
[1186] | 608 | vertex *new_point2 = new vertex(origin_coord2); |
---|
| 609 | |
---|
| 610 | //********************************************************************************************************* |
---|
| 611 | // The algorithm for recursive build of a new Hasse diagram representing the space structure from the old |
---|
| 612 | // diagram works so that you create two copies of the old Hasse diagram, you shift them up one level (points |
---|
| 613 | // will be segments, segments will be areas etc.) and you connect each one of the original copied polyhedrons |
---|
| 614 | // with its offspring by a parent-child relation. Also each of the segments in the first (second) copy is |
---|
| 615 | // connected to the first (second) newly created vertex by a parent-child relation. |
---|
| 616 | //********************************************************************************************************* |
---|
[1172] | 617 | |
---|
[1186] | 618 | |
---|
| 619 | // Create the vectors of vectors of pointers to polyhedrons to hold the copies of the old Hasse diagram |
---|
[1172] | 620 | vector<vector<polyhedron*>> new_statistic1; |
---|
| 621 | vector<vector<polyhedron*>> new_statistic2; |
---|
| 622 | |
---|
[1186] | 623 | // Copy the statistic by rows |
---|
[1172] | 624 | for(int j=0;j<statistic.size();j++) |
---|
[1171] | 625 | { |
---|
[1186] | 626 | // an element counter |
---|
[1171] | 627 | int element_number = 0; |
---|
| 628 | |
---|
[1186] | 629 | vector<polyhedron*> supportnew_1; |
---|
| 630 | vector<polyhedron*> supportnew_2; |
---|
| 631 | |
---|
| 632 | new_statistic1.push_back(supportnew_1); |
---|
| 633 | new_statistic2.push_back(supportnew_2); |
---|
| 634 | |
---|
| 635 | // for each polyhedron in the given row |
---|
[1172] | 636 | for(vector<polyhedron*>::iterator horiz_ref = statistic[j].begin();horiz_ref<statistic[j].end();horiz_ref++) |
---|
[1171] | 637 | { |
---|
[1186] | 638 | // Append an extra zero coordinate to each of the vertices for the new dimension |
---|
| 639 | // If j==0 => we loop through vertices |
---|
| 640 | if(j == 0) |
---|
| 641 | { |
---|
| 642 | // cast the polyhedron pointer to a vertex pointer and push a zero to its vector of coordinates |
---|
| 643 | ((vertex*) (*horiz_ref))->push_coordinate(0); |
---|
| 644 | } |
---|
| 645 | |
---|
| 646 | // if it has parents |
---|
[1172] | 647 | if(!(*horiz_ref)->parents.empty()) |
---|
[1171] | 648 | { |
---|
[1186] | 649 | // save the relative address of this child in a vector kids_rel_addresses of all its parents. |
---|
| 650 | // This information will later be used for copying the whole Hasse diagram with each of the |
---|
| 651 | // relations contained within. |
---|
[1172] | 652 | for(vector<polyhedron*>::iterator parent_ref = (*horiz_ref)->parents.begin();parent_ref < (*horiz_ref)->parents.end();parent_ref++) |
---|
[1171] | 653 | { |
---|
[1172] | 654 | (*parent_ref)->kids_rel_addresses.push_back(element_number); |
---|
| 655 | } |
---|
[1171] | 656 | } |
---|
| 657 | |
---|
[1186] | 658 | // ************************************************************************************************** |
---|
| 659 | // Here we begin creating a new polyhedron, which will be a copy of the old one. Each such polyhedron |
---|
| 660 | // will be created as a toprow, but this information will be later forgotten and only the polyhedrons |
---|
| 661 | // in the top row of the Hasse diagram will be considered toprow for later use. |
---|
| 662 | // ************************************************************************************************** |
---|
| 663 | |
---|
| 664 | // First we create vectors specifying a toprow condition. In the case of a preconstructed statistic |
---|
| 665 | // this condition will be a vector of zeros. There are two vectors, because we need two copies of |
---|
| 666 | // the original Hasse diagram. |
---|
[1204] | 667 | vec vec1(i+2); |
---|
| 668 | vec1.zeros(); |
---|
[1171] | 669 | |
---|
[1204] | 670 | vec vec2(i+2); |
---|
| 671 | vec2.zeros(); |
---|
| 672 | |
---|
[1186] | 673 | // We create a new toprow with the previously specified condition. |
---|
[1172] | 674 | toprow *current_copy1 = new toprow(vec1); |
---|
| 675 | toprow *current_copy2 = new toprow(vec2); |
---|
[1171] | 676 | |
---|
[1186] | 677 | // The vertices of the copies will be inherited, because there will be a parent/child relation |
---|
| 678 | // between each polyhedron and its offspring (comming from the copy) and a parent has all the |
---|
| 679 | // vertices of its child plus more. |
---|
[1172] | 680 | for(vector<vertex*>::iterator vert_ref = (*horiz_ref)->vertices.begin();vert_ref<(*horiz_ref)->vertices.end();vert_ref++) |
---|
[1171] | 681 | { |
---|
[1172] | 682 | current_copy1->vertices.push_back(*vert_ref); |
---|
| 683 | current_copy2->vertices.push_back(*vert_ref); |
---|
[1171] | 684 | } |
---|
[1172] | 685 | |
---|
[1186] | 686 | // The only new vertex of the offspring should be the newly created point. |
---|
[1172] | 687 | current_copy1->vertices.push_back(new_point1); |
---|
| 688 | current_copy2->vertices.push_back(new_point2); |
---|
[1186] | 689 | |
---|
| 690 | // This method guarantees that each polyhedron is already triangulated, therefore its triangulation |
---|
| 691 | // is only one set of vertices and it is the set of all its vertices. |
---|
[1172] | 692 | current_copy1->triangulations.push_back(current_copy1->vertices); |
---|
| 693 | current_copy2->triangulations.push_back(current_copy2->vertices); |
---|
| 694 | |
---|
[1186] | 695 | // Now we have copied the polyhedron and we have to copy all of its relations. Because we are copying |
---|
| 696 | // in the Hasse diagram from bottom up, we always have to copy the parent/child relations to all the |
---|
| 697 | // kids and when we do that and know the child, in the child we will remember the parent we came from. |
---|
| 698 | // This way all the parents/children relations are saved in both the parent and the child. |
---|
[1172] | 699 | if(!(*horiz_ref)->kids_rel_addresses.empty()) |
---|
[1171] | 700 | { |
---|
[1172] | 701 | for(vector<int>::iterator kid_ref = (*horiz_ref)->kids_rel_addresses.begin();kid_ref<(*horiz_ref)->kids_rel_addresses.end();kid_ref++) |
---|
[1186] | 702 | { |
---|
| 703 | // find the child and save the relation to the parent |
---|
| 704 | current_copy1->children.push_back(new_statistic1[j-1][(*kid_ref)]); |
---|
| 705 | current_copy2->children.push_back(new_statistic2[j-1][(*kid_ref)]); |
---|
[1171] | 706 | |
---|
[1186] | 707 | // in the child save the parents' address |
---|
| 708 | new_statistic1[j-1][(*kid_ref)]->parents.push_back(current_copy1); |
---|
| 709 | new_statistic2[j-1][(*kid_ref)]->parents.push_back(current_copy2); |
---|
[1172] | 710 | } |
---|
[1171] | 711 | |
---|
[1186] | 712 | // Here we clear the parents kids_rel_addresses vector for later use (when we need to widen the |
---|
| 713 | // Hasse diagram again) |
---|
[1172] | 714 | (*horiz_ref)->kids_rel_addresses.clear(); |
---|
[1171] | 715 | } |
---|
[1186] | 716 | // If there were no children previously, we are copying a polyhedron that has been a vertex before. |
---|
| 717 | // In this case it is a segment now and it will have a relation to its mother (copywise) and to the |
---|
| 718 | // newly created point. Here we create the connection to the new point, again from both sides. |
---|
[1171] | 719 | else |
---|
| 720 | { |
---|
[1186] | 721 | // Add the address of the new point in the former vertex |
---|
[1172] | 722 | current_copy1->children.push_back(new_point1); |
---|
| 723 | current_copy2->children.push_back(new_point2); |
---|
[1171] | 724 | |
---|
[1186] | 725 | // Add the address of the former vertex in the new point |
---|
[1172] | 726 | new_point1->parents.push_back(current_copy1); |
---|
| 727 | new_point2->parents.push_back(current_copy2); |
---|
[1171] | 728 | } |
---|
| 729 | |
---|
[1186] | 730 | // Save the mother in its offspring |
---|
[1172] | 731 | current_copy1->children.push_back(*horiz_ref); |
---|
| 732 | current_copy2->children.push_back(*horiz_ref); |
---|
[1171] | 733 | |
---|
[1186] | 734 | // Save the offspring in its mother |
---|
| 735 | (*horiz_ref)->parents.push_back(current_copy1); |
---|
| 736 | (*horiz_ref)->parents.push_back(current_copy2); |
---|
| 737 | |
---|
[1171] | 738 | |
---|
[1186] | 739 | // Add the copies into the relevant statistic. The statistic will later be appended to the previous |
---|
| 740 | // Hasse diagram |
---|
| 741 | new_statistic1[j].push_back(current_copy1); |
---|
| 742 | new_statistic2[j].push_back(current_copy2); |
---|
| 743 | |
---|
| 744 | // Raise the count in the vector of polyhedrons |
---|
[1171] | 745 | element_number++; |
---|
[1186] | 746 | |
---|
[1172] | 747 | } |
---|
[1186] | 748 | } |
---|
| 749 | |
---|
| 750 | statistic[0].push_back(new_point1); |
---|
| 751 | statistic[0].push_back(new_point2); |
---|
| 752 | |
---|
| 753 | // Merge the new statistics into the old one. This will either be the final statistic or we will |
---|
| 754 | // reenter the widening loop. |
---|
| 755 | for(int j=0;j<new_statistic1.size();j++) |
---|
| 756 | { |
---|
| 757 | if(j+1==statistic.size()) |
---|
| 758 | { |
---|
| 759 | vector<polyhedron*> support; |
---|
| 760 | statistic.push_back(support); |
---|
| 761 | } |
---|
| 762 | |
---|
| 763 | statistic[j+1].insert(statistic[j+1].end(),new_statistic1[j].begin(),new_statistic1[j].end()); |
---|
| 764 | statistic[j+1].insert(statistic[j+1].end(),new_statistic2[j].begin(),new_statistic2[j].end()); |
---|
| 765 | } |
---|
[1171] | 766 | } |
---|
| 767 | } |
---|
| 768 | |
---|
| 769 | |
---|
[976] | 770 | |
---|
[1171] | 771 | |
---|
[976] | 772 | }; |
---|
| 773 | |
---|
[1208] | 774 | /* |
---|
| 775 | |
---|
[976] | 776 | //! Robust Bayesian AR model for Multicriteria-Laplace-Inverse-Gamma density |
---|
[1204] | 777 | class RARX : public BM |
---|
| 778 | { |
---|
| 779 | private: |
---|
| 780 | |
---|
| 781 | emlig posterior; |
---|
| 782 | |
---|
| 783 | public: |
---|
| 784 | RARX():BM() |
---|
| 785 | { |
---|
| 786 | }; |
---|
| 787 | |
---|
| 788 | void bayes(const itpp::vec &yt, const itpp::vec &cond = empty_vec) |
---|
| 789 | { |
---|
| 790 | |
---|
| 791 | } |
---|
| 792 | |
---|
[1208] | 793 | };*/ |
---|
[976] | 794 | |
---|
| 795 | |
---|
[1208] | 796 | |
---|
[976] | 797 | #endif //TRAGE_H |
---|