Changeset 1320 for applications/robust

Timestamp:

03/31/11 18:46:06 (14 years ago)

Author:

sindj

Message:

Rozdelano samplovani z emlig. Zbyva dodelat vypocet v jednotlivych simplexech triangulace. JS

Location:

applications/robust

Files:

• robustlib.cpp (modified) (2 diffs)
• robustlib.h (modified) (7 diffs)

applications/robust/robustlib.cpp

@@ -16 +16 @@
                         vert_simplex.insert((vertex*)this);
 
-                        triangulation.push_back(vert_simplex);
+                        triangulation.insert(pair<double,set<vertex*>>(0,vert_simplex));
                 }
 
                 for(list<polyhedron*>::iterator child_ref = children.begin();child_ref!=children.end();child_ref++)
                 {                       
-                        for(list<set<vertex*>>::iterator t_ref = (*child_ref)->triangulation.begin();t_ref!=(*child_ref)->triangulation.end();t_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).begin(),(*t_ref).end());                            
+                                new_simplex.insert((*t_ref).second.begin(),(*t_ref).second.end());                              
 
                                 pair<set<vertex*>::iterator,bool> ret_val = new_simplex.insert(*vertices.begin());
@@ -30 +30 @@
                                 if(ret_val.second == true)
                                 {
-                                        triangulation.push_back(new_simplex);
+                                        double cur_prob = 0;
 
                                         if(should_integrate)
                                         {
-                                                ((toprow *)this)->probability += ((toprow *)this)->integrate_simplex(new_simplex, 'S');
-                                        }
+                                                cur_prob = ((toprow *)this)->integrate_simplex(new_simplex, 'S');
+
+                                                ((toprow *)this)->probability += cur_prob;
+                                        }
+
+                                        triangulation.insert(pair<double,set<vertex*>>(cur_prob,new_simplex));
                                 } 
                         }       

applications/robust/robustlib.h

@@ -164 +164 @@
 
         /// List of triangulation polyhedrons of the polyhedron given by their relative vertices. 
-        list<set<vertex*>> triangulation;
+        map<double,set<vertex*>> triangulation;
 
         /// A list of relative addresses serving for Hasse diagram construction.
@@ -279 +279 @@
                 vert_simplex.insert(this);
 
-                triangulation.push_back(vert_simplex);
+                triangulation.insert(pair<double,set<vertex*>>(0,vert_simplex));
         }
 
@@ -773 +773 @@
                                                                 toprow* cur_par_toprow = ((toprow*)current_parent);
                                                                 cur_par_toprow->probability = 0.0;
-                                                                                
-                                                                for(list<set<vertex*>>::iterator t_ref = current_parent->triangulation.begin();t_ref!=current_parent->triangulation.end();t_ref++)
+                                                                
+                                                                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++)
                                                                 {
-                                                                        cur_par_toprow->probability += cur_par_toprow->integrate_simplex(*t_ref,'C');
-                                                                }                                                                       
+                                                                        double cur_prob = cur_par_toprow->integrate_simplex((*t_ref).second,'C');
+                                                                        
+                                                                        cur_par_toprow->probability += cur_prob;
+
+                                                                        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());
                                                         }
                                                 }
@@ -870 +879 @@
                                                                 cur_par_toprow->probability = 0.0;
                                                                         
-                                                                for(list<set<vertex*>>::iterator t_ref = current_parent->triangulation.begin();t_ref!=current_parent->triangulation.end();t_ref++)
+                                                                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++)
                                                                 {
-                                                                        cur_par_toprow->probability += cur_par_toprow->integrate_simplex(*t_ref,'C');
-                                                                }                                                                       
+                                                                        double cur_prob = cur_par_toprow->integrate_simplex((*t_ref).second,'C');
+                                                                        
+                                                                        cur_par_toprow->probability += cur_prob;
+
+                                                                        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());
                                                         }
 
@@ -1690 +1708 @@
                 }
 
+
+
+        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)
+                        {
+                                toprow* current_top = (toprow*)top_ref;
+
+                                sum+=current_top->probability;
+                                ordered_toprows.insert(pair<double,toprow*>(sum,current_top));
+                        }
+                }
+                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;
+        }
+
 protected:
 
@@ -1858 +1926 @@
                                         t_simplex2.insert(current_copy2->vertices.begin(),current_copy2->vertices.end());
                                         
-                                        current_copy1->triangulation.push_back(t_simplex1);
-                                        current_copy2->triangulation.push_back(t_simplex2);                                     
+                                        current_copy1->triangulation.insert(pair<double,set<vertex*>>(0,t_simplex1));
+                                        current_copy2->triangulation.insert(pair<double,set<vertex*>>(0,t_simplex2));                                   
                                         
                                         // Now we have copied the polyhedron and we have to copy all of its relations. Because we are copying
@@ -2021 +2089 @@
                 //posterior->step_me(0);
                 
+                /// \TODO tohle je spatne, tady musi byt jiny vypocet poctu podminek, kdyby nejaka byla multiplicitni, tak tohle bude spatne
                 if(conditions.size()>window_size && window_size!=0)
                 {