smp_stage_quad.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_QUAD_H
#define LINE_API_MC_SMP_STAGE_QUAD_H
 
/**
 * @file
 * @ingroup api_mc
 * Embedded jump probabilities and mean holding time of a semi-Markov stage.
 *
 * Port of `matlab/src/api/mc/smp_stage_quad.m`. For KERNEL arcs both are in
 * closed form: P = probs and m = sum_h p_h E[F_h]. For COMPETING arcs they are
 * P_h = CIF_h(Inf) and m = int_0^Inf (1 - Fhold(t)) dt (see smp_stage_cif.h),
 * evaluated on [0, U], where U starts at the largest arc mean and doubles
 * until the joint survival is below 1e-12 (or U reaches 1e8 times that mean).
 * The grid is 20001 uniform points on [0, U] merged with 4001 points on
 * [0, min(U, 40 E[F_h])] and the point min(E[F_h], U) per arc, so a fast clock
 * beside a slow one is still resolved. P is renormalized to sum to one, and
 * the mean is Simpson's rule with the midpoint survival from a second
 * smp_stage_cif call.
 *
 * `U` is returned as well: the horizon the quadrature covered.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdio>
#include <string>
#include <vector>
 
#include "line/api/mc/smp_stage_cif.h"
#include "line/lang/distribution.h"
#include "line/lang/lang_types.h"
#include "line/util/error.h"
 
namespace line {
namespace mc {
 
/** What `smp_stage_quad` returns. */
struct SmpStageQuad {
    std::vector<double> P;  ///< the embedded jump probabilities, one per arc
    double m = 0.0;         ///< the mean holding time
    double U = 0.0;         ///< the quadrature horizon
};
 
template <class T>
SmpStageQuad smp_stage_quad(const std::vector<lang::Distrib<T>>& arcs,
                            const std::vector<double>& probs) {
    const std::size_t E = arcs.size();
    SmpStageQuad out;
    out.P.assign(E, 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 (live.empty()) return out;
    std::vector<double> mu(live.size(), 0.0);
    for (std::size_t k = 0; k < live.size(); ++k) {
        mu[k] = num_traits<T>::to_double(arcs[live[k]].mean);
        if (!(std::isfinite(mu[k]) && mu[k] > 0.0)) {
            char buf[64];
            std::snprintf(buf, sizeof buf, "%g", mu[k]);
            throw InputError(std::string("smp_stage_quad: the ") +
                             lang::process_to_text(arcs[live[k]].type) +
                             " law of the transition to stage " + std::to_string(live[k] + 1) +
                             " has mean " + buf +
                             ": a semi-Markov stage needs every transition law to have a finite "
                             "positive mean.");
        }
    }
    if (!probs.empty()) {
        for (std::size_t k = 0; k < live.size(); ++k) {
            out.P[live[k]] = probs[live[k]];
            out.m += probs[live[k]] * mu[k];
        }
        out.U = *std::max_element(mu.begin(), mu.end());
        return out;
    }
    const double mumax = *std::max_element(mu.begin(), mu.end());
    const auto survival = [&](double t) {
        double s = 1.0;
        for (std::size_t h : live)
            s *= 1.0 - num_traits<T>::to_double(
                           lang::dist_cdf(arcs[h], num_traits<T>::from_double(t)));
        return s;
    };
    double U = mumax;
    while (survival(U) > 1e-12 && U < 1e8 * mumax) U = 2.0 * U;
    out.U = U;
    std::vector<double> g = smp_detail::linspace(0.0, U, 20001);
    for (std::size_t k = 0; k < live.size(); ++k) {
        const std::vector<double> gk = smp_detail::linspace(0.0, std::min(U, 40.0 * mu[k]), 4001);
        g.insert(g.end(), gk.begin(), gk.end());
        g.push_back(std::min(mu[k], U));
    }
    std::sort(g.begin(), g.end());
    g.erase(std::unique(g.begin(), g.end()), g.end());
    const SmpStageCif c = smp_stage_cif(arcs, std::vector<double>(), g);
    const std::size_t ng = g.size();
    double tot = 0.0;
    for (std::size_t h = 0; h < E; ++h) {
        out.P[h] = c.C(h, ng - 1);
        tot += out.P[h];
    }
    if (tot > 0.0)
        for (std::size_t h = 0; h < E; ++h) out.P[h] /= tot;
    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;
    const SmpStageCif cm = smp_stage_cif(arcs, std::vector<double>(), gm);
    double acc = 0.0;
    for (std::size_t j = 0; j + 1 < ng; ++j)
        acc += ((1.0 - c.Fhold[j]) + 4.0 * (1.0 - cm.Fhold[j]) + (1.0 - c.Fhold[j + 1])) *
               (g[j + 1] - g[j]);
    out.m = acc / 6.0;
    return out;
}
 
}  // namespace mc
}  // namespace line
 
#endif  // LINE_API_MC_SMP_STAGE_QUAD_H