src/OpenFOAM/containers/Lists/Distribution/Distribution.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2009-2011 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 #include "ListOps.H"
  29
  30 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  31
  32 template<class Type>
  33 Foam::Distribution<Type>::Distribution()
  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>
  42 Foam::Distribution<Type>::Distribution(const Type& binWidth)
  43 :
  44     List< List<scalar> >(pTraits<Type>::nComponents),
  45     binWidth_(binWidth),
  46     listStarts_(pTraits<Type>::nComponents, 0)
  47 {}
  48
  49
  50 template<class Type>
  51 Foam::Distribution<Type>::Distribution(const Distribution<Type>& d)
  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>
  62 Foam::Distribution<Type>::~Distribution()
  63 {}
  64
  65
  66 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
  67
  68 template<class Type>
  69 Foam::scalar Foam::Distribution<Type>::totalWeight(direction cmpt) const
  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>
  85 Foam::List<Foam::label> Foam::Distribution<Type>::keys(direction cmpt) const
  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>
  99 Foam::label Foam::Distribution<Type>::index
 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>
 165 Foam::Pair<Foam::label> Foam::Distribution<Type>::validLimits
 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>
 197 Type Foam::Distribution<Type>::mean() const
 198 {
 199     Type meanValue(pTraits<Type>::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>
 226 Type Foam::Distribution<Type>::median() const
 227 {
 228     Type medianValue(pTraits<Type>::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>
 312 void Foam::Distribution<Type>::add
 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>
 334 Foam::List< Foam::List< Foam::Pair<Foam::scalar> > >Foam::
 335 Distribution<Type>::normalised() const
 336 {
 337     List< List < Pair<scalar> > > normDistribution(pTraits<Type>::nComponents);
 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>
 382 Foam::List< Foam::List< Foam::Pair<Foam::scalar> > >Foam::
 383 Distribution<Type>::raw() const
 384 {
 385     List< List < Pair<scalar> > > rawDistribution(pTraits<Type>::nComponents);
 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>
 424 Foam::List< Foam::List< Foam::Pair<Foam::scalar> > >Foam::
 425 Distribution<Type>::cumulativeNormalised() const
 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>
 460 Foam::List< Foam::List< Foam::Pair<Foam::scalar> > >Foam::
 461 Distribution<Type>::cumulativeRaw() const
 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>
 492 void Foam::Distribution<Type>::clear()
 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
 531     List< List < Pair<scalar> > > normCumDist = cumulativeNormalised();
 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     {
 568         FatalErrorIn
 569         (
 570             "Foam::Distribution<Type>::operator="
 571             "(const Foam::Distribution<Type>&)"
 572         )   << "Attempted assignment to self"
 573             << abort(FatalError);
 574     }
 575
 576     List< List<scalar> >::operator=(rhs);
 577
 578     binWidth_ = rhs.binWidth();
 579
 580     listStarts_ = rhs.listStarts();
 581 }
 582
 583
 584 // * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
 585
 586 template <class Type>
 587 Foam::Istream& Foam::operator>>
 588 (
 589     Istream& is,
 590     Distribution<Type>& d
 591 )
 592 {
 593     is  >> static_cast<List< List<scalar> >&>(d)
 594         >> d.binWidth_
 595         >> d.listStarts_;
 596
 597     // Check state of Istream
 598     is.check("Istream& operator>>(Istream&, Distribution<Type>&)");
 599
 600     return is;
 601 }
 602
 603
 604 template<class Type>
 605 Foam::Ostream& Foam::operator<<
 606 (
 607     Ostream& os,
 608     const Distribution<Type>& d
 609 )
 610 {
 611     os  <<  static_cast<const List< List<scalar> >& >(d)
 612         << d.binWidth_ << token::SPACE
 613         << d.listStarts_;
 614
 615     // Check state of Ostream
 616     os.check("Ostream& operator<<(Ostream&, " "const Distribution&)");
 617
 618     return os;
 619 }
 620
 621
 622 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 623
 624 template <class Type>
 625 Foam::Distribution<Type> Foam::operator+
 626 (
 627     const Distribution<Type>& d1,
 628     const Distribution<Type>& d2
 629 )
 630 {
 631     // The coarsest binWidth is the sensible choice
 632     Distribution<Type> d(max(d1.binWidth(), d2.binWidth()));
 633
 634     List< List< List < Pair<scalar> > > > rawDists(2);
 635
 636     rawDists[0] = d1.raw();
 637     rawDists[1] = d2.raw();
 638
 639     forAll(rawDists, rDI)
 640     {
 641         for (direction cmpt = 0; cmpt < pTraits<Type>::nComponents; cmpt++)
 642         {
 643             List<scalar>& cmptDistribution = d[cmpt];
 644
 645             const List < Pair<scalar> >& cmptRaw = rawDists[rDI][cmpt];
 646
 647             forAll(cmptRaw, rI)
 648             {
 649                 scalar valueToAdd = cmptRaw[rI].first();
 650                 scalar cmptWeight = cmptRaw[rI].second();
 651
 652                 label n =
 653                 label
 654                 (
 655                     component(valueToAdd, cmpt)
 656                    /component(d.binWidth(), cmpt)
 657                 )
 658                 - label
 659                 (
 660                     neg(component(valueToAdd, cmpt)
 661                    /component(d.binWidth(), cmpt))
 662                 );
 663
 664                 label listIndex = d.index(cmpt, n);
 665
 666                 cmptDistribution[listIndex] += cmptWeight;
 667             }
 668         }
 669     }
 670
 671     return Distribution<Type>(d);
 672 }
 673
 674
 675 // ************************************************************************* //