Index: /applications/robust/main.cpp
===================================================================
--- /applications/robust/main.cpp (revision 1343)
+++ /applications/robust/main.cpp (revision 1346)
@@ -18,5 +18,5 @@
 
 const int emlig_size = 2;
-const int utility_constant = 2;
+const int utility_constant = 3;
 
 
@@ -139,5 +139,5 @@
 	vector<vector<string>> strings;
 
-	char* file_strings[3] = {"c:\\dataADClosePercDiff","c:\\ar_student_single","c:\\ar_cauchy_single"};
+	char* file_strings[3] = {"c:\\Users\\Hontik\\Desktop\\dataGCClosePercDiff","c:\\ar_student_single","c:\\ar_cauchy_single"};
 
 	for(int i = 0;i<3;i++)
@@ -171,5 +171,5 @@
 		vector<vec> conditions;
 		//emlig* emliga = new emlig(2);
-		RARX* my_rarx = new RARX(3,30,true);
+		RARX* my_rarx = new RARX(2,30,false);
 		
 		
@@ -211,23 +211,23 @@
 					
 				
-				if(k>10)
+				if(k>8)
 				{
 					//my_rarx->posterior->step_me(0);
 
-					mat samples = my_rarx->posterior->sample_mat(10);
-
-					pair<vec,mat> imp_samples = my_rarx->posterior->importance_sample(20);
-
-					cout << imp_samples.first << endl;
+					//mat samples = my_rarx->posterior->sample_mat(10);
+
+					pair<vec,mat> imp_samples = my_rarx->posterior->importance_sample(1000);
+
+					//cout << imp_samples.first << endl;
 					
 					vec sample_prediction;
-					for(int t = 0;t<samples.cols();t++)
+					for(int t = 0;t<imp_samples.first.size();t++)
 					{
 						vec lap_sample = conditions[k-3].left(2);
-						lap_sample.ins(lap_sample.size(),1.0);
+						//lap_sample.ins(lap_sample.size(),1.0);
 						
 						lap_sample.ins(0,LapRNG());
 
-						sample_prediction.ins(0,lap_sample*samples.get_col(t));
+						sample_prediction.ins(0,lap_sample*imp_samples.second.get_col(t));
 					}
 
@@ -237,9 +237,15 @@
 					// cout << sample_prediction << endl;
 					vec poly_coefs;
+					double prediction;
 					bool stop_iteration = false;
 					int en = 1;
 					do
 					{
-						double poly_coef = ones(sample_pow.size())*sample_pow/sample_pow.size();
+						double poly_coef = imp_samples.first*sample_pow/(imp_samples.first*ones(imp_samples.first.size()));
+
+						if(en==1)
+						{
+							prediction = poly_coef;
+						}
 
 						poly_coef = poly_coef*en*fact(utility_constant-2+en)/fact(utility_constant-2);
@@ -278,5 +284,5 @@
 					*/
 
-					// cout << "Coefficients: " << poly_coefs << endl;
+					cout << "Coefficients: " << poly_coefs << endl;
 										
 					/*
@@ -287,4 +293,6 @@
 					*/
 					
+
+
 					cvec actions = roots(poly_coefs);
 					
@@ -295,6 +303,4 @@
 						if(actions[t].imag() == 0)
 						{
-							
-
 							double second_derivative = 0;
 							for(int q = 1;q<poly_coefs.size();q++)
@@ -319,9 +325,10 @@
 					// cout << "MaxLik coords:" << my_rarx->posterior->minimal_vertex->get_coordinates() << endl;
 
+					/*
 					double prediction = 0;
 					for(int s = 1;s<samples.rows();s++)
 					{
 						
-						double avg_parameter = samples.get_row(s)*ones(samples.cols())/samples.cols();
+						double avg_parameter = imp_samples.get_row(s)*ones(samples.cols())/samples.cols();
 
 						prediction += avg_parameter*conditions[k-3][s-1];
@@ -355,7 +362,9 @@
 						myfile.close();
 						*/
-					}
-
-					cout << "Prediction: "<< prediction << endl;
+
+					
+					//}
+
+					// cout << "Prediction: "<< prediction << endl;
 					
 					enorm<ldmat>* pred_mat = my_arx->epredictor(conditions[k-3].left(2));
@@ -401,5 +410,5 @@
 					}
 					myfile.close();
-					
+					//*/
 
 				}					
Index: /applications/robust/robustlib.h
===================================================================
--- /applications/robust/robustlib.h (revision 1343)
+++ /applications/robust/robustlib.h (revision 1346)
@@ -847,5 +847,5 @@
 							//current_parent->set_state(0,MERGE);	
 
-							if(level == number_of_parameters - 1)
+							if((level == number_of_parameters - 1) && (!shouldsplit))
 							{
 								toprow* cur_par_toprow = ((toprow*)current_parent);
@@ -863,5 +863,5 @@
 								}
 
-								cur_par_toprow->my_emlig->normalization_factor += cur_par_toprow->probability;
+								normalization_factor += cur_par_toprow->probability;
 
 								//current_parent->triangulation.clear();
@@ -873,5 +873,5 @@
 						{
 							for_merging[level+1].push_back(current_parent);
-							//current_parent->parentconditions.erase(toremove);
+							//current_parent->parentconditions.erase(toremove);							
 						}						
 
@@ -953,5 +953,5 @@
 							((toprow*)current_parent)->condition_order++;
 
-							if(level == number_of_parameters - 1)
+							if(level == number_of_parameters - 1 && current_parent->mergechild == NULL)
 							{
 								toprow* cur_par_toprow = ((toprow*)current_parent);
@@ -969,5 +969,5 @@
 								}
 
-								cur_par_toprow->my_emlig->normalization_factor += cur_par_toprow->probability;
+								normalization_factor += cur_par_toprow->probability;
 
 								//current_parent->triangulation.clear();
@@ -1154,6 +1154,5 @@
 		vec null_vector = "";
 
-		add_and_remove_condition(null_vector, toremove);
-	
+		add_and_remove_condition(null_vector, toremove);	
 	}
 
@@ -1338,4 +1337,6 @@
 					else
 					{
+						bool will_be_split = false;
+						
 						toprow* current_positive = (toprow*)(*merge_ref)->positiveparent;
 						toprow* current_negative = (toprow*)(*merge_ref)->negativeparent;
@@ -1395,4 +1396,5 @@
 							current_positive->set_state(0,SPLIT);
 							for_splitting[k].push_back(current_positive);
+							will_be_split = true;
 						}
 						
@@ -1428,4 +1430,6 @@
 							current_positive->negativeneutralvertices.insert(current_negative->negativeneutralvertices.begin(),current_negative->negativeneutralvertices.end());
 							current_positive->positiveneutralvertices.insert(current_negative->positiveneutralvertices.begin(),current_negative->positiveneutralvertices.end());
+
+							will_be_split = true;
 						}
 						else
@@ -1505,5 +1509,5 @@
 						current_positive->totally_neutral = (current_positive->totally_neutral && current_negative->totally_neutral);
 
-						current_positive->my_emlig->normalization_factor += current_positive->triangulate(k==for_splitting.size()-1);
+						normalization_factor += current_positive->triangulate(k==for_splitting.size()-1 && !will_be_split);
 						
 						statistic.delete_polyhedron(k,current_negative);
@@ -1722,6 +1726,6 @@
 					new_totally_neutral_child->triangulate(false);
 
-					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);
+					normalization_factor += positive_poly->triangulate(k==for_splitting.size()-1);
+					normalization_factor += negative_poly->triangulate(k==for_splitting.size()-1);
 					
 					statistic.append_polyhedron(k-1, new_totally_neutral_child);					
@@ -1822,5 +1826,5 @@
 		double  prob_sum     = 0;	
 		toprow* sampled_toprow;				
-		for(polyhedron* top_ref = statistic.rows[statistic.size()-1];top_ref!=statistic.row_ends[statistic.size()-1];top_ref=top_ref->next_poly)
+		for(polyhedron* top_ref = statistic.rows[number_of_parameters];top_ref!=statistic.end_poly;top_ref=top_ref->next_poly)
 		{
 			// cout << "CDF:"<< (*top_ref).first << endl;
@@ -1834,5 +1838,12 @@
 				sampled_toprow = (toprow*)top_ref;
 				break;
-			}						
+			}
+			else
+			{
+				if(top_ref->next_poly==statistic.end_poly)
+				{
+					cout << "Error.";
+				}
+			}
 		}				
 
@@ -1857,54 +1868,61 @@
 	pair<double,double> choose_sigma(simplex* sampled_simplex)
 	{
-		double rnumber = randu();
-		double sigma;
-					
-		double sum_g = 0;
-		for(int i = 0;i<sampled_simplex->positive_gamma_parameters.size();i++)
-		{
-			for(multimap<double,double>::iterator g_ref = sampled_simplex->positive_gamma_parameters[i].begin();g_ref != sampled_simplex->positive_gamma_parameters[i].end();g_ref++)
-			{
-				sum_g += (*g_ref).first/sampled_simplex->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)();
+		double sigma = 0;
+		double pg_sum;
+		double ng_sum;
+		do
+		{			
+			double rnumber = randu();
+			
+						
+			double sum_g = 0;
+			for(int i = 0;i<sampled_simplex->positive_gamma_parameters.size();i++)
+			{
+				for(multimap<double,double>::iterator g_ref = sampled_simplex->positive_gamma_parameters[i].begin();g_ref != sampled_simplex->positive_gamma_parameters[i].end();g_ref++)
+				{
+					sum_g += (*g_ref).first/sampled_simplex->positive_gamma_sum;
+
 								
-					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;								
+					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;
-				}							
-			}
-
-			if(sigma!=0)
-			{
-				break;
-			}
-		}
-
-		rnumber = randu();
-
-		double pg_sum = 0;
-		for(vector<multimap<double,double>>::iterator v_ref = sampled_simplex->positive_gamma_parameters.begin();v_ref!=sampled_simplex->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((sampled_simplex->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 = sampled_simplex->negative_gamma_parameters.begin();v_ref!=sampled_simplex->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((sampled_simplex->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;
-			}					
-		}
+				}
+			}
+
+			rnumber = randu();
+
+			pg_sum = 0;
+			for(vector<multimap<double,double>>::iterator v_ref = sampled_simplex->positive_gamma_parameters.begin();v_ref!=sampled_simplex->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((sampled_simplex->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;
+				}					
+			}
+
+			ng_sum = 0;
+			for(vector<multimap<double,double>>::iterator v_ref = sampled_simplex->negative_gamma_parameters.begin();v_ref!=sampled_simplex->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((sampled_simplex->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;
+				}					
+			}
+		}
+		while(pg_sum-ng_sum<0);
 
 		return pair<double,double>((pg_sum-ng_sum)/pg_sum,sigma);
@@ -2277,5 +2295,6 @@
 			extended_coords.ins(0,-1.0);
 
-			double sample_prob = 1/condition_and_simplex.second->probability/pow(probability_and_sigma.second,(int)conditions.size()-number_of_parameters)*exp(1/probability_and_sigma.second*(extended_coords*condition_and_simplex.first));
+			double exponent = extended_coords*condition_and_simplex.first;
+			double sample_prob = 1/condition_and_simplex.second->probability/pow(probability_and_sigma.second,(int)conditions.size()-number_of_parameters)*exp((-1)/probability_and_sigma.second*exponent);
 			sample_prob *= probability_and_sigma.first;
 
