diff --git a/src/runner_pt.cpp b/src/runner_pt.cpp index c264c4ae42860dc4b266436d9494569dfabd6e1c..168904f1894dbeff5791313a099454178dccc25a 100644 --- a/src/runner_pt.cpp +++ b/src/runner_pt.cpp @@ -110,32 +110,30 @@ std::tuple 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; } diff --git a/src/util.cpp b/src/util.cpp index e01fb944be24ffdbcfeffd756e33289ce5b4421a..1ac1897ad7bbbf8d6e9e1a2654214bbdba0cc3e7 100644 --- a/src/util.cpp +++ b/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]; diff --git a/test/meson.build b/test/meson.build index 0b4807bc9f0e8c5e4c51d3fc10b4c57bf71c96d8..843812c5897a7581cb50c3eaba181bf3404ff264 100644 --- a/test/meson.build +++ b/test/meson.build @@ -1,8 +1,9 @@ -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) diff --git a/test/monotone_interpolator.cpp b/test/monotone_interpolator.cpp new file mode 100644 index 0000000000000000000000000000000000000000..4dddcc09571da05b7df29557dbe6d4459f9d1d6b --- /dev/null +++ b/test/monotone_interpolator.cpp @@ -0,0 +1,47 @@ +#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 x(npoints); + std::vector 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); +}