Index: /applications/robust/robustlib.cpp
===================================================================
--- /applications/robust/robustlib.cpp (revision 1331)
+++ /applications/robust/robustlib.cpp (revision 1335)
@@ -2,5 +2,5 @@
 
 
-void polyhedron::triangulate(bool should_integrate)
+double polyhedron::triangulate(bool should_integrate)
 	{
 		for(set<simplex*>::iterator t_ref = triangulation.begin();t_ref!=triangulation.end();t_ref++)
@@ -48,5 +48,14 @@
 				}
 			}	
-		}		
+		}	
+
+		if(should_integrate)
+		{
+			return ((toprow *)this)->probability;
+		}
+		else
+		{
+			return 0.0;
+		}
 		
 		/*
@@ -165,15 +174,15 @@
 						
 
-						if(a_0<current_emlig->likelihood_value)
+						if(a_0<current_emlig->min_ll)
 						{
 							current_emlig->minimal_vertex = base_vertex;
-							current_emlig->likelihood_value = a_0;
+							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->likelihood_value)
+						if(a_m<current_emlig->min_ll)
 						{
 							current_emlig->minimal_vertex = (*vert_ref);
-							current_emlig->likelihood_value = a_m;						
+							current_emlig->min_ll = a_m;						
 						}
 
@@ -212,5 +221,5 @@
 			}*/
 
-			double int_value = 0;
+			double int_value = 0;			
 
 			// cout << jacobian << endl;
Index: /applications/robust/robustlib.h
===================================================================
--- /applications/robust/robustlib.h (revision 1334)
+++ /applications/robust/robustlib.h (revision 1335)
@@ -82,4 +82,6 @@
 
 	double negative_gamma_sum;
+
+	double min_beta;
 	
 
@@ -99,9 +101,10 @@
 	{
 		positive_gamma_parameters.clear();
-		negative_gamma_parameters.clear();
-		
+		negative_gamma_parameters.clear();		
 		
 		positive_gamma_sum = 0;
 		negative_gamma_sum = 0;
+
+		min_beta = numeric_limits<double>::max();
 	}
 
@@ -131,4 +134,9 @@
 
 			negative_gamma_parameters[order].insert(pair<double,double>(-weight,beta));
+		}
+
+		if(beta < min_beta)
+		{
+			min_beta = beta;
 		}
 	}
@@ -321,5 +329,5 @@
 
 	/// A method for triangulation of given polyhedron.
-	void triangulate(bool should_integrate);	
+	double triangulate(bool should_integrate);	
 };
 
@@ -854,4 +862,6 @@
 									//new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second));
 								}
+
+								cur_par_toprow->my_emlig->normalization_factor += cur_par_toprow->probability;
 
 								//current_parent->triangulation.clear();
@@ -959,4 +969,6 @@
 									//new_triangulation.insert(pair<double,set<vertex*>>(cur_prob,(*t_ref).second));
 								}
+
+								cur_par_toprow->my_emlig->normalization_factor += cur_par_toprow->probability;
 
 								//current_parent->triangulation.clear();
@@ -999,5 +1011,5 @@
 	vertex* minimal_vertex;
 
-	double likelihood_value;
+	double min_ll;
 
 	vector<multiset<my_ivec>> correction_factors;
@@ -1013,5 +1025,5 @@
 		create_statistic(number_of_parameters);
 
-		likelihood_value = numeric_limits<double>::max();
+		min_ll = numeric_limits<double>::max();
 	}
 
@@ -1022,5 +1034,5 @@
 		this->statistic = statistic;	
 
-		likelihood_value = numeric_limits<double>::max();
+		min_ll = numeric_limits<double>::max();
 	}
 
@@ -1146,6 +1158,6 @@
 	{
 		//step_me(0);
-		
-		likelihood_value = numeric_limits<double>::max();
+		normalization_factor = 0;
+		min_ll = numeric_limits<double>::max();
 
 		bool should_remove = (toremove.size() != 0);
@@ -1480,5 +1492,5 @@
 						current_positive->totally_neutral = (current_positive->totally_neutral && current_negative->totally_neutral);
 
-						current_positive->triangulate(k==for_splitting.size()-1);
+						current_positive->my_emlig->normalization_factor += current_positive->triangulate(k==for_splitting.size()-1);
 						
 						statistic.delete_polyhedron(k,current_negative);
@@ -1693,6 +1705,6 @@
 					new_totally_neutral_child->triangulate(false);
 
-					positive_poly->triangulate(k==for_splitting.size()-1);
-					negative_poly->triangulate(k==for_splitting.size()-1);
+					positive_poly->my_emlig->normalization_factor += positive_poly->triangulate(k==for_splitting.size()-1);
+					negative_poly->my_emlig->normalization_factor += negative_poly->triangulate(k==for_splitting.size()-1);
 					
 					statistic.append_polyhedron(k-1, new_totally_neutral_child);					
@@ -1790,10 +1802,10 @@
 		/// \TODO tady je to spatne, tady nesmi byt conditions.size(), viz RARX.bayes()
 		if(conditions.size()-2-number_of_parameters>=0)
-		{
-			
-			
+		{			
 			mat sample_mat;
 			map<double,toprow*> ordered_toprows;			
 			double sum_a = 0;
+			
+			//cout << "Likelihoods of toprows:" << endl;
 
 			for(polyhedron* top_ref = statistic.rows[number_of_parameters];top_ref!=statistic.end_poly;top_ref=top_ref->next_poly)
@@ -1802,6 +1814,18 @@
 
 				sum_a+=current_top->probability;
+				/*
+				cout << current_top->probability << "   ";
+
+				for(set<vertex*>::iterator vert_ref = (*top_ref).vertices.begin();vert_ref!=(*top_ref).vertices.end();vert_ref++)
+				{
+					cout << round(100*(*vert_ref)->get_coordinates())/100 << " ; ";
+				}
+				*/
+
+				// cout << endl;
 				ordered_toprows.insert(pair<double,toprow*>(sum_a,current_top));
-			}
+			}			
+			
+			// cout << "Sum N: " << normalization_factor << endl;
 
 			while(sample_mat.cols()<n)
@@ -1847,6 +1871,10 @@
 
 				cout << "Simplex/Count: " << simplex_count << "/" << sampled_toprow->triangulation.size() << endl;
+				cout << "Simplex factor: " << (*s_ref)->probability << endl;
+				cout << "Toprow factor:  " << sampled_toprow->probability << endl;
+				cout << "Emlig factor:   " << normalization_factor << endl;
 				// cout << &(*tri_ref) << endl;
 
+				int number_of_runs = 0;
 				bool have_sigma = false;
 				double sigma = 0;
@@ -1862,10 +1890,13 @@
 							sum_g += (*g_ref).first/(*s_ref)->positive_gamma_sum;
 
+							
 							if(sum_g>rnumber)
 							{
 								itpp::Gamma_RNG* gamma = new itpp::Gamma_RNG(conditions.size()-number_of_parameters,1/(*g_ref).second);
 								sigma = 1/(*gamma)();
+								
+								cout << "Sigma mean:   " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl;								
 								break;
-							}						
+							}							
 						}
 
@@ -1883,5 +1914,5 @@
 						for(multimap<double,double>::iterator pg_ref = (*v_ref).begin();pg_ref!=(*v_ref).end();pg_ref++)
 						{
-							pg_sum += exp(-(*pg_ref).second/sigma)*pow((*pg_ref).second/sigma,(int)conditions.size()-number_of_parameters-1)*(*pg_ref).second/fact(conditions.size()-number_of_parameters-1)*(*pg_ref).first;
+							pg_sum += exp(((*s_ref)->min_beta-(*pg_ref).second)/sigma)*pow((*pg_ref).second/sigma,(int)conditions.size()-number_of_parameters-1)*(*pg_ref).second/fact(conditions.size()-number_of_parameters-1)*(*pg_ref).first;
 						}					
 					}
@@ -1892,5 +1923,5 @@
 						for(multimap<double,double>::iterator ng_ref = (*v_ref).begin();ng_ref!=(*v_ref).end();ng_ref++)
 						{
-							ng_sum += exp(-(*ng_ref).second/sigma)*pow((*ng_ref).second/sigma,(int)conditions.size()-number_of_parameters-1)*(*ng_ref).second/fact(conditions.size()-number_of_parameters-1)*(*ng_ref).first;
+							ng_sum += exp(((*s_ref)->min_beta-(*ng_ref).second)/sigma)*pow((*ng_ref).second/sigma,(int)conditions.size()-number_of_parameters-1)*(*ng_ref).second/fact(conditions.size()-number_of_parameters-1)*(*ng_ref).first;
 						}					
 					}
@@ -1900,6 +1931,11 @@
 						have_sigma = true;
 					}
+
+					number_of_runs++;
 				}
 				while(!have_sigma);
+
+				cout << "Sigma: " << sigma << endl;
+				cout << "Nr. of runs: " << number_of_runs << endl;
 
 				int dimension = (*s_ref)->vertices.size()-1;
@@ -1909,4 +1945,6 @@
 
 				vertex* base_vert = *(*s_ref)->vertices.begin();
+
+				cout << "Base vertex coords(should be close to est. param.): " << base_vert->get_coordinates() << endl;
 				
 				int row_count = 0;
@@ -1914,10 +1952,14 @@
 				for(set<vertex*>::iterator vert_ref = ++(*s_ref)->vertices.begin();vert_ref!=(*s_ref)->vertices.end();vert_ref++)
 				{
-					vec relative_coords = (*vert_ref)->get_coordinates()-base_vert->get_coordinates();				
+					vec current_coords = (*vert_ref)->get_coordinates();
+
+					cout << "Coords of vertex[" << row_count << "]: " << current_coords << endl; 
+					
+					vec relative_coords = current_coords-base_vert->get_coordinates();				
 
 					jacobian.set_row(row_count,relative_coords);
 
 					row_count++;
-				}
+				}				
 				
 				cout << "Jacobian: " << jacobian << endl;
@@ -1926,4 +1968,6 @@
 								
 				gradient = jacobian*gradient;	
+
+				cout << "Gradient after trafo:" << gradient << endl;
 
 				// vec normal_gradient = gradient/sqrt(gradient*gradient);
@@ -1943,5 +1987,5 @@
 					x_axis = rotation_matrix*x_axis;
 
-					double t = gradient[i]/gradient*x_axis;
+					double t = abs(gradient[i]/gradient*x_axis);
 
 					double sin_theta = sign(gradient[i])*t/sqrt(1+pow(t,2));
@@ -1961,43 +2005,21 @@
 
 				// cout << rotation_matrix << endl;
-				cout << "Gradient after trafo:" << gradient << endl;
-
-				vec minima = numeric_limits<double>::max()*ones(number_of_parameters);
-				vec maxima = -numeric_limits<double>::max()*ones(number_of_parameters);
-
-				vertex* min_grad;
-				vertex* max_grad;
-								
-				for(set<vertex*>::iterator vert_ref = (*s_ref)->vertices.begin();vert_ref!=(*s_ref)->vertices.end();vert_ref++)
-				{
-					vec new_coords = rotation_matrix*jacobian*(*vert_ref)->get_coordinates();
-
-					// cout << new_coords << endl;
-
-					for(int j = 0;j<new_coords.size();j++)
-					{
-						if(new_coords[j]>maxima[j])
-						{
-							maxima[j] = new_coords[j];
-
-							if(j==0)
-							{
-								max_grad = (*vert_ref);
-							}
-						}
-
-						if(new_coords[j]<minima[j])
-						{
-							minima[j] = new_coords[j];
-
-							if(j==0)
-							{
-								min_grad = (*vert_ref);
-							}
-						}
-					}					
-				}
-
-				vec sample_coordinates;				
+				
+				mat extended_rotation = rotation_matrix;
+				extended_rotation.ins_col(0,zeros(extended_rotation.rows()));
+
+				cout << "Extended rotation: " << extended_rotation << endl;
+				
+				vec minima = itpp::min(extended_rotation,2);
+				vec maxima = itpp::max(extended_rotation,2);
+
+				cout << "Minima: " << minima << endl;
+				cout << "Maxima: " << maxima << endl;
+
+				vec sample_coordinates;		
+				bool is_inside = true;
+				
+				vec new_sample;
+				sample_coordinates = new_sample;
 
 				for(int j = 0;j<number_of_parameters;j++)
@@ -2008,17 +2030,54 @@
 
 					if(j==0)
-					{
-						double pos_value = sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0];
+					{						
+						vec new_gradient = rotation_matrix*gradient;
 						
-						vec new_gradient = rotation_matrix*gradient;
-
 						cout << "New gradient(should have only first component nonzero):" << new_gradient << endl;
 
-						cout << "Max: " << maxima[0] << "  Min: " << minima[0] << "  Grad:" << new_gradient[0] << endl;
+						// cout << "Max: " << maxima[0] << "  Min: " << minima[0] << "  Grad:" << new_gradient[0] << endl;
 						
-						double log_bracket = rnumber*exp(new_gradient[0]*maxima[0]/sigma)+(1-rnumber)*exp(new_gradient[0]*minima[0]/sigma);
+						double log_bracket = 1-rnumber*(1-exp(new_gradient[0]/sigma*(minima[0]-maxima[0])));
 						
-						coordinate = log(log_bracket);
-					}
+						coordinate = minima[0]-sigma/new_gradient[0]*log(log_bracket);
+					}
+					else
+					{
+						coordinate = minima[j]+rnumber*(maxima[j]-minima[j]);
+					}
+
+					sample_coordinates.ins(j,coordinate);
+				}
+
+				cout << "Sampled coordinates(gradient direction): " << sample_coordinates << endl;
+
+				sample_coordinates = rotation_matrix.T()*sample_coordinates;
+
+				cout << "Sampled coordinates(backrotated direction):" << sample_coordinates << endl;
+
+				
+				for(int j = 0;j<sample_coordinates.size();j++)
+				{
+					if(sample_coordinates[j]<0)
+					{
+						is_inside = false;
+					}
+				}
+
+				double above_criterion = ones(sample_coordinates.size())*sample_coordinates;
+
+				if(above_criterion>1)
+				{
+					is_inside = false;
+				}
+
+				if(is_inside)
+				{					
+					sample_coordinates = jacobian.T()*sample_coordinates+(*base_vert).get_coordinates();
+					
+					sample_coordinates.ins(0,sigma);
+					
+					cout << "Sampled coordinates(parameter space):" << sample_coordinates << endl;
+
+					sample_mat.ins_col(0,sample_coordinates);
 				}
 
@@ -2026,4 +2085,5 @@
 				// cout << sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0] << endl;
 
+				
 			}
 
