FreeFOAM The Cross-Platform CFD Toolkit
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) 2008-2010 OpenCFD Ltd.
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 
28 namespace Foam
29 {
30 
31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 
34 :
35  Map<label>(),
36  binWidth_(1)
37 {}
38 
39 
40 distribution::distribution(const scalar binWidth)
41 :
42  Map<label>(),
43  binWidth_(binWidth)
44 {}
45 
46 
48 :
49  Map<label>(static_cast< Map<label> >(d)),
50  binWidth_(d.binWidth())
51 {}
52 
53 
54 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
55 
57 {}
58 
59 
60 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
61 
63 {
64  label sumOfEntries = 0;
65 
66  forAllConstIter(Map<label>, *this, iter)
67  {
68  sumOfEntries += iter();
69 
70  if (sumOfEntries < 0)
71  {
72  WarningIn("label distribution::totalEntries()")
73  << "Accumulated distribution values total has become negative: "
74  << "sumOfEntries = " << sumOfEntries
75  << ". This is most likely to be because too many samples "
76  << "have been added to the bins and the label has 'rolled "
77  << "round'. Try distribution::approxTotalEntries which "
78  << "returns a scalar." << endl;
79 
80  sumOfEntries = -1;
81 
82  break;
83  }
84  }
85 
86  return sumOfEntries;
87 }
88 
89 
91 {
92  scalar sumOfEntries = 0;
93 
94  forAllConstIter(Map<label>, *this, iter)
95  {
96  sumOfEntries += scalar(iter());
97  }
98 
99  return sumOfEntries;
100 }
101 
102 
103 scalar distribution::mean() const
104 {
105  scalar runningSum = 0;
106 
107  scalar totEnt = approxTotalEntries();
108 
109  List<label> keys = toc();
110 
111  forAll(keys,k)
112  {
113  label key = keys[k];
114 
115  runningSum +=
116  (0.5 + scalar(key))
117  *binWidth_
118  *scalar((*this)[key])
119  /totEnt;
120  }
121 
122  return runningSum;
123 }
124 
125 
127 {
128  // From:
129  // http://mathworld.wolfram.com/StatisticalMedian.html
130  // The statistical median is the value of the distribution variable
131  // where the cumulative distribution = 0.5.
132 
133  scalar median = 0.0;
134 
135  scalar runningSum = 0.0;
136 
137  List<Pair<scalar> > normDist(normalised());
138 
139  if (normDist.size())
140  {
141  if (normDist.size() == 1)
142  {
143  median = normDist[0].first();
144  }
145  else if
146  (
147  normDist.size() > 1
148  && normDist[0].second()*binWidth_ > 0.5
149  )
150  {
151  scalar xk = normDist[1].first();
152  scalar xkm1 = normDist[0].first();
153  scalar Sk =
154  (normDist[0].second() + normDist[1].second())*binWidth_;
155  scalar Skm1 = normDist[0].second()*binWidth_;
156 
157  median = (0.5 - Skm1)*(xk - xkm1)/(Sk - Skm1) + xkm1;
158  }
159  else
160  {
161  label lastNonZeroIndex = 0;
162 
163  forAll(normDist,nD)
164  {
165  if (runningSum + (normDist[nD].second()*binWidth_) > 0.5)
166  {
167  scalar xk = normDist[nD].first();
168  scalar xkm1 = normDist[lastNonZeroIndex].first();
169  scalar Sk = runningSum + (normDist[nD].second()*binWidth_);
170  scalar Skm1 = runningSum;
171 
172  median = (0.5 - Skm1)*(xk - xkm1)/(Sk - Skm1) + xkm1;
173 
174  break;
175  }
176  else if (normDist[nD].second() > 0.0)
177  {
178  runningSum += normDist[nD].second()*binWidth_;
179 
180  lastNonZeroIndex = nD;
181  }
182  }
183  }
184  }
185 
186  return median;
187 }
188 
189 
190 void distribution::add(const scalar valueToAdd)
191 {
192  iterator iter(this->begin());
193 
194  label n = label(valueToAdd/binWidth_) - label(neg(valueToAdd/binWidth_));
195 
196  iter = find(n);
197 
198  if (iter == this->end())
199  {
200  this->insert(n,1);
201  }
202  else
203  {
204  (*this)[n]++;
205  }
206 
207  if ((*this)[n] < 0)
208  {
209  FatalErrorIn("distribution::add(const scalar valueToAdd)")
210  << "Accumulated distribution value has become negative: "
211  << "bin = " << (0.5 + scalar(n)) * binWidth_
212  << ", value = " << (*this)[n]
213  << ". This is most likely to be because too many samples "
214  << "have been added to a bin and the label has 'rolled round'"
215  << abort(FatalError);
216  }
217 }
218 
219 
220 void distribution::add(const label valueToAdd)
221 {
222  add(scalar(valueToAdd));
223 }
224 
225 
227 {
228  iterator iter(this->begin());
229 
230  List<label> keys = toc();
231 
232  sort(keys);
233 
234  if (keys.size())
235  {
236  for (label k = keys[0]; k < keys[keys.size()-1]; k++)
237  {
238  iter = find(k);
239 
240  if (iter == this->end())
241  {
242  this->insert(k,0);
243  }
244  }
245  }
246 }
247 
248 
250 {
251  scalar totEnt = approxTotalEntries();
252 
254 
255  List<label> keys = toc();
256 
257  sort(keys);
258 
259  List<Pair<scalar> > normDist(size());
260 
261  forAll(keys,k)
262  {
263  label key = keys[k];
264 
265  normDist[k].first() = (0.5 + scalar(key))*binWidth_;
266 
267  normDist[k].second() = scalar((*this)[key])/totEnt/binWidth_;
268  }
269 
270  return normDist;
271 }
272 
273 
275 {
276  return normalisedShifted(mean());
277 }
278 
279 
281 {
282  List<Pair<scalar> > oldDist(normalised());
283 
284  List<Pair<scalar> > newDist(oldDist.size());
285 
286  forAll(oldDist,u)
287  {
288  oldDist[u].first() -= shiftValue;
289  }
290 
291  scalar lowestOldBin = oldDist[0].first()/binWidth_ - 0.5;
292 
293  label lowestNewKey = label
294  (
295  lowestOldBin + 0.5*sign(lowestOldBin)
296  );
297 
298  scalar interpolationStartDirection =
299  sign(scalar(lowestNewKey) - lowestOldBin);
300 
301  label newKey = lowestNewKey;
302 
303 // Info << shiftValue
304 // << nl << lowestOldBin
305 // << nl << lowestNewKey
306 // << nl << interpolationStartDirection
307 // << endl;
308 
309 // scalar checkNormalisation = 0;
310 
311 // forAll (oldDist, oD)
312 // {
313 // checkNormalisation += oldDist[oD].second()*binWidth_;
314 // }
315 
316 // Info << "Initial normalisation = " << checkNormalisation << endl;
317 
318  forAll(oldDist,u)
319  {
320  newDist[u].first() = (0.5 + scalar(newKey)) * binWidth_;
321 
322  if (interpolationStartDirection < 0)
323  {
324  if (u == 0)
325  {
326  newDist[u].second() =
327  (0.5 + scalar(newKey))*oldDist[u].second()
328  - oldDist[u].second()
329  *(oldDist[u].first() - binWidth_)/ binWidth_;
330  }
331  else
332  {
333  newDist[u].second() =
334  (0.5 + scalar(newKey))
335  *(oldDist[u].second() - oldDist[u-1].second())
336  +
337  (
338  oldDist[u-1].second()*oldDist[u].first()
339  - oldDist[u].second()*oldDist[u-1].first()
340  )
341  /binWidth_;
342  }
343  }
344  else
345  {
346  if (u == oldDist.size() - 1)
347  {
348  newDist[u].second() =
349  (0.5 + scalar(newKey))*-oldDist[u].second()
350  + oldDist[u].second()*(oldDist[u].first() + binWidth_)
351  /binWidth_;
352  }
353  else
354  {
355  newDist[u].second() =
356  (0.5 + scalar(newKey))
357  *(oldDist[u+1].second() - oldDist[u].second())
358  +
359  (
360  oldDist[u].second()*oldDist[u+1].first()
361  - oldDist[u+1].second()*oldDist[u].first()
362  )
363  /binWidth_;
364  }
365  }
366 
367  newKey++;
368  }
369 
370 // checkNormalisation = 0;
371 
372 // forAll (newDist, nD)
373 // {
374 // checkNormalisation += newDist[nD].second()*binWidth_;
375 // }
376 
377 // Info << "Shifted normalisation = " << checkNormalisation << endl;
378 
379  return newDist;
380 }
381 
382 
384 {
386 
387  List<label> keys = toc();
388 
389  sort(keys);
390 
391  List<Pair<scalar> > rawDist(size());
392 
393  forAll(keys,k)
394  {
395  label key = keys[k];
396 
397  rawDist[k].first() = (0.5 + scalar(key))*binWidth_;
398 
399  rawDist[k].second() = scalar((*this)[key]);
400  }
401 
402  return rawDist;
403 }
404 
405 
406 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
407 
409 {
410  // Check for assignment to self
411  if (this == &rhs)
412  {
413  FatalErrorIn("distribution::operator=(const distribution&)")
414  << "Attempted assignment to self"
415  << abort(FatalError);
416  }
417 
419 
420  binWidth_ = rhs.binWidth();
421 }
422 
423 
424 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
425 
427 {
428  os << d.binWidth_
429  << static_cast<const Map<label>&>(d);
430 
431  // Check state of Ostream
432  os.check
433  (
434  "Ostream& operator<<(Ostream&, "
435  "const distribution&)"
436  );
437 
438  return os;
439 }
440 
441 
442 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
443 
444 } // End namespace Foam
445 
446 // ************************ vim: set sw=4 sts=4 et: ************************ //