root/applications/robust/robustlib.cpp @ 1257

#include "robustlib.h"

void polyhedron::triangulate(bool should_integrate)
        {               
                for(list<polyhedron*>::iterator child_ref = children.begin();child_ref!=children.end();child_ref++)
                {
                        set<double> simplex_integrals;
                        
                        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());                            

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

                                if(ret_val.second == true)
                                {
                                        triangulation.push_back(new_simplex);

                                        if(should_integrate)
                                        {
                                                toprow* as_toprow = (toprow*)this;

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

                                                vertex* base_vertex = (*new_simplex.begin());

                                                int dimension = new_simplex.size()-1;

                                                mat jacobian(dimension,dimension);

                                                double a_0 = base_vertex->get_coordinates()*cur_condition+as_toprow->condition[0];
                                                vec as;

                                                int row_count = 0;

                                                for(set<vertex*>::iterator vert_ref = (++new_simplex.begin()); vert_ref!=new_simplex.end();vert_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;                                                   
                                                                                                                
                                                        as.ins(as.size(),new_a);                                                        
                                                        
                                                        row_count++;
                                                }

                                                double int_value = 0;
                                                for(int a_count = 0;a_count<as.size();a_count++)
                                                {
                                                        vec reduced_as = (as.get(a_count)-as);

                                                        reduced_as.del(a_count);

                                                        int fac_order = ((toprow*)this)->condition_order-reduced_as.size()-1;

                                                        double fac_value = ((double)tgamma(fac_order)*det(jacobian))/as[a_count]/pow(as[a_count]-a_0,fac_order); 

                                                        for(int b_count = 0;b_count<reduced_as.size();b_count++)
                                                        {
                                                                fac_value /= reduced_as[b_count];                                                               
                                                        }

                                                        int_value += fac_value;                                                 
                                                }

                                                simplex_integrals.insert(int_value);
                                        }
                                } 
                        }

                        if(should_integrate)
                        {
                                ((toprow*)this)->probability = 0.0;

                                for(set<double>::iterator integ_ref = simplex_integrals.begin();integ_ref!=simplex_integrals.end();integ_ref++)
                                {
                                        ((toprow*)this)->probability += (*integ_ref);
                                }
                        }
                }               
        }