VectorRng.h

Go to the documentation of this file.

#pragma once
 
 
#include "math/rng/Rng.h" // not really needed, but it doesn't make sense to use VectorRng without some Rng object
#include "objects/wrappers/Function.h"
 
NAMESPACE_SPH_BEGIN
 
template <typename TScalarRng>
class VectorRng : public Noncopyable {
private:
    TScalarRng rngImpl;
 
public:
    VectorRng() {
        static_assert(!std::is_reference<TScalarRng>::value, "Cannot be used for references");
    }
 
    template <typename... TArgs>
    explicit VectorRng(TArgs&&... args)
        : rngImpl(std::forward<TArgs>(args)...) {}
 
    VectorRng(VectorRng&& other)
        : rngImpl(std::move(other.rngImpl)) {
        static_assert(!std::is_reference<TScalarRng>::value, "Cannot be used for references");
    }
 
    VectorRng& operator=(VectorRng&& other) {
        static_assert(!std::is_reference<TScalarRng>::value, "Cannot be used for references");
        rngImpl = std::move(other.rngImpl);
        return *this;
    }
 
    Vector operator()() {
        return Vector(rngImpl(0), rngImpl(1), rngImpl(2));
    }
 
    Float getAdditional(const Size i) {
        SPH_ASSERT(i >= 3);
        return rngImpl(i);
    }
};
 
template <typename TScalarRng>
class VectorPdfRng : public Noncopyable {
private:
    Box box;
    VectorRng<TScalarRng> vectorRng;
    Function<Float(const Vector&)> pdf;
    Function<Float(const Vector&)> jacobian;
    Float maxPdf;
 
public:
    VectorPdfRng(
        const Box& box,
        const Function<Float(const Vector&)>& pdf = [](const Vector&) { return 1._f; },
        const Function<Float(const Vector&)>& jacobian = [](const Vector&) { return 1._f; },
        const Float maximalPdf = 0._f)
        : box(box)
        , pdf(pdf)
        , jacobian(jacobian)
        , maxPdf(maximalPdf) {
        if (maxPdf == 0._f) {
            // step for finding pdf maximum
            const Vector delta = 0.05_f * box.size();
            box.iterate(delta, [&](const Vector& v) { //
                this->maxPdf = max(this->maxPdf, this->pdf(v) * this->jacobian(v));
            });
            SPH_ASSERT(maxPdf > 0._f);
        }
    }
 
    Vector operator()() {
        SPH_ASSERT(getLength(box.size()) > 0._f);
        for (Size i = 0; i < 10000; ++i) {
            Vector v = vectorRng() * box.size() + box.lower();
            const Float p = vectorRng.getAdditional(4) * maxPdf;
            if (p < pdf(v) * jacobian(v)) {
                return v;
            }
        }
        SPH_ASSERT(false, "Couldn't generate vector in 1000 iterations");
        return Vector(0._f);
    }
};
 
 
NAMESPACE_SPH_END