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#include "ListOps.H"
32
33// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
34
35template<class Type>
37:
38 List<List<scalar>>(pTraits<Type>::nComponents),
39 binWidth_(pTraits<Type>::one),
40 listStarts_(pTraits<Type>::nComponents, 0)
41{}
42
43
44template<class Type>
46:
47 List<List<scalar>>(pTraits<Type>::nComponents),
48 binWidth_(binWidth),
49 listStarts_(pTraits<Type>::nComponents, 0)
50{}
51
52
53template<class Type>
55:
56 List<List<scalar>>(static_cast<const List<List<scalar>>&>(d)),
57 binWidth_(d.binWidth()),
58 listStarts_(d.listStarts())
59{}
60
61
62// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
63
64template<class Type>
66{}
67
68
69// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
70
71template<class Type>
73{
74 const List<scalar>& cmptDistribution = (*this)[cmpt];
75
76 scalar sumOfWeights = 0.0;
77
78 forAll(cmptDistribution, i)
79 {
80 sumOfWeights += cmptDistribution[i];
81 }
82
83 return sumOfWeights;
84}
85
86
87template<class Type>
89{
90 return identity((*this)[cmpt].size(), listStarts_[cmpt]);
91}
92
93
94template<class Type>
96(
97 direction cmpt,
98 label n
99)
100{
101 List<scalar>& cmptDistribution = (*this)[cmpt];
102
103 if (cmptDistribution.empty())
104 {
105 // Initialise this list with this value
106
107 cmptDistribution.setSize(2, 0.0);
108
109 listStarts_[cmpt] = n;
110
111 return 0;
112 }
113
114 label listIndex = -1;
115
116 label& listStart = listStarts_[cmpt];
117
118 label testIndex = n - listStart;
119
120 if (testIndex < 0)
121 {
122 // Underflow of this List, storage increase and remapping
123 // required
124
125 List<scalar> newCmptDistribution(2*cmptDistribution.size(), Zero);
126
127 label sOld = cmptDistribution.size();
128
129 forAll(cmptDistribution, i)
130 {
131 newCmptDistribution[i + sOld] = cmptDistribution[i];
132 }
133
134 cmptDistribution = newCmptDistribution;
135
136 listStart -= sOld;
137
138 // Recursively call this function in case another remap is required.
139 listIndex = index(cmpt, n);
140 }
141 else if (testIndex > cmptDistribution.size() - 1)
142 {
143 // Overflow of this List, storage increase required
144
145 cmptDistribution.setSize(2*cmptDistribution.size(), 0.0);
146
147 // Recursively call this function in case another storage
148 // alteration is required.
149 listIndex = index(cmpt, n);
150 }
151 else
152 {
153 listIndex = n - listStart;
154 }
155
156 return listIndex;
157}
158
159
160template<class Type>
162(
163 direction cmpt
164) const
165{
166 const List<scalar>& cmptDistribution = (*this)[cmpt];
167
168 // limits.first(): lower bound, i.e. the first non-zero entry
169 // limits.second(): upper bound, i.e. the last non-zero entry
170 Pair<label> limits(-1, -1);
171
172 forAll(cmptDistribution, i)
173 {
174 if (cmptDistribution[i] > 0.0)
175 {
176 if (limits.first() == -1)
177 {
178 limits.first() = i;
179 limits.second() = i;
180 }
181 else
182 {
183 limits.second() = i;
184 }
185 }
186 }
187
188 return limits;
189}
190
191
192template<class Type>
194{
195 Type meanValue(Zero);
196
197 for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
198 {
199 const List<scalar>& cmptDistribution = (*this)[cmpt];
200
201 scalar totalCmptWeight = totalWeight(cmpt);
202
203 List<label> theKeys = keys(cmpt);
204
205 forAll(theKeys, k)
206 {
207 label key = theKeys[k];
208
209 setComponent(meanValue, cmpt) +=
210 (0.5 + scalar(key))
211 *component(binWidth_, cmpt)
212 *cmptDistribution[k]
213 /totalCmptWeight;
214 }
215 }
216
217 return meanValue;
218}
219
220
221template<class Type>
223{
224 Type medianValue(Zero);
225
226 List<List<Pair<scalar>>> normDistribution = normalised();
227
228 for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
229 {
230 List<Pair<scalar>>& normDist = normDistribution[cmpt];
231
232 if (normDist.size())
233 {
234 if (normDist.size() == 1)
235 {
236 setComponent(medianValue, cmpt) = normDist[0].first();
237 }
238 else if
239 (
240 normDist.size() > 1
241 && normDist[0].second()*component(binWidth_, cmpt) > 0.5
242 )
243 {
244 scalar xk =
245 normDist[0].first()
246 + 0.5*component(binWidth_, cmpt);
247
248 scalar xkm1 =
249 normDist[0].first()
250 - 0.5*component(binWidth_, cmpt);
251
252 scalar Sk = (normDist[0].second())*component(binWidth_, cmpt);
253
254 setComponent(medianValue, cmpt) = 0.5*(xk - xkm1)/(Sk) + xkm1;
255 }
256 else
257 {
258 label previousNonZeroIndex = 0;
259
260 scalar cumulative = 0.0;
261
262 forAll(normDist, nD)
263 {
264 if
265 (
266 cumulative
267 + (normDist[nD].second()*component(binWidth_, cmpt))
268 > 0.5
269 )
270 {
271 scalar xk =
272 normDist[nD].first()
273 + 0.5*component(binWidth_, cmpt);
274
275 scalar xkm1 =
276 normDist[previousNonZeroIndex].first()
277 + 0.5*component(binWidth_, cmpt);
278
279 scalar Sk =
280 cumulative
281 + (normDist[nD].second()*component(binWidth_, cmpt));
282
283 scalar Skm1 = cumulative;
284
285 setComponent(medianValue, cmpt) =
286 (0.5 - Skm1)*(xk - xkm1)/(Sk - Skm1) + xkm1;
287
288 break;
289 }
290 else if (mag(normDist[nD].second()) > VSMALL)
291 {
292 cumulative +=
293 normDist[nD].second()*component(binWidth_, cmpt);
294
295 previousNonZeroIndex = nD;
296 }
297 }
298 }
299 }
300
301 }
302
303 return medianValue;
304}
305
306
307template<class Type>
309(
310 const Type& valueToAdd,
311 const Type& weight
312)
313{
314 for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
315 {
316 List<scalar>& cmptDistribution = (*this)[cmpt];
317
318 label n =
319 label(component(valueToAdd, cmpt)/component(binWidth_, cmpt))
320 - label(neg(component(valueToAdd, cmpt)/component(binWidth_, cmpt)));
321
322 label listIndex = index(cmpt, n);
323
324 cmptDistribution[listIndex] += component(weight, cmpt);
325 }
326}
327
328
329template<class Type>
332{
334
335 for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
336 {
337 const List<scalar>& cmptDistribution = (*this)[cmpt];
338
339 if (cmptDistribution.empty())
340 {
341 continue;
342 }
343
344 scalar totalCmptWeight = totalWeight(cmpt);
345
346 List<label> cmptKeys = keys(cmpt);
347
348 List<Pair<scalar>>& normDist = normDistribution[cmpt];
349
350 Pair<label> limits = validLimits(cmpt);
351
352 normDist.setSize(limits.second() - limits.first() + 1);
353
354 for
355 (
356 label k = limits.first(), i = 0;
357 k <= limits.second();
358 k++, i++
359 )
360 {
361 label key = cmptKeys[k];
362
363 normDist[i].first() =
364 (0.5 + scalar(key))*component(binWidth_, cmpt);
365
366 normDist[i].second() =
367 cmptDistribution[k]
368 /totalCmptWeight
369 /component(binWidth_, cmpt);
370 }
371 }
372
373 return normDistribution;
374}
375
376
377template<class Type>
380{
382
383 for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
384 {
385 const List<scalar>& cmptDistribution = (*this)[cmpt];
386
387 if (cmptDistribution.empty())
388 {
389 continue;
390 }
391
392 List<label> cmptKeys = keys(cmpt);
393
394 List<Pair<scalar>>& rawDist = rawDistribution[cmpt];
395
396 Pair<label> limits = validLimits(cmpt);
397
398 rawDist.setSize(limits.second() - limits.first() + 1);
399
400 for
401 (
402 label k = limits.first(), i = 0;
403 k <= limits.second();
404 k++, i++
405 )
406 {
407 label key = cmptKeys[k];
408
409 rawDist[i].first() = (0.5 + scalar(key))*component(binWidth_, cmpt);
410
411 rawDist[i].second() = cmptDistribution[k];
412 }
413 }
414
415 return rawDistribution;
416}
417
418
419template<class Type>
422{
423 List<List<Pair<scalar>>> normalisedDistribution = normalised();
424
425 List<List<Pair<scalar>>> cumulativeNormalisedDistribution =
426 normalisedDistribution;
427
428 for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
429 {
430 const List<Pair<scalar>>& normalisedCmpt =
431 normalisedDistribution[cmpt];
432
433 List<Pair<scalar>>& cumNormalisedCmpt =
434 cumulativeNormalisedDistribution[cmpt];
435
436 scalar sum = 0.0;
437
438 forAll(normalisedCmpt, i)
439 {
440 cumNormalisedCmpt[i].first() =
441 normalisedCmpt[i].first()
442 + 0.5*component(binWidth_, cmpt);
443
444 cumNormalisedCmpt[i].second() =
445 normalisedCmpt[i].second()*component(binWidth_, cmpt) + sum;
446
447 sum = cumNormalisedCmpt[i].second();
448 }
449 }
450
451 return cumulativeNormalisedDistribution;
452}
453
454
455template<class Type>
458{
459 List<List<Pair<scalar>>> rawDistribution = raw();
460
461 List<List<Pair<scalar>>> cumulativeRawDistribution = rawDistribution;
462
463 for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
464 {
465 const List<Pair<scalar>>& rawCmpt = rawDistribution[cmpt];
466
467 List<Pair<scalar>>& cumRawCmpt = cumulativeRawDistribution[cmpt];
468
469 scalar sum = 0.0;
470
471 forAll(rawCmpt, i)
472 {
473 cumRawCmpt[i].first() =
474 rawCmpt[i].first()
475 + 0.5*component(binWidth_, cmpt);
476
477 cumRawCmpt[i].second() = rawCmpt[i].second() + sum;
478
479 sum = cumRawCmpt[i].second();
480 }
481 }
482
483 return cumulativeRawDistribution;
484}
485
486
487template<class Type>
489{
490 for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
491 {
492 (*this)[cmpt].clear();
493
494 listStarts_[cmpt] = 0;
495 }
496}
497
498
499template<class Type>
500void Foam::Distribution<Type>::write(const fileName& filePrefix) const
501{
502 List<List<Pair<scalar>>> rawDistribution = raw();
503
504 List<List<Pair<scalar>>> normDistribution = normalised();
505
506 for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
507 {
508 const List<Pair<scalar>>& rawPairs = rawDistribution[cmpt];
509
510 const List<Pair<scalar>>& normPairs = normDistribution[cmpt];
511
512 OFstream os(filePrefix + '_' + pTraits<Type>::componentNames[cmpt]);
513
514 os << "# key normalised raw" << endl;
515
516 forAll(normPairs, i)
517 {
518 os << normPairs[i].first()
519 << ' ' << normPairs[i].second()
520 << ' ' << rawPairs[i].second()
521 << nl;
522 }
523 }
524
525 List<List<Pair<scalar>>> rawCumDist = cumulativeRaw();
526
527 List<List<Pair<scalar>>> normCumDist = cumulativeNormalised();
528
529 for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
530 {
531 const List<Pair<scalar>>& rawPairs = rawCumDist[cmpt];
532
533 const List<Pair<scalar>>& normPairs = normCumDist[cmpt];
534
536 (
537 filePrefix + "_cumulative_" + pTraits<Type>::componentNames[cmpt]
538 );
539
540 os << "# key normalised raw" << endl;
541
542 forAll(normPairs, i)
543 {
544 os << normPairs[i].first()
545 << ' ' << normPairs[i].second()
546 << ' ' << rawPairs[i].second()
547 << nl;
548 }
549 }
550}
551
552
553// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * //
554
555template<class Type>
557(
558 const Distribution<Type>& rhs
559)
560{
561 if (this == &rhs)
562 {
563 return; // Self-assignment is a no-op
564 }
565
566 List<List<scalar>>::operator=(rhs);
567
568 binWidth_ = rhs.binWidth();
569
570 listStarts_ = rhs.listStarts();
571}
572
573
574// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
575
576template<class Type>
577Foam::Istream& Foam::operator>>
578(
579 Istream& is,
581)
582{
583 is >> static_cast<List<List<scalar>>&>(d)
584 >> d.binWidth_
585 >> d.listStarts_;
586
587 is.check(FUNCTION_NAME);
588 return is;
589}
590
591
592template<class Type>
593Foam::Ostream& Foam::operator<<
594(
595 Ostream& os,
596 const Distribution<Type>& d
597)
598{
599 os << static_cast<const List<List<scalar>>&>(d)
600 << d.binWidth_ << token::SPACE
601 << d.listStarts_;
602
604 return os;
605}
606
607
608// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
609
610template<class Type>
611Foam::Distribution<Type> Foam::operator+
612(
613 const Distribution<Type>& d1,
614 const Distribution<Type>& d2
615)
616{
617 // The coarsest binWidth is the sensible choice
618 Distribution<Type> d(max(d1.binWidth(), d2.binWidth()));
619
620 List<List<List<Pair<scalar>>>> rawDists(2);
621
622 rawDists[0] = d1.raw();
623 rawDists[1] = d2.raw();
624
625 forAll(rawDists, rDI)
626 {
627 for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
628 {
629 List<scalar>& cmptDistribution = d[cmpt];
630
631 const List<Pair<scalar>>& cmptRaw = rawDists[rDI][cmpt];
632
633 forAll(cmptRaw, rI)
634 {
635 scalar valueToAdd = cmptRaw[rI].first();
636 scalar cmptWeight = cmptRaw[rI].second();
637
638 label n =
639 label
640 (
641 component(valueToAdd, cmpt)
642 /component(d.binWidth(), cmpt)
643 )
644 - label
645 (
646 neg(component(valueToAdd, cmpt)
647 /component(d.binWidth(), cmpt))
648 );
649
650 label listIndex = d.index(cmpt, n);
651
652 cmptDistribution[listIndex] += cmptWeight;
653 }
654 }
655 }
656
657 return Distribution<Type>(d);
658}
659
660
661// ************************************************************************* //
label k
Various functions to operate on Lists.
label n
Accumulating histogram of component values. Specified bin resolution, automatic generation of bins.
Definition: Distribution.H:67
List< List< Pair< scalar > > > cumulativeNormalised() const
Return the cumulative normalised distribution and.
Definition: Distribution.C:421
Type mean() const
Definition: Distribution.C:193
~Distribution()
Destructor.
Definition: Distribution.C:65
Type median() const
Definition: Distribution.C:222
Distribution()
Construct null.
Definition: Distribution.C:36
scalar totalWeight(direction cmpt) const
Sum the total weight added to the component in the.
Definition: Distribution.C:72
void clear()
Resets the Distribution by clearing the stored lists.
Definition: Distribution.C:488
Pair< label > validLimits(direction cmpt) const
Returns the indices of the first and last non-zero entries.
Definition: Distribution.C:162
const Type & binWidth() const
Return the bin width.
Definition: DistributionI.H:32
List< List< Pair< scalar > > > cumulativeRaw() const
Return the cumulative total bin weights and integration.
Definition: Distribution.C:457
List< List< Pair< scalar > > > raw() const
Return the distribution of the total bin weights.
Definition: Distribution.C:379
List< List< Pair< scalar > > > normalised() const
Return the normalised distribution (probability density)
Definition: Distribution.C:331
T & first() noexcept
The first element of the list, position [0].
Definition: FixedListI.H:207
const_iterator_pair< const_key_iterator, this_type > keys() const
A const iterator begin/end pair for iterating over keys.
Definition: HashTable.H:924
virtual bool check(const char *operation) const
Check IOstream status for given operation.
Definition: IOstream.C:58
An Istream is an abstract base class for all input systems (streams, files, token lists etc)....
Definition: Istream.H:64
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
void setSize(const label n)
Alias for resize()
Definition: List.H:218
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
const T & second() const noexcept
Return second element, which is also the last element.
Definition: PairI.H:120
T & first()
Return the first element of the list.
Definition: UListI.H:202
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
Definition: UListI.H:427
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
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
A class representing the concept of 1 (one) that can be used to avoid manipulating objects known to b...
Definition: one.H:62
A traits class, which is primarily used for primitives.
Definition: pTraits.H:59
@ SPACE
Space [isspace].
Definition: token.H:125
OBJstream os(runTime.globalPath()/outputName)
#define FUNCTION_NAME
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:680
void component(FieldField< Field, typename FieldField< Field, Type >::cmptType > &sf, const FieldField< Field, Type > &f, const direction d)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
label & setComponent(label &val, const direction) noexcept
Non-const access to integer-type (has no components)
Definition: label.H:123
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
uint8_t direction
Definition: direction.H:56
dimensionedScalar neg(const dimensionedScalar &ds)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333