CategoricalHistory.h

Go to the documentation of this file.

#pragma once

 
#include <cstddef>
#include <iterator>
#include <list>
#include <optional>
#include <utility>
 
#include <EmrysFilter/Core/Allocation.h>
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Core/Math.h>
#include <EmrysFilter/Core/Multiply.h>
#include <EmrysFilter/Core/Weights.h>
#include <EmrysFilter/Distributions/Categorical.h>
#include <EmrysFilter/Filters/HiddenMarkovFilter.h>
#include <EmrysFilter/Smoothers/Lag.h>
#include <EmrysFilter/Spaces/Discrete.h>
#include <EmrysFilter/Spaces/Space.h>
 
namespace Emrys
{
 
namespace internal
{

template <typename Space, typename Scalar>
[[nodiscard]] Scalar
smoothingRatioScale(
    CategoricalPredictedJoint<Space, Scalar> const & joint,
    CategoricalParams<Space, Scalar> const & smoothedNext
)
{
    EMRYS_FILTER_ASSERT(
        smoothedNext.logTotalWeight() && "the smoothed distribution at the step's end must have valid weights"
    );
    Scalar const minusInfinity = -Eigen::NumTraits<Scalar>::infinity();
    auto const logNext = logProbabilitiesOf(smoothedNext);
    auto const logPredicted = logProbabilitiesOf(joint.predicted);
    Scalar scale = minusInfinity;
    for (Eigen::Index j = 0; j < logNext.size(); ++j)
    {
        if (logNext(j) != minusInfinity)
        {
            EMRYS_FILTER_ASSERT(
                logPredicted(j) != minusInfinity && "a value the prediction rules out has smoothed probability"
            );
            scale = Eigen::numext::maxi(scale, Scalar(logNext(j) - logPredicted(j)));
        }
    }
    return scale;
}

template <typename Space, typename Scalar>
[[nodiscard]] typename CategoricalParams<Space, Scalar>::ArrayType
smoothingRatios(
    CategoricalPredictedJoint<Space, Scalar> const & joint,
    CategoricalParams<Space, Scalar> const & smoothedNext
)
{
    EMRYS_FILTER_ASSERT(joint.filtered.space() == smoothedNext.space() && "the distributions are on different spaces");
    Scalar const minusInfinity = -Eigen::NumTraits<Scalar>::infinity();
    Scalar const scale = smoothingRatioScale(joint, smoothedNext);
    auto const logNext = logProbabilitiesOf(smoothedNext);
    auto const logPredicted = logProbabilitiesOf(joint.predicted);
    auto ratios = typename CategoricalParams<Space, Scalar>::ArrayType{};
    ratios.resize(logNext.size());
    for (Eigen::Index j = 0; j < logNext.size(); ++j)
    {
        ratios(j) =
            logNext(j) == minusInfinity ? Scalar(0) : Eigen::numext::exp(Scalar(logNext(j) - logPredicted(j) - scale));
    }
    return ratios;
}
 
}  // namespace internal

template <typename Space, typename Scalar>
[[nodiscard]] CategoricalParams<Space, Scalar>
smoothStep(
    CategoricalPredictedJoint<Space, Scalar> const & joint,
    CategoricalParams<Space, Scalar> const & smoothedNext
)
{
    using ArrayType = typename CategoricalParams<Space, Scalar>::ArrayType;
    ArrayType const ratios = internal::smoothingRatios(joint, smoothedNext);
    ArrayType const logSmoothed = internal::logProbabilitiesOf(joint.filtered)
                                + internal::multiply(joint.transitions, ratios.matrix()).array().log();
    auto smoothed = CategoricalParams<Space, Scalar>{joint.filtered.space(), logSmoothed};
    // A value with probability at the end was predicted from one with probability at the start.
    [[maybe_unused]] auto const normalized = smoothed.normalize();
    EMRYS_FILTER_ASSERT(normalized && "no value smoothed at the step's end came from a value at its start");
    return smoothed;
}

template <typename Space, typename Scalar>
[[nodiscard]] typename CategoricalPredictedJoint<Space, Scalar>::MatrixType
transitionPosterior(
    CategoricalPredictedJoint<Space, Scalar> const & joint,
    CategoricalParams<Space, Scalar> const & smoothedNext
)
{
    using ArrayType = typename CategoricalParams<Space, Scalar>::ArrayType;
    ArrayType const ratios = internal::smoothingRatios(joint, smoothedNext);
    ArrayType const filtered = internal::logProbabilitiesOf(joint.filtered).exp();
    typename CategoricalPredictedJoint<Space, Scalar>::MatrixType posterior =
        filtered.matrix().asDiagonal() * joint.transitions * ratios.matrix().asDiagonal();
    posterior /= posterior.sum();
    return posterior;
}

template <DiscreteSpace Space_, RealScalar Scalar_ = double, typename Allocator_ = DefaultAllocator>
class CategoricalHistory
{
public:
 
    using Scalar = Scalar_;                                        
    using SpaceType = Space_;                                      
    using CategoricalType = CategoricalParams<Space_, Scalar_>;    
    using JointType = CategoricalPredictedJoint<Space_, Scalar_>;  
    using allocator_type = Allocator_;                             
    using Index = Eigen::Index;                                    

    struct Step
    {
        JointType joint;  

        CategoricalType smoothed;
    };
 
private:
 
    using StepAllocator = internal::ReBound<Allocator_, Step>;  
    using ListType = std::list<Step, StepAllocator>;            
 
public:
 
    using const_iterator = typename ListType::const_iterator;  

    CategoricalHistory() = default;

    explicit CategoricalHistory(Allocator_ const & allocator)
        : m_steps(StepAllocator(allocator))
    {}

    explicit CategoricalHistory(Lag lag, Allocator_ const & allocator = {})
        : m_steps(StepAllocator(allocator))
        , m_lag{lag.steps}
    {
        EMRYS_FILTER_ASSERT(lag.steps > 0 && "a fixed-lag history keeps at least one step");
    }

    CategoricalHistory(CategoricalHistory const & other)
        : CategoricalHistory{other, other.get_allocator()}
    {}

    CategoricalHistory(CategoricalHistory const & other, Allocator_ const & allocator)
        : m_steps(other.m_steps, StepAllocator(allocator))
        , m_lag{other.m_lag}
    {}

    CategoricalHistory(CategoricalHistory && other) = default;

    CategoricalHistory(CategoricalHistory && other, Allocator_ const & allocator)
        : m_steps(std::move(other.m_steps), StepAllocator(allocator))
        , m_lag{other.m_lag}
    {}

    CategoricalHistory &
    operator=(CategoricalHistory const & other) = default;

    CategoricalHistory &
    operator=(CategoricalHistory && other) = default;
 
    ~CategoricalHistory() = default;

    [[nodiscard]] allocator_type
    get_allocator() const
    {
        return allocator_type(m_steps.get_allocator());
    }

    [[nodiscard]] std::optional<Index>
    lag() const
    {
        return m_lag > 0 ? std::optional<Index>{m_lag} : std::nullopt;
    }

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

    [[nodiscard]] bool
    empty() const
    {
        return m_steps.empty();
    }

    [[nodiscard]] Step const &
    front() const
    {
        EMRYS_FILTER_ASSERT(!empty() && "the history is empty");
        return m_steps.front();
    }

    [[nodiscard]] Step const &
    back() const
    {
        EMRYS_FILTER_ASSERT(!empty() && "the history is empty");
        return m_steps.back();
    }

    [[nodiscard]] const_iterator
    begin() const
    {
        return m_steps.begin();
    }

    [[nodiscard]] const_iterator
    end() const
    {
        return m_steps.end();
    }

    void
    push(JointType joint)
    {
        if (m_lag > 0 && size() == m_lag)
        {
            m_steps.splice(m_steps.end(), m_steps, m_steps.begin());
            Step & step = m_steps.back();
            step.joint = std::move(joint);
            step.smoothed = step.joint.filtered;
        }
        else
        {
            CategoricalType smoothed = joint.filtered;
            m_steps.push_back(Step{std::move(joint), std::move(smoothed)});
        }
    }

    void
    smooth(CategoricalType const & latest)
    {
        CategoricalType const * next{&latest};
        for (auto step = m_steps.rbegin(); step != m_steps.rend(); ++step)
        {
            step->smoothed = smoothStep(step->joint, *next);
            next = &step->smoothed;
        }
    }

    void
    clear()
    {
        m_steps.clear();
    }
 
private:
 
    ListType m_steps;  
    Index m_lag = 0;   
};

template <DiscreteSpace Space, RealScalar Scalar, typename Allocator = DefaultAllocator>
struct ValueSequence
{
    internal::Buffer<Eigen::Index, Allocator> values;

    Scalar logProbability;
};

template <typename Space, typename Scalar, typename Allocator>
[[nodiscard]] ValueSequence<Space, Scalar, Allocator>
mostProbableSequence(
    CategoricalHistory<Space, Scalar, Allocator> const & history,
    CategoricalParams<Space, Scalar> const & latest
)
{
    using ArrayType = typename CategoricalParams<Space, Scalar>::ArrayType;
    Scalar const minusInfinity = -Eigen::NumTraits<Scalar>::infinity();
    Eigen::Index const count = latest.cardinality();
    std::size_t const steps = std::size_t(history.size());
    auto backPointers = internal::bufferOf<Eigen::Index>(history.get_allocator());
    backPointers.resize(steps * std::size_t(count));
 
    ArrayType delta = internal::logProbabilitiesOf(history.empty() ? latest : history.front().joint.filtered);
    auto next = ArrayType{};
    next.resize(count);
    std::size_t k = 0;
    for (auto step = history.begin(); step != history.end(); ++step, ++k)
    {
        auto const following = std::next(step);
        auto const filtered =
            internal::logProbabilitiesOf(following == history.end() ? latest : following->joint.filtered);
        auto const predicted = internal::logProbabilitiesOf(step->joint.predicted);
        for (Eigen::Index j = 0; j < count; ++j)
        {
            Scalar best = minusInfinity;
            Eigen::Index from = 0;
            for (Eigen::Index i = 0; i < count; ++i)
            {
                Scalar const candidate = delta(i) + Eigen::numext::log(step->joint.transitions(i, j));
                if (candidate > best)
                {
                    best = candidate;
                    from = i;
                }
            }
            backPointers[k * std::size_t(count) + std::size_t(j)] = from;
            bool const reachable = best != minusInfinity && predicted(j) != minusInfinity;
            next(j) = reachable ? Scalar(best + filtered(j) - predicted(j)) : minusInfinity;
        }
        delta.swap(next);
    }
 
    auto result = ValueSequence<Space, Scalar, Allocator>{
        internal::bufferOf<Eigen::Index>(history.get_allocator()),
        minusInfinity
    };
    result.values.resize(steps + 1);
    Eigen::Index value = 0;
    result.logProbability = delta.maxCoeff(&value);
    result.values[steps] = value;
    for (std::size_t s = steps; s > 0; --s)
    {
        value = backPointers[(s - 1) * std::size_t(count) + std::size_t(value)];
        result.values[s - 1] = value;
    }
    return result;
}
 
}  // namespace Emrys