root/applications/robust/robustlib.cpp @ 1376

#include "robustlib.h"


double polyhedron::triangulate(bool should_integrate)
        {
                for(set<simplex*>::iterator t_ref = triangulation.begin();t_ref!=triangulation.end();t_ref++)
                {
                        delete (*t_ref);
                }               
                triangulation.clear();

                if(should_integrate)
                {
                        ((toprow *)this)->probability = 0.0;
                }               
                
                if(vertices.size()==1)
                {
                        simplex* vert_simplex = new simplex((vertex*)this);                     

                        triangulation.insert(vert_simplex);
                }

                for(list<polyhedron*>::iterator child_ref = children.begin();child_ref!=children.end();child_ref++)
                {                       
                        for(set<simplex*>::iterator s_ref = (*child_ref)->triangulation.begin();s_ref!=(*child_ref)->triangulation.end();s_ref++)
                        {
                                simplex* new_simplex = new simplex((*s_ref)->vertices);                                                         

                                pair<set<vertex*>::iterator,bool> ret_val = new_simplex->vertices.insert(*vertices.begin());

                                if(ret_val.second == true)
                                {
                                        double cur_prob = 0;

                                        /*
                                        if(should_integrate&&new_simplex->vertices.size()!=3)
                                        {
                                                cout << "Error: Wrong vertex count for integration!";
                                        }
                                        */

                                        if(should_integrate)
                                        {
                                                cur_prob = ((toprow *)this)->integrate_simplex(new_simplex, 'S');

                                                ((toprow *)this)->probability += cur_prob;
                                        }

                                        triangulation.insert(new_simplex);
                                }
                                else
                                {
                                        delete new_simplex;
                                }
                        }       
                }       

                if(should_integrate)
                {
                        return ((toprow *)this)->probability;
                }
                else
                {
                        return 0.0;
                }
                
                /*
                for(list<polyhedron*>::iterator child_ref = children.begin();child_ref!=children.end();child_ref++)
                {                       
                        for(list<set<vertex*>>::iterator t_ref = (*child_ref)->triangulation.begin();t_ref!=(*child_ref)->triangulation.end();t_ref++)
                        {
                                set<vertex*> new_simplex;
                                new_simplex.insert((*t_ref).begin(),(*t_ref).end());                    

                                for(set<vertex*>::iterator suitable_vert_ref = vertices.begin();*suitable_vert_ref<*(*t_ref).begin();suitable_vert_ref++)
                                {
                                        set<vertex*> suitable_simplex;
                                        suitable_simplex.insert(new_simplex.begin(),new_simplex.end());
                                        
                                        suitable_simplex.insert(*suitable_vert_ref);

                                        triangulation.push_back(suitable_simplex);

                                        if(should_integrate)
                                        {
                                                ((toprow *)this)->probability += ((toprow *)this)->integrate_simplex(suitable_simplex, 'S');
                                        }
                                }                               
                        }       
                }*/             
        }


        double toprow::integrate_simplex(simplex* simplex, char c)
        {
                // cout << ((toprow*)this)->condition << endl;          
                int sigma_order = ((toprow*)this)->condition_order-simplex->vertices.size()-1;

                // cout << "C:" << condition_order << "  N:" << my_emlig->number_of_parameters << "  C+N:" << condition_order-my_emlig->number_of_parameters << endl;
                // pause(0.1);

                if(sigma_order >= 0)
                {

                        //cout << endl;
                        //cout << ((toprow*)this)->condition << 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;
                        simplex->clear_gammas();

                        if(this->my_emlig!=NULL)
                        {
                                current_emlig = this->my_emlig;
                        }
                        else
                        {
                                throw exception("The statistic of the polyhedron you are trying to integrate over doesn't belong to any emlig!");
                        }                                               

                        toprow* as_toprow = (toprow*)this;                                              

                        int dimension = simplex->vertices.size()-1;

                        mat jacobian(dimension,dimension);                      

                        map<double,int> as;
                        vertex* base_vertex = simplex->vertices.begin();
                                                        
                                
                        for(set<vertex*>::iterator vert_ref = simplex->vertices.begin(); vert_ref!=simplex->vertices.end();vert_ref++)
                        {
                                if(vert_ref!=simplex->vertices.begin())
                                {
                                        vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates();
                                        jacobian.set_row(row_count,relative_coords);
                                }

                                double a = -((*vert_ref)->get_coordinates().ins(0,-1)*cur_condition);
                                if(a<current_emlig->min_ll)
                                {
                                        current_emlig->minimal_vertex = (*vert_ref);
                                        current_emlig->min_ll = a;                                              
                                }

                                //cout << "a0:" << a_0 << " a0 coords:" << base_vertex->get_coordinates() << " am:" << a_m << " am coords:" << (*vert_ref)->get_coordinates() << 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>(a,1));
                                if(returned.second == false)
                                {
                                        (*returned.first).second++;
                                }                                               
                        }
                        
                        
                        /*
                        cout << "a_0: " << a_0 << "    ";
                        int as_count = 1;
                        for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++)
                        {
                                cout << "a_" << as_count << ": " << (*as_ref).first << "    ";
                                as_count++;
                        }
                        */

                        double int_value = 0;                   

                        // cout << jacobian << endl;

                        double det_jacobian    = abs(det(jacobian));
                        double correction_term;
                        for(map<double,int>::iterator a_ref = as.begin();as_ref!=as.end();as_ref++)
                        {
                                double multiplier = det_jacobian/fact(jacobian.rows());                         
                                
                                int a_order = (*a_ref).second;                  
                                current_emlig->set_correction_factors(a_order);

                                vector<double> factors;
                                int number_of_factors = 0;                                                              
                                for(map<double,int>::iterator a2_ref = as.begin();a2_ref!=as.end();a2_ref++)
                                {                                       
                                        if(a2_ref!=a_ref)
                                        {                                                                               
                                                for(int k = 0;k<(*a2_ref).second;k++)
                                                {
                                                        factors.push_back((*a_ref).first-(*a2_ref).first);
                                                }

                                                multiplier        /= pow((*a_ref).first-(*a2_ref).first,(*a2_ref).second);
                                                number_of_factors += (*a2_ref).second;
                                        }                                                                               
                                }       

                                double cur_as_ref = (*a_ref).first;                             

                                double gamma_multiplier = -cur_as_ref-a_0;

                                double bracket = fact(as1_order-1)/pow(gamma_multiplier,sigma_order);
                                                                
                                simplex->insert_gamma(0,bracket*multiplier,gamma_multiplier);                                                           

                                // bracket *= fact(sigma_order);                                

                                for(int k = 0;k < as1_order-1;k++)
                                {
                                        double local_bracket = 0;
                                        double bracket_factor = 1/fact(as1_order-1-k)/pow(gamma_multiplier,sigma_order-k-1);//pow((double)-1,k+1)
                                        
                                        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++)
                                                {
                                                        //cout << "Factor vector:" << (*combi_ref) << endl;
                                                        
                                                        bracket_combination /= factors[(*combi_ref)[j]-1];
                                                }                                               
                                                                                                
                                                local_bracket += bracket_factor*bracket_combination;                                                                    
                                        }

                                        simplex->insert_gamma(k+1,local_bracket*multiplier,gamma_multiplier);

                                        int division_factor = 1;
                                        for(int s = 0;s<k;s++)
                                        {
                                                division_factor *= k-s;
                                        }
                                        
                                        bracket += local_bracket/division_factor; //local_bracket*fact(sigma_order-k);
                                }                               

                                int_value += multiplier*bracket;

                                
                                                                                                                                                
                        }

                        double correction_term_base = correction_term/pow(-a_0,sigma_order);

                        simplex->insert_gamma(0,correction_term_base,-a_0);                     

                        correction_term = correction_term_base;//fact(sigma_order)*correction_term_base;

                        //cout << c << int_value << endl;

                        int_value += correction_term;                   

                        //cout << "Probability:" << int_value << endl;
                        //pause(0.100);

                        simplex->probability = int_value;
                        
                        return int_value;       
                        
                }
                else
                {
                        cout << "Improper probability density." << endl;
                        return 0.0;
                }
        }