root/applications/robust/robustlib.h @ 1281

Revision 1281, 40.5 kB (checked in by sindj, 13 years ago)

Tak to vypada, ze to konecne pocita spravne, tak to radsi commitnu, protoze vlastne nevim, jak jsem to udelal :-) JS

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