ContinuousTimeMarkovChain.h

Go to the documentation of this file.

#pragma once

 
#include <tuple>
#include <type_traits>
#include <utility>
 
#include <unsupported/Eigen/MatrixFunctions>
 
#include <EmrysFilter/Core/Forwarding.h>
#include <EmrysFilter/Core/Macros.h>
#include <EmrysFilter/Core/Math.h>
#include <EmrysFilter/Core/Time.h>
#include <EmrysFilter/Distributions/Categorical.h>
#include <EmrysFilter/Spaces/Discrete.h>
#include <EmrysFilter/Spaces/Space.h>
 
namespace Emrys
{

template <DiscreteSpace Space_, RealScalar Scalar_ = double>
class ContinuousTimeMarkovChain
{
public:
 
    using Scalar = Scalar_;                                               
    using SpaceType = Space_;                                             
    static constexpr int CardinalityAtCompileTime = Space_::Cardinality;  
    static constexpr int MaxCardinalityAtCompileTime = Space_::MaxCardinality;  

    using MatrixType = typename MarkovChain<Space_, Scalar_>::MatrixType;

    explicit ContinuousTimeMarkovChain(Space_ space = {})
        : m_space{std::move(space)}
        , m_generator{MatrixType::Zero(m_space.cardinality(), m_space.cardinality())}
    {}

    template <MatrixExpression Generator>
    ContinuousTimeMarkovChain(Space_ space, Generator && generator)
        : m_space{std::move(space)}
        , m_generator{internal::castScalar<Scalar>(std::forward<Generator>(generator))}
    {
        EMRYS_FILTER_ASSERT(
            m_generator.rows() == m_space.cardinality() && m_generator.cols() == m_space.cardinality()
            && "the generator has a row and a column per value"
        );
        [[maybe_unused]] Scalar scale = 0;
        for (Eigen::Index i = 0; i < m_generator.rows(); ++i)
        {
            for (Eigen::Index j = 0; j < m_generator.cols(); ++j)
            {
                EMRYS_FILTER_ASSERT((i == j || m_generator(i, j) >= Scalar(0)) && "rates must be nonnegative");
            }
            scale = Eigen::numext::maxi(scale, Scalar(-m_generator(i, i)));
        }
        EMRYS_FILTER_ASSERT(
            (m_generator.rowwise().sum().array().abs()
             <= Eigen::NumTraits<Scalar>::dummy_precision() * Eigen::numext::maxi(scale, Scalar(1)))
                .all()
            && "each row of the generator must sum to zero"
        );
    }

    template <MatrixExpression Generator>
        requires FixedSizeSpace<Space_>
    explicit ContinuousTimeMarkovChain(Generator && generator)
        : ContinuousTimeMarkovChain{Space_{}, std::forward<Generator>(generator)}
    {}

    template <DiscreteSpace... Axes>
        requires(
            FixedSizeSpace<Space_> && internal::IsDiscreteProduct<Space_>
            && sizeof...(Axes) == std::tuple_size_v<internal::DiscreteProductAxes<Space_>>
        )
    explicit ContinuousTimeMarkovChain(ContinuousTimeMarkovChain<Axes, Scalar_> const &... chains)
        : ContinuousTimeMarkovChain{Space_{}, chains...}
    {}

    template <DiscreteSpace... Axes>
        requires(
            internal::IsDiscreteProduct<Space_>
            && sizeof...(Axes) == std::tuple_size_v<internal::DiscreteProductAxes<Space_>>
        )
    ContinuousTimeMarkovChain(Space_ space, ContinuousTimeMarkovChain<Axes, Scalar_> const &... chains)
        : m_space{std::move(space)}
    {
        static_assert(
            std::is_same_v<std::tuple<Axes...>, internal::DiscreteProductAxes<Space_>>,
            "give one chain per axis of the product, in order"
        );
        Eigen::Index const n = m_space.cardinality();
        m_generator = MatrixType::Zero(n, n);
        Eigen::Index stride = n;
        (
            [&] {
                Eigen::Index const count = chains.space().cardinality();
                stride /= count;
                // The axis jumps from i to j while every other coordinate stays: i and j differ
                // only in this axis's coordinate, which a multiple of the stride carries.
                for (Eigen::Index i = 0; i < n; ++i)
                {
                    Eigen::Index const from = (i / stride) % count;
                    Eigen::Index const base = i - from * stride;
                    for (Eigen::Index to = 0; to < count; ++to)
                    {
                        m_generator(i, base + to * stride) += chains.generator()(from, to);
                    }
                }
            }(),
            ...
        );
        EMRYS_FILTER_ASSERT(stride == 1 && "each chain must have its axis's cardinality");
    }

    [[nodiscard]] static ContinuousTimeMarkovChain
    persistent(Scalar rate, Space_ space = {})
    {
        EMRYS_FILTER_ASSERT(rate >= Scalar(0) && "the rate must be nonnegative");
        Eigen::Index const n = space.cardinality();
        auto result = ContinuousTimeMarkovChain{std::move(space)};
        if (n > 1)
        {
            result.m_generator.setConstant(rate / Scalar(n - 1));
            result.m_generator.diagonal().setConstant(-rate);
        }
        return result;
    }

    template <typename DwellTimes, typename ExitWeights>
    [[nodiscard]] static ContinuousTimeMarkovChain
    fromDwellTimes(
        Eigen::MatrixBase<DwellTimes> const & dwellTimes,
        Eigen::MatrixBase<ExitWeights> const & exitWeights,
        Space_ space = {}
    )
    {
        Eigen::Index const n = space.cardinality();
        MatrixType generator = internal::dwellGenerator<Scalar, MatrixType>(dwellTimes, exitWeights);
        EMRYS_FILTER_ASSERT(generator.rows() == n && "give a dwell time per value");
        return ContinuousTimeMarkovChain{std::move(space), std::move(generator)};
    }

    template <typename DwellTimes>
    [[nodiscard]] static ContinuousTimeMarkovChain
    fromDwellTimes(Eigen::MatrixBase<DwellTimes> const & dwellTimes, Space_ space = {})
    {
        Eigen::Index const n = space.cardinality();
        return fromDwellTimes(dwellTimes, MatrixType::Ones(n, n) - MatrixType::Identity(n, n), std::move(space));
    }

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

    [[nodiscard]] MatrixType const &
    generator() const
    {
        return m_generator;
    }

    [[nodiscard]] Scalar
    rate(Eigen::Index from, Eigen::Index to) const
    {
        return m_generator(from, to);
    }

    [[nodiscard]] MarkovChain<Space_, Scalar_>
    transition(TimeStep<Scalar_> auto const & time) const
    {
        Scalar const dt = internal::toScalarTime<Scalar_>(time);
        MatrixType const scaled = m_generator * dt;
        MatrixType transitions = scaled.exp();
        // The exponential of a generator is stochastic up to rounding: clear the negative
        // rounding off the zeros and rescale each row to sum to one.
        transitions = transitions.cwiseMax(Scalar(0));
        transitions.array().colwise() /= transitions.rowwise().sum().array();
        return MarkovChain<Space_, Scalar_>{m_space, std::move(transitions)};
    }
 
private:
 
    Space_ m_space;          
    MatrixType m_generator;  
};

template <DiscreteSpace Space, MatrixExpression Generator>
ContinuousTimeMarkovChain(Space, Generator &&)
    -> ContinuousTimeMarkovChain<Space, typename std::remove_cvref_t<Generator>::Scalar>;


template <typename Scalar, typename Space>
[[nodiscard]] CategoricalParams<Space, Scalar>
predict(
    CategoricalParams<Space, Scalar> const & prior,
    ContinuousTimeMarkovChain<Space, Scalar> const & chain,
    TimeStep<Scalar> auto const & dt
)
{
    return predict(prior, chain.transition(dt));
}

template <typename Scalar, typename Space>
[[nodiscard]] std::expected<CategoricalParams<Space, Scalar>, ChainError>
stationary(ContinuousTimeMarkovChain<Space, Scalar> const & chain)
{
    using ArrayType = typename CategoricalParams<Space, Scalar>::ArrayType;
    std::expected<ArrayType, ChainError> const probabilities =
        internal::stationaryProbabilities<Scalar, ArrayType>(chain.generator());
    if (!probabilities)
    {
        return std::unexpected{probabilities.error()};
    }
    return CategoricalParams<Space, Scalar>{chain.space(), probabilities->log()};
}
 
}  // namespace Emrys