really actually fix the new optimizer

src/runner_pt.cpp

@@ -110,32 +110,30 @@ std::tuple<double, double> pt_chain::optimize_params() {
 
 
 	double sum{};
+	double efficiency{};
 	for(size_t i = 0; i < comm_barrier.size()-1; i++) {
 		comm_barrier[i] = sum;
 		sum += rejection_est[i]; 
+		efficiency += rejection_est[i]/(1-rejection_est[i]);
 	}
 	comm_barrier[comm_barrier.size()-1] = sum;
 
 	monotonic_interpolator lambda{comm_barrier, params};
 	double convergence{};
-	double nonlinearity = 0;
 
 	for(size_t i = 1; i < params.size()-1; i++) {
 		double new_param = lambda(sum*i/(params.size()-1));
 		double d = (new_param - params[i]);
 		
-		double ideal_comm_barrier = sum*i/(params.size()-1);
-		nonlinearity += pow(comm_barrier[i]-ideal_comm_barrier,2);
 		convergence += d*d;
 		params[i] = new_param;
 	}
 
-	int n = params.size();
-	double nonlinearity_worst = sqrt(n/3. + 1./(6*(n-1))-5./6.);
-	nonlinearity = sqrt(nonlinearity)/nonlinearity_worst;
 	convergence = sqrt(convergence)/params.size();
 
-	return std::tie(nonlinearity, convergence);
+	double round_trip_rate = (1+sum)/(1+efficiency);
+
+	return std::tie(round_trip_rate, convergence);
 }
 
 bool pt_chain::is_done() {
@@ -539,7 +537,6 @@ void runner_pt_master::react() {
 			if(chain_run.weight_ratios[pos] < 0) {
 				double weight;
 				MPI_Recv(&weight, 1, MPI_DOUBLE, target_rank, 0, MPI_COMM_WORLD, &stat);
-				assert(weight >= 0);
 				chain_run.weight_ratios[pos] = weight;
 			}
 		}
@@ -571,11 +568,11 @@ void runner_pt_master::pt_global_update(pt_chain &chain, pt_chain_run &chain_run
 		double w1 = chain_run.weight_ratios[i];
 		double w2 = chain_run.weight_ratios[i + 1];
 
+		double p = std::min(exp(w1+w2),1.);
 		double r = rng_->random_double();
 
-		chain.rejection_rates[i] += 1 - std::min(w1 * w2, 1.);
-		chain.rejection_rate_entries[chain_run.swap_odd]++;
-		if(r < w1 * w2) {
+		chain.rejection_rates[i] += 1 - p;
+		if(r < p) {
 			int rank0{};
 			int rank1{};
 
@@ -596,6 +593,7 @@ void runner_pt_master::pt_global_update(pt_chain &chain, pt_chain_run &chain_run
 		}
 	}
 
+	chain.rejection_rate_entries[chain_run.swap_odd]++;
 	chain_run.swap_odd = !chain_run.swap_odd;
 }
 

src/util.cpp

@@ -54,11 +54,11 @@ double monotonic_interpolator::operator()(double x0) {
 	// replace by binary search if necessary!
 	size_t idx = 0;
 	assert(x0 < x_[x_.size()-1]);
-	while(x0 >= x_[idx]) {
-		idx++;
+	for(; idx < x_.size()-1; idx++) {
+		if((x0-x_[idx])*(x0-x_[idx+1]) <= 0) {
+			break;
+		}
 	}
-	assert(idx >= 1);
-	idx--;
 	assert(idx < x_.size()-1);
 
 	double del = x_[idx+1]-x_[idx];

test/meson.build

 
-t1 = executable('test_duration_parser', 'duration_parser.cpp',
+t1 = executable('tests',
+  ['duration_parser.cpp', 'monotone_interpolator.cpp'],
   dependencies : loadleveller_dep,
   include_directories : include_directories('../src')
 )
 
-test('test duration parser', t1)
+test('tests', t1)
         

test/monotone_interpolator.cpp

+#include "catch.hpp"
+#include "util.h"
+
+using namespace loadl;
+
+double testfun1(double x) {
+	return x+sin(x);
+}
+
+double testfun2(double x) {
+	return -x*x-sin(x*x);
+}
+
+double testfun3(double x) {
+	return x*x*x;
+}
+
+double interpolate(double (*f)(double), double xmax, int npoints, int ncheck) {
+	double dx = xmax/(npoints-1);
+	std::vector<double> x(npoints);
+	std::vector<double> y(npoints);
+
+	for(int i = 0; i < npoints; i++) {
+		x[i] = dx*i;
+		y[i] = f(x[i]);
+	}
+
+	monotonic_interpolator inter(x, y);
+
+	double rms = 0;
+	double dx_check = xmax/(ncheck+1);
+	for(int i = 0; i < ncheck; i++) {
+		double x = i*dx_check;
+		double y = inter(x);
+		rms += (f(x)-y)*(f(x)-y);
+	}
+
+	rms = sqrt(rms)/ncheck;
+	return rms;
+}
+
+TEST_CASE("monotonic interpolator") {
+	REQUIRE(interpolate(testfun1, 10, 30, 100) < 0.0001);
+	REQUIRE(interpolate(testfun2, 10, 150, 200) < 0.001);
+	REQUIRE(interpolate(testfun2, 3, 60, 100) < 0.0001);
+	REQUIRE(interpolate(testfun3, 1, 10, 100) < 0.0005);
+}