Skip to content

Commit 613a815

Browse files
sethrjsethrj
authored
Merge pull request #8 from swig-fortran/random
Extend and improve flc_random
2 parents f425106 + aca714a commit 613a815

13 files changed

+729
-299
lines changed

‎CMakeLists.txt

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66
#---------------------------------------------------------------------------#
77

88
cmake_minimum_required(VERSION 3.8)
9-
project(Flibcpp VERSION 0.2.1 LANGUAGES CXX Fortran)
9+
project(Flibcpp VERSION 0.3.0 LANGUAGES CXX Fortran)
1010

1111
#---------------------------------------------------------------------------#
1212
# OPTIONS

‎doc/modules/algorithm.rst

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -206,12 +206,12 @@ shuffle
206206
-------
207207

208208
The "shuffle" subroutine depends on the :ref:`modules_random` module so that it
209-
can use the supported random number generator to randomly reorder an array.
209+
can use the default random number generator to randomly reorder an array.
210210

211211
Example::
212212

213213
use flc_algorithm, only : shuffle
214-
use flc_random, only : Engine
214+
use flc_random, only : Engine => MersenneEngine4
215215
implicit none
216216
integer :: i
217217
integer, dimension(8) :: iarr = (/ ((i), i = -4, 3) /)

‎doc/modules/random.rst

Lines changed: 58 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -9,73 +9,106 @@ Random
99
******
1010

1111
The ``flc_random`` module contains advanced pseudo-random number generation
12-
from the `<random>`_ C++ header. Currently
13-
random number generation is hardwired to the ``std::mt19937_64`` 64-bit
14-
Mersenne Twister by Matsumoto and Nishimura.
12+
from the `<random>`_ C++ header.
1513

1614
The C++11 way of generating random numbers takes some getting used to. Rather
1715
than having a single global function that returns a random real number in the
1816
range :math:`[0,1)`, C++11 has independent *engine* objects that generate
1917
streams of random bits. Different *distribution* objects convert those bits
2018
into samples of a distribution.
2119

20+
Although C++11 defines a dizzying variety of random number engines,
21+
``flc_random`` wraps just two: the 32- and 64-bit `Mersenne Twister`
22+
algorithms. The 32-bit version (``MersenneEngine4``) is currently the only
23+
engine type that can be used with distributions and algorithms.
24+
2225
Flibcpp wraps distribution objects as independent subroutines. Each subroutine
2326
accepts the constructor parameters of the distribution, the random engine, and
2427
a target Fortran array to be filled with random numbers.
2528

2629
Generating an array with 10 normally-distributed reals with a mean of 8 and a
2730
standard deviation of 2 is done as such::
2831

29-
use flc_random, only : Engine, normal_distribution
32+
use flc_random, only : Engine => MersenneEngine4, normal_distribution
3033
real(C_DOUBLE), dimension(20) :: arr
3134
type(Engine) :: rng
3235

3336
rng = Engine()
3437
call normal_distribution(8.0d0, 2.0d0, rng, arr)
38+
call rng%release()
3539

3640
.. _<random> : https://en.cppreference.com/w/cpp/numeric/random
41+
.. _Mersenne Twister : https://en.wikipedia.org/wiki/Mersenne_Twister
42+
43+
Engines
44+
=======
45+
46+
The two Mersenne twister engines in ``flc_random`` return different-sized
47+
integers per call:
48+
49+
- ``MersenneEngine4``: each invocation returns a 32-bit ``integer(4)``
50+
- ``MersenneEngine8``: each invocation returns a 64-bit ``integer(8)``
51+
52+
Engines can be constructed using one of two interface functions: the
53+
argument-free ``MersenneEngine4()`` uses the default seed, and the engine takes
54+
a single argument ``MersenneEngine4(1234567)`` with the seed. Alternatively,
55+
the seed can be set (or reset) using the ``seed()`` type-bound procedure.
3756

38-
Engine
39-
======
57+
Generally, engines are used with distributions (described below). However, if
58+
necessary, individual randomly generated values can be obtained by calling
59+
the ``next()`` type-bound procedure.
4060

41-
The ``Engine`` is a derived type that wraps the Mersenne Twister PRNG engine.
42-
Its interface is::
61+
.. warning:: C++ generates *unsigned* integers with entropy in every bit. This
62+
means that the integers obtained from ``engine%next()``, reinterpreted as
63+
signed Fortran integers, may be negative.
4364

44-
type, public :: Engine
45-
contains
46-
procedure :: seed => swigf_Engine_seed
47-
procedure :: discard => swigf_Engine_discard
48-
procedure :: release => delete_Engine
49-
end type Engine
50-
interface Engine
51-
module procedure new_Engine
52-
end interface
65+
In most cases, the default distribution-compatible ``MersenneEngine4`` should
66+
be used, since the distributions described below require it.
5367

5468
Distributions
5569
=============
5670

5771
Distributions produce numerical values from the random bitstream provided by
58-
the Engine.
72+
an RNG engine. For efficiency, each distribution subroutine accepts an *array*
73+
of values that are filled with samples of the distribution.
5974

6075
normal_distribution
6176
-------------------
6277

63-
Fill the given array with Gaussian-distributed numbers generated using the
64-
given random number engine about the given mean value with the given standard
65-
deviation, which defaults to 1.
78+
Each element of the sampled array is distributed according to a Gaussian
79+
function with the given mean and standard deviation.
6680

6781
uniform_int_distribution
6882
------------------------
6983

70-
Fill the given array with uniformly distributed integers generated using the
71-
given random number engine, between the two bounds (inclusive on both sides).
84+
Each element is uniformly sampled between the two provided bounds, inclusive on
85+
both sides.
7286

7387
uniform_real_distribution
7488
-------------------------
7589

76-
Fill the given array with uniformly distributed real numbers generated using the
77-
given random number engine, between the two bounds (inclusive on left side only).
90+
Each element is a sample of a uniform distribution between the two bounds,
91+
inclusive on left side only.
92+
93+
discrete_distribution
94+
---------------------
95+
96+
The discrete distribution is constructed with an array of :math:`N` weights:
97+
the probability that an index in the range :math:`[1, N]` will be selected.
98+
::
99+
100+
integer(C_INT), dimension(4), parameter :: weights = [1, 1, 2, 4]
101+
integer(C_INT), dimension(1024) :: sampled
102+
call discrete_distribution(weights, Engine(), sampled)
103+
104+
In the above example, ``1`` and ``2`` will be present in the ``sampled`` array
105+
about the same number of times since those indices have equal weight; ``3``
106+
will be present about twice as often as ``1`` and ``4`` will be present about
107+
four times as often.
78108

109+
.. note:: The C++ distribution returns values in :math:`[0, N)`, so in
110+
accordance with Flibcpp's :ref:`indexing convention <conventions_indexing>`
111+
the result is transformed when provided to Fortran users.
79112

80113
.. ############################################################################
81114
.. end of doc/modules/random.rst

‎example/sort.f90

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ program sort_example
88
use, intrinsic :: ISO_C_BINDING
99
use flc
1010
use flc_algorithm, only : sort, shuffle
11-
use flc_random, only : Engine, normal_distribution
11+
use flc_random, only : Engine => MersenneEngine4, normal_distribution
1212
use example_utils, only : write_version, read_positive_int, STDOUT
1313
implicit none
1414
integer :: arr_size

‎include/flc_algorithm.i

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -212,7 +212,7 @@ static bool includes_impl(const T *data1, size_t size1,
212212

213213
%inline {
214214
template<class T>
215-
static void shuffle(std::SWIG_MERSENNE_TWISTER& g, T *DATA, size_t DATASIZE) {
215+
static void shuffle(std::FLC_DEFAULT_ENGINE& g, T *DATA, size_t DATASIZE) {
216216
std::shuffle(DATA, DATA + DATASIZE, g);
217217
}
218218
}

‎include/flc_random.i

Lines changed: 74 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -9,79 +9,102 @@
99
%include "import_flc.i"
1010
%flc_add_header
1111

12-
/* -------------------------------------------------------------------------
13-
* Generator class definition
14-
* ------------------------------------------------------------------------- */
15-
1612
%{
1713
#include <random>
1814
%}
1915

20-
// TODO: define a CMake-configurable option for selecting the 32-bit twister
21-
#if 0
22-
#define SWIG_MERSENNE_TWISTER mt19937
23-
#define SWIG_MERSENNE_RESULT_TYPE int32_t
24-
#else
25-
#define SWIG_MERSENNE_TWISTER mt19937_64
26-
#define SWIG_MERSENNE_RESULT_TYPE int64_t
27-
#endif
28-
29-
%rename(Engine) std::SWIG_MERSENNE_TWISTER;
30-
%fortran_autofree_rvalue(std::SWIG_MERSENNE_TWISTER);
16+
/* -------------------------------------------------------------------------
17+
* Macros
18+
* ------------------------------------------------------------------------- */
3119

20+
%define %flc_random_engine(NAME, GENERATOR, RESULT_TYPE)
21+
%fortran_autofree_rvalue(std::GENERATOR);
3222
namespace std {
33-
class SWIG_MERSENNE_TWISTER
23+
24+
%rename(NAME) GENERATOR;
25+
%rename("next") GENERATOR::operator();
26+
27+
class GENERATOR
3428
{
3529
public:
36-
typedef SWIG_MERSENNE_RESULT_TYPE result_type;
30+
typedef RESULT_TYPE result_type;
3731

38-
SWIG_MERSENNE_TWISTER();
39-
explicit SWIG_MERSENNE_TWISTER(result_type seed_value);
32+
GENERATOR();
33+
explicit GENERATOR(result_type seed_value);
4034
void seed(result_type seed_value);
4135
void discard(unsigned long long count);
36+
result_type operator()();
4237
};
38+
4339
} // namespace std
40+
%enddef
4441

4542
/* -------------------------------------------------------------------------
4643
* RNG distribution routines
47-
*
48-
* The generated subroutines will be called from Fortran like:
49-
*
50-
* call uniform_real_distribution(gen, -10, 10, fill_array)
5144
* ------------------------------------------------------------------------- */
5245

53-
%define %flc_random_distribution1(DISTNAME, TYPE, ARG1)
54-
%inline {
55-
static void DISTNAME(TYPE ARG1,
56-
std::SWIG_MERSENNE_TWISTER& g,
57-
TYPE *DATA, size_t DATASIZE) {
58-
std::DISTNAME<TYPE> dist(ARG1);
59-
TYPE *end = DATA + DATASIZE;
60-
while (DATA != end) {
61-
*DATA++ = dist(g);
62-
}
46+
%{
47+
template<class D, class G, class T>
48+
static inline void flc_generate(D dist, G& g, T* data, size_t size) {
49+
T* const end = data + size;
50+
while (data != end) {
51+
*data++ = dist(g);
52+
}
6353
}
54+
%}
55+
56+
%apply (const SWIGTYPE *DATA, size_t SIZE) {
57+
(const int32_t *WEIGHTS, size_t WEIGHTSIZE),
58+
(const int64_t *WEIGHTS, size_t WEIGHTSIZE) };
59+
60+
%inline %{
61+
template<class T, class G>
62+
static void uniform_int_distribution(T left, T right,
63+
G& engine, T* DATA, size_t DATASIZE) {
64+
flc_generate(std::uniform_int_distribution<T>(left, right),
65+
engine, DATA, DATASIZE);
6466
}
65-
%enddef
66-
%define %flc_random_distribution2(DISTNAME, TYPE, ARG1, ARG2)
67-
%inline {
68-
static void DISTNAME(TYPE ARG1, TYPE ARG2,
69-
std::SWIG_MERSENNE_TWISTER& g,
70-
TYPE *DATA, size_t DATASIZE) {
71-
std::DISTNAME<TYPE> dist(ARG1, ARG2);
72-
TYPE *end = DATA + DATASIZE;
73-
while (DATA != end) {
74-
*DATA++ = dist(g);
75-
}
67+
68+
template<class T, class G>
69+
static void uniform_real_distribution(T left, T right,
70+
G& engine, T* DATA, size_t DATASIZE) {
71+
flc_generate(std::uniform_real_distribution<T>(left, right),
72+
engine, DATA, DATASIZE);
73+
}
74+
75+
template<class T, class G>
76+
static void normal_distribution(T mean, T stddev,
77+
G& engine, T* DATA, size_t DATASIZE) {
78+
flc_generate(std::normal_distribution<T>(mean, stddev),
79+
engine, DATA, DATASIZE);
7680
}
81+
82+
template<class T, class G>
83+
static void discrete_distribution(const T* WEIGHTS, size_t WEIGHTSIZE,
84+
G& engine, T* DATA, size_t DATASIZE) {
85+
std::discrete_distribution<T> dist(WEIGHTS, WEIGHTS + WEIGHTSIZE);
86+
T* const end = DATA + DATASIZE;
87+
while (DATA != end) {
88+
*DATA++ = dist(engine) + 1; // Note: transform to Fortran 1-offset
89+
}
7790
}
91+
%}
92+
93+
%define %flc_distribution(NAME, STDENGINE, TYPE)
94+
%template(NAME##_distribution) NAME##_distribution< TYPE, std::STDENGINE >;
7895
%enddef
7996

80-
// Uniform distributions
81-
%flc_random_distribution2(uniform_int_distribution, int32_t, left, right)
82-
%flc_random_distribution2(uniform_int_distribution, int64_t, left, right)
83-
%flc_random_distribution2(uniform_real_distribution, double, left, right)
97+
// Engines
98+
%flc_random_engine(MersenneEngine4, mt19937, int32_t)
99+
%flc_random_engine(MersenneEngine8, mt19937_64, int64_t)
100+
101+
#define FLC_DEFAULT_ENGINE mt19937
102+
%flc_distribution(uniform_int, FLC_DEFAULT_ENGINE, int32_t)
103+
%flc_distribution(uniform_int, FLC_DEFAULT_ENGINE, int64_t)
104+
%flc_distribution(uniform_real, FLC_DEFAULT_ENGINE, double)
105+
106+
%flc_distribution(normal, FLC_DEFAULT_ENGINE, double)
84107

85-
// Gaussian distribution
86-
%flc_random_distribution1(normal_distribution, double, mean)
87-
%flc_random_distribution2(normal_distribution, double, mean, stddev)
108+
// Discrete sampling distribution
109+
%flc_distribution(discrete, FLC_DEFAULT_ENGINE, int32_t)
110+
%flc_distribution(discrete, FLC_DEFAULT_ENGINE, int64_t)

‎include/flc_vector.i

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@
1515
* Macro definitions
1616
* ------------------------------------------------------------------------- */
1717

18-
%define EXTEND_STD_VECTOR_POD(CTYPE)
18+
%define %flc_extend_vector_pod(CTYPE)
1919
%apply (const SWIGTYPE *DATA, ::size_t SIZE)
2020
{ (const CTYPE* DATA, size_type SIZE) };
2121

@@ -50,6 +50,9 @@
5050
* These provide an efficient constructor from a Fortan array view. It also
5151
* offers a "view" functionality for getting an array pointer to the
5252
* vector-owned data.
53+
*
54+
* This definition is considered part of the \em public API so that downstream
55+
* apps that generate FLC-based bindings can instantiate their own POD vectors.
5356
*/
5457
%define %specialize_std_vector_pod(T)
5558

@@ -62,7 +65,7 @@ namespace std {
6265
SWIG_STD_VECTOR_COMMON(T, const T&)
6366
SWIG_STD_VECTOR_REF(T)
6467
%extend {
65-
EXTEND_STD_VECTOR_POD(T)
68+
%flc_extend_vector_pod(T)
6669
}
6770
};
6871
}

‎src/flc_algorithm.f90

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1508,7 +1508,7 @@ function swigf_includes__SWIG_7(data1, data2, cmp) &
15081508

15091509
subroutine swigf_shuffle__SWIG_1(g, data)
15101510
use, intrinsic :: ISO_C_BINDING
1511-
class(Engine), intent(in) :: g
1511+
class(MersenneEngine4), intent(in) :: g
15121512
integer(C_INT32_T), dimension(:), target :: data
15131513
type(SwigClassWrapper) :: farg1
15141514
type(SwigArrayWrapper) :: farg2
@@ -1520,7 +1520,7 @@ subroutine swigf_shuffle__SWIG_1(g, data)
15201520

15211521
subroutine swigf_shuffle__SWIG_2(g, data)
15221522
use, intrinsic :: ISO_C_BINDING
1523-
class(Engine), intent(in) :: g
1523+
class(MersenneEngine4), intent(in) :: g
15241524
integer(C_INT64_T), dimension(:), target :: data
15251525
type(SwigClassWrapper) :: farg1
15261526
type(SwigArrayWrapper) :: farg2
@@ -1532,7 +1532,7 @@ subroutine swigf_shuffle__SWIG_2(g, data)
15321532

15331533
subroutine swigf_shuffle__SWIG_3(g, data)
15341534
use, intrinsic :: ISO_C_BINDING
1535-
class(Engine), intent(in) :: g
1535+
class(MersenneEngine4), intent(in) :: g
15361536
real(C_DOUBLE), dimension(:), target :: data
15371537
type(SwigClassWrapper) :: farg1
15381538
type(SwigArrayWrapper) :: farg2
@@ -1544,7 +1544,7 @@ subroutine swigf_shuffle__SWIG_3(g, data)
15441544

15451545
subroutine swigf_shuffle__SWIG_4(g, data)
15461546
use, intrinsic :: ISO_C_BINDING
1547-
class(Engine), intent(in) :: g
1547+
class(MersenneEngine4), intent(in) :: g
15481548
type(C_PTR), dimension(:), target :: data
15491549
type(SwigClassWrapper) :: farg1
15501550
type(SwigArrayWrapper) :: farg2

0 commit comments

Comments
 (0)