Solved the cjloss deviation problem

.gitignore

@@ -44,3 +44,4 @@ doc/
 *.lat
 *.bnd
 *.plt
+*.pdf

GH

@@ -1,3 +1,13 @@
+Utána nézni, hogy eddig milyen eredmények vannak a gauss heurisztikáról?
+Alapos kisérleteket végezni, hogy mitől függ, hogy mekkora sugár kell, hogy működjön:
+- random lattice-kkel full enumerationt csinálni (olyan alacsony dimenzióban, ahol még működik a dolog)
+- fix determináns, véletlen lattice
+- kisérletek:
+	- fix sugár, változó lattice, mekkora a hiba (hiba előjelét is nézni)
+	- több különböző sugárral elvégezni és hisztogrammokat késziteni
+	- Hermite faktorok mennyire különböznek? Mennyire függ a hiba a hermite faktortól? 
+
+-------------------------------------------
 Érvek amik amellett szólnak, hogy a heurisztika számitással van a gond:
 - Mindkét enumeration algoritmus ugyanazt az eredményt adja.
 - Már a pruned esetben is ugyanazt adták és hasonló volt az eltérés mint a teljes esetben

TODO

@@ -1,28 +1,8 @@
-+ Bounding function inputot implementálni, hogy gyorsabban lehessen tesztelni (kapcsoló, hogy milyen legyen a kimenet NTL vagy gnuplot)
-+ A korlát gyökvonásánál biztositani, hogy felfelé kerekitse double-ra (RR-ben gyököt vonni és azt kerekiteni felfelé double-ra)
-	+ nem, illetve nem úgy működik az RR felfelé kerekités, kézzel kell implementálni double-el is rendesen működik dolog, nem kell kinlódni
-+ tesztelni a bounding function-t (BKZ 30, dimenzió 80)
-	+ tesztelni jól kerekitett R-el és kitalálni, hogy miért nem működik 
-		+ az egyik hiba az volt, hogy az algoritmus az R_i^2-eket várja én meg az R_i-ket adom meg neki.
-		+ most viszont sokkal tovább fut mint igérte. Lehetséges okok:
-			+ Végtelen ciklusba esett amiatt amit az NTL megpróbál kiküszübülni -nem ez a gond, az NTL-es változat is sokáig fut
-			+ Rosszul számolja a node-ok számát
-				+ paramétereket keresni, ahol az LLL lineáris pruninggal mérhető ideig fut
-				+ tesztadatokat generálni vele és más alacsony success rate-űvel (hogy végig kelljen mennie) leszámolni a node-okat
-				+ összevetni a gyakorlatban tapasztalt node- számot az előrejelzettel
-				+ R_k val kell normalizálni az egyes körökben, de én R_n-el normalizáltam
-					- először egy robosztus referenciaimplementációt csinálni (pvolandscale függvénybe kitenni a skálázást és a volume számitást)	
-					- Nagyon lassit vele, meg kellene próbálni skálázásal kiváltani
-			+ Valami el van cseszve az enumeration algoritmusban -nem volt
-	- Implementálni a bázis randomizálást
-+ debugolni a node predikálót:
-	+ olyan cjloss lattice-t konstruálni ahol nagyjából egybeesik a gauss heurisztika és a tényleges legrövidebb vektor és azzal letesztelni - továbbra is alulbecsül függetlenül a gaus heuristika pontosságától. - a cjloss lattice-kre valamiért nem megy
-	+ letesztelni random lattice-okkal - azokat is alulbecsli a leszámlálás vége meg elszáll
-		+ tévedés, rosszul implementáltam a random lattice-ket
-	+ utánanézni a gauss heurisztikának, hogy mit írnak róla és annak alapján finomítani a teszteket
-	+ Alaposabban megérteni az enumeration algoritmusát és debugolni
-	+ Ellenőrizni és debugolni a gauss heurisztika számítását (a cjloss lattice-re a grafikonom laposabb volt mint Phongé)
-+ a cjlossos tesztprogramból automatikus tesztet irni,  - a cjlossossal nem működik és nem is tudom, hogyan lehetne működésre birni (lásd GH)
+- Minimális pruning function keresését kikisérletezni:
+	- tesztelni a gyakorlatban a képletemet állitható kezdőértékkel (dupla és sima random lattice-el is)
+	- konstans pruninggal megkeresni a megfelelő arányokat az egyes dimenziókban (dupla és sima random lattice-el is)
+		- több különböző lattice-ra és összehasonlitani
+
 - a random lattice-s tesztprogramból automatikus teszteket irni
 	- node számlálás full (double)
 	- extreme pruning tesztet külön lefuttatni, hogy ha van még valami gond, akkor az itt kiderüljön

ratio.tmp

@@ -1 +1 @@
-3 22 144 700 3005 11883 40502 126184 362750 988043 2381045 5510916 11939329 23671558 45555319 79560409 126235491 192195797 298829270 450540580 637102013 839021320 1094276433 1269535717 1380349197 1469832912 1555795329 1549768249 1453589791 1291319860 1096733132 883994409 727279166 542925664 388830593 279242934 193941217 123415192 75655355 46379468 25107579 13496579 7105497 3606957 1786140 812953 367352 158658 67212 26555 10025 3598 1341 493 159 54 17 2 1 0 
+0 0 1 1 3 6 11 18 29 39 64 87 143 182 284 361 543 682 1001 1247 1727 1867 2477 2449 2998 2880 3275 3049 3480 3033 3390 2982 3040 2470 2316 1697 1566 1128 915 629 535 384 297 190 145 93 71 38 23 15 11 6 3 1 0 0 0 0 0 

src/boundary.cpp

@@ -105,6 +105,38 @@ double ball_vol(int k, double r) {
 		return 2 * pow(4*M_PI, k/2) * fact(k/2) / fact(k) * pow(r, k);
 }
 
+double integral_odd(int ltilde, int l, double tvec[], double vvec[]) {
+	double ret= 0;
+
+	if(ltilde==0) 
+		return 1;
+
+	vvec[ltilde-1]= integral_odd(ltilde-1, l, tvec, vvec);
+
+	ret-= pow(1-tvec[l-ltilde], (2*ltilde+1)/2.0)*pow(2,2*ltilde+1)*fact(ltilde+1)/fact(2*ltilde+2);
+	for(int i= 1; i < ltilde; i++)
+		ret+= pow(tvec[l-ltilde],ltilde-i)*vvec[i]/fact(ltilde-i)*((ltilde-i-1)%2==0?1:-1);
+
+	if(ltilde!=l)		
+		return ret;
+
+	return ret + pow(2,2*ltilde+1)*fact(ltilde+1)/fact(2*ltilde+2);
+}
+
+double integral_even(int ltilde, int l, double tvec[], double vvec[]) {
+	double ret= 0;
+
+	if(ltilde==0) 
+		return 1;
+
+	vvec[ltilde-1]= integral_even(ltilde-1, l, tvec, vvec);
+
+	for(int i= 0; i < ltilde; i++)
+		ret+= pow(tvec[l-ltilde],ltilde-i)*vvec[i]/fact(ltilde-i)*((ltilde-i-1)%2==0?1:-1);
+
+	return ret;
+}
+
 double polytope_volume(double vols[], double bounds[], int dim) {
 	double ret= 0;
 
@@ -277,31 +309,34 @@ double n_full_gsa(double R, double scale, int bsize, int n) {
 	return N; 
 }
 
-double pvol_and_scale(double Rvec[], int k) {
+double ci_prob(double Rvec[], int k) {
 	double ret;
-	int l= (k+1)/2-1;
+	int l= k/2;
 	double* bounds= new double[l];
-	double* polytope_vols= new double[l];
+	double* vvec= new double[l];
 
 
-/*	cout << "khalf: " << khalf << endl;
-	cout << "Rvec: [ ";
-	for(int i= 0; i< 2*khalf; i++)
+	/*cout << "# khalf: " << l << endl;
+	cout << "# Rvec: [ ";
+	for(int i= 0; i< 2*l; i++)
 		cout << Rvec[i] << " ";	
 	cout << "] "<< endl;*/
 
 	for(int i= 0; i< l; i++)
 		bounds[i]= Rvec[2*i]*Rvec[2*i]/(Rvec[k-1]*Rvec[k-1]);
 
-	/*cout << "bounds: [ ";
-	for(int i= 0; i< khalf; i++)
+	/*cout << "# bounds: [ ";
+	for(int i= 0; i< l; i++)
 		cout << bounds[i] << " ";	
-	cout << "] "<< endl << endl;*/
+	cout << "] "<< endl << "# " << endl;*/
 
-	ret= polytope_volume(polytope_vols, bounds, l);
+	if(k%2==0)
+		ret= integral_even(l, l, bounds, vvec)*fact(l);
+	else
+		ret= integral_odd(l, l, bounds, vvec)*fact(2*l+2)/(pow(2,2*l+1)*fact(l+1));
 
 	delete [] bounds;
-	delete [] polytope_vols;
+	delete [] vvec;
 
 	return ret;
 }
@@ -359,22 +394,28 @@ double t_extreme_reference(double Rvec[], double b_star_norm[], double t_node, d
 		for(int i= n-k; i < n; i++) 
 			denom*= b_star_norm[i]; 
 
-		int psd= (k+1)/2-1;
-		V_act= ball_vol(k, Rvec[k-1]) * pvol_and_scale(Rvec, k) * fact(psd); 
+		V_act= ball_vol(k, Rvec[k-1]) * ci_prob(Rvec, k); 
 
 		N+= V_act/denom/2;
 
-		cout << "# Simvol_" << psd << " = " <<  fact(psd) << endl;
-		cout << "# Polvol_" << psd << " = " <<  pvol_and_scale(Rvec, k) << endl;
+		//cout << "# Simvol_" << psd << " = " <<  fact(psd) << endl;
+		//cout << "# Polvol_" << psd << " = " <<  pvol_and_scale(Rvec, k) << endl;
+		//cout << "# Polprob_" << k/2 << " = " <<  polprob_even(Rvec, k) << endl;
 		cout << "# V_ball(" << k << ", " << Rvec[k-1] << ") = " << ball_vol(k, Rvec[k-1]) << endl;
-		cout << "# |b^*_" << n-k << "| = " << b_star_norm[n-k] << endl;
-		cout << "# V_" << k << " = " << V_act << endl;
+		//cout << "# |b^*_" << n-k << "| = " << b_star_norm[n-k] << endl;
+		//cout << "# V_" << k << " = " << V_act << endl;
 		cout << "# denom_" << k << " = " << denom << endl;
-		cout << "# N_" << k << " = " << V_act/denom/2 << endl; 
+		cout << "# ratio_" << k << " = " << pow(ball_vol(k, Rvec[k-1])/denom,1.0/k) << endl;
+		//cout << "# N_" << k << " = " << V_act/denom/2 << endl; 
 		myfile >> nodes; sum+= nodes; 
 		cout << "# Measured nodes: " << nodes << endl;
-		cout << "# Difference: " << nodes - V_act/denom/2 << endl;		
-		cout << "# Ratio: " << nodes/(V_act/denom/2) << endl;
+		//cout << "# Difference: " << nodes - V_act/denom/2 << endl;		
+		//cout << "# Ratio: " << nodes/(V_act/denom/2) << endl;
+		//cout << "# GH_" << k << " = " << pow(denom/ball_vol(k, 1),1.0/k) << endl;
+		cout << "# R_" << k << " = " << Rvec[k-1] << endl;
+		if( Rvec[k-1] < sqrt(k))
+			cout << "# WARNING! Gaussian Heuristics violation: " << k << "^(1/2) = " << sqrt(k) << endl;
+		//cout << "# R/GH in dim " << k << " = " << Rvec[k-1]/ pow(denom/ball_vol(k, 1),1.0/k) << endl;
 
 		//cout << k-1 << " " << nodes/(V_act/denom/2) << endl;
 		//cout << k-1 << " " << nodes - V_act/denom/2 << endl;		
@@ -386,7 +427,7 @@ double t_extreme_reference(double Rvec[], double b_star_norm[], double t_node, d
 	//N/= 2;
 	cout << "# Predicted nodes: " << N << endl;
 	cout << "# Measured nodes: " << sum << endl;
-	return (t_reduc+t_node*N)/(pvol_and_scale(Rvec, n/2)*fact(n/2)); 
+	return (t_reduc+t_node*N)/ci_prob(Rvec, n/2); 
 }
 
 double t_extreme(double Rvec[], double b_star_norm[], double t_node, double t_reduc, int n) {

src/cjloss.cpp

@@ -46,12 +46,12 @@ cjloss::cjloss(long dimension, double density, long seed) {
 
 	ZZ N;
 	conv(N,ceil(sqrt(dimension-1)));
-	for(int i= 1; i<dimension; i++){
-		basis[i-1][i]= 2;
-		basis[i-1][0]= 2*values[i-1]*N;
+	for(int i= 0; i<dimension-1; i++){
+		basis[i][i]= 2;
+		basis[i][dimension-1]= 2*values[i]*N;
 		basis[dimension-1][i]= 1;
 		}
-	basis[dimension-1][0]= 2*sum*N;
+	basis[dimension-1][dimension-1]= 2*sum*N;
 	}
 
 /* generates a random knapsack problem with equally many zeroes and ones in the solution (condition applied in the extreme pruning article) */ 

test/CMakeLists.txt

@@ -1,6 +1,7 @@
 SET(TARGET_NAME1 unit_tests)
 SET(TARGET_NAME2 act_test)
 SET(TARGET_NAME3 randtest)
+SET(TARGET_NAME4 detprofile)
 INCLUDE_DIRECTORIES(${PROJECT_SOURCE_DIR}/include)
 LINK_DIRECTORIES(${LIBRARY_OUTPUT_PATH})
 
@@ -12,3 +13,6 @@ TARGET_LINK_LIBRARIES(${TARGET_NAME2} cleanbkz ntl)
 
 ADD_EXECUTABLE(${TARGET_NAME3} randtest.cpp)
 TARGET_LINK_LIBRARIES(${TARGET_NAME3} cleanbkz ntl)
+
+ADD_EXECUTABLE(${TARGET_NAME4} detprofile.cpp)
+TARGET_LINK_LIBRARIES(${TARGET_NAME4} cleanbkz ntl)