HiddenMarkovFilter.h

Go to the documentation of this file.

#pragma once

 
#include <expected>
#include <type_traits>
#include <utility>
 
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Core/Math.h>
#include <EmrysFilter/Core/Print.h>
#include <EmrysFilter/Core/Weights.h>
#include <EmrysFilter/Distributions/Categorical.h>
#include <EmrysFilter/Distributions/ContinuousTimeMarkovChain.h>
#include <EmrysFilter/Filters/FilterUpkeep.h>
#include <EmrysFilter/Filters/Recursion.h>
#include <EmrysFilter/Models/Emission.h>
#include <EmrysFilter/Spaces/Discrete.h>
 
namespace Emrys
{
 
namespace internal
{

template <typename Chain>
inline constexpr bool IsContinuousTimeChain = false;

template <typename Space, typename Scalar>
inline constexpr bool IsContinuousTimeChain<ContinuousTimeMarkovChain<Space, Scalar>> = true;
 
}  // namespace internal

template <DiscreteSpace Space, RealScalar Scalar>
struct CategoricalPredictedJoint
{
    using MatrixType = typename MarkovChain<Space, Scalar>::MatrixType;  
 
    CategoricalParams<Space, Scalar> filtered;   
    CategoricalParams<Space, Scalar> predicted;  
    MatrixType transitions;                      
};

template <typename Scalar = double>
class HiddenMarkovFilterUpdate
{
public:

    [[nodiscard]] Scalar
    logLikelihood() const
    {
        return m_logLikelihood;
    }
 
    bool accepted = false;  
 
private:
 
    friend struct internal::UpdateAccess;
 
    Scalar m_logLikelihood = 0;  
};

template <
    DiscreteSpace Space_,
    typename Chain_ = MarkovChain<Space_>,
    RealScalar Scalar_ = double,
    typename Upkeep_ = upkeep::None>
class HiddenMarkovFilter : public internal::HeldUpkeep<HiddenMarkovFilter<Space_, Chain_, Scalar_, Upkeep_>, Upkeep_>
{
    using UpkeepBase = internal::HeldUpkeep<HiddenMarkovFilter, Upkeep_>;  
 
    static_assert(
        std::is_same_v<typename Chain_::SpaceType, Space_> && std::is_same_v<typename Chain_::Scalar, Scalar_>,
        "the chain must be on the filter's space, in its scalar"
    );
 
public:
 
    using Scalar = Scalar_;                                
    using SpaceType = Space_;                              
    using ChainType = Chain_;                              
    using StateType = CategoricalParams<Space_, Scalar_>;  

    HiddenMarkovFilter(Chain_ chain, StateType state, std::type_identity_t<Upkeep_> upkeep = {})
        : UpkeepBase{std::move(upkeep)}
        , m_chain{std::move(chain)}
        , m_state{std::move(state)}
    {
        EMRYS_FILTER_ASSERT(m_chain.space() == m_state.space() && "the chain is on a different space");
    }

    explicit HiddenMarkovFilter(Chain_ chain)
        : UpkeepBase{Upkeep_{}}
        , m_chain{std::move(chain)}
        , m_state{m_chain.space()}
    {}

    [[nodiscard]] StateType const &
    state() const
    {
        return m_state;
    }

    [[nodiscard]] StateType &
    state()
    {
        return m_state;
    }

    void
    setState(StateType state)
    {
        EMRYS_FILTER_ASSERT(m_chain.space() == state.space() && "the state is on a different space");
        m_state = std::move(state);
    }

    [[nodiscard]] Chain_ const &
    chain() const
    {
        return m_chain;
    }

    [[nodiscard]] Chain_ &
    chain()
    {
        return m_chain;
    }

    [[nodiscard]] Scalar
    logEvidence() const
    {
        return m_logEvidence;
    }

    void
    resetLogEvidence()
    {
        m_logEvidence = Scalar(0);
    }

    void
    predict()
        requires(!internal::IsContinuousTimeChain<Chain_>)
    {
        m_state = Emrys::predict(m_state, m_chain);
    }

    void
    predict(TimeStep<Scalar> auto const & dt)
        requires internal::IsContinuousTimeChain<Chain_>
    {
        m_state = Emrys::predict(m_state, m_chain, dt);
    }

    [[nodiscard]] std::expected<Scalar, WeightError>
    predict(TimeStep<Scalar> auto const & dt, ArrivalInterval<Space_, Scalar_> const & interval)
        requires internal::IsContinuousTimeChain<Chain_>
    {
        return apply(Emrys::predict(m_state, m_chain, dt, interval));
    }

    [[nodiscard]] CategoricalPredictedJoint<Space_, Scalar_>
    predictJoint()
        requires(!internal::IsContinuousTimeChain<Chain_>)
    {
        auto predicted = Emrys::predict(m_state, m_chain);
        auto joint = CategoricalPredictedJoint<Space_, Scalar_>{std::move(m_state), predicted, m_chain.matrix()};
        m_state = std::move(predicted);
        return joint;
    }

    [[nodiscard]] CategoricalPredictedJoint<Space_, Scalar_>
    predictJoint(TimeStep<Scalar> auto const & dt)
        requires internal::IsContinuousTimeChain<Chain_>
    {
        auto step = m_chain.transition(dt);
        auto predicted = Emrys::predict(m_state, step);
        auto joint = CategoricalPredictedJoint<Space_, Scalar_>{std::move(m_state), predicted, step.matrix()};
        m_state = std::move(predicted);
        return joint;
    }

    [[nodiscard]] std::expected<CategoricalPredictedJoint<Space_, Scalar_>, WeightError>
    predictJoint(TimeStep<Scalar> auto const & dt, ArrivalInterval<Space_, Scalar_> const & interval)
        requires internal::IsContinuousTimeChain<Chain_>
    {
        auto step = internal::arrivalKernel(m_chain, internal::toScalarTime<Scalar>(dt), interval);
        auto moved = internal::moveByKernel(m_state, step.kernel, step.logScale);
        if (!moved)
        {
            return std::unexpected{moved.error()};
        }
        // Divide the kernel by the step's total, so that it moves the filtered state onto the
        // normalized prediction exactly.
        step.kernel /= Eigen::numext::exp(Scalar(moved->logEvidence - step.logScale));
        auto joint =
            CategoricalPredictedJoint<Space_, Scalar_>{std::move(m_state), moved->posterior, std::move(step.kernel)};
        m_state = std::move(moved->posterior);
        m_logEvidence += moved->logEvidence;
        return joint;
    }

    template <EmissionModel<Space_> E>
    [[nodiscard]] std::expected<HiddenMarkovFilterUpdate<Scalar_>, WeightError>
    update(E const & emission)
    {
        auto const prior = m_state.logTotalWeight();
        if (!prior)
        {
            return std::unexpected{prior.error()};
        }
        auto const evidence = m_state.update(internal::emissionLogLikelihoods(m_state, emission));
        auto result = HiddenMarkovFilterUpdate<Scalar_>{};
        if (!evidence)
        {
            // The prior is valid, so a sum of zero is the measurement's: no value explains it.
            if (evidence.error() != WeightError::ZeroSum)
            {
                return std::unexpected{evidence.error()};
            }
            internal::UpdateAccess::logLikelihood(result) = -Eigen::NumTraits<Scalar_>::infinity();
            m_logEvidence = -Eigen::NumTraits<Scalar_>::infinity();
            return result;
        }
        internal::UpdateAccess::logLikelihood(result) = *evidence;
        result.accepted = true;
        m_logEvidence += *evidence;
        return result;
    }
 
private:

    [[nodiscard]] std::expected<Scalar, WeightError>
    apply(std::expected<CategoricalUpdate<Space_, Scalar_>, WeightError> result)
    {
        if (!result)
        {
            return std::unexpected{result.error()};
        }
        m_state = std::move(result->posterior);
        m_logEvidence += result->logEvidence;
        return result->logEvidence;
    }
 
    Chain_ m_chain;            
    StateType m_state;         
    Scalar m_logEvidence = 0;  
};

template <typename Chain>
HiddenMarkovFilter(Chain, CategoricalParams<typename Chain::SpaceType, typename Chain::Scalar>)
    -> HiddenMarkovFilter<typename Chain::SpaceType, Chain, typename Chain::Scalar>;
 
template <typename Chain>
HiddenMarkovFilter(Chain) -> HiddenMarkovFilter<typename Chain::SpaceType, Chain, typename Chain::Scalar>;
 
template <typename Chain, typename U>
HiddenMarkovFilter(Chain, CategoricalParams<typename Chain::SpaceType, typename Chain::Scalar>, U)
    -> HiddenMarkovFilter<typename Chain::SpaceType, Chain, typename Chain::Scalar, U>;

 
}  // namespace Emrys
 
#ifndef EMRYS_FILTER_NO_IO
 
namespace Emrys
{
 
namespace internal
{

template <typename Space_, typename Chain_, typename Scalar_, typename Upkeep_>
void
describe(std::ostream & stream, HiddenMarkovFilter<Space_, Chain_, Scalar_, Upkeep_> const & filter, Indent indent)
{
    describeFilter(stream, "hidden Markov filter", filter.state(), indent);
}
 
}  // namespace internal
 
}  // namespace Emrys
 
#endif  // EMRYS_FILTER_NO_IO