Changeset 1262 for applications/robust/robustlib.cpp

Timestamp:

12/13/10 19:25:45 (14 years ago)

Author:

sindj

Message:

Dodelavani korekcnich faktoru pri stjene hodnote funkce ve dvou bodech polyhedronu. Rozdelano. JS

Files:

• applications/robust/robustlib.cpp (modified) (3 diffs)

applications/robust/robustlib.cpp

@@ -31 +31 @@
 
                                                 double a_0 = base_vertex->get_coordinates()*cur_condition+as_toprow->condition[0];
-                                                vec as;
+                                                map<double,int> as;
 
                                                 int row_count = 0;
@@ -41 +41 @@
 
                                                         double new_a = relative_coords*cur_condition;                                                   
-                                                                                                                
-                                                        as.ins(as.size(),new_a);                                                        
+                                                        
+                                                        pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(new_a,1));
+                                                        if(returned.second == false)
+                                                        {
+                                                                (*returned.first).second++;
+                                                        }
+                                                        //as.ins(as.size(),new_a);                                                      
                                                         
                                                         row_count++;
@@ -48 +53 @@
 
                                                 double int_value = 0;
-                                                for(int a_count = 0;a_count<as.size();a_count++)
+                                                for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++)
                                                 {
-                                                        vec reduced_as = (as.get(a_count)-as);
+                                                        int fac_order = ((toprow*)this)->condition_order-as.size()-1;
 
-                                                        reduced_as.del(a_count);
+                                                        double fac_value = 1;
+                                                        if(a_0!=(*as_ref).first)
+                                                        {
+                                                                //fac_value = ((double)tgamma(fac_order)*det(jacobian))/(*as_ref).first/pow((*as_ref).first-a_0,fac_order); 
+                                                                
+                                                                fac_value /= (*as_ref).first;
 
-                                                        int fac_order = ((toprow*)this)->condition_order-reduced_as.size()-1;
+                                                                for(map<double,int>::iterator as2_ref = as.begin();as2_ref!=as.end();as2_ref++)
+                                                                {
+                                                                        if(as2_ref!=as_ref)
+                                                                        {
+                                                                                double current_factor = (*as_ref).first-(*as2_ref).first;
 
-                                                        double fac_value = ((double)tgamma(fac_order)*det(jacobian))/as[a_count]/pow(as[a_count]-a_0,fac_order); 
+                                                                                fac_value /= current_factor;
 
-                                                        for(int b_count = 0;b_count<reduced_as.size();b_count++)
+                                                                                                                                                                
+                                                                        }
+                                                                }
+                                                        }
+                                                        else
                                                         {
-                                                                fac_value /= reduced_as[b_count];                                                               
+                                                                throw new exception("Should this happen? a_i=a_0 in the formula for integrating over a simplex! I think it can happen with 0 probability(when phi*2*v_0=phi*v_i)");
                                                         }
 
+                                                        
                                                         int_value += fac_value;                                                 
                                                 }