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