random number generation from a distribution defined by a histogram

cmake/doctest.cmake

@@ -1,17 +1,35 @@
 # Downloads and compiles DocTest unit testing framework
 include(FetchContent)
-set(FETCHCONTENT_QUIET ON)
-set(FETCHCONTENT_UPDATES_DISCONNECTED ON)
+#set(FETCHCONTENT_QUIET ON)
+#set(FETCHCONTENT_UPDATES_DISCONNECTED ON)
 
-set(DOCTEST_WITH_TESTS OFF CACHE BOOL "doctest tests and examples")
-set(DOCTEST_WITH_MAIN_IN_STATIC_LIB ON CACHE BOOL "enable doctest_with_main")
-set(DOCTEST_NO_INSTALL OFF CACHE BOOL "Skip the installation process")
-#set(DOCTEST_USE_STD_HEADERS OFF CACHE BOOL "Use std headers")
+set(CMAKE_POLICY_VERSION_MINIMUM 3.10) # turn off the CMake warnings
 
 FetchContent_Declare(doctest
         GIT_REPOSITORY https://github.com/doctest/doctest.git
-        GIT_TAG v2.4.11    # "main" for latest
+        GIT_TAG v2.4.12   # "main" for latest
         GIT_SHALLOW TRUE  # download specific revision only (git clone --depth 1)
         GIT_PROGRESS TRUE # show download progress in Ninja
-        USES_TERMINAL_DOWNLOAD TRUE)
+        EXCLUDE_FROM_ALL ON # don't build if not used
+        FIND_PACKAGE_ARGS 2.4.12)
+
+set(DOCTEST_WITH_TESTS OFF CACHE BOOL "Build tests/examples")
+set(DOCTEST_WITH_MAIN_IN_STATIC_LIB ON CACHE BOOL "Build a static lib for doctest::doctest_with_main")
+set(DOCTEST_NO_INSTALL OFF CACHE BOOL "Skip the installation process")
+set(DOCTEST_USE_STD_HEADERS OFF CACHE BOOL "Use std headers")
+
 FetchContent_MakeAvailable(doctest)
+
+if(doctest_FOUND) # find_package
+    message(STATUS "Found doctest: ${doctest_DIR}")
+else(doctest_FOUND) # FetchContent
+    message(STATUS "Fetched doctest: ${doctest_SOURCE_DIR}")
+endif(doctest_FOUND)
+
+if (TARGET doctest::doctest)
+    message(STATUS "    Available target: doctest::doctest")
+endif ()
+
+if (TARGET doctest::doctest_with_main)
+    message(STATUS "    Available target: doctest::doctest_with_main")
+endif ()

examples/histogram.csv

@@ -0,0 +1,8 @@
+#offset,width,height
+2000, 1000,  10
+3000, 1000, 150
+4000, 1000,   0
+5000, 1000, 220
+6000, 1000,  40
+7000, 1000, 300
+9000, 1000, 280

src/CMakeLists.txt

@@ -1,14 +1,18 @@
 add_custom_target(data
         COMMAND ${CMAKE_COMMAND} -E copy ${PROJECT_SOURCE_DIR}/table_input.csv ${CMAKE_CURRENT_BINARY_DIR}/table_input.csv
-        DEPENDS ${PROJECT_SOURCE_DIR}/table_input.csv
-        BYPRODUCTS table_input.csv)
+        COMMAND ${CMAKE_COMMAND} -E copy ${PROJECT_SOURCE_DIR}/examples/histogram.csv ${CMAKE_CURRENT_BINARY_DIR}/histogram.csv
+        DEPENDS ${PROJECT_SOURCE_DIR}/table_input.csv ${PROJECT_SOURCE_DIR}/examples/histogram.csv
+        BYPRODUCTS table_input.csv histogram.csv)
 
 add_library(errors OBJECT errors.cpp)
 
 add_library(table SHARED table.cpp)
 target_link_libraries(table PRIVATE errors)
 add_dependencies(table data)
 
+add_library(generator SHARED generator.cpp)
+add_dependencies(generator data)
+
 if (UPPAALLIBS_WITH_TESTS)
     add_executable(test_table test_table.cpp)
     target_link_libraries(test_table PRIVATE doctest::doctest_with_main)
@@ -20,4 +24,9 @@ if (UPPAALLIBS_WITH_TESTS)
         target_link_libraries(test_errors PRIVATE errors doctest::doctest_with_main)
         add_test(NAME test_errors COMMAND test_errors)
     endif()
+
+    add_executable(test_generator test_generator.cpp)
+    target_link_libraries(test_generator PRIVATE generator doctest::doctest_with_main)
+    add_test(NAME test_generator COMMAND test_generator)
+
 endif (UPPAALLIBS_WITH_TESTS)

src/dynlib.h

@@ -1,6 +1,8 @@
 #ifndef UPPAAL_LIBS_DYNLIB_H
 #define UPPAAL_LIBS_DYNLIB_H
 
+/// Macros for exporting library symbols
+
 #if defined(_WIN32)
 #define C_PUBLIC extern "C" __declspec(dllexport)
 #elif defined(__linux__)

src/errors.cpp

@@ -10,17 +10,26 @@
 #include <cstdarg>	// va_list
 #include <cerrno>	// errno
 
-static auto error_own = false;	// do we own the error file?
-static auto error_path = std::string{"error.log"};
+static bool& error_own()
+{
+	static auto own = false;	// do we own the error file?
+	return own;
+}
+
+static std::string& error_path()
+{
+	static auto path = std::string{"error.log"};
+	return path;
+}
 
 C_PUBLIC int set_error_path(const char* path)
 {
-	error_path = path;
-	error_own = false;
+	error_path() = path;
+	error_own() = false;
 	return 0;
 }
 
-C_PUBLIC const char* get_error_path() { return error_path.c_str(); }
+C_PUBLIC const char* get_error_path() { return error_path().c_str(); }
 
 FILE* open_error_file()
 {
@@ -39,11 +48,11 @@ FILE* open_error_file()
 		file = nullptr;
 	}
 #else
-	if (error_own) {
+	if (error_own()) {
 		file = std::fopen(get_error_path(), "a");
 	} else {
 		file = std::fopen(get_error_path(), "w");
-		error_own = true;
+		error_own() = true;
 	}
 	if (file == nullptr) {
 		fprintf(stderr, "error while opening %s: %d\n", path, errno);

src/generator.cpp

@@ -0,0 +1,154 @@
+#include "generator.hpp"
+
+#include <vector>
+#include <fstream>
+#include <random>
+#include <cassert>
+
+static std::istream& skip_comments(std::istream& is)
+{
+	auto line = std::string{};
+	while (is) {
+		if (is.peek() == '#')
+			std::getline(is, line);
+		else if (std::isspace(is.peek()))
+			is.get();
+		else
+			break;
+	}
+	return is;
+}
+
+static std::istream& skip_separator(std::istream& is)
+{
+	while (is && std::isspace(is.peek()))
+		is.get();
+	if (is) {
+		switch (is.peek()) {
+		case ',':
+		case ';':
+			is.get(); break;
+		default: is.setstate(std::ios::badbit);
+		}
+	} else {
+		is.setstate(std::ios::badbit);
+	}
+	return is;
+}
+
+
+struct Bar
+{
+	double offset;
+	double width; // remove if distribution is discrete
+	double height;
+};
+
+static std::istream& operator>>(std::istream& is, Bar& bar)
+{
+	is >> skip_comments >> bar.offset >> skip_separator >> bar.width >> skip_separator >> bar.height;
+	if (bar.width < 0 || bar.height < 0)
+		is.setstate(std::ios::badbit);
+	return is;
+}
+
+static double sum(const std::vector<Bar>& bars)
+{
+	auto res = 0.0;
+	for (const auto& bar: bars)
+		res += bar.width * bar.height;
+	return res;
+}
+
+static std::mt19937& generator()
+{
+	static auto gen = std::mt19937{std::random_device{}()};
+	return gen;
+}
+
+class Histogram
+{
+	std::vector<Bar> bars;
+	double total;
+public:
+	Histogram(std::vector<Bar> bars): bars{std::move(bars)}, total{sum(this->bars)} {}
+	Histogram(): Histogram{{}} {}
+	void clear() { bars.clear(); total = 0; }
+	void add(double offset, double width, double height) {
+		if (not bars.empty() and offset < bars.back().offset + bars.back().width)
+			throw std::logic_error{"Offsets must be increasing monotonically and cannot overlap"};
+		bars.emplace_back(offset, width, height);
+		total += width * height;
+	}
+	unsigned int size() const { return bars.size(); }
+	double sample() {
+		if (bars.empty())
+			return 0;
+		auto dist = std::uniform_real_distribution<double>{0, total};
+		auto pick = dist(generator());
+		auto i = 0u;
+		auto weight = bars[i].width * bars[i].height;
+		while (weight < pick) {
+			pick -= weight;
+			++i;
+			assert(i < bars.size());
+			weight = bars[i].width * bars[i].height;
+		}
+		return bars[i].offset + pick / bars[i].height;
+	}
+};
+
+static Histogram& histogram()
+{
+	static auto hist = Histogram{};
+	return hist;
+}
+
+C_PUBLIC int init_generator(unsigned int seed)
+{
+	generator() = std::mt19937{seed};
+	return 1;
+}
+
+/// Loads a histogram from a CSV file
+C_PUBLIC int load_histogram(const char* path)
+{
+	auto last_offset = -std::numeric_limits<double>::infinity();
+	auto bars = std::vector<Bar>{};
+	auto is = std::ifstream{path};
+	auto bar = Bar{};
+	while (is >> bar) {
+		if (last_offset <= bar.offset) {
+			bars.push_back(bar);
+			last_offset = bar.offset + bar.width;
+		} else
+			is.setstate(std::ios::badbit);
+	}
+	histogram() = Histogram{std::move(bars)};
+	return histogram().size();
+}
+
+/// Clears the histogram
+C_PUBLIC int clear_histogram()
+{
+	histogram().clear();
+	return 1;
+}
+
+/// Adds a histogram bar
+C_PUBLIC int add_bar(double offset, double width, double height)
+{
+	try {
+		histogram().add(offset, width, height);
+		return 1;
+	} catch (std::exception& e) {
+		return 0;
+	}
+}
+
+
+/// generates random values from bars
+C_PUBLIC double generate()
+{
+	return histogram().sample();
+}

src/generator.hpp

@@ -0,0 +1,24 @@
+#ifndef UPPAAL_LIBS_GENERATOR_HPP
+#define UPPAAL_LIBS_GENERATOR_HPP
+
+#include "dynlib.h"
+
+/** Initializes the pseudo random number generator with specific seed
+ * @returns 1 upon success */
+C_PUBLIC int init_generator(unsigned int seed);
+
+/** Loads a histogram from a CSV file
+ * @returns the number of successfully loaded data rows
+ */
+C_PUBLIC int load_histogram(const char* path);
+
+/// Clears the histogram
+C_PUBLIC int clear_histogram();
+
+/// Adds a histogram bar
+C_PUBLIC int add_bar(double offset, double width, double height);
+
+/// Generates random values from a previously loaded histogram
+C_PUBLIC double generate();
+
+#endif	// UPPAAL_LIBS_GENERATOR_HPP

src/test_generator.cpp

@@ -0,0 +1,56 @@
+#include "generator.hpp"
+
+#include <doctest/doctest.h>
+
+#include <vector>
+#include <cmath> // round
+
+TEST_CASE("Blank histogram")
+{
+	for (auto i = 0; i < 10; ++i) {
+		const auto n = generate();
+		CHECK(n == 0);
+	}
+}
+
+static const auto approx = doctest::Approx{0}.epsilon(0.05); // within 5% margin
+
+TEST_CASE("Histogram")
+{
+	const auto init = init_generator(42*42*42); // make a predictable sequence
+	REQUIRE(init == 1);
+	SUBCASE("From a file")
+	{
+		const auto rows = load_histogram("histogram.csv");
+		REQUIRE(rows == 7);
+	}
+	SUBCASE("From scratch")
+	{   // same as in "histogram_input.csv"
+		CHECK(clear_histogram() == 1);
+		CHECK(add_bar(2000, 1000,  10) == 1);
+		CHECK(add_bar(3000, 1000, 150) == 1);
+		CHECK(add_bar(4000, 1000,   0) == 1);
+		CHECK(add_bar(5000, 1000, 220) == 1);
+		CHECK(add_bar(6000, 1000,  40) == 1);
+		CHECK(add_bar(7000, 1000, 300) == 1);
+		CHECK(add_bar(9000, 1000, 280) == 1);
+	}
+	auto bars = std::vector<size_t>(8, 0);
+	constexpr auto N = 1000;
+	for (auto i = 0; i < 1000*N; ++i) {
+		const auto n = generate();
+		CHECK(2000 <= n);
+		CHECK(n < 10000);
+		const auto bar = static_cast<size_t>((n-2000) / 1000);
+		REQUIRE(bar < bars.size());
+		++bars[bar];
+	}
+	CHECK(bars[0] == approx(10.*N));
+	CHECK(bars[1] == approx(150.*N));
+	CHECK(bars[2] == 0);
+	CHECK(bars[3] == approx(220.*N));
+	CHECK(bars[4] == approx(40.*N));
+	CHECK(bars[5] == approx(300.*N));
+	CHECK(bars[6] == 0);
+	CHECK(bars[7] == approx(280.*N));
+}