Index: applications/robust/robustlib.cpp
===================================================================
--- applications/robust/robustlib.cpp (revision 1270)
+++ applications/robust/robustlib.cpp (revision 1271)
@@ -3,8 +3,11 @@
 void polyhedron::triangulate(bool should_integrate)
 	{		
+		if(should_integrate)
+		{
+			((toprow *)this)->probability = 0.0;
+		}		
+		
 		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++)
 			{
@@ -18,136 +21,147 @@
 					triangulation.push_back(new_simplex);
 
-					int condition_order = ((toprow*)this)->condition_order-1;
+					if(should_integrate)
+					{
+						((toprow *)this)->probability += ((toprow *)this)->integrate_simplex(new_simplex, 'S');
+					}
+				} 
+			}	
+		}		
+	}
 
-					// cout << "C:" << condition_order << "  N:" << my_emlig->number_of_parameters << "  C+N:" << condition_order-my_emlig->number_of_parameters << endl;
-					// pause(0.1);
-					if(should_integrate && condition_order-my_emlig->number_of_parameters >= 0)
+
+	double toprow::integrate_simplex(set<vertex*> simplex, char c)
+	{
+		int condition_order = ((toprow*)this)->condition_order-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(condition_order-my_emlig->number_of_parameters >= 0)
+		{
+			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);
+
+			vertex* base_vertex = (*simplex.begin());
+
+			int dimension = 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 = (++simplex.begin()); vert_ref!=simplex.end();vert_ref++)
+			{
+				vec relative_coords = (*vert_ref)->get_coordinates()-base_vertex->get_coordinates();
+
+				// cout << relative_coords << endl;
+
+				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++;
+				}
+				else
+				{
+					cout << "a[" << row_count << "] = " << new_a << endl;
+				}
+				//as.ins(as.size(),new_a);							
+				
+				row_count++;
+			}
+
+			// cout << endl;
+
+			double int_value = 0;
+
+			// cout << jacobian << endl;
+
+			double det_jacobian = det(jacobian);
+			for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++)
+			{
+				double multiplier = det_jacobian;
+				if(a_0!=(*as_ref).first)
+				{								
+					int as1_order = (*as_ref).second;
+
+					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++)
 					{
-						emlig* current_emlig;
+						
+						if(as2_ref!=as_ref)
+						{										
+							for(int k = 0;k<(*as2_ref).second;k++)
+							{
+								factors.push_back((*as_ref).first-(*as2_ref).first);
+							}
 
-						if(this->my_emlig!=NULL)
-						{
-							current_emlig = this->my_emlig;
+							multiplier        /= pow((*as_ref).first-(*as2_ref).first,(*as2_ref).second);
+							number_of_factors += (*as2_ref).second;
 						}
 						else
 						{
-							throw exception("The statistic of the polyhedron you are trying to integrate over doesn't belong to any emlig!");
-						}						
+							factors.push_back((*as_ref).first);
 
-						toprow* as_toprow = (toprow*)this;
+							multiplier        /= (*as_ref).first;
+							number_of_factors += 1;
+						}									
+					}					
 
-						vec cur_condition = as_toprow->condition.get(1,as_toprow->condition.size()-1);
+					double bracket = fact(condition_order-number_of_factors)/fact(as1_order-1)/pow((*as_ref).first-a_0,condition_order-number_of_factors);
+					for(int k = 0;k < as1_order-1;k++)
+					{
+						double bracket_factor = pow((double)-1,k+1)*fact(condition_order-1-number_of_factors-k)/fact(as1_order-2-k)/pow((*as_ref).first-a_0,condition_order-1-number_of_factors-k);
+						
+						for(set<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[k].begin();combi_ref!=this->my_emlig->correction_factors[k].end();combi_ref++)
+						{
+							double bracket_combination = 1;
+							for(int j = 0;j<(*combi_ref).size();j++)
+							{
+								bracket_combination /= factors[(*combi_ref)[j]];
+							}
 
-						vertex* base_vertex = (*new_simplex.begin());
+							bracket+=bracket_factor*bracket_combination;									
+						}									
+					}								
 
-						int dimension = new_simplex.size()-1;
+					int_value += multiplier*bracket;	
 
-						mat jacobian(dimension,dimension);
+				}
+				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)");
+				}
 
-						double a_0 = base_vertex->get_coordinates()*cur_condition+as_toprow->condition[0];
-						map<double,int> as;
+				// cout << c << int_value << endl;
 
-						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;
-						double det_jacobian = det(jacobian);
-						for(map<double,int>::iterator as_ref = as.begin();as_ref!=as.end();as_ref++)
-						{
-							double multiplier = det_jacobian;
-							if(a_0!=(*as_ref).first)
-							{								
-								int as1_order = (*as_ref).second;
-
-								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 facti = fact(0);
-
-								double bracket = fact(condition_order-number_of_factors)/fact(as1_order-1)/pow((*as_ref).first-a_0,condition_order-number_of_factors);
-								for(int k = 0;k < as1_order-1;k++)
-								{
-									double bracket_factor = pow((double)-1,k+1)*fact(condition_order-1-number_of_factors-k)/fact(as1_order-2-k)/pow((*as_ref).first-a_0,condition_order-1-number_of_factors-k);
-									
-									for(set<my_ivec>::iterator combi_ref = this->my_emlig->correction_factors[k].begin();combi_ref!=this->my_emlig->correction_factors[k].end();combi_ref++)
-									{
-										double bracket_combination = 1;
-										for(int j = 0;j<(*combi_ref).size();j++)
-										{
-											bracket_combination /= factors[(*combi_ref)[j]];
-										}
-
-										bracket+=bracket_factor*bracket_combination;									
-									}									
-								}								
-
-								int_value += multiplier*bracket;	
-
-							}
-							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)");
-							}						
-													
-						}
-
-						cout << int_value << endl;
-						pause(0.1);
-
-
-						simplex_integrals.insert(int_value);
-					}
-				} 
+				return 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);
-				}
-			}
-		}		
+			
+		}
+		else
+		{
+			return 0.0;
+		}
 	}
 
