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> |
---|
11 | #include <itpp/itbase.h> |
---|
12 | #include <map> |
---|
13 | #include <limits> |
---|
14 | #include <vector> |
---|
15 | #include <list> |
---|
16 | #include <set> |
---|
17 | #include <algorithm> |
---|
18 | |
---|
19 | //using namespace bdm; |
---|
20 | using namespace std; |
---|
21 | using namespace itpp; |
---|
22 | |
---|
23 | const double max_range = 10;//numeric_limits<double>::max()/10e-10; |
---|
24 | |
---|
25 | enum actions {MERGE, SPLIT}; |
---|
26 | |
---|
27 | |
---|
28 | |
---|
29 | class polyhedron; |
---|
30 | class vertex; |
---|
31 | |
---|
32 | /* |
---|
33 | class t_simplex |
---|
34 | { |
---|
35 | public: |
---|
36 | set<vertex*> minima; |
---|
37 | |
---|
38 | set<vertex*> simplex; |
---|
39 | |
---|
40 | t_simplex(vertex* origin_vertex) |
---|
41 | { |
---|
42 | simplex.insert(origin_vertex); |
---|
43 | minima.insert(origin_vertex); |
---|
44 | } |
---|
45 | };*/ |
---|
46 | |
---|
47 | class emlig; |
---|
48 | class polyhedron; |
---|
49 | |
---|
50 | class condition |
---|
51 | { |
---|
52 | public: |
---|
53 | vec value; |
---|
54 | |
---|
55 | int multiplicity; |
---|
56 | |
---|
57 | condition(vec value) |
---|
58 | { |
---|
59 | this->value = value; |
---|
60 | multiplicity = 1; |
---|
61 | } |
---|
62 | }; |
---|
63 | |
---|
64 | |
---|
65 | /// A class describing a single polyhedron of the split complex. From a collection of such classes a Hasse diagram |
---|
66 | /// of the structure in the exponent of a Laplace-Inverse-Gamma density will be created. |
---|
67 | class polyhedron |
---|
68 | { |
---|
69 | /// A property having a value of 1 usually, with higher value only if the polyhedron arises as a coincidence of |
---|
70 | /// more than just the necessary number of conditions. For example if a newly created line passes through an already |
---|
71 | /// existing point, the points multiplicity will rise by 1. |
---|
72 | int multiplicity; |
---|
73 | |
---|
74 | int split_state; |
---|
75 | |
---|
76 | int merge_state; |
---|
77 | |
---|
78 | |
---|
79 | |
---|
80 | public: |
---|
81 | emlig* my_emlig; |
---|
82 | |
---|
83 | /// A list of polyhedrons parents within the Hasse diagram. |
---|
84 | list<polyhedron*> parents; |
---|
85 | |
---|
86 | /// A list of polyhedrons children withing the Hasse diagram. |
---|
87 | list<polyhedron*> children; |
---|
88 | |
---|
89 | /// All the vertices of the given polyhedron |
---|
90 | set<vertex*> vertices; |
---|
91 | |
---|
92 | set<condition*> parentconditions; |
---|
93 | |
---|
94 | /// A list used for storing children that lie in the positive region related to a certain condition |
---|
95 | list<polyhedron*> positivechildren; |
---|
96 | |
---|
97 | /// A list used for storing children that lie in the negative region related to a certain condition |
---|
98 | list<polyhedron*> negativechildren; |
---|
99 | |
---|
100 | /// Children intersecting the condition |
---|
101 | list<polyhedron*> neutralchildren; |
---|
102 | |
---|
103 | set<polyhedron*> totallyneutralgrandchildren; |
---|
104 | |
---|
105 | set<polyhedron*> totallyneutralchildren; |
---|
106 | |
---|
107 | set<polyhedron*> grandparents; |
---|
108 | |
---|
109 | set<vertex*> positiveneutralvertices; |
---|
110 | |
---|
111 | set<vertex*> negativeneutralvertices; |
---|
112 | |
---|
113 | bool totally_neutral; |
---|
114 | |
---|
115 | polyhedron* mergechild; |
---|
116 | |
---|
117 | polyhedron* positiveparent; |
---|
118 | |
---|
119 | polyhedron* negativeparent; |
---|
120 | |
---|
121 | polyhedron* next_poly; |
---|
122 | |
---|
123 | polyhedron* prev_poly; |
---|
124 | |
---|
125 | int message_counter; |
---|
126 | |
---|
127 | /// List of triangulation polyhedrons of the polyhedron given by their relative vertices. |
---|
128 | list<set<vertex*>> triangulation; |
---|
129 | |
---|
130 | /// A list of relative addresses serving for Hasse diagram construction. |
---|
131 | list<int> kids_rel_addresses; |
---|
132 | |
---|
133 | /// Default constructor |
---|
134 | polyhedron() |
---|
135 | { |
---|
136 | multiplicity = 1; |
---|
137 | |
---|
138 | message_counter = 0; |
---|
139 | |
---|
140 | totally_neutral = NULL; |
---|
141 | |
---|
142 | mergechild = NULL; |
---|
143 | } |
---|
144 | |
---|
145 | /// Setter for raising multiplicity |
---|
146 | void raise_multiplicity() |
---|
147 | { |
---|
148 | multiplicity++; |
---|
149 | } |
---|
150 | |
---|
151 | /// Setter for lowering multiplicity |
---|
152 | void lower_multiplicity() |
---|
153 | { |
---|
154 | multiplicity--; |
---|
155 | } |
---|
156 | |
---|
157 | int get_multiplicity() |
---|
158 | { |
---|
159 | return multiplicity; |
---|
160 | } |
---|
161 | |
---|
162 | /// An obligatory operator, when the class is used within a C++ STL structure like a vector |
---|
163 | int operator==(polyhedron polyhedron2) |
---|
164 | { |
---|
165 | return true; |
---|
166 | } |
---|
167 | |
---|
168 | /// An obligatory operator, when the class is used within a C++ STL structure like a vector |
---|
169 | int operator<(polyhedron polyhedron2) |
---|
170 | { |
---|
171 | return false; |
---|
172 | } |
---|
173 | |
---|
174 | |
---|
175 | |
---|
176 | int set_state(double state_indicator, actions action) |
---|
177 | { |
---|
178 | switch(action) |
---|
179 | { |
---|
180 | case MERGE: |
---|
181 | merge_state = (int)sign(state_indicator); |
---|
182 | return merge_state; |
---|
183 | case SPLIT: |
---|
184 | split_state = (int)sign(state_indicator); |
---|
185 | return split_state; |
---|
186 | } |
---|
187 | } |
---|
188 | |
---|
189 | int get_state(actions action) |
---|
190 | { |
---|
191 | switch(action) |
---|
192 | { |
---|
193 | case MERGE: |
---|
194 | return merge_state; |
---|
195 | break; |
---|
196 | case SPLIT: |
---|
197 | return split_state; |
---|
198 | break; |
---|
199 | } |
---|
200 | } |
---|
201 | |
---|
202 | int number_of_children() |
---|
203 | { |
---|
204 | return children.size(); |
---|
205 | } |
---|
206 | |
---|
207 | |
---|
208 | void triangulate(bool should_integrate); |
---|
209 | }; |
---|
210 | |
---|
211 | |
---|
212 | /// A class for representing 0-dimensional polyhedron - a vertex. It will be located in the bottom row of the Hasse |
---|
213 | /// diagram representing a complex of polyhedrons. It has its coordinates in the parameter space. |
---|
214 | class vertex : public polyhedron |
---|
215 | { |
---|
216 | /// A dynamic array representing coordinates of the vertex |
---|
217 | vec coordinates; |
---|
218 | |
---|
219 | public: |
---|
220 | |
---|
221 | double function_value; |
---|
222 | |
---|
223 | /// Default constructor |
---|
224 | vertex(); |
---|
225 | |
---|
226 | /// Constructor of a vertex from a set of coordinates |
---|
227 | vertex(vec coordinates) |
---|
228 | { |
---|
229 | this->coordinates = coordinates; |
---|
230 | |
---|
231 | vertices.insert(this); |
---|
232 | |
---|
233 | set<vertex*> vert_simplex; |
---|
234 | |
---|
235 | vert_simplex.insert(this); |
---|
236 | |
---|
237 | triangulation.push_back(vert_simplex); |
---|
238 | } |
---|
239 | |
---|
240 | /// A method that widens the set of coordinates of given vertex. It is used when a complex in a parameter |
---|
241 | /// space of certain dimension is established, but the dimension is not known when the vertex is created. |
---|
242 | void push_coordinate(double coordinate) |
---|
243 | { |
---|
244 | coordinates = concat(coordinates,coordinate); |
---|
245 | } |
---|
246 | |
---|
247 | /// A method obtaining the set of coordinates of a vertex. These coordinates are not obtained as a pointer |
---|
248 | /// (not given by reference), but a new copy is created (they are given by value). |
---|
249 | vec get_coordinates() |
---|
250 | { |
---|
251 | return coordinates; |
---|
252 | } |
---|
253 | |
---|
254 | }; |
---|
255 | |
---|
256 | |
---|
257 | /// A class representing a polyhedron in a top row of the complex. Such polyhedron has a condition that differitiates |
---|
258 | /// it from polyhedrons in other rows. |
---|
259 | class toprow : public polyhedron |
---|
260 | { |
---|
261 | |
---|
262 | public: |
---|
263 | double probability; |
---|
264 | |
---|
265 | vertex* minimal_vertex; |
---|
266 | |
---|
267 | /// A condition used for determining the function of a Laplace-Inverse-Gamma density resulting from Bayesian estimation |
---|
268 | vec condition_sum; |
---|
269 | |
---|
270 | int condition_order; |
---|
271 | |
---|
272 | /// Default constructor |
---|
273 | toprow(){}; |
---|
274 | |
---|
275 | /// Constructor creating a toprow from the condition |
---|
276 | toprow(condition *condition, int condition_order) |
---|
277 | { |
---|
278 | this->condition_sum = condition->value; |
---|
279 | this->condition_order = condition_order; |
---|
280 | } |
---|
281 | |
---|
282 | toprow(vec condition_sum, int condition_order) |
---|
283 | { |
---|
284 | this->condition_sum = condition_sum; |
---|
285 | } |
---|
286 | |
---|
287 | double integrate_simplex(set<vertex*> simplex, char c); |
---|
288 | |
---|
289 | }; |
---|
290 | |
---|
291 | |
---|
292 | |
---|
293 | |
---|
294 | |
---|
295 | |
---|
296 | class c_statistic |
---|
297 | { |
---|
298 | |
---|
299 | public: |
---|
300 | polyhedron* end_poly; |
---|
301 | polyhedron* start_poly; |
---|
302 | |
---|
303 | vector<polyhedron*> rows; |
---|
304 | |
---|
305 | vector<polyhedron*> row_ends; |
---|
306 | |
---|
307 | c_statistic() |
---|
308 | { |
---|
309 | end_poly = new polyhedron(); |
---|
310 | start_poly = new polyhedron(); |
---|
311 | }; |
---|
312 | |
---|
313 | void append_polyhedron(int row, polyhedron* appended_start, polyhedron* appended_end) |
---|
314 | { |
---|
315 | if(row>((int)rows.size())-1) |
---|
316 | { |
---|
317 | if(row>rows.size()) |
---|
318 | { |
---|
319 | throw new exception("You are trying to append a polyhedron whose children are not in the statistic yet!"); |
---|
320 | return; |
---|
321 | } |
---|
322 | |
---|
323 | rows.push_back(end_poly); |
---|
324 | row_ends.push_back(end_poly); |
---|
325 | } |
---|
326 | |
---|
327 | // POSSIBLE FAILURE: the function is not checking if start and end are connected |
---|
328 | |
---|
329 | if(rows[row] != end_poly) |
---|
330 | { |
---|
331 | appended_start->prev_poly = row_ends[row]; |
---|
332 | row_ends[row]->next_poly = appended_start; |
---|
333 | |
---|
334 | } |
---|
335 | else if((row>0 && rows[row-1]!=end_poly)||row==0) |
---|
336 | { |
---|
337 | appended_start->prev_poly = start_poly; |
---|
338 | rows[row]= appended_start; |
---|
339 | } |
---|
340 | else |
---|
341 | { |
---|
342 | throw new exception("Wrong polyhedron insertion into statistic: missing intermediary polyhedron!"); |
---|
343 | } |
---|
344 | |
---|
345 | appended_end->next_poly = end_poly; |
---|
346 | row_ends[row] = appended_end; |
---|
347 | } |
---|
348 | |
---|
349 | void append_polyhedron(int row, polyhedron* appended_poly) |
---|
350 | { |
---|
351 | append_polyhedron(row,appended_poly,appended_poly); |
---|
352 | } |
---|
353 | |
---|
354 | void insert_polyhedron(int row, polyhedron* inserted_poly, polyhedron* following_poly) |
---|
355 | { |
---|
356 | if(following_poly != end_poly) |
---|
357 | { |
---|
358 | inserted_poly->next_poly = following_poly; |
---|
359 | inserted_poly->prev_poly = following_poly->prev_poly; |
---|
360 | |
---|
361 | if(following_poly->prev_poly == start_poly) |
---|
362 | { |
---|
363 | rows[row] = inserted_poly; |
---|
364 | } |
---|
365 | else |
---|
366 | { |
---|
367 | inserted_poly->prev_poly->next_poly = inserted_poly; |
---|
368 | } |
---|
369 | |
---|
370 | following_poly->prev_poly = inserted_poly; |
---|
371 | } |
---|
372 | else |
---|
373 | { |
---|
374 | this->append_polyhedron(row, inserted_poly); |
---|
375 | } |
---|
376 | |
---|
377 | } |
---|
378 | |
---|
379 | |
---|
380 | |
---|
381 | |
---|
382 | void delete_polyhedron(int row, polyhedron* deleted_poly) |
---|
383 | { |
---|
384 | if(deleted_poly->prev_poly != start_poly) |
---|
385 | { |
---|
386 | deleted_poly->prev_poly->next_poly = deleted_poly->next_poly; |
---|
387 | } |
---|
388 | else |
---|
389 | { |
---|
390 | rows[row] = deleted_poly->next_poly; |
---|
391 | } |
---|
392 | |
---|
393 | if(deleted_poly->next_poly!=end_poly) |
---|
394 | { |
---|
395 | deleted_poly->next_poly->prev_poly = deleted_poly->prev_poly; |
---|
396 | } |
---|
397 | else |
---|
398 | { |
---|
399 | row_ends[row] = deleted_poly->prev_poly; |
---|
400 | } |
---|
401 | |
---|
402 | |
---|
403 | |
---|
404 | deleted_poly->next_poly = NULL; |
---|
405 | deleted_poly->prev_poly = NULL; |
---|
406 | } |
---|
407 | |
---|
408 | int size() |
---|
409 | { |
---|
410 | return rows.size(); |
---|
411 | } |
---|
412 | |
---|
413 | polyhedron* get_end() |
---|
414 | { |
---|
415 | return end_poly; |
---|
416 | } |
---|
417 | |
---|
418 | polyhedron* get_start() |
---|
419 | { |
---|
420 | return start_poly; |
---|
421 | } |
---|
422 | |
---|
423 | int row_size(int row) |
---|
424 | { |
---|
425 | if(this->size()>row && row>=0) |
---|
426 | { |
---|
427 | int row_size = 0; |
---|
428 | |
---|
429 | for(polyhedron* row_poly = rows[row]; row_poly!=end_poly; row_poly=row_poly->next_poly) |
---|
430 | { |
---|
431 | row_size++; |
---|
432 | } |
---|
433 | |
---|
434 | return row_size; |
---|
435 | } |
---|
436 | else |
---|
437 | { |
---|
438 | throw new exception("There is no row to obtain size from!"); |
---|
439 | } |
---|
440 | } |
---|
441 | }; |
---|
442 | |
---|
443 | |
---|
444 | class my_ivec : public ivec |
---|
445 | { |
---|
446 | public: |
---|
447 | my_ivec():ivec(){}; |
---|
448 | |
---|
449 | my_ivec(ivec origin):ivec() |
---|
450 | { |
---|
451 | this->ins(0,origin); |
---|
452 | } |
---|
453 | |
---|
454 | bool operator>(const my_ivec &second) const |
---|
455 | { |
---|
456 | return max(*this)>max(second); |
---|
457 | |
---|
458 | /* |
---|
459 | int size1 = this->size(); |
---|
460 | int size2 = second.size(); |
---|
461 | |
---|
462 | int counter1 = 0; |
---|
463 | while(0==0) |
---|
464 | { |
---|
465 | if((*this)[counter1]==0) |
---|
466 | { |
---|
467 | size1--; |
---|
468 | } |
---|
469 | |
---|
470 | if((*this)[counter1]!=0) |
---|
471 | break; |
---|
472 | |
---|
473 | counter1++; |
---|
474 | } |
---|
475 | |
---|
476 | int counter2 = 0; |
---|
477 | while(0==0) |
---|
478 | { |
---|
479 | if(second[counter2]==0) |
---|
480 | { |
---|
481 | size2--; |
---|
482 | } |
---|
483 | |
---|
484 | if(second[counter2]!=0) |
---|
485 | break; |
---|
486 | |
---|
487 | counter2++; |
---|
488 | } |
---|
489 | |
---|
490 | if(size1!=size2) |
---|
491 | { |
---|
492 | return size1>size2; |
---|
493 | } |
---|
494 | else |
---|
495 | { |
---|
496 | for(int i = 0;i<size1;i++) |
---|
497 | { |
---|
498 | if((*this)[counter1+i]!=second[counter2+i]) |
---|
499 | { |
---|
500 | return (*this)[counter1+i]>second[counter2+i]; |
---|
501 | } |
---|
502 | } |
---|
503 | |
---|
504 | return false; |
---|
505 | }*/ |
---|
506 | } |
---|
507 | |
---|
508 | |
---|
509 | bool operator==(const my_ivec &second) const |
---|
510 | { |
---|
511 | return max(*this)==max(second); |
---|
512 | |
---|
513 | /* |
---|
514 | int size1 = this->size(); |
---|
515 | int size2 = second.size(); |
---|
516 | |
---|
517 | int counter = 0; |
---|
518 | while(0==0) |
---|
519 | { |
---|
520 | if((*this)[counter]==0) |
---|
521 | { |
---|
522 | size1--; |
---|
523 | } |
---|
524 | |
---|
525 | if((*this)[counter]!=0) |
---|
526 | break; |
---|
527 | |
---|
528 | counter++; |
---|
529 | } |
---|
530 | |
---|
531 | counter = 0; |
---|
532 | while(0==0) |
---|
533 | { |
---|
534 | if(second[counter]==0) |
---|
535 | { |
---|
536 | size2--; |
---|
537 | } |
---|
538 | |
---|
539 | if(second[counter]!=0) |
---|
540 | break; |
---|
541 | |
---|
542 | counter++; |
---|
543 | } |
---|
544 | |
---|
545 | if(size1!=size2) |
---|
546 | { |
---|
547 | return false; |
---|
548 | } |
---|
549 | else |
---|
550 | { |
---|
551 | for(int i=0;i<size1;i++) |
---|
552 | { |
---|
553 | if((*this)[size()-1-i]!=second[second.size()-1-i]) |
---|
554 | { |
---|
555 | return false; |
---|
556 | } |
---|
557 | } |
---|
558 | |
---|
559 | return true; |
---|
560 | }*/ |
---|
561 | } |
---|
562 | |
---|
563 | bool operator<(const my_ivec &second) const |
---|
564 | { |
---|
565 | return !(((*this)>second)||((*this)==second)); |
---|
566 | } |
---|
567 | |
---|
568 | bool operator!=(const my_ivec &second) const |
---|
569 | { |
---|
570 | return !((*this)==second); |
---|
571 | } |
---|
572 | |
---|
573 | bool operator<=(const my_ivec &second) const |
---|
574 | { |
---|
575 | return !((*this)>second); |
---|
576 | } |
---|
577 | |
---|
578 | bool operator>=(const my_ivec &second) const |
---|
579 | { |
---|
580 | return !((*this)<second); |
---|
581 | } |
---|
582 | |
---|
583 | my_ivec right(my_ivec original) |
---|
584 | { |
---|
585 | |
---|
586 | } |
---|
587 | }; |
---|
588 | |
---|
589 | |
---|
590 | |
---|
591 | |
---|
592 | |
---|
593 | |
---|
594 | |
---|
595 | //! Conditional(e) Multicriteria-Laplace-Inverse-Gamma distribution density |
---|
596 | class emlig // : eEF |
---|
597 | { |
---|
598 | |
---|
599 | /// A statistic in a form of a Hasse diagram representing a complex of convex polyhedrons obtained as a result |
---|
600 | /// of data update from Bayesian estimation or set by the user if this emlig is a prior density |
---|
601 | |
---|
602 | |
---|
603 | vector<list<polyhedron*>> for_splitting; |
---|
604 | |
---|
605 | vector<list<polyhedron*>> for_merging; |
---|
606 | |
---|
607 | list<condition*> conditions; |
---|
608 | |
---|
609 | double normalization_factor; |
---|
610 | |
---|
611 | |
---|
612 | |
---|
613 | void alter_toprow_conditions(condition *condition, bool should_be_added) |
---|
614 | { |
---|
615 | for(polyhedron* horiz_ref = statistic.rows[statistic.size()-1];horiz_ref!=statistic.get_end();horiz_ref=horiz_ref->next_poly) |
---|
616 | { |
---|
617 | set<vertex*>::iterator vertex_ref = horiz_ref->vertices.begin(); |
---|
618 | |
---|
619 | do |
---|
620 | { |
---|
621 | vertex_ref++; |
---|
622 | } |
---|
623 | while((*vertex_ref)->parentconditions.find(condition)==(*vertex_ref)->parentconditions.end()); |
---|
624 | |
---|
625 | double product = (*vertex_ref)->get_coordinates()*condition->value; |
---|
626 | |
---|
627 | if((product>0 && should_be_added)||(product<0 && !should_be_added)) |
---|
628 | { |
---|
629 | ((toprow*) horiz_ref)->condition_sum += condition->value; |
---|
630 | } |
---|
631 | else |
---|
632 | { |
---|
633 | ((toprow*) horiz_ref)->condition_sum -= condition->value; |
---|
634 | } |
---|
635 | } |
---|
636 | } |
---|
637 | |
---|
638 | |
---|
639 | |
---|
640 | void send_state_message(polyhedron* sender, condition *toadd, condition *toremove, int level) |
---|
641 | { |
---|
642 | |
---|
643 | bool shouldmerge = (toremove == NULL); |
---|
644 | bool shouldsplit = (toadd == NULL); |
---|
645 | |
---|
646 | if(shouldsplit||shouldmerge) |
---|
647 | { |
---|
648 | for(list<polyhedron*>::iterator parent_iterator = sender->parents.begin();parent_iterator!=sender->parents.end();parent_iterator++) |
---|
649 | { |
---|
650 | polyhedron* current_parent = *parent_iterator; |
---|
651 | |
---|
652 | current_parent->message_counter++; |
---|
653 | |
---|
654 | bool is_last = (current_parent->message_counter == current_parent->number_of_children()); |
---|
655 | bool is_first = (current_parent->message_counter == 1); |
---|
656 | |
---|
657 | if(shouldmerge) |
---|
658 | { |
---|
659 | int child_state = sender->get_state(MERGE); |
---|
660 | int parent_state = current_parent->get_state(MERGE); |
---|
661 | |
---|
662 | if(parent_state == 0||is_first) |
---|
663 | { |
---|
664 | parent_state = current_parent->set_state(child_state, MERGE); |
---|
665 | } |
---|
666 | |
---|
667 | if(child_state == 0) |
---|
668 | { |
---|
669 | if(current_parent->mergechild == NULL) |
---|
670 | { |
---|
671 | current_parent->mergechild = sender; |
---|
672 | } |
---|
673 | } |
---|
674 | |
---|
675 | if(is_last) |
---|
676 | { |
---|
677 | if(current_parent->mergechild!=NULL) |
---|
678 | { |
---|
679 | if(current_parent->mergechild->get_multiplicity()==1) |
---|
680 | { |
---|
681 | if(parent_state > 0) |
---|
682 | { |
---|
683 | current_parent->mergechild->positiveparent = current_parent; |
---|
684 | } |
---|
685 | |
---|
686 | if(parent_state < 0) |
---|
687 | { |
---|
688 | current_parent->mergechild->negativeparent = current_parent; |
---|
689 | } |
---|
690 | } |
---|
691 | } |
---|
692 | |
---|
693 | |
---|
694 | if(parent_state == 0) |
---|
695 | { |
---|
696 | for_merging[level+1].push_back(current_parent); |
---|
697 | } |
---|
698 | else |
---|
699 | { |
---|
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) |
---|
728 | { |
---|
729 | current_parent->totally_neutral = sender->totally_neutral; |
---|
730 | } |
---|
731 | else |
---|
732 | { |
---|
733 | current_parent->totally_neutral = current_parent->totally_neutral && sender->totally_neutral; |
---|
734 | } |
---|
735 | |
---|
736 | if(sender->totally_neutral) |
---|
737 | { |
---|
738 | current_parent->totallyneutralchildren.insert(sender); |
---|
739 | } |
---|
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 | } |
---|
812 | |
---|
813 | } |
---|
814 | |
---|
815 | } |
---|
816 | } |
---|
817 | |
---|
818 | public: |
---|
819 | c_statistic statistic; |
---|
820 | |
---|
821 | vertex* minimal_vertex; |
---|
822 | |
---|
823 | double likelihood_value; |
---|
824 | |
---|
825 | vector<multiset<my_ivec>> correction_factors; |
---|
826 | |
---|
827 | int number_of_parameters; |
---|
828 | |
---|
829 | /// A default constructor creates an emlig with predefined statistic representing only the range of the given |
---|
830 | /// parametric space, where the number of parameters of the needed model is given as a parameter to the constructor. |
---|
831 | emlig(int number_of_parameters) |
---|
832 | { |
---|
833 | this->number_of_parameters = number_of_parameters; |
---|
834 | |
---|
835 | create_statistic(number_of_parameters); |
---|
836 | |
---|
837 | likelihood_value = numeric_limits<double>::max(); |
---|
838 | } |
---|
839 | |
---|
840 | /// A constructor for creating an emlig when the user wants to create the statistic by himself. The creation of a |
---|
841 | /// statistic is needed outside the constructor. Used for a user defined prior distribution on the parameters. |
---|
842 | emlig(c_statistic statistic) |
---|
843 | { |
---|
844 | this->statistic = statistic; |
---|
845 | |
---|
846 | likelihood_value = numeric_limits<double>::max(); |
---|
847 | } |
---|
848 | |
---|
849 | void step_me(int marker) |
---|
850 | { |
---|
851 | /* |
---|
852 | for(int i = 0;i<statistic.size();i++) |
---|
853 | { |
---|
854 | for(polyhedron* horiz_ref = statistic.rows[i];horiz_ref!=statistic.get_end();horiz_ref=horiz_ref->next_poly) |
---|
855 | { |
---|
856 | if(i==statistic.size()-1) |
---|
857 | { |
---|
858 | cout << ((toprow*)horiz_ref)->condition << " " << ((toprow*)horiz_ref)->probability << endl; |
---|
859 | cout << "Order:" << ((toprow*)horiz_ref)->condition_order << endl; |
---|
860 | } |
---|
861 | char* string = "Checkpoint"; |
---|
862 | } |
---|
863 | } |
---|
864 | */ |
---|
865 | |
---|
866 | /* |
---|
867 | list<vec> table_entries; |
---|
868 | for(polyhedron* horiz_ref = statistic.rows[statistic.size()-1];horiz_ref!=statistic.row_ends[statistic.size()-1];horiz_ref=horiz_ref->next_poly) |
---|
869 | { |
---|
870 | toprow *current_toprow = (toprow*)(horiz_ref); |
---|
871 | for(list<set<vertex*>>::iterator tri_ref = current_toprow->triangulation.begin();tri_ref!=current_toprow->triangulation.end();tri_ref++) |
---|
872 | { |
---|
873 | for(set<vertex*>::iterator vert_ref = (*tri_ref).begin();vert_ref!=(*tri_ref).end();vert_ref++) |
---|
874 | { |
---|
875 | vec table_entry = vec(); |
---|
876 | |
---|
877 | table_entry.ins(0,(*vert_ref)->get_coordinates()*current_toprow->condition.get(1,current_toprow->condition.size()-1)-current_toprow->condition.get(0,0)); |
---|
878 | |
---|
879 | table_entry.ins(0,(*vert_ref)->get_coordinates()); |
---|
880 | |
---|
881 | table_entries.push_back(table_entry); |
---|
882 | } |
---|
883 | } |
---|
884 | } |
---|
885 | |
---|
886 | unique(table_entries.begin(),table_entries.end()); |
---|
887 | |
---|
888 | |
---|
889 | |
---|
890 | for(list<vec>::iterator entry_ref = table_entries.begin();entry_ref!=table_entries.end();entry_ref++) |
---|
891 | { |
---|
892 | ofstream myfile; |
---|
893 | myfile.open("robust_data.txt", ios::out | ios::app); |
---|
894 | if (myfile.is_open()) |
---|
895 | { |
---|
896 | for(int i = 0;i<(*entry_ref).size();i++) |
---|
897 | { |
---|
898 | myfile << (*entry_ref)[i] << ";"; |
---|
899 | } |
---|
900 | myfile << endl; |
---|
901 | |
---|
902 | myfile.close(); |
---|
903 | } |
---|
904 | else |
---|
905 | { |
---|
906 | cout << "File problem." << endl; |
---|
907 | } |
---|
908 | } |
---|
909 | */ |
---|
910 | |
---|
911 | |
---|
912 | return; |
---|
913 | } |
---|
914 | |
---|
915 | int statistic_rowsize(int row) |
---|
916 | { |
---|
917 | return statistic.row_size(row); |
---|
918 | } |
---|
919 | |
---|
920 | void add_condition(vec toadd) |
---|
921 | { |
---|
922 | vec null_vector = ""; |
---|
923 | |
---|
924 | add_and_remove_condition(toadd, null_vector); |
---|
925 | } |
---|
926 | |
---|
927 | |
---|
928 | void remove_condition(vec toremove) |
---|
929 | { |
---|
930 | vec null_vector = ""; |
---|
931 | |
---|
932 | add_and_remove_condition(null_vector, toremove); |
---|
933 | |
---|
934 | } |
---|
935 | |
---|
936 | void add_and_remove_condition(vec toadd, vec toremove) |
---|
937 | { |
---|
938 | likelihood_value = numeric_limits<double>::max(); |
---|
939 | |
---|
940 | bool should_remove = (toremove.size() != 0); |
---|
941 | bool should_add = (toadd.size() != 0); |
---|
942 | |
---|
943 | for_splitting.clear(); |
---|
944 | for_merging.clear(); |
---|
945 | |
---|
946 | for(int i = 0;i<statistic.size();i++) |
---|
947 | { |
---|
948 | list<polyhedron*> empty_split; |
---|
949 | list<polyhedron*> empty_merge; |
---|
950 | |
---|
951 | for_splitting.push_back(empty_split); |
---|
952 | for_merging.push_back(empty_merge); |
---|
953 | } |
---|
954 | |
---|
955 | list<condition*>::iterator toremove_ref = conditions.end(); |
---|
956 | bool condition_should_be_added = false; |
---|
957 | |
---|
958 | for(list<condition*>::iterator ref = conditions.begin();ref!=conditions.end();ref++) |
---|
959 | { |
---|
960 | if(should_remove) |
---|
961 | { |
---|
962 | if((*ref)->value == toremove) |
---|
963 | { |
---|
964 | if((*ref)->multiplicity>1) |
---|
965 | { |
---|
966 | (*ref)->multiplicity--; |
---|
967 | |
---|
968 | alter_toprow_conditions(*ref,false); |
---|
969 | |
---|
970 | should_remove = false; |
---|
971 | } |
---|
972 | else |
---|
973 | { |
---|
974 | toremove_ref = ref; |
---|
975 | } |
---|
976 | } |
---|
977 | } |
---|
978 | |
---|
979 | if(should_add) |
---|
980 | { |
---|
981 | if((*ref)->value == toadd) |
---|
982 | { |
---|
983 | (*ref)->multiplicity++; |
---|
984 | |
---|
985 | alter_toprow_conditions(*ref,true); |
---|
986 | |
---|
987 | should_add = false; |
---|
988 | } |
---|
989 | else |
---|
990 | { |
---|
991 | condition_should_be_added = true; |
---|
992 | } |
---|
993 | } |
---|
994 | } |
---|
995 | |
---|
996 | condition* condition_to_remove = NULL; |
---|
997 | |
---|
998 | if(toremove_ref!=conditions.end()) |
---|
999 | { |
---|
1000 | conditions.erase(toremove_ref); |
---|
1001 | condition_to_remove = *toremove_ref; |
---|
1002 | } |
---|
1003 | |
---|
1004 | condition* condition_to_add = NULL; |
---|
1005 | |
---|
1006 | if(condition_should_be_added) |
---|
1007 | { |
---|
1008 | conditions.push_back(new condition(toadd)); |
---|
1009 | } |
---|
1010 | |
---|
1011 | |
---|
1012 | |
---|
1013 | for(polyhedron* horizontal_position = statistic.rows[0];horizontal_position!=statistic.get_end();horizontal_position=horizontal_position->next_poly) |
---|
1014 | { |
---|
1015 | vertex* current_vertex = (vertex*)horizontal_position; |
---|
1016 | |
---|
1017 | if(should_add||should_remove) |
---|
1018 | { |
---|
1019 | vec appended_coords = current_vertex->get_coordinates(); |
---|
1020 | appended_coords.ins(0,-1.0); |
---|
1021 | |
---|
1022 | if(should_add) |
---|
1023 | { |
---|
1024 | double local_condition = 0;// = toadd*(appended_coords.first/=appended_coords.second); |
---|
1025 | |
---|
1026 | local_condition = appended_coords*toadd; |
---|
1027 | |
---|
1028 | current_vertex->set_state(local_condition,SPLIT); |
---|
1029 | |
---|
1030 | /// \TODO There should be a rounding error tolerance used here to insure we are not having too many points because of rounding error. |
---|
1031 | if(local_condition == 0) |
---|
1032 | { |
---|
1033 | current_vertex->totally_neutral = true; |
---|
1034 | |
---|
1035 | current_vertex->raise_multiplicity(); |
---|
1036 | |
---|
1037 | current_vertex->negativeneutralvertices.insert(current_vertex); |
---|
1038 | current_vertex->positiveneutralvertices.insert(current_vertex); |
---|
1039 | } |
---|
1040 | } |
---|
1041 | |
---|
1042 | if(should_remove) |
---|
1043 | { |
---|
1044 | set<condition*>::iterator cond_ref; |
---|
1045 | |
---|
1046 | for(cond_ref = current_vertex->parentconditions.begin();cond_ref!=current_vertex->parentconditions.end();cond_ref++) |
---|
1047 | { |
---|
1048 | if(*cond_ref == condition_to_remove) |
---|
1049 | { |
---|
1050 | break; |
---|
1051 | } |
---|
1052 | } |
---|
1053 | |
---|
1054 | if(cond_ref!=current_vertex->parentconditions.end()) |
---|
1055 | { |
---|
1056 | current_vertex->parentconditions.erase(cond_ref); |
---|
1057 | current_vertex->set_state(0,MERGE); |
---|
1058 | for_merging[0].push_back(current_vertex); |
---|
1059 | } |
---|
1060 | else |
---|
1061 | { |
---|
1062 | double local_condition = toremove*appended_coords; |
---|
1063 | current_vertex->set_state(local_condition,MERGE); |
---|
1064 | } |
---|
1065 | } |
---|
1066 | } |
---|
1067 | |
---|
1068 | send_state_message(current_vertex, condition_to_add, condition_to_remove, 0); |
---|
1069 | |
---|
1070 | } |
---|
1071 | |
---|
1072 | |
---|
1073 | |
---|
1074 | if(should_remove) |
---|
1075 | { |
---|
1076 | for(int i = 0;i<for_merging.size();i++) |
---|
1077 | { |
---|
1078 | for(list<polyhedron*>::iterator merge_ref = for_merging[i].begin();merge_ref!=for_merging[i].end();merge_ref++) |
---|
1079 | { |
---|
1080 | cout << (*merge_ref)->get_state(MERGE) << ","; |
---|
1081 | } |
---|
1082 | |
---|
1083 | cout << endl; |
---|
1084 | } |
---|
1085 | |
---|
1086 | set<vertex*> vertices_to_be_reduced; |
---|
1087 | |
---|
1088 | int k = 1; |
---|
1089 | |
---|
1090 | for(vector<list<polyhedron*>>::iterator vert_ref = for_merging.begin();vert_ref<for_merging.end();vert_ref++) |
---|
1091 | { |
---|
1092 | for(list<polyhedron*>::reverse_iterator merge_ref = vert_ref->rbegin();merge_ref!=vert_ref->rend();merge_ref++) |
---|
1093 | { |
---|
1094 | if((*merge_ref)->get_multiplicity()>1) |
---|
1095 | { |
---|
1096 | if(k==1) |
---|
1097 | { |
---|
1098 | vertices_to_be_reduced.insert((vertex*)(*merge_ref)); |
---|
1099 | } |
---|
1100 | else |
---|
1101 | { |
---|
1102 | (*merge_ref)->lower_multiplicity(); |
---|
1103 | } |
---|
1104 | } |
---|
1105 | else |
---|
1106 | { |
---|
1107 | toprow* current_positive = (toprow*)(*merge_ref)->positiveparent; |
---|
1108 | toprow* current_negative = (toprow*)(*merge_ref)->negativeparent; |
---|
1109 | |
---|
1110 | current_positive->condition_sum -= toremove; |
---|
1111 | current_positive->condition_order--; |
---|
1112 | |
---|
1113 | current_positive->children.insert(current_positive->children.end(),current_negative->children.begin(),current_negative->children.end()); |
---|
1114 | current_positive->children.remove(*merge_ref); |
---|
1115 | |
---|
1116 | for(list<polyhedron*>::iterator child_ref = current_negative->children.begin();child_ref!=current_negative->children.end();child_ref++) |
---|
1117 | { |
---|
1118 | (*child_ref)->parents.remove(current_negative); |
---|
1119 | (*child_ref)->parents.push_back(current_positive); |
---|
1120 | } |
---|
1121 | |
---|
1122 | current_positive->negativechildren.insert(current_positive->negativechildren.end(),current_negative->negativechildren.begin(),current_negative->negativechildren.end()); |
---|
1123 | |
---|
1124 | current_positive->positivechildren.insert(current_positive->positivechildren.end(),current_negative->positivechildren.begin(),current_negative->positivechildren.end()); |
---|
1125 | |
---|
1126 | current_positive->neutralchildren.insert(current_positive->neutralchildren.end(),current_negative->positivechildren.begin(),current_negative->positivechildren.end()); |
---|
1127 | |
---|
1128 | switch((*merge_ref)->get_state(SPLIT)) |
---|
1129 | { |
---|
1130 | case -1: |
---|
1131 | current_positive->negativechildren.remove(*merge_ref); |
---|
1132 | break; |
---|
1133 | case 0: |
---|
1134 | current_positive->neutralchildren.remove(*merge_ref); |
---|
1135 | break; |
---|
1136 | case 1: |
---|
1137 | current_positive->positivechildren.remove(*merge_ref); |
---|
1138 | break; |
---|
1139 | } |
---|
1140 | |
---|
1141 | 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()); |
---|
1143 | |
---|
1144 | for(list<polyhedron*>::iterator parent_ref = current_negative->parents.begin();parent_ref!=current_negative->parents.end();parent_ref++) |
---|
1145 | { |
---|
1146 | (*parent_ref)->children.remove(current_negative); |
---|
1147 | (*parent_ref)->children.push_back(current_positive); |
---|
1148 | } |
---|
1149 | |
---|
1150 | current_positive->totallyneutralchildren.insert(current_negative->totallyneutralchildren.begin(),current_negative->totallyneutralchildren.end()); |
---|
1151 | current_positive->totallyneutralchildren.erase(*merge_ref); |
---|
1152 | |
---|
1153 | current_positive->totallyneutralgrandchildren.insert(current_negative->totallyneutralgrandchildren.begin(),current_negative->totallyneutralgrandchildren.end()); |
---|
1154 | |
---|
1155 | current_positive->vertices.insert(current_negative->vertices.begin(),current_negative->vertices.end()); |
---|
1156 | current_positive->negativeneutralvertices.insert(current_negative->negativeneutralvertices.begin(),current_negative->negativeneutralvertices.end()); |
---|
1157 | current_positive->positiveneutralvertices.insert(current_negative->positiveneutralvertices.begin(),current_negative->positiveneutralvertices.end()); |
---|
1158 | |
---|
1159 | for(set<vertex*>::iterator vert_ref = (*merge_ref)->vertices.begin();vert_ref!=(*merge_ref)->vertices.end();vert_ref++) |
---|
1160 | { |
---|
1161 | if((*vert_ref)->get_multiplicity()==1) |
---|
1162 | { |
---|
1163 | current_positive->vertices.erase(*vert_ref); |
---|
1164 | current_positive->negativeneutralvertices.erase(*vert_ref); |
---|
1165 | current_positive->positiveneutralvertices.erase(*vert_ref); |
---|
1166 | } |
---|
1167 | } |
---|
1168 | |
---|
1169 | if(current_negative->get_state(SPLIT)==0&&!current_negative->totally_neutral) |
---|
1170 | { |
---|
1171 | for_splitting[k].remove(current_negative); |
---|
1172 | |
---|
1173 | if(current_positive->get_state(SPLIT)!=0||current_positive->totally_neutral) |
---|
1174 | { |
---|
1175 | for_splitting[k].push_back(current_positive); |
---|
1176 | } |
---|
1177 | } |
---|
1178 | |
---|
1179 | if(current_positive->totally_neutral) |
---|
1180 | { |
---|
1181 | if(!current_negative->totally_neutral) |
---|
1182 | { |
---|
1183 | for(set<polyhedron*>::iterator grand_ref = current_positive->grandparents.begin();grand_ref!=current_positive->grandparents.end();grand_ref++) |
---|
1184 | { |
---|
1185 | (*grand_ref)->totallyneutralgrandchildren.erase(current_positive); |
---|
1186 | } |
---|
1187 | |
---|
1188 | current_positive->grandparents.clear(); |
---|
1189 | } |
---|
1190 | else |
---|
1191 | { |
---|
1192 | for(set<polyhedron*>::iterator grand_ref = current_negative->grandparents.begin();grand_ref!=current_negative->grandparents.end();grand_ref++) |
---|
1193 | { |
---|
1194 | (*grand_ref)->totallyneutralgrandchildren.erase(current_negative); |
---|
1195 | (*grand_ref)->totallyneutralgrandchildren.insert(current_positive); |
---|
1196 | } |
---|
1197 | |
---|
1198 | current_positive->grandparents.insert(current_negative->grandparents.begin(),current_negative->grandparents.end()); |
---|
1199 | } |
---|
1200 | } |
---|
1201 | else |
---|
1202 | { |
---|
1203 | if(current_negative->totally_neutral) |
---|
1204 | { |
---|
1205 | for(set<polyhedron*>::iterator grand_ref = current_negative->grandparents.begin();grand_ref!=current_negative->grandparents.end();grand_ref++) |
---|
1206 | { |
---|
1207 | (*grand_ref)->totallyneutralgrandchildren.erase(current_negative); |
---|
1208 | } |
---|
1209 | } |
---|
1210 | } |
---|
1211 | |
---|
1212 | current_positive->totally_neutral = (current_positive->totally_neutral && current_negative->totally_neutral); |
---|
1213 | |
---|
1214 | current_positive->triangulate(k==for_splitting.size()-1); |
---|
1215 | |
---|
1216 | statistic.delete_polyhedron(k,current_negative); |
---|
1217 | |
---|
1218 | delete current_negative; |
---|
1219 | |
---|
1220 | for(list<polyhedron*>::iterator child_ref = (*merge_ref)->children.begin();child_ref!=(*merge_ref)->children.end();child_ref++) |
---|
1221 | { |
---|
1222 | (*child_ref)->parents.remove(*merge_ref); |
---|
1223 | } |
---|
1224 | |
---|
1225 | for(list<polyhedron*>::iterator parent_ref = (*merge_ref)->parents.begin();parent_ref!=(*merge_ref)->parents.end();parent_ref++) |
---|
1226 | { |
---|
1227 | (*parent_ref)->positivechildren.remove(*merge_ref); |
---|
1228 | (*parent_ref)->negativechildren.remove(*merge_ref); |
---|
1229 | (*parent_ref)->neutralchildren.remove(*merge_ref); |
---|
1230 | (*parent_ref)->children.remove(*merge_ref); |
---|
1231 | } |
---|
1232 | |
---|
1233 | for(set<polyhedron*>::iterator grand_ch_ref = (*merge_ref)->totallyneutralgrandchildren.begin();grand_ch_ref!=(*merge_ref)->totallyneutralgrandchildren.end();grand_ch_ref++) |
---|
1234 | { |
---|
1235 | (*grand_ch_ref)->grandparents.erase(*merge_ref); |
---|
1236 | } |
---|
1237 | |
---|
1238 | for(set<polyhedron*>::iterator grand_p_ref = (*merge_ref)->grandparents.begin();grand_p_ref!=(*merge_ref)->grandparents.end();grand_p_ref++) |
---|
1239 | { |
---|
1240 | (*grand_p_ref)->totallyneutralgrandchildren.erase(*merge_ref); |
---|
1241 | } |
---|
1242 | |
---|
1243 | statistic.delete_polyhedron(k-1,*merge_ref); |
---|
1244 | |
---|
1245 | if(k==1) |
---|
1246 | { |
---|
1247 | vertices_to_be_reduced.insert((vertex*)(*merge_ref)); |
---|
1248 | } |
---|
1249 | else |
---|
1250 | { |
---|
1251 | delete *merge_ref; |
---|
1252 | } |
---|
1253 | } |
---|
1254 | } |
---|
1255 | |
---|
1256 | k++; |
---|
1257 | |
---|
1258 | } |
---|
1259 | |
---|
1260 | for(set<vertex*>::iterator vert_ref = vertices_to_be_reduced.begin();vert_ref!=vertices_to_be_reduced.end();vert_ref++) |
---|
1261 | { |
---|
1262 | if((*vert_ref)->get_multiplicity()>1) |
---|
1263 | { |
---|
1264 | (*vert_ref)->lower_multiplicity(); |
---|
1265 | } |
---|
1266 | else |
---|
1267 | { |
---|
1268 | delete *vert_ref; |
---|
1269 | } |
---|
1270 | } |
---|
1271 | } |
---|
1272 | |
---|
1273 | |
---|
1274 | if(should_add) |
---|
1275 | { |
---|
1276 | int k = 1; |
---|
1277 | |
---|
1278 | vector<list<polyhedron*>>::iterator beginning_ref = ++for_splitting.begin(); |
---|
1279 | |
---|
1280 | for(vector<list<polyhedron*>>::iterator vert_ref = beginning_ref;vert_ref<for_splitting.end();vert_ref++) |
---|
1281 | { |
---|
1282 | |
---|
1283 | for(list<polyhedron*>::reverse_iterator split_ref = vert_ref->rbegin();split_ref != vert_ref->rend();split_ref++) |
---|
1284 | { |
---|
1285 | polyhedron* new_totally_neutral_child; |
---|
1286 | |
---|
1287 | polyhedron* current_polyhedron = (*split_ref); |
---|
1288 | |
---|
1289 | if(vert_ref == beginning_ref) |
---|
1290 | { |
---|
1291 | vec coordinates1 = ((vertex*)(*(current_polyhedron->children.begin())))->get_coordinates(); |
---|
1292 | vec coordinates2 = ((vertex*)(*(++current_polyhedron->children.begin())))->get_coordinates(); |
---|
1293 | |
---|
1294 | vec extended_coord2 = coordinates2; |
---|
1295 | extended_coord2.ins(0,-1.0); |
---|
1296 | |
---|
1297 | double t = (-toadd*extended_coord2)/(toadd(1,toadd.size()-1)*(coordinates1-coordinates2)); |
---|
1298 | |
---|
1299 | vec new_coordinates = (1-t)*coordinates2+t*coordinates1; |
---|
1300 | |
---|
1301 | // cout << "c1:" << coordinates1 << endl << "c2:" << coordinates2 << endl << "nc:" << new_coordinates << endl; |
---|
1302 | |
---|
1303 | vertex* neutral_vertex = new vertex(new_coordinates); |
---|
1304 | |
---|
1305 | new_totally_neutral_child = neutral_vertex; |
---|
1306 | } |
---|
1307 | else |
---|
1308 | { |
---|
1309 | toprow* neutral_toprow = new toprow(); |
---|
1310 | |
---|
1311 | neutral_toprow->condition_sum = ((toprow*)current_polyhedron)->condition_sum; // tohle tu bylo driv: zeros(number_of_parameters+1); |
---|
1312 | neutral_toprow->condition_order = ((toprow*)current_polyhedron)->condition_order+1; |
---|
1313 | |
---|
1314 | new_totally_neutral_child = neutral_toprow; |
---|
1315 | } |
---|
1316 | |
---|
1317 | new_totally_neutral_child->my_emlig = this; |
---|
1318 | |
---|
1319 | new_totally_neutral_child->children.insert(new_totally_neutral_child->children.end(), |
---|
1320 | current_polyhedron->totallyneutralgrandchildren.begin(), |
---|
1321 | current_polyhedron->totallyneutralgrandchildren.end()); |
---|
1322 | |
---|
1323 | |
---|
1324 | |
---|
1325 | // cout << ((toprow*)current_polyhedron)->condition << endl << toadd << endl; |
---|
1326 | |
---|
1327 | toprow* positive_poly = new toprow(((toprow*)current_polyhedron)->condition_sum+toadd, ((toprow*)current_polyhedron)->condition_order+1); |
---|
1328 | toprow* negative_poly = new toprow(((toprow*)current_polyhedron)->condition_sum-toadd, ((toprow*)current_polyhedron)->condition_order+1); |
---|
1329 | |
---|
1330 | positive_poly->my_emlig = this; |
---|
1331 | negative_poly->my_emlig = this; |
---|
1332 | |
---|
1333 | for(set<polyhedron*>::iterator grand_ref = current_polyhedron->totallyneutralgrandchildren.begin(); grand_ref != current_polyhedron->totallyneutralgrandchildren.end();grand_ref++) |
---|
1334 | { |
---|
1335 | (*grand_ref)->parents.push_back(new_totally_neutral_child); |
---|
1336 | (*grand_ref)->grandparents.insert(positive_poly); |
---|
1337 | (*grand_ref)->grandparents.insert(negative_poly); |
---|
1338 | |
---|
1339 | new_totally_neutral_child->vertices.insert((*grand_ref)->vertices.begin(),(*grand_ref)->vertices.end()); |
---|
1340 | } |
---|
1341 | |
---|
1342 | positive_poly->children.push_back(new_totally_neutral_child); |
---|
1343 | negative_poly->children.push_back(new_totally_neutral_child); |
---|
1344 | |
---|
1345 | |
---|
1346 | for(list<polyhedron*>::iterator parent_ref = current_polyhedron->parents.begin();parent_ref!=current_polyhedron->parents.end();parent_ref++) |
---|
1347 | { |
---|
1348 | (*parent_ref)->totallyneutralgrandchildren.insert(new_totally_neutral_child); |
---|
1349 | new_totally_neutral_child->grandparents.insert(*parent_ref); |
---|
1350 | |
---|
1351 | (*parent_ref)->neutralchildren.remove(current_polyhedron); |
---|
1352 | (*parent_ref)->children.remove(current_polyhedron); |
---|
1353 | |
---|
1354 | (*parent_ref)->children.push_back(positive_poly); |
---|
1355 | (*parent_ref)->children.push_back(negative_poly); |
---|
1356 | (*parent_ref)->positivechildren.push_back(positive_poly); |
---|
1357 | (*parent_ref)->negativechildren.push_back(negative_poly); |
---|
1358 | } |
---|
1359 | |
---|
1360 | positive_poly->parents.insert(positive_poly->parents.end(), |
---|
1361 | current_polyhedron->parents.begin(), |
---|
1362 | current_polyhedron->parents.end()); |
---|
1363 | |
---|
1364 | negative_poly->parents.insert(negative_poly->parents.end(), |
---|
1365 | current_polyhedron->parents.begin(), |
---|
1366 | current_polyhedron->parents.end()); |
---|
1367 | |
---|
1368 | |
---|
1369 | |
---|
1370 | new_totally_neutral_child->parents.push_back(positive_poly); |
---|
1371 | new_totally_neutral_child->parents.push_back(negative_poly); |
---|
1372 | |
---|
1373 | for(list<polyhedron*>::iterator child_ref = current_polyhedron->positivechildren.begin();child_ref!=current_polyhedron->positivechildren.end();child_ref++) |
---|
1374 | { |
---|
1375 | (*child_ref)->parents.remove(current_polyhedron); |
---|
1376 | (*child_ref)->parents.push_back(positive_poly); |
---|
1377 | } |
---|
1378 | |
---|
1379 | positive_poly->children.insert(positive_poly->children.end(), |
---|
1380 | current_polyhedron->positivechildren.begin(), |
---|
1381 | current_polyhedron->positivechildren.end()); |
---|
1382 | |
---|
1383 | for(list<polyhedron*>::iterator child_ref = current_polyhedron->negativechildren.begin();child_ref!=current_polyhedron->negativechildren.end();child_ref++) |
---|
1384 | { |
---|
1385 | (*child_ref)->parents.remove(current_polyhedron); |
---|
1386 | (*child_ref)->parents.push_back(negative_poly); |
---|
1387 | } |
---|
1388 | |
---|
1389 | negative_poly->children.insert(negative_poly->children.end(), |
---|
1390 | current_polyhedron->negativechildren.begin(), |
---|
1391 | current_polyhedron->negativechildren.end()); |
---|
1392 | |
---|
1393 | positive_poly->vertices.insert(current_polyhedron->positiveneutralvertices.begin(),current_polyhedron->positiveneutralvertices.end()); |
---|
1394 | positive_poly->vertices.insert(new_totally_neutral_child->vertices.begin(),new_totally_neutral_child->vertices.end()); |
---|
1395 | |
---|
1396 | negative_poly->vertices.insert(current_polyhedron->negativeneutralvertices.begin(),current_polyhedron->negativeneutralvertices.end()); |
---|
1397 | negative_poly->vertices.insert(new_totally_neutral_child->vertices.begin(),new_totally_neutral_child->vertices.end()); |
---|
1398 | |
---|
1399 | new_totally_neutral_child->triangulate(false); |
---|
1400 | |
---|
1401 | positive_poly->triangulate(k==for_splitting.size()-1); |
---|
1402 | negative_poly->triangulate(k==for_splitting.size()-1); |
---|
1403 | |
---|
1404 | statistic.append_polyhedron(k-1, new_totally_neutral_child); |
---|
1405 | |
---|
1406 | statistic.insert_polyhedron(k, positive_poly, current_polyhedron); |
---|
1407 | statistic.insert_polyhedron(k, negative_poly, current_polyhedron); |
---|
1408 | |
---|
1409 | statistic.delete_polyhedron(k, current_polyhedron); |
---|
1410 | |
---|
1411 | delete current_polyhedron; |
---|
1412 | } |
---|
1413 | |
---|
1414 | k++; |
---|
1415 | } |
---|
1416 | } |
---|
1417 | |
---|
1418 | |
---|
1419 | vector<int> sizevector; |
---|
1420 | for(int s = 0;s<statistic.size();s++) |
---|
1421 | { |
---|
1422 | sizevector.push_back(statistic.row_size(s)); |
---|
1423 | } |
---|
1424 | |
---|
1425 | /* |
---|
1426 | for(polyhedron* topr_ref = statistic.rows[statistic.size()-1];topr_ref!=statistic.row_ends[statistic.size()-1]->next_poly;topr_ref=topr_ref->next_poly) |
---|
1427 | { |
---|
1428 | cout << ((toprow*)topr_ref)->condition << endl; |
---|
1429 | } |
---|
1430 | */ |
---|
1431 | |
---|
1432 | } |
---|
1433 | |
---|
1434 | void set_correction_factors(int order) |
---|
1435 | { |
---|
1436 | for(int remaining_order = correction_factors.size();remaining_order<order;remaining_order++) |
---|
1437 | { |
---|
1438 | multiset<my_ivec> factor_templates; |
---|
1439 | multiset<my_ivec> final_factors; |
---|
1440 | |
---|
1441 | my_ivec orig_template = my_ivec(); |
---|
1442 | |
---|
1443 | for(int i = 1;i<number_of_parameters-remaining_order+1;i++) |
---|
1444 | { |
---|
1445 | bool in_cycle = false; |
---|
1446 | for(int j = 0;j<=remaining_order;j++) { |
---|
1447 | |
---|
1448 | multiset<my_ivec>::iterator fac_ref = factor_templates.begin(); |
---|
1449 | |
---|
1450 | do |
---|
1451 | { |
---|
1452 | my_ivec current_template; |
---|
1453 | if(!in_cycle) |
---|
1454 | { |
---|
1455 | current_template = orig_template; |
---|
1456 | in_cycle = true; |
---|
1457 | } |
---|
1458 | else |
---|
1459 | { |
---|
1460 | current_template = (*fac_ref); |
---|
1461 | fac_ref++; |
---|
1462 | } |
---|
1463 | |
---|
1464 | current_template.ins(current_template.size(),i); |
---|
1465 | |
---|
1466 | // cout << "template:" << current_template << endl; |
---|
1467 | |
---|
1468 | if(current_template.size()==remaining_order+1) |
---|
1469 | { |
---|
1470 | final_factors.insert(current_template); |
---|
1471 | } |
---|
1472 | else |
---|
1473 | { |
---|
1474 | factor_templates.insert(current_template); |
---|
1475 | } |
---|
1476 | } |
---|
1477 | while(fac_ref!=factor_templates.end()); |
---|
1478 | } |
---|
1479 | } |
---|
1480 | |
---|
1481 | correction_factors.push_back(final_factors); |
---|
1482 | |
---|
1483 | } |
---|
1484 | } |
---|
1485 | |
---|
1486 | protected: |
---|
1487 | |
---|
1488 | /// A method for creating plain default statistic representing only the range of the parameter space. |
---|
1489 | void create_statistic(int number_of_parameters) |
---|
1490 | { |
---|
1491 | for(int i = 0;i<number_of_parameters;i++) |
---|
1492 | { |
---|
1493 | vec condition_vec = zeros(number_of_parameters+1); |
---|
1494 | condition_vec[i+1] = 1; |
---|
1495 | |
---|
1496 | condition* new_condition = new condition(condition_vec); |
---|
1497 | |
---|
1498 | conditions.push_back(new_condition); |
---|
1499 | } |
---|
1500 | |
---|
1501 | // An empty vector of coordinates. |
---|
1502 | vec origin_coord; |
---|
1503 | |
---|
1504 | // We create an origin - this point will have all the coordinates zero, but now it has an empty vector of coords. |
---|
1505 | vertex *origin = new vertex(origin_coord); |
---|
1506 | |
---|
1507 | origin->my_emlig = this; |
---|
1508 | |
---|
1509 | /* |
---|
1510 | // As a statistic, we have to create a vector of vectors of polyhedron pointers. It will then represent the Hasse |
---|
1511 | // diagram. First we create a vector of polyhedrons.. |
---|
1512 | list<polyhedron*> origin_vec; |
---|
1513 | |
---|
1514 | // ..we fill it with the origin.. |
---|
1515 | origin_vec.push_back(origin); |
---|
1516 | |
---|
1517 | // ..and we fill the statistic with the created vector. |
---|
1518 | statistic.push_back(origin_vec); |
---|
1519 | */ |
---|
1520 | |
---|
1521 | statistic = *(new c_statistic()); |
---|
1522 | |
---|
1523 | statistic.append_polyhedron(0, origin); |
---|
1524 | |
---|
1525 | // Now we have a statistic for a zero dimensional space. Regarding to how many dimensional space we need to |
---|
1526 | // describe, we have to widen the descriptional default statistic. We use an iterative procedure as follows: |
---|
1527 | for(int i=0;i<number_of_parameters;i++) |
---|
1528 | { |
---|
1529 | // We first will create two new vertices. These will be the borders of the parameter space in the dimension |
---|
1530 | // of newly added parameter. Therefore they will have all coordinates except the last one zero. We get the |
---|
1531 | // right amount of zero cooridnates by reading them from the origin |
---|
1532 | vec origin_coord = origin->get_coordinates(); |
---|
1533 | |
---|
1534 | // And we incorporate the nonzero coordinates into the new cooordinate vectors |
---|
1535 | vec origin_coord1 = concat(origin_coord,-max_range); |
---|
1536 | vec origin_coord2 = concat(origin_coord,max_range); |
---|
1537 | |
---|
1538 | |
---|
1539 | // Now we create the points |
---|
1540 | vertex* new_point1 = new vertex(origin_coord1); |
---|
1541 | vertex* new_point2 = new vertex(origin_coord2); |
---|
1542 | |
---|
1543 | new_point1->my_emlig = this; |
---|
1544 | new_point2->my_emlig = this; |
---|
1545 | |
---|
1546 | //********************************************************************************************************* |
---|
1547 | // The algorithm for recursive build of a new Hasse diagram representing the space structure from the old |
---|
1548 | // diagram works so that you create two copies of the old Hasse diagram, you shift them up one level (points |
---|
1549 | // will be segments, segments will be areas etc.) and you connect each one of the original copied polyhedrons |
---|
1550 | // with its offspring by a parent-child relation. Also each of the segments in the first (second) copy is |
---|
1551 | // connected to the first (second) newly created vertex by a parent-child relation. |
---|
1552 | //********************************************************************************************************* |
---|
1553 | |
---|
1554 | |
---|
1555 | /* |
---|
1556 | // Create the vectors of vectors of pointers to polyhedrons to hold the copies of the old Hasse diagram |
---|
1557 | vector<vector<polyhedron*>> new_statistic1; |
---|
1558 | vector<vector<polyhedron*>> new_statistic2; |
---|
1559 | */ |
---|
1560 | |
---|
1561 | c_statistic* new_statistic1 = new c_statistic(); |
---|
1562 | c_statistic* new_statistic2 = new c_statistic(); |
---|
1563 | |
---|
1564 | |
---|
1565 | // Copy the statistic by rows |
---|
1566 | for(int j=0;j<statistic.size();j++) |
---|
1567 | { |
---|
1568 | |
---|
1569 | |
---|
1570 | // an element counter |
---|
1571 | int element_number = 0; |
---|
1572 | |
---|
1573 | /* |
---|
1574 | vector<polyhedron*> supportnew_1; |
---|
1575 | vector<polyhedron*> supportnew_2; |
---|
1576 | |
---|
1577 | new_statistic1.push_back(supportnew_1); |
---|
1578 | new_statistic2.push_back(supportnew_2); |
---|
1579 | */ |
---|
1580 | |
---|
1581 | // for each polyhedron in the given row |
---|
1582 | for(polyhedron* horiz_ref = statistic.rows[j];horiz_ref!=statistic.get_end();horiz_ref=horiz_ref->next_poly) |
---|
1583 | { |
---|
1584 | // Append an extra zero coordinate to each of the vertices for the new dimension |
---|
1585 | // If vert_ref is at the first index => we loop through vertices |
---|
1586 | if(j == 0) |
---|
1587 | { |
---|
1588 | // cast the polyhedron pointer to a vertex pointer and push a zero to its vector of coordinates |
---|
1589 | ((vertex*) horiz_ref)->push_coordinate(0); |
---|
1590 | } |
---|
1591 | /* |
---|
1592 | else |
---|
1593 | { |
---|
1594 | ((toprow*) (*horiz_ref))->condition.ins(0,0); |
---|
1595 | }*/ |
---|
1596 | |
---|
1597 | // if it has parents |
---|
1598 | if(!horiz_ref->parents.empty()) |
---|
1599 | { |
---|
1600 | // save the relative address of this child in a vector kids_rel_addresses of all its parents. |
---|
1601 | // This information will later be used for copying the whole Hasse diagram with each of the |
---|
1602 | // relations contained within. |
---|
1603 | for(list<polyhedron*>::iterator parent_ref = horiz_ref->parents.begin();parent_ref != horiz_ref->parents.end();parent_ref++) |
---|
1604 | { |
---|
1605 | (*parent_ref)->kids_rel_addresses.push_back(element_number); |
---|
1606 | } |
---|
1607 | } |
---|
1608 | |
---|
1609 | // ************************************************************************************************** |
---|
1610 | // Here we begin creating a new polyhedron, which will be a copy of the old one. Each such polyhedron |
---|
1611 | // will be created as a toprow, but this information will be later forgotten and only the polyhedrons |
---|
1612 | // in the top row of the Hasse diagram will be considered toprow for later use. |
---|
1613 | // ************************************************************************************************** |
---|
1614 | |
---|
1615 | // First we create vectors specifying a toprow condition. In the case of a preconstructed statistic |
---|
1616 | // this condition will be a vector of zeros. There are two vectors, because we need two copies of |
---|
1617 | // the original Hasse diagram. |
---|
1618 | vec vec1(number_of_parameters+1); |
---|
1619 | vec1.zeros(); |
---|
1620 | |
---|
1621 | vec vec2(number_of_parameters+1); |
---|
1622 | vec2.zeros(); |
---|
1623 | |
---|
1624 | // We create a new toprow with the previously specified condition. |
---|
1625 | toprow* current_copy1 = new toprow(vec1, 0); |
---|
1626 | toprow* current_copy2 = new toprow(vec2, 0); |
---|
1627 | |
---|
1628 | current_copy1->my_emlig = this; |
---|
1629 | current_copy2->my_emlig = this; |
---|
1630 | |
---|
1631 | // The vertices of the copies will be inherited, because there will be a parent/child relation |
---|
1632 | // between each polyhedron and its offspring (comming from the copy) and a parent has all the |
---|
1633 | // vertices of its child plus more. |
---|
1634 | for(set<vertex*>::iterator vertex_ref = horiz_ref->vertices.begin();vertex_ref!=horiz_ref->vertices.end();vertex_ref++) |
---|
1635 | { |
---|
1636 | current_copy1->vertices.insert(*vertex_ref); |
---|
1637 | current_copy2->vertices.insert(*vertex_ref); |
---|
1638 | } |
---|
1639 | |
---|
1640 | // The only new vertex of the offspring should be the newly created point. |
---|
1641 | current_copy1->vertices.insert(new_point1); |
---|
1642 | current_copy2->vertices.insert(new_point2); |
---|
1643 | |
---|
1644 | // This method guarantees that each polyhedron is already triangulated, therefore its triangulation |
---|
1645 | // is only one set of vertices and it is the set of all its vertices. |
---|
1646 | set<vertex*> t_simplex1; |
---|
1647 | set<vertex*> t_simplex2; |
---|
1648 | |
---|
1649 | t_simplex1.insert(current_copy1->vertices.begin(),current_copy1->vertices.end()); |
---|
1650 | t_simplex2.insert(current_copy2->vertices.begin(),current_copy2->vertices.end()); |
---|
1651 | |
---|
1652 | current_copy1->triangulation.push_back(t_simplex1); |
---|
1653 | current_copy2->triangulation.push_back(t_simplex2); |
---|
1654 | |
---|
1655 | // Now we have copied the polyhedron and we have to copy all of its relations. Because we are copying |
---|
1656 | // in the Hasse diagram from bottom up, we always have to copy the parent/child relations to all the |
---|
1657 | // kids and when we do that and know the child, in the child we will remember the parent we came from. |
---|
1658 | // This way all the parents/children relations are saved in both the parent and the child. |
---|
1659 | if(!horiz_ref->kids_rel_addresses.empty()) |
---|
1660 | { |
---|
1661 | for(list<int>::iterator kid_ref = horiz_ref->kids_rel_addresses.begin();kid_ref!=horiz_ref->kids_rel_addresses.end();kid_ref++) |
---|
1662 | { |
---|
1663 | polyhedron* new_kid1 = new_statistic1->rows[j-1]; |
---|
1664 | polyhedron* new_kid2 = new_statistic2->rows[j-1]; |
---|
1665 | |
---|
1666 | // THIS IS NOT EFFECTIVE: It could be improved by having the list indexed for new_statistic, but |
---|
1667 | // not indexed for statistic. Hopefully this will not cause a big slowdown - happens only offline. |
---|
1668 | if(*kid_ref) |
---|
1669 | { |
---|
1670 | for(int k = 1;k<=(*kid_ref);k++) |
---|
1671 | { |
---|
1672 | new_kid1=new_kid1->next_poly; |
---|
1673 | new_kid2=new_kid2->next_poly; |
---|
1674 | } |
---|
1675 | } |
---|
1676 | |
---|
1677 | // find the child and save the relation to the parent |
---|
1678 | current_copy1->children.push_back(new_kid1); |
---|
1679 | current_copy2->children.push_back(new_kid2); |
---|
1680 | |
---|
1681 | // in the child save the parents' address |
---|
1682 | new_kid1->parents.push_back(current_copy1); |
---|
1683 | new_kid2->parents.push_back(current_copy2); |
---|
1684 | } |
---|
1685 | |
---|
1686 | // Here we clear the parents kids_rel_addresses vector for later use (when we need to widen the |
---|
1687 | // Hasse diagram again) |
---|
1688 | horiz_ref->kids_rel_addresses.clear(); |
---|
1689 | } |
---|
1690 | // If there were no children previously, we are copying a polyhedron that has been a vertex before. |
---|
1691 | // In this case it is a segment now and it will have a relation to its mother (copywise) and to the |
---|
1692 | // newly created point. Here we create the connection to the new point, again from both sides. |
---|
1693 | else |
---|
1694 | { |
---|
1695 | // Add the address of the new point in the former vertex |
---|
1696 | current_copy1->children.push_back(new_point1); |
---|
1697 | current_copy2->children.push_back(new_point2); |
---|
1698 | |
---|
1699 | // Add the address of the former vertex in the new point |
---|
1700 | new_point1->parents.push_back(current_copy1); |
---|
1701 | new_point2->parents.push_back(current_copy2); |
---|
1702 | } |
---|
1703 | |
---|
1704 | // Save the mother in its offspring |
---|
1705 | current_copy1->children.push_back(horiz_ref); |
---|
1706 | current_copy2->children.push_back(horiz_ref); |
---|
1707 | |
---|
1708 | // Save the offspring in its mother |
---|
1709 | horiz_ref->parents.push_back(current_copy1); |
---|
1710 | horiz_ref->parents.push_back(current_copy2); |
---|
1711 | |
---|
1712 | |
---|
1713 | // Add the copies into the relevant statistic. The statistic will later be appended to the previous |
---|
1714 | // Hasse diagram |
---|
1715 | new_statistic1->append_polyhedron(j,current_copy1); |
---|
1716 | new_statistic2->append_polyhedron(j,current_copy2); |
---|
1717 | |
---|
1718 | // Raise the count in the vector of polyhedrons |
---|
1719 | element_number++; |
---|
1720 | |
---|
1721 | } |
---|
1722 | |
---|
1723 | } |
---|
1724 | |
---|
1725 | /* |
---|
1726 | statistic.begin()->push_back(new_point1); |
---|
1727 | statistic.begin()->push_back(new_point2); |
---|
1728 | */ |
---|
1729 | |
---|
1730 | statistic.append_polyhedron(0, new_point1); |
---|
1731 | statistic.append_polyhedron(0, new_point2); |
---|
1732 | |
---|
1733 | // Merge the new statistics into the old one. This will either be the final statistic or we will |
---|
1734 | // reenter the widening loop. |
---|
1735 | for(int j=0;j<new_statistic1->size();j++) |
---|
1736 | { |
---|
1737 | /* |
---|
1738 | if(j+1==statistic.size()) |
---|
1739 | { |
---|
1740 | list<polyhedron*> support; |
---|
1741 | statistic.push_back(support); |
---|
1742 | } |
---|
1743 | |
---|
1744 | (statistic.begin()+j+1)->insert((statistic.begin()+j+1)->end(),new_statistic1[j].begin(),new_statistic1[j].end()); |
---|
1745 | (statistic.begin()+j+1)->insert((statistic.begin()+j+1)->end(),new_statistic2[j].begin(),new_statistic2[j].end()); |
---|
1746 | */ |
---|
1747 | statistic.append_polyhedron(j+1,new_statistic1->rows[j],new_statistic1->row_ends[j]); |
---|
1748 | statistic.append_polyhedron(j+1,new_statistic2->rows[j],new_statistic2->row_ends[j]); |
---|
1749 | } |
---|
1750 | } |
---|
1751 | |
---|
1752 | /* |
---|
1753 | vector<list<toprow*>> toprow_statistic; |
---|
1754 | int line_count = 0; |
---|
1755 | |
---|
1756 | for(vector<list<polyhedron*>>::iterator polyhedron_ref = ++statistic.begin(); polyhedron_ref!=statistic.end();polyhedron_ref++) |
---|
1757 | { |
---|
1758 | list<toprow*> support_list; |
---|
1759 | toprow_statistic.push_back(support_list); |
---|
1760 | |
---|
1761 | for(list<polyhedron*>::iterator polyhedron_ref2 = polyhedron_ref->begin(); polyhedron_ref2 != polyhedron_ref->end(); polyhedron_ref2++) |
---|
1762 | { |
---|
1763 | toprow* support_top = (toprow*)(*polyhedron_ref2); |
---|
1764 | |
---|
1765 | toprow_statistic[line_count].push_back(support_top); |
---|
1766 | } |
---|
1767 | |
---|
1768 | line_count++; |
---|
1769 | }*/ |
---|
1770 | |
---|
1771 | /* |
---|
1772 | vector<int> sizevector; |
---|
1773 | for(int s = 0;s<statistic.size();s++) |
---|
1774 | { |
---|
1775 | sizevector.push_back(statistic.row_size(s)); |
---|
1776 | } |
---|
1777 | */ |
---|
1778 | |
---|
1779 | } |
---|
1780 | |
---|
1781 | |
---|
1782 | |
---|
1783 | |
---|
1784 | }; |
---|
1785 | |
---|
1786 | |
---|
1787 | |
---|
1788 | //! Robust Bayesian AR model for Multicriteria-Laplace-Inverse-Gamma density |
---|
1789 | class RARX //: public BM |
---|
1790 | { |
---|
1791 | private: |
---|
1792 | |
---|
1793 | emlig* posterior; |
---|
1794 | |
---|
1795 | int window_size; |
---|
1796 | |
---|
1797 | public: |
---|
1798 | RARX(int number_of_parameters, const int window_size)//:BM() |
---|
1799 | { |
---|
1800 | posterior = new emlig(number_of_parameters); |
---|
1801 | |
---|
1802 | this->window_size = window_size; |
---|
1803 | }; |
---|
1804 | |
---|
1805 | void bayes(const itpp::vec &yt, const itpp::vec &cond = "") |
---|
1806 | { |
---|
1807 | |
---|
1808 | } |
---|
1809 | |
---|
1810 | }; |
---|
1811 | |
---|
1812 | |
---|
1813 | |
---|
1814 | #endif //TRAGE_H |
---|