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 | /// An enumeration of possible actions performed on the polyhedrons. We can merge them or split them. |
---|
26 | enum actions {MERGE, SPLIT}; |
---|
27 | |
---|
28 | // Forward declaration of polyhedron, vertex and emlig |
---|
29 | class polyhedron; |
---|
30 | class vertex; |
---|
31 | class emlig; |
---|
32 | |
---|
33 | |
---|
34 | /* |
---|
35 | class t_simplex |
---|
36 | { |
---|
37 | public: |
---|
38 | set<vertex*> minima; |
---|
39 | |
---|
40 | set<vertex*> simplex; |
---|
41 | |
---|
42 | t_simplex(vertex* origin_vertex) |
---|
43 | { |
---|
44 | simplex.insert(origin_vertex); |
---|
45 | minima.insert(origin_vertex); |
---|
46 | } |
---|
47 | };*/ |
---|
48 | |
---|
49 | /// A class representing a single condition that can be added to the emlig. A condition represents data entries in a statistical model. |
---|
50 | class condition |
---|
51 | { |
---|
52 | public: |
---|
53 | /// Value of the condition representing the data |
---|
54 | vec value; |
---|
55 | |
---|
56 | /// Mulitplicity of the given condition may represent multiple occurences of same data entry. |
---|
57 | int multiplicity; |
---|
58 | |
---|
59 | /// Default constructor of condition class takes the value of data entry and creates a condition with multiplicity 1 (first occurence of the data). |
---|
60 | condition(vec value) |
---|
61 | { |
---|
62 | this->value = value; |
---|
63 | multiplicity = 1; |
---|
64 | } |
---|
65 | }; |
---|
66 | |
---|
67 | class simplex |
---|
68 | { |
---|
69 | public: |
---|
70 | |
---|
71 | set<vertex*> vertices; |
---|
72 | |
---|
73 | double probability; |
---|
74 | |
---|
75 | vector<multimap<double,double>> gamma_parameters; |
---|
76 | |
---|
77 | double gamma_sum; |
---|
78 | |
---|
79 | |
---|
80 | simplex(set<vertex*> vertices) |
---|
81 | { |
---|
82 | this->vertices.insert(vertices.begin(),vertices.end()); |
---|
83 | probability = 0; |
---|
84 | } |
---|
85 | |
---|
86 | simplex(vertex* vertex) |
---|
87 | { |
---|
88 | this->vertices.insert(vertex); |
---|
89 | probability = 0; |
---|
90 | } |
---|
91 | }; |
---|
92 | |
---|
93 | |
---|
94 | /// A class describing a single polyhedron of the split complex. From a collection of such classes a Hasse diagram |
---|
95 | /// of the structure in the exponent of a Laplace-Inverse-Gamma density will be created. |
---|
96 | class polyhedron |
---|
97 | { |
---|
98 | /// A property having a value of 1 usually, with higher value only if the polyhedron arises as a coincidence of |
---|
99 | /// more than just the necessary number of conditions. For example if a newly created line passes through an already |
---|
100 | /// existing point, the points multiplicity will rise by 1. |
---|
101 | int multiplicity; |
---|
102 | |
---|
103 | /// A property representing the position of the polyhedron related to current condition with relation to which we |
---|
104 | /// are splitting the parameter space (new data has arrived). This property is setup within a classification procedure and |
---|
105 | /// is only valid while the new condition is being added. It has to be reset when new condition is added and new classification |
---|
106 | /// has to be performed. |
---|
107 | int split_state; |
---|
108 | |
---|
109 | /// A property representing the position of the polyhedron related to current condition with relation to which we |
---|
110 | /// are merging the parameter space (data is being deleted usually due to a moving window model which is more adaptive and |
---|
111 | /// steps in for the forgetting in a classical Gaussian AR model). This property is setup within a classification procedure and |
---|
112 | /// is only valid while the new condition is being removed. It has to be reset when new condition is removed and new classification |
---|
113 | /// has to be performed. |
---|
114 | int merge_state; |
---|
115 | |
---|
116 | |
---|
117 | |
---|
118 | public: |
---|
119 | /// A pointer to the multi-Laplace inverse gamma distribution this polyhedron belongs to. |
---|
120 | emlig* my_emlig; |
---|
121 | |
---|
122 | /// A list of polyhedrons parents within the Hasse diagram. |
---|
123 | list<polyhedron*> parents; |
---|
124 | |
---|
125 | /// A list of polyhedrons children withing the Hasse diagram. |
---|
126 | list<polyhedron*> children; |
---|
127 | |
---|
128 | /// All the vertices of the given polyhedron |
---|
129 | set<vertex*> vertices; |
---|
130 | |
---|
131 | /// The conditions that gave birth to the polyhedron. If some of them is removed, the polyhedron ceases to exist. |
---|
132 | set<condition*> parentconditions; |
---|
133 | |
---|
134 | /// A list used for storing children that lie in the positive region related to a certain condition |
---|
135 | list<polyhedron*> positivechildren; |
---|
136 | |
---|
137 | /// A list used for storing children that lie in the negative region related to a certain condition |
---|
138 | list<polyhedron*> negativechildren; |
---|
139 | |
---|
140 | /// Children intersecting the condition |
---|
141 | list<polyhedron*> neutralchildren; |
---|
142 | |
---|
143 | /// A set of grandchildren of the polyhedron that when new condition is added lie exactly on the condition hyperplane. These grandchildren |
---|
144 | /// behave differently from other grandchildren, when the polyhedron is split. New grandchild is not necessarily created on the crossection of |
---|
145 | /// the polyhedron and new condition. |
---|
146 | set<polyhedron*> totallyneutralgrandchildren; |
---|
147 | |
---|
148 | /// A set of children of the polyhedron that when new condition is added lie exactly on the condition hyperplane. These children |
---|
149 | /// behave differently from other children, when the polyhedron is split. New child is not necessarily created on the crossection of |
---|
150 | /// the polyhedron and new condition. |
---|
151 | set<polyhedron*> totallyneutralchildren; |
---|
152 | |
---|
153 | /// Reverse relation to the totallyneutralgrandchildren set is needed for merging of already existing polyhedrons to keep |
---|
154 | /// totallyneutralgrandchildren list up to date. |
---|
155 | set<polyhedron*> grandparents; |
---|
156 | |
---|
157 | /// Vertices of the polyhedron classified as positive related to an added condition. When the polyhderon is split by the new condition, |
---|
158 | /// these vertices will belong to the positive part of the splitted polyhedron. |
---|
159 | set<vertex*> positiveneutralvertices; |
---|
160 | |
---|
161 | /// Vertices of the polyhedron classified as negative related to an added condition. When the polyhderon is split by the new condition, |
---|
162 | /// these vertices will belong to the negative part of the splitted polyhedron. |
---|
163 | set<vertex*> negativeneutralvertices; |
---|
164 | |
---|
165 | /// A bool specifying if the polyhedron lies exactly on the newly added condition or not. |
---|
166 | bool totally_neutral; |
---|
167 | |
---|
168 | /// When two polyhedrons are merged, there always exists a child lying on the former border of the polyhedrons. This child manages the merge |
---|
169 | /// of the two polyhedrons. This property gives us the address of the mediator child. |
---|
170 | polyhedron* mergechild; |
---|
171 | |
---|
172 | /// If the polyhedron serves as a mergechild for two of its parents, we need to have the address of the parents to access them. This |
---|
173 | /// is the pointer to the positive parent being merged. |
---|
174 | polyhedron* positiveparent; |
---|
175 | |
---|
176 | /// If the polyhedron serves as a mergechild for two of its parents, we need to have the address of the parents to access them. This |
---|
177 | /// is the pointer to the negative parent being merged. |
---|
178 | polyhedron* negativeparent; |
---|
179 | |
---|
180 | /// Adressing withing the statistic. Next_poly is a pointer to the next polyhedron in the statistic on the same level (if this is a point, |
---|
181 | /// next_poly will be a point etc.). |
---|
182 | polyhedron* next_poly; |
---|
183 | |
---|
184 | /// Adressing withing the statistic. Prev_poly is a pointer to the previous polyhedron in the statistic on the same level (if this is a point, |
---|
185 | /// next_poly will be a point etc.). |
---|
186 | polyhedron* prev_poly; |
---|
187 | |
---|
188 | /// A property counting the number of messages obtained from children within a classification procedure of position of the polyhedron related |
---|
189 | /// an added/removed condition. If the message counter reaches the number of children, we know the polyhedrons' position has been fully classified. |
---|
190 | int message_counter; |
---|
191 | |
---|
192 | /// List of triangulation polyhedrons of the polyhedron given by their relative vertices. |
---|
193 | set<simplex*> triangulation; |
---|
194 | |
---|
195 | /// A list of relative addresses serving for Hasse diagram construction. |
---|
196 | list<int> kids_rel_addresses; |
---|
197 | |
---|
198 | /// Default constructor |
---|
199 | polyhedron() |
---|
200 | { |
---|
201 | multiplicity = 1; |
---|
202 | |
---|
203 | message_counter = 0; |
---|
204 | |
---|
205 | totally_neutral = NULL; |
---|
206 | |
---|
207 | mergechild = NULL; |
---|
208 | } |
---|
209 | |
---|
210 | /// Setter for raising multiplicity |
---|
211 | void raise_multiplicity() |
---|
212 | { |
---|
213 | multiplicity++; |
---|
214 | } |
---|
215 | |
---|
216 | /// Setter for lowering multiplicity |
---|
217 | void lower_multiplicity() |
---|
218 | { |
---|
219 | multiplicity--; |
---|
220 | } |
---|
221 | |
---|
222 | int get_multiplicity() |
---|
223 | { |
---|
224 | return multiplicity; |
---|
225 | } |
---|
226 | |
---|
227 | /// An obligatory operator, when the class is used within a C++ STL structure like a vector |
---|
228 | int operator==(polyhedron polyhedron2) |
---|
229 | { |
---|
230 | return true; |
---|
231 | } |
---|
232 | |
---|
233 | /// An obligatory operator, when the class is used within a C++ STL structure like a vector |
---|
234 | int operator<(polyhedron polyhedron2) |
---|
235 | { |
---|
236 | return false; |
---|
237 | } |
---|
238 | |
---|
239 | |
---|
240 | /// A setter of state of current polyhedron relative to the action specified in the argument. The three possible states of the |
---|
241 | /// polyhedron are -1 - NEGATIVE, 0 - NEUTRAL, 1 - POSITIVE. Neutral state means that either the state has been reset or the polyhedron is |
---|
242 | /// ready to be split/merged. |
---|
243 | int set_state(double state_indicator, actions action) |
---|
244 | { |
---|
245 | switch(action) |
---|
246 | { |
---|
247 | case MERGE: |
---|
248 | merge_state = (int)sign(state_indicator); |
---|
249 | return merge_state; |
---|
250 | case SPLIT: |
---|
251 | split_state = (int)sign(state_indicator); |
---|
252 | return split_state; |
---|
253 | } |
---|
254 | } |
---|
255 | |
---|
256 | /// A getter of state of current polyhedron relative to the action specified in the argument. The three possible states of the |
---|
257 | /// polyhedron are -1 - NEGATIVE, 0 - NEUTRAL, 1 - POSITIVE. Neutral state means that either the state has been reset or the polyhedron is |
---|
258 | /// ready to be split/merged. |
---|
259 | int get_state(actions action) |
---|
260 | { |
---|
261 | switch(action) |
---|
262 | { |
---|
263 | case MERGE: |
---|
264 | return merge_state; |
---|
265 | break; |
---|
266 | case SPLIT: |
---|
267 | return split_state; |
---|
268 | break; |
---|
269 | } |
---|
270 | } |
---|
271 | |
---|
272 | /// Method for obtaining the number of children of given polyhedron. |
---|
273 | int number_of_children() |
---|
274 | { |
---|
275 | return children.size(); |
---|
276 | } |
---|
277 | |
---|
278 | /// A method for triangulation of given polyhedron. |
---|
279 | void triangulate(bool should_integrate); |
---|
280 | }; |
---|
281 | |
---|
282 | |
---|
283 | /// A class for representing 0-dimensional polyhedron - a vertex. It will be located in the bottom row of the Hasse |
---|
284 | /// diagram representing a complex of polyhedrons. It has its coordinates in the parameter space. |
---|
285 | class vertex : public polyhedron |
---|
286 | { |
---|
287 | /// A dynamic array representing coordinates of the vertex |
---|
288 | vec coordinates; |
---|
289 | |
---|
290 | public: |
---|
291 | /// A property specifying the value of the density (ted nevim, jestli je to jakoby log nebo ne) above the vertex. |
---|
292 | double function_value; |
---|
293 | |
---|
294 | /// Default constructor |
---|
295 | vertex(); |
---|
296 | |
---|
297 | /// Constructor of a vertex from a set of coordinates |
---|
298 | vertex(vec coordinates) |
---|
299 | { |
---|
300 | this->coordinates = coordinates; |
---|
301 | |
---|
302 | vertices.insert(this); |
---|
303 | |
---|
304 | simplex* vert_simplex = new simplex(vertices); |
---|
305 | |
---|
306 | triangulation.insert(vert_simplex); |
---|
307 | } |
---|
308 | |
---|
309 | /// A method that widens the set of coordinates of given vertex. It is used when a complex in a parameter |
---|
310 | /// space of certain dimension is established, but the dimension is not known when the vertex is created. |
---|
311 | void push_coordinate(double coordinate) |
---|
312 | { |
---|
313 | coordinates = concat(coordinates,coordinate); |
---|
314 | } |
---|
315 | |
---|
316 | /// A method obtaining the set of coordinates of a vertex. These coordinates are not obtained as a pointer |
---|
317 | /// (not given by reference), but a new copy is created (they are given by value). |
---|
318 | vec get_coordinates() |
---|
319 | { |
---|
320 | return coordinates; |
---|
321 | } |
---|
322 | |
---|
323 | }; |
---|
324 | |
---|
325 | |
---|
326 | /// A class representing a polyhedron in a top row of the complex. Such polyhedron has a condition that differen tiates |
---|
327 | /// it from polyhedrons in other rows. |
---|
328 | class toprow : public polyhedron |
---|
329 | { |
---|
330 | |
---|
331 | public: |
---|
332 | double probability; |
---|
333 | |
---|
334 | vertex* minimal_vertex; |
---|
335 | |
---|
336 | /// A condition used for determining the function of a Laplace-Inverse-Gamma density resulting from Bayesian estimation |
---|
337 | vec condition_sum; |
---|
338 | |
---|
339 | int condition_order; |
---|
340 | |
---|
341 | /// Default constructor |
---|
342 | toprow(){}; |
---|
343 | |
---|
344 | /// Constructor creating a toprow from the condition |
---|
345 | toprow(condition *condition, int condition_order) |
---|
346 | { |
---|
347 | this->condition_sum = condition->value; |
---|
348 | this->condition_order = condition_order; |
---|
349 | } |
---|
350 | |
---|
351 | toprow(vec condition_sum, int condition_order) |
---|
352 | { |
---|
353 | this->condition_sum = condition_sum; |
---|
354 | this->condition_order = condition_order; |
---|
355 | } |
---|
356 | |
---|
357 | double integrate_simplex(simplex* simplex, char c); |
---|
358 | |
---|
359 | }; |
---|
360 | |
---|
361 | |
---|
362 | |
---|
363 | |
---|
364 | |
---|
365 | |
---|
366 | |
---|
367 | class c_statistic |
---|
368 | { |
---|
369 | |
---|
370 | public: |
---|
371 | polyhedron* end_poly; |
---|
372 | polyhedron* start_poly; |
---|
373 | |
---|
374 | vector<polyhedron*> rows; |
---|
375 | |
---|
376 | vector<polyhedron*> row_ends; |
---|
377 | |
---|
378 | c_statistic() |
---|
379 | { |
---|
380 | end_poly = new polyhedron(); |
---|
381 | start_poly = new polyhedron(); |
---|
382 | }; |
---|
383 | |
---|
384 | void append_polyhedron(int row, polyhedron* appended_start, polyhedron* appended_end) |
---|
385 | { |
---|
386 | if(row>((int)rows.size())-1) |
---|
387 | { |
---|
388 | if(row>rows.size()) |
---|
389 | { |
---|
390 | throw new exception("You are trying to append a polyhedron whose children are not in the statistic yet!"); |
---|
391 | return; |
---|
392 | } |
---|
393 | |
---|
394 | rows.push_back(end_poly); |
---|
395 | row_ends.push_back(end_poly); |
---|
396 | } |
---|
397 | |
---|
398 | // POSSIBLE FAILURE: the function is not checking if start and end are connected |
---|
399 | |
---|
400 | if(rows[row] != end_poly) |
---|
401 | { |
---|
402 | appended_start->prev_poly = row_ends[row]; |
---|
403 | row_ends[row]->next_poly = appended_start; |
---|
404 | |
---|
405 | } |
---|
406 | else if((row>0 && rows[row-1]!=end_poly)||row==0) |
---|
407 | { |
---|
408 | appended_start->prev_poly = start_poly; |
---|
409 | rows[row]= appended_start; |
---|
410 | } |
---|
411 | else |
---|
412 | { |
---|
413 | throw new exception("Wrong polyhedron insertion into statistic: missing intermediary polyhedron!"); |
---|
414 | } |
---|
415 | |
---|
416 | appended_end->next_poly = end_poly; |
---|
417 | row_ends[row] = appended_end; |
---|
418 | } |
---|
419 | |
---|
420 | void append_polyhedron(int row, polyhedron* appended_poly) |
---|
421 | { |
---|
422 | append_polyhedron(row,appended_poly,appended_poly); |
---|
423 | } |
---|
424 | |
---|
425 | void insert_polyhedron(int row, polyhedron* inserted_poly, polyhedron* following_poly) |
---|
426 | { |
---|
427 | if(following_poly != end_poly) |
---|
428 | { |
---|
429 | inserted_poly->next_poly = following_poly; |
---|
430 | inserted_poly->prev_poly = following_poly->prev_poly; |
---|
431 | |
---|
432 | if(following_poly->prev_poly == start_poly) |
---|
433 | { |
---|
434 | rows[row] = inserted_poly; |
---|
435 | } |
---|
436 | else |
---|
437 | { |
---|
438 | inserted_poly->prev_poly->next_poly = inserted_poly; |
---|
439 | } |
---|
440 | |
---|
441 | following_poly->prev_poly = inserted_poly; |
---|
442 | } |
---|
443 | else |
---|
444 | { |
---|
445 | this->append_polyhedron(row, inserted_poly); |
---|
446 | } |
---|
447 | |
---|
448 | } |
---|
449 | |
---|
450 | |
---|
451 | |
---|
452 | |
---|
453 | void delete_polyhedron(int row, polyhedron* deleted_poly) |
---|
454 | { |
---|
455 | if(deleted_poly->prev_poly != start_poly) |
---|
456 | { |
---|
457 | deleted_poly->prev_poly->next_poly = deleted_poly->next_poly; |
---|
458 | } |
---|
459 | else |
---|
460 | { |
---|
461 | rows[row] = deleted_poly->next_poly; |
---|
462 | } |
---|
463 | |
---|
464 | if(deleted_poly->next_poly!=end_poly) |
---|
465 | { |
---|
466 | deleted_poly->next_poly->prev_poly = deleted_poly->prev_poly; |
---|
467 | } |
---|
468 | else |
---|
469 | { |
---|
470 | row_ends[row] = deleted_poly->prev_poly; |
---|
471 | } |
---|
472 | |
---|
473 | |
---|
474 | |
---|
475 | deleted_poly->next_poly = NULL; |
---|
476 | deleted_poly->prev_poly = NULL; |
---|
477 | } |
---|
478 | |
---|
479 | int size() |
---|
480 | { |
---|
481 | return rows.size(); |
---|
482 | } |
---|
483 | |
---|
484 | polyhedron* get_end() |
---|
485 | { |
---|
486 | return end_poly; |
---|
487 | } |
---|
488 | |
---|
489 | polyhedron* get_start() |
---|
490 | { |
---|
491 | return start_poly; |
---|
492 | } |
---|
493 | |
---|
494 | int row_size(int row) |
---|
495 | { |
---|
496 | if(this->size()>row && row>=0) |
---|
497 | { |
---|
498 | int row_size = 0; |
---|
499 | |
---|
500 | for(polyhedron* row_poly = rows[row]; row_poly!=end_poly; row_poly=row_poly->next_poly) |
---|
501 | { |
---|
502 | row_size++; |
---|
503 | } |
---|
504 | |
---|
505 | return row_size; |
---|
506 | } |
---|
507 | else |
---|
508 | { |
---|
509 | throw new exception("There is no row to obtain size from!"); |
---|
510 | } |
---|
511 | } |
---|
512 | }; |
---|
513 | |
---|
514 | |
---|
515 | class my_ivec : public ivec |
---|
516 | { |
---|
517 | public: |
---|
518 | my_ivec():ivec(){}; |
---|
519 | |
---|
520 | my_ivec(ivec origin):ivec() |
---|
521 | { |
---|
522 | this->ins(0,origin); |
---|
523 | } |
---|
524 | |
---|
525 | bool operator>(const my_ivec &second) const |
---|
526 | { |
---|
527 | return max(*this)>max(second); |
---|
528 | |
---|
529 | /* |
---|
530 | int size1 = this->size(); |
---|
531 | int size2 = second.size(); |
---|
532 | |
---|
533 | int counter1 = 0; |
---|
534 | while(0==0) |
---|
535 | { |
---|
536 | if((*this)[counter1]==0) |
---|
537 | { |
---|
538 | size1--; |
---|
539 | } |
---|
540 | |
---|
541 | if((*this)[counter1]!=0) |
---|
542 | break; |
---|
543 | |
---|
544 | counter1++; |
---|
545 | } |
---|
546 | |
---|
547 | int counter2 = 0; |
---|
548 | while(0==0) |
---|
549 | { |
---|
550 | if(second[counter2]==0) |
---|
551 | { |
---|
552 | size2--; |
---|
553 | } |
---|
554 | |
---|
555 | if(second[counter2]!=0) |
---|
556 | break; |
---|
557 | |
---|
558 | counter2++; |
---|
559 | } |
---|
560 | |
---|
561 | if(size1!=size2) |
---|
562 | { |
---|
563 | return size1>size2; |
---|
564 | } |
---|
565 | else |
---|
566 | { |
---|
567 | for(int i = 0;i<size1;i++) |
---|
568 | { |
---|
569 | if((*this)[counter1+i]!=second[counter2+i]) |
---|
570 | { |
---|
571 | return (*this)[counter1+i]>second[counter2+i]; |
---|
572 | } |
---|
573 | } |
---|
574 | |
---|
575 | return false; |
---|
576 | }*/ |
---|
577 | } |
---|
578 | |
---|
579 | |
---|
580 | bool operator==(const my_ivec &second) const |
---|
581 | { |
---|
582 | return max(*this)==max(second); |
---|
583 | |
---|
584 | /* |
---|
585 | int size1 = this->size(); |
---|
586 | int size2 = second.size(); |
---|
587 | |
---|
588 | int counter = 0; |
---|
589 | while(0==0) |
---|
590 | { |
---|
591 | if((*this)[counter]==0) |
---|
592 | { |
---|
593 | size1--; |
---|
594 | } |
---|
595 | |
---|
596 | if((*this)[counter]!=0) |
---|
597 | break; |
---|
598 | |
---|
599 | counter++; |
---|
600 | } |
---|
601 | |
---|
602 | counter = 0; |
---|
603 | while(0==0) |
---|
604 | { |
---|
605 | if(second[counter]==0) |
---|
606 | { |
---|
607 | size2--; |
---|
608 | } |
---|
609 | |
---|
610 | if(second[counter]!=0) |
---|
611 | break; |
---|
612 | |
---|
613 | counter++; |
---|
614 | } |
---|
615 | |
---|
616 | if(size1!=size2) |
---|
617 | { |
---|
618 | return false; |
---|
619 | } |
---|
620 | else |
---|
621 | { |
---|
622 | for(int i=0;i<size1;i++) |
---|
623 | { |
---|
624 | if((*this)[size()-1-i]!=second[second.size()-1-i]) |
---|
625 | { |
---|
626 | return false; |
---|
627 | } |
---|
628 | } |
---|
629 | |
---|
630 | return true; |
---|
631 | }*/ |
---|
632 | } |
---|
633 | |
---|
634 | bool operator<(const my_ivec &second) const |
---|
635 | { |
---|
636 | return !(((*this)>second)||((*this)==second)); |
---|
637 | } |
---|
638 | |
---|
639 | bool operator!=(const my_ivec &second) const |
---|
640 | { |
---|
641 | return !((*this)==second); |
---|
642 | } |
---|
643 | |
---|
644 | bool operator<=(const my_ivec &second) const |
---|
645 | { |
---|
646 | return !((*this)>second); |
---|
647 | } |
---|
648 | |
---|
649 | bool operator>=(const my_ivec &second) const |
---|
650 | { |
---|
651 | return !((*this)<second); |
---|
652 | } |
---|
653 | |
---|
654 | my_ivec right(my_ivec original) |
---|
655 | { |
---|
656 | |
---|
657 | } |
---|
658 | }; |
---|
659 | |
---|
660 | |
---|
661 | |
---|
662 | |
---|
663 | |
---|
664 | |
---|
665 | |
---|
666 | //! Conditional(e) Multicriteria-Laplace-Inverse-Gamma distribution density |
---|
667 | class emlig // : eEF |
---|
668 | { |
---|
669 | |
---|
670 | /// A statistic in a form of a Hasse diagram representing a complex of convex polyhedrons obtained as a result |
---|
671 | /// of data update from Bayesian estimation or set by the user if this emlig is a prior density |
---|
672 | |
---|
673 | |
---|
674 | vector<list<polyhedron*>> for_splitting; |
---|
675 | |
---|
676 | vector<list<polyhedron*>> for_merging; |
---|
677 | |
---|
678 | list<condition*> conditions; |
---|
679 | |
---|
680 | double normalization_factor; |
---|
681 | |
---|
682 | |
---|
683 | |
---|
684 | void alter_toprow_conditions(condition *condition, bool should_be_added) |
---|
685 | { |
---|
686 | for(polyhedron* horiz_ref = statistic.rows[statistic.size()-1];horiz_ref!=statistic.get_end();horiz_ref=horiz_ref->next_poly) |
---|
687 | { |
---|
688 | set<vertex*>::iterator vertex_ref = horiz_ref->vertices.begin(); |
---|
689 | |
---|
690 | do |
---|
691 | { |
---|
692 | vertex_ref++; |
---|
693 | } |
---|
694 | while((*vertex_ref)->parentconditions.find(condition)==(*vertex_ref)->parentconditions.end()); |
---|
695 | |
---|
696 | double product = (*vertex_ref)->get_coordinates()*condition->value; |
---|
697 | |
---|
698 | if(should_be_added) |
---|
699 | { |
---|
700 | ((toprow*) horiz_ref)->condition_order++; |
---|
701 | |
---|
702 | if(product>0) |
---|
703 | { |
---|
704 | ((toprow*) horiz_ref)->condition_sum += condition->value; |
---|
705 | } |
---|
706 | else |
---|
707 | { |
---|
708 | ((toprow*) horiz_ref)->condition_sum -= condition->value; |
---|
709 | } |
---|
710 | } |
---|
711 | else |
---|
712 | { |
---|
713 | ((toprow*) horiz_ref)->condition_order--; |
---|
714 | |
---|
715 | if(product<0) |
---|
716 | { |
---|
717 | ((toprow*) horiz_ref)->condition_sum += condition->value; |
---|
718 | } |
---|
719 | else |
---|
720 | { |
---|
721 | ((toprow*) horiz_ref)->condition_sum -= condition->value; |
---|
722 | } |
---|
723 | } |
---|
724 | } |
---|
725 | } |
---|
726 | |
---|
727 | |
---|
728 | |
---|
729 | void send_state_message(polyhedron* sender, condition *toadd, condition *toremove, int level) |
---|
730 | { |
---|
731 | |
---|
732 | bool shouldmerge = (toremove != NULL); |
---|
733 | bool shouldsplit = (toadd != NULL); |
---|
734 | |
---|
735 | if(shouldsplit||shouldmerge) |
---|
736 | { |
---|
737 | for(list<polyhedron*>::iterator parent_iterator = sender->parents.begin();parent_iterator!=sender->parents.end();parent_iterator++) |
---|
738 | { |
---|
739 | polyhedron* current_parent = *parent_iterator; |
---|
740 | |
---|
741 | current_parent->message_counter++; |
---|
742 | |
---|
743 | bool is_last = (current_parent->message_counter == current_parent->number_of_children()); |
---|
744 | bool is_first = (current_parent->message_counter == 1); |
---|
745 | |
---|
746 | if(shouldmerge) |
---|
747 | { |
---|
748 | int child_state = sender->get_state(MERGE); |
---|
749 | int parent_state = current_parent->get_state(MERGE); |
---|
750 | |
---|
751 | if(parent_state == 0||is_first) |
---|
752 | { |
---|
753 | parent_state = current_parent->set_state(child_state, MERGE); |
---|
754 | } |
---|
755 | |
---|
756 | if(child_state == 0) |
---|
757 | { |
---|
758 | if(current_parent->mergechild == NULL) |
---|
759 | { |
---|
760 | current_parent->mergechild = sender; |
---|
761 | } |
---|
762 | } |
---|
763 | |
---|
764 | if(is_last) |
---|
765 | { |
---|
766 | if(parent_state == 1) |
---|
767 | { |
---|
768 | ((toprow*)current_parent)->condition_sum-=toremove->value; |
---|
769 | ((toprow*)current_parent)->condition_order--; |
---|
770 | } |
---|
771 | |
---|
772 | if(parent_state == -1) |
---|
773 | { |
---|
774 | ((toprow*)current_parent)->condition_sum+=toremove->value; |
---|
775 | ((toprow*)current_parent)->condition_order--; |
---|
776 | } |
---|
777 | |
---|
778 | if(current_parent->mergechild != NULL) |
---|
779 | { |
---|
780 | if(current_parent->mergechild->get_multiplicity()==1) |
---|
781 | { |
---|
782 | if(parent_state > 0) |
---|
783 | { |
---|
784 | current_parent->mergechild->positiveparent = current_parent; |
---|
785 | } |
---|
786 | |
---|
787 | if(parent_state < 0) |
---|
788 | { |
---|
789 | current_parent->mergechild->negativeparent = current_parent; |
---|
790 | } |
---|
791 | } |
---|
792 | } |
---|
793 | else |
---|
794 | { |
---|
795 | //current_parent->set_state(0,MERGE); |
---|
796 | |
---|
797 | if(level == number_of_parameters - 1) |
---|
798 | { |
---|
799 | toprow* cur_par_toprow = ((toprow*)current_parent); |
---|
800 | cur_par_toprow->probability = 0.0; |
---|
801 | |
---|
802 | //set<simplex*> new_triangulation; |
---|
803 | |
---|
804 | for(set<simplex*>::iterator s_ref = current_parent->triangulation.begin();s_ref!=current_parent->triangulation.end();s_ref++) |
---|
805 | { |
---|
806 | double cur_prob = cur_par_toprow->integrate_simplex((*s_ref),'C'); |
---|
807 | |
---|
808 | cur_par_toprow->probability += cur_prob; |
---|
809 | |
---|
810 | //new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second)); |
---|
811 | } |
---|
812 | |
---|
813 | //current_parent->triangulation.clear(); |
---|
814 | //current_parent->triangulation.insert(new_triangulation.begin(),new_triangulation.end()); |
---|
815 | } |
---|
816 | } |
---|
817 | |
---|
818 | if(parent_state == 0) |
---|
819 | { |
---|
820 | for_merging[level+1].push_back(current_parent); |
---|
821 | //current_parent->parentconditions.erase(toremove); |
---|
822 | } |
---|
823 | |
---|
824 | |
---|
825 | } |
---|
826 | } |
---|
827 | |
---|
828 | if(shouldsplit) |
---|
829 | { |
---|
830 | current_parent->totallyneutralgrandchildren.insert(sender->totallyneutralchildren.begin(),sender->totallyneutralchildren.end()); |
---|
831 | |
---|
832 | for(set<polyhedron*>::iterator tot_child_ref = sender->totallyneutralchildren.begin();tot_child_ref!=sender->totallyneutralchildren.end();tot_child_ref++) |
---|
833 | { |
---|
834 | (*tot_child_ref)->grandparents.insert(current_parent); |
---|
835 | } |
---|
836 | |
---|
837 | switch(sender->get_state(SPLIT)) |
---|
838 | { |
---|
839 | case 1: |
---|
840 | current_parent->positivechildren.push_back(sender); |
---|
841 | current_parent->positiveneutralvertices.insert(sender->vertices.begin(),sender->vertices.end()); |
---|
842 | break; |
---|
843 | case 0: |
---|
844 | current_parent->neutralchildren.push_back(sender); |
---|
845 | current_parent->positiveneutralvertices.insert(sender->positiveneutralvertices.begin(),sender->positiveneutralvertices.end()); |
---|
846 | current_parent->negativeneutralvertices.insert(sender->negativeneutralvertices.begin(),sender->negativeneutralvertices.end()); |
---|
847 | |
---|
848 | if(current_parent->totally_neutral == NULL) |
---|
849 | { |
---|
850 | current_parent->totally_neutral = sender->totally_neutral; |
---|
851 | } |
---|
852 | else |
---|
853 | { |
---|
854 | current_parent->totally_neutral = current_parent->totally_neutral && sender->totally_neutral; |
---|
855 | } |
---|
856 | |
---|
857 | if(sender->totally_neutral) |
---|
858 | { |
---|
859 | current_parent->totallyneutralchildren.insert(sender); |
---|
860 | } |
---|
861 | |
---|
862 | break; |
---|
863 | case -1: |
---|
864 | current_parent->negativechildren.push_back(sender); |
---|
865 | current_parent->negativeneutralvertices.insert(sender->vertices.begin(),sender->vertices.end()); |
---|
866 | break; |
---|
867 | } |
---|
868 | |
---|
869 | if(is_last) |
---|
870 | { |
---|
871 | |
---|
872 | /// \TODO Nechapu druhou podminku, zda se mi ze je to spatne.. Nemela by byt jen prvni? Nebo se jedna o nastaveni totalni neutrality? |
---|
873 | if((current_parent->negativechildren.size()>0&¤t_parent->positivechildren.size()>0)|| |
---|
874 | (current_parent->neutralchildren.size()>0&¤t_parent->totally_neutral==false)) |
---|
875 | { |
---|
876 | for_splitting[level+1].push_back(current_parent); |
---|
877 | |
---|
878 | current_parent->set_state(0, SPLIT); |
---|
879 | } |
---|
880 | else |
---|
881 | { |
---|
882 | if(current_parent->negativechildren.size()>0) |
---|
883 | { |
---|
884 | current_parent->set_state(-1, SPLIT); |
---|
885 | |
---|
886 | ((toprow*)current_parent)->condition_sum-=toadd->value; |
---|
887 | |
---|
888 | } |
---|
889 | else if(current_parent->positivechildren.size()>0) |
---|
890 | { |
---|
891 | current_parent->set_state(1, SPLIT); |
---|
892 | |
---|
893 | ((toprow*)current_parent)->condition_sum+=toadd->value; |
---|
894 | } |
---|
895 | else |
---|
896 | { |
---|
897 | current_parent->raise_multiplicity(); |
---|
898 | } |
---|
899 | |
---|
900 | ((toprow*)current_parent)->condition_order++; |
---|
901 | |
---|
902 | if(level == number_of_parameters - 1) |
---|
903 | { |
---|
904 | toprow* cur_par_toprow = ((toprow*)current_parent); |
---|
905 | cur_par_toprow->probability = 0.0; |
---|
906 | |
---|
907 | //map<double,set<vertex*>> new_triangulation; |
---|
908 | |
---|
909 | for(set<simplex*>::iterator s_ref = current_parent->triangulation.begin();s_ref!=current_parent->triangulation.end();s_ref++) |
---|
910 | { |
---|
911 | double cur_prob = cur_par_toprow->integrate_simplex((*s_ref),'C'); |
---|
912 | |
---|
913 | cur_par_toprow->probability += cur_prob; |
---|
914 | |
---|
915 | //new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second)); |
---|
916 | } |
---|
917 | |
---|
918 | //current_parent->triangulation.clear(); |
---|
919 | //current_parent->triangulation.insert(new_triangulation.begin(),new_triangulation.end()); |
---|
920 | } |
---|
921 | |
---|
922 | if(current_parent->mergechild == NULL) |
---|
923 | { |
---|
924 | current_parent->positivechildren.clear(); |
---|
925 | current_parent->negativechildren.clear(); |
---|
926 | current_parent->neutralchildren.clear(); |
---|
927 | current_parent->totallyneutralchildren.clear(); |
---|
928 | current_parent->totallyneutralgrandchildren.clear(); |
---|
929 | // current_parent->grandparents.clear(); |
---|
930 | current_parent->positiveneutralvertices.clear(); |
---|
931 | current_parent->negativeneutralvertices.clear(); |
---|
932 | current_parent->totally_neutral = NULL; |
---|
933 | current_parent->kids_rel_addresses.clear(); |
---|
934 | } |
---|
935 | } |
---|
936 | } |
---|
937 | } |
---|
938 | |
---|
939 | if(is_last) |
---|
940 | { |
---|
941 | current_parent->mergechild = NULL; |
---|
942 | current_parent->message_counter = 0; |
---|
943 | |
---|
944 | send_state_message(current_parent,toadd,toremove,level+1); |
---|
945 | } |
---|
946 | |
---|
947 | } |
---|
948 | |
---|
949 | } |
---|
950 | } |
---|
951 | |
---|
952 | public: |
---|
953 | c_statistic statistic; |
---|
954 | |
---|
955 | vertex* minimal_vertex; |
---|
956 | |
---|
957 | double likelihood_value; |
---|
958 | |
---|
959 | vector<multiset<my_ivec>> correction_factors; |
---|
960 | |
---|
961 | int number_of_parameters; |
---|
962 | |
---|
963 | /// A default constructor creates an emlig with predefined statistic representing only the range of the given |
---|
964 | /// parametric space, where the number of parameters of the needed model is given as a parameter to the constructor. |
---|
965 | emlig(int number_of_parameters) |
---|
966 | { |
---|
967 | this->number_of_parameters = number_of_parameters; |
---|
968 | |
---|
969 | create_statistic(number_of_parameters); |
---|
970 | |
---|
971 | likelihood_value = numeric_limits<double>::max(); |
---|
972 | } |
---|
973 | |
---|
974 | /// A constructor for creating an emlig when the user wants to create the statistic by himself. The creation of a |
---|
975 | /// statistic is needed outside the constructor. Used for a user defined prior distribution on the parameters. |
---|
976 | emlig(c_statistic statistic) |
---|
977 | { |
---|
978 | this->statistic = statistic; |
---|
979 | |
---|
980 | likelihood_value = numeric_limits<double>::max(); |
---|
981 | } |
---|
982 | |
---|
983 | void step_me(int marker) |
---|
984 | { |
---|
985 | |
---|
986 | for(int i = 0;i<statistic.size();i++) |
---|
987 | { |
---|
988 | //int zero = 0; |
---|
989 | //int one = 0; |
---|
990 | //int two = 0; |
---|
991 | |
---|
992 | for(polyhedron* horiz_ref = statistic.rows[i];horiz_ref!=statistic.get_end();horiz_ref=horiz_ref->next_poly) |
---|
993 | { |
---|
994 | |
---|
995 | |
---|
996 | if(i==statistic.size()-1) |
---|
997 | { |
---|
998 | cout << ((toprow*)horiz_ref)->condition_sum << " " << ((toprow*)horiz_ref)->probability << endl; |
---|
999 | cout << "Order:" << ((toprow*)horiz_ref)->condition_order << endl; |
---|
1000 | } |
---|
1001 | |
---|
1002 | /* |
---|
1003 | if(i==0) |
---|
1004 | { |
---|
1005 | cout << ((vertex*)horiz_ref)->get_coordinates() << endl; |
---|
1006 | } |
---|
1007 | */ |
---|
1008 | |
---|
1009 | /* |
---|
1010 | char* string = "Checkpoint"; |
---|
1011 | |
---|
1012 | |
---|
1013 | if((*horiz_ref).parentconditions.size()==0) |
---|
1014 | { |
---|
1015 | zero++; |
---|
1016 | } |
---|
1017 | else if((*horiz_ref).parentconditions.size()==1) |
---|
1018 | { |
---|
1019 | one++; |
---|
1020 | } |
---|
1021 | else |
---|
1022 | { |
---|
1023 | two++; |
---|
1024 | } |
---|
1025 | */ |
---|
1026 | |
---|
1027 | } |
---|
1028 | } |
---|
1029 | |
---|
1030 | |
---|
1031 | /* |
---|
1032 | list<vec> table_entries; |
---|
1033 | for(polyhedron* horiz_ref = statistic.rows[statistic.size()-1];horiz_ref!=statistic.row_ends[statistic.size()-1];horiz_ref=horiz_ref->next_poly) |
---|
1034 | { |
---|
1035 | toprow *current_toprow = (toprow*)(horiz_ref); |
---|
1036 | for(list<set<vertex*>>::iterator tri_ref = current_toprow->triangulation.begin();tri_ref!=current_toprow->triangulation.end();tri_ref++) |
---|
1037 | { |
---|
1038 | for(set<vertex*>::iterator vert_ref = (*tri_ref).begin();vert_ref!=(*tri_ref).end();vert_ref++) |
---|
1039 | { |
---|
1040 | vec table_entry = vec(); |
---|
1041 | |
---|
1042 | table_entry.ins(0,(*vert_ref)->get_coordinates()*current_toprow->condition.get(1,current_toprow->condition.size()-1)-current_toprow->condition.get(0,0)); |
---|
1043 | |
---|
1044 | table_entry.ins(0,(*vert_ref)->get_coordinates()); |
---|
1045 | |
---|
1046 | table_entries.push_back(table_entry); |
---|
1047 | } |
---|
1048 | } |
---|
1049 | } |
---|
1050 | |
---|
1051 | unique(table_entries.begin(),table_entries.end()); |
---|
1052 | |
---|
1053 | |
---|
1054 | |
---|
1055 | for(list<vec>::iterator entry_ref = table_entries.begin();entry_ref!=table_entries.end();entry_ref++) |
---|
1056 | { |
---|
1057 | ofstream myfile; |
---|
1058 | myfile.open("robust_data.txt", ios::out | ios::app); |
---|
1059 | if (myfile.is_open()) |
---|
1060 | { |
---|
1061 | for(int i = 0;i<(*entry_ref).size();i++) |
---|
1062 | { |
---|
1063 | myfile << (*entry_ref)[i] << ";"; |
---|
1064 | } |
---|
1065 | myfile << endl; |
---|
1066 | |
---|
1067 | myfile.close(); |
---|
1068 | } |
---|
1069 | else |
---|
1070 | { |
---|
1071 | cout << "File problem." << endl; |
---|
1072 | } |
---|
1073 | } |
---|
1074 | */ |
---|
1075 | |
---|
1076 | |
---|
1077 | return; |
---|
1078 | } |
---|
1079 | |
---|
1080 | int statistic_rowsize(int row) |
---|
1081 | { |
---|
1082 | return statistic.row_size(row); |
---|
1083 | } |
---|
1084 | |
---|
1085 | void add_condition(vec toadd) |
---|
1086 | { |
---|
1087 | vec null_vector = ""; |
---|
1088 | |
---|
1089 | add_and_remove_condition(toadd, null_vector); |
---|
1090 | } |
---|
1091 | |
---|
1092 | |
---|
1093 | void remove_condition(vec toremove) |
---|
1094 | { |
---|
1095 | vec null_vector = ""; |
---|
1096 | |
---|
1097 | add_and_remove_condition(null_vector, toremove); |
---|
1098 | |
---|
1099 | } |
---|
1100 | |
---|
1101 | void add_and_remove_condition(vec toadd, vec toremove) |
---|
1102 | { |
---|
1103 | //step_me(0); |
---|
1104 | |
---|
1105 | likelihood_value = numeric_limits<double>::max(); |
---|
1106 | |
---|
1107 | bool should_remove = (toremove.size() != 0); |
---|
1108 | bool should_add = (toadd.size() != 0); |
---|
1109 | |
---|
1110 | for_splitting.clear(); |
---|
1111 | for_merging.clear(); |
---|
1112 | |
---|
1113 | for(int i = 0;i<statistic.size();i++) |
---|
1114 | { |
---|
1115 | list<polyhedron*> empty_split; |
---|
1116 | list<polyhedron*> empty_merge; |
---|
1117 | |
---|
1118 | for_splitting.push_back(empty_split); |
---|
1119 | for_merging.push_back(empty_merge); |
---|
1120 | } |
---|
1121 | |
---|
1122 | list<condition*>::iterator toremove_ref = conditions.end(); |
---|
1123 | bool condition_should_be_added = should_add; |
---|
1124 | |
---|
1125 | for(list<condition*>::iterator ref = conditions.begin();ref!=conditions.end();ref++) |
---|
1126 | { |
---|
1127 | if(should_remove) |
---|
1128 | { |
---|
1129 | if((*ref)->value == toremove) |
---|
1130 | { |
---|
1131 | if((*ref)->multiplicity>1) |
---|
1132 | { |
---|
1133 | (*ref)->multiplicity--; |
---|
1134 | |
---|
1135 | alter_toprow_conditions(*ref,false); |
---|
1136 | |
---|
1137 | should_remove = false; |
---|
1138 | } |
---|
1139 | else |
---|
1140 | { |
---|
1141 | toremove_ref = ref; |
---|
1142 | } |
---|
1143 | } |
---|
1144 | } |
---|
1145 | |
---|
1146 | if(should_add) |
---|
1147 | { |
---|
1148 | if((*ref)->value == toadd) |
---|
1149 | { |
---|
1150 | (*ref)->multiplicity++; |
---|
1151 | |
---|
1152 | alter_toprow_conditions(*ref,true); |
---|
1153 | |
---|
1154 | should_add = false; |
---|
1155 | |
---|
1156 | condition_should_be_added = false; |
---|
1157 | } |
---|
1158 | } |
---|
1159 | } |
---|
1160 | |
---|
1161 | condition* condition_to_remove = NULL; |
---|
1162 | |
---|
1163 | if(toremove_ref!=conditions.end()) |
---|
1164 | { |
---|
1165 | condition_to_remove = *toremove_ref; |
---|
1166 | conditions.erase(toremove_ref); |
---|
1167 | } |
---|
1168 | |
---|
1169 | condition* condition_to_add = NULL; |
---|
1170 | |
---|
1171 | if(condition_should_be_added) |
---|
1172 | { |
---|
1173 | condition* new_condition = new condition(toadd); |
---|
1174 | |
---|
1175 | conditions.push_back(new_condition); |
---|
1176 | condition_to_add = new_condition; |
---|
1177 | } |
---|
1178 | |
---|
1179 | for(polyhedron* horizontal_position = statistic.rows[0];horizontal_position!=statistic.get_end();horizontal_position=horizontal_position->next_poly) |
---|
1180 | { |
---|
1181 | vertex* current_vertex = (vertex*)horizontal_position; |
---|
1182 | |
---|
1183 | if(should_add||should_remove) |
---|
1184 | { |
---|
1185 | vec appended_coords = current_vertex->get_coordinates(); |
---|
1186 | appended_coords.ins(0,-1.0); |
---|
1187 | |
---|
1188 | if(should_add) |
---|
1189 | { |
---|
1190 | double local_condition = 0;// = toadd*(appended_coords.first/=appended_coords.second); |
---|
1191 | |
---|
1192 | local_condition = appended_coords*toadd; |
---|
1193 | |
---|
1194 | current_vertex->set_state(local_condition,SPLIT); |
---|
1195 | |
---|
1196 | /// \TODO There should be a rounding error tolerance used here to insure we are not having too many points because of rounding error. |
---|
1197 | if(local_condition == 0) |
---|
1198 | { |
---|
1199 | current_vertex->totally_neutral = true; |
---|
1200 | |
---|
1201 | current_vertex->raise_multiplicity(); |
---|
1202 | |
---|
1203 | current_vertex->negativeneutralvertices.insert(current_vertex); |
---|
1204 | current_vertex->positiveneutralvertices.insert(current_vertex); |
---|
1205 | } |
---|
1206 | } |
---|
1207 | |
---|
1208 | if(should_remove) |
---|
1209 | { |
---|
1210 | set<condition*>::iterator cond_ref; |
---|
1211 | |
---|
1212 | for(cond_ref = current_vertex->parentconditions.begin();cond_ref!=current_vertex->parentconditions.end();cond_ref++) |
---|
1213 | { |
---|
1214 | if(*cond_ref == condition_to_remove) |
---|
1215 | { |
---|
1216 | break; |
---|
1217 | } |
---|
1218 | } |
---|
1219 | |
---|
1220 | if(cond_ref!=current_vertex->parentconditions.end()) |
---|
1221 | { |
---|
1222 | current_vertex->parentconditions.erase(cond_ref); |
---|
1223 | current_vertex->set_state(0,MERGE); |
---|
1224 | for_merging[0].push_back(current_vertex); |
---|
1225 | } |
---|
1226 | else |
---|
1227 | { |
---|
1228 | double local_condition = toremove*appended_coords; |
---|
1229 | current_vertex->set_state(local_condition,MERGE); |
---|
1230 | } |
---|
1231 | } |
---|
1232 | } |
---|
1233 | |
---|
1234 | send_state_message(current_vertex, condition_to_add, condition_to_remove, 0); |
---|
1235 | |
---|
1236 | } |
---|
1237 | |
---|
1238 | |
---|
1239 | |
---|
1240 | if(should_remove) |
---|
1241 | { |
---|
1242 | for(int i = 0;i<for_merging.size();i++) |
---|
1243 | { |
---|
1244 | for(list<polyhedron*>::iterator merge_ref = for_merging[i].begin();merge_ref!=for_merging[i].end();merge_ref++) |
---|
1245 | { |
---|
1246 | cout << (*merge_ref)->get_state(MERGE) << ","; |
---|
1247 | } |
---|
1248 | |
---|
1249 | cout << endl; |
---|
1250 | } |
---|
1251 | |
---|
1252 | set<vertex*> vertices_to_be_reduced; |
---|
1253 | |
---|
1254 | int k = 1; |
---|
1255 | |
---|
1256 | for(vector<list<polyhedron*>>::iterator vert_ref = for_merging.begin();vert_ref<for_merging.end();vert_ref++) |
---|
1257 | { |
---|
1258 | for(list<polyhedron*>::reverse_iterator merge_ref = vert_ref->rbegin();merge_ref!=vert_ref->rend();merge_ref++) |
---|
1259 | { |
---|
1260 | if((*merge_ref)->get_multiplicity()>1) |
---|
1261 | { |
---|
1262 | if(k==1) |
---|
1263 | { |
---|
1264 | vertices_to_be_reduced.insert((vertex*)(*merge_ref)); |
---|
1265 | } |
---|
1266 | else |
---|
1267 | { |
---|
1268 | (*merge_ref)->lower_multiplicity(); |
---|
1269 | } |
---|
1270 | } |
---|
1271 | else |
---|
1272 | { |
---|
1273 | toprow* current_positive = (toprow*)(*merge_ref)->positiveparent; |
---|
1274 | toprow* current_negative = (toprow*)(*merge_ref)->negativeparent; |
---|
1275 | |
---|
1276 | //current_positive->condition_sum -= toremove; |
---|
1277 | //current_positive->condition_order--; |
---|
1278 | |
---|
1279 | current_positive->parentconditions.erase(condition_to_remove); |
---|
1280 | |
---|
1281 | current_positive->children.insert(current_positive->children.end(),current_negative->children.begin(),current_negative->children.end()); |
---|
1282 | current_positive->children.remove(*merge_ref); |
---|
1283 | |
---|
1284 | for(list<polyhedron*>::iterator child_ref = current_negative->children.begin();child_ref!=current_negative->children.end();child_ref++) |
---|
1285 | { |
---|
1286 | (*child_ref)->parents.remove(current_negative); |
---|
1287 | (*child_ref)->parents.push_back(current_positive); |
---|
1288 | } |
---|
1289 | |
---|
1290 | // current_positive->parents.insert(current_positive->parents.begin(),current_negative->parents.begin(),current_negative->parents.end()); |
---|
1291 | // unique(current_positive->parents.begin(),current_positive->parents.end()); |
---|
1292 | |
---|
1293 | for(list<polyhedron*>::iterator parent_ref = current_negative->parents.begin();parent_ref!=current_negative->parents.end();parent_ref++) |
---|
1294 | { |
---|
1295 | (*parent_ref)->children.remove(current_negative); |
---|
1296 | |
---|
1297 | switch(current_negative->get_state(SPLIT)) |
---|
1298 | { |
---|
1299 | case -1: |
---|
1300 | (*parent_ref)->negativechildren.remove(current_negative); |
---|
1301 | break; |
---|
1302 | case 0: |
---|
1303 | (*parent_ref)->neutralchildren.remove(current_negative); |
---|
1304 | break; |
---|
1305 | case 1: |
---|
1306 | (*parent_ref)->positivechildren.remove(current_negative); |
---|
1307 | break; |
---|
1308 | } |
---|
1309 | //(*parent_ref)->children.push_back(current_positive); |
---|
1310 | } |
---|
1311 | |
---|
1312 | if(current_positive->get_state(SPLIT)!=0&¤t_negative->get_state(SPLIT)==0) |
---|
1313 | { |
---|
1314 | for(list<polyhedron*>::iterator parent_ref = current_positive->parents.begin();parent_ref!=current_positive->parents.end();parent_ref++) |
---|
1315 | { |
---|
1316 | if(current_positive->get_state(SPLIT)==1) |
---|
1317 | { |
---|
1318 | (*parent_ref)->positivechildren.remove(current_positive); |
---|
1319 | } |
---|
1320 | else |
---|
1321 | { |
---|
1322 | (*parent_ref)->negativechildren.remove(current_positive); |
---|
1323 | } |
---|
1324 | |
---|
1325 | (*parent_ref)->neutralchildren.push_back(current_positive); |
---|
1326 | } |
---|
1327 | |
---|
1328 | current_positive->set_state(0,SPLIT); |
---|
1329 | for_splitting[k].push_back(current_positive); |
---|
1330 | } |
---|
1331 | |
---|
1332 | if((current_positive->get_state(SPLIT)==0&&!current_positive->totally_neutral)||(current_negative->get_state(SPLIT)==0&&!current_negative->totally_neutral)) |
---|
1333 | { |
---|
1334 | current_positive->negativechildren.insert(current_positive->negativechildren.end(),current_negative->negativechildren.begin(),current_negative->negativechildren.end()); |
---|
1335 | |
---|
1336 | current_positive->positivechildren.insert(current_positive->positivechildren.end(),current_negative->positivechildren.begin(),current_negative->positivechildren.end()); |
---|
1337 | |
---|
1338 | current_positive->neutralchildren.insert(current_positive->neutralchildren.end(),current_negative->neutralchildren.begin(),current_negative->neutralchildren.end()); |
---|
1339 | |
---|
1340 | switch((*merge_ref)->get_state(SPLIT)) |
---|
1341 | { |
---|
1342 | case -1: |
---|
1343 | current_positive->negativechildren.remove(*merge_ref); |
---|
1344 | break; |
---|
1345 | case 0: |
---|
1346 | current_positive->neutralchildren.remove(*merge_ref); |
---|
1347 | break; |
---|
1348 | case 1: |
---|
1349 | current_positive->positivechildren.remove(*merge_ref); |
---|
1350 | break; |
---|
1351 | } |
---|
1352 | |
---|
1353 | current_positive->totallyneutralchildren.insert(current_negative->totallyneutralchildren.begin(),current_negative->totallyneutralchildren.end()); |
---|
1354 | |
---|
1355 | current_positive->totallyneutralchildren.erase(*merge_ref); |
---|
1356 | |
---|
1357 | current_positive->totallyneutralgrandchildren.insert(current_negative->totallyneutralgrandchildren.begin(),current_negative->totallyneutralgrandchildren.end()); |
---|
1358 | |
---|
1359 | current_positive->negativeneutralvertices.insert(current_negative->negativeneutralvertices.begin(),current_negative->negativeneutralvertices.end()); |
---|
1360 | current_positive->positiveneutralvertices.insert(current_negative->positiveneutralvertices.begin(),current_negative->positiveneutralvertices.end()); |
---|
1361 | } |
---|
1362 | else |
---|
1363 | { |
---|
1364 | if(!current_positive->totally_neutral) |
---|
1365 | { |
---|
1366 | current_positive->positivechildren.clear(); |
---|
1367 | current_positive->negativechildren.clear(); |
---|
1368 | current_positive->neutralchildren.clear(); |
---|
1369 | current_positive->totallyneutralchildren.clear(); |
---|
1370 | current_positive->totallyneutralgrandchildren.clear(); |
---|
1371 | current_positive->positiveneutralvertices.clear(); |
---|
1372 | current_positive->negativeneutralvertices.clear(); |
---|
1373 | current_positive->totally_neutral = NULL; |
---|
1374 | current_positive->kids_rel_addresses.clear(); |
---|
1375 | } |
---|
1376 | |
---|
1377 | } |
---|
1378 | |
---|
1379 | |
---|
1380 | |
---|
1381 | current_positive->vertices.insert(current_negative->vertices.begin(),current_negative->vertices.end()); |
---|
1382 | |
---|
1383 | |
---|
1384 | for(set<vertex*>::iterator vert_ref = (*merge_ref)->vertices.begin();vert_ref!=(*merge_ref)->vertices.end();vert_ref++) |
---|
1385 | { |
---|
1386 | if((*vert_ref)->get_multiplicity()==1) |
---|
1387 | { |
---|
1388 | current_positive->vertices.erase(*vert_ref); |
---|
1389 | |
---|
1390 | if((current_positive->get_state(SPLIT)==0&&!current_positive->totally_neutral)||(current_negative->get_state(SPLIT)==0&&!current_negative->totally_neutral)) |
---|
1391 | { |
---|
1392 | current_positive->negativeneutralvertices.erase(*vert_ref); |
---|
1393 | current_positive->positiveneutralvertices.erase(*vert_ref); |
---|
1394 | } |
---|
1395 | } |
---|
1396 | } |
---|
1397 | |
---|
1398 | if(current_negative->get_state(SPLIT)==0&&!current_negative->totally_neutral) |
---|
1399 | { |
---|
1400 | for_splitting[k].remove(current_negative); |
---|
1401 | } |
---|
1402 | |
---|
1403 | if(current_positive->totally_neutral) |
---|
1404 | { |
---|
1405 | if(!current_negative->totally_neutral) |
---|
1406 | { |
---|
1407 | for(set<polyhedron*>::iterator grand_ref = current_positive->grandparents.begin();grand_ref!=current_positive->grandparents.end();grand_ref++) |
---|
1408 | { |
---|
1409 | (*grand_ref)->totallyneutralgrandchildren.erase(current_positive); |
---|
1410 | } |
---|
1411 | } |
---|
1412 | else |
---|
1413 | { |
---|
1414 | for(set<polyhedron*>::iterator grand_ref = current_negative->grandparents.begin();grand_ref!=current_negative->grandparents.end();grand_ref++) |
---|
1415 | { |
---|
1416 | (*grand_ref)->totallyneutralgrandchildren.erase(current_negative); |
---|
1417 | (*grand_ref)->totallyneutralgrandchildren.insert(current_positive); |
---|
1418 | } |
---|
1419 | } |
---|
1420 | } |
---|
1421 | else |
---|
1422 | { |
---|
1423 | if(current_negative->totally_neutral) |
---|
1424 | { |
---|
1425 | for(set<polyhedron*>::iterator grand_ref = current_negative->grandparents.begin();grand_ref!=current_negative->grandparents.end();grand_ref++) |
---|
1426 | { |
---|
1427 | (*grand_ref)->totallyneutralgrandchildren.erase(current_negative); |
---|
1428 | } |
---|
1429 | } |
---|
1430 | } |
---|
1431 | |
---|
1432 | current_positive->grandparents.clear(); |
---|
1433 | |
---|
1434 | |
---|
1435 | |
---|
1436 | current_positive->totally_neutral = (current_positive->totally_neutral && current_negative->totally_neutral); |
---|
1437 | |
---|
1438 | current_positive->triangulate(k==for_splitting.size()-1); |
---|
1439 | |
---|
1440 | statistic.delete_polyhedron(k,current_negative); |
---|
1441 | |
---|
1442 | delete current_negative; |
---|
1443 | |
---|
1444 | for(list<polyhedron*>::iterator child_ref = (*merge_ref)->children.begin();child_ref!=(*merge_ref)->children.end();child_ref++) |
---|
1445 | { |
---|
1446 | (*child_ref)->parents.remove(*merge_ref); |
---|
1447 | } |
---|
1448 | |
---|
1449 | /* |
---|
1450 | for(list<polyhedron*>::iterator parent_ref = (*merge_ref)->parents.begin();parent_ref!=(*merge_ref)->parents.end();parent_ref++) |
---|
1451 | { |
---|
1452 | (*parent_ref)->positivechildren.remove(*merge_ref); |
---|
1453 | (*parent_ref)->negativechildren.remove(*merge_ref); |
---|
1454 | (*parent_ref)->neutralchildren.remove(*merge_ref); |
---|
1455 | (*parent_ref)->children.remove(*merge_ref); |
---|
1456 | } |
---|
1457 | */ |
---|
1458 | |
---|
1459 | for(set<polyhedron*>::iterator grand_ch_ref = (*merge_ref)->totallyneutralgrandchildren.begin();grand_ch_ref!=(*merge_ref)->totallyneutralgrandchildren.end();grand_ch_ref++) |
---|
1460 | { |
---|
1461 | (*grand_ch_ref)->grandparents.erase(*merge_ref); |
---|
1462 | } |
---|
1463 | |
---|
1464 | |
---|
1465 | for(set<polyhedron*>::iterator grand_p_ref = (*merge_ref)->grandparents.begin();grand_p_ref!=(*merge_ref)->grandparents.end();grand_p_ref++) |
---|
1466 | { |
---|
1467 | (*grand_p_ref)->totallyneutralgrandchildren.erase(*merge_ref); |
---|
1468 | } |
---|
1469 | |
---|
1470 | for_splitting[k-1].remove(*merge_ref); |
---|
1471 | |
---|
1472 | statistic.delete_polyhedron(k-1,*merge_ref); |
---|
1473 | |
---|
1474 | if(k==1) |
---|
1475 | { |
---|
1476 | vertices_to_be_reduced.insert((vertex*)(*merge_ref)); |
---|
1477 | } |
---|
1478 | else |
---|
1479 | { |
---|
1480 | delete *merge_ref; |
---|
1481 | } |
---|
1482 | } |
---|
1483 | } |
---|
1484 | |
---|
1485 | k++; |
---|
1486 | |
---|
1487 | } |
---|
1488 | |
---|
1489 | for(set<vertex*>::iterator vert_ref = vertices_to_be_reduced.begin();vert_ref!=vertices_to_be_reduced.end();vert_ref++) |
---|
1490 | { |
---|
1491 | if((*vert_ref)->get_multiplicity()>1) |
---|
1492 | { |
---|
1493 | (*vert_ref)->lower_multiplicity(); |
---|
1494 | } |
---|
1495 | else |
---|
1496 | { |
---|
1497 | delete *vert_ref; |
---|
1498 | } |
---|
1499 | } |
---|
1500 | |
---|
1501 | delete condition_to_remove; |
---|
1502 | } |
---|
1503 | |
---|
1504 | vector<int> sizevector; |
---|
1505 | for(int s = 0;s<statistic.size();s++) |
---|
1506 | { |
---|
1507 | sizevector.push_back(statistic.row_size(s)); |
---|
1508 | cout << statistic.row_size(s) << ", "; |
---|
1509 | } |
---|
1510 | |
---|
1511 | cout << endl; |
---|
1512 | |
---|
1513 | if(should_add) |
---|
1514 | { |
---|
1515 | int k = 1; |
---|
1516 | |
---|
1517 | vector<list<polyhedron*>>::iterator beginning_ref = ++for_splitting.begin(); |
---|
1518 | |
---|
1519 | for(vector<list<polyhedron*>>::iterator vert_ref = beginning_ref;vert_ref<for_splitting.end();vert_ref++) |
---|
1520 | { |
---|
1521 | |
---|
1522 | for(list<polyhedron*>::reverse_iterator split_ref = vert_ref->rbegin();split_ref != vert_ref->rend();split_ref++) |
---|
1523 | { |
---|
1524 | polyhedron* new_totally_neutral_child; |
---|
1525 | |
---|
1526 | polyhedron* current_polyhedron = (*split_ref); |
---|
1527 | |
---|
1528 | if(vert_ref == beginning_ref) |
---|
1529 | { |
---|
1530 | vec coordinates1 = ((vertex*)(*(current_polyhedron->children.begin())))->get_coordinates(); |
---|
1531 | vec coordinates2 = ((vertex*)(*(++current_polyhedron->children.begin())))->get_coordinates(); |
---|
1532 | |
---|
1533 | vec extended_coord2 = coordinates2; |
---|
1534 | extended_coord2.ins(0,-1.0); |
---|
1535 | |
---|
1536 | double t = (-toadd*extended_coord2)/(toadd(1,toadd.size()-1)*(coordinates1-coordinates2)); |
---|
1537 | |
---|
1538 | vec new_coordinates = (1-t)*coordinates2+t*coordinates1; |
---|
1539 | |
---|
1540 | // cout << "c1:" << coordinates1 << endl << "c2:" << coordinates2 << endl << "nc:" << new_coordinates << endl; |
---|
1541 | |
---|
1542 | vertex* neutral_vertex = new vertex(new_coordinates); |
---|
1543 | |
---|
1544 | new_totally_neutral_child = neutral_vertex; |
---|
1545 | } |
---|
1546 | else |
---|
1547 | { |
---|
1548 | toprow* neutral_toprow = new toprow(); |
---|
1549 | |
---|
1550 | neutral_toprow->condition_sum = ((toprow*)current_polyhedron)->condition_sum; // tohle tu bylo driv: zeros(number_of_parameters+1); |
---|
1551 | neutral_toprow->condition_order = ((toprow*)current_polyhedron)->condition_order+1; |
---|
1552 | |
---|
1553 | new_totally_neutral_child = neutral_toprow; |
---|
1554 | } |
---|
1555 | |
---|
1556 | new_totally_neutral_child->parentconditions.insert(current_polyhedron->parentconditions.begin(),current_polyhedron->parentconditions.end()); |
---|
1557 | new_totally_neutral_child->parentconditions.insert(condition_to_add); |
---|
1558 | |
---|
1559 | new_totally_neutral_child->my_emlig = this; |
---|
1560 | |
---|
1561 | new_totally_neutral_child->children.insert(new_totally_neutral_child->children.end(), |
---|
1562 | current_polyhedron->totallyneutralgrandchildren.begin(), |
---|
1563 | current_polyhedron->totallyneutralgrandchildren.end()); |
---|
1564 | |
---|
1565 | |
---|
1566 | |
---|
1567 | // cout << ((toprow*)current_polyhedron)->condition << endl << toadd << endl; |
---|
1568 | |
---|
1569 | toprow* positive_poly = new toprow(((toprow*)current_polyhedron)->condition_sum+toadd, ((toprow*)current_polyhedron)->condition_order+1); |
---|
1570 | toprow* negative_poly = new toprow(((toprow*)current_polyhedron)->condition_sum-toadd, ((toprow*)current_polyhedron)->condition_order+1); |
---|
1571 | |
---|
1572 | positive_poly->my_emlig = this; |
---|
1573 | negative_poly->my_emlig = this; |
---|
1574 | |
---|
1575 | positive_poly->parentconditions.insert(current_polyhedron->parentconditions.begin(),current_polyhedron->parentconditions.end()); |
---|
1576 | negative_poly->parentconditions.insert(current_polyhedron->parentconditions.begin(),current_polyhedron->parentconditions.end()); |
---|
1577 | |
---|
1578 | for(set<polyhedron*>::iterator grand_ref = current_polyhedron->totallyneutralgrandchildren.begin(); grand_ref != current_polyhedron->totallyneutralgrandchildren.end();grand_ref++) |
---|
1579 | { |
---|
1580 | (*grand_ref)->parents.push_back(new_totally_neutral_child); |
---|
1581 | |
---|
1582 | // tohle tu nebylo. ma to tu byt? |
---|
1583 | //positive_poly->totallyneutralgrandchildren.insert(*grand_ref); |
---|
1584 | //negative_poly->totallyneutralgrandchildren.insert(*grand_ref); |
---|
1585 | |
---|
1586 | //(*grand_ref)->grandparents.insert(positive_poly); |
---|
1587 | //(*grand_ref)->grandparents.insert(negative_poly); |
---|
1588 | |
---|
1589 | new_totally_neutral_child->vertices.insert((*grand_ref)->vertices.begin(),(*grand_ref)->vertices.end()); |
---|
1590 | } |
---|
1591 | |
---|
1592 | positive_poly->children.push_back(new_totally_neutral_child); |
---|
1593 | negative_poly->children.push_back(new_totally_neutral_child); |
---|
1594 | |
---|
1595 | |
---|
1596 | for(list<polyhedron*>::iterator parent_ref = current_polyhedron->parents.begin();parent_ref!=current_polyhedron->parents.end();parent_ref++) |
---|
1597 | { |
---|
1598 | (*parent_ref)->totallyneutralgrandchildren.insert(new_totally_neutral_child); |
---|
1599 | // new_totally_neutral_child->grandparents.insert(*parent_ref); |
---|
1600 | |
---|
1601 | (*parent_ref)->neutralchildren.remove(current_polyhedron); |
---|
1602 | (*parent_ref)->children.remove(current_polyhedron); |
---|
1603 | |
---|
1604 | (*parent_ref)->children.push_back(positive_poly); |
---|
1605 | (*parent_ref)->children.push_back(negative_poly); |
---|
1606 | (*parent_ref)->positivechildren.push_back(positive_poly); |
---|
1607 | (*parent_ref)->negativechildren.push_back(negative_poly); |
---|
1608 | } |
---|
1609 | |
---|
1610 | positive_poly->parents.insert(positive_poly->parents.end(), |
---|
1611 | current_polyhedron->parents.begin(), |
---|
1612 | current_polyhedron->parents.end()); |
---|
1613 | |
---|
1614 | negative_poly->parents.insert(negative_poly->parents.end(), |
---|
1615 | current_polyhedron->parents.begin(), |
---|
1616 | current_polyhedron->parents.end()); |
---|
1617 | |
---|
1618 | |
---|
1619 | |
---|
1620 | new_totally_neutral_child->parents.push_back(positive_poly); |
---|
1621 | new_totally_neutral_child->parents.push_back(negative_poly); |
---|
1622 | |
---|
1623 | for(list<polyhedron*>::iterator child_ref = current_polyhedron->positivechildren.begin();child_ref!=current_polyhedron->positivechildren.end();child_ref++) |
---|
1624 | { |
---|
1625 | (*child_ref)->parents.remove(current_polyhedron); |
---|
1626 | (*child_ref)->parents.push_back(positive_poly); |
---|
1627 | } |
---|
1628 | |
---|
1629 | positive_poly->children.insert(positive_poly->children.end(), |
---|
1630 | current_polyhedron->positivechildren.begin(), |
---|
1631 | current_polyhedron->positivechildren.end()); |
---|
1632 | |
---|
1633 | for(list<polyhedron*>::iterator child_ref = current_polyhedron->negativechildren.begin();child_ref!=current_polyhedron->negativechildren.end();child_ref++) |
---|
1634 | { |
---|
1635 | (*child_ref)->parents.remove(current_polyhedron); |
---|
1636 | (*child_ref)->parents.push_back(negative_poly); |
---|
1637 | } |
---|
1638 | |
---|
1639 | negative_poly->children.insert(negative_poly->children.end(), |
---|
1640 | current_polyhedron->negativechildren.begin(), |
---|
1641 | current_polyhedron->negativechildren.end()); |
---|
1642 | |
---|
1643 | positive_poly->vertices.insert(current_polyhedron->positiveneutralvertices.begin(),current_polyhedron->positiveneutralvertices.end()); |
---|
1644 | positive_poly->vertices.insert(new_totally_neutral_child->vertices.begin(),new_totally_neutral_child->vertices.end()); |
---|
1645 | |
---|
1646 | negative_poly->vertices.insert(current_polyhedron->negativeneutralvertices.begin(),current_polyhedron->negativeneutralvertices.end()); |
---|
1647 | negative_poly->vertices.insert(new_totally_neutral_child->vertices.begin(),new_totally_neutral_child->vertices.end()); |
---|
1648 | |
---|
1649 | new_totally_neutral_child->triangulate(false); |
---|
1650 | |
---|
1651 | positive_poly->triangulate(k==for_splitting.size()-1); |
---|
1652 | negative_poly->triangulate(k==for_splitting.size()-1); |
---|
1653 | |
---|
1654 | statistic.append_polyhedron(k-1, new_totally_neutral_child); |
---|
1655 | |
---|
1656 | statistic.insert_polyhedron(k, positive_poly, current_polyhedron); |
---|
1657 | statistic.insert_polyhedron(k, negative_poly, current_polyhedron); |
---|
1658 | |
---|
1659 | statistic.delete_polyhedron(k, current_polyhedron); |
---|
1660 | |
---|
1661 | delete current_polyhedron; |
---|
1662 | } |
---|
1663 | |
---|
1664 | k++; |
---|
1665 | } |
---|
1666 | } |
---|
1667 | |
---|
1668 | |
---|
1669 | sizevector.clear(); |
---|
1670 | for(int s = 0;s<statistic.size();s++) |
---|
1671 | { |
---|
1672 | sizevector.push_back(statistic.row_size(s)); |
---|
1673 | cout << statistic.row_size(s) << ", "; |
---|
1674 | } |
---|
1675 | |
---|
1676 | cout << endl; |
---|
1677 | |
---|
1678 | /* |
---|
1679 | 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) |
---|
1680 | { |
---|
1681 | cout << ((toprow*)topr_ref)->condition << endl; |
---|
1682 | } |
---|
1683 | */ |
---|
1684 | |
---|
1685 | // step_me(0); |
---|
1686 | |
---|
1687 | } |
---|
1688 | |
---|
1689 | void set_correction_factors(int order) |
---|
1690 | { |
---|
1691 | for(int remaining_order = correction_factors.size();remaining_order<order;remaining_order++) |
---|
1692 | { |
---|
1693 | multiset<my_ivec> factor_templates; |
---|
1694 | multiset<my_ivec> final_factors; |
---|
1695 | |
---|
1696 | my_ivec orig_template = my_ivec(); |
---|
1697 | |
---|
1698 | for(int i = 1;i<number_of_parameters-remaining_order+1;i++) |
---|
1699 | { |
---|
1700 | bool in_cycle = false; |
---|
1701 | for(int j = 0;j<=remaining_order;j++) { |
---|
1702 | |
---|
1703 | multiset<my_ivec>::iterator fac_ref = factor_templates.begin(); |
---|
1704 | |
---|
1705 | do |
---|
1706 | { |
---|
1707 | my_ivec current_template; |
---|
1708 | if(!in_cycle) |
---|
1709 | { |
---|
1710 | current_template = orig_template; |
---|
1711 | in_cycle = true; |
---|
1712 | } |
---|
1713 | else |
---|
1714 | { |
---|
1715 | current_template = (*fac_ref); |
---|
1716 | fac_ref++; |
---|
1717 | } |
---|
1718 | |
---|
1719 | current_template.ins(current_template.size(),i); |
---|
1720 | |
---|
1721 | // cout << "template:" << current_template << endl; |
---|
1722 | |
---|
1723 | if(current_template.size()==remaining_order+1) |
---|
1724 | { |
---|
1725 | final_factors.insert(current_template); |
---|
1726 | } |
---|
1727 | else |
---|
1728 | { |
---|
1729 | factor_templates.insert(current_template); |
---|
1730 | } |
---|
1731 | } |
---|
1732 | while(fac_ref!=factor_templates.end()); |
---|
1733 | } |
---|
1734 | } |
---|
1735 | |
---|
1736 | correction_factors.push_back(final_factors); |
---|
1737 | |
---|
1738 | } |
---|
1739 | } |
---|
1740 | |
---|
1741 | |
---|
1742 | |
---|
1743 | mat sample_mat(int n) |
---|
1744 | { |
---|
1745 | |
---|
1746 | |
---|
1747 | |
---|
1748 | /// \TODO tady je to spatne, tady nesmi byt conditions.size(), viz RARX.bayes() |
---|
1749 | if(conditions.size()-2-number_of_parameters>=0) |
---|
1750 | { |
---|
1751 | mat sample_mat; |
---|
1752 | map<double,toprow*> ordered_toprows; |
---|
1753 | double sum_a = 0; |
---|
1754 | |
---|
1755 | for(polyhedron* top_ref = statistic.rows[number_of_parameters];top_ref!=statistic.end_poly;top_ref=top_ref->next_poly) |
---|
1756 | { |
---|
1757 | toprow* current_top = (toprow*)top_ref; |
---|
1758 | |
---|
1759 | sum_a+=current_top->probability; |
---|
1760 | ordered_toprows.insert(pair<double,toprow*>(sum_a,current_top)); |
---|
1761 | } |
---|
1762 | |
---|
1763 | while(sample_mat.cols()<n) |
---|
1764 | { |
---|
1765 | double rnumber = randu()*sum_a; |
---|
1766 | |
---|
1767 | // cout << "RND:" << rnumber << endl; |
---|
1768 | |
---|
1769 | // This could be more efficient (log n), but map::upper_bound() doesn't let me dereference returned iterator |
---|
1770 | toprow* sampled_toprow; |
---|
1771 | for(map<double,toprow*>::iterator top_ref = ordered_toprows.begin();top_ref!=ordered_toprows.end();top_ref++) |
---|
1772 | { |
---|
1773 | // cout << "CDF:"<< (*top_ref).first << endl; |
---|
1774 | |
---|
1775 | if((*top_ref).first >= rnumber) |
---|
1776 | { |
---|
1777 | sampled_toprow = (*top_ref).second; |
---|
1778 | break; |
---|
1779 | } |
---|
1780 | } |
---|
1781 | |
---|
1782 | // cout << &sampled_toprow << ";"; |
---|
1783 | |
---|
1784 | rnumber = randu(); |
---|
1785 | |
---|
1786 | set<simplex*>::iterator s_ref; |
---|
1787 | double sum_b = 0; |
---|
1788 | for(s_ref = sampled_toprow->triangulation.begin();s_ref!=sampled_toprow->triangulation.end();s_ref++) |
---|
1789 | { |
---|
1790 | sum_b += (*s_ref)->probability; |
---|
1791 | |
---|
1792 | if(sum_b/sampled_toprow->probability >= rnumber) |
---|
1793 | break; |
---|
1794 | } |
---|
1795 | |
---|
1796 | // cout << &(*tri_ref) << endl; |
---|
1797 | |
---|
1798 | rnumber = randu(); |
---|
1799 | |
---|
1800 | double sigma = 0; |
---|
1801 | double sum_g = 0; |
---|
1802 | for(int i = 0;i<(*s_ref)->gamma_parameters.size();i++) |
---|
1803 | { |
---|
1804 | for(multimap<double,double>::iterator g_ref = (*s_ref)->gamma_parameters[i].begin();g_ref != (*s_ref)->gamma_parameters[i].end();g_ref++) |
---|
1805 | { |
---|
1806 | sum_g += (*g_ref).first/(*s_ref)->gamma_sum; |
---|
1807 | |
---|
1808 | if(sum_g>rnumber) |
---|
1809 | { |
---|
1810 | itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,(*g_ref).second); |
---|
1811 | sigma = (*gamma)(); |
---|
1812 | break; |
---|
1813 | } |
---|
1814 | } |
---|
1815 | |
---|
1816 | if(sigma!=0) |
---|
1817 | { |
---|
1818 | break; |
---|
1819 | } |
---|
1820 | } |
---|
1821 | |
---|
1822 | int dimension = (*s_ref)->vertices.size()-1; |
---|
1823 | |
---|
1824 | mat jacobian(dimension,dimension); |
---|
1825 | vec gradient = sampled_toprow->condition_sum.right(dimension); |
---|
1826 | |
---|
1827 | vertex* base_vert = *(*s_ref)->vertices.begin(); |
---|
1828 | |
---|
1829 | int row_count = 0; |
---|
1830 | |
---|
1831 | for(set<vertex*>::iterator vert_ref = ++(*s_ref)->vertices.begin();vert_ref!=(*s_ref)->vertices.end();vert_ref++) |
---|
1832 | { |
---|
1833 | jacobian.set_row(row_count,(*vert_ref)->get_coordinates()-base_vert->get_coordinates()); |
---|
1834 | |
---|
1835 | row_count++; |
---|
1836 | } |
---|
1837 | |
---|
1838 | // cout << gradient << endl; |
---|
1839 | // cout << jacobian << endl; |
---|
1840 | |
---|
1841 | gradient = jacobian*gradient; |
---|
1842 | |
---|
1843 | //cout << gradient << endl; |
---|
1844 | |
---|
1845 | mat rotation_matrix = eye(dimension); |
---|
1846 | |
---|
1847 | double squared_length = 1; |
---|
1848 | |
---|
1849 | for(int i = 1;i<dimension;i++) |
---|
1850 | { |
---|
1851 | double t = gradient[i]/sqrt(squared_length); |
---|
1852 | |
---|
1853 | double sin_theta = t/sqrt(1+pow(t,2)); |
---|
1854 | double cos_theta = 1/sqrt(1+pow(t,2)); |
---|
1855 | |
---|
1856 | mat partial_rotation = eye(dimension); |
---|
1857 | |
---|
1858 | partial_rotation.set(0,0,cos_theta); |
---|
1859 | partial_rotation.set(i,i,cos_theta); |
---|
1860 | |
---|
1861 | partial_rotation.set(0,i,sin_theta); |
---|
1862 | partial_rotation.set(i,0,-sin_theta); |
---|
1863 | |
---|
1864 | rotation_matrix = rotation_matrix*partial_rotation; |
---|
1865 | |
---|
1866 | squared_length += pow(gradient[i],2); |
---|
1867 | } |
---|
1868 | |
---|
1869 | // cout << rotation_matrix << endl; |
---|
1870 | // cout << gradient << endl; |
---|
1871 | |
---|
1872 | vec minima = numeric_limits<double>::max()*ones(number_of_parameters); |
---|
1873 | vec maxima = -numeric_limits<double>::max()*ones(number_of_parameters); |
---|
1874 | |
---|
1875 | vertex* min_grad; |
---|
1876 | vertex* max_grad; |
---|
1877 | |
---|
1878 | for(set<vertex*>::iterator vert_ref = (*s_ref)->vertices.begin();vert_ref!=(*s_ref)->vertices.end();vert_ref++) |
---|
1879 | { |
---|
1880 | vec new_coords = rotation_matrix*jacobian*(*vert_ref)->get_coordinates(); |
---|
1881 | |
---|
1882 | // cout << new_coords << endl; |
---|
1883 | |
---|
1884 | for(int j = 0;j<new_coords.size();j++) |
---|
1885 | { |
---|
1886 | if(new_coords[j]>maxima[j]) |
---|
1887 | { |
---|
1888 | maxima[j] = new_coords[j]; |
---|
1889 | |
---|
1890 | if(j==0) |
---|
1891 | { |
---|
1892 | max_grad = (*vert_ref); |
---|
1893 | } |
---|
1894 | } |
---|
1895 | |
---|
1896 | if(new_coords[j]<minima[j]) |
---|
1897 | { |
---|
1898 | minima[j] = new_coords[j]; |
---|
1899 | |
---|
1900 | if(j==0) |
---|
1901 | { |
---|
1902 | min_grad = (*vert_ref); |
---|
1903 | } |
---|
1904 | } |
---|
1905 | } |
---|
1906 | } |
---|
1907 | |
---|
1908 | vec sample_coordinates; |
---|
1909 | |
---|
1910 | for(int j = 0;j<number_of_parameters;j++) |
---|
1911 | { |
---|
1912 | rnumber = randu(); |
---|
1913 | |
---|
1914 | double coordinate; |
---|
1915 | |
---|
1916 | if(j==0) |
---|
1917 | { |
---|
1918 | double pos_value = sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0]; |
---|
1919 | |
---|
1920 | |
---|
1921 | //coordinate = log(rnumber*(exp(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0])-exp(sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*min_grad->get_coordinates()-sampled_toprow->condition_sum[0]))+ |
---|
1922 | } |
---|
1923 | } |
---|
1924 | |
---|
1925 | // cout << sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*min_grad->get_coordinates()-sampled_toprow->condition_sum[0] << endl; |
---|
1926 | // cout << sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0] << endl; |
---|
1927 | |
---|
1928 | } |
---|
1929 | |
---|
1930 | return sample_mat; |
---|
1931 | } |
---|
1932 | else |
---|
1933 | { |
---|
1934 | throw new exception("You are trying to sample from density that is not determined (parameters can't be integrated out)!"); |
---|
1935 | |
---|
1936 | return 0; |
---|
1937 | } |
---|
1938 | |
---|
1939 | |
---|
1940 | } |
---|
1941 | |
---|
1942 | protected: |
---|
1943 | |
---|
1944 | /// A method for creating plain default statistic representing only the range of the parameter space. |
---|
1945 | void create_statistic(int number_of_parameters) |
---|
1946 | { |
---|
1947 | /* |
---|
1948 | for(int i = 0;i<number_of_parameters;i++) |
---|
1949 | { |
---|
1950 | vec condition_vec = zeros(number_of_parameters+1); |
---|
1951 | condition_vec[i+1] = 1; |
---|
1952 | |
---|
1953 | condition* new_condition = new condition(condition_vec); |
---|
1954 | |
---|
1955 | conditions.push_back(new_condition); |
---|
1956 | } |
---|
1957 | */ |
---|
1958 | |
---|
1959 | // An empty vector of coordinates. |
---|
1960 | vec origin_coord; |
---|
1961 | |
---|
1962 | // We create an origin - this point will have all the coordinates zero, but now it has an empty vector of coords. |
---|
1963 | vertex *origin = new vertex(origin_coord); |
---|
1964 | |
---|
1965 | origin->my_emlig = this; |
---|
1966 | |
---|
1967 | /* |
---|
1968 | // As a statistic, we have to create a vector of vectors of polyhedron pointers. It will then represent the Hasse |
---|
1969 | // diagram. First we create a vector of polyhedrons.. |
---|
1970 | list<polyhedron*> origin_vec; |
---|
1971 | |
---|
1972 | // ..we fill it with the origin.. |
---|
1973 | origin_vec.push_back(origin); |
---|
1974 | |
---|
1975 | // ..and we fill the statistic with the created vector. |
---|
1976 | statistic.push_back(origin_vec); |
---|
1977 | */ |
---|
1978 | |
---|
1979 | statistic = *(new c_statistic()); |
---|
1980 | |
---|
1981 | statistic.append_polyhedron(0, origin); |
---|
1982 | |
---|
1983 | // Now we have a statistic for a zero dimensional space. Regarding to how many dimensional space we need to |
---|
1984 | // describe, we have to widen the descriptional default statistic. We use an iterative procedure as follows: |
---|
1985 | for(int i=0;i<number_of_parameters;i++) |
---|
1986 | { |
---|
1987 | // We first will create two new vertices. These will be the borders of the parameter space in the dimension |
---|
1988 | // of newly added parameter. Therefore they will have all coordinates except the last one zero. We get the |
---|
1989 | // right amount of zero cooridnates by reading them from the origin |
---|
1990 | vec origin_coord = origin->get_coordinates(); |
---|
1991 | |
---|
1992 | // And we incorporate the nonzero coordinates into the new cooordinate vectors |
---|
1993 | vec origin_coord1 = concat(origin_coord,-max_range); |
---|
1994 | vec origin_coord2 = concat(origin_coord,max_range); |
---|
1995 | |
---|
1996 | |
---|
1997 | // Now we create the points |
---|
1998 | vertex* new_point1 = new vertex(origin_coord1); |
---|
1999 | vertex* new_point2 = new vertex(origin_coord2); |
---|
2000 | |
---|
2001 | new_point1->my_emlig = this; |
---|
2002 | new_point2->my_emlig = this; |
---|
2003 | |
---|
2004 | //********************************************************************************************************* |
---|
2005 | // The algorithm for recursive build of a new Hasse diagram representing the space structure from the old |
---|
2006 | // diagram works so that you create two copies of the old Hasse diagram, you shift them up one level (points |
---|
2007 | // will be segments, segments will be areas etc.) and you connect each one of the original copied polyhedrons |
---|
2008 | // with its offspring by a parent-child relation. Also each of the segments in the first (second) copy is |
---|
2009 | // connected to the first (second) newly created vertex by a parent-child relation. |
---|
2010 | //********************************************************************************************************* |
---|
2011 | |
---|
2012 | |
---|
2013 | /* |
---|
2014 | // Create the vectors of vectors of pointers to polyhedrons to hold the copies of the old Hasse diagram |
---|
2015 | vector<vector<polyhedron*>> new_statistic1; |
---|
2016 | vector<vector<polyhedron*>> new_statistic2; |
---|
2017 | */ |
---|
2018 | |
---|
2019 | c_statistic* new_statistic1 = new c_statistic(); |
---|
2020 | c_statistic* new_statistic2 = new c_statistic(); |
---|
2021 | |
---|
2022 | |
---|
2023 | // Copy the statistic by rows |
---|
2024 | for(int j=0;j<statistic.size();j++) |
---|
2025 | { |
---|
2026 | |
---|
2027 | |
---|
2028 | // an element counter |
---|
2029 | int element_number = 0; |
---|
2030 | |
---|
2031 | /* |
---|
2032 | vector<polyhedron*> supportnew_1; |
---|
2033 | vector<polyhedron*> supportnew_2; |
---|
2034 | |
---|
2035 | new_statistic1.push_back(supportnew_1); |
---|
2036 | new_statistic2.push_back(supportnew_2); |
---|
2037 | */ |
---|
2038 | |
---|
2039 | // for each polyhedron in the given row |
---|
2040 | for(polyhedron* horiz_ref = statistic.rows[j];horiz_ref!=statistic.get_end();horiz_ref=horiz_ref->next_poly) |
---|
2041 | { |
---|
2042 | // Append an extra zero coordinate to each of the vertices for the new dimension |
---|
2043 | // If vert_ref is at the first index => we loop through vertices |
---|
2044 | if(j == 0) |
---|
2045 | { |
---|
2046 | // cast the polyhedron pointer to a vertex pointer and push a zero to its vector of coordinates |
---|
2047 | ((vertex*) horiz_ref)->push_coordinate(0); |
---|
2048 | } |
---|
2049 | /* |
---|
2050 | else |
---|
2051 | { |
---|
2052 | ((toprow*) (*horiz_ref))->condition.ins(0,0); |
---|
2053 | }*/ |
---|
2054 | |
---|
2055 | // if it has parents |
---|
2056 | if(!horiz_ref->parents.empty()) |
---|
2057 | { |
---|
2058 | // save the relative address of this child in a vector kids_rel_addresses of all its parents. |
---|
2059 | // This information will later be used for copying the whole Hasse diagram with each of the |
---|
2060 | // relations contained within. |
---|
2061 | for(list<polyhedron*>::iterator parent_ref = horiz_ref->parents.begin();parent_ref != horiz_ref->parents.end();parent_ref++) |
---|
2062 | { |
---|
2063 | (*parent_ref)->kids_rel_addresses.push_back(element_number); |
---|
2064 | } |
---|
2065 | } |
---|
2066 | |
---|
2067 | // ************************************************************************************************** |
---|
2068 | // Here we begin creating a new polyhedron, which will be a copy of the old one. Each such polyhedron |
---|
2069 | // will be created as a toprow, but this information will be later forgotten and only the polyhedrons |
---|
2070 | // in the top row of the Hasse diagram will be considered toprow for later use. |
---|
2071 | // ************************************************************************************************** |
---|
2072 | |
---|
2073 | // First we create vectors specifying a toprow condition. In the case of a preconstructed statistic |
---|
2074 | // this condition will be a vector of zeros. There are two vectors, because we need two copies of |
---|
2075 | // the original Hasse diagram. |
---|
2076 | vec vec1(number_of_parameters+1); |
---|
2077 | vec1.zeros(); |
---|
2078 | |
---|
2079 | vec vec2(number_of_parameters+1); |
---|
2080 | vec2.zeros(); |
---|
2081 | |
---|
2082 | // We create a new toprow with the previously specified condition. |
---|
2083 | toprow* current_copy1 = new toprow(vec1, 0); |
---|
2084 | toprow* current_copy2 = new toprow(vec2, 0); |
---|
2085 | |
---|
2086 | current_copy1->my_emlig = this; |
---|
2087 | current_copy2->my_emlig = this; |
---|
2088 | |
---|
2089 | // The vertices of the copies will be inherited, because there will be a parent/child relation |
---|
2090 | // between each polyhedron and its offspring (comming from the copy) and a parent has all the |
---|
2091 | // vertices of its child plus more. |
---|
2092 | for(set<vertex*>::iterator vertex_ref = horiz_ref->vertices.begin();vertex_ref!=horiz_ref->vertices.end();vertex_ref++) |
---|
2093 | { |
---|
2094 | current_copy1->vertices.insert(*vertex_ref); |
---|
2095 | current_copy2->vertices.insert(*vertex_ref); |
---|
2096 | } |
---|
2097 | |
---|
2098 | // The only new vertex of the offspring should be the newly created point. |
---|
2099 | current_copy1->vertices.insert(new_point1); |
---|
2100 | current_copy2->vertices.insert(new_point2); |
---|
2101 | |
---|
2102 | // This method guarantees that each polyhedron is already triangulated, therefore its triangulation |
---|
2103 | // is only one set of vertices and it is the set of all its vertices. |
---|
2104 | simplex* t_simplex1 = new simplex(current_copy1->vertices); |
---|
2105 | simplex* t_simplex2 = new simplex(current_copy2->vertices); |
---|
2106 | |
---|
2107 | current_copy1->triangulation.insert(t_simplex1); |
---|
2108 | current_copy2->triangulation.insert(t_simplex2); |
---|
2109 | |
---|
2110 | // Now we have copied the polyhedron and we have to copy all of its relations. Because we are copying |
---|
2111 | // in the Hasse diagram from bottom up, we always have to copy the parent/child relations to all the |
---|
2112 | // kids and when we do that and know the child, in the child we will remember the parent we came from. |
---|
2113 | // This way all the parents/children relations are saved in both the parent and the child. |
---|
2114 | if(!horiz_ref->kids_rel_addresses.empty()) |
---|
2115 | { |
---|
2116 | for(list<int>::iterator kid_ref = horiz_ref->kids_rel_addresses.begin();kid_ref!=horiz_ref->kids_rel_addresses.end();kid_ref++) |
---|
2117 | { |
---|
2118 | polyhedron* new_kid1 = new_statistic1->rows[j-1]; |
---|
2119 | polyhedron* new_kid2 = new_statistic2->rows[j-1]; |
---|
2120 | |
---|
2121 | // THIS IS NOT EFFECTIVE: It could be improved by having the list indexed for new_statistic, but |
---|
2122 | // not indexed for statistic. Hopefully this will not cause a big slowdown - happens only offline. |
---|
2123 | if(*kid_ref) |
---|
2124 | { |
---|
2125 | for(int k = 1;k<=(*kid_ref);k++) |
---|
2126 | { |
---|
2127 | new_kid1=new_kid1->next_poly; |
---|
2128 | new_kid2=new_kid2->next_poly; |
---|
2129 | } |
---|
2130 | } |
---|
2131 | |
---|
2132 | // find the child and save the relation to the parent |
---|
2133 | current_copy1->children.push_back(new_kid1); |
---|
2134 | current_copy2->children.push_back(new_kid2); |
---|
2135 | |
---|
2136 | // in the child save the parents' address |
---|
2137 | new_kid1->parents.push_back(current_copy1); |
---|
2138 | new_kid2->parents.push_back(current_copy2); |
---|
2139 | } |
---|
2140 | |
---|
2141 | // Here we clear the parents kids_rel_addresses vector for later use (when we need to widen the |
---|
2142 | // Hasse diagram again) |
---|
2143 | horiz_ref->kids_rel_addresses.clear(); |
---|
2144 | } |
---|
2145 | // If there were no children previously, we are copying a polyhedron that has been a vertex before. |
---|
2146 | // In this case it is a segment now and it will have a relation to its mother (copywise) and to the |
---|
2147 | // newly created point. Here we create the connection to the new point, again from both sides. |
---|
2148 | else |
---|
2149 | { |
---|
2150 | // Add the address of the new point in the former vertex |
---|
2151 | current_copy1->children.push_back(new_point1); |
---|
2152 | current_copy2->children.push_back(new_point2); |
---|
2153 | |
---|
2154 | // Add the address of the former vertex in the new point |
---|
2155 | new_point1->parents.push_back(current_copy1); |
---|
2156 | new_point2->parents.push_back(current_copy2); |
---|
2157 | } |
---|
2158 | |
---|
2159 | // Save the mother in its offspring |
---|
2160 | current_copy1->children.push_back(horiz_ref); |
---|
2161 | current_copy2->children.push_back(horiz_ref); |
---|
2162 | |
---|
2163 | // Save the offspring in its mother |
---|
2164 | horiz_ref->parents.push_back(current_copy1); |
---|
2165 | horiz_ref->parents.push_back(current_copy2); |
---|
2166 | |
---|
2167 | |
---|
2168 | // Add the copies into the relevant statistic. The statistic will later be appended to the previous |
---|
2169 | // Hasse diagram |
---|
2170 | new_statistic1->append_polyhedron(j,current_copy1); |
---|
2171 | new_statistic2->append_polyhedron(j,current_copy2); |
---|
2172 | |
---|
2173 | // Raise the count in the vector of polyhedrons |
---|
2174 | element_number++; |
---|
2175 | |
---|
2176 | } |
---|
2177 | |
---|
2178 | } |
---|
2179 | |
---|
2180 | /* |
---|
2181 | statistic.begin()->push_back(new_point1); |
---|
2182 | statistic.begin()->push_back(new_point2); |
---|
2183 | */ |
---|
2184 | |
---|
2185 | statistic.append_polyhedron(0, new_point1); |
---|
2186 | statistic.append_polyhedron(0, new_point2); |
---|
2187 | |
---|
2188 | // Merge the new statistics into the old one. This will either be the final statistic or we will |
---|
2189 | // reenter the widening loop. |
---|
2190 | for(int j=0;j<new_statistic1->size();j++) |
---|
2191 | { |
---|
2192 | /* |
---|
2193 | if(j+1==statistic.size()) |
---|
2194 | { |
---|
2195 | list<polyhedron*> support; |
---|
2196 | statistic.push_back(support); |
---|
2197 | } |
---|
2198 | |
---|
2199 | (statistic.begin()+j+1)->insert((statistic.begin()+j+1)->end(),new_statistic1[j].begin(),new_statistic1[j].end()); |
---|
2200 | (statistic.begin()+j+1)->insert((statistic.begin()+j+1)->end(),new_statistic2[j].begin(),new_statistic2[j].end()); |
---|
2201 | */ |
---|
2202 | statistic.append_polyhedron(j+1,new_statistic1->rows[j],new_statistic1->row_ends[j]); |
---|
2203 | statistic.append_polyhedron(j+1,new_statistic2->rows[j],new_statistic2->row_ends[j]); |
---|
2204 | } |
---|
2205 | } |
---|
2206 | |
---|
2207 | /* |
---|
2208 | vector<list<toprow*>> toprow_statistic; |
---|
2209 | int line_count = 0; |
---|
2210 | |
---|
2211 | for(vector<list<polyhedron*>>::iterator polyhedron_ref = ++statistic.begin(); polyhedron_ref!=statistic.end();polyhedron_ref++) |
---|
2212 | { |
---|
2213 | list<toprow*> support_list; |
---|
2214 | toprow_statistic.push_back(support_list); |
---|
2215 | |
---|
2216 | for(list<polyhedron*>::iterator polyhedron_ref2 = polyhedron_ref->begin(); polyhedron_ref2 != polyhedron_ref->end(); polyhedron_ref2++) |
---|
2217 | { |
---|
2218 | toprow* support_top = (toprow*)(*polyhedron_ref2); |
---|
2219 | |
---|
2220 | toprow_statistic[line_count].push_back(support_top); |
---|
2221 | } |
---|
2222 | |
---|
2223 | line_count++; |
---|
2224 | }*/ |
---|
2225 | |
---|
2226 | /* |
---|
2227 | vector<int> sizevector; |
---|
2228 | for(int s = 0;s<statistic.size();s++) |
---|
2229 | { |
---|
2230 | sizevector.push_back(statistic.row_size(s)); |
---|
2231 | } |
---|
2232 | */ |
---|
2233 | |
---|
2234 | } |
---|
2235 | |
---|
2236 | |
---|
2237 | |
---|
2238 | |
---|
2239 | }; |
---|
2240 | |
---|
2241 | |
---|
2242 | |
---|
2243 | //! Robust Bayesian AR model for Multicriteria-Laplace-Inverse-Gamma density |
---|
2244 | class RARX //: public BM |
---|
2245 | { |
---|
2246 | private: |
---|
2247 | |
---|
2248 | |
---|
2249 | |
---|
2250 | int window_size; |
---|
2251 | |
---|
2252 | list<vec> conditions; |
---|
2253 | |
---|
2254 | public: |
---|
2255 | emlig* posterior; |
---|
2256 | |
---|
2257 | RARX(int number_of_parameters, const int window_size)//:BM() |
---|
2258 | { |
---|
2259 | posterior = new emlig(number_of_parameters); |
---|
2260 | |
---|
2261 | this->window_size = window_size; |
---|
2262 | }; |
---|
2263 | |
---|
2264 | void bayes(const itpp::vec &yt, const itpp::vec &cond = "") |
---|
2265 | { |
---|
2266 | conditions.push_back(yt); |
---|
2267 | |
---|
2268 | //posterior->step_me(0); |
---|
2269 | |
---|
2270 | /// \TODO tohle je spatne, tady musi byt jiny vypocet poctu podminek, kdyby nejaka byla multiplicitni, tak tohle bude spatne |
---|
2271 | if(conditions.size()>window_size && window_size!=0) |
---|
2272 | { |
---|
2273 | posterior->add_and_remove_condition(yt,conditions.front()); |
---|
2274 | conditions.pop_front(); |
---|
2275 | |
---|
2276 | //posterior->step_me(1); |
---|
2277 | } |
---|
2278 | else |
---|
2279 | { |
---|
2280 | posterior->add_condition(yt); |
---|
2281 | } |
---|
2282 | |
---|
2283 | } |
---|
2284 | |
---|
2285 | }; |
---|
2286 | |
---|
2287 | |
---|
2288 | |
---|
2289 | #endif //TRAGE_H |
---|