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