src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C

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