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-------------------------------------------------------------------------------
11License
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
34namespace Foam
35{
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
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
128Foam::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
215void 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
245void 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// ************************************************************************* //
label k
label n
Forward iterator with non-const access.
Definition: HashTable.H:720
virtual bool check(const char *operation) const
Check IOstream status for given operation.
Definition: IOstream.C:58
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
void operator=(const this_type &rhs)
Copy assignment.
Definition: Map.H:117
Output to file stream, using an OSstream.
Definition: OFstream.H:57
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: Pair.H:69
T & first()
Return the first element of the list.
Definition: UListI.H:202
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T & last()
Return the last element of the list.
Definition: UListI.H:216
Accumulating histogram of values. Specified bin resolution automatic generation of bins.
Definition: distribution.H:64
void operator=(const distribution &)
Definition: distribution.C:433
List< Pair< scalar > > normalisedMinusMean()
Definition: distribution.C:293
List< Pair< scalar > > raw()
Definition: distribution.C:411
virtual ~distribution()
Destructor.
Definition: distribution.C:81
List< Pair< scalar > > normalised()
Definition: distribution.C:266
List< Pair< scalar > > normalisedShifted(scalar shiftValue)
Definition: distribution.C:300
label totalEntries() const
Definition: distribution.C:87
scalar approxTotalEntries() const
Definition: distribution.C:115
distribution()
Construct null.
Definition: distribution.C:58
scalar binWidth() const
Definition: distributionI.H:30
scalar mean() const
Definition: distribution.C:128
A class for handling file names.
Definition: fileName.H:76
virtual bool write()
Write the output fields.
Sums a given list of (at least two or more) fields and outputs the result into a new field,...
Definition: add.H:161
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
OBJstream os(runTime.globalPath()/outputName)
#define WarningInFunction
Report a warning using Foam::Warning.
#define FUNCTION_NAME
Namespace for OpenFOAM.
dimensionedScalar sign(const dimensionedScalar &ds)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:680
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & operator<<(Ostream &, const boundaryPatch &p)
Write boundaryPatch as dictionary entries (without surrounding braces)
Definition: boundaryPatch.C:83
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
void add(FieldField< Field1, typename typeOfSum< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
errorManip< error > abort(error &err)
Definition: errorManip.H:144
dimensionedScalar neg(const dimensionedScalar &ds)
error FatalError
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
srcOptions insert("case", fileName(rootDirSource/caseDirSource))
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278