From a1169711b8ede816bfa6f3539fdfba6da5400734 Mon Sep 17 00:00:00 2001
From: Kristofer Munsterhjelm <kristofer@munsterhjelm.no>
Date: Mon, 23 Sep 2024 18:30:13 +0200
Subject: [PATCH] Gaussian integral: Show percentile where the optimum is
 located.

---
 src/multiwinner/auxiliary/normal_integral.cc | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/multiwinner/auxiliary/normal_integral.cc b/src/multiwinner/auxiliary/normal_integral.cc
index 0cbc5fc..a26e805 100644
--- a/src/multiwinner/auxiliary/normal_integral.cc
+++ b/src/multiwinner/auxiliary/normal_integral.cc
@@ -52,7 +52,6 @@ double harmonic_integral_monte_carlo(double x_1, double delta) {
 	double sumIIA = 0, sumIIB = 0;
 	double sumIIIA = 0, sumIIIB = 0;
 	double sumIVA = 0, sumIVB = 0;
-	double debug = 0;
 
 	if (x_1 > 0) {
 		throw std::invalid_argument("Harmonic integral MC: x_1 must be "
@@ -61,7 +60,7 @@ double harmonic_integral_monte_carlo(double x_1, double delta) {
 
 	double x_2 = -x_1;
 
-	for (int i = 0; i < maxiter; ++i) {
+	for (size_t i = 0; i < maxiter; ++i) {
 		double x = rnorm().first; // wasteful but wth
 		double contrib = 0;
 		if (dist(x, x_1) < dist(x, x_2)) {
@@ -123,10 +122,10 @@ double inverf(double x) {
 
 int main() {
 
-	double x_1 = -0.2, x_2 = 0.2;
 	double delta = 0.5;
 
 	double recordholder = -1, record = -1;
+	double x_1;
 
 	for (x_1 = -1; x_1 < 0; x_1 += 0.01) {
 		if (harmonic_integral(x_1, delta) > record) {
@@ -146,5 +145,8 @@ int main() {
 	std::cout << "The quality function is then " << record << "\n";
 
 	std::cout << "Compare semi-analytical optimum: " <<
-		sqrt(2) * inverf(1 - (2*delta+3)/(2.0*delta+2)) << "\n";
+		sqrt(2) * inverf(-1/(2*delta+2)) << "\n";
+
+	std::cout << "x_1 at the " << 100 * (2 * delta + 1) /
+		(4 * delta + 4) << " percentile.\n";
 }
\ No newline at end of file