Changeset 1335 for applications/robust

Timestamp:

04/19/11 18:15:22 (14 years ago)

Author:

sindj

Message:

Dokoncovani samplingu, uz to skoro funguje, jen to spatne generuje sigma. JS

Location:

applications/robust

Files:

• robustlib.cpp (modified) (4 diffs)
• robustlib.h (modified) (26 diffs)

applications/robust/robustlib.cpp

@@ -2 +2 @@
 
 
-void polyhedron::triangulate(bool should_integrate)
+double polyhedron::triangulate(bool should_integrate)
         {
                 for(set<simplex*>::iterator t_ref = triangulation.begin();t_ref!=triangulation.end();t_ref++)
@@ -48 +48 @@
                                 }
                         }       
-                }               
+                }       
+
+                if(should_integrate)
+                {
+                        return ((toprow *)this)->probability;
+                }
+                else
+                {
+                        return 0.0;
+                }
                 
                 /*
@@ -165 +174 @@
                                                 
 
-                                                if(a_0<current_emlig->likelihood_value)
+                                                if(a_0<current_emlig->min_ll)
                                                 {
                                                         current_emlig->minimal_vertex = base_vertex;
-                                                        current_emlig->likelihood_value = a_0;
+                                                        current_emlig->min_ll = a_0;
                                                 }
 
                                                 double a_m = (*vert_ref)->get_coordinates()*cur_condition-as_toprow->condition_sum[0];
-                                                if(a_m<current_emlig->likelihood_value)
+                                                if(a_m<current_emlig->min_ll)
                                                 {
                                                         current_emlig->minimal_vertex = (*vert_ref);
-                                                        current_emlig->likelihood_value = a_m;                                          
+                                                        current_emlig->min_ll = a_m;                                            
                                                 }
 
@@ -212 +221 @@
                         }*/
 
-                        double int_value = 0;
+                        double int_value = 0;                   
 
                         // cout << jacobian << endl;

applications/robust/robustlib.h

@@ -82 +82 @@
 
         double negative_gamma_sum;
+
+        double min_beta;
         
 
@@ -99 +101 @@
         {
                 positive_gamma_parameters.clear();
-                negative_gamma_parameters.clear();
-                
+                negative_gamma_parameters.clear();              
                 
                 positive_gamma_sum = 0;
                 negative_gamma_sum = 0;
+
+                min_beta = numeric_limits<double>::max();
         }
 
@@ -131 +134 @@
 
                         negative_gamma_parameters[order].insert(pair<double,double>(-weight,beta));
+                }
+
+                if(beta < min_beta)
+                {
+                        min_beta = beta;
                 }
         }
@@ -321 +329 @@
 
         /// A method for triangulation of given polyhedron.
-        void triangulate(bool should_integrate);        
+        double triangulate(bool should_integrate);      
 };
 
@@ -854 +862 @@
                                                                         //new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second));
                                                                 }
+
+                                                                cur_par_toprow->my_emlig->normalization_factor += cur_par_toprow->probability;
 
                                                                 //current_parent->triangulation.clear();
@@ -959 +969 @@
                                                                         //new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second));
                                                                 }
+
+                                                                cur_par_toprow->my_emlig->normalization_factor += cur_par_toprow->probability;
 
                                                                 //current_parent->triangulation.clear();
@@ -999 +1011 @@
         vertex* minimal_vertex;
 
-        double likelihood_value;
+        double min_ll;
 
         vector<multiset<my_ivec>> correction_factors;
@@ -1013 +1025 @@
                 create_statistic(number_of_parameters);
 
-                likelihood_value = numeric_limits<double>::max();
+                min_ll = numeric_limits<double>::max();
         }
 
@@ -1022 +1034 @@
                 this->statistic = statistic;    
 
-                likelihood_value = numeric_limits<double>::max();
+                min_ll = numeric_limits<double>::max();
         }
 
@@ -1146 +1158 @@
         {
                 //step_me(0);
-                
-                likelihood_value = numeric_limits<double>::max();
+                normalization_factor = 0;
+                min_ll = numeric_limits<double>::max();
 
                 bool should_remove = (toremove.size() != 0);
@@ -1480 +1492 @@
                                                 current_positive->totally_neutral = (current_positive->totally_neutral && current_negative->totally_neutral);
 
-                                                current_positive->triangulate(k==for_splitting.size()-1);
+                                                current_positive->my_emlig->normalization_factor += current_positive->triangulate(k==for_splitting.size()-1);
                                                 
                                                 statistic.delete_polyhedron(k,current_negative);
@@ -1693 +1705 @@
                                         new_totally_neutral_child->triangulate(false);
 
-                                        positive_poly->triangulate(k==for_splitting.size()-1);
-                                        negative_poly->triangulate(k==for_splitting.size()-1);
+                                        positive_poly->my_emlig->normalization_factor += positive_poly->triangulate(k==for_splitting.size()-1);
+                                        negative_poly->my_emlig->normalization_factor += negative_poly->triangulate(k==for_splitting.size()-1);
                                         
                                         statistic.append_polyhedron(k-1, new_totally_neutral_child);                                    
@@ -1790 +1802 @@
                 /// \TODO tady je to spatne, tady nesmi byt conditions.size(), viz RARX.bayes()
                 if(conditions.size()-2-number_of_parameters>=0)
-                {
-                        
-                        
+                {                       
                         mat sample_mat;
                         map<double,toprow*> ordered_toprows;                    
                         double sum_a = 0;
+                        
+                        //cout << "Likelihoods of toprows:" << endl;
 
                         for(polyhedron* top_ref = statistic.rows[number_of_parameters];top_ref!=statistic.end_poly;top_ref=top_ref->next_poly)
@@ -1802 +1814 @@
 
                                 sum_a+=current_top->probability;
+                                /*
+                                cout << current_top->probability << "   ";
+
+                                for(set<vertex*>::iterator vert_ref = (*top_ref).vertices.begin();vert_ref!=(*top_ref).vertices.end();vert_ref++)
+                                {
+                                        cout << round(100*(*vert_ref)->get_coordinates())/100 << " ; ";
+                                }
+                                */
+
+                                // cout << endl;
                                 ordered_toprows.insert(pair<double,toprow*>(sum_a,current_top));
-                        }
+                        }                       
+                        
+                        // cout << "Sum N: " << normalization_factor << endl;
 
                         while(sample_mat.cols()<n)
@@ -1847 +1871 @@
 
                                 cout << "Simplex/Count: " << simplex_count << "/" << sampled_toprow->triangulation.size() << endl;
+                                cout << "Simplex factor: " << (*s_ref)->probability << endl;
+                                cout << "Toprow factor:  " << sampled_toprow->probability << endl;
+                                cout << "Emlig factor:   " << normalization_factor << endl;
                                 // cout << &(*tri_ref) << endl;
 
+                                int number_of_runs = 0;
                                 bool have_sigma = false;
                                 double sigma = 0;
@@ -1862 +1890 @@
                                                         sum_g += (*g_ref).first/(*s_ref)->positive_gamma_sum;
 
+                                                        
                                                         if(sum_g>rnumber)
                                                         {
                                                                 itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,1/(*g_ref).second);
                                                                 sigma = 1/(*gamma)();
+                                                                
+                                                                cout << "Sigma mean:   " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl;                                                         
                                                                 break;
-                                                        }                                               
+                                                        }                                                       
                                                 }
 
@@ -1883 +1914 @@
                                                 for(multimap<double,double>::iterator pg_ref = (*v_ref).begin();pg_ref!=(*v_ref).end();pg_ref++)
                                                 {
-                                                        pg_sum += exp(-(*pg_ref).second/sigma)*pow((*pg_ref).second/sigma,(int)conditions.size()-number_of_parameters-1)*(*pg_ref).second/fact(conditions.size()-number_of_parameters-1)*(*pg_ref).first;
+                                                        pg_sum += exp(((*s_ref)->min_beta-(*pg_ref).second)/sigma)*pow((*pg_ref).second/sigma,(int)conditions.size()-number_of_parameters-1)*(*pg_ref).second/fact(conditions.size()-number_of_parameters-1)*(*pg_ref).first;
                                                 }                                       
                                         }
@@ -1892 +1923 @@
                                                 for(multimap<double,double>::iterator ng_ref = (*v_ref).begin();ng_ref!=(*v_ref).end();ng_ref++)
                                                 {
-                                                        ng_sum += exp(-(*ng_ref).second/sigma)*pow((*ng_ref).second/sigma,(int)conditions.size()-number_of_parameters-1)*(*ng_ref).second/fact(conditions.size()-number_of_parameters-1)*(*ng_ref).first;
+                                                        ng_sum += exp(((*s_ref)->min_beta-(*ng_ref).second)/sigma)*pow((*ng_ref).second/sigma,(int)conditions.size()-number_of_parameters-1)*(*ng_ref).second/fact(conditions.size()-number_of_parameters-1)*(*ng_ref).first;
                                                 }                                       
                                         }
@@ -1900 +1931 @@
                                                 have_sigma = true;
                                         }
+
+                                        number_of_runs++;
                                 }
                                 while(!have_sigma);
+
+                                cout << "Sigma: " << sigma << endl;
+                                cout << "Nr. of runs: " << number_of_runs << endl;
 
                                 int dimension = (*s_ref)->vertices.size()-1;
@@ -1909 +1945 @@
 
                                 vertex* base_vert = *(*s_ref)->vertices.begin();
+
+                                cout << "Base vertex coords(should be close to est. param.): " << base_vert->get_coordinates() << endl;
                                 
                                 int row_count = 0;
@@ -1914 +1952 @@
                                 for(set<vertex*>::iterator vert_ref = ++(*s_ref)->vertices.begin();vert_ref!=(*s_ref)->vertices.end();vert_ref++)
                                 {
-                                        vec relative_coords = (*vert_ref)->get_coordinates()-base_vert->get_coordinates();                              
+                                        vec current_coords = (*vert_ref)->get_coordinates();
+
+                                        cout << "Coords of vertex[" << row_count << "]: " << current_coords << endl; 
+                                        
+                                        vec relative_coords = current_coords-base_vert->get_coordinates();                              
 
                                         jacobian.set_row(row_count,relative_coords);
 
                                         row_count++;
-                                }
+                                }                               
                                 
                                 cout << "Jacobian: " << jacobian << endl;
@@ -1926 +1968 @@
                                                                 
                                 gradient = jacobian*gradient;   
+
+                                cout << "Gradient after trafo:" << gradient << endl;
 
                                 // vec normal_gradient = gradient/sqrt(gradient*gradient);
@@ -1943 +1987 @@
                                         x_axis = rotation_matrix*x_axis;
 
-                                        double t = gradient[i]/gradient*x_axis;
+                                        double t = abs(gradient[i]/gradient*x_axis);
 
                                         double sin_theta = sign(gradient[i])*t/sqrt(1+pow(t,2));
@@ -1961 +2005 @@
 
                                 // cout << rotation_matrix << endl;
-                                cout << "Gradient after trafo:" << gradient << endl;
-
-                                vec minima = numeric_limits<double>::max()*ones(number_of_parameters);
-                                vec maxima = -numeric_limits<double>::max()*ones(number_of_parameters);
-
-                                vertex* min_grad;
-                                vertex* max_grad;
-                                                                
-                                for(set<vertex*>::iterator vert_ref = (*s_ref)->vertices.begin();vert_ref!=(*s_ref)->vertices.end();vert_ref++)
-                                {
-                                        vec new_coords = rotation_matrix*jacobian*(*vert_ref)->get_coordinates();
-
-                                        // cout << new_coords << endl;
-
-                                        for(int j = 0;j<new_coords.size();j++)
-                                        {
-                                                if(new_coords[j]>maxima[j])
-                                                {
-                                                        maxima[j] = new_coords[j];
-
-                                                        if(j==0)
-                                                        {
-                                                                max_grad = (*vert_ref);
-                                                        }
-                                                }
-
-                                                if(new_coords[j]<minima[j])
-                                                {
-                                                        minima[j] = new_coords[j];
-
-                                                        if(j==0)
-                                                        {
-                                                                min_grad = (*vert_ref);
-                                                        }
-                                                }
-                                        }                                       
-                                }
-
-                                vec sample_coordinates;                         
+                                
+                                mat extended_rotation = rotation_matrix;
+                                extended_rotation.ins_col(0,zeros(extended_rotation.rows()));
+
+                                cout << "Extended rotation: " << extended_rotation << endl;
+                                
+                                vec minima = itpp::min(extended_rotation,2);
+                                vec maxima = itpp::max(extended_rotation,2);
+
+                                cout << "Minima: " << minima << endl;
+                                cout << "Maxima: " << maxima << endl;
+
+                                vec sample_coordinates;         
+                                bool is_inside = true;
+                                
+                                vec new_sample;
+                                sample_coordinates = new_sample;
 
                                 for(int j = 0;j<number_of_parameters;j++)
@@ -2008 +2030 @@
 
                                         if(j==0)
-                                        {
-                                                double pos_value = sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0];
+                                        {                                               
+                                                vec new_gradient = rotation_matrix*gradient;
                                                 
-                                                vec new_gradient = rotation_matrix*gradient;
-
                                                 cout << "New gradient(should have only first component nonzero):" << new_gradient << endl;
 
-                                                cout << "Max: " << maxima[0] << "  Min: " << minima[0] << "  Grad:" << new_gradient[0] << endl;
+                                                // cout << "Max: " << maxima[0] << "  Min: " << minima[0] << "  Grad:" << new_gradient[0] << endl;
                                                 
-                                                double log_bracket = rnumber*exp(new_gradient[0]*maxima[0]/sigma)+(1-rnumber)*exp(new_gradient[0]*minima[0]/sigma);
+                                                double log_bracket = 1-rnumber*(1-exp(new_gradient[0]/sigma*(minima[0]-maxima[0])));
                                                 
-                                                coordinate = log(log_bracket);
-                                        }
+                                                coordinate = minima[0]-sigma/new_gradient[0]*log(log_bracket);
+                                        }
+                                        else
+                                        {
+                                                coordinate = minima[j]+rnumber*(maxima[j]-minima[j]);
+                                        }
+
+                                        sample_coordinates.ins(j,coordinate);
+                                }
+
+                                cout << "Sampled coordinates(gradient direction): " << sample_coordinates << endl;
+
+                                sample_coordinates = rotation_matrix.T()*sample_coordinates;
+
+                                cout << "Sampled coordinates(backrotated direction):" << sample_coordinates << endl;
+
+                                
+                                for(int j = 0;j<sample_coordinates.size();j++)
+                                {
+                                        if(sample_coordinates[j]<0)
+                                        {
+                                                is_inside = false;
+                                        }
+                                }
+
+                                double above_criterion = ones(sample_coordinates.size())*sample_coordinates;
+
+                                if(above_criterion>1)
+                                {
+                                        is_inside = false;
+                                }
+
+                                if(is_inside)
+                                {                                       
+                                        sample_coordinates = jacobian.T()*sample_coordinates+(*base_vert).get_coordinates();
+                                        
+                                        sample_coordinates.ins(0,sigma);
+                                        
+                                        cout << "Sampled coordinates(parameter space):" << sample_coordinates << endl;
+
+                                        sample_mat.ins_col(0,sample_coordinates);
                                 }
 
@@ -2026 +2085 @@
                                 // cout << sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0] << endl;
 
+                                
                         }