Program Listing for File PDFSegmentPowerlaw.hpp

↰ Return to documentation for file (src/pdfs/PDFSegmentPowerlaw.hpp)

#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