From 844685d29d3fc8389fb1b477a178c1909c752ec5 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Mon, 23 Nov 2009 12:22:13 +0000
Subject: [PATCH] Finished working version of templated Distribution class.

Changed the underlying map to be a Map<scalar> to allow the addition
of a value to the Distribution to be weighted (by timestep or by cell
size for example) component by-component - default weight is one.

Adding normalisation, write, mean snd median functions.

Test application DistributionTest draws numbers from the Random class
and creates distributions of them, writing to disk, for the main
types: scalar, vector, symmTensior, sphericalTensor and tensor.

Creating a labelVector distibution, but need to check what is
happening to the mean and median values and what happens when an odd
number is used as a binWidth, as truncation to label will be
occurring.
---
 .../test/Distribution/DistributionTest.C      | 193 ++++++++++--
 .../Map/Distribution/Distribution.C           | 277 +++++++++++++++---
 .../Map/Distribution/Distribution.H           |  53 ++--
 3 files changed, 435 insertions(+), 88 deletions(-)

diff --git a/applications/test/Distribution/DistributionTest.C b/applications/test/Distribution/DistributionTest.C
index cdbbc4e469..687d5746ac 100644
--- a/applications/test/Distribution/DistributionTest.C
+++ b/applications/test/Distribution/DistributionTest.C
@@ -28,12 +28,24 @@ Application
 Description
     Test the Distribution class
 
+    Plot normal distribution test in gnuplot using:
+
+    @verbatim
+    normalDistribution(mean, sigma, x) = \
+        sqrt(1.0/(2.0*pi*sigma**2))*exp(-(x - mean)**2.0/(2.0*sigma**2))
+
+    plot normalDistribution(8.5, 2.5, x)
+    @endverbatim
+
 \*---------------------------------------------------------------------------*/
 
 #include "vector.H"
 #include "labelVector.H"
 #include "tensor.H"
 #include "Distribution.H"
+#include "Random.H"
+#include "dimensionedTypes.H"
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,36 +53,173 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
-    Distribution<scalar> dS(scalar(3.0));
-    Distribution<vector> dV(vector(0.1, 0.24, 0.18));
-    Distribution<labelVector> dLV(labelVector(2,3,4));
+    Random R(918273);
 
-    dS.add(1.0);
-    dS.add(2.0);
-    dS.add(12.0);
-    dS.add(1.3);
+    {
+        // scalar
+        label randomDistributionTestSize = 50000000;
 
-    Info<< dS << nl << dS.raw() << endl;
+        Distribution<scalar> dS(scalar(5e-2));
 
-    vector vA(1.2, 1.3, 1.1);
-    vector vB(1.3, 1.5, 1.6);
-    vector vC(0.5, 5.3, 1.1);
+        Info<< nl << "Distribution<scalar>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from GaussNormal distribution."
+            << endl;
 
-    dV.add(vA);
-    dV.add(vB);
-    dV.add(vC);
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dS.add(2.5*R.GaussNormal() + 8.5);
+        }
 
-    Info<< dV << nl << dV.raw() << endl;
+        Info<< "Mean " << dS.mean() << nl
+            << "Median " << dS.median()
+            << endl;
 
-    labelVector lVA(6, 8, 11);
-    labelVector lVB(-12, -3, 6);
-    labelVector lVC(-4, -2, 5);
+        dS.write("Distribution_scalar_test", dS.normalised());
+    }
 
-    dLV.add(lVA);
-    dLV.add(lVB);
-    dLV.add(lVC);
+    {
+        // vector
+        Distribution<vector> dV(vector(0.1, 0.05, 0.15));
 
-    Info<< dLV << nl << dLV.raw() << endl;
+        label randomDistributionTestSize = 1000000;
+
+        Info<< nl << "Distribution<vector>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from uniform and GaussNormal distribution."
+            << endl;
+
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dV.add(R.vector01());
+
+            // Adding separate GaussNormal components with component
+            // weights
+
+            dV.add
+            (
+                vector
+                (
+                    R.GaussNormal()*3.0 + 1.5,
+                    R.GaussNormal()*0.25 + 4.0,
+                    R.GaussNormal()*3.0 - 1.5
+                ),
+                vector(1.0, 2.0, 5.0)
+            );
+        }
+
+        Info<< "Mean " << dV.mean() << nl
+            << "Median " << dV.median()
+            << endl;
+
+        dV.write("Distribution_vector_test", dV.normalised());
+    }
+
+    {
+        // labelVector
+        Distribution<labelVector> dLV(labelVector::one*10);
+
+        label randomDistributionTestSize = 2000000;
+
+        Info<< nl << "Distribution<labelVector>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from uniform distribution."
+            << endl;
+
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dLV.add
+            (
+                labelVector
+                (
+                    R.integer(-1000, 1000),
+                    R.integer(-5000, 5000),
+                    R.integer(-2000, 7000)
+                )
+            );
+        }
+
+        Info<< "Mean " << dLV.mean() << nl
+            << "Median " << dLV.median()
+            << endl;
+
+        dLV.write("Distribution_labelVector_test", dLV.normalised());
+    }
+
+    {
+        // tensor
+        Distribution<tensor> dT(tensor::one*1e-2);
+
+        label randomDistributionTestSize = 2000000;
+
+        Info<< nl << "Distribution<tensor>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from uniform distribution."
+            << endl;
+
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dT.add(R.tensor01());
+        }
+
+        Info<< "Mean " << dT.mean() << nl
+            << "Median " << dT.median()
+            << endl;
+
+        dT.write("Distribution_tensor_test", dT.normalised());
+    }
+
+    {
+        // symmTensor
+        Distribution<symmTensor> dSyT(symmTensor::one*1e-2);
+
+        label randomDistributionTestSize = 2000000;
+
+        Info<< nl << "Distribution<symmTensor>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from uniform distribution."
+            << endl;
+
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dSyT.add(R.symmTensor01());
+        }
+
+        Info<< "Mean " << dSyT.mean() << nl
+            << "Median " << dSyT.median()
+            << endl;
+
+        dSyT.write("Distribution_symmTensor_test", dSyT.normalised());
+    }
+
+    {
+        // sphericalTensor
+        Distribution<sphericalTensor> dSpT(sphericalTensor::one*1e-2);
+
+        label randomDistributionTestSize = 50000000;
+
+        Info<< nl << "Distribution<sphericalTensor>" << nl
+            << "Sampling "
+            << randomDistributionTestSize
+            << " times from uniform distribution."
+            << endl;
+
+        for (label i = 0; i < randomDistributionTestSize; i++)
+        {
+            dSpT.add(R.sphericalTensor01());
+        }
+
+        Info<< "Mean " << dSpT.mean() << nl
+            << "Median " << dSpT.median()
+            << endl;
+
+        dSpT.write("Distribution_sphericalTensor_test", dSpT.normalised());
+    }
 
     Info<< nl << "End" << nl << endl;
 
diff --git a/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.C b/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.C
index 3a0defe22c..cee8014f09 100644
--- a/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.C
+++ b/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.C
@@ -25,16 +25,14 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "Distribution.H"
-
-// * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * * //
-
+#include "OFstream.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class Type>
 Foam::Distribution<Type>::Distribution()
 :
-    List< Map<label> >(pTraits<Type>::nComponents),
+    List< Map<scalar> >(pTraits<Type>::nComponents),
     binWidth_(pTraits<Type>::one)
 {}
 
@@ -42,7 +40,7 @@ Foam::Distribution<Type>::Distribution()
 template<class Type>
 Foam::Distribution<Type>::Distribution(const Type& binWidth)
 :
-    List< Map<label> >(pTraits<Type>::nComponents),
+    List< Map<scalar> >(pTraits<Type>::nComponents),
     binWidth_(binWidth)
 {}
 
@@ -53,7 +51,7 @@ Foam::Distribution<Type>::Distribution(const Type& binWidth)
 //     const cmptType& binWidth
 // )
 // :
-//     List< Map<label> >(pTraits<Type>::nComponents),
+//     List< Map<scalar> >(pTraits<Type>::nComponents),
 //     binWidth_(binWidth*pTraits<Type>::one)
 // {}
 
@@ -61,7 +59,7 @@ Foam::Distribution<Type>::Distribution(const Type& binWidth)
 template<class Type>
 Foam::Distribution<Type>::Distribution(const Distribution<Type>& d)
 :
-    List< Map<label> >(static_cast< const List< Map<label> >& >(d)),
+    List< Map<scalar> >(static_cast< const List< Map<scalar> >& >(d)),
     binWidth_(d.binWidth())
 {}
 
@@ -80,13 +78,145 @@ Foam::Distribution<Type>::~Distribution()
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 template<class Type>
-void Foam::Distribution<Type>::add(const Type& valueToAdd)
+Foam::scalar Foam::Distribution<Type>::totalWeight(direction cmpt) const
 {
-    for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; cmpt++)
-    {
-        Map<label>& cmptDistribution = (*this)[cmpt];
+    const Map<scalar>& cmptDistribution = (*this)[cmpt];
 
-        Map<label>::iterator iter(cmptDistribution.begin());
+    scalar sumOfWeights = 0.0;
+
+    forAllConstIter(Map<scalar>, cmptDistribution, iter)
+    {
+        sumOfWeights += iter();
+    }
+
+    return sumOfWeights;
+}
+
+
+template<class Type>
+inline Type Foam::Distribution<Type>::mean() const
+{
+    Type meanValue(pTraits<Type>::zero);
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        const Map<scalar>& cmptDistribution = (*this)[cmpt];
+
+        scalar totalCmptWeight = totalWeight(cmpt);
+
+        List<label> keys = cmptDistribution.sortedToc();
+
+        forAll(keys,k)
+        {
+            label key = keys[k];
+
+            setComponent(meanValue, cmpt) +=
+                (0.5 + scalar(key))
+               *component(binWidth_, cmpt)
+               *cmptDistribution[key]
+               /totalCmptWeight;
+        }
+    }
+
+    return meanValue;
+}
+
+
+template<class Type>
+inline Type Foam::Distribution<Type>::median()
+{
+    Type medianValue(pTraits<Type>::zero);
+
+    List< List < Pair<scalar> > > normDistribution = normalised();
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        List<Pair<scalar> >& normDist = normDistribution[cmpt];
+
+        scalar cumulative = 0.0;
+
+        if (normDist.size())
+        {
+            if (normDist.size() == 1)
+            {
+                setComponent(medianValue, cmpt) = normDist[0].first();
+            }
+            else if
+            (
+                normDist.size() > 1
+             && normDist[0].second()*component(binWidth_, cmpt) > 0.5
+            )
+            {
+                scalar xk = normDist[1].first();
+
+                scalar xkm1 = normDist[0].first();
+
+                scalar Sk =
+                    (normDist[0].second() + normDist[1].second())
+                   *component(binWidth_, cmpt);
+
+                scalar Skm1 = normDist[0].second()*component(binWidth_, cmpt);
+
+                setComponent(medianValue, cmpt) =
+                    (0.5 - Skm1)*(xk - xkm1)/(Sk - Skm1) + xkm1;
+            }
+            else
+            {
+                label lastNonZeroIndex = 0;
+
+                forAll(normDist,nD)
+                {
+                    if
+                    (
+                        cumulative
+                      + (normDist[nD].second()*component(binWidth_, cmpt))
+                      > 0.5
+                    )
+                    {
+                        scalar xk = normDist[nD].first();
+
+                        scalar xkm1 = normDist[lastNonZeroIndex].first();
+
+                        scalar Sk =
+                            cumulative
+                          + (normDist[nD].second()*component(binWidth_, cmpt));
+
+                        scalar Skm1 = cumulative;
+
+                        setComponent(medianValue, cmpt) =
+                            (0.5 - Skm1)*(xk - xkm1)/(Sk - Skm1) + xkm1;
+
+                        break;
+                    }
+                    else if (normDist[nD].second() > 0.0)
+                    {
+                        cumulative +=
+                            normDist[nD].second()*component(binWidth_, cmpt);
+
+                        lastNonZeroIndex = nD;
+                    }
+                }
+            }
+        }
+
+    }
+
+    return medianValue;
+}
+
+
+template<class Type>
+void Foam::Distribution<Type>::add
+(
+    const Type& valueToAdd,
+    const Type& weight
+)
+{
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        Map<scalar>& cmptDistribution = (*this)[cmpt];
+
+        Map<scalar>::iterator iter(cmptDistribution.begin());
 
         label n =
             label(component(valueToAdd, cmpt)/component(binWidth_, cmpt))
@@ -96,23 +226,23 @@ void Foam::Distribution<Type>::add(const Type& valueToAdd)
 
         if (iter == cmptDistribution.end())
         {
-            cmptDistribution.insert(n,1);
+            cmptDistribution.insert(n, component(weight, cmpt));
         }
         else
         {
-            cmptDistribution[n]++;
+            cmptDistribution[n] += component(weight, cmpt);
         }
 
-        if (cmptDistribution[n] < 0)
-        {
-            FatalErrorIn("Distribution::add(const scalar valueToAdd)")
-                << "Accumulated Distribution value has become negative: "
-                << "bin = " << (0.5 + scalar(n))*component(binWidth_, cmpt)
-                << ", value = " << cmptDistribution[n]
-                << ". This is most likely to be because too many samples "
-                << "have been added to a bin and the label has 'rolled round'"
-                << abort(FatalError);
-        }
+        // if (cmptDistribution[n] < 0)
+        // {
+        //     FatalErrorIn("Distribution::add(const scalar valueToAdd)")
+        //         << "Accumulated Distribution value has become negative: "
+        //         << "bin = " << (0.5 + scalar(n))*component(binWidth_, cmpt)
+        //         << ", value = " << cmptDistribution[n]
+        //         << ". This is most likely to be because too many samples "
+        //         << "have been added to a bin and the weight has 'rolled round'"
+        //         << abort(FatalError);
+        // }
     }
 }
 
@@ -120,11 +250,11 @@ void Foam::Distribution<Type>::add(const Type& valueToAdd)
 template<class Type>
 void Foam::Distribution<Type>::insertMissingKeys()
 {
-    for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
     {
-        Map<label>& cmptDistribution = (*this)[cmpt];
+        Map<scalar>& cmptDistribution = (*this)[cmpt];
 
-        Map<label>::iterator iter(cmptDistribution.begin());
+        Map<scalar>::iterator iter(cmptDistribution.begin());
 
         List<label> keys = cmptDistribution.sortedToc();
 
@@ -146,19 +276,57 @@ void Foam::Distribution<Type>::insertMissingKeys()
 
 template<class Type>
 Foam::List< Foam::List< Foam::Pair<Foam::scalar> > >Foam::
-Distribution<Type>::raw()
+Distribution<Type>::normalised()
 {
-    List< List < Pair<scalar> > > rawDistributions(pTraits<Type>::nComponents);
+    List< List < Pair<scalar> > > normDistribution(pTraits<Type>::nComponents);
 
     insertMissingKeys();
 
-    for (direction cmpt=0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
     {
-        Map<label>& cmptDistribution = (*this)[cmpt];
+        Map<scalar>& cmptDistribution = (*this)[cmpt];
+
+        scalar totalCmptWeight = totalWeight(cmpt);
 
         List<label> keys = cmptDistribution.sortedToc();
 
-        List<Pair<scalar> >& rawDist = rawDistributions[cmpt];
+        List<Pair<scalar> >& normDist = normDistribution[cmpt];
+
+        normDist.setSize(keys.size());
+
+        forAll(keys,k)
+        {
+            label key = keys[k];
+
+            normDist[k].first() =
+                (0.5 + scalar(key))*component(binWidth_, cmpt);
+
+            normDist[k].second() =
+                cmptDistribution[key]
+               /totalCmptWeight
+               /component(binWidth_, cmpt);
+        }
+    }
+
+    return normDistribution;
+}
+
+
+template<class Type>
+Foam::List< Foam::List< Foam::Pair<Foam::scalar> > >Foam::
+Distribution<Type>::raw()
+{
+    List< List < Pair<scalar> > > rawDistribution(pTraits<Type>::nComponents);
+
+    insertMissingKeys();
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        Map<scalar>& cmptDistribution = (*this)[cmpt];
+
+        List<label> keys = cmptDistribution.sortedToc();
+
+        List<Pair<scalar> >& rawDist = rawDistribution[cmpt];
 
         rawDist.setSize(keys.size());
 
@@ -168,11 +336,48 @@ Distribution<Type>::raw()
 
             rawDist[k].first() = (0.5 + scalar(key))*component(binWidth_, cmpt);
 
-            rawDist[k].second() = scalar(cmptDistribution[key]);
+            rawDist[k].second() = cmptDistribution[key];
         }
     }
 
-    return rawDistributions;
+    return rawDistribution;
+}
+
+
+template<class Type>
+void Foam::Distribution<Type>::write
+(
+    const fileName& filePrefix,
+    const List< List<Pair<scalar> > >& pairs
+) const
+{
+    if (pairs.size() != pTraits<Type>::nComponents)
+    {
+        FatalErrorIn
+        (
+            "Distribution::write"
+            "("
+                "const fileName& filePrefix,"
+                "const List< List<Pair<scalar> > >& pairs"
+            ")"
+        )
+            << "List of pairs (" << pairs.size()
+            << ") is not the same size as the number of components ("
+            << pTraits<Type>::nComponents << ")." << nl
+            << abort(FatalError);
+    }
+
+    for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+    {
+        const List<Pair<scalar> >& cmptPairs = pairs[cmpt];
+
+        OFstream os(filePrefix + '_' + pTraits<Type>::componentNames[cmpt]);
+
+        forAll(cmptPairs, i)
+        {
+            os  << cmptPairs[i].first() << ' ' << cmptPairs[i].second() << nl;
+        }
+    }
 }
 
 
@@ -195,7 +400,7 @@ void Foam::Distribution<Type>::operator=
             << abort(FatalError);
     }
 
-    List< Map<label> >::operator=(rhs);
+    List< Map<scalar> >::operator=(rhs);
 
     binWidth_ = rhs.binWidth();
 }
@@ -211,7 +416,7 @@ Foam::Ostream& Foam::operator<<
 )
 {
     os  << d.binWidth_
-        << static_cast<const List< Map<label> >& >(d);
+        << static_cast<const List< Map<scalar> >& >(d);
 
     // Check state of Ostream
     os.check
diff --git a/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.H b/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.H
index 0f60ee5511..377ffae830 100644
--- a/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.H
+++ b/src/OpenFOAM/containers/HashTables/Map/Distribution/Distribution.H
@@ -60,37 +60,19 @@ Ostream& operator<<(Ostream&, const Distribution<Type>&);
 template<class Type>
 class Distribution
 :
-    public List< Map<label> >
+    public List< Map<scalar> >
 {
     // Private data
 
         Type binWidth_;
 
 
-    // Private Member Functions
-
-        // //- Disallow default bitwise copy construct
-        // Distribution(const Distribution<Type>&);
-
-        // //- Disallow default bitwise assignment
-        // void operator=(const Distribution<Type>&);
-
-
 public:
 
     //- Component type
     typedef typename pTraits<Type>::cmptType cmptType;
 
 
-    // Static data members
-
-        // static void write
-        // (
-        //     const fileName& file,
-        //     const List<Pair<scalar> >& pairs
-        // );
-
-
     // Constructors
 
         //- Construct null
@@ -99,7 +81,7 @@ public:
         //- Construct from separate binWidth for each component
         Distribution(const Type& binWidth);
 
-        // //- Construct from single binWidth for each component
+        //- Construct from single binWidth for each component
         // Distribution(const cmptType& binWidth);
 
         //- Construct as copy
@@ -111,17 +93,24 @@ public:
 
     // Member Functions
 
-        // Type totalEntries() const;
+        scalar totalWeight(direction cmpt) const;
 
-        // Type mean() const;
+        Type mean() const;
 
-        // Type median();
+        // From http://mathworld.wolfram.com/StatisticalMedian.html
+        // The statistical median is the value of the Distribution
+        // variable where the cumulative Distribution = 0.5.
+        Type median();
 
-        void add(const Type& valueToAdd);
+        void add
+        (
+            const Type& valueToAdd,
+            const Type& weight = pTraits<Type>::one
+        );
 
         void insertMissingKeys();
 
-        // List<Pair<scalar> > normalised();
+        List< List<Pair<scalar> > > normalised();
 
         // List<Pair<scalar> > normalisedMinusMean();
 
@@ -133,15 +122,19 @@ public:
 
         inline const Type& binWidth() const;
 
+    // Write
+
+        void write
+        (
+            const fileName& filePrefix,
+            const List< List<Pair<scalar> > >& pairs
+        ) const;
+
+
     // Member Operators
 
         void operator=(const Distribution<Type>&);
 
-
-    // Friend Functions
-
-    // Friend Operators
-
     // IOstream Operators
 
         friend Ostream& operator<< <Type>