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