#include "robustlib.h"


void polyhedron::triangulate(bool should_integrate)
	{
		this->triangulation.clear();

		if(should_integrate)
		{
			((toprow *)this)->probability = 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++)
				{
					new_simplex.insert(*suitable_vert_ref);

					triangulation.push_back(new_simplex);

					if(should_integrate)
					{
						((toprow *)this)->probability += ((toprow *)this)->integrate_simplex(new_simplex, 'S');
					}
				}				
			}	
		}		
	}


	double toprow::integrate_simplex(set<vertex*> simplex, char c)
	{
		// cout << ((toprow*)this)->condition << endl;		

		int condition_order = ((toprow*)this)->condition_order-2; // -2 by bylo, pokud chceme uniformni apriorno

		int sigma_order = condition_order-my_emlig->number_of_parameters;

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

			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.get(1,as_toprow->condition.size()-1);

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

			int dimension = simplex.size()-1;

			mat jacobian(dimension,dimension);			

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

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

				base_vertex = (*base_ref);

			
				cout << endl << "Base coords:" << base_vertex->get_coordinates() << endl;

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

				int row_count = 0;
				for(set<vertex*>::iterator vert_ref = simplex.begin(); vert_ref!=simplex.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.end())
							{
								throw new exception("Equal local conditions are paired. If this ever occurs, the software has to be modified to include multiplied a_0!!");
							}
							
							order_correct = false;						

							break;
						}						
						
						

						if(a_0<current_emlig->likelihood_value)
						{
							current_emlig->minimal_vertex = base_vertex;
							current_emlig->likelihood_value = a_0;
						}

						double a_m = (*vert_ref)->get_coordinates()*cur_condition-as_toprow->condition[0];
						if(a_m<current_emlig->likelihood_value)
						{
							current_emlig->minimal_vertex = (*vert_ref);
							current_emlig->likelihood_value = 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 bracket = fact(sigma_order)/fact(as1_order-1)/pow(a_0-(*as_ref).first,sigma_order+1);
				for(int k = 0;k < as1_order-1;k++)
				{
 					double bracket_factor = -pow((double)-1,k+1)*fact(sigma_order-k)/fact(as1_order-1-k)/pow(a_0-(*as_ref).first,sigma_order-k);
					
					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];
						}

						bracket+=bracket_factor*bracket_combination;									
					}									
				}

				

				int_value += multiplier*bracket;
																		
			}

			
			correction_term *= fact(sigma_order)/abs(pow(a_0,sigma_order+1));

			// cout << c << int_value << endl;

			int_value += correction_term;
		

			//cout << "Probability:" << int_value << endl;
			
			return int_value;

			
			
		}
		else
		{
			return 0.0;
		}
	}