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