make rng backend modular

meson_options.txt

+option('rng_backend', type : 'combo', choices : ['stl_mt19937', 'stl_mt19937_64', 'internal_mersenne'], description : 'The backend used for the random number generator provided by loadleveller')

src/config.h.in

+#pragma once
+
+#define RNG_BACKEND @rng_backend@

src/mc.cpp

@@ -11,14 +11,14 @@ void mc::write_output(const std::string &) {}
 
 void mc::random_init() {
 	if(param.defined("seed")) {
-		rng.reset(new randomnumbergenerator(param.get<uint64_t>("seed")));
+		rng.reset(new random_number_generator(param.get<uint64_t>("seed")));
 	} else {
-		rng.reset(new randomnumbergenerator());
+		rng.reset(new random_number_generator());
 	}
 }
 
 double mc::random01() {
-	return rng->d();
+	return rng->random_double();
 }
 
 int mc::sweep() const {
@@ -82,7 +82,6 @@ void mc::_write(const std::string &dir) {
 static bool file_exists(const std::string &path) {
 	struct stat buf;
 	return stat(path.c_str(), &buf) == 0;
-
 }
 
 bool mc::_read(const std::string &dir) {

src/mc.h

@@ -19,7 +19,7 @@ private:
 
 protected:
 	parser param;
-	std::unique_ptr<randomnumbergenerator> rng;
+	std::unique_ptr<random_number_generator> rng;
 
 	virtual void init() = 0;
 	virtual void checkpoint_write(const iodump::group &out) = 0;

src/meson.build

+rng_names = {
+  'stl_mt19937' : 'rng_stl<std::mt19937>',
+  'stl_mt19937_64' : 'rng_stl<std::mt19937_64>',
+  'internal_mersenne' : 'rng_internal_mersenne',
+}
+
+conf_data = configuration_data()
+conf_data.set('rng_backend', rng_names[get_option('rng_backend')])
+configure_file(input : 'config.h.in',
+  output : 'config.h',
+  configuration : conf_data)
+
 loadleveller_sources = files([
   'dump.cpp',
   'evalable.cpp',

src/random.cpp

@@ -2,114 +2,28 @@
 
 namespace loadl {
 
-#ifdef MCL_RNG_MT
-randomnumbergenerator::randomnumbergenerator() {
-	seed_ = mtrand_.get_seed();
+void rng_internal_mersenne::set_seed(uint64_t seed) {
+	mtrand_.seed(seed);
 }
 
-randomnumbergenerator::randomnumbergenerator(uint64_t seed) : mtrand_{seed}, seed_{seed} {}
-
-void randomnumbergenerator::checkpoint_write(const iodump::group &d) {
-	std::vector<uint64_t> randState;
-	mtrand_.save(randState);
-	// if you implement your own random number generator backend,
-	// make sure to keep this file layout for scripts that like to read the seed.
-	d.write("state", randState);
-	d.write("seed", seed_);
-
-	// This is for future compatibility. B) You can add RNG implementations without
-	// breaking compatibility to old runs.
-	d.write("type", static_cast<uint32_t>(RNG_MersenneTwister));
-}
-
-void randomnumbergenerator::checkpoint_read(const iodump::group &d) {
-	std::vector<uint64_t> randState;
-	d.read("state", randState);
-	d.read("seed", seed_);
-	mtrand_.load(randState);
-}
-
-double randomnumbergenerator::norm() { // adjusted from https://de.wikipedia.org/wiki/Polar-Methode
-	/* for now, I do not need the extra speed coming from this.
-	 * If you add hidden state like this you have to change the dump format
-	 * (breaking compatibility with your old runs)
-	   if(norm_is_stored) {
-	    norm_is_stored=false;
-	    return norm_stored;
+void rng_internal_mersenne::backend_checkpoint_write(const iodump::group &d) {
+	std::vector<uint64_t> rand_state;
+	mtrand_.save(rand_state);
+	d.write("state", rand_state);
 }
-	*/
-	double u, v, q, p;
-
-	do {
-		u = 2.0 * d() - 1;
-		v = 2.0 * d() - 1;
-		q = u * u + v * v;
-	} while(q <= 0.0 || q >= 1.0);
-
-	p = sqrt(-2 * log(q) / q);
-	/*
-	norm_stored = u * p;
-	norm_is_stored = true;
-	*/
-	return v * p;
-}
-
-#endif
 
-#ifdef MCL_RNG_SPRNG_4
-randomnumbergenerator::randomnumbergenerator() {
-	int gtype = 4; // Multipl. Lagg. Fib.
-	ptr = SelectType(gtype);
-	seed_ = make_sprng_seed();
-	ptr->init_sprng(0, 1, seed_, SPRNG_DEFAULT);
+void rng_internal_mersenne::backend_checkpoint_read(const iodump::group &d) {
+	std::vector<uint64_t> rand_state;
+	d.read("state", rand_state);
+	mtrand_.load(rand_state);
 }
 
-randomnumbergenerator::randomnumbergenerator(luint seed) {
-	int gtype = 4; // Multipl. Lagg. Fib.
-	ptr = SelectType(gtype);
-	seed_ = seed;
-	ptr->init_sprng(0, 1, seed_, SPRNG_DEFAULT);
+double rng_internal_mersenne::random_double() {
+	return mtrand_.randDblExc(1);
 }
 
-void randomnumbergenerator::write(odump &d) {
-	char *buffer;
-	int buffersize = ptr->pack_sprng(&buffer);
-	d.write(buffersize);
-	d.write(buffer, buffersize);
-	d.write(seed_);
-	delete buffer;
+int rng_internal_mersenne::random_integer(int bound) {
+	return mtrand_.randInt(bound - 1);
 }
 
-void randomnumbergenerator::read(idump &d) {
-	char *buffer;
-	int buffersize;
-	d.read(buffersize);
-	buffer = new char[buffersize];
-	d.read(buffer, buffersize);
-	d.read(seed_);
-	ptr->unpack_sprng(buffer);
-	delete buffer;
-}
-#endif
-
-#ifdef MCL_RNG_BOOST
-void randomnumbergenerator::write(odump &d) {}
-
-void randomnumbergenerator::read(idump &d) {}
-
-randomnumbergenerator::randomnumbergenerator() {
-	rng = new boost::mt19937;
-	val = new boost::uniform_real<>(0.0, 1.0);
-
-	seed_ = 5489;
-}
-
-randomnumbergenerator::randomnumbergenerator(int seed) {
-	rng = new boost::mt19937;
-	rng->seed(seed);
-	val = new boost::uniform_real<>(0.0, 1.0);
-
-	seed_ = seed;
-}
-#endif
 }

src/random.h

 #pragma once
 
-#include "dump.h"
+#include "config.h"
 
-#define MCL_RNG_MT
+#include "MersenneTwister.h"
+#include "dump.h"
+#include <random>
+#include <sstream>
+#include <typeinfo>
 
-// everything but MCL_RNG_MT is old and needs to be ported to the new api if you want to use it.
-// in that case, also get rid of the macros and use templates to switch between different generators
-// or so.
+// A whole bunch of template magic to make the random backend modular.
+// We don’t want a vtable on such a performance critical object, otherwise I
+// would make it settable at runtime.
 
 namespace loadl {
-enum rng_type {
-	RNG_MersenneTwister = 0,
+
+template<class base>
+class rng_base : public base {
+private:
+	static uint32_t get_rank() {
+		int initialized;
+		MPI_Initialized(&initialized);
+		if(initialized) {
+			int rank;
+			MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+			return rank;
+		}
+		return 0;
+	}
+
+	uint64_t seed_;
+
+public:
+	rng_base() {
+		std::seed_seq seq = {static_cast<uint64_t>(time(0)), static_cast<uint64_t>(clock()),
+		                     static_cast<uint64_t>(get_rank())};
+		std::vector<uint64_t> seed(1);
+		seq.generate(seed.begin(), seed.end());
+		static_cast<base *>(this)->set_seed(seed[0]);
 	};
+
+	rng_base(uint64_t seed) {
+		static_cast<base *>(this)->set_seed(seed);
 	}
 
-#ifdef MCL_RNG_MT
+	std::string type() const {
+		return typeid(base).name();
+	}
 
-#include "MersenneTwister.h"
-namespace loadl {
+	void checkpoint_write(const iodump::group &d) {
+		static_cast<base *>(this)->backend_checkpoint_write(d);
+		d.write("seed", seed_);
+		d.write("type", type());
+	}
 
-class randomnumbergenerator {
-public:
-	randomnumbergenerator();
-	randomnumbergenerator(uint64_t seed);
-	void checkpoint_write(const iodump::group &dump_file);
-	void checkpoint_read(const iodump::group &dump_file);
-	uint64_t seed() {
-		return seed_;
+	void checkpoint_read(const iodump::group &d) {
+		static_cast<base *>(this)->backend_checkpoint_read(d);
+		d.read("seed", seed_);
+		std::string read_type;
+		d.read("type", read_type);
+		if(read_type != type()) {
+			throw std::runtime_error(
+			    fmt::format("dump file was created using rng backend '{}' while this version of "
+			                "loadleveller was built with backend '{}'",
+			                read_type, type()));
+		}
 	}
-	double d(double supp = 1) {
-		return mtrand_.randDblExc(supp);
-	} // in (0,supp)
-	int i(int bound) {
-		return mtrand_.randInt(bound - 1);
-	} // in [0,bound)
 
-	double norm(); // normal distribution, mean = 0, std = 1
+	// implemented by base:
+	//
+
+	uint32_t rand_integer(uint32_t bound); // in [0,bound)
+	double rand_double();                  // in [0,1]
+};
+
+// based on a dinosaur code in the MersenneTwister.h header
+class rng_internal_mersenne {
+public:
+	void backend_checkpoint_write(const iodump::group &dump_file);
+	void backend_checkpoint_read(const iodump::group &dump_file);
+	void set_seed(uint64_t seed);
+
+	double random_double();
+	int random_integer(int bound);
 
 private:
 	MTRand mtrand_;
-	uint64_t seed_;
 };
-}
 
-#endif
-
-#ifdef MCL_RNG_SPRNG_4
-#include "sprng_cpp.h"
-class randomnumbergenerator {
+// based on the c++ stl implementation
+template<typename stl_engine>
+class rng_stl {
 public:
-	randomnumbergenerator();
-	randomnumbergenerator(uint64_t seed);
-	~randomnumbergenerator() {
-		delete ptr;
+	void backend_checkpoint_write(const iodump::group &dump_file) {
+		std::stringstream buf;
+		buf << generator_;
+		buf << real_dist_;
+		dump_file.write("rng_state", buf.str());
 	}
-	void write(odump &);
-	void read(idump &);
-	uint64_t seed() {
-		return my_seed;
+	void backend_checkpoint_read(const iodump::group &dump_file) {
+		std::string state;
+		dump_file.read("rng_state", state);
+		std::stringstream buf(state);
+		buf >> generator_;
+		buf >> real_dist_;
 	}
-	double d() {
-		return ptr->sprng();
-	} // in (0,1)
-	double d(double supp) {
-		return supp * d();
-	} // in (0,supp)
-	int i(int bound) {
-		return int(bound * d());
-	} // in [0,bound)
 
-private:
-	Sprng *ptr;
-	uint64_t my_seed;
-};
-#endif
+	void set_seed(uint64_t seed) {
+		generator_.seed(seed);
+	}
 
-#ifdef MCL_RNG_BOOST
-#include <boost/random.hpp>
-class randomnumbergenerator {
-public:
-	randomnumbergenerator();
-	randomnumbergenerator(int seed);
-	void write(odump &);
-	void read(idump &);
-	uint64_t seed() {
-		return my_seed;
+	double random_double() {
+		return real_dist_(generator_);
+	}
+
+	uint32_t random_integer(uint32_t bound) {
+		std::uniform_int_distribution<uint32_t> int_dist(0, bound);
+		return int_dist(generator_);
 	}
-	double d() {
-		return (*val)(*rng);
-	} // returns a value between 0 (excl) and 1 (excl).
-	double d(double supp) {
-		return supp * d();
-	} // returns a value between 0 (excl) and supp (excl.)
-	int i(int bound) {
-		return int(bound * d());
-	} // returns a value between 0 (incl) and bound (excl.)
 
 private:
-	boost::mt19937 *rng;
-	boost::uniform_real<> *val;
-	uint64_t my_seed;
+	stl_engine generator_;
+	std::uniform_real_distribution<double> real_dist_{0, 1};
 };
-#endif
+
+// RNG_BACKEND is a macro set by the build system. If you add backends and you can help it,
+// avoid using huge blocks of #ifdefs as it will lead to dead code nobody uses for 10 years.
+// using random_number_generator = rng_base<RNG_BACKEND, #RNG_BACKEND>;
+
+using random_number_generator = rng_base<RNG_BACKEND>;
+}