Changeset 1254

Timestamp:

11/19/10 11:15:48 (14 years ago)

Author:

sindj

Message:

Rozdelane pocitani integralu. JS

Location:

applications/robust

Files:

• main.cpp (modified) (1 diff)
• robustlib.cpp (modified) (1 diff)
• robustlib.h (modified) (3 diffs)

applications/robust/main.cpp

@@ -23 +23 @@
                 cout << "Step:" << i << endl;
 
-                double condition[9];            
+                double condition[4];            
 
                 condition[0] = rand()/1000.0;

applications/robust/robustlib.cpp

@@ -1 +1 @@
 #include "robustlib.h"
 
+void polyhedron::triangulate(bool should_integrate)
+        {               
+                for(list<polyhedron*>::iterator child_ref = children.begin();child_ref!=children.end();child_ref++)
+                {
+                        set<double> simplex_integrals;
+                        
+                        for(list<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());                            
+
+                                pair<set<vertex*>::iterator,bool> ret_val = new_simplex.insert(*vertices.begin());
+
+                                if(ret_val.second == true)
+                                {
+                                        triangulation.push_back(new_simplex);
+
+                                        if(should_integrate)
+                                        {
+                                                toprow* as_toprow = (toprow*)this;
+
+                                                vec cur_condition = as_toprow->condition.get(1,as_toprow->condition.size()-1);
+
+                                                vertex* base_vertex = (*new_simplex.begin());
+
+                                                int dimension = new_simplex.size()-1;
+
+                                                mat jacobian(dimension,dimension);
+
+                                                double a_0 = base_vertex->get_coordinates()*cur_condition+as_toprow->condition[0];
+                                                vec as;
+
+                                                int row_count = 0;
+
+                                                for(set<vertex*>::iterator vert_ref = (++new_simplex.begin()); vert_ref!=new_simplex.end();vert_ref++)
+                                                {
+                                                        vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates();
+                                                        jacobian.set_row(row_count,relative_coords);
+
+                                                        double new_a = relative_coords*cur_condition;                                                   
+                                                                                                                
+                                                        as.ins(as.size(),new_a);                                                        
+                                                        
+                                                        row_count++;
+                                                }
+
+                                                double int_value = 0;
+                                                for(int a_count = 0;a_count<as.size();a_count++)
+                                                {
+                                                        vec reduced_as = (as.get(a_count)-as);
+
+                                                        reduced_as.del(a_count);
+
+                                                        int fac_order = ((toprow*)this)->condition_order-reduced_as.size()-1;
+
+                                                        double fac_value = ((double)tgamma(fac_order)*det(jacobian))/as[a_count]/pow(as[a_count]-a_0,fac_order); 
+
+                                                        for(int b_count = 0;b_count<reduced_as.size();b_count++)
+                                                        {
+                                                                fac_value /= reduced_as[b_count];                                                               
+                                                        }
+
+                                                        int_value += fac_value;                                                 
+                                                }
+
+                                                simplex_integrals.insert(int_value);
+                                        }
+                                } 
+                        }
+
+                        if(should_integrate)
+                        {
+                                ((toprow*)this)->probability = 0.0;
+
+                                for(set<double>::iterator integ_ref = simplex_integrals.begin();integ_ref!=simplex_integrals.end();integ_ref++)
+                                {
+                                        ((toprow*)this)->probability += (*integ_ref);
+                                }
+                        }
+                }               
+        }
+

applications/robust/robustlib.h

@@ -19 +19 @@
 using namespace itpp;
 
-const double max_range = 99999999999999999999999.0;//numeric_limits<double>::max()/10e-10;
+const double max_range = 1000.0;//numeric_limits<double>::max()/10e-10;
 
 enum actions {MERGE, SPLIT};
@@ -25 +25 @@
 class polyhedron;
 class vertex;
+class toprow;
 
 /*
@@ -170 +171 @@
         }
 
-        void triangulate(bool should_integrate)
-        {               
-                for(list<polyhedron*>::iterator child_ref = children.begin();child_ref!=children.end();child_ref++)
-                {
-                        set<double> simplex_integrals;
-                        
-                        for(list<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());                            
-
-                                pair<set<vertex*>::iterator,bool> ret_val = new_simplex.insert(*vertices.begin());
-
-                                if(ret_val.second == true)
-                                {
-                                        triangulation.push_back(new_simplex);
-
-                                        if(should_integrate)
-                                        {
-                                                toprow* as_toprow = (toprow*)this;                                              
-
-                                                vertex* base_vertex = (*new_simplex.begin());
-
-                                                double a_0 = base_vertex->get_coordinates()*as_toprow->condition(1,as_toprow->condition.size()-1)+as_toprow->condition[0];
-                                                list<double> as;
-
-                                                for(set<vertex*>::iterator vert_ref = (++new_simplex.begin()); vert_ref!=new_simplex.end();vert_ref++)
-                                                {
-                                                        vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates();
-
-                                                        double new_a = relative_coords*as_toprow->condition(1,as_toprow->condition.size()-1);
-
-                                                        as.push_back(new_a);
-                                                }
-
-                                                for(list<double>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++)
-                                                {
-                                                        //TODO DODELAT VYPOCET!
-                                                }
-                                        }
-                                } 
-                        }
-
-                        if(should_integrate)
-                        {
-                                ((toprow*)this)->probability = 0.0;
-
-                                for(set<double>::iterator integ_ref = simplex_integrals.begin();integ_ref!=simplex_integrals.end();integ_ref++)
-                                {
-                                        ((toprow*)this)->probability += (*integ_ref);
-                                }
-                        }
-                }               
-        }
+        void triangulate(bool should_integrate);
+