Program Listing for File PDFSegmentExponential.hpp

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

#ifndef PDFSEGMENTEXPONENTIAL_HPP
#define PDFSEGMENTEXPONENTIAL_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 PDFSegmentExponential : public PDFSegment {
    public:

        // Constructor and destructor
        PDFSegmentExponential(double sMin, double sMax, double scale) :
            PDFSegment(sMin, sMax), scale_(scale)
        {
            // Calculate normalization constant for the PDF segment
            const double exMin = std::exp(-sMin_ / scale_);
            const double exMax = std::exp(-sMax_ / scale_);
            norm_ = 1.0 / (scale_ * (exMin - exMax));
        }
        PDFSegmentExponential(std::ifstream& file,
            FileFormats fmt,
            double &sMin,
            double &sMax,
            double &wgt);
        PDFSegmentExponential(toml::node_view<const toml::node> node,
            FileFormats fmt,
            double &sMin,
            double &sMax,
            double &wgt);
        PDFSegmentExponential(const PDFSegmentExponential&) = default;
        auto operator=(const PDFSegmentExponential&) -> PDFSegmentExponential& = default;
        PDFSegmentExponential(PDFSegmentExponential&&) = default;
        auto operator=(PDFSegmentExponential&&) -> PDFSegmentExponential& = default;
        ~PDFSegmentExponential() 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::exp(-x / scale_); // 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_);
            const double num1 = aClamped * std::exp(-aClamped / scale_);
            const double num2 = bClamped * std::exp(-bClamped / scale_);
            const double denom1 = std::exp(-aClamped / scale_);
            const double denom2 = std::exp(-bClamped / scale_);
            return scale_ + ((num1 - num2) / (denom1 - denom2));
        }
        [[nodiscard]] auto expectationValue() const -> double override {
            const double num1 = sMin_ * std::exp(-sMin_ / scale_);
            const double num2 = sMax_ * std::exp(-sMax_ / scale_);
            const double denom1 = std::exp(-sMin_ / scale_);
            const double denom2 = std::exp(-sMax_ / scale_);
            return scale_ + ((num1 - num2) / (denom1 - denom2));
        }
        [[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_);
            return norm_ * scale_ * (std::exp(-aClamped / scale_) - std::exp(-bClamped / scale_));
        }

        // 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_);
            const double u = utils::canonical(utils::rng()());
            const double expA = std::exp(-aClamped / scale_);
            const double expB= std::exp(-bClamped / scale_);
            return -scale_ * std::log(expA - (u * (expA - expB))); // Inverse transform sampling for exponential distribution
        }
        [[nodiscard]] auto draw() const -> double override {
            return draw(sMin_, sMax_);
        }

        // Serialization
        void toToml(toml::table& tbl) const override;

    private:
        double scale_;
        double norm_;
    };

} // namespace pdfs

#endif // PDFSEGMENTEXPONENTIAL_HPP