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