From 3bcbb1615b907c8b81ff5fcf274ff2ac83137316 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 12 Jun 2018 08:16:06 +0200
Subject: [PATCH] ENH: improve uniformity of Random::position<label>()
(closes #865)
- start/end values were underrepresented due to rounding.
Now extend the range to include -0.5 and +0.5 beyond the usual
range to ensure the same number density.
---
applications/test/Random/Test-Random.C | 16 ++++++++++-
.../primitives/random/Random/Random.C | 28 ++++++++++---------
2 files changed, 30 insertions(+), 14 deletions(-)
diff --git a/applications/test/Random/Test-Random.C b/applications/test/Random/Test-Random.C
index b4bd2f25522..a724fe83577 100644
--- a/applications/test/Random/Test-Random.C
+++ b/applications/test/Random/Test-Random.C
@@ -174,7 +174,21 @@ int main(int argc, char *argv[])
}
}
- cout<< nl << "Done." << endl;
+ // Test uniformity of random
+ {
+ List<label> samples(20, Zero);
+
+ Random rnd(123456);
+ for (label i=0; i < 1000*samples.size(); ++i)
+ {
+ ++samples[rnd.position<label>(0,19)];
+ }
+
+ Info<< nl << "uniform [0,20)" << nl << " "
+ << flatOutput(samples) << nl;
+ }
+
+ Info<< nl << "Done." << endl;
return 0;
}
diff --git a/src/OpenFOAM/primitives/random/Random/Random.C b/src/OpenFOAM/primitives/random/Random/Random.C
index 307fcb0eb1b..ed09267b126 100644
--- a/src/OpenFOAM/primitives/random/Random/Random.C
+++ b/src/OpenFOAM/primitives/random/Random/Random.C
@@ -121,14 +121,16 @@ Foam::scalar Foam::Random::position
template<>
Foam::label Foam::Random::position(const label& start, const label& end)
{
- return start + round(scalar01()*(end - start));
+ // Extend range from [0, N-1] to (-0.5, N-0.5) to ensure that round()
+ // results in the same number density at the ends.
+ return start + round(scalar01()*((end - start) + 0.998) - 0.499);
}
template<>
Foam::scalar Foam::Random::globalSample01()
{
- scalar value = -GREAT;
+ scalar value(-GREAT);
if (Pstream::master())
{
@@ -144,7 +146,7 @@ Foam::scalar Foam::Random::globalSample01()
template<>
Foam::label Foam::Random::globalSample01()
{
- label value = labelMin;
+ label value(labelMin);
if (Pstream::master())
{
@@ -160,7 +162,7 @@ Foam::label Foam::Random::globalSample01()
template<>
Foam::scalar Foam::Random::globalGaussNormal()
{
- scalar value = -GREAT;
+ scalar value(-GREAT);
if (Pstream::master())
{
@@ -176,16 +178,16 @@ Foam::scalar Foam::Random::globalGaussNormal()
template<>
Foam::label Foam::Random::globalGaussNormal()
{
- scalar value = -GREAT;
+ label value(labelMin);
if (Pstream::master())
{
- value = GaussNormal<scalar>();
+ value = GaussNormal<label>();
}
Pstream::scatter(value);
- return round(value);
+ return value;
}
@@ -196,16 +198,16 @@ Foam::scalar Foam::Random::globalPosition
const scalar& end
)
{
- scalar value = -GREAT;
+ scalar value(-GREAT);
if (Pstream::master())
{
- value = scalar01()*(end - start);
+ value = position<scalar>(start, end);
}
Pstream::scatter(value);
- return start + value;
+ return value;
}
@@ -216,16 +218,16 @@ Foam::label Foam::Random::globalPosition
const label& end
)
{
- label value = labelMin;
+ label value(labelMin);
if (Pstream::master())
{
- value = round(scalar01()*(end - start));
+ value = position<label>(start, end);
}
Pstream::scatter(value);
- return start + value;
+ return value;
}
--
GitLab