#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];
						map<double,int> 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;							
							
							pair<map<double,int>::iterator,bool> returned = as.insert(pair<double,int>(new_a,1));
							if(returned.second == false)
							{
								(*returned.first).second++;
							}
							//as.ins(as.size(),new_a);							
							
							row_count++;
						}

						double int_value = 0;
						for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++)
						{
							int fac_order = ((toprow*)this)->condition_order-as.size()-1;

							double fac_value = 1;
							if(a_0!=(*as_ref).first)
							{
								//fac_value = ((double)tgamma(fac_order)*det(jacobian))/(*as_ref).first/pow((*as_ref).first-a_0,fac_order); 
								
								fac_value /= (*as_ref).first;

								for(map<double,int>::iterator as2_ref = as.begin();as2_ref!=as.end();as2_ref++)
								{
									if(as2_ref!=as_ref)
									{
										double current_factor = (*as_ref).first-(*as2_ref).first;

										fac_value /= current_factor;

																				
									}
								}
							}
							else
							{
								throw new exception("Should this happen? a_i=a_0 in the formula for integrating over a simplex! I think it can happen with 0 probability(when phi*2*v_0=phi*v_i)");
							}

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