Changeset 1262 for applications/robust

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

Location:

applications/robust

Files:

• robustlib.cpp (modified) (3 diffs)
• robustlib.h (modified) (5 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;                                                 
                                                 }

applications/robust/robustlib.h

@@ -12 +12 @@
 #include <vector>
 #include <list>
+#include <map>
 #include <set>
 #include <algorithm>
+#include <string>
         
 using namespace bdm;
@@ -413 +415 @@
 class emlig // : eEF
 {
+        vector<set<ivec>> correction_factors;
 
         /// A statistic in a form of a Hasse diagram representing a complex of convex polyhedrons obtained as a result
@@ -425 +428 @@
 
         double normalization_factor;
+
+        
 
         void alter_toprow_conditions(vec condition, bool should_be_added)
@@ -915 +920 @@
                         sizevector.push_back(statistic.row_size(s));
                 }*/
-        }
+
+                
+        }
+
+        vector<list<ivec>> get_correction_factors(int order)
+                {
+                        for(int remaining_order = correction_factors.size();remaining_order>order;remaining_order++)
+                        {
+                                set<ivec> factor_templates;
+                                set<ivec> final_factors;
+
+                                for(int i = 1;i==number_of_parameters-order+1;i++)
+                                {
+                                        for(int j = 1;j==remaining_order;j++)
+                                        {
+                                                factor_templates.insert(zeros(number_of_parameters-order+2));
+
+                                                for(set<ivec>::iterator fac_ref = factor_templates.begin();fac_ref!=factor_templates.end();fac_ref++)
+                                                {
+                                                        ivec current_template = (*fac_ref);
+
+                                                        current_template[0]+=1;
+                                                        current_template[i]+=1;
+
+                                                        if(current_template[0]==remaining_order)
+                                                        {
+                                                                final_factors.insert(current_template.right(current_template.size()-1);
+                                                        }
+                                                        else
+                                                        {
+                                                                factor_templates.insert(current_template);
+                                                        }
+                                                }                                       
+                                        }
+                                }       
+
+                                correction_factors.push_back(final_factors);
+
+                        }
+                }
 
 protected:
@@ -1234 +1278 @@
 
 
+
 #endif //TRAGE_H