smp_stage_cif.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_MC_SMP_STAGE_CIF_H
#define LINE_API_MC_SMP_STAGE_CIF_H
 
/**
 * @file
 * @ingroup api_mc
 * Cumulative incidence functions of one semi-Markov stage.
 *
 * Port of `matlab/src/api/mc/smp_stage_cif.m`. A stage left by E possible
 * transitions is described in one of two ways, which are the two
 * parameterizations of a semi-Markov kernel compared by Asanjarani, Liquet and
 * Nazarathy (arXiv:2005.14462, 2020):
 *
 *   COMPETING arcs (transition intensities, `probs` empty): every arc h is a
 *   latent clock with CDF F_h started on entry, the first to ring picks the
 *   destination, and
 *       CIF_h(t) = int_0^t prod_{k!=h} (1 - F_k(s)) dF_h(s),
 *       Fhold(t) = 1 - prod_k (1 - F_k(t)).
 *
 *   KERNEL arcs (sojourn times, `probs` the jump probabilities): the
 *   destination is drawn on entry and the sojourn then has the law F_h, so
 *       CIF_h(t) = p_h F_h(t),   Fhold(t) = sum_h p_h F_h(t).
 *
 * The competing integral is the Riemann-Stieltjes sum on the grid `t` with the
 * other clocks' survival read at the interval midpoints, which needs no
 * density and so admits a deterministic clock. When `t(1) > 0` the sum starts
 * at zero on 200 uniform intervals below `t(1)`, whose values are dropped.
 * Every CDF is clamped to [0, 1] before use.
 */
 
#include <algorithm>
#include <cstddef>
#include <vector>
 
#include "line/lang/distribution.h"
#include "line/lang/lang_types.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace mc {
 
/** What `smp_stage_cif` returns: `C(h, j)` = CIF_h(t_j) and `Fhold[j]`. */
struct SmpStageCif {
    Matrix<double> C;
    std::vector<double> Fhold;
};
 
namespace smp_detail {
 
/** `linspace(a, b, n)`, with the last point exactly `b`. */
inline std::vector<double> linspace(double a, double b, std::size_t n) {
    std::vector<double> y(n, a);
    if (n == 0) return y;
    if (n == 1) {
        y[0] = b;
        return y;
    }
    const double n1 = static_cast<double>(n - 1);
    for (std::size_t k = 0; k < n; ++k) y[k] = a + static_cast<double>(k) * (b - a) / n1;
    y[n - 1] = b;
    return y;
}
 
/** The CDF of `d` at every instant of `t`, clamped to [0, 1] (`cdfRow`). */
template <class T>
std::vector<double> cdf_row(const lang::Distrib<T>& d, const std::vector<double>& t) {
    std::vector<double> F(t.size(), 0.0);
    for (std::size_t j = 0; j < t.size(); ++j) {
        const double v =
            num_traits<T>::to_double(lang::dist_cdf(d, num_traits<T>::from_double(t[j])));
        F[j] = std::min(std::max(v, 0.0), 1.0);
    }
    return F;
}
 
}  // namespace smp_detail
 
/**
 * CIF and holding-time CDF of a stage at the ascending instants `t`.
 *
 * `arcs[h]` is the law of the transition to stage h, DISABLED where there is
 * none; `probs` is empty for competing arcs, or holds the E kernel
 * probabilities.
 */
template <class T>
SmpStageCif smp_stage_cif(const std::vector<lang::Distrib<T>>& arcs,
                          const std::vector<double>& probs, const std::vector<double>& t) {
    const std::size_t E = arcs.size(), nt = t.size();
    SmpStageCif out;
    out.C = Matrix<double>(E, nt, 0.0);
    out.Fhold.assign(nt, 0.0);
    std::vector<std::size_t> live;
    for (std::size_t h = 0; h < E; ++h)
        if (!arcs[h].disabled) live.push_back(h);
    if (nt == 0 || live.empty()) return out;
    if (!probs.empty()) {
        if (probs.size() != E)
            throw InputError("smp_stage_cif: the kernel probabilities need one entry per arc");
        for (std::size_t h : live) {
            const std::vector<double> Fh = smp_detail::cdf_row(arcs[h], t);
            for (std::size_t j = 0; j < nt; ++j) {
                out.C(h, j) = probs[h] * Fh[j];
                out.Fhold[j] += probs[h] * Fh[j];
            }
        }
        return out;
    }
    std::size_t npre = 0;
    std::vector<double> g;
    if (t[0] > 0.0) {
        const std::vector<double> pre = smp_detail::linspace(0.0, t[0], 201);
        g.assign(pre.begin(), pre.end() - 1);
        npre = 200;
    }
    g.insert(g.end(), t.begin(), t.end());
    const std::size_t ng = g.size(), nl = live.size();
    std::vector<double> gm(ng - 1);
    for (std::size_t j = 0; j + 1 < ng; ++j) gm[j] = (g[j] + g[j + 1]) / 2.0;
    std::vector<std::vector<double>> F(nl), Fm(nl);
    for (std::size_t k = 0; k < nl; ++k) {
        F[k] = smp_detail::cdf_row(arcs[live[k]], g);
        Fm[k] = smp_detail::cdf_row(arcs[live[k]], gm);
    }
    for (std::size_t k = 0; k < nl; ++k) {
        double acc = 0.0;
        for (std::size_t j = 1; j < ng; ++j) {
            double Sother = 1.0;
            for (std::size_t q = 0; q < nl; ++q)
                if (q != k) Sother *= 1.0 - Fm[q][j - 1];
            acc += (F[k][j] - F[k][j - 1]) * Sother;
            if (j >= npre) out.C(live[k], j - npre) = acc;
        }
    }
    for (std::size_t j = 0; j < nt; ++j) {
        double S = 1.0;
        for (std::size_t k = 0; k < nl; ++k) S *= 1.0 - F[k][npre + j];
        out.Fhold[j] = 1.0 - S;
    }
    return out;
}
 
}  // namespace mc
}  // namespace line
 
#endif  // LINE_API_MC_SMP_STAGE_CIF_H