Program Listing for File PDFSegmentLognormal.hpp

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

#ifndef PDFSEGMENTLOGNORMAL_HPP
#define PDFSEGMENTLOGNORMAL_HPP

#include "../utils/RandomDist.hpp"
#include "../utils/RngThread.hpp"
#include "PDFCommons.hpp"
#include "PDFSegment.hpp"
#include <algorithm>
#include <cmath>
#include <fstream>
#include <numbers>
#include <toml.hpp>

namespace pdfs {

    class PDFSegmentLognormal : public PDFSegment {
    public:

        // Constructor and destructor
        PDFSegmentLognormal(double sMin, double sMax, double mean, double stddev) :
            PDFSegment(sMin, sMax),
            mean_(mean),
            stddev_(stddev),
            logMean_(std::log(mean_)),
            dispDex_(stddev / std::numbers::ln10),
            root2dev_(std::numbers::sqrt2 * stddev_)
            {
                norm_ = std::numbers::sqrt2 * std::numbers::inv_sqrtpi /
                    stddev_ / (
                    std::erf( -std::log(sMin_/mean_) / root2dev_ ) -
                    std::erf( -std::log(sMax_/mean_) / root2dev_ )
                    );
            }
        PDFSegmentLognormal(std::ifstream& file,
            FileFormats fmt,
            double &sMin,
            double &sMax,
            double &wgt);
        PDFSegmentLognormal(toml::node_view<const toml::node> node,
            FileFormats fmt,
            double &sMin,
            double &sMax,
            double &wgt);
        PDFSegmentLognormal(const PDFSegmentLognormal&) = default;
        auto operator=(const PDFSegmentLognormal&) -> PDFSegmentLognormal& = default;
        PDFSegmentLognormal(PDFSegmentLognormal&&) = default;
        auto operator=(PDFSegmentLognormal&&) -> PDFSegmentLognormal& = default;
       ~PDFSegmentLognormal() 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( -0.5 * std::pow( std::log(x / mean_) / stddev_, 2) ) / 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 logA = std::log(aClamped);
            const double logB = std::log(bClamped);
            return mean_ * std::exp(std::pow(stddev_, 2) / 2) *
                (
                    std::erf( ((logMean_ - logA + std::pow(stddev_, 2) ) / root2dev_ ) ) -
                    std::erf( ((logMean_ - logB + std::pow(stddev_, 2) ) / root2dev_ ) )
                ) / (
                    std::erf( (logB - logMean_) / root2dev_ ) -
                    std::erf( (logA - logMean_) / root2dev_ )
                );
        }
        [[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 expectation value calculation
            }
            const double aClamped = std::max(a, sMin_);
            const double bClamped = std::min(b, sMax_);
            const double logA = std::log(aClamped);
            const double logB = std::log(bClamped);
            return norm_ /
                (std::numbers::inv_sqrtpi * std::numbers::sqrt2) *
                stddev_ * (
                std::erf( (logMean_ - logA) / root2dev_ ) -
                std::erf( (logMean_ - logB) / root2dev_ )
                );
        }

        // 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
            }
            utils::LognormalDist dist(logMean_, stddev_);
            double sample = NAN;
            while (true) {
                sample = dist(utils::rng()());
                if (sample >= aClamped && sample <= bClamped) {break; } // Rejection sampling
            }
            return sample;
        }
        [[nodiscard]] auto draw() const -> double override {
            return draw(sMin_, sMax_);
        }

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

    private:
        double mean_;
        double stddev_;
        double norm_;
        double logMean_;
        double dispDex_;
        double root2dev_;
    };

} // namespace pdfs

#endif // PDFSEGMENTLOGNORMAL_HPP