root/applications/robust/robustlib.h @ 1269

Revision 1269, 38.4 kB (checked in by sindj, 13 years ago)

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