Uniform.h

Go to the documentation of this file.

#pragma once

 
#include <cmath>
#include <expected>
#include <limits>
#include <optional>
#include <random>
#include <utility>
 
#include <EmrysFilter/Core/Forwarding.h>
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Core/Math.h>
#include <EmrysFilter/Core/Print.h>
#include <EmrysFilter/Distributions/Moments.h>
#include <EmrysFilter/Distributions/SampleSet.h>
#include <EmrysFilter/Distributions/Sampling.h>
#include <EmrysFilter/Spaces/Composite.h>
#include <EmrysFilter/Spaces/Indexing.h>
#include <EmrysFilter/Spaces/Space.h>
 
namespace Emrys
{

enum class UniformError
{
    Empty,            
    NonFinite,        
    WiderThanPeriod,  
};

template <Space Space_, RealScalar Scalar_ = double>
class Uniform
{
public:
 
    using Scalar = Scalar_;                               
    using Index = Eigen::Index;                           
    using SpaceType = Space_;                             
    using MomentsType = Moments<Space_, Scalar_>;         
    using VectorType = typename MomentsType::VectorType;  
    using MatrixType = typename MomentsType::MatrixType;  

    static constexpr int ValueDimensionAtCompileTime = MomentsType::ValueDimensionAtCompileTime;

    static constexpr int MaxValueDimensionAtCompileTime = MomentsType::MaxValueDimensionAtCompileTime;

    static constexpr int DifferenceDimensionAtCompileTime = MomentsType::DifferenceDimensionAtCompileTime;

    static constexpr int MaxDifferenceDimensionAtCompileTime = MomentsType::MaxDifferenceDimensionAtCompileTime;
 
    static_assert(
        internal::hasEqualSizes<Space_>,
        "a uniform box needs a space whose values and differences have the same coordinates"
    );

    explicit Uniform(Space_ space = {})
        requires(!isLinear<Space_>)
        : m_space{std::move(space)}
    {
        Index const n = Emrys::valueDimension(m_space);
        m_lower.resize(n);
        m_width.resize(n);
        m_period.resize(n);
        internal::forEachLeaf(m_space, [&](Index offset, auto const & leaf) {
            using Leaf = std::remove_cvref_t<decltype(leaf)>;
            static_assert(hasDefaultInterval<Leaf>(), "a coordinate has no default interval; give the box explicitly");
            constexpr auto found = internal::leafSupport<Leaf>();
            constexpr internal::LeafSupport support = *found;
            constexpr double width = support.upper - support.lower;
            auto const rows = Eigen::seqN(offset, leaf.valueDimension());
            m_lower(rows).setConstant(Scalar(support.lower));
            m_width(rows).setConstant(Scalar(width));
            // Only a Periodic coordinate has a default interval.
            static_assert(support.periodic);
            m_period(rows).setConstant(Scalar(width));
        });
        m_logVolume = m_width.array().log().sum();
    }

    template <MatrixExpression Lower, typename UpperDerived>
    [[nodiscard]] static std::expected<Uniform, UniformError>
    create(Space_ space, Lower && lower, Eigen::MatrixBase<UpperDerived> const & upper)
    {
        VectorType start = std::forward<Lower>(lower);
        EMRYS_FILTER_ASSERT(upper.size() == start.size() && "give one upper bound per lower bound");
        VectorType width = upper - start;
        internal::fitValueDimension(space, start.size());
        if (!allFinite(start) || !allFinite(width))
        {
            return std::unexpected{UniformError::NonFinite};
        }
        if (!(width.array() > Scalar(0)).all())
        {
            return std::unexpected{UniformError::Empty};
        }
        auto period = VectorType{VectorType::Zero(start.size())};
        bool fits = true;
        internal::forEachLeaf(space, [&](Index offset, auto const & leaf) {
            constexpr auto support = internal::leafSupport<std::remove_cvref_t<decltype(leaf)>>();
            if constexpr (support.has_value() && support->periodic)
            {
                constexpr double leafPeriod = support->upper - support->lower;
                for (Index i = offset; i < offset + leaf.valueDimension(); ++i)
                {
                    fits = fits && width(i) <= Scalar(leafPeriod);
                    period(i) = Scalar(leafPeriod);
                }
            }
        });
        if (!fits)
        {
            return std::unexpected{UniformError::WiderThanPeriod};
        }
        return Uniform{std::move(space), std::move(start), std::move(width), std::move(period)};
    }

    template <MatrixExpression Lower, typename UpperDerived>
        requires internal::SizedByLength<Space_>
    [[nodiscard]] static std::expected<Uniform, UniformError>
    create(Lower && lower, Eigen::MatrixBase<UpperDerived> const & upper)
    {
        return create(Space_{}, std::forward<Lower>(lower), upper);
    }

    [[nodiscard]] Index
    valueDimension() const
    {
        return m_lower.size();
    }

    [[nodiscard]] Index
    differenceDimension() const
    {
        return m_lower.size();
    }

    [[nodiscard]] Space_ const &
    space() const
    {
        return m_space;
    }

    [[nodiscard]] VectorType const &
    lower() const
    {
        return m_lower;
    }

    [[nodiscard]] VectorType
    upper() const
    {
        return m_lower + m_width;
    }

    [[nodiscard]] VectorType const &
    width() const
    {
        return m_width;
    }

    [[nodiscard]] Scalar
    logVolume() const
    {
        return m_logVolume;
    }

    template <typename Derived>
    [[nodiscard]] bool
    contains(Eigen::MatrixBase<Derived> const & x) const
    {
        static_assert(Derived::IsVectorAtCompileTime, "contains() expects a vector");
        EMRYS_FILTER_ASSERT(x.size() == valueDimension());
        for (Index i = 0; i < valueDimension(); ++i)
        {
            Scalar offset = Scalar(x(i)) - m_lower(i);
            if (m_period(i) > Scalar(0))
            {
                offset -= m_period(i) * floorOf<Scalar>(offset / m_period(i));
            }
            if (!(offset >= Scalar(0) && offset <= m_width(i)))
            {
                return false;
            }
        }
        return true;
    }

    template <typename Derived>
    [[nodiscard]] Scalar
    logDensity(Eigen::MatrixBase<Derived> const & x) const
    {
        return contains(x) ? -m_logVolume : -std::numeric_limits<Scalar>::infinity();
    }

    template <typename Derived>
    [[nodiscard]] Scalar
    density(Eigen::MatrixBase<Derived> const & x) const
    {
        return Eigen::numext::exp(logDensity(x));
    }

    [[nodiscard]] MomentsType
    moments() const
    {
        auto result = MomentsType{valueDimension(), differenceDimension()};
        result.mean() = m_lower + Scalar(0.5) * m_width;
        normalizeValues(m_space, result.mean());
        result.covarianceLower() = (m_width.array().square() / Scalar(12)).matrix().asDiagonal();
        return result;
    }

    template <typename Path>
    [[nodiscard]] auto
    marginal() const
    {
        return internal::marginalOverPath<Path>(m_space, [&](auto const & target, auto space) {
            return restricted(target.values, std::move(space));
        });
    }

    template <typename Path>
    [[nodiscard]] auto
    marginal(Index element) const
    {
        return internal::marginalOverElement<Path>(m_space, element, [&](auto const & target, auto space) {
            return restricted(target.values, std::move(space));
        });
    }

    template <typename GivenPath, typename Derived>
    [[nodiscard]] auto
    conditional(Eigen::MatrixBase<Derived> const & values) const
    {
        return internal::tagConditional<GivenPath>(
            m_space,
            [&](auto const & t, auto const & g, auto space, auto const &) {
                EMRYS_FILTER_ASSERT(values.size() == g.values.size());
                return restricted(t.values, std::move(space));
            }
        );
    }
 
private:

    template <Space, RealScalar>
    friend class Uniform;

    template <typename Leaf>
    static consteval bool
    hasDefaultInterval()
    {
        constexpr auto support = internal::leafSupport<Leaf>();
        return support && support->lower > -std::numeric_limits<double>::infinity()
            && support->upper < std::numeric_limits<double>::infinity();
    }

    Uniform(Space_ space, VectorType lower, VectorType width, VectorType period)
        : m_space{std::move(space)}
        , m_lower{std::move(lower)}
        , m_width{std::move(width)}
        , m_period{std::move(period)}
        , m_logVolume{m_width.array().log().sum()}
    {}

    template <typename Target, typename SubSpace>
    [[nodiscard]] Uniform<SubSpace, Scalar>
    restricted(Target const & target, SubSpace space) const
    {
        using Result = Uniform<SubSpace, Scalar>;
        using ResultVector = typename Result::VectorType;
        return Result{
            std::move(space),
            ResultVector{m_lower(target)},
            ResultVector{m_width(target)},
            ResultVector{m_period(target)}
        };
    }
 
    EMRYS_FILTER_NO_UNIQUE_ADDRESS Space_ m_space;  
    VectorType m_lower;                             
    VectorType m_width;                             
    VectorType m_period;                            
    Scalar m_logVolume{};                           
};

template <typename Scalar, typename Space, std::uniform_random_bit_generator G>
[[nodiscard]] SampleSet<Space, Scalar>
drawSamples(Uniform<Space, Scalar> const & uniform, Eigen::Index n, G & rng)
{
    EMRYS_FILTER_ASSERT(n >= 0);
    using SampleSetType = SampleSet<Space, Scalar>;
    SampleSetType samples =
        uniform.width().asDiagonal() * internal::standardUniform<SampleSetType>(uniform.valueDimension(), n, rng);
    samples.colwise() += uniform.lower();
    normalizeValues(uniform.space(), samples);
    return samples;
}

template <typename Scalar, typename Space, std::uniform_random_bit_generator G>
[[nodiscard]] typename Moments<Space, Scalar>::VectorType
drawSample(Uniform<Space, Scalar> const & uniform, G & rng)
{
    static_assert(std::floating_point<Scalar>, "sampling requires a standard floating-point scalar");
    using Vector = typename Moments<Space, Scalar>::VectorType;
    auto draw = std::uniform_real_distribution<Scalar>{};
    auto sample = Vector{};
    sample.resize(uniform.valueDimension());
    for (Eigen::Index i = 0; i < sample.size(); ++i)
    {
        sample(i) = uniform.lower()(i) + uniform.width()(i) * draw(rng);
    }
    normalizeValues(uniform.space(), sample);
    return sample;
}

template <typename Scalar, typename Space, std::uniform_random_bit_generator G>
[[nodiscard]] SamplesWithLogDensity<Space, Scalar>
drawSamplesWithLogDensity(Uniform<Space, Scalar> const & uniform, Eigen::Index n, G & rng)
{
    return {drawSamples(uniform, n, rng), Eigen::Array<Scalar, Eigen::Dynamic, 1>::Constant(n, -uniform.logVolume())};
}

template <typename Scalar, typename Space, std::uniform_random_bit_generator G>
[[nodiscard]] SampleWithLogDensity<Space, Scalar>
drawSampleWithLogDensity(Uniform<Space, Scalar> const & uniform, G & rng)
{
    return {drawSample(uniform, rng), -uniform.logVolume()};
}
 
}  // namespace Emrys
 
#ifndef EMRYS_FILTER_NO_IO
 
namespace Emrys
{
 
namespace internal
{

inline void
describe(std::ostream & stream, UniformError error, Indent indent)
{
    stream << indent;
    switch (error)
    {
        case UniformError::Empty:
            stream << "an upper bound is not above its lower bound";
            break;
        case UniformError::NonFinite:
            stream << "a bound is not finite";
            break;
        case UniformError::WiderThanPeriod:
            stream << "an interval is wider than its coordinate's period";
            break;
    }
}
 
}  // namespace internal
 
}  // namespace Emrys
 
#endif  // EMRYS_FILTER_NO_IO