Branch data     Line data    Source code

       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/Statistics/UniPDF.hpp
       4                 :            :   \copyright 2012-2015 J. Bakosi,
       5                 :            :              2016-2018 Los Alamos National Security, LLC.,
       6                 :            :              2019-2021 Triad National Security, LLC.,
       7                 :            :              2022-2024 J. Bakosi
       8                 :            :              All rights reserved. See the LICENSE file for details.
       9                 :            :   \brief     Univariate PDF estimator
      10                 :            :   \details   Univariate PDF estimator. This class can be used to estimate a
      11                 :            :     probability density function of (PDF) a scalar variable from an ensemble.
      12                 :            :     The implementation uses the standard container std::unordered_map, which is
      13                 :            :     a hash-based associative container with linear algorithmic complexity for
      14                 :            :     insertion of a new sample.
      15                 :            : */
      16                 :            : // *****************************************************************************
      17                 :            : #ifndef UniPDF_h
      18                 :            : #define UniPDF_h
      19                 :            : 
      20                 :            : #include <array>
      21                 :            : #include <unordered_map>
      22                 :            : #include <algorithm>
      23                 :            : #include <cfenv>
      24                 :            : 
      25                 :            : #include "Types.hpp"
      26                 :            : #include "Exception.hpp"
      27                 :            : #include "PUPUtil.hpp"
      28                 :            : 
      29                 :            : namespace tk {
      30                 :            : 
      31                 :            : //! Univariate PDF estimator
      32         [ +  - ]:      34745 : class UniPDF {
      33                 :            : 
      34                 :            :   public:
      35                 :            :     //! Number of sample space dimensions
      36                 :            :     static const std::size_t dim = 1;
      37                 :            : 
      38                 :            :     //! Key type
      39                 :            :     using key_type = long;
      40                 :            : 
      41                 :            :     //! Pair type
      42                 :            :     using pair_type = std::pair< const key_type, tk::real >;
      43                 :            : 
      44                 :            :     //! \brief Univariate PDF
      45                 :            :     //! \details The underlying container type is an unordered_map where the key
      46                 :            :     //!   is one bin id corresponding to the single sample space dimension, and
      47                 :            :     //!   the mapped value is the sample counter. The hasher functor used here
      48                 :            :     //!   is the default for the key type provided by the standard library.
      49                 :            :     using map_type = std::unordered_map< key_type, tk::real >;
      50                 :            : 
      51                 :            :     //! Empty constructor for Charm++
      52         [ +  - ]:       8766 :     explicit UniPDF() : m_binsize( 0 ), m_nsample( 0 ), m_pdf() {}
      53                 :            : 
      54                 :            :     //! Constructor: Initialize univariate PDF container
      55                 :            :     //! \param[in] bs Sample space bin size
      56                 :            :     explicit UniPDF( tk::real bs ) :
      57         [ +  - ]:       2459 :       m_binsize( bs ), m_nsample( 0 ), m_pdf() {}
      58                 :            : 
      59                 :            :     //! Accessor to number of samples
      60                 :            :     //! \return Number of samples collected
      61         [ +  - ]:     168996 :     std::size_t nsample() const noexcept { return m_nsample; }
      62                 :            : 
      63                 :            :     //! Add sample to univariate PDF
      64                 :            :     //! \param[in] sample Sample to insert
      65                 :    3828504 :     void add( tk::real sample ) {
      66                 :            :       Assert( m_binsize > 0, "Bin size must be positive" );
      67                 :    3828504 :       ++m_nsample;
      68                 :            :       fenv_t fe;
      69                 :    3828504 :       feholdexcept( &fe );
      70                 :    3828504 :       ++m_pdf[ std::lround( sample / m_binsize ) ];
      71                 :    3828504 :       feclearexcept( FE_UNDERFLOW );
      72                 :    3828504 :       feupdateenv( &fe );
      73                 :    3828504 :     }
      74                 :            : 
      75                 :            :     //! Add multiple samples from a PDF
      76                 :            :     //! \param[in] p PDF whose samples to add
      77                 :       6651 :     void addPDF( const UniPDF& p ) {
      78                 :       6651 :       m_binsize = p.binsize();
      79                 :       6651 :       m_nsample += p.nsample();
      80         [ +  + ]:     360689 :       for (const auto& e : p.map()) m_pdf[ e.first ] += e.second;
      81                 :       6651 :     }
      82                 :            : 
      83                 :            :     //! Zero bins
      84                 :            :     void zero() noexcept { m_nsample = 0; m_pdf.clear(); }
      85                 :            : 
      86                 :            :     //! Constant accessor to underlying PDF map
      87                 :            :     //! \return Constant reference to underlying map
      88                 :            :     const map_type& map() const noexcept { return m_pdf; }
      89                 :            : 
      90                 :            :     //! Constant accessor to bin size
      91                 :            :     //! \return Sample space bin size
      92         [ -  + ]:       7377 :     tk::real binsize() const noexcept { return m_binsize; }
      93                 :            : 
      94                 :            :     //! Return minimum and maximum bin ids of sample space
      95                 :            :     //! \return {min,max} Minimum and maximum of the bin ids
      96                 :        726 :     std::array< long, 2*dim > extents() const {
      97                 :            :       Assert( !m_pdf.empty(), "PDF empty" );
      98                 :        726 :       auto x = std::minmax_element( begin(m_pdf), end(m_pdf),
      99                 :            :                  []( const pair_type& a, const pair_type& b )
     100 [ +  + ][ +  + ]:      83414 :                  { return a.first < b.first; } );
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
         [ +  + ][ +  + ]
     101                 :        726 :       return {{ x.first->first, x.second->first }};
     102                 :            :     }
     103                 :            : 
     104                 :            :     //! Compute integral of the distribution across the whole sample space
     105                 :            :     //! \return Integral of the distribution
     106                 :            :     tk::real integral() const {
     107                 :            :       return std::accumulate( m_pdf.cbegin(), m_pdf.cend(), 0.0,
     108                 :            :         [&]( tk::real i, const pair_type& p ){
     109         [ +  - ]:     163071 :           return i + p.second; } ) / static_cast< tk::real >( m_nsample );
     110                 :            :     }
     111                 :            : 
     112                 :            :     /** @name Pack/Unpack: Serialize UniPDF object for Charm++ */
     113                 :            :     ///@{
     114                 :            :     //! Pack/Unpack serialize member function
     115                 :            :     //! \param[in,out] p Charm++'s PUP::er serializer object reference
     116                 :            :     // cppcheck-suppress constParameter
     117                 :      26298 :     void pup( PUP::er& p ) {
     118                 :      26298 :       p | m_binsize;
     119                 :      26298 :       p | m_nsample;
     120                 :      26298 :       p | m_pdf;
     121                 :      26298 :     }
     122                 :            :     //! \brief Pack/Unpack serialize operator|
     123                 :            :     //! \param[in,out] p Charm++'s PUP::er serializer object reference
     124                 :            :     //! \param[in,out] c UniPDF object reference
     125         [ +  - ]:      26298 :     friend void operator|( PUP::er& p, UniPDF& c ) { c.pup(p); }
     126                 :            :     ///@}
     127                 :            : 
     128                 :            :   private:
     129                 :            :     tk::real m_binsize;         //!< Sample space bin size
     130                 :            :     std::size_t m_nsample;      //!< Number of samples collected
     131                 :            :     map_type m_pdf;             //!< Probability density function
     132                 :            : };
     133                 :            : 
     134                 :            : //! Output univariate PDF to output stream
     135                 :            : //! \param[in,out] os Stream to output to
     136                 :            : //! \param[in] p PDF to output
     137                 :            : //! \return Updated stream
     138                 :            : //! \note Used for debugging.
     139                 :            : static inline
     140                 :            : std::ostream& operator<< ( std::ostream& os, const tk::UniPDF& p ) {
     141                 :            :   os << p.binsize() << ", " << p.nsample() << ": ";
     142                 :            :   std::map< typename tk::UniPDF::key_type, tk::real >
     143                 :            :     sorted( p.map().begin(), p.map().end() );
     144                 :            :   for (const auto& b : sorted) os << '(' << b.first << ',' << b.second << ") ";
     145                 :            :   return os;
     146                 :            : }
     147                 :            : 
     148                 :            : } // tk::
     149                 :            : 
     150                 :            : #endif // UniPDF_h