Import GCC's std::normal_distribution to make things portable.

CMakeLists.txt

@@ -21,6 +21,9 @@ PROJECT (SampleFlow
 #########################################
 ### Set up compiler flags and input paths
 
+message(STATUS "Using compiler ${CMAKE_CXX_COMPILER_ID}, version ${CMAKE_CXX_COMPILER_VERSION}")
+
+
 #########################################
 # Also make sure we link with the threads library in question.
 # We need this when linking C++ stuff that uses threads.
@@ -114,13 +117,19 @@ INSTALL(FILES "${PROJECT_BINARY_DIR}/${PROJECT_NAME}config.cmake"
 
 # Now also define a library that consists of the C++20 modules,
 # assuming the compiler and generator supports this:
-if (${CMAKE_CXX_COMPILER_ID} STREQUAL "Clang"
-    AND
-    ${CMAKE_CXX_COMPILER_VERSION} VERSION_GREATER_EQUAL 16
+if (((${CMAKE_CXX_COMPILER_ID} STREQUAL "Clang"
+      AND
+      ${CMAKE_CXX_COMPILER_VERSION} VERSION_GREATER_EQUAL 16)
+     OR
+     (${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU"
+      AND
+      ${CMAKE_CXX_COMPILER_VERSION} VERSION_GREATER_EQUAL 14))   # does not actually work in practice
     AND
     ${CMAKE_GENERATOR} STREQUAL "Ninja")
   message(STATUS "Enabling the use of C++20-style modules")
   set(SAMPLEFLOW_BUILD_MODULE "ON")
+else()
+  message(STATUS "Not using C++20-style modules")
 endif()
 
 if (SAMPLEFLOW_BUILD_MODULE)

tests/adaptive_mh_01.cc

@@ -25,6 +25,8 @@
 
 #include <eigen3/Eigen/Dense>
 
+#include "tests.h"
+
 #ifndef SAMPLEFLOW_TEST_WITH_MODULE
 #  include <sampleflow/producers/metropolis_hastings.h>
 #  include <sampleflow/consumers/covariance_matrix.h>
@@ -78,7 +80,7 @@ std::pair<SampleType,double> perturb_adaptive (const SampleType &x,
   SampleType random_vector;
   for (unsigned int i=0; i<random_vector.size(); ++i)
     random_vector(i) = 2.4/std::sqrt(1.*x.size()) *
-                       std::normal_distribution<double>(0,1)(rng);
+                       SampleFlow::Testing::NormalDistribution<double>(0,1)(rng);
 
   const SampleType y = x + LLt.matrixL() * random_vector;
 

tests/auto_covariance_matrix_05.cc

@@ -29,6 +29,8 @@
 
 #include <eigen3/Eigen/Dense>
 
+#include "tests.h"
+
 #ifndef SAMPLEFLOW_TEST_WITH_MODULE
 #  include <sampleflow/producers/metropolis_hastings.h>
 #  include <sampleflow/consumers/acceptance_ratio.h>
@@ -64,7 +66,7 @@ std::pair<SampleType,double> perturb (const SampleType &x)
   // Perturb the current sample using a Gaussian distribution
   // around the current point with standard deviation 1.5.
   static std::mt19937 rng;
-  std::normal_distribution<double> distribution(0., 1.5);
+  SampleFlow::Testing::NormalDistribution<double> distribution(0., 1.5);
 
   SampleType y = x;
   for (auto &el : y)

tests/auto_covariance_trace_05.cc

@@ -29,6 +29,8 @@
 
 #include <eigen3/Eigen/Dense>
 
+#include "tests.h"
+
 #ifndef SAMPLEFLOW_TEST_WITH_MODULE
 #  include <sampleflow/producers/metropolis_hastings.h>
 #  include <sampleflow/consumers/acceptance_ratio.h>
@@ -53,7 +55,7 @@ std::pair<SampleType,double> perturb (const SampleType &x)
   // Perturb the current sample using a Gaussian distribution
   // around the current point with standard deviation 1.5.
   static std::mt19937 rng;
-  std::normal_distribution<double> distribution(0., 1.5);
+  SampleFlow::Testing::NormalDistribution<double> distribution(0., 1.5);
 
   SampleType y = x;
   for (auto &el : y)

tests/covariance_matrix_11.cc

@@ -31,6 +31,8 @@
 
 #include <eigen3/Eigen/Dense>
 
+#include "tests.h"
+
 #ifndef SAMPLEFLOW_TEST_WITH_MODULE
 #  include <sampleflow/producers/metropolis_hastings.h>
 #  include <sampleflow/consumers/covariance_matrix.h>
@@ -65,7 +67,7 @@ std::pair<SampleType,double> perturb (const SampleType &x)
   static std::mt19937 rng;
   SampleType random_vector;
   for (unsigned int i=0; i<random_vector.size(); ++i)
-    random_vector(i) = std::normal_distribution<double>(0,1)(rng);
+    random_vector(i) = SampleFlow::Testing::NormalDistribution<double>(0,1)(rng);
 
   const SampleType y = (LLt.matrixL()) * random_vector + x;
 

tests/delayed_rejection_mh_producer_03.cc

@@ -22,6 +22,8 @@
 #include <random>
 #include <cmath>
 
+#include "tests.h"
+
 #ifndef SAMPLEFLOW_TEST_WITH_MODULE
 #  include <sampleflow/producers/delayed_rejection_mh.h>
 #  include <sampleflow/filters/conversion.h>
@@ -43,7 +45,7 @@ double log_likelihood (const SampleType &x)
 std::pair<SampleType,double> perturb (const SampleType &x, const std::vector<SampleType> &)
 {
   static std::mt19937 rng;
-  std::normal_distribution<double> distribution(0, 1);
+  SampleFlow::Testing::NormalDistribution<double> distribution(0, 1);
   const double perturbation = distribution(rng);
   const SampleType x_tilde = x + perturbation;
   return {x_tilde, 1.0};

tests/my_triangle.h

@@ -20,6 +20,7 @@
 #include <limits>
 #include <valarray>
 
+#include "tests.h"
 
 /**
  * A class that serves as a custom sample class that is not a vector
@@ -54,11 +55,12 @@ std::ostream &operator<<(std::ostream &os, const MyTriangle &tri)
 }
 
 
+
 inline
 std::pair<MyTriangle, double> perturb(const MyTriangle &sample)
 {
   static std::mt19937 gen;
-  std::normal_distribution<> d {0, 1};
+  SampleFlow::Testing::NormalDistribution<> d {0, 1};
   double side_a = sample.side_lengths[0] + d(gen);
   double side_b = sample.side_lengths[1] + d(gen);
   double side_c = sample.side_lengths[2] + d(gen);

tests/tests.h

@@ -0,0 +1,106 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the SampleFlow authors.
+//
+// This file is part of the SampleFlow library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef SAMPLEFLOW_TESTS_H
+#define SAMPLEFLOW_TESTS_H
+
+
+namespace SampleFlow
+{
+  namespace Testing
+  {
+    /**
+     * Whereas the random number generators are portable between
+     * compilers, the random number distributions are not. This class
+     * makes sure we remain compatible at least between GCC and Clang
+     * by creating a class that, given a random number generator,
+     * returns normally distributed random floating point numbers.
+     *
+     * The class is a drop-in replacement for std::normal_distribution.
+     */
+    template<typename RealType = double>
+    class NormalDistribution
+    {
+        static_assert(std::is_floating_point_v<RealType>,
+                      "result_type must be a floating point type");
+
+      public:
+        using result_type = RealType;
+
+
+      public:
+        explicit
+        NormalDistribution(result_type mu,
+                           result_type sigma)
+          : mean(mu), stddev(sigma)
+        { }
+
+        template<typename RNG>
+        result_type
+        operator()(RNG &rng);
+
+      private:
+        result_type mean, stddev;
+        result_type saved_value = 0;
+        bool        saved_value_available = false;
+    };
+
+
+
+    template<typename RealType>
+    template<typename RNG>
+    typename NormalDistribution<RealType>::result_type
+    NormalDistribution<RealType>::
+    operator()(RNG &rng)
+    {
+      result_type ret;
+
+      auto get_real = [&rng]()
+      {
+        return std::generate_canonical<RealType,
+               std::numeric_limits<RealType>::digits,
+               RNG>(rng);
+      };
+
+      if (saved_value_available)
+        {
+          saved_value_available = false;
+          ret = saved_value;
+        }
+      else
+        {
+          result_type x, y, r2;
+          do
+            {
+              x = result_type(2.0) * get_real() - 1.0;
+              y = result_type(2.0) * get_real() - 1.0;
+              r2 = x * x + y * y;
+            }
+          while (r2 > 1.0 || r2 == 0.0);
+
+          const result_type mult = std::sqrt(-2 * std::log(r2) / r2);
+          ret = y * mult;
+          saved_value = x * mult;
+          saved_value_available = true;
+
+        }
+      return ret * stddev + mean;
+    }
+  }
+}
+
+
+#endif