#include "robustlib.h" double polyhedron::triangulate(bool should_integrate) { for(set::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::iterator child_ref = children.begin();child_ref!=children.end();child_ref++) { for(set::iterator s_ref = (*child_ref)->triangulation.begin();s_ref!=(*child_ref)->triangulation.end();s_ref++) { simplex* new_simplex = new simplex((*s_ref)->vertices); pair::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::iterator child_ref = children.begin();child_ref!=children.end();child_ref++) { for(list>::iterator t_ref = (*child_ref)->triangulation.begin();t_ref!=(*child_ref)->triangulation.end();t_ref++) { set new_simplex; new_simplex.insert((*t_ref).begin(),(*t_ref).end()); for(set::iterator suitable_vert_ref = vertices.begin();*suitable_vert_ref<*(*t_ref).begin();suitable_vert_ref++) { set 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 as; vertex* base_vertex = simplex->vertices.begin(); for(set::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(amin_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::iterator,bool> returned = as.insert(pair(a,1)); if(returned.second == false) { (*returned.first).second++; } } /* cout << "a_0: " << a_0 << " "; int as_count = 1; for(map::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::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 factors; int number_of_factors = 0; for(map::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::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;sinsert_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; } }