diff --git a/Framework/Random/CMakeLists.txt b/Framework/Random/CMakeLists.txt index 3c1947d2508b318eafe6ab426bceac9af7b9e752..230f84f67f4ed6eb02652c5cb10853955f7de776 100644 --- a/Framework/Random/CMakeLists.txt +++ b/Framework/Random/CMakeLists.txt @@ -7,6 +7,7 @@ set ( set ( CORSIKArandom_HEADERS RNGManager.h + UniformRealDistribution.h ) set ( diff --git a/Framework/Random/UniformRealDistribution.h b/Framework/Random/UniformRealDistribution.h new file mode 100644 index 0000000000000000000000000000000000000000..c90ab674ea41976647320a64095a5e8b68d6484e --- /dev/null +++ b/Framework/Random/UniformRealDistribution.h @@ -0,0 +1,32 @@ +#ifndef _include_random_distributions_h +#define _include_random_distributions_h + +#include <corsika/units/PhysicalUnits.h> +#include <random> + +namespace corsika::random { + + template <class TQuantity> + class UniformRealDistribution { + using RealType = typename TQuantity::value_type; + std::uniform_real_distribution<RealType> dist{RealType(0.), RealType(1.)}; + + TQuantity const a, b; + + public: + UniformRealDistribution(TQuantity b) + : a{TQuantity(phys::units::detail::magnitude_tag, 0)} + , b(b) {} + UniformRealDistribution(TQuantity a, TQuantity b) + : a(a) + , b(b) {} + + template <class Generator> + TQuantity operator()(Generator& g) { + return a + dist(g) * (b - a); + } + }; + +} // namespace corsika::random + +#endif diff --git a/Framework/Random/testRandom.cc b/Framework/Random/testRandom.cc index 041a96c21ff241e3a86a903a887e521ff17c9e29..56d8157a74530d404f5d8ba1c8cf8aefe9f75f99 100644 --- a/Framework/Random/testRandom.cc +++ b/Framework/Random/testRandom.cc @@ -1,4 +1,3 @@ - /** * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * @@ -14,7 +13,11 @@ #include <catch2/catch.hpp> #include <corsika/random/RNGManager.h> +#include <corsika/random/UniformRealDistribution.h> +#include <corsika/units/PhysicalUnits.h> #include <iostream> +#include <limits> +#include <random> using namespace corsika::random; @@ -37,3 +40,46 @@ SCENARIO("random-number streams can be registered and retrieved") { } } } + +TEST_CASE("UniformRealDistribution") { + using namespace corsika::units::si; + std::mt19937 rng; + + corsika::random::UniformRealDistribution<LengthType> dist(1_m, 2_m); + + SECTION("range") { + corsika::random::UniformRealDistribution<LengthType> dist(1_m, 2_m); + + LengthType min = + +1_m * std::numeric_limits<typename LengthType::value_type>::infinity(); + LengthType max = + -1_m * std::numeric_limits<typename LengthType::value_type>::infinity(); + + for (int i{0}; i < 1'000'000; ++i) { + LengthType x = dist(rng); + min = std::min(min, x); + max = std::max(max, x); + } + + CHECK(min / 1_m == Approx(1.)); + CHECK(max / 2_m == Approx(1.)); + } + + SECTION("range") { + corsika::random::UniformRealDistribution<LengthType> dist(18_cm); + + LengthType min = + +1_m * std::numeric_limits<typename LengthType::value_type>::infinity(); + LengthType max = + -1_m * std::numeric_limits<typename LengthType::value_type>::infinity(); + + for (int i{0}; i < 1'000'000; ++i) { + LengthType x = dist(rng); + min = std::min(min, x); + max = std::max(max, x); + } + + CHECK(min / 1_m == Approx(0.).margin(1e-3)); + CHECK(max / 18_cm == Approx(1.)); + } +}