distribution.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2019 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "distribution.H"
30 #include "OFstream.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36  defineTypeNameAndDebug(distribution, 0);
37 }
38 
39 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
40 
42 (
43  const fileName& file,
44  const List<Pair<scalar>>& pairs
45 )
46 {
47  OFstream os(file);
48 
49  forAll(pairs, i)
50  {
51  os << pairs[i].first() << ' ' << pairs[i].second() << nl;
52  }
53 }
54 
55 
56 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
57 
59 :
60  Map<label>(),
61  binWidth_(1)
62 {}
63 
64 
65 Foam::distribution::distribution(const scalar binWidth)
66 :
67  Map<label>(),
68  binWidth_(binWidth)
69 {}
70 
71 
73 :
74  Map<label>(static_cast<Map<label>>(d)),
75  binWidth_(d.binWidth())
76 {}
77 
78 
79 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
80 
82 {}
83 
84 
85 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
86 
88 {
89  label sumOfEntries = 0;
90 
91  forAllConstIters(*this, iter)
92  {
93  sumOfEntries += iter.val();
94 
95  if (sumOfEntries < 0)
96  {
98  << "Accumulated distribution values total has become negative: "
99  << "sumOfEntries = " << sumOfEntries
100  << ". This is most likely to be because too many samples "
101  << "have been added to the bins and the label has 'rolled "
102  << "round'. Try distribution::approxTotalEntries which "
103  << "returns a scalar." << endl;
104 
105  sumOfEntries = -1;
106 
107  break;
108  }
109  }
110 
111  return sumOfEntries;
112 }
113 
114 
116 {
117  scalar sumOfEntries = 0;
118 
119  forAllConstIters(*this, iter)
120  {
121  sumOfEntries += scalar(iter.val());
122  }
123 
124  return sumOfEntries;
125 }
126 
127 
128 Foam::scalar Foam::distribution::mean() const
129 {
130  scalar runningSum = 0;
131 
132  scalar totEnt = approxTotalEntries();
133 
134  List<label> keys = toc();
135 
136  forAll(keys,k)
137  {
138  label key = keys[k];
139 
140  runningSum +=
141  (0.5 + scalar(key))
142  *binWidth_
143  *scalar((*this)[key])
144  /totEnt;
145  }
146 
147  return runningSum;
148 }
149 
150 
152 {
153  // Reference:
154  // http://mathworld.wolfram.com/StatisticalMedian.html
155  // The statistical median is the value of the distribution variable
156  // where the cumulative distribution = 0.5.
157 
158  scalar median = 0.0;
159 
160  scalar runningSum = 0.0;
161 
162  List<Pair<scalar>> normDist(normalised());
163 
164  if (normDist.size())
165  {
166  if (normDist.size() == 1)
167  {
168  median = normDist[0].first();
169  }
170  else if
171  (
172  normDist.size() > 1
173  && normDist[0].second()*binWidth_ > 0.5
174  )
175  {
176  scalar xk = normDist[1].first();
177  scalar xkm1 = normDist[0].first();
178  scalar Sk =
179  (normDist[0].second() + normDist[1].second())*binWidth_;
180  scalar Skm1 = normDist[0].second()*binWidth_;
181 
182  median = (0.5 - Skm1)*(xk - xkm1)/(Sk - Skm1) + xkm1;
183  }
184  else
185  {
186  label lastNonZeroIndex = 0;
187 
188  forAll(normDist,nD)
189  {
190  if (runningSum + (normDist[nD].second()*binWidth_) > 0.5)
191  {
192  scalar xk = normDist[nD].first();
193  scalar xkm1 = normDist[lastNonZeroIndex].first();
194  scalar Sk = runningSum + (normDist[nD].second()*binWidth_);
195  scalar Skm1 = runningSum;
196 
197  median = (0.5 - Skm1)*(xk - xkm1)/(Sk - Skm1) + xkm1;
198 
199  break;
200  }
201  else if (normDist[nD].second() > 0.0)
202  {
203  runningSum += normDist[nD].second()*binWidth_;
204 
205  lastNonZeroIndex = nD;
206  }
207  }
208  }
209  }
210 
211  return median;
212 }
213 
214 
215 void Foam::distribution::add(const scalar valueToAdd)
216 {
217  iterator iter(this->begin());
218 
219  label n = label(valueToAdd/binWidth_) - label(neg(valueToAdd/binWidth_));
220 
221  iter = find(n);
222 
223  if (iter == this->end())
224  {
225  this->insert(n,1);
226  }
227  else
228  {
229  (*this)[n]++;
230  }
231 
232  if ((*this)[n] < 0)
233  {
235  << "Accumulated distribution value has become negative: "
236  << "bin = " << (0.5 + scalar(n)) * binWidth_
237  << ", value = " << (*this)[n]
238  << ". This is most likely to be because too many samples "
239  << "have been added to a bin and the label has 'rolled round'"
240  << abort(FatalError);
241  }
242 }
243 
244 
245 void Foam::distribution::add(const label valueToAdd)
246 {
247  add(scalar(valueToAdd));
248 }
249 
250 
252 {
253  List<label> keys = sortedToc();
254 
255  if (keys.size() > 2)
256  {
257  for (label k = keys[1]; k < keys.last(); k++)
258  {
259  // Insert 0, if the entry does not already exist
260  this->insert(k,0);
261  }
262  }
263 }
264 
265 
267 {
268  scalar totEnt = approxTotalEntries();
269 
270  insertMissingKeys();
271 
272  List<label> keys = sortedToc();
273  List<Pair<scalar>> normDist(size());
274 
275  forAll(keys,k)
276  {
277  label key = keys[k];
278 
279  normDist[k].first() = (0.5 + scalar(key))*binWidth_;
280 
281  normDist[k].second() = scalar((*this)[key])/totEnt/binWidth_;
282  }
283 
284  if (debug)
285  {
286  Info<< "totEnt: " << totEnt << endl;
287  }
288 
289  return normDist;
290 }
291 
292 
294 {
295  return normalisedShifted(mean());
296 }
297 
298 
300 (
301  scalar shiftValue
302 )
303 {
304  List<Pair<scalar>> oldDist(normalised());
305 
306  List<Pair<scalar>> newDist(oldDist.size());
307 
308  forAll(oldDist,u)
309  {
310  oldDist[u].first() -= shiftValue;
311  }
312 
313  scalar lowestOldBin = oldDist[0].first()/binWidth_ - 0.5;
314 
315  label lowestNewKey = label
316  (
317  lowestOldBin + 0.5*sign(lowestOldBin)
318  );
319 
320  scalar interpolationStartDirection =
321  sign(scalar(lowestNewKey) - lowestOldBin);
322 
323  label newKey = lowestNewKey;
324 
325  if (debug)
326  {
327  Info<< shiftValue
328  << nl << lowestOldBin
329  << nl << lowestNewKey
330  << nl << interpolationStartDirection
331  << endl;
332 
333  scalar checkNormalisation = 0;
334 
335  forAll(oldDist, oD)
336  {
337  checkNormalisation += oldDist[oD].second()*binWidth_;
338  }
339 
340  Info<< "Initial normalisation = " << checkNormalisation << endl;
341  }
342 
343  forAll(oldDist,u)
344  {
345  newDist[u].first() = (0.5 + scalar(newKey)) * binWidth_;
346 
347  if (interpolationStartDirection < 0)
348  {
349  if (u == 0)
350  {
351  newDist[u].second() =
352  (0.5 + scalar(newKey))*oldDist[u].second()
353  - oldDist[u].second()
354  *(oldDist[u].first() - binWidth_)/ binWidth_;
355  }
356  else
357  {
358  newDist[u].second() =
359  (0.5 + scalar(newKey))
360  *(oldDist[u].second() - oldDist[u-1].second())
361  +
362  (
363  oldDist[u-1].second()*oldDist[u].first()
364  - oldDist[u].second()*oldDist[u-1].first()
365  )
366  /binWidth_;
367  }
368  }
369  else
370  {
371  if (u == oldDist.size() - 1)
372  {
373  newDist[u].second() =
374  (0.5 + scalar(newKey))*-oldDist[u].second()
375  + oldDist[u].second()*(oldDist[u].first() + binWidth_)
376  /binWidth_;
377  }
378  else
379  {
380  newDist[u].second() =
381  (0.5 + scalar(newKey))
382  *(oldDist[u+1].second() - oldDist[u].second())
383  +
384  (
385  oldDist[u].second()*oldDist[u+1].first()
386  - oldDist[u+1].second()*oldDist[u].first()
387  )
388  /binWidth_;
389  }
390  }
391 
392  newKey++;
393  }
394 
395  if (debug)
396  {
397  scalar checkNormalisation = 0;
398 
399  forAll(newDist, nD)
400  {
401  checkNormalisation += newDist[nD].second()*binWidth_;
402  }
403 
404  Info<< "Shifted normalisation = " << checkNormalisation << endl;
405  }
406 
407  return newDist;
408 }
409 
410 
412 {
413  insertMissingKeys();
414 
415  List<label> keys = sortedToc();
416  List<Pair<scalar>> rawDist(size());
417 
418  forAll(keys,k)
419  {
420  label key = keys[k];
421 
422  rawDist[k].first() = (0.5 + scalar(key))*binWidth_;
423 
424  rawDist[k].second() = scalar((*this)[key]);
425  }
426 
427  return rawDist;
428 }
429 
430 
431 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
432 
434 {
435  if (this == &rhs)
436  {
437  return; // Self-assignment is a no-op
438  }
439 
441 
442  binWidth_ = rhs.binWidth();
443 }
444 
445 
446 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
447 
449 {
450  os << d.binWidth_
451  << static_cast<const Map<label>&>(d);
452 
454  return os;
455 }
456 
457 
458 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
insert
srcOptions insert("case", fileName(rootDirSource/caseDirSource))
Foam::distribution::binWidth
scalar binWidth() const
Definition: distributionI.H:30
Foam::distribution::mean
scalar mean() const
Definition: distribution.C:128
Foam::distribution
Accumulating histogram of values. Specified bin resolution automatic generation of bins.
Definition: distribution.H:61
Foam::distribution::add
void add(const scalar valueToAdd)
Add a value to the appropriate bin of the distribution.
Definition: distribution.C:215
stdFoam::begin
constexpr auto begin(C &c) -> decltype(c.begin())
Return iterator to the beginning of the container c.
Definition: stdFoam.H:97
Foam::fileName
A class for handling file names.
Definition: fileName.H:73
Foam::distribution::operator=
void operator=(const distribution &)
Definition: distribution.C:433
Foam::Map::operator=
void operator=(const this_type &rhs)
Copy assignment.
Definition: Map.H:117
Foam::glTF::key
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
Definition: foamGltfBase.H:108
Foam::Map
A HashTable to objects of type <T> with a label key.
Definition: lumpedPointController.H:69
Foam::distribution::normalisedMinusMean
List< Pair< scalar > > normalisedMinusMean()
Definition: distribution.C:293
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::sign
dimensionedScalar sign(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:166
Foam::Map< label >::iterator
typename parent_type::iterator iterator
Definition: Map.H:69
Foam::distribution::median
scalar median()
Definition: distribution.C:151
Foam::distribution::totalEntries
label totalEntries() const
Definition: distribution.C:87
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
OFstream.H
Foam::operator<<
Ostream & operator<<(Ostream &, const boundaryPatch &p)
Write boundaryPatch as dictionary entries (without surrounding braces)
Definition: boundaryPatch.C:83
n
label n
Definition: TABSMDCalcMethod2.H:31
distribution.H
Foam::distribution::insertMissingKeys
void insertMissingKeys()
Definition: distribution.C:251
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::distribution::approxTotalEntries
scalar approxTotalEntries() const
Definition: distribution.C:115
Foam::distribution::raw
List< Pair< scalar > > raw()
Definition: distribution.C:411
Foam::IOstream::check
virtual bool check(const char *operation) const
Check IOstream status for given operation.
Definition: IOstream.C:58
Foam::FatalError
error FatalError
Foam::add
void add(FieldField< Field1, typename typeOfSum< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
Definition: FieldFieldFunctions.C:939
os
OBJstream os(runTime.globalPath()/outputName)
Foam::distribution::write
static void write(const fileName &file, const List< Pair< scalar >> &pairs)
Write to file.
Definition: distribution.C:42
stdFoam::end
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:121
Foam::ListOps::find
label find(const ListType &input, const UnaryPredicate &pred, const label start=0)
Find index of the first occurrence that satisfies the predicate.
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::distribution::~distribution
virtual ~distribution()
Destructor.
Definition: distribution.C:81
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Foam::OFstream
Output to file stream, using an OSstream.
Definition: OFstream.H:53
Foam::distribution::distribution
distribution()
Construct null.
Definition: distribution.C:58
Foam::normalised
VectorSpace< Form, Cmpt, Ncmpts > normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:487
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::Pair< scalar >
forAllConstIters
forAllConstIters(mixture.phases(), phase)
Definition: pEqn.H:28
Foam::List
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: BitOps.H:63
k
label k
Boltzmann constant.
Definition: LISASMDCalcMethod2.H:41
FUNCTION_NAME
#define FUNCTION_NAME
Definition: messageStream.H:295
Foam::distribution::normalised
List< Pair< scalar > > normalised()
Definition: distribution.C:266
Foam::Ostream
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:56
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
Foam::neg
dimensionedScalar neg(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:199
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
Foam::distribution::normalisedShifted
List< Pair< scalar > > normalisedShifted(scalar shiftValue)
Definition: distribution.C:300