Changeset 1272 for applications/robust/robustlib.cpp

Timestamp:

01/14/11 18:28:55 (13 years ago)

Author:

sindj

Message:

Dodelano pocitani, odstranena chyba se zapornymi pravdepodobnostmi. JS

Files:

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

applications/robust/robustlib.cpp

@@ -55 +55 @@
                         vec cur_condition = as_toprow->condition.get(1,as_toprow->condition.size()-1);
 
+                        // cout << as_toprow->condition << endl;
+
                         vertex* base_vertex = (*simplex.begin());
 
@@ -61 +63 @@
                         mat jacobian(dimension,dimension);
 
-                        double a_0 = base_vertex->get_coordinates()*cur_condition+as_toprow->condition[0];
+                        // cout << base_vertex->get_coordinates() << endl;
+
+                        double a_0 = base_vertex->get_coordinates()*cur_condition-as_toprow->condition[0];
                         map<double,int> as;
 
@@ -69 +73 @@
                                 vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates();
 
+                                // cout << (*vert_ref)->get_coordinates() << endl;
                                 // cout << relative_coords << endl;
 
@@ -80 +85 @@
                                         (*returned.first).second++;
                                 }
+                                /*
                                 else
                                 {
                                         cout << "a[" << row_count << "] = " << new_a << endl;
                                 }
+                                */
+                                
                                 //as.ins(as.size(),new_a);                                                      
                                 
@@ -89 +97 @@
                         }
 
-                        // cout << endl;
+                        
 
                         double int_value = 0;
@@ -95 +103 @@
                         // cout << jacobian << endl;
 
-                        double det_jacobian = det(jacobian);
+                        double det_jacobian    = abs(det(jacobian));
+                        double correction_term = det_jacobian;                  
                         for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++)
                         {
                                 double multiplier = det_jacobian;
+                                
                                 if(a_0!=(*as_ref).first)
-                                {                                                               
+                                {
                                         int as1_order = (*as_ref).second;
-
+                                        
+                                        correction_term /= -pow((*as_ref).first,as1_order);                                     
+                                        
                                         current_emlig->set_correction_factors(as1_order);
 
@@ -126 +138 @@
                                                         multiplier        /= (*as_ref).first;
                                                         number_of_factors += 1;
-                                                }                                                                       
+                                                }
+                                                
                                         }                                       
 
-                                        double bracket = fact(condition_order-number_of_factors)/fact(as1_order-1)/pow((*as_ref).first-a_0,condition_order-number_of_factors);
+                                        double bracket = fact(condition_order-number_of_factors)/fact(as1_order-1)/pow(a_0+(*as_ref).first,condition_order-number_of_factors+1);
                                         for(int k = 0;k < as1_order-1;k++)
                                         {
-                                                double bracket_factor = pow((double)-1,k+1)*fact(condition_order-1-number_of_factors-k)/fact(as1_order-2-k)/pow((*as_ref).first-a_0,condition_order-1-number_of_factors-k);
+                                                double bracket_factor = pow((double)-1,k+1)*fact(condition_order-1-number_of_factors-k)/fact(as1_order-2-k)/pow(a_0+(*as_ref).first,condition_order-number_of_factors-k);
                                                 
                                                 for(set<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[k].begin();combi_ref!=this->my_emlig->correction_factors[k].end();combi_ref++)
@@ -144 +157 @@
                                                         bracket+=bracket_factor*bracket_combination;                                                                    
                                                 }                                                                       
-                                        }                                                               
+                                        }
 
-                                        int_value += multiplier*bracket;        
+                                        
 
+                                        int_value += multiplier*bracket;
                                 }
                                 else
                                 {
                                         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)");
-                                }
+                                }                                                                                                               
+                        }
 
-                                // cout << c << int_value << endl;
+                        
+                        correction_term *= fact(condition_order-my_emlig->number_of_parameters)/pow(a_0,condition_order-my_emlig->number_of_parameters+1);
 
-                                return int_value;                                                                               
-                        }
+                        // cout << c << int_value << endl;
+
+                        int_value += correction_term;
+
+                        // cout << int_value << endl;
+
+                        // cout  << "***************************"  << endl << endl;
+
+                        return int_value;
                         
                 }