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.
This commit is contained in:
graham
2009-11-23 12:22:13 +00:00
parent a74f31aace
commit 844685d29d
3 changed files with 435 additions and 88 deletions

View File

@ -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;

View File

@ -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

View File

@ -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>