root/applications/robust/robustlib.h @ 1381

Revision 1379, 84.1 kB (checked in by sindj, 13 years ago)

Prepracovani integrace v robustlib.cpp, zjednoduseni, zesymetrizovani. Nevim jestli funguje, nutno overit a vyzkouset testovani hypotez a integraci na jednicku (u normalizacnich faktoru).JS

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