Changeset 1324

Timestamp:

04/07/11 18:33:51 (14 years ago)

Author:

sindj

Message:

Temer dodelane samplovani, rozdelany samling sigma, sampling alpha temer dokoncen. JS

Location:

applications/robust

Files:

• main.cpp (modified) (3 diffs)
• robustlib.cpp (modified) (9 diffs)
• robustlib.h (modified) (13 diffs)

applications/robust/main.cpp

@@ -134 +134 @@
         vector<vector<string>> strings;
 
-        char* file_strings[3] = {"c:\\ar_student_single.txt", "c:\\ar_cauchy_single.txt","c:\\ar_normal_single.txt"};
+        char* file_strings[3] = {"c:\\ar_cauchy_single.txt", "c:\\ar_normal_single.txt","c:\\ar_student_single.txt"};
 
         for(int i = 0;i<3;i++)
@@ -162 +162 @@
                 vector<vec> conditions;
                 //emlig* emliga = new emlig(2);
-                RARX* my_rarx = new RARX(2,30);
+                RARX* my_rarx = new RARX(2,70);
 
                 for(int k = 1;k<170;k++)
@@ -191 +191 @@
                                 if(k>5)
                                 {
+                                        //my_rarx->posterior->step_me(0);
+
+                                        //my_rarx->posterior->sample_mat(1);
+
                                         cout << "MaxLik coords:" << my_rarx->posterior->minimal_vertex->get_coordinates() << endl;
 

applications/robust/robustlib.cpp

@@ -4 +4 @@
 void polyhedron::triangulate(bool should_integrate)
         {
+                for(set<simplex*>::iterator t_ref = triangulation.begin();t_ref!=triangulation.end();t_ref++)
+                {
+                        delete (*t_ref);
+                }               
                 triangulation.clear();
 
@@ -13 +17 @@
                 if(vertices.size()==1)
                 {
-                        set<vertex*> vert_simplex;
-                        vert_simplex.insert((vertex*)this);
-
-                        triangulation.insert(pair<double,set<vertex*>>(0,vert_simplex));
+                        simplex* vert_simplex = new simplex((vertex*)this);                     
+
+                        triangulation.insert(vert_simplex);
                 }
 
                 for(list<polyhedron*>::iterator child_ref = children.begin();child_ref!=children.end();child_ref++)
                 {                       
-                        for(map<double,set<vertex*>>::iterator t_ref = (*child_ref)->triangulation.begin();t_ref!=(*child_ref)->triangulation.end();t_ref++)
-                        {
-                                set<vertex*> new_simplex;
-                                new_simplex.insert((*t_ref).second.begin(),(*t_ref).second.end());                              
-
-                                pair<set<vertex*>::iterator,bool> ret_val = new_simplex.insert(*vertices.begin());
+                        for(set<simplex*>::iterator s_ref = (*child_ref)->triangulation.begin();s_ref!=(*child_ref)->triangulation.end();s_ref++)
+                        {
+                                simplex* new_simplex = new simplex((*s_ref)->vertices);                                                         
+
+                                pair<set<vertex*>::iterator,bool> ret_val = new_simplex->vertices.insert(*vertices.begin());
 
                                 if(ret_val.second == true)
@@ -39 +41 @@
                                         }
 
-                                        triangulation.insert(pair<double,set<vertex*>>(cur_prob,new_simplex));
-                                } 
+                                        triangulation.insert(new_simplex);
+                                }
+                                else
+                                {
+                                        delete new_simplex;
+                                }
                         }       
                 }               
@@ -71 +77 @@
 
 
-        double toprow::integrate_simplex(set<vertex*> simplex, char c)
+        double toprow::integrate_simplex(simplex* simplex, char c)
         {
                 // cout << ((toprow*)this)->condition << endl;          
@@ -107 +113 @@
                         // cout << as_toprow->condition << endl;                        
 
-                        int dimension = simplex.size()-1;
+                        int dimension = simplex->vertices.size()-1;
 
                         mat jacobian(dimension,dimension);                      
@@ -114 +120 @@
                         map<double,int> as;
                         vertex* base_vertex;
-                        set<vertex*>::iterator base_ref = simplex.begin();
+                        set<vertex*>::iterator base_ref = simplex->vertices.begin();
                         bool order_correct;
                                                 
@@ -132 +138 @@
 
                                 int row_count = 0;
-                                for(set<vertex*>::iterator vert_ref = simplex.begin(); vert_ref!=simplex.end();vert_ref++)
+                                for(set<vertex*>::iterator vert_ref = simplex->vertices.begin(); vert_ref!=simplex->vertices.end();vert_ref++)
                                 {
                                         if(vert_ref != base_ref)
@@ -146 +152 @@
                                                         base_ref++;
 
-                                                        if(base_ref == simplex.end())
+                                                        if(base_ref == simplex->vertices.end())
                                                         {
                                                                 throw new exception("Equal local conditions are paired. If this ever occurs, the software has to be modified to include multiplied a_0!!");
@@ -285 +291 @@
                 
 
-                        //cout << "Probability:" << int_value << endl;
-                        
-                        return int_value;
-
-                        
+                        // cout << "Probability:" << int_value << endl;
+                        // pause(0.100);
+
+                        simplex->probability = int_value;
+                        
+                        return int_value;       
                         
                 }

applications/robust/robustlib.h

@@ -31 +31 @@
 class emlig;
 
+
 /*
 class t_simplex
@@ -61 +62 @@
                 this->value = value;
                 multiplicity = 1;
+        }
+};
+
+class simplex
+{
+public:
+
+        set<vertex*> vertices;
+
+        double probability;
+
+        vector<list<pair<double,double>>> gamma_parameters;
+
+        simplex(set<vertex*> vertices)
+        {
+                this->vertices.insert(vertices.begin(),vertices.end());
+                probability = 0;
+        }
+
+        simplex(vertex* vertex)
+        {
+                this->vertices.insert(vertex);
+                probability = 0;
         }
 };
@@ -164 +188 @@
 
         /// List of triangulation polyhedrons of the polyhedron given by their relative vertices. 
-        map<double,set<vertex*>> triangulation;
+        set<simplex*> triangulation;
 
         /// A list of relative addresses serving for Hasse diagram construction.
@@ -275 +299 @@
                 vertices.insert(this);
 
-                set<vertex*> vert_simplex;
-
-                vert_simplex.insert(this);
-
-                triangulation.insert(pair<double,set<vertex*>>(0,vert_simplex));
+                simplex* vert_simplex = new simplex(vertices);          
+
+                triangulation.insert(vert_simplex);
         }
 
@@ -330 +352 @@
         }
 
-        double integrate_simplex(set<vertex*> simplex, char c);
+        double integrate_simplex(simplex* simplex, char c);
 
 };
+
 
 
@@ -774 +797 @@
                                                                 cur_par_toprow->probability = 0.0;
                                                                 
-                                                                map<double,set<vertex*>> new_triangulation;
-
-                                                                for(map<double,set<vertex*>>::iterator t_ref = current_parent->triangulation.begin();t_ref!=current_parent->triangulation.end();t_ref++)
+                                                                //set<simplex*> new_triangulation;
+
+                                                                for(set<simplex*>::iterator s_ref = current_parent->triangulation.begin();s_ref!=current_parent->triangulation.end();s_ref++)
                                                                 {
-                                                                        double cur_prob = cur_par_toprow->integrate_simplex((*t_ref).second,'C');
+                                                                        double cur_prob = cur_par_toprow->integrate_simplex((*s_ref),'C');
                                                                         
                                                                         cur_par_toprow->probability += cur_prob;
 
-                                                                        new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second));
+                                                                        //new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second));
                                                                 }
 
-                                                                current_parent->triangulation.clear();
-                                                                current_parent->triangulation.insert(new_triangulation.begin(),new_triangulation.end());
+                                                                //current_parent->triangulation.clear();
+                                                                //current_parent->triangulation.insert(new_triangulation.begin(),new_triangulation.end());
                                                         }
                                                 }
@@ -879 +902 @@
                                                                 cur_par_toprow->probability = 0.0;
                                                                         
-                                                                map<double,set<vertex*>> new_triangulation;
+                                                                //map<double,set<vertex*>> new_triangulation;
                                                                 
-                                                                for(map<double,set<vertex*>>::iterator t_ref = current_parent->triangulation.begin();t_ref!=current_parent->triangulation.end();t_ref++)
+                                                                for(set<simplex*>::iterator s_ref = current_parent->triangulation.begin();s_ref!=current_parent->triangulation.end();s_ref++)
                                                                 {
-                                                                        double cur_prob = cur_par_toprow->integrate_simplex((*t_ref).second,'C');
+                                                                        double cur_prob = cur_par_toprow->integrate_simplex((*s_ref),'C');
                                                                         
                                                                         cur_par_toprow->probability += cur_prob;
 
-                                                                        new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second));
+                                                                        //new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second));
                                                                 }
 
-                                                                current_parent->triangulation.clear();
-                                                                current_parent->triangulation.insert(new_triangulation.begin(),new_triangulation.end());
+                                                                //current_parent->triangulation.clear();
+                                                                //current_parent->triangulation.insert(new_triangulation.begin(),new_triangulation.end());
                                                         }
 
@@ -960 +983 @@
                 for(int i = 0;i<statistic.size();i++)
                 {
-                        int zero = 0;
-                        int one  = 0;
-                        int two  = 0;
+                        //int zero = 0;
+                        //int one  = 0;
+                        //int two  = 0;
 
                         for(polyhedron* horiz_ref = statistic.rows[i];horiz_ref!=statistic.get_end();horiz_ref=horiz_ref->next_poly)
                         {
+                                
+                                
+                                if(i==statistic.size()-1)
+                                {
+                                        cout << ((toprow*)horiz_ref)->condition_sum << "   " << ((toprow*)horiz_ref)->probability << endl;
+                                        cout << "Order:" << ((toprow*)horiz_ref)->condition_order << endl;
+                                }
                                 
                                 /*
-                                if(i==statistic.size()-1)
-                                {
-                                        //cout << ((toprow*)horiz_ref)->condition_sum << "   " << ((toprow*)horiz_ref)->probability << endl;
-                                        cout << "Order:" << ((toprow*)horiz_ref)->condition_order << endl;
-                                }
                                 if(i==0)
                                 {
@@ -979 +1004 @@
                                 */
 
+                                /*
                                 char* string = "Checkpoint";
 
+
                                 if((*horiz_ref).parentconditions.size()==0)
                                 {
@@ -993 +1020 @@
                                         two++;
                                 }
+                                */
                                 
                         }
@@ -1070 +1098 @@
         void add_and_remove_condition(vec toadd, vec toremove)
         {
-                step_me(0);
+                //step_me(0);
                 
                 likelihood_value = numeric_limits<double>::max();
@@ -1712 +1740 @@
         mat sample_mat(int n)
         {
-                mat sample_mat;
-                map<double,toprow*> ordered_toprows;
+                
+                
 
                 /// \TODO tady je to spatne, tady nesmi byt conditions.size(), viz RARX.bayes()
                 if(conditions.size()-2-number_of_parameters>=0)
                 {
-                        double sum = 0;
-
-                        for(polyhedron* top_ref = statistic.rows[number_of_parameters-1];top_ref!=statistic.end_poly;top_ref=top_ref->next_poly)
+                        mat sample_mat;
+                        map<double,toprow*> ordered_toprows;                    
+                        double sum_a = 0;
+                        
+
+                        for(polyhedron* top_ref = statistic.rows[number_of_parameters];top_ref!=statistic.end_poly;top_ref=top_ref->next_poly)
                         {
                                 toprow* current_top = (toprow*)top_ref;
 
-                                sum+=current_top->probability;
-                                ordered_toprows.insert(pair<double,toprow*>(sum,current_top));
-                        }
+                                sum_a+=current_top->probability;
+                                ordered_toprows.insert(pair<double,toprow*>(sum_a,current_top));
+                        }
+
+                        while(sample_mat.cols()<n)
+                        {
+                                double rnumber = randu()*sum_a;
+
+                                // cout << "RND:" << rnumber << endl;
+
+                                // This could be more efficient (log n), but map::upper_bound() doesn't let me dereference returned iterator
+                                toprow* sampled_toprow;                         
+                                for(map<double,toprow*>::iterator top_ref = ordered_toprows.begin();top_ref!=ordered_toprows.end();top_ref++)
+                                {
+                                        // cout << "CDF:"<< (*top_ref).first << endl;
+
+                                        if((*top_ref).first >= rnumber)
+                                        {
+                                                sampled_toprow = (*top_ref).second;
+                                                break;
+                                        }                                               
+                                }                               
+
+                                // cout << &sampled_toprow << ";";
+
+                                rnumber = randu();                              
+
+                                set<simplex*>::iterator s_ref;
+                                double sum_b = 0;
+
+                                for(s_ref = sampled_toprow->triangulation.begin();s_ref!=sampled_toprow->triangulation.end();s_ref++)
+                                {
+                                        sum_b += (*s_ref)->probability;
+
+                                        if(sum_b/sampled_toprow->probability >= rnumber)
+                                                break;
+                                }
+
+                                // cout << &(*tri_ref) << endl;
+
+                                int dimension = (*s_ref)->vertices.size()-1;
+
+                                mat jacobian(dimension,dimension);
+                                vec gradient = sampled_toprow->condition_sum.right(dimension);
+
+                                vertex* base_vert = *(*s_ref)->vertices.begin();
+                                
+                                int row_count = 0;
+
+                                for(set<vertex*>::iterator vert_ref = ++(*s_ref)->vertices.begin();vert_ref!=(*s_ref)->vertices.end();vert_ref++)
+                                {
+                                        jacobian.set_row(row_count,(*vert_ref)->get_coordinates()-base_vert->get_coordinates());
+
+                                        row_count++;
+                                }
+
+                                // cout << gradient << endl;
+                                // cout << jacobian << endl;
+
+                                gradient = jacobian*gradient;   
+
+                                //cout << gradient << endl;
+
+                                mat rotation_matrix = eye(dimension);
+
+                                double squared_length = 1;
+
+                                for(int i = 1;i<dimension;i++)
+                                {
+                                        double t = gradient[i]/sqrt(squared_length);
+
+                                        double sin_theta = t/sqrt(1+pow(t,2));
+                                        double cos_theta = 1/sqrt(1+pow(t,2));
+
+                                        mat partial_rotation = eye(dimension);
+
+                                        partial_rotation.set(0,0,cos_theta);
+                                        partial_rotation.set(i,i,cos_theta);
+                                        
+                                        partial_rotation.set(0,i,sin_theta);
+                                        partial_rotation.set(i,0,-sin_theta);
+                                        
+                                        rotation_matrix = rotation_matrix*partial_rotation;
+                                        
+                                        squared_length += pow(gradient[i],2);
+                                }
+
+                                // cout << rotation_matrix << endl;
+                                // cout << gradient << endl;
+
+                                vec minima = numeric_limits<double>::max()*ones(number_of_parameters);
+                                vec maxima = -numeric_limits<double>::max()*ones(number_of_parameters);
+
+                                vertex* min_grad;
+                                vertex* max_grad;
+                                
+                                for(set<vertex*>::iterator vert_ref = (*s_ref)->vertices.begin();vert_ref!=(*s_ref)->vertices.end();vert_ref++)
+                                {
+                                        vec new_coords = rotation_matrix*jacobian*(*vert_ref)->get_coordinates();
+
+                                        // cout << new_coords << endl;
+
+                                        for(int j = 0;j<new_coords.size();j++)
+                                        {
+                                                if(new_coords[j]>maxima[j])
+                                                {
+                                                        maxima[j] = new_coords[j];
+
+                                                        if(j==0)
+                                                        {
+                                                                max_grad = (*vert_ref);
+                                                        }
+                                                }
+
+                                                if(new_coords[j]<minima[j])
+                                                {
+                                                        minima[j] = new_coords[j];
+
+                                                        if(j==0)
+                                                        {
+                                                                min_grad = (*vert_ref);
+                                                        }
+                                                }
+                                        }                                       
+                                }
+
+                                vec sample_coordinates;
+
+                                for(int j = 0;j<number_of_parameters;j++)
+                                {
+                                        rnumber = randu();
+                                        
+                                        double coordinate;
+
+                                        if(j==0)
+                                        {
+                                                double pos_value = sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0];
+                                                
+                                                
+                                                //coordinate = log(rnumber*(exp(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0])-exp(sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*min_grad->get_coordinates()-sampled_toprow->condition_sum[0]))+
+                                        }
+                                }
+
+                                // cout << sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*min_grad->get_coordinates()-sampled_toprow->condition_sum[0] << endl;
+                                // cout << sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0] << endl;
+
+                        }
+
+                        return sample_mat;
                 }
                 else
                 {
                         throw new exception("You are trying to sample from density that is not determined (parameters can't be integrated out)!");
-                }
-
-                while(sample_mat.cols()<n)
-                {
-                        double rnumber = randu();
-
-                        toprow* sampled_toprow = ordered_toprows.upper_bound(rnumber)->second;
-
-                        rnumber = randu();
-
-                        map<double,set<vertex*>>::iterator tri_ref = sampled_toprow->triangulation.begin()
-                        double sum = 0;
-
-                        for(tri_ref = sampled_toprow->triangulation.begin();tri_ref!=sampled_toprow->triangulation.end();tri_ref++)
-                        {
-                                sum += (*tri_ref).first;
-
-                                if(sum/sampled_toprow->probability >= rnumber)
-                                        break;
-                        }
-
-
-                }
-
-                return sample_mat;
+                
+                        return 0;
+                }
+
+                
         }
 
@@ -1920 +2077 @@
                                         // This method guarantees that each polyhedron is already triangulated, therefore its triangulation
                                         // is only one set of vertices and it is the set of all its vertices.
-                                        set<vertex*> t_simplex1;
-                                        set<vertex*> t_simplex2;
-
-                                        t_simplex1.insert(current_copy1->vertices.begin(),current_copy1->vertices.end());
-                                        t_simplex2.insert(current_copy2->vertices.begin(),current_copy2->vertices.end());
+                                        simplex* t_simplex1 = new simplex(current_copy1->vertices);
+                                        simplex* t_simplex2 = new simplex(current_copy2->vertices);                                     
                                         
-                                        current_copy1->triangulation.insert(pair<double,set<vertex*>>(0,t_simplex1));
-                                        current_copy2->triangulation.insert(pair<double,set<vertex*>>(0,t_simplex2));                                   
+                                        current_copy1->triangulation.insert(t_simplex1);
+                                        current_copy2->triangulation.insert(t_simplex2);                                        
                                         
                                         // Now we have copied the polyhedron and we have to copy all of its relations. Because we are copying