Merge pull request #8 from swig-fortran/random

Extend and improve flc_random

Author: sethrj

CMakeLists.txt

@@ -6,7 +6,7 @@
 #---------------------------------------------------------------------------#
 
 cmake_minimum_required(VERSION 3.8)
-project(Flibcpp VERSION 0.2.1 LANGUAGES CXX Fortran)
+project(Flibcpp VERSION 0.3.0 LANGUAGES CXX Fortran)
 
 #---------------------------------------------------------------------------#
 # OPTIONS

doc/modules/algorithm.rst

@@ -206,12 +206,12 @@ shuffle
 -------
 
 The "shuffle" subroutine depends on the :ref:`modules_random` module so that it
-can use the supported random number generator to randomly reorder an array.
+can use the default random number generator to randomly reorder an array.
 
 Example::
 
   use flc_algorithm, only : shuffle
-  use flc_random, only : Engine
+  use flc_random, only : Engine => MersenneEngine4
   implicit none
   integer :: i
   integer, dimension(8) :: iarr = (/ ((i), i = -4, 3) /)

doc/modules/random.rst

@@ -9,73 +9,106 @@ Random
 ******
 
 The ``flc_random`` module contains advanced pseudo-random number generation
-from the `<random>`_ C++ header. Currently
-random number generation is hardwired to the ``std::mt19937_64`` 64-bit
-Mersenne Twister by Matsumoto and Nishimura.
+from the `<random>`_ C++ header.
 
 The C++11 way of generating random numbers takes some getting used to. Rather
 than having a single global function that returns a random real number in the
 range :math:`[0,1)`, C++11 has independent *engine* objects that generate
 streams of random bits. Different *distribution* objects convert those bits
 into samples of a distribution.
 
+Although C++11 defines a dizzying variety of random number engines,
+``flc_random`` wraps just two: the 32- and 64-bit `Mersenne Twister`
+algorithms. The 32-bit version (``MersenneEngine4``) is currently the only
+engine type that can be used with distributions and algorithms.
+
 Flibcpp wraps distribution objects as independent subroutines. Each subroutine
 accepts the constructor parameters of the distribution, the random engine, and
 a target Fortran array to be filled with random numbers.
 
 Generating an array with 10 normally-distributed reals with a mean of 8 and a
 standard deviation of 2 is done as such::
 
-    use flc_random, only : Engine, normal_distribution
+    use flc_random, only : Engine => MersenneEngine4, normal_distribution
     real(C_DOUBLE), dimension(20) :: arr
     type(Engine) :: rng
 
     rng = Engine()
     call normal_distribution(8.0d0, 2.0d0, rng, arr)
+    call rng%release()
 
 .. _<random> : https://en.cppreference.com/w/cpp/numeric/random
+.. _Mersenne Twister : https://en.wikipedia.org/wiki/Mersenne_Twister
+
+Engines
+=======
+
+The two Mersenne twister engines in ``flc_random`` return different-sized
+integers per call:
+
+ - ``MersenneEngine4``: each invocation returns a 32-bit ``integer(4)``
+ - ``MersenneEngine8``: each invocation returns a 64-bit ``integer(8)``
+
+Engines can be constructed using one of two interface functions: the
+argument-free ``MersenneEngine4()`` uses the default seed, and the engine takes
+a single argument ``MersenneEngine4(1234567)`` with the seed. Alternatively,
+the seed can be set (or reset) using the ``seed()`` type-bound procedure.
 
-Engine
-======
+Generally, engines are used with distributions (described below). However, if
+necessary, individual randomly generated values can be obtained by calling
+the ``next()`` type-bound procedure.
 
-The ``Engine`` is a derived type that wraps the Mersenne Twister PRNG engine.
-Its interface is::
+.. warning:: C++ generates *unsigned* integers with entropy in every bit. This
+   means that the integers obtained from ``engine%next()``, reinterpreted as
+   signed Fortran integers, may be negative.
 
-    type, public :: Engine
-    contains
-     procedure :: seed => swigf_Engine_seed
-     procedure :: discard => swigf_Engine_discard
-     procedure :: release => delete_Engine
-    end type Engine
-    interface Engine
-     module procedure new_Engine
-    end interface
+In most cases, the default distribution-compatible ``MersenneEngine4`` should
+be used, since the distributions described below require it.
 
 Distributions
 =============
 
 Distributions produce numerical values from the random bitstream provided by
-the Engine.
+an RNG engine. For efficiency, each distribution subroutine accepts an *array*
+of values that are filled with samples of the distribution.
 
 normal_distribution
 -------------------
 
-Fill the given array with Gaussian-distributed numbers generated using the
-given random number engine about the given mean value with the given standard
-deviation, which defaults to 1.
+Each element of the sampled array is distributed according to a Gaussian
+function with the given mean and standard deviation.
 
 uniform_int_distribution
 ------------------------
 
-Fill the given array with uniformly distributed integers generated using the
-given random number engine, between the two bounds (inclusive on both sides).
+Each element is uniformly sampled between the two provided bounds, inclusive on
+both sides.
 
 uniform_real_distribution
 -------------------------
 
-Fill the given array with uniformly distributed real numbers generated using the
-given random number engine, between the two bounds (inclusive on left side only).
+Each element is a sample of a uniform distribution between the two bounds,
+inclusive on left side only.
+
+discrete_distribution
+---------------------
+
+The discrete distribution is constructed with an array of :math:`N` weights:
+the probability that an index in the range :math:`[1, N]` will be selected.
+::
+
+   integer(C_INT), dimension(4), parameter :: weights = [1, 1, 2, 4]
+   integer(C_INT), dimension(1024) :: sampled
+   call discrete_distribution(weights, Engine(), sampled)
+
+In the above example, ``1`` and ``2`` will be present in the ``sampled`` array
+about the same number of times since those indices have equal weight; ``3``
+will be present about twice as often as ``1`` and ``4`` will be present about
+four times as often.
 
+.. note:: The C++ distribution returns values in :math:`[0, N)`, so in
+   accordance with Flibcpp's :ref:`indexing convention <conventions_indexing>`
+   the result is transformed when provided to Fortran users.
 
 .. ############################################################################
 .. end of doc/modules/random.rst

example/sort.f90

@@ -8,7 +8,7 @@ program sort_example
   use, intrinsic :: ISO_C_BINDING
   use flc
   use flc_algorithm, only : sort, shuffle
-  use flc_random, only : Engine, normal_distribution
+  use flc_random, only : Engine => MersenneEngine4, normal_distribution
   use example_utils, only : write_version, read_positive_int, STDOUT
   implicit none
   integer :: arr_size

include/flc_algorithm.i

@@ -212,7 +212,7 @@ static bool includes_impl(const T *data1, size_t size1,
 
 %inline {
 template<class T>
-static void shuffle(std::SWIG_MERSENNE_TWISTER& g, T *DATA, size_t DATASIZE) {
+static void shuffle(std::FLC_DEFAULT_ENGINE& g, T *DATA, size_t DATASIZE) {
     std::shuffle(DATA, DATA + DATASIZE, g);
 }
 }

include/flc_random.i

@@ -9,79 +9,102 @@
 %include "import_flc.i"
 %flc_add_header
 
-/* -------------------------------------------------------------------------
- * Generator class definition
- * ------------------------------------------------------------------------- */
-
 %{
 #include <random>
 %}
 
-// TODO: define a CMake-configurable option for selecting the 32-bit twister
-#if 0
-#define SWIG_MERSENNE_TWISTER mt19937
-#define SWIG_MERSENNE_RESULT_TYPE int32_t
-#else
-#define SWIG_MERSENNE_TWISTER mt19937_64
-#define SWIG_MERSENNE_RESULT_TYPE int64_t
-#endif
-
-%rename(Engine) std::SWIG_MERSENNE_TWISTER;
-%fortran_autofree_rvalue(std::SWIG_MERSENNE_TWISTER);
+/* -------------------------------------------------------------------------
+ * Macros
+ * ------------------------------------------------------------------------- */
 
+%define %flc_random_engine(NAME, GENERATOR, RESULT_TYPE)
+%fortran_autofree_rvalue(std::GENERATOR);
 namespace std {
-class SWIG_MERSENNE_TWISTER
+
+%rename(NAME) GENERATOR;
+%rename("next") GENERATOR::operator();
+
+class GENERATOR
 {
   public:
-    typedef SWIG_MERSENNE_RESULT_TYPE result_type;
+    typedef RESULT_TYPE result_type;
 
-    SWIG_MERSENNE_TWISTER();
-    explicit SWIG_MERSENNE_TWISTER(result_type seed_value);
+    GENERATOR();
+    explicit GENERATOR(result_type seed_value);
     void seed(result_type seed_value);
     void discard(unsigned long long count);
+    result_type operator()();
 };
+
 } // namespace std
+%enddef
 
 /* -------------------------------------------------------------------------
  * RNG distribution routines
- *
- * The generated subroutines will be called from Fortran like:
- *
- *     call uniform_real_distribution(gen, -10, 10, fill_array)
  * ------------------------------------------------------------------------- */
 
-%define %flc_random_distribution1(DISTNAME, TYPE, ARG1)
-%inline {
-static void DISTNAME(TYPE ARG1,
-                     std::SWIG_MERSENNE_TWISTER& g,
-                     TYPE *DATA, size_t DATASIZE) {
-    std::DISTNAME<TYPE> dist(ARG1);
-    TYPE *end = DATA + DATASIZE;
-    while (DATA != end) {
-        *DATA++ = dist(g);
-    }
+%{
+template<class D, class G, class T>
+static inline void flc_generate(D dist, G& g, T* data, size_t size) {
+  T* const end = data + size;
+  while (data != end) {
+    *data++ = dist(g);
+  }
 }
+%}
+
+%apply (const SWIGTYPE *DATA, size_t SIZE) {
+       (const int32_t *WEIGHTS, size_t WEIGHTSIZE),
+       (const int64_t *WEIGHTS, size_t WEIGHTSIZE) };
+
+%inline %{
+template<class T, class G>
+static void uniform_int_distribution(T left, T right,
+                                     G& engine, T* DATA, size_t DATASIZE) {
+  flc_generate(std::uniform_int_distribution<T>(left, right),
+               engine, DATA, DATASIZE);
 }
-%enddef
-%define %flc_random_distribution2(DISTNAME, TYPE, ARG1, ARG2)
-%inline {
-static void DISTNAME(TYPE ARG1, TYPE ARG2,
-                     std::SWIG_MERSENNE_TWISTER& g,
-                     TYPE *DATA, size_t DATASIZE) {
-    std::DISTNAME<TYPE> dist(ARG1, ARG2);
-    TYPE *end = DATA + DATASIZE;
-    while (DATA != end) {
-        *DATA++ = dist(g);
-    }
+
+template<class T, class G>
+static void uniform_real_distribution(T left, T right,
+                                      G& engine, T* DATA, size_t DATASIZE) {
+  flc_generate(std::uniform_real_distribution<T>(left, right),
+               engine, DATA, DATASIZE);
+}
+
+template<class T, class G>
+static void normal_distribution(T mean, T stddev,
+                                G& engine, T* DATA, size_t DATASIZE) {
+  flc_generate(std::normal_distribution<T>(mean, stddev),
+               engine, DATA, DATASIZE);
 }
+
+template<class T, class G>
+static void discrete_distribution(const T* WEIGHTS, size_t WEIGHTSIZE,
+                                  G& engine, T* DATA, size_t DATASIZE) {
+  std::discrete_distribution<T> dist(WEIGHTS, WEIGHTS + WEIGHTSIZE);
+  T* const end = DATA + DATASIZE;
+  while (DATA != end) {
+    *DATA++ = dist(engine) + 1; // Note: transform to Fortran 1-offset
+  }
 }
+%}
+
+%define %flc_distribution(NAME, STDENGINE, TYPE)
+%template(NAME##_distribution) NAME##_distribution< TYPE, std::STDENGINE >;
 %enddef
 
-// Uniform distributions
-%flc_random_distribution2(uniform_int_distribution, int32_t, left, right)
-%flc_random_distribution2(uniform_int_distribution, int64_t, left, right)
-%flc_random_distribution2(uniform_real_distribution, double, left, right)
+// Engines
+%flc_random_engine(MersenneEngine4, mt19937,    int32_t)
+%flc_random_engine(MersenneEngine8, mt19937_64, int64_t)
+
+#define FLC_DEFAULT_ENGINE mt19937
+%flc_distribution(uniform_int,  FLC_DEFAULT_ENGINE, int32_t)
+%flc_distribution(uniform_int,  FLC_DEFAULT_ENGINE, int64_t)
+%flc_distribution(uniform_real, FLC_DEFAULT_ENGINE, double)
+
+%flc_distribution(normal, FLC_DEFAULT_ENGINE, double)
 
-// Gaussian distribution
-%flc_random_distribution1(normal_distribution, double, mean)
-%flc_random_distribution2(normal_distribution, double, mean, stddev)
+// Discrete sampling distribution
+%flc_distribution(discrete, FLC_DEFAULT_ENGINE, int32_t)
+%flc_distribution(discrete, FLC_DEFAULT_ENGINE, int64_t)

include/flc_vector.i

@@ -15,7 +15,7 @@
  * Macro definitions
  * ------------------------------------------------------------------------- */
 
-%define EXTEND_STD_VECTOR_POD(CTYPE)
+%define %flc_extend_vector_pod(CTYPE)
   %apply (const SWIGTYPE *DATA, ::size_t SIZE)
     { (const CTYPE* DATA, size_type SIZE) };
 
@@ -50,6 +50,9 @@
  * These provide an efficient constructor from a Fortan array view. It also
  * offers a "view" functionality for getting an array pointer to the
  * vector-owned data.
+ *
+ * This definition is considered part of the \em public API so that downstream
+ * apps that generate FLC-based bindings can instantiate their own POD vectors.
  */
 %define %specialize_std_vector_pod(T)
 
@@ -62,7 +65,7 @@ namespace std {
     SWIG_STD_VECTOR_COMMON(T, const T&)
     SWIG_STD_VECTOR_REF(T)
     %extend {
-      EXTEND_STD_VECTOR_POD(T)
+      %flc_extend_vector_pod(T)
     }
   };
 }

src/flc_algorithm.f90

@@ -1508,7 +1508,7 @@ function swigf_includes__SWIG_7(data1, data2, cmp) &
 
 subroutine swigf_shuffle__SWIG_1(g, data)
 use, intrinsic :: ISO_C_BINDING
-class(Engine), intent(in) :: g
+class(MersenneEngine4), intent(in) :: g
 integer(C_INT32_T), dimension(:), target :: data
 type(SwigClassWrapper) :: farg1 
 type(SwigArrayWrapper) :: farg2 
@@ -1520,7 +1520,7 @@ subroutine swigf_shuffle__SWIG_1(g, data)
 
 subroutine swigf_shuffle__SWIG_2(g, data)
 use, intrinsic :: ISO_C_BINDING
-class(Engine), intent(in) :: g
+class(MersenneEngine4), intent(in) :: g
 integer(C_INT64_T), dimension(:), target :: data
 type(SwigClassWrapper) :: farg1 
 type(SwigArrayWrapper) :: farg2 
@@ -1532,7 +1532,7 @@ subroutine swigf_shuffle__SWIG_2(g, data)
 
 subroutine swigf_shuffle__SWIG_3(g, data)
 use, intrinsic :: ISO_C_BINDING
-class(Engine), intent(in) :: g
+class(MersenneEngine4), intent(in) :: g
 real(C_DOUBLE), dimension(:), target :: data
 type(SwigClassWrapper) :: farg1 
 type(SwigArrayWrapper) :: farg2 
@@ -1544,7 +1544,7 @@ subroutine swigf_shuffle__SWIG_3(g, data)
 
 subroutine swigf_shuffle__SWIG_4(g, data)
 use, intrinsic :: ISO_C_BINDING
-class(Engine), intent(in) :: g
+class(MersenneEngine4), intent(in) :: g
 type(C_PTR), dimension(:), target :: data
 type(SwigClassWrapper) :: farg1 
 type(SwigArrayWrapper) :: farg2