root/applications/robust/robustlib.cpp @ 1365

#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;

                        vec cur_condition = as_toprow->condition_sum.get(1,as_toprow->condition_sum.size()-1);

                        // cout << as_toprow->condition << endl;                        

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

                        mat jacobian(dimension,dimension);                      

                        double a_0;
                        map<double,int> as;
                        vertex* base_vertex;
                        set<vertex*>::iterator base_ref = simplex->vertices.begin();
                        bool order_correct;
                                                

                        do
                        {
                                order_correct = true;
                                as.clear();

                                base_vertex = (*base_ref);
                        
                                //cout << "Base coords: " << base_vertex->get_coordinates() << endl;
                                //cout << "Condition: " << as_toprow->condition_sum << endl;
                                //pause(0.1);

                                a_0 = -base_vertex->get_coordinates()*cur_condition+as_toprow->condition_sum[0];                
                                
                                int row_count = 0;
                                for(set<vertex*>::iterator vert_ref = simplex->vertices.begin(); vert_ref!=simplex->vertices.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 || new_a == 0)
                                                {
                                                        base_ref++;

                                                        if(base_ref == simplex->vertices.end()&&simplex->vertices.size()!=2)
                                                        {
                                                                throw new exception("Equal local conditions are paired. If this ever occurs, the software has to be modified to include multiplied a_0!!");
                                                        }
                                                        /*
                                                        else if(implex->vertices.size()==2)
                                                        {
                                                                return relative_coords[0]*exp
                                                        }
                                                        */
                                                        
                                                        order_correct = false;                                          

                                                        break;
                                                }                                               
                                                
                                                

                                                if(a_0<current_emlig->min_ll)
                                                {
                                                        current_emlig->minimal_vertex = base_vertex;
                                                        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->min_ll)
                                                {
                                                        current_emlig->minimal_vertex = (*vert_ref);
                                                        current_emlig->min_ll = a_m;                                            
                                                }

                                                //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>(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++;
                                        }
                                }
                        }
                        while(!order_correct);
                        
                        /*
                        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 = det_jacobian;                  
                        for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++)
                        {
                                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++)
                                {
                                        
                                        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 gamma_multiplier = -(*as_ref).first-a_0;

                                double bracket = fact(as1_order-1)/pow(gamma_multiplier,sigma_order+1);
                                                                
                                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);//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+1);

                        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;
                }
        }