From 42a9a6ae5a29e2a75d4466fc690b8accf7610d4b Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Thu, 21 Nov 2019 09:40:00 +0100
Subject: [PATCH] ENH: add generator class for uniform/gaussian random numbers

- can be used in combination with std::generate, or as a substitute
  unary operator to supply random numbers for std::transform.
---
 applications/test/Random/Test-Random.C        |  54 ++++++-
 .../primitives/random/Random/Random.H         | 134 +++++++++++++++++-
 2 files changed, 185 insertions(+), 3 deletions(-)

diff --git a/applications/test/Random/Test-Random.C b/applications/test/Random/Test-Random.C
index 8664bf8266..75b286755c 100644
--- a/applications/test/Random/Test-Random.C
+++ b/applications/test/Random/Test-Random.C
@@ -31,7 +31,6 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "Rand48.H"
 #include "Random.H"
 
 #include <cstdlib>
@@ -218,6 +217,59 @@ int main(int argc, char *argv[])
     testPosition(20);
     testPosition(3);
 
+
+    // Generators
+    {
+        const label n = 20;
+        List<label> samples(n, Zero);
+
+        Random::uniformGeneratorOp<label> gen(0, n-1);
+
+        // Test uniformity of random
+        samples = Zero;
+        for (label i=0; i < 100000*n; ++i)
+        {
+            // Calling with/without parameter is the same
+            if (i % 3)
+            {
+                ++samples[gen()];
+            }
+            else
+            {
+                // Unary ignores any parameter
+                ++samples[gen(3.14159)];
+            }
+        }
+
+        Info<< nl << "Uniform generator [0," << n << ")\n  "
+            << flatOutput(samples) << nl;
+    }
+
+    {
+        Random::gaussianGeneratorOp<scalar> gen;
+
+        Info<< "Some gaussian generated values" << nl;
+
+        for (label i=0; i < 20; ++i)
+        {
+            Info<< ' ';
+
+            // Calling with/without parameter is the same
+            if (i % 3)
+            {
+                Info<< gen();
+            }
+            else
+            {
+                // Unary ignores any parameter
+                Info<< gen(3.14159);
+            }
+        }
+
+        Info<< nl;
+    }
+
+
     // This should fail (in FULLDEBUG)
     const bool throwingError = FatalError.throwExceptions();
     try
diff --git a/src/OpenFOAM/primitives/random/Random/Random.H b/src/OpenFOAM/primitives/random/Random/Random.H
index 81d68c3171..c09fd3ac58 100644
--- a/src/OpenFOAM/primitives/random/Random/Random.H
+++ b/src/OpenFOAM/primitives/random/Random/Random.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2017-2018 OpenCFD Ltd.
+    Copyright (C) 2017-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,7 +59,7 @@ template<class T> class UList;
 
 class Random
 {
-    // Private data
+    // Private Data
 
         //- Initial random number seed
         label seed_;
@@ -85,6 +85,11 @@ class Random
 
 public:
 
+        // Forward declarations - generator classes
+        template<class T> class uniformGeneratorOp;
+        template<class T> class gaussianGeneratorOp;
+
+
     // Public Static Data
 
         //- The default seed value (name may change in the future)
@@ -203,6 +208,131 @@ template<>
 label Random::globalPosition<label>(const label& start, const label& end);
 
 
+/*---------------------------------------------------------------------------*\
+                 Class Random::uniformGeneratorOp Declaration
+\*---------------------------------------------------------------------------*/
+
+//- A generator class returning a uniformly distributed random number
+//- on the given interval.
+//
+//  \sa std::generate()
+template<class T>
+class Random::uniformGeneratorOp
+:
+    public Random
+{
+    // Private data
+
+        //- The interval
+        const T min_;
+        const T max_;
+
+
+    // Private Member Functions
+
+        //- Generate a random number. Treat as mutable.
+        T generate() const
+        {
+            return const_cast<uniformGeneratorOp&>(*this)
+                .position<T>(min_, max_);
+        }
+
+
+public:
+
+    // Constructors
+
+        //- Construct null or with seed value. Uses default [0,1] interval
+        explicit uniformGeneratorOp(const label seed = Random::defaultSeed)
+        :
+            uniformGeneratorOp(seed, pTraits<T>::zero, pTraits<T>::one)
+        {}
+
+        //- Construct with specified interval, using default seed
+        uniformGeneratorOp(const T& minval, const T& maxval)
+        :
+            uniformGeneratorOp(Random::defaultSeed, minval, maxval)
+        {}
+
+        //- Construct with seed value and specified interval
+        uniformGeneratorOp(const label seed, const T& minval, const T& maxval)
+        :
+            Random(seed),
+            min_(minval),
+            max_(maxval)
+        {}
+
+
+    // Member Operators
+
+        //- Generate a random number
+        T operator()()
+        {
+            return generate();
+        }
+
+        //- Ignore parameter and generate a random number, which allows it
+        //- to be used as a replacement for general unary operators.
+        template<class U>
+        T operator()(const U&) const
+        {
+            return generate();
+        }
+};
+
+
+/*---------------------------------------------------------------------------*\
+                 Class Random::gaussianGeneratorOp Declaration
+\*---------------------------------------------------------------------------*/
+
+//- A generator class returning a gaussian distributed random number.
+//
+//  \sa std::generate()
+template<class T>
+class Random::gaussianGeneratorOp
+:
+    public Random
+{
+
+    // Private Member Functions
+
+        //- Generate a random number. Treat as mutable.
+        T generate() const
+        {
+            return const_cast<gaussianGeneratorOp&>(*this)
+                .GaussNormal<T>();
+        }
+
+
+public:
+
+    // Constructors
+
+        //- Construct null or with seed value. Uses default [0,1] interval
+        explicit gaussianGeneratorOp(const label seed = Random::defaultSeed)
+        :
+            Random(seed)
+        {}
+
+
+    // Member Operators
+
+        //- Generate a random number
+        T operator()()
+        {
+            return generate();
+        }
+
+        //- Ignore parameter and generate a random number, which allows it
+        //- to be used as a replacement for general unary operators.
+        template<class U>
+        T operator()(const U&) const
+        {
+            return generate();
+        }
+};
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam