Changeset 1275 for applications/robust

Timestamp:

02/16/11 20:38:48 (13 years ago)

Author:

sindj

Message:

Pokracovani spojovani polyhedronu, ladeni vypoctu. JS

Location:

applications/robust

Files:

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

applications/robust/main.cpp

@@ -11 +11 @@
 using namespace itpp;
 
-using namespace bdm;
+//using namespace bdm;
 
-const int emlig_size = 4;
+const int emlig_size = 2;
 
 
@@ -39 +39 @@
         }*/
         
+        vec condition1 = "0.3 1.0 1.0";
+        emlig1->add_condition(condition1);
 
-        //emlig1->step_me(0);
+        vec condition2 = "-0.5 -1.0 1.0";
+        emlig1->add_condition(condition2);
 
+        //vec condition3 = "-0.3 0.5 0.5";
+        //emlig1->add_condition(condition3);
+
+        // emlig1->step_me(0);
+
+        //emlig1->remove_condition(condition1);
+
+        
+
+        
+
+        /*
         for(int i = 0;i<500;i++)
         {
@@ -65 +80 @@
                 }
                 */
-
+                /*
                 cout << emlig1->statistic_rowsize(emlig_size) << endl << endl;
         
@@ -84 +99 @@
                 }
                 */
-        }
+        //}
+    
 
 

applications/robust/robustlib.cpp

@@ -33 +33 @@
         double toprow::integrate_simplex(set<vertex*> simplex, char c)
         {
-                int condition_order = ((toprow*)this)->condition_order-1;
+                // cout << ((toprow*)this)->condition << endl;
+                
+                int condition_order = ((toprow*)this)->condition_order+1; // -2 by bylo, pokud chceme uniformni apriorno
 
                 // cout << "C:" << condition_order << "  N:" << my_emlig->number_of_parameters << "  C+N:" << condition_order-my_emlig->number_of_parameters << endl;
                 // pause(0.1);
 
-                if(condition_order-my_emlig->number_of_parameters >= 0)
+                if(condition_order-my_emlig->number_of_parameters > 0)
                 {
+
+                        cout << endl;
+                        cout << ((toprow*)this)->condition << endl;
+                        cout << "C:" << condition_order << "  N:" << my_emlig->number_of_parameters << "  C+N:" << condition_order-my_emlig->number_of_parameters << endl;
+
                         emlig* current_emlig;
 
@@ -55 +62 @@
                         vec cur_condition = as_toprow->condition.get(1,as_toprow->condition.size()-1);
 
-                        // cout << as_toprow->condition << endl;
-
-                        vertex* base_vertex = (*simplex.begin());
+                        // cout << as_toprow->condition << endl;                        
 
                         int dimension = simplex.size()-1;
 
-                        mat jacobian(dimension,dimension);
-
-                        // cout << base_vertex->get_coordinates() << endl;
-
-                        double a_0 = base_vertex->get_coordinates()*cur_condition-as_toprow->condition[0];
+                        mat jacobian(dimension,dimension);                      
+
+                        double a_0;
                         map<double,int> as;
-
-                        int row_count = 0;
-                        for(set<vertex*>::iterator vert_ref = (++simplex.begin()); vert_ref!=simplex.end();vert_ref++)
-                        {
-                                vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates();
-
-                                // cout << (*vert_ref)->get_coordinates() << endl;
-                                // cout << relative_coords << endl;
-
-                                jacobian.set_row(row_count,relative_coords);
-
-                                double new_a = relative_coords*cur_condition;                                                   
-                                
-                                pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(new_a,1));
-                                if(returned.second == false)
-                                {
-                                        (*returned.first).second++;
+                        vertex* base_vertex;
+                        set<vertex*>::iterator base_ref = simplex.begin();
+                        bool order_correct;
+                                                
+
+                        do
+                        {
+                                order_correct = true;
+                                as.clear();
+
+                                base_vertex = (*base_ref);
+
+                                cout << "Base coords:" << base_vertex->get_coordinates() << endl;
+
+                                a_0 = base_vertex->get_coordinates()*cur_condition-as_toprow->condition[0];
+                                
+                                int row_count = 0;
+                                for(set<vertex*>::iterator vert_ref = simplex.begin(); vert_ref!=simplex.end();vert_ref++)
+                                {
+                                        if(vert_ref != base_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;
+
+                                                if(new_a + a_0 == 0)
+                                                {
+                                                        base_ref++;
+
+                                                        if(base_ref == simplex.end())
+                                                        {
+                                                                throw new exception("Equal local conditions are paired. If this ever occurs, the software has to be modified to include multiplied a_0!!");
+                                                        }
+                                                        
+                                                        order_correct = false;                                          
+
+                                                        break;
+                                                }
+                                                
+                                                cout << "Absolute coords:(V"  << row_count << ")" << (*vert_ref)->get_coordinates() << endl;
+                                                //cout << "Relative coords:(V"  << row_count << ")" << relative_coords << endl;
+
+                                                pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(new_a,1));
+                                                if(returned.second == false)
+                                                {
+                                                        (*returned.first).second++;
+                                                }
+                                                /*
+                                                else
+                                                {
+                                                        cout << "a[" << row_count << "] = " << new_a << endl;
+                                                }
+                                                */
+                                                
+                                                //as.ins(as.size(),new_a);                                                      
+                                                
+                                                row_count++;
+                                        }
                                 }
-                                /*
-                                else
-                                {
-                                        cout << "a[" << row_count << "] = " << new_a << endl;
-                                }
-                                */
-                                
-                                //as.ins(as.size(),new_a);                                                      
-                                
-                                row_count++;
                         }
-
+                        while(!order_correct);
                         
 
@@ -107 +145 @@
                         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;
+                                double multiplier = det_jacobian;                               
+                                
+                                int as1_order = (*as_ref).second;
+                                
+                                correction_term /= -pow((*as_ref).first,as1_order);                                     
+                                
+                                current_emlig->set_correction_factors(as1_order);
+
+                                vector<double> factors;
+                                int number_of_factors = 0;                                                              
+                                for(map<double,int>::iterator as2_ref = as.begin();as2_ref!=as.end();as2_ref++)
+                                {
                                         
-                                        correction_term /= -pow((*as_ref).first,as1_order);                                     
+                                        if(as2_ref!=as_ref)
+                                        {                                                                               
+                                                for(int k = 0;k<(*as2_ref).second;k++)
+                                                {
+                                                        factors.push_back((*as_ref).first-(*as2_ref).first);
+                                                }
+
+                                                multiplier        /= pow((*as_ref).first-(*as2_ref).first,(*as2_ref).second);
+                                                number_of_factors += (*as2_ref).second;
+                                        }
+                                        else
+                                        {
+                                                factors.push_back((*as_ref).first);
+
+                                                multiplier        /= (*as_ref).first;
+                                                number_of_factors += 1;
+                                        }
                                         
-                                        current_emlig->set_correction_factors(as1_order);
-
-                                        vector<double> factors;
-                                        int number_of_factors = 0;                                                              
-                                        for(map<double,int>::iterator as2_ref = as.begin();as2_ref!=as.end();as2_ref++)
-                                        {
-                                                
-                                                if(as2_ref!=as_ref)
-                                                {                                                                               
-                                                        for(int k = 0;k<(*as2_ref).second;k++)
-                                                        {
-                                                                factors.push_back((*as_ref).first-(*as2_ref).first);
-                                                        }
-
-                                                        multiplier        /= pow((*as_ref).first-(*as2_ref).first,(*as2_ref).second);
-                                                        number_of_factors += (*as2_ref).second;
-                                                }
-                                                else
-                                                {
-                                                        factors.push_back((*as_ref).first);
-
-                                                        multiplier        /= (*as_ref).first;
-                                                        number_of_factors += 1;
-                                                }
-                                                
-                                        }                                       
-
-                                        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(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++)
-                                                {
-                                                        double bracket_combination = 1;
-                                                        for(int j = 0;j<(*combi_ref).size();j++)
-                                                        {
-                                                                bracket_combination /= factors[(*combi_ref)[j]];
-                                                        }
-
-                                                        bracket+=bracket_factor*bracket_combination;                                                                    
-                                                }                                                                       
-                                        }
-
+                                }                                       
+
+                                double bracket = fact(condition_order-number_of_factors)/fact(as1_order-1)/pow(a_0+(*as_ref).first,condition_order-number_of_factors);
+                                for(int k = 0;k < as1_order-1;k++)
+                                {
+                                        double bracket_factor = pow((double)-1,k+1)*fact(condition_order-number_of_factors-k)/fact(as1_order-2-k)/pow(a_0+(*as_ref).first,condition_order-number_of_factors-k);
                                         
-
-                                        int_value += multiplier*bracket;
+                                        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++)
+                                        {
+                                                double bracket_combination = 1;
+                                                for(int j = 0;j<(*combi_ref).size();j++)
+                                                {
+                                                        bracket_combination /= factors[(*combi_ref)[j]];
+                                                }
+
+                                                bracket+=bracket_factor*bracket_combination;                                                                    
+                                        }                                                                       
                                 }
-                                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)");
-                                }                                                                                                               
+
+                                
+
+                                int_value += multiplier*bracket;
+                                                                                                                                                
                         }
 
@@ -175 +207 @@
 
                         int_value += correction_term;
-
-                        // cout << int_value << endl;
-
-                        // cout  << "***************************"  << endl << endl;
-
+                
+
+                        cout << "Probability:" << int_value << endl;
+                        
                         return int_value;
+
+                        
                         
                 }

applications/robust/robustlib.h

@@ -822 +822 @@
                         for(polyhedron* horiz_ref = statistic.rows[i];horiz_ref!=statistic.get_end();horiz_ref=horiz_ref->next_poly)
                         {
+                                if(i==statistic.size()-1)
+                                {
+                                        cout << ((toprow*)horiz_ref)->condition << "   " << ((toprow*)horiz_ref)->probability << endl;
+                                        cout << "Order:" << ((toprow*)horiz_ref)->condition_order << endl;
+                                }
                                 char* string = "Checkpoint";
                         }