.. _program_listing_file_src_pdfs_PDF.hpp:

Program Listing for File PDF.hpp
================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_pdfs_PDF.hpp>` (``src/pdfs/PDF.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   
   #ifndef PDF_HPP
   #define PDF_HPP
   
   #include "../utils/RandomDist.hpp"
   #include "../utils/RngThread.hpp"
   #include "PDFCommons.hpp"
   #include "PDFSegment.hpp"
   #include <algorithm>
   #include <array>
   #include <cmath>
   #include <concepts>
   #include <format>
   #include <limits>
   #include <memory>
   #include <ranges>
   #include <stdexcept>
   #include <string>
   #include <toml.hpp>
   #include <utility>
   #include <valarray>
   #include <vector>
   
   namespace pdfs {
   
       class PDF {
       public:
   
           // Constructors and destructor
           template <typename SC, typename WC>
           PDF(SC seg,
               WC wgt,
               SamplingMethods method = SamplingMethods::stopNearest,
               bool normalize = true) :
               seg_(std::move(seg)),
               method_(method),
               normalized_(normalize)
           {
               const auto n = std::size(wgt);
               const auto* ptr = std::data(wgt);
               if (n == 0 || ptr == nullptr) {
                   throw std::runtime_error("PDF: must have at least one segment");
               }
               wgt_ = std::valarray<double>(ptr, n);
               if (wgt_.min() <= 0) { // Safety check
                   throw std::runtime_error("PDF: elements of wgt must be positive");
               }
               if (seg_.size() != wgt_.size()) {
                   throw std::runtime_error("PDF: seg and wgt must be of equal size");
               }
               if (normalize) {
                   normalizePDF();
               }
               sMin_ = seg_.front()->getMin();
               sMax_ = seg_.front()->getMax();
               for (auto& s : seg_) {
                   sMin_ = std::min(sMin_, s->getMin());
                   sMax_ = std::max(sMax_, s->getMax());
               }
           }
   
           template <typename SegT>
               requires std::derived_from<SegT, PDFSegment>
           explicit PDF(std::unique_ptr<SegT> seg,
               const double wgt = 1.0,
               SamplingMethods method = SamplingMethods::stopNearest) :
               wgt_(wgt, 1),
               sMin_(seg->getMin()),
               sMax_(seg->getMax()),
               method_(method),
               normalized_(wgt == 1.0)
           {
               seg_.push_back(std::move(seg));
           }
   
           PDF() :
               sMin_(std::numeric_limits<double>::quiet_NaN()),
               sMax_(std::numeric_limits<double>::quiet_NaN()),
               method_(pdfs::SamplingMethods::none),
               normalized_(false)
           { }
           virtual ~PDF() = default;
   
           // Disallow copying, since we use unique_ptr
           // objects to hold segments
           PDF(const PDF&) = delete;
           auto operator=(const PDF&) -> PDF& = delete;
   
           // Default moves
           PDF(PDF&&) = default; 
           auto operator=(PDF&&) -> PDF& = default; 
   
           // Getters for internal state
           [[nodiscard]] auto valid() const { return !seg_.empty(); }
           [[nodiscard]] auto getMin() const -> double { return sMin_; }
           [[nodiscard]] auto getMax() const -> double { return sMax_; }
           [[nodiscard]] auto getWeights() const -> const std::valarray<double>& { return wgt_; }
   
           // Sampling policy getter and setter
           [[nodiscard]] auto getSampling() const -> SamplingMethods
           {
               return method_;
           }
           void setSampling(SamplingMethods method)
           {
               method_ = method;
           }
   
           // Normalization method
           auto normalizePDF() -> void
           {
               wgt_ /= wgt_.sum();
               normalized_ = true;
           }
           [[nodiscard]] auto normalized() const -> bool { return normalized_; }
           auto markNormalized() -> void { normalized_ = true; }
   
           // Serialization
           [[nodiscard]] virtual auto toToml() const -> std::string
           {
               if (!valid())
               {
                   throw std::logic_error("PDF::toToml: cannot serialize an invalid (empty) PDF");
               }
               if (method_ == SamplingMethods::none)
               {
                   throw std::logic_error("PDF::toToml: cannot serialize a PDF whose sampling "
                       "method is 'none', which no PDF descriptor can specify");
               }
               std::string out;
               out += "format = \"advanced\"\n";
               out += "method = \"" + std::string(samplingMethodName(method_)) + "\"\n";
               out += std::string("normalized = ") + (normalized_ ? "true" : "false") + "\n";
               int n = 1;
               for (auto const& [s,w] : std::views::zip(seg_,wgt_)) {
                   toml::table segTbl;
                   s->toToml(segTbl);
                   out += "\n[segment" + std::to_string(n++) + "]\n";
                   out += "type = " + tomlValueString(*segTbl.get("type")) + "\n";
                   out += "weight = " + tomlValueString(toml::value<double>(w)) + "\n";
                   for (auto const& [key, val] : segTbl) {
                       if (key.str() == "type") { continue; }
                       out += std::string(key.str()) + " = " + tomlValueString(val) + "\n";
                   }
               }
               return out;
           }
   
           // Evaluation functions
           virtual auto operator()(double x) const -> double
           {
               double p = 0.0;
               for (auto const& [s,w] : std::views::zip(seg_,wgt_)) {
                   if (s->getMin() < x && s->getMax() >= x) {
                       p += w * (*s)(x);
                   }
               }
               return p;
           }
           [[nodiscard]] virtual auto expectationValue(const double a, const double b) const -> double
           {
               double e = 0.0;
               double wSum = 0.0;
               for (auto const& [s,w] : std::views::zip(seg_,wgt_)) {
                   const double wi = w * s->integral(a, b);
                   e += wi * s->expectationValue(a, b);
                   wSum += wi;
               }
               return e / wSum;
           }
           [[nodiscard]] virtual auto expectationValue() const -> double
           {
               return expectationValue(sMin_, sMax_);
           }
           [[nodiscard]] virtual auto integral(const double a, const double b) const -> double {
               if (a >= b) {
                   return 0.0; // Interval is empty
               }
               double sum = 0.0;
               for (auto const& [s,w] : std::views::zip(seg_,wgt_)) {
                   sum += w * s->integral(a,b);  // Segment overlaps range
               }
               return sum;
           }
           [[nodiscard]] auto integral() const -> double { return wgt_.sum(); }
           [[nodiscard]] virtual auto deltas() const -> std::vector<std::pair<double, double>>
           {
               std::vector<std::pair<double, double>> result;
               for (auto const& [s,w] : std::views::zip(seg_,wgt_)) {
                   if (s->getMin() == s->getMax()) { result.emplace_back(s->getMin(), w); }
               }
               return result;
           }
   
           // Drawing functions
           [[nodiscard]]
           virtual auto draw(const double a = std::numeric_limits<double>::lowest(),
            const double b = std::numeric_limits<double>::max()) const -> double {
               std::vector<double> wLim(wgt_.size());
               for (auto const& [wL,s,w] : std::views::zip(wLim,seg_,wgt_)) {
                   wL = w * s->integral(a,b);
               }
               const utils::DiscreteDist dist(wLim);
               return seg_.at(dist(utils::rng()()))->draw(a,b);
           }
           [[nodiscard]] virtual auto draw(unsigned int nDraw,
               double a = std::numeric_limits<double>::lowest(),
               double b = std::numeric_limits<double>::max())
               const -> std::vector<double>
           {
               std::vector<double> result(nDraw);
               for (auto &r : result) { r = draw(a, b); }
               return result;
           }
           [[nodiscard]] virtual auto drawTarget(double target,
               double a = std::numeric_limits<double>::lowest(),
               double b = std::numeric_limits<double>::max()) const
               -> std::vector<double>
           {
               // A range containing no probability -- either degenerate
               // (a >= b) or lying entirely outside every segment's
               // support -- has no samples to draw, regardless of target.
               // Falling through would build a utils::DiscreteDist from
               // all-zero per-segment weights (which throws), or divide
               // by a zero integral in expectationValue(a, b), giving a
               // NaN sample count or Poisson mean (or, for the stop
               // methods, an endless loop of zero-valued draws).
               // integral() consumes no rng state, so this check does not
               // perturb the draw sequence.
               if (a >= b || !(integral(a, b) > 0.0)) { return {}; }
   
               std::vector<double> sample;
               switch (method_)
               {
                   case SamplingMethods::stopNearest:
                   case SamplingMethods::stopBefore:
                   case SamplingMethods::stopAfter:
                   case SamplingMethods::stop50:
                   {
                       // Stop methods
                       double sum = 0.0;
                       while (true)
                       {
                           auto s = draw(a, b);
                           if (sum + s <= target)
                           {
                               // Taget not yet reached; add to sample and continue
                               sample.push_back(s);
                               sum += s;
                           } else {
                               // Target reached; decide whether to keep last number based on policy
                               const bool keep =
                                   (
                                       method_ == SamplingMethods::stopNearest &&
                                       s + sum - target < target - sum
                                   ) || (
                                       method_ == SamplingMethods::stopAfter
                                   ) || (
                                       method_ == SamplingMethods::stop50 &&
                                       utils::canonical(utils::rng()()) > 0.5
                                   );
                               if (keep) { sample.push_back(s); }
                               break;
                           }
                       }
                       return sample;
                   }
                   case SamplingMethods::number:
                   {
                       // Number method: compute number of samples, then draw
                       const auto nSamp = 
                           static_cast<unsigned int>(
                           std::round(target / expectationValue(a,b))
                           );
                       return draw(nSamp, a, b);
                   }
                   case SamplingMethods::poisson:
                   {
                       // Poisson method: draw number of samples, then draw
                       utils::PoissonDist dist(target / expectationValue(a,b));
                       const auto nSamp = dist(utils::rng()());
                       return draw(nSamp, a, b);
                   }
                   case SamplingMethods::sorted:
                   {
                       // Sorted sampling: first draw repeatedly using number
                       // method, until target is exceeded, then sort list and
                       // keep or discard last entry following stop nearest method.
                       std::vector<double> samples;
                       double sum = 0.0;
                       while (sum < target) {
                           const auto nSamp = 
                               static_cast<unsigned int>(std::round((target - sum) / expectationValue(a,b)));
                           auto sampleSet = draw(nSamp, a, b);
                           samples.resize(samples.size() + nSamp);
                           samples.insert(samples.end(), sampleSet.begin(), sampleSet.end());
                           for (auto s : sampleSet) { sum += s; }
                       }
                       std::ranges::sort(samples);
                       if (std::fabs(target - (sum - samples.back())) <
                           std::fabs(target - sum))
                       {
                           samples.pop_back();
                       }
                       return samples;
                   }
                   case SamplingMethods::none:
                   {
                       throw std::runtime_error("PDF: invalid sampling method");
                   }
               }
               // Every enumerator of SamplingMethods is handled above, but
               // method_'s underlying type can in principle hold a value
               // outside the declared enumerators, so the compiler cannot
               // prove the switch is exhaustive; make that explicit rather
               // than falling off the end of a non-void function.
               throw std::runtime_error("PDF: invalid sampling method");
           }
   
       protected:
   
           [[nodiscard]] static auto tomlValueString(const toml::node& val) -> std::string
           {
               if (const auto* str = val.as_string()) { return "\"" + str->get() + "\""; }
               if (const auto* b = val.as_boolean()) { return b->get() ? "true" : "false"; }
               if (const auto* i = val.as_integer()) { return std::to_string(i->get()); }
               if (const auto* d = val.as_floating_point())
               {
                   std::string s = std::format("{}", d->get());
                   if (s.find_first_of(".eni") == std::string::npos) { s += ".0"; }
                   return s;
               }
               throw std::logic_error("PDF::toToml: unsupported toml value type");
           }
   
   
           std::vector<std::unique_ptr<PDFSegment> > seg_; 
           std::valarray<double> wgt_;     
           double sMin_;                   
           double sMax_;                   
           SamplingMethods method_;        
           bool normalized_;               
       };
   
   } // namespace pdfs
   
   #endif // PDF_HPP