From 2bdfcc4fbb425d9ad6af0da3de5aa1562abf1d5c Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Tue, 15 Jan 2019 13:23:18 +0100
Subject: [PATCH] added ExponentialDistribution
---
Framework/Random/CMakeLists.txt | 1 +
Framework/Random/ExponentialDistribution.h | 38 ++++++++++++++++++++++
Framework/Random/UniformRealDistribution.h | 5 ++-
Framework/Random/testRandom.cc | 24 ++++++++++++++
4 files changed, 65 insertions(+), 3 deletions(-)
create mode 100644 Framework/Random/ExponentialDistribution.h
diff --git a/Framework/Random/CMakeLists.txt b/Framework/Random/CMakeLists.txt
index 3607b585..84479e12 100644
--- a/Framework/Random/CMakeLists.txt
+++ b/Framework/Random/CMakeLists.txt
@@ -8,6 +8,7 @@ set (
CORSIKArandom_HEADERS
RNGManager.h
UniformRealDistribution.h
+ ExponentialDistribution.h
)
set (
diff --git a/Framework/Random/ExponentialDistribution.h b/Framework/Random/ExponentialDistribution.h
new file mode 100644
index 00000000..248d785a
--- /dev/null
+++ b/Framework/Random/ExponentialDistribution.h
@@ -0,0 +1,38 @@
+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#ifndef _include_ExponentialDistribution_h
+#define _include_ExponentialDistribution_h
+
+#include <corsika/units/PhysicalUnits.h>
+#include <random>
+
+namespace corsika::random {
+
+ template <class TQuantity>
+ class ExponentialDistribution {
+ using RealType = typename TQuantity::value_type;
+ std::exponential_distribution<RealType> dist{1.};
+
+ TQuantity const fBeta;
+
+ public:
+ ExponentialDistribution(TQuantity beta)
+ : fBeta(beta) {}
+
+ template <class Generator>
+ TQuantity operator()(Generator& g) {
+ return fBeta * dist(g);
+ }
+ };
+
+} // namespace corsika::random
+
+#endif
diff --git a/Framework/Random/UniformRealDistribution.h b/Framework/Random/UniformRealDistribution.h
index 73872928..fd501ab3 100644
--- a/Framework/Random/UniformRealDistribution.h
+++ b/Framework/Random/UniformRealDistribution.h
@@ -1,4 +1,3 @@
-
/**
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
@@ -9,8 +8,8 @@
* the license.
*/
-#ifndef _include_random_distributions_h
-#define _include_random_distributions_h
+#ifndef _include_UniformRealDistribution_h
+#define _include_UniformRealDistribution_h
#include <corsika/units/PhysicalUnits.h>
#include <random>
diff --git a/Framework/Random/testRandom.cc b/Framework/Random/testRandom.cc
index 56d8157a..97923eb0 100644
--- a/Framework/Random/testRandom.cc
+++ b/Framework/Random/testRandom.cc
@@ -14,10 +14,12 @@
#include <corsika/random/RNGManager.h>
#include <corsika/random/UniformRealDistribution.h>
+#include <corsika/random/ExponentialDistribution.h>
#include <corsika/units/PhysicalUnits.h>
#include <iostream>
#include <limits>
#include <random>
+#include <type_traits>
using namespace corsika::random;
@@ -83,3 +85,25 @@ TEST_CASE("UniformRealDistribution") {
CHECK(max / 18_cm == Approx(1.));
}
}
+
+TEST_CASE("ExponentialDistribution") {
+ using namespace corsika::units::si;
+ std::mt19937 rng;
+
+ auto const beta = 15_m;
+
+ corsika::random::ExponentialDistribution dist(beta);
+
+ SECTION("mean") {
+ std::remove_const<decltype(beta)>::type mean = beta * 0;
+
+ int constexpr N = 1'000'000;
+
+ for (int i{0}; i < N; ++i) {
+ decltype(beta) x = dist(rng);
+ mean += x / N;
+ }
+
+ CHECK(mean / beta == Approx(1).margin(1e-2));
+ }
+}
--
GitLab