Program Listing for File PDF.hpp
↰ Return to documentation for file (src/pdfs/PDF.hpp)
#ifndef PDF_HPP
#define PDF_HPP
#include "../utils/RandomDist.hpp"
#include "../utils/RngThread.hpp"
#include "PDFCommons.hpp"
#include "PDFSegment.hpp"
#include <algorithm>
#include <array>
#include <cmath>
#include <concepts>
#include <format>
#include <limits>
#include <memory>
#include <ranges>
#include <stdexcept>
#include <string>
#include <toml.hpp>
#include <utility>
#include <valarray>
#include <vector>
namespace pdfs {
class PDF {
public:
// Constructors and destructor
template <typename SC, typename WC>
PDF(SC seg,
WC wgt,
SamplingMethods method = SamplingMethods::stopNearest,
bool normalize = true) :
seg_(std::move(seg)),
method_(method),
normalized_(normalize)
{
const auto n = std::size(wgt);
const auto* ptr = std::data(wgt);
if (n == 0 || ptr == nullptr) {
throw std::runtime_error("PDF: must have at least one segment");
}
wgt_ = std::valarray<double>(ptr, n);
if (wgt_.min() <= 0) { // Safety check
throw std::runtime_error("PDF: elements of wgt must be positive");
}
if (seg_.size() != wgt_.size()) {
throw std::runtime_error("PDF: seg and wgt must be of equal size");
}
if (normalize) {
normalizePDF();
}
sMin_ = seg_.front()->getMin();
sMax_ = seg_.front()->getMax();
for (auto& s : seg_) {
sMin_ = std::min(sMin_, s->getMin());
sMax_ = std::max(sMax_, s->getMax());
}
}
template <typename SegT>
requires std::derived_from<SegT, PDFSegment>
explicit PDF(std::unique_ptr<SegT> seg,
const double wgt = 1.0,
SamplingMethods method = SamplingMethods::stopNearest) :
wgt_(wgt, 1),
sMin_(seg->getMin()),
sMax_(seg->getMax()),
method_(method),
normalized_(wgt == 1.0)
{
seg_.push_back(std::move(seg));
}
PDF() :
sMin_(std::numeric_limits<double>::quiet_NaN()),
sMax_(std::numeric_limits<double>::quiet_NaN()),
method_(pdfs::SamplingMethods::none),
normalized_(false)
{ }
virtual ~PDF() = default;
// Disallow copying, since we use unique_ptr
// objects to hold segments
PDF(const PDF&) = delete;
auto operator=(const PDF&) -> PDF& = delete;
// Default moves
PDF(PDF&&) = default;
auto operator=(PDF&&) -> PDF& = default;
// Getters for internal state
[[nodiscard]] auto valid() const { return !seg_.empty(); }
[[nodiscard]] auto getMin() const -> double { return sMin_; }
[[nodiscard]] auto getMax() const -> double { return sMax_; }
[[nodiscard]] auto getWeights() const -> const std::valarray<double>& { return wgt_; }
// Sampling policy getter and setter
[[nodiscard]] auto getSampling() const -> SamplingMethods
{
return method_;
}
void setSampling(SamplingMethods method)
{
method_ = method;
}
// Normalization method
auto normalizePDF() -> void
{
wgt_ /= wgt_.sum();
normalized_ = true;
}
[[nodiscard]] auto normalized() const -> bool { return normalized_; }
auto markNormalized() -> void { normalized_ = true; }
// Serialization
[[nodiscard]] virtual auto toToml() const -> std::string
{
if (!valid())
{
throw std::logic_error("PDF::toToml: cannot serialize an invalid (empty) PDF");
}
if (method_ == SamplingMethods::none)
{
throw std::logic_error("PDF::toToml: cannot serialize a PDF whose sampling "
"method is 'none', which no PDF descriptor can specify");
}
std::string out;
out += "format = \"advanced\"\n";
out += "method = \"" + std::string(samplingMethodName(method_)) + "\"\n";
out += std::string("normalized = ") + (normalized_ ? "true" : "false") + "\n";
int n = 1;
for (auto const& [s,w] : std::views::zip(seg_,wgt_)) {
toml::table segTbl;
s->toToml(segTbl);
out += "\n[segment" + std::to_string(n++) + "]\n";
out += "type = " + tomlValueString(*segTbl.get("type")) + "\n";
out += "weight = " + tomlValueString(toml::value<double>(w)) + "\n";
for (auto const& [key, val] : segTbl) {
if (key.str() == "type") { continue; }
out += std::string(key.str()) + " = " + tomlValueString(val) + "\n";
}
}
return out;
}
// Evaluation functions
virtual auto operator()(double x) const -> double
{
double p = 0.0;
for (auto const& [s,w] : std::views::zip(seg_,wgt_)) {
if (s->getMin() < x && s->getMax() >= x) {
p += w * (*s)(x);
}
}
return p;
}
[[nodiscard]] virtual auto expectationValue(const double a, const double b) const -> double
{
double e = 0.0;
double wSum = 0.0;
for (auto const& [s,w] : std::views::zip(seg_,wgt_)) {
const double wi = w * s->integral(a, b);
e += wi * s->expectationValue(a, b);
wSum += wi;
}
return e / wSum;
}
[[nodiscard]] virtual auto expectationValue() const -> double
{
return expectationValue(sMin_, sMax_);
}
[[nodiscard]] virtual auto integral(const double a, const double b) const -> double {
if (a >= b) {
return 0.0; // Interval is empty
}
double sum = 0.0;
for (auto const& [s,w] : std::views::zip(seg_,wgt_)) {
sum += w * s->integral(a,b); // Segment overlaps range
}
return sum;
}
[[nodiscard]] auto integral() const -> double { return wgt_.sum(); }
[[nodiscard]] virtual auto deltas() const -> std::vector<std::pair<double, double>>
{
std::vector<std::pair<double, double>> result;
for (auto const& [s,w] : std::views::zip(seg_,wgt_)) {
if (s->getMin() == s->getMax()) { result.emplace_back(s->getMin(), w); }
}
return result;
}
// Drawing functions
[[nodiscard]]
virtual auto draw(const double a = std::numeric_limits<double>::lowest(),
const double b = std::numeric_limits<double>::max()) const -> double {
std::vector<double> wLim(wgt_.size());
for (auto const& [wL,s,w] : std::views::zip(wLim,seg_,wgt_)) {
wL = w * s->integral(a,b);
}
const utils::DiscreteDist dist(wLim);
return seg_.at(dist(utils::rng()()))->draw(a,b);
}
[[nodiscard]] virtual auto draw(unsigned int nDraw,
double a = std::numeric_limits<double>::lowest(),
double b = std::numeric_limits<double>::max())
const -> std::vector<double>
{
std::vector<double> result(nDraw);
for (auto &r : result) { r = draw(a, b); }
return result;
}
[[nodiscard]] virtual auto drawTarget(double target,
double a = std::numeric_limits<double>::lowest(),
double b = std::numeric_limits<double>::max()) const
-> std::vector<double>
{
// A range containing no probability -- either degenerate
// (a >= b) or lying entirely outside every segment's
// support -- has no samples to draw, regardless of target.
// Falling through would build a utils::DiscreteDist from
// all-zero per-segment weights (which throws), or divide
// by a zero integral in expectationValue(a, b), giving a
// NaN sample count or Poisson mean (or, for the stop
// methods, an endless loop of zero-valued draws).
// integral() consumes no rng state, so this check does not
// perturb the draw sequence.
if (a >= b || !(integral(a, b) > 0.0)) { return {}; }
std::vector<double> sample;
switch (method_)
{
case SamplingMethods::stopNearest:
case SamplingMethods::stopBefore:
case SamplingMethods::stopAfter:
case SamplingMethods::stop50:
{
// Stop methods
double sum = 0.0;
while (true)
{
auto s = draw(a, b);
if (sum + s <= target)
{
// Taget not yet reached; add to sample and continue
sample.push_back(s);
sum += s;
} else {
// Target reached; decide whether to keep last number based on policy
const bool keep =
(
method_ == SamplingMethods::stopNearest &&
s + sum - target < target - sum
) || (
method_ == SamplingMethods::stopAfter
) || (
method_ == SamplingMethods::stop50 &&
utils::canonical(utils::rng()()) > 0.5
);
if (keep) { sample.push_back(s); }
break;
}
}
return sample;
}
case SamplingMethods::number:
{
// Number method: compute number of samples, then draw
const auto nSamp =
static_cast<unsigned int>(
std::round(target / expectationValue(a,b))
);
return draw(nSamp, a, b);
}
case SamplingMethods::poisson:
{
// Poisson method: draw number of samples, then draw
utils::PoissonDist dist(target / expectationValue(a,b));
const auto nSamp = dist(utils::rng()());
return draw(nSamp, a, b);
}
case SamplingMethods::sorted:
{
// Sorted sampling: first draw repeatedly using number
// method, until target is exceeded, then sort list and
// keep or discard last entry following stop nearest method.
std::vector<double> samples;
double sum = 0.0;
while (sum < target) {
const auto nSamp =
static_cast<unsigned int>(std::round((target - sum) / expectationValue(a,b)));
auto sampleSet = draw(nSamp, a, b);
samples.resize(samples.size() + nSamp);
samples.insert(samples.end(), sampleSet.begin(), sampleSet.end());
for (auto s : sampleSet) { sum += s; }
}
std::ranges::sort(samples);
if (std::fabs(target - (sum - samples.back())) <
std::fabs(target - sum))
{
samples.pop_back();
}
return samples;
}
case SamplingMethods::none:
{
throw std::runtime_error("PDF: invalid sampling method");
}
}
// Every enumerator of SamplingMethods is handled above, but
// method_'s underlying type can in principle hold a value
// outside the declared enumerators, so the compiler cannot
// prove the switch is exhaustive; make that explicit rather
// than falling off the end of a non-void function.
throw std::runtime_error("PDF: invalid sampling method");
}
protected:
[[nodiscard]] static auto tomlValueString(const toml::node& val) -> std::string
{
if (const auto* str = val.as_string()) { return "\"" + str->get() + "\""; }
if (const auto* b = val.as_boolean()) { return b->get() ? "true" : "false"; }
if (const auto* i = val.as_integer()) { return std::to_string(i->get()); }
if (const auto* d = val.as_floating_point())
{
std::string s = std::format("{}", d->get());
if (s.find_first_of(".eni") == std::string::npos) { s += ".0"; }
return s;
}
throw std::logic_error("PDF::toToml: unsupported toml value type");
}
std::vector<std::unique_ptr<PDFSegment> > seg_;
std::valarray<double> wgt_;
double sMin_;
double sMax_;
SamplingMethods method_;
bool normalized_;
};
} // namespace pdfs
#endif // PDF_HPP