root/applications/robust/robustlib.h @ 1267

Revision 1267, 38.5 kB (checked in by sindj, 13 years ago)

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