Changeset 1280 for applications/robust

Timestamp:

02/18/11 21:50:33 (13 years ago)

Author:

sindj

Message:

Oprava chyb pri vypoctech integralu - nedokonceno, ale blizim se. JS

Location:

applications/robust

Files:

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

applications/robust/main.cpp

@@ -39 +39 @@
         }*/
         
-        vec condition1 = "0.3 1.0 1.0";
+        vec condition1 = "1.0 1.0 1.0";
         emlig1->add_condition(condition1);
 
-        vec condition2 = "-0.5 -1.0 1.0";
+        vec condition2 = "-1.0 1.0 1.0";
         emlig1->add_condition(condition2);
+
+        vec condition3 = "0.5 -1.1 1.0";
+        emlig1->add_condition(condition3);
+
+        vec condition4 = "-0.5 -1.0 1.0";
+        emlig1->add_condition(condition4);
 
         //vec condition3 = "-0.3 0.5 0.5";

applications/robust/robustlib.cpp

@@ -35 +35 @@
                 // cout << ((toprow*)this)->condition << endl;
                 
-                int condition_order = ((toprow*)this)->condition_order+1; // -2 by bylo, pokud chceme uniformni apriorno
+                int condition_order = ((toprow*)this)->condition_order-2; // -2 by bylo, pokud chceme uniformni apriorno
+
+                int sigma_order = condition_order-my_emlig->number_of_parameters;
 
                 // 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(sigma_order >= 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;
+                        cout << "C:" << condition_order+2 << "  N:" << my_emlig->number_of_parameters << "  C+N:" << condition_order-my_emlig->number_of_parameters+2 << endl;
 
                         emlig* current_emlig;
@@ -111 +113 @@
                                                 }
                                                 
-                                                // cout << "Absolute coords:(V"  << row_count << ")" << (*vert_ref)->get_coordinates() << endl;
-                                                cout << "Relative coords:(V"  << row_count << ")" << relative_coords << endl;
+                                                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));
@@ -147 +149 @@
                                 int as1_order = (*as_ref).second;
                                 
-                                correction_term /= -pow((*as_ref).first,as1_order);                                     
+                                correction_term /= pow(-(*as_ref).first,as1_order);                                     
                                 
                                 current_emlig->set_correction_factors(as1_order);
@@ -174 +176 @@
                                         }
                                         
-                                }                                       
-
-                                double bracket = fact(condition_order-my_emlig->number_of_parameters)/fact(as1_order-1)/pow(a_0-(*as_ref).first,condition_order-my_emlig->number_of_parameters+2);
+                                }       
+
+                                
+
+                                double bracket = fact(sigma_order)/fact(as1_order-1)/pow(a_0-(*as_ref).first,sigma_order+2);
                                 for(int k = 0;k < as1_order-1;k++)
                                 {
-                                        double bracket_factor = pow((double)-1,k+1)*fact(condition_order-my_emlig->number_of_parameters-k)/fact(as1_order-1-k)/pow(a_0-(*as_ref).first,condition_order-my_emlig->number_of_parameters+1-k);
-                                        
-                                        // TODO TADY NEKDE JE CHYBA, NEDOJDE KE SPRAVNEMU NAPLNENI CORRECTION FAKTORU!!!
-                                        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_factor = pow((double)-1,k+1)*fact(sigma_order-k)/fact(as1_order-1-k)/pow(a_0-(*as_ref).first,sigma_order+2-k);
+                                        
+                                        ivec control_vec = ivec();
+                                        control_vec.ins(0,my_emlig->number_of_parameters-as1_order+1);          
+                                        
+                                        
+                                        for(multiset<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[k].begin();combi_ref!=this->my_emlig->correction_factors[k].upper_bound(my_ivec(control_vec));combi_ref++)
                                         {
                                                 double bracket_combination = 1;
-                                                for(int j = 0;j<=(*combi_ref).size();j++)
-                                                {
-                                                        bracket_combination /= factors[(*combi_ref)[j]];
+                                                for(int j = 0;j<(*combi_ref).size();j++)
+                                                {
+                                                        cout << "Factor vector:" << (*combi_ref) << endl;
+                                                        
+                                                        bracket_combination /= factors[(*combi_ref)[j]-1];
                                                 }
 
@@ -201 +210 @@
 
                         
-                        correction_term *= fact(condition_order-my_emlig->number_of_parameters)/pow(a_0,condition_order-my_emlig->number_of_parameters+2);
+                        correction_term *= fact(sigma_order)/pow(a_0,sigma_order+2);
 
                         // cout << c << int_value << endl;

applications/robust/robustlib.h

@@ -21 +21 @@
 using namespace itpp;
 
-const double max_range = 1000.0;//numeric_limits<double>::max()/10e-10;
+const double max_range = 10.0;//numeric_limits<double>::max()/10e-10;
 
 enum actions {MERGE, SPLIT};
@@ -435 +435 @@
         bool operator>(const my_ivec &second) const
         {
+                return max(*this)>max(second);
+                
+                /*
                 int size1 = this->size();
-                int size2 = second.size();
+                int size2 = second.size();              
                  
                 int counter1 = 0;
@@ -481 +484 @@
 
                         return false;
-                }
+                }*/
         }
 
@@ -487 +490 @@
         bool operator==(const my_ivec &second) const
         {
+                return max(*this)==max(second);
+                
+                /*
                 int size1 = this->size();
-                int size2 = second.size();
+                int size2 = second.size();              
                  
                 int counter = 0;
@@ -533 +539 @@
 
                         return true;
-                }
+                }*/
         }
 
@@ -796 +802 @@
         c_statistic statistic;
 
-        vector<set<my_ivec>> correction_factors;
+        vector<multiset<my_ivec>> correction_factors;
 
         int number_of_parameters;
@@ -1121 +1127 @@
         void set_correction_factors(int order)
                 {
-                        for(int remaining_order = correction_factors.size();!(remaining_order>order);remaining_order++)
-                        {
-                                set<my_ivec> factor_templates;
-                                set<my_ivec> final_factors;
-
-                                
-                                for(int i = 1;i!=number_of_parameters-order+1;i++)
+                        for(int remaining_order = correction_factors.size();remaining_order<order;remaining_order++)
+                        {
+                                multiset<my_ivec> factor_templates;
+                                multiset<my_ivec> final_factors;                                
+
+                                my_ivec orig_template = my_ivec();                              
+
+                                for(int i = 1;i<number_of_parameters-remaining_order+1;i++)
                                 {                                       
-                                        my_ivec new_template = my_ivec();
-                                        new_template.ins(0,1);
-                                        new_template.ins(1,i);
-                                        factor_templates.insert(new_template);
-
-                                        
-                                        for(int j = 1;j<remaining_order;j++)
-                                        {
+                                        bool in_cycle = false;
+                                        for(int j = 0;j<=remaining_order;j++)                                   {
                                                 
-                                                for(set<my_ivec>::iterator fac_ref = factor_templates.begin();fac_ref!=factor_templates.end();fac_ref++)
+                                                multiset<my_ivec>::iterator fac_ref = factor_templates.begin();
+
+                                                do
                                                 {
-                                                        ivec current_template = (*fac_ref);
+                                                        my_ivec current_template;
+                                                        if(!in_cycle)
+                                                        {
+                                                                current_template = orig_template;
+                                                                in_cycle = true;
+                                                        }
+                                                        else
+                                                        {
+                                                                current_template = (*fac_ref);
+                                                                fac_ref++;
+                                                        }                                                       
                                                         
-                                                        current_template[0]+=1;
                                                         current_template.ins(current_template.size(),i);
 
+                                                        cout << "template:" << current_template << endl;
                                                         
-                                                        if(current_template[0]==remaining_order)
+                                                        if(current_template.size()==remaining_order+1)
                                                         {
-                                                                final_factors.insert(current_template.right(current_template.size()-1));
+                                                                final_factors.insert(current_template);
                                                         }
                                                         else
@@ -1155 +1168 @@
                                                         }
                                                 }
+                                                while(fac_ref!=factor_templates.end());
                                         }
                                 }