#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 cur_as_ref = (*as_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;
		}
	}