Index: applications/robust/robustlib.h
===================================================================
--- applications/robust/robustlib.h (revision 1395)
+++ applications/robust/robustlib.h (revision 1396)
@@ -2252,5 +2252,6 @@
 				if(sum_g>rnumber)
 				{					
-					GamRNG.setup(condition_order-number_of_parameters-1+i,1/(*g_ref).second);
+					// tady je mozna chyba ve vaskove kodu
+					GamRNG.setup(condition_order-number_of_parameters-1+i,(*g_ref).second);
 																
 					sigma = 1/GamRNG();
@@ -2273,5 +2274,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,condition_order-number_of_parameters-1+i)/sigma/fact(condition_order-number_of_parameters-1+i)*(*pg_ref).first;  // pg_sum += exp((sampled_simplex->min_beta-(*pg_ref).second)/sigma)*pow((*pg_ref).second/sigma,(int)conditions.size())*(*pg_ref).second/fact(conditions.size())*(*pg_ref).first;
+				pg_sum += exp(-(*pg_ref).second/sigma)*pow((*pg_ref).second/sigma,condition_order-number_of_parameters-1+i)/sigma/fact(condition_order-number_of_parameters-2+i)*(*pg_ref).first;  // pg_sum += exp((sampled_simplex->min_beta-(*pg_ref).second)/sigma)*pow((*pg_ref).second/sigma,(int)conditions.size())*(*pg_ref).second/fact(conditions.size())*(*pg_ref).first;
 			}
 
@@ -2279,325 +2280,13 @@
 		}	
 
-		return pair<double,double>(sum_g/pg_sum,sigma);
-	}
-
-	mat sample_mat(int n)
-	{		
-
-		/// \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)
-			{
-				toprow* current_top = (toprow*)top_ref;
-
-				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)
-			{
-				//// cout << "*************************************" << endl;
-
-				
-				
-				double rnumber = randu()*sum_a;
-
-				// cout << "RND:" << rnumber << endl;
-
-				// This could be more efficient (log n), but map::upper_bound() doesn't let me dereference returned iterator
-				int toprow_count = 0;
-				toprow* sampled_toprow;				
-				for(map<double,toprow*>::iterator top_ref = ordered_toprows.begin();top_ref!=ordered_toprows.end();top_ref++)
-				{
-					// cout << "CDF:"<< (*top_ref).first << endl;
-					toprow_count++;
-
-					if((*top_ref).first >= rnumber)
-					{
-						sampled_toprow = (*top_ref).second;
-						break;
-					}						
-				}				
-
-				//// cout << "Toprow/Count: " << toprow_count << "/" << ordered_toprows.size() << endl;
-				// cout << &sampled_toprow << ";";
-
-				rnumber = randu();				
-
-				set<simplex*>::iterator s_ref;
-				double sum_b = 0;
-				int simplex_count = 0;
-				for(s_ref = sampled_toprow->triangulation.begin();s_ref!=sampled_toprow->triangulation.end();s_ref++)
-				{
-					simplex_count++;
-					
-					sum_b += (*s_ref)->probability;
-
-					if(sum_b/sampled_toprow->probability >= rnumber)
-						break;
-				}
-
-				//// 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;
-				do
-				{
-					rnumber = randu();
-					
-					double sum_g = 0;
-					for(int i = 0;i<(*s_ref)->positive_gamma_parameters.size();i++)
-					{
-						for(multimap<double,double>::iterator g_ref = (*s_ref)->positive_gamma_parameters[i].begin();g_ref != (*s_ref)->positive_gamma_parameters[i].end();g_ref++)
-						{
-							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)();
-								
-								GamRNG.setup(conditions.size()-number_of_parameters,(*g_ref).second);
-																
-								sigma = 1/GamRNG();
-
-								// cout << "Sigma mean:   " << (*g_ref).second/(conditions.size()-number_of_parameters-1) << endl;								
-								break;
-							}							
-						}
-
-						if(sigma!=0)
-						{
-							break;
-						}
-					}
-
-					rnumber = randu();
-
-					double pg_sum = 0;
-					for(vector<multimap<double,double>>::iterator v_ref = (*s_ref)->positive_gamma_parameters.begin();v_ref!=(*s_ref)->positive_gamma_parameters.end();v_ref++)
-					{
-						for(multimap<double,double>::iterator pg_ref = (*v_ref).begin();pg_ref!=(*v_ref).end();pg_ref++)
-						{
-							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;
-						}					
-					}
-
-					double ng_sum = 0;
-					for(vector<multimap<double,double>>::iterator v_ref = (*s_ref)->negative_gamma_parameters.begin();v_ref!=(*s_ref)->negative_gamma_parameters.end();v_ref++)
-					{
-						for(multimap<double,double>::iterator ng_ref = (*v_ref).begin();ng_ref!=(*v_ref).end();ng_ref++)
-						{
-							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;
-						}					
-					}
-					
-					if((pg_sum-ng_sum)/pg_sum>rnumber)
-					{
-						have_sigma = true;
-					}
-
-					number_of_runs++;
-				}
-				while(!have_sigma);
-
-				//// cout << "Sigma: " << sigma << endl;
-				//// cout << "Nr. of sigma runs: " << number_of_runs << endl;
-
-				int dimension = (*s_ref)->vertices.size()-1;
-
-				mat jacobian(dimension,dimension);
-				vec gradient = sampled_toprow->condition_sum.right(dimension);
-
-				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;
-
-				for(set<vertex*>::iterator vert_ref = ++(*s_ref)->vertices.begin();vert_ref!=(*s_ref)->vertices.end();vert_ref++)
-				{
-					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;
-
-				//// cout << "Gradient before trafo:" << gradient << endl;
-								
-				gradient = jacobian*gradient;	
-
-				//// cout << "Gradient after trafo:" << gradient << endl;
-
-				// vec normal_gradient = gradient/sqrt(gradient*gradient);
-				// cout << gradient << endl;
-				// cout << normal_gradient << endl;
-				// cout << sqrt(gradient*gradient) << endl;
-
-				mat rotation_matrix = eye(dimension);				
-
-								
-
-				for(int i = 1;i<dimension;i++)
-				{
-					vec x_axis = zeros(dimension);
-					x_axis.set(0,1);
-
-					x_axis = rotation_matrix*x_axis;
-
-					double t = abs(gradient[i]/gradient*x_axis);
-
-					double sin_theta = sign(gradient[i])*t/sqrt(1+pow(t,2));
-					double cos_theta = sign(gradient*x_axis)/sqrt(1+pow(t,2));
-
-					mat partial_rotation = eye(dimension);
-
-					partial_rotation.set(0,0,cos_theta);
-					partial_rotation.set(i,i,cos_theta);
-					
-					partial_rotation.set(0,i,sin_theta);
-					partial_rotation.set(i,0,-sin_theta);
-					
-					rotation_matrix = rotation_matrix*partial_rotation;				
-					
-				}
-
-				// cout << rotation_matrix << endl;
-				
-				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++)
-				{
-					rnumber = randu();
-					
-					double coordinate;
-
-					if(j==0)
-					{						
-						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;
-						
-						double log_bracket = 1-rnumber*(1-exp(new_gradient[0]/sigma*(minima[0]-maxima[0])));
-						
-						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);
-
-					// cout << sample_mat.cols() << ",";
-				}
-
-				// cout << sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*min_grad->get_coordinates()-sampled_toprow->condition_sum[0] << endl;
-				// cout << sampled_toprow->condition_sum.right(sampled_toprow->condition_sum.size()-1)*max_grad->get_coordinates()-sampled_toprow->condition_sum[0] << endl;
-
-				
-			}
-
-			cout << endl;
-			return sample_mat;
-		}
-		else
-		{
-			throw new exception("You are trying to sample from density that is not determined (parameters can't be integrated out)!");
+		return pair<double,double>(sampled_simplex->positive_gamma_sum/pg_sum,sigma);
+	}
+
+	pair<vec,mat> sample(int n, bool rejection)
+	{
+		vec probabilities;
+		mat samples;		
 		
-			return 0;
-		}
-
-		
-	}
-
-	pair<vec,mat> importance_sample(int n)
-	{
-		vec probabilities;
-		mat samples;
-		
-		for(int i = 0;i<n;i++)
+		while(samples.cols()<n)
 		{
 			pair<vec,simplex*> condition_and_simplex = choose_simplex();
@@ -2639,16 +2328,43 @@
 
 			double exponent = extended_coords*condition_and_simplex.first;
-			double sample_prob = 1*condition_and_simplex.second->volume/condition_and_simplex.second->probability/pow(2*probability_and_sigma.second,condition_order)*exp(-exponent/probability_and_sigma.second)*probability_and_sigma.first;			
+			double sample_prob = 1*condition_and_simplex.second->volume/condition_and_simplex.second->probability/pow(2*probability_and_sigma.second,condition_order)*exp(-exponent/probability_and_sigma.second);//*probability_and_sigma.first;			
 
 			sample_coords.ins(sample_coords.size(),probability_and_sigma.second);
-
+			
 			samples.ins_col(0,sample_coords);
 			probabilities.ins(0,sample_prob);
 
-			cout << "C:" << sample_coords << "   p:" << sample_prob << endl;
-			pause(1 );
-		}
-	
-		return pair<vec,mat>(probabilities,samples);
+			
+
+			//cout << "C:" << sample_coords << "   p:" << sample_prob << endl;
+			//pause(1);
+		}
+
+		if(rejection)
+		{
+			double max_prob = max(probabilities);
+
+			set<int> indices;
+			for(int i = 0;i<n;i++)
+			{
+				double randn = randu();
+
+				if(probabilities.get(i)<randn*max_prob)
+				{
+					indices.insert(i);
+				}
+			}
+
+			for(set<int>::reverse_iterator ind_ref = indices.rbegin();ind_ref!=indices.rend();ind_ref++)
+			{
+				samples.del_col(*ind_ref);				
+			}
+
+			return pair<vec,mat>(ones(samples.cols()),samples);
+		}
+		else
+		{	
+			return pair<vec,mat>(probabilities,samples);
+		}
 	}
 
@@ -3032,5 +2748,5 @@
 		//posterior->step_me(0);
 		
-		cout << "Current condition:" << yt << endl;
+		cout << "*************************************" << endl << "Current condition:" << yt << endl << "*************************************" << endl;
 
 		/// \TODO tohle je spatne, tady musi byt jiny vypocet poctu podminek, kdyby nejaka byla multiplicitni, tak tohle bude spatne
