root/applications/robust/robustlib.h @ 1318

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