From abd1ee0f1db6bc844b9053634952cf8d5f933853 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Thu, 11 Feb 2010 17:21:16 +0000
Subject: [PATCH] ENH: Distribution. Adding operator+ to allow
reduce(dist, sumOp< Distribution<scalar> >());
operations for parallel data. Combines distributions using the
coarsest binWidth.
Changed write function to write normalised and raw to the same file.
BUG: Distribution. Corrected inaccurate median calculation.
---
.../test/Distribution/DistributionTest.C | 137 +++++---
.../Lists/Distribution/Distribution.C | 294 +++++++++++-------
.../Lists/Distribution/Distribution.H | 51 +--
3 files changed, 309 insertions(+), 173 deletions(-)
diff --git a/applications/test/Distribution/DistributionTest.C b/applications/test/Distribution/DistributionTest.C
index e59cb1af1ea..3d181b22fbb 100644
--- a/applications/test/Distribution/DistributionTest.C
+++ b/applications/test/Distribution/DistributionTest.C
@@ -45,7 +45,8 @@ Description
#include "Distribution.H"
#include "Random.H"
#include "dimensionedTypes.H"
-
+#include "argList.H"
+#include "PstreamReduceOps.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -53,6 +54,8 @@ using namespace Foam;
int main(int argc, char *argv[])
{
+ # include "setRootCase.H"
+
Random R(918273);
{
@@ -76,7 +79,67 @@ int main(int argc, char *argv[])
<< "Median " << dS.median()
<< endl;
- dS.write("Distribution_scalar_test", dS.normalised());
+ dS.write("Distribution_scalar_test_1");
+
+ Distribution<scalar> dS2(scalar(1e-2));
+
+ Info<< nl << "Distribution<scalar>" << nl
+ << "Sampling "
+ << randomDistributionTestSize
+ << " times from GaussNormal distribution."
+ << endl;
+
+ for (label i = 0; i < randomDistributionTestSize; i++)
+ {
+ dS2.add(1.5*R.GaussNormal() -6.0);
+ }
+
+ Info<< "Mean " << dS2.mean() << nl
+ << "Median " << dS2.median()
+ << endl;
+
+ dS2.write("Distribution_scalar_test_2");
+
+ Info<< nl << "Adding previous two Distribution<scalar>" << endl;
+
+ dS = dS + dS2;
+
+ dS.write("Distribution_scalar_test_1+2");
+ }
+
+ if (Pstream::parRun())
+ {
+ // scalar in parallel
+ label randomDistributionTestSize = 100000000;
+
+ Distribution<scalar> dS(scalar(1e-1));
+
+ Pout<< "Distribution<scalar>" << nl
+ << "Sampling "
+ << randomDistributionTestSize
+ << " times from uniform distribution."
+ << endl;
+
+ for (label i = 0; i < randomDistributionTestSize; i++)
+ {
+ dS.add(R.scalar01() + 10*Pstream::myProcNo());
+ }
+
+ Pout<< "Mean " << dS.mean() << nl
+ << "Median " << dS.median()
+ << endl;
+
+ reduce(dS, sumOp< Distribution<scalar> >());
+
+ if (Pstream::master())
+ {
+ Info<< "Reducing parallel Distribution<scalar>" << nl
+ << "Mean " << dS.mean() << nl
+ << "Median " << dS.median()
+ << endl;
+
+ dS.write("Distribution_scalar_test_parallel_reduced");
+ }
}
{
@@ -114,40 +177,40 @@ int main(int argc, char *argv[])
<< "Median " << dV.median()
<< endl;
- dV.write("Distribution_vector_test", dV.normalised());
+ dV.write("Distribution_vector_test");
}
- {
- // 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());
- }
+ // {
+ // // 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");
+ // }
{
// tensor
@@ -170,7 +233,7 @@ int main(int argc, char *argv[])
<< "Median " << dT.median()
<< endl;
- dT.write("Distribution_tensor_test", dT.normalised());
+ dT.write("Distribution_tensor_test");
}
{
@@ -194,7 +257,7 @@ int main(int argc, char *argv[])
<< "Median " << dSyT.median()
<< endl;
- dSyT.write("Distribution_symmTensor_test", dSyT.normalised());
+ dSyT.write("Distribution_symmTensor_test");
}
{
@@ -218,7 +281,7 @@ int main(int argc, char *argv[])
<< "Median " << dSpT.median()
<< endl;
- dSpT.write("Distribution_sphericalTensor_test", dSpT.normalised());
+ dSpT.write("Distribution_sphericalTensor_test");
}
Info<< nl << "End" << nl << endl;
diff --git a/src/OpenFOAM/containers/Lists/Distribution/Distribution.C b/src/OpenFOAM/containers/Lists/Distribution/Distribution.C
index 263853428dd..d242feb11a1 100644
--- a/src/OpenFOAM/containers/Lists/Distribution/Distribution.C
+++ b/src/OpenFOAM/containers/Lists/Distribution/Distribution.C
@@ -28,7 +28,73 @@ License
#include "OFstream.H"
#include "ListOps.H"
-// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Distribution<Type>::Distribution()
+:
+ List< List<scalar> >(pTraits<Type>::nComponents),
+ binWidth_(pTraits<Type>::one),
+ listStarts_(pTraits<Type>::nComponents, 0)
+{}
+
+
+template<class Type>
+Foam::Distribution<Type>::Distribution(const Type& binWidth)
+:
+ List< List<scalar> >(pTraits<Type>::nComponents),
+ binWidth_(binWidth),
+ listStarts_(pTraits<Type>::nComponents, 0)
+{}
+
+
+template<class Type>
+Foam::Distribution<Type>::Distribution(const Distribution<Type>& d)
+:
+ List< List<scalar> >(static_cast< const List< List<scalar> >& >(d)),
+ binWidth_(d.binWidth()),
+ listStarts_(d.listStarts())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Distribution<Type>::~Distribution()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::scalar Foam::Distribution<Type>::totalWeight(direction cmpt) const
+{
+ const List<scalar>& cmptDistribution = (*this)[cmpt];
+
+ scalar sumOfWeights = 0.0;
+
+ forAll(cmptDistribution, i)
+ {
+ sumOfWeights += cmptDistribution[i];
+ }
+
+ return sumOfWeights;
+}
+
+
+template<class Type>
+Foam::List<Foam::label> Foam::Distribution<Type>::keys(direction cmpt) const
+{
+ List<label> keys = identity((*this)[cmpt].size());
+
+ forAll(keys, k)
+ {
+ keys[k] += listStarts_[cmpt];
+ }
+
+ return keys;
+}
+
template<class Type>
Foam::label Foam::Distribution<Type>::index
@@ -128,73 +194,6 @@ Foam::Pair<Foam::label> Foam::Distribution<Type>::validLimits
}
-// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::Distribution<Type>::Distribution()
-:
- List< List<scalar> >(pTraits<Type>::nComponents),
- binWidth_(pTraits<Type>::one),
- listStarts_(pTraits<Type>::nComponents, 0)
-{}
-
-
-template<class Type>
-Foam::Distribution<Type>::Distribution(const Type& binWidth)
-:
- List< List<scalar> >(pTraits<Type>::nComponents),
- binWidth_(binWidth),
- listStarts_(pTraits<Type>::nComponents, 0)
-{}
-
-
-template<class Type>
-Foam::Distribution<Type>::Distribution(const Distribution<Type>& d)
-:
- List< List<scalar> >(static_cast< const List< List<scalar> >& >(d)),
- binWidth_(d.binWidth()),
- listStarts_(d.listStarts())
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::Distribution<Type>::~Distribution()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
-
-template<class Type>
-Foam::scalar Foam::Distribution<Type>::totalWeight(direction cmpt) const
-{
- const List<scalar>& cmptDistribution = (*this)[cmpt];
-
- scalar sumOfWeights = 0.0;
-
- forAll(cmptDistribution, i)
- {
- sumOfWeights += cmptDistribution[i];
- }
-
- return sumOfWeights;
-}
-
-template<class Type>
-Foam::List<Foam::label> Foam::Distribution<Type>::keys(direction cmpt) const
-{
- List<label> keys = identity((*this)[cmpt].size());
-
- forAll(keys, k)
- {
- keys[k] += listStarts_[cmpt];
- }
-
- return keys;
-}
-
-
template<class Type>
Type Foam::Distribution<Type>::mean() const
{
@@ -225,7 +224,7 @@ Type Foam::Distribution<Type>::mean() const
template<class Type>
-Type Foam::Distribution<Type>::median()
+Type Foam::Distribution<Type>::median() const
{
Type medianValue(pTraits<Type>::zero);
@@ -247,18 +246,17 @@ Type Foam::Distribution<Type>::median()
&& normDist[0].second()*component(binWidth_, cmpt) > 0.5
)
{
- scalar xk = normDist[1].first();
-
- scalar xkm1 = normDist[0].first();
+ scalar xk =
+ normDist[0].first()
+ + 0.5*component(binWidth_, cmpt);
- scalar Sk =
- (normDist[0].second() + normDist[1].second())
- *component(binWidth_, cmpt);
+ scalar xkm1 =
+ normDist[0].first()
+ - 0.5*component(binWidth_, cmpt);
- scalar Skm1 = normDist[0].second()*component(binWidth_, cmpt);
+ scalar Sk = (normDist[0].second())*component(binWidth_, cmpt);
- setComponent(medianValue, cmpt) =
- (0.5 - Skm1)*(xk - xkm1)/(Sk - Skm1) + xkm1;
+ setComponent(medianValue, cmpt) = 0.5*(xk - xkm1)/(Sk) + xkm1;
}
else
{
@@ -266,7 +264,7 @@ Type Foam::Distribution<Type>::median()
scalar cumulative = 0.0;
- forAll(normDist,nD)
+ forAll(normDist, nD)
{
if
(
@@ -275,9 +273,13 @@ Type Foam::Distribution<Type>::median()
> 0.5
)
{
- scalar xk = normDist[nD].first();
+ scalar xk =
+ normDist[nD].first()
+ + 0.5*component(binWidth_, cmpt);
- scalar xkm1 = normDist[previousNonZeroIndex].first();
+ scalar xkm1 =
+ normDist[previousNonZeroIndex].first()
+ + 0.5*component(binWidth_, cmpt);
scalar Sk =
cumulative
@@ -314,7 +316,6 @@ void Foam::Distribution<Type>::add
const Type& weight
)
{
-
for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
{
List<scalar>& cmptDistribution = (*this)[cmpt];
@@ -332,13 +333,13 @@ void Foam::Distribution<Type>::add
template<class Type>
Foam::List< Foam::List< Foam::Pair<Foam::scalar> > >Foam::
-Distribution<Type>::normalised()
+Distribution<Type>::normalised() const
{
List< List < Pair<scalar> > > normDistribution(pTraits<Type>::nComponents);
for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
{
- List<scalar>& cmptDistribution = (*this)[cmpt];
+ const List<scalar>& cmptDistribution = (*this)[cmpt];
if (cmptDistribution.empty())
{
@@ -380,13 +381,13 @@ Distribution<Type>::normalised()
template<class Type>
Foam::List< Foam::List< Foam::Pair<Foam::scalar> > >Foam::
-Distribution<Type>::raw()
+Distribution<Type>::raw() const
{
List< List < Pair<scalar> > > rawDistribution(pTraits<Type>::nComponents);
for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
{
- List<scalar>& cmptDistribution = (*this)[cmpt];
+ const List<scalar>& cmptDistribution = (*this)[cmpt];
if (cmptDistribution.empty())
{
@@ -433,37 +434,28 @@ void Foam::Distribution<Type>::clear()
template<class Type>
-void Foam::Distribution<Type>::write
-(
- const fileName& filePrefix,
- const List< List< Pair<scalar> > >& pairs
-) const
+void Foam::Distribution<Type>::write(const fileName& filePrefix) 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);
- }
+ List< List<Pair<scalar> > > rawDistribution = raw();
+
+ List< List < Pair<scalar> > > normDistribution = normalised();
for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
{
- const List< Pair<scalar> >& cmptPairs = pairs[cmpt];
+ const List< Pair<scalar> >& rawPairs = rawDistribution[cmpt];
+
+ const List< Pair<scalar> >& normPairs = normDistribution[cmpt];
OFstream os(filePrefix + '_' + pTraits<Type>::componentNames[cmpt]);
- forAll(cmptPairs, i)
+ os << "# key normalised raw" << endl;
+
+ forAll(normPairs, i)
{
- os << cmptPairs[i].first() << ' ' << cmptPairs[i].second() << nl;
+ os << normPairs[i].first()
+ << ' ' << normPairs[i].second()
+ << ' ' << rawPairs[i].second()
+ << nl;
}
}
}
@@ -491,11 +483,31 @@ void Foam::Distribution<Type>::operator=
List< List<scalar> >::operator=(rhs);
binWidth_ = rhs.binWidth();
+
+ listStarts_ = rhs.listStarts();
}
// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
+template <class Type>
+Foam::Istream& Foam::operator>>
+(
+ Istream& is,
+ Distribution<Type>& d
+)
+{
+ is >> static_cast<List< List<scalar> >&>(d)
+ >> d.binWidth_
+ >> d.listStarts_;
+
+ // Check state of Istream
+ is.check("Istream& operator>>(Istream&, Distribution<Type>&)");
+
+ return is;
+}
+
+
template<class Type>
Foam::Ostream& Foam::operator<<
(
@@ -503,18 +515,68 @@ Foam::Ostream& Foam::operator<<
const Distribution<Type>& d
)
{
- os << d.binWidth_
- << static_cast<const List< List<scalar> >& >(d);
+ os << static_cast<const List< List<scalar> >& >(d)
+ << d.binWidth_ << token::SPACE
+ << d.listStarts_;
// Check state of Ostream
- os.check
- (
- "Ostream& operator<<(Ostream&, "
- "const Distribution&)"
- );
+ os.check("Ostream& operator<<(Ostream&, " "const Distribution&)");
return os;
}
+// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
+
+template <class Type>
+Foam::Distribution<Type> Foam::operator+
+(
+ const Distribution<Type>& d1,
+ const Distribution<Type>& d2
+)
+{
+ // The coarsest binWidth is the sensible choice
+ Distribution<Type> d(max(d1.binWidth(), d2.binWidth()));
+
+ List< List< List < Pair<scalar> > > > rawDists(2);
+
+ rawDists[0] = d1.raw();
+ rawDists[1] = d2.raw();
+
+ forAll(rawDists, rDI)
+ {
+ for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
+ {
+ List<scalar>& cmptDistribution = d[cmpt];
+
+ const List < Pair<scalar> >& cmptRaw = rawDists[rDI][cmpt];
+
+ forAll(cmptRaw, rI)
+ {
+ scalar valueToAdd = cmptRaw[rI].first();
+ scalar cmptWeight = cmptRaw[rI].second();
+
+ label n =
+ label
+ (
+ component(valueToAdd, cmpt)
+ /component(d.binWidth(), cmpt)
+ )
+ - label
+ (
+ neg(component(valueToAdd, cmpt)
+ /component(d.binWidth(), cmpt))
+ );
+
+ label listIndex = d.index(cmpt, n);
+
+ cmptDistribution[listIndex] += cmptWeight;
+ }
+ }
+ }
+
+ return Distribution<Type>(d);
+}
+
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/containers/Lists/Distribution/Distribution.H b/src/OpenFOAM/containers/Lists/Distribution/Distribution.H
index 5dbdb444c46..ef6ef25676c 100644
--- a/src/OpenFOAM/containers/Lists/Distribution/Distribution.H
+++ b/src/OpenFOAM/containers/Lists/Distribution/Distribution.H
@@ -50,6 +50,9 @@ namespace Foam
template<class Type>
class Distribution;
+template<class Type>
+Istream& operator>>(Istream&, Distribution<Type>&);
+
template<class Type>
Ostream& operator<<(Ostream&, const Distribution<Type>&);
@@ -70,15 +73,6 @@ class Distribution
//- The start bin index of each component
List<label> listStarts_;
- // Private Member Functions
-
- //- Return the appropriate List index for the given bin index.
- // Resizes the List if required
- label index(direction cmpt, label n);
-
- //- Returns the indices of the first and last non-zero entries
- Pair<label> validLimits(direction cmpt) const;
-
public:
@@ -109,12 +103,19 @@ public:
List<label> keys(direction cmpt) const;
+ //- Return the appropriate List index for the given bin index.
+ // Resizes the List if required
+ label index(direction cmpt, label n);
+
+ //- Returns the indices of the first and last non-zero entries
+ Pair<label> validLimits(direction cmpt) const;
+
Type mean() const;
// 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();
+ Type median() const;
//- Add a value to the distribution, optionally specifying a weight
void add
@@ -124,10 +125,10 @@ public:
);
//- Return the normalised distribution and bins
- List< List<Pair<scalar> > > normalised();
+ List< List<Pair<scalar> > > normalised() const;
//- Return the distribution of the total bin weights
- List< List < Pair<scalar> > > raw();
+ List< List < Pair<scalar> > > raw() const;
//- Resets the Distribution by clearing the stored lists.
// Leaves the same number of them and the same binWidth.
@@ -144,13 +145,9 @@ public:
// Write
- //- Write the distribution to file. Produces a separate file
- // for each component.
- void write
- (
- const fileName& filePrefix,
- const List< List<Pair<scalar> > >& pairs
- ) const;
+ //- Write the distribution to file: key normalised raw.
+ // Produces a separate file for each component.
+ void write(const fileName& filePrefix) const;
// Member Operators
@@ -159,6 +156,12 @@ public:
// IOstream Operators
+ friend Istream& operator>> <Type>
+ (
+ Istream&,
+ Distribution<Type>&
+ );
+
friend Ostream& operator<< <Type>
(
Ostream&,
@@ -167,7 +170,15 @@ public:
};
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
+
+template<class Type>
+Distribution<Type> operator+
+(
+ const Distribution<Type>&,
+ const Distribution<Type>&
+);
+
} // End namespace Foam
--
GitLab