Skip to content
GitLab
Commit f666f549 authored by Andrew Heather's avatar Andrew Heather
Browse files

ENH: distributionModels - refactored binned and general models

parent c0f1a87d
Hide whitespace changes
Inline Side-by-side
src/lagrangian/distributionModels/binned/binned.C
View file @ f666f549
......@@ -39,7 +39,38 @@ namespace Foam
const char* Foam::distributionModels::binned::header =
"minValue maxValue (binBoundaries) (values)";
"(bin probability)";
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::distributionModels::binned::initialise()
{
const label nSample(xy_.size());
forAll(xy_, bini)
{
xy_[bini][1] /= scalar(nSample);
}
// Convert values to integral values
for (label bini = 1; bini < nSample; ++bini)
{
xy_[bini][1] += xy_[bini - 1][1];
}
// Calculate the mean value
label bini = 0;
forAll(xy_, i)
{
if (xy_[i][1] > 0.5)
{
bini = i;
break;
}
}
meanValue_ = xy_[bini][1];
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
......@@ -51,121 +82,93 @@ Foam::distributionModels::binned::binned
)
:
distributionModel(typeName, dict, rndGen),
minValue_(readScalar(distributionModelDict_.lookup("minValue"))),
maxValue_(readScalar(distributionModelDict_.lookup("maxValue"))),
meanValue_(0),
binBoundaries_(distributionModelDict_.lookup("binBoundaries")),
values_(distributionModelDict_.lookup("values"))
xy_(distributionModelDict_.lookup("distribution")),
meanValue_(0)
{
if (minValue_ < 0)
{
FatalErrorInFunction
<< "Minimum value must be greater than zero. "
<< "Supplied minValue = " << minValue_
<< exit(FatalError);
}
if (maxValue_ < minValue_)
if (maxValue() < minValue())
{
FatalErrorInFunction
<< "Maximum value is smaller than the minimum value:"
<< " maxValue = " << maxValue_ << ", minValue = " << minValue_
<< " maxValue = " << maxValue()
<< ", minValue = " << minValue()
<< exit(FatalError);
}
if (binBoundaries_.size() != values_.size() + 1)
{
FatalErrorInFunction
<< "Number of bin boundaries must be number of values + 1"
<< " number of binBoundaries : " << binBoundaries_.size() << nl
<< " number of values : " << values_.size() << nl
<< exit(FatalError);
}
initialise();
}
Foam::distributionModels::binned::binned
(
const scalarField& values,
const UList<scalar>& sampleData,
const scalar binWidth,
cachedRandom& rndGen
)
:
distributionModel(typeName, dictionary::null, rndGen),
minValue_(0),
maxValue_(0),
meanValue_(0),
binBoundaries_(),
values_()
xy_(),
meanValue_(0)
{
label bin0 = floor(min(values)/binWidth);
label bin1 = ceil(max(values)/binWidth);
label nBin = bin1 - bin0;
scalar minValue = GREAT;
scalar maxValue = -GREAT;
forAll(sampleData, i)
{
minValue = min(minValue, sampleData[i]);
maxValue = max(maxValue, sampleData[i]);
}
minValue_ = bin0*binWidth;
maxValue_ = bin1*binWidth;
const label bin0 = floor(minValue/binWidth);
const label bin1 = ceil(maxValue/binWidth);
const label nBin = bin1 - bin0;
if (nBin == 0)
{
WarningInFunction
<< "Data cannot be binned - zero bins generated" << nl
<< " Bin width : " << binWidth << nl
<< " Values : " << values
<< " Bin width : " << binWidth << nl
<< " Sample data : " << sampleData
<< endl;
return;
}
// Create the bins
binBoundaries_.setSize(nBin + 1);
forAll(binBoundaries_, binI)
// Populate bin boundaries and initialise occurrences
xy_.setSize(nBin);
forAll(xy_, bini)
{
binBoundaries_[binI] = (bin0 + binI)*binWidth;
xy_[bini][0] = (bin0 + bini)*binWidth;
xy_[bini][1] = 0;
}
// Bin and normalise the data
values_.setSize(nBin, 0);
forAll(values, i)
// Bin the data
forAll(sampleData, i)
{
label binI = floor(values[i]/binWidth) - bin0;
values_[binI]++;
}
values_ /= sum(values_);
// Choose the nearest bin
label bini = floor(sampleData[i]/binWidth) - bin0;
label binii = min(bini + 1, nBin - 1);
// Calculate the mean value
scalar sum = 0;
label valueI = 0;
forAll(values_, i)
{
sum += values_[i];
if (sum > 0.5)
scalar d1 = mag(sampleData[i] - xy_[bini][0]);
scalar d2 = mag(xy_[binii][0] - sampleData[i]);
if (d1 < d2)
{
valueI = i;
break;
xy_[bini][1]++;
}
else
{
xy_[binii][1]++;
}
}
if (valueI == 0)
{
meanValue_ = binBoundaries_[valueI];
}
else
{
meanValue_ =
binBoundaries_[valueI]
+ (0.5 - binBoundaries_[valueI])
/ (binBoundaries_[valueI + 1] - binBoundaries_[valueI]);
}
initialise();
}
Foam::distributionModels::binned::binned(const binned& p)
:
distributionModel(p),
minValue_(p.minValue_),
maxValue_(p.maxValue_),
binBoundaries_(p.binBoundaries_),
values_(p.values_)
xy_(p.xy_),
meanValue_(p.meanValue_)
{}
......@@ -181,31 +184,36 @@ Foam::scalar Foam::distributionModels::binned::sample() const
{
scalar y = rndGen_.sample01<scalar>();
scalar sum = 0;
forAll(values_, i)
for (label i = 0; i < xy_.size() - 1; ++i)
{
sum += values_[i];
if (sum > y)
if (xy_[i][1] > y)
{
scalar binMin = binBoundaries_[i];
scalar binMax = binBoundaries_[i + 1];
return binMin + (1 - (sum - y)/values_[i])*(binMax - binMin);
scalar d1 = y - xy_[i][1];
scalar d2 = xy_[i+1][1] - y;
if (d1 < d2)
{
return xy_[i][0];
}
else
{
return xy_[i+1][0];
}
}
}
return maxValue_;
return xy_.last()[0];
}
Foam::scalar Foam::distributionModels::binned::minValue() const
{
return minValue_;
return xy_.first()[0];
}
Foam::scalar Foam::distributionModels::binned::maxValue() const
{
return maxValue_;
return xy_.last()[0];
}
......@@ -218,20 +226,14 @@ Foam::scalar Foam::distributionModels::binned::meanValue() const
void Foam::distributionModels::binned::readData(Istream& is)
{
// distributionModel::readData(is);
is >> minValue_
>> maxValue_
>> binBoundaries_
>> values_;
is >> xy_;
}
void Foam::distributionModels::binned::writeData(Ostream& os) const
{
// distributionModel::writeData(os);
os << minValue_ << token::SPACE
<< maxValue_ << token::SPACE
<< binBoundaries_ << token::SPACE
<< values_;
os << xy_ ;
}
......@@ -242,11 +244,7 @@ Foam::dictionary Foam::distributionModels::binned::writeDict
{
// dictionary dict = distributionModel::writeDict(dictName);
dictionary dict(dictName);
dict.add("minValue", minValue_);
dict.add("maxValue", maxValue_);
dict.add("binBoundaries", binBoundaries_);
dict.add("values", values_);
dict.add("distribution", xy_);
return dict;
}
......@@ -255,26 +253,26 @@ Foam::dictionary Foam::distributionModels::binned::writeDict
void Foam::distributionModels::binned::readDict(const dictionary& dict)
{
// distributionModel::readDict(dict);
dict.lookup("minValue") >> minValue_;
dict.lookup("maxValue") >> maxValue_;
dict.lookup("binBoundaries") >> binBoundaries_;
dict.lookup("values") >> values_;
dict.lookup("distribution") >> xy_;
}
Foam::Ostream& Foam::operator<<(Ostream& os, const distributionModels::binned& b)
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const distributionModels::binned& b
)
{
os.check("Foam::Ostream& Foam::operator<<(Ostream&, const binned&)");
os.check(FUNCTION_NAME);
b.writeData(os);
return os;
}
Foam::Istream& Foam::operator>>(Istream& is, distributionModels::binned& b)
Foam::Istream& Foam::operator>>(Istream& is, distributionModels::binned& b)
{
is.check("Foam::Istream& Foam::operator>>(Istream&, binned&)");
is.check(FUNCTION_NAME);
b.readData(is);
return is;
......
src/lagrangian/distributionModels/binned/binned.H
View file @ f666f549
......@@ -22,10 +22,11 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::binned
Foam::distributionModels::binned
Description
A binned distribution model
A binned distribution model where the random sample is taken from a
discrete set of bins
SourceFiles
binned.C
......@@ -69,23 +70,19 @@ class binned
{
// Private data
//- Distribution minimum
scalar minValue_;
typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
//- Distribution maximum
scalar maxValue_;
// List of (bin probability)
List<pair> xy_;
//- Distribution mean value
scalar meanValue_;
// Model coefficients
// Private member functions
//- List of bin boundaries
Field<scalar> binBoundaries_;
//- List of normalised bin values
Field<scalar> values_;
//- Initialise the distribution parameters
void initialise();
public:
......@@ -104,7 +101,7 @@ public:
//- Construct from components
binned
(
const scalarField& values,
const UList<scalar>& sampleData,
const scalar binWidth,
cachedRandom& rndGen
);
......
src/lagrangian/distributionModels/general/general.C
View file @ f666f549
......@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -37,6 +37,40 @@ namespace Foam
}
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::distributionModels::general::initialise()
{
const label nEntries = xy_.size();
integral_.setSize(nEntries);
// Normalize the cumulative distribution
integral_[0] = 0.0;
for (label i = 1; i < nEntries; i++)
{
scalar k = (xy_[i][1] - xy_[i-1][1])/(xy_[i][0] - xy_[i-1][0]);
scalar d = xy_[i-1][1] - k*xy_[i-1][0];
scalar y1 = xy_[i][0]*(0.5*k*xy_[i][0] + d);
scalar y0 = xy_[i-1][0]*(0.5*k*xy_[i-1][0] + d);
scalar area = y1 - y0;
integral_[i] = area + integral_[i-1];
}
scalar sumArea = integral_.last();
meanValue_ = sumArea/(maxValue() - minValue());
for (label i=0; i < nEntries; i++)
{
xy_[i][1] /= sumArea;
integral_[i] /= sumArea;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::distributionModels::general::general
......@@ -47,39 +81,67 @@ Foam::distributionModels::general::general
:
distributionModel(typeName, dict, rndGen),
xy_(distributionModelDict_.lookup("distribution")),
nEntries_(xy_.size()),
minValue_(xy_[0][0]),
maxValue_(xy_[nEntries_-1][0]),
meanValue_(0.0),
integral_(nEntries_)
integral_()
{
check();
// normalize the cumulative distributionModel
initialise();
}
integral_[0] = 0.0;
for (label i=1; i<nEntries_; i++)
Foam::distributionModels::general::general
(
const UList<scalar>& sampleData,
const scalar binWidth,
cachedRandom& rndGen
)
:
distributionModel(typeName, dictionary::null, rndGen),
xy_(),
meanValue_(0.0),
integral_()
{
scalar minValue = GREAT;
scalar maxValue = -GREAT;
forAll(sampleData, i)
{
minValue = min(minValue, sampleData[i]);
maxValue = max(maxValue, sampleData[i]);
}
scalar k = (xy_[i][1] - xy_[i-1][1])/(xy_[i][0] - xy_[i-1][0]);
scalar d = xy_[i-1][1] - k*xy_[i-1][0];
scalar y1 = xy_[i][0]*(0.5*k*xy_[i][0] + d);
scalar y0 = xy_[i-1][0]*(0.5*k*xy_[i-1][0] + d);
scalar area = y1 - y0;
label bin0 = floor(minValue/binWidth);
label bin1 = ceil(maxValue/binWidth);
label nEntries = bin1 - bin0;
integral_[i] = area + integral_[i-1];
if (nEntries == 0)
{
WarningInFunction
<< "Data cannot be binned - zero bins generated" << nl
<< " Bin width : " << binWidth << nl
<< " Sample data : " << sampleData
<< endl;
return;
}
scalar sumArea = integral_.last();
xy_.setSize(nEntries);
meanValue_ = sumArea/(maxValue_ - minValue_);
// Populate bin boundaries and initialise occurrences
for (label bini = 0; bini < nEntries; ++bini)
{
xy_[bini][0] = (bin0 + bini)*binWidth;
xy_[bini][1] = 0;
}
for (label i=0; i<nEntries_; i++)
// Populate occurrences
forAll(sampleData, i)
{
xy_[i][1] /= sumArea;
integral_[i] /= sumArea;
label bini = floor(sampleData[i]/binWidth) - bin0;
xy_[bini][1]++;
}
initialise();
}
......@@ -87,9 +149,6 @@ Foam::distributionModels::general::general(const general& p)
:
distributionModel(p),
xy_(p.xy_),
nEntries_(p.nEntries_),
minValue_(p.minValue_),
maxValue_(p.maxValue_),
integral_(p.integral_)
{}
......@@ -106,8 +165,8 @@ Foam::scalar Foam::distributionModels::general::sample() const
{
scalar y = rndGen_.sample01<scalar>();
// find the interval where y is in the table
label n=1;
// Find the interval where y is in the table
label n = 1;
while (integral_[n] <= y)
{
n++;
......@@ -119,7 +178,7 @@ Foam::scalar Foam::distributionModels::general::sample() const
scalar alpha = y + xy_[n-1][0]*(0.5*k*xy_[n-1][0] + d) - integral_[n-1];
scalar x = 0.0;
// if k is small it is a linear equation, otherwise it is of second order
// If k is small it is a linear equation, otherwise it is of second order
if (mag(k) > SMALL)
{
scalar p = 2.0*d/k;
......@@ -148,13 +207,13 @@ Foam::scalar Foam::distributionModels::general::sample() const
Foam::scalar Foam::distributionModels::general::minValue() const
{
return minValue_;
return xy_.first()[0];
}
Foam::scalar Foam::distributionModels::general::maxValue() const
{
return maxValue_;
return xy_.last()[0];
}
......@@ -164,4 +223,62 @@ Foam::scalar Foam::distributionModels::general::meanValue() const
}
void Foam::distributionModels::general::readData(Istream& is)
{
// distributionModel::readData(is);
is >> xy_;
initialise();
}
void Foam::distributionModels::general::writeData(Ostream& os) const
{
// distributionModel::writeData(os);
os << xy_;
}
Foam::dictionary Foam::distributionModels::general::writeDict
(
const word& dictName
) const
{
// dictionary dict = distributionModel::writeDict(dictName);
dictionary dict(dictName);
dict.add("distribution", xy_);
return dict;
}
void Foam::distributionModels::general::readDict(const dictionary& dict)
{
// distributionModel::readDict(dict);
dict.lookup("distribution") >> xy_;
initialise();
}
Foam::Ostream& Foam::operator<<
(
Ostream& os,
const distributionModels::general& b
)
{
os.check(FUNCTION_NAME);
b.writeData(os);
return os;
}
Foam::Istream& Foam::operator>>(Istream