Program Listing for File PDFSegmentSchechter.hpp

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

#ifndef PDFSEGMENTSCHECHTER_HPP
#define PDFSEGMENTSCHECHTER_HPP

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

namespace pdfs {

    class PDFSegmentSchechter : public PDFSegment {
    public:

        PDFSegmentSchechter(double sMin, double sMax, double sStar, double alpha) :
            PDFSegment(sMin, sMax),
            sStar_(sStar),
            alpha_(alpha),
            norm_(1.0 / (
                    std::pow(sStar_, alpha_ + 1) *
                    (gsl_sf_gamma_inc(alpha_ + 1, sMin_ / sStar_) -
                    gsl_sf_gamma_inc(alpha_ + 1, sMax_ / sStar_)
                    ))
                )
            { }
        PDFSegmentSchechter(std::ifstream& file,
            FileFormats fmt,
            double &sMin,
            double &sMax,
            double &wgt);
        PDFSegmentSchechter(toml::node_view<const toml::node> node,
            FileFormats fmt,
            double &sMin,
            double &sMax,
            double &wgt);
        PDFSegmentSchechter(const PDFSegmentSchechter&) = default;
        auto operator=(const PDFSegmentSchechter&) -> PDFSegmentSchechter& = default;
        PDFSegmentSchechter(PDFSegmentSchechter&&) = default;
        auto operator=(PDFSegmentSchechter&&) -> PDFSegmentSchechter& = default;
        ~PDFSegmentSchechter() override = default;

        // Evaluation functions
        [[nodiscard]] auto operator()(double x) const -> double override {
            if (x < sMin_ || x > sMax_) {
                return 0.0;
            }
            return norm_ * std::pow(x, alpha_) * std::exp(-x / sStar_);
        }
        [[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_);
            return sStar_ * (
                gsl_sf_gamma_inc(alpha_ + 2, aClamped / sStar_) -
                gsl_sf_gamma_inc(alpha_ + 2, bClamped / sStar_)
            ) / (
                gsl_sf_gamma_inc(alpha_ + 1, aClamped / sStar_) -
                gsl_sf_gamma_inc(alpha_ + 1, bClamped / sStar_)
            );
        }
        [[nodiscard]] auto expectationValue() const -> double override {
            return expectationValue(sMin_, sMax_); // Expectation value over the full range of the segment
        }
        [[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_ * std::pow(sStar_, alpha_ + 1) * (
                gsl_sf_gamma_inc(alpha_ + 1, aClamped / sStar_) -
                gsl_sf_gamma_inc(alpha_ + 1, bClamped / sStar_)
            );
        }

        // Drawing function
        [[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)
            // Find the value of the deviate y by numerically solving
            // u = \int_a^y x^alpha exp(-x / sStar) dx /
            //     \int_a^b x^alpha exp(-x / sStar) dx
            // using a simple bisection search.
            const double gammaA = gsl_sf_gamma_inc(alpha_ + 1, aClamped / sStar_);
            const double gammaB = gsl_sf_gamma_inc(alpha_ + 1, bClamped / sStar_);
            const double denom = gammaA - gammaB;
            double y = aClamped;
            double step = (bClamped - aClamped) / 2.0;
            while (true) {
                y += step;
                const double gammaY = gsl_sf_gamma_inc(alpha_ + 1, y / sStar_);
                const double resid = u - ((gammaA - gammaY) / denom);
                if (fabs(resid) < 1e-8) {
                    break;  // Converged
                }
                if (resid < 0) {
                    y -= step;
                    step /= 2.0;
                    if (step < 1e-8) {
                        break;
                    }
                }
            }
            return y;
        }
        [[nodiscard]] auto draw() const -> double override {
            return draw(sMin_, sMax_);
        }

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

    private:
        double sStar_;
        double alpha_;
        double norm_;
    };

} // namespace pdfs

#endif // PDFSEGMENTSCHECHTER_HPP