.. _program_listing_file_src_pdfs_PDFSegmentPowerlaw.hpp:

Program Listing for File PDFSegmentPowerlaw.hpp
===============================================

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

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

.. code-block:: cpp

   
   #ifndef PDFSEGMENTPOWERLAW_HPP
   #define PDFSEGMENTPOWERLAW_HPP
   
   #include "../utils/RandomDist.hpp"
   #include "../utils/RngThread.hpp"
   #include "PDFCommons.hpp"
   #include "PDFSegment.hpp"
   #include <algorithm>
   #include <cmath>
   #include <fstream>
   #include <toml.hpp>
   
   namespace pdfs {
   
       class PDFSegmentPowerlaw : public PDFSegment {
       public:
   
           // Constructor and destructor
           PDFSegmentPowerlaw(double sMin, double sMax, double alpha) :
               PDFSegment(sMin, sMax), alpha_(alpha) {
                   // Calculate normalization constant for the PDF segment
                   if (alpha_ != -1) {
                       norm_ = (alpha_ + 1) / (std::pow(sMax_, alpha_ + 1) - std::pow(sMin_, alpha_ + 1));
                   } else {
                       norm_ = 1.0 / std::log(sMax_ / sMin_);
                   }
               }
           PDFSegmentPowerlaw(std::ifstream& file,
               FileFormats fmt,
               double &sMin,
               double &sMax,
               double &wgt);
           PDFSegmentPowerlaw(toml::node_view<const toml::node> node,
               FileFormats fmt,
               double &sMin,
               double &sMax,
               double &wgt);
           PDFSegmentPowerlaw(const PDFSegmentPowerlaw&) = default;
           auto operator=(const PDFSegmentPowerlaw&) -> PDFSegmentPowerlaw& = default;
           PDFSegmentPowerlaw(PDFSegmentPowerlaw&&) = default;
           auto operator=(PDFSegmentPowerlaw&&) -> PDFSegmentPowerlaw& = default;
           ~PDFSegmentPowerlaw() override = default;
   
           // Evaluation functions
           [[nodiscard]] auto operator()(double x) const -> double override {
               if (x < sMin_ || x > sMax_) {
                   return 0.0; // PDF is zero outside the segment
               }
               return norm_ * std::pow(x, alpha_); // PDF value at x
           }
   
           [[nodiscard]] auto expectationValue(const double a, const double b) const -> double override {
               if (a >= b || a > sMax_ || b < sMin_) {
                   return 0.0; // Invalid range for expectation value calculation
               }
               if (a == sMax_) {
                   return sMax_; // Handle edge cases
               }
               if (b == sMin_) {
                   return sMin_; // Handle edge cases
               }
               const double aClamped = std::max(a, sMin_);
               const double bClamped = std::min(b, sMax_);
               if (alpha_ != -1 && alpha_ != -2) {
                   return (alpha_ + 1) / (alpha_ + 2) *
                       (std::pow(bClamped, alpha_ + 2) - std::pow(aClamped, alpha_ + 2)) /
                       (std::pow(bClamped, alpha_ + 1) - std::pow(aClamped, alpha_ + 1));
               }
               if (alpha_ == -1) {
                   return (bClamped - aClamped) / (std::log(bClamped / aClamped));
               }
               return aClamped * bClamped * std::log(bClamped / aClamped) / (bClamped - aClamped);
           }
   
           [[nodiscard]] auto expectationValue() const -> double override {
               return expectationValue(sMin_, sMax_);
           }
   
           [[nodiscard]] auto integral(const double a, const double b) const -> double override {
               if (a >= b || a >= sMax_ || b <= sMin_) {
                   return 0.0; // Invalid range for integral calculation
               }
               const double aClamped = std::max(a, sMin_);
               const double bClamped = std::min(b, sMax_);
               if (alpha_ != -1) {
                   // Special case
                   return norm_ / (alpha_ + 1) *
                       (std::pow(bClamped, alpha_ + 1) - std::pow(aClamped, alpha_ + 1));
               }
               // General case
               return norm_ * std::log(bClamped / aClamped);
           }
   
           // Drawing functions
           [[nodiscard]] auto draw(const double a, const double b) const -> double override {
               const double aClamped = std::max(a, sMin_);
               const double bClamped = std::min(b, sMax_);
               if (aClamped >= bClamped) {
                   return 0.0; // Invalid range for drawing
               }
               const double u = utils::canonical(utils::rng()()); // Uniform random number in [0, 1)
               if (alpha_ != -1)
               {
                   // Special case
                   return std::pow(
                       (u * std::pow(bClamped, alpha_ + 1)) +
                       ((1 - u) * std::pow(aClamped, alpha_ + 1)),
                       1.0 / (alpha_ + 1)
                   );
               }
               // General case
               return aClamped * std::pow(bClamped / aClamped, u);
           }
           [[nodiscard]] auto draw() const -> double override {
               return draw(sMin_, sMax_);
           }
   
           // Serialization
           void toToml(toml::table& tbl) const override;
   
       private:
           double alpha_;  
           double norm_;   
       };
   
   } // namespace pdfs
   
   #endif // PDFSEGMENTPOWERLAW_HPP