infer_smp.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_INFER_INFER_SMP_H
#define LINE_API_INFER_INFER_SMP_H
 
/**
 * @file
 * @ingroup api_infer
 * Maximum-likelihood fit of a parametric semi-Markov multi-state model.
 *
 * Port of `matlab/src/api/infer/infer_smp.m` and its four companions,
 * `infer_smp_family.m`, `infer_smp_env.m`, `infer_smp_hazard.m` and
 * `infer_smp_nelson_aalen.m`, which this header carries together.
 *
 * DATA. One row per observed sojourn: the process entered state `from[k]`,
 * stayed `dur[k] > 0`, and then moved to state `to[k]`, or was still in
 * `from[k]` when observation stopped, coded `to[k] = 0` (right censoring).
 * States are 1..E, as in the reference, since 0 is the censoring code.
 *
 * THE TWO APPROACHES of Asanjarani, Liquet and Nazarathy (arXiv:2005.14462,
 * 2020). `intensity` (Approach II, the default) fits each transition i -> j as
 * a latent clock with its own law, so the likelihood splits into one censored
 * two-state likelihood per transition, every sojourn out of i that did not end
 * in j counting as censored for i -> j. `sojourn` (Approach I) fits the
 * embedded chain p_ij and the law F_ij of a sojourn that ends in j; a censored
 * sojourn contributes S_i(U) = sum_j p_ij S_ij(U), which couples the
 * parameters of a state, so a state with censored rows is fitted jointly (p by
 * a softmax over the logits of destinations 2..m) and one without in closed
 * form for p.
 *
 * FAMILIES, in unconstrained coordinates THETA (`infer_smp_family_*`):
 *   exp        THETA = log(rate)
 *   weibull    THETA = [log(shape), log(scale)]
 *   gamma      THETA = [log(shape), log(scale)]
 *   lognormal  THETA = [mu, log(sigma)]
 * Covariates enter as a proportional-hazards factor exp(beta' z).
 *
 * THE MAXIMIZER is the reference's deterministic Nelder-Mead (the Lagarias et
 * al. variant of fminsearch, restarted from its own optimum until it stops
 * improving), step for step, so the four codebases walk the same simplex; an
 * exponential transition without covariates has its closed-form MLE. Standard
 * errors come from the central-difference Hessian of the negative
 * log-likelihood in THETA coordinates, NaN when it is not positive definite.
 */
 
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstddef>
#include <functional>
#include <limits>
#include <string>
#include <vector>
 
#include "line/api/mc/smp_stage_cif.h"
#include "line/api/mc/smp_stage_quad.h"
#include "line/api/pfqn/pfqn_sjn.h"
#include "line/lang/lang_types.h"
#include "line/lang/qn/environment.h"
#include "line/lang/qn/network_struct.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace infer {
 
// ---- infer_smp_family -------------------------------------------------------
 
/** Lower-case copy of a family name. */
inline std::string smp_family_name(const std::string& f) {
    std::string s = f;
    std::transform(s.begin(), s.end(), s.begin(),
                   [](unsigned char c) { return static_cast<char>(std::tolower(c)); });
    return s;
}
 
/** `infer_smp_family(F, 'npar')`: the number of parameters, validating the name. */
inline std::size_t infer_smp_family_npar(const std::string& family) {
    const std::string f = smp_family_name(family);
    if (f == "exp") return 1;
    if (f == "weibull" || f == "gamma" || f == "lognormal") return 2;
    throw InputError("infer_smp_family: Unknown family '" + f +
                     "': use 'exp', 'weibull', 'gamma' or 'lognormal'.");
}
 
/**
 * `infer_smp_family(F, 'logfs', THETA, T, ETA)`: log density and log survival
 * at `t` under the proportional-hazards factor exp(eta): H = H0 exp(eta),
 * log h = log h0 + eta, log f = log h - H. `eta` holds one value per instant,
 * or a single value for all of them, or nothing for the baseline. The survival
 * of 'gamma' and 'lognormal' is floored at 1e-300 before the log.
 */
inline void infer_smp_family_logfs(const std::string& family, const std::vector<double>& th,
                                   const std::vector<double>& t, const std::vector<double>& eta,
                                   std::vector<double>& logf, std::vector<double>& logS) {
    const std::string f = smp_family_name(family);
    const std::size_t n = t.size();
    logf.assign(n, 0.0);
    logS.assign(n, 0.0);
    if (!(eta.empty() || eta.size() == 1 || eta.size() == n))
        throw InputError("infer_smp_family: ETA must be a scalar or hold one value per instant");
    const double log2pi = std::log(2.0 * 3.14159265358979323846);
    for (std::size_t k = 0; k < n; ++k) {
        const double tk = t[k];
        double H0 = 0.0, logh0 = 0.0;
        if (f == "exp") {
            const double lam = std::exp(th[0]);
            logh0 = std::log(lam);
            H0 = lam * tk;
        } else if (f == "weibull") {
            const double kk = std::exp(th[0]), s = std::exp(th[1]);
            H0 = std::pow(tk / s, kk);
            logh0 = std::log(kk) - std::log(s) + (kk - 1.0) * (std::log(tk) - std::log(s));
        } else if (f == "gamma") {
            const double a = std::exp(th[0]), b = std::exp(th[1]);
            const double lf0 = (a - 1.0) * std::log(tk) - tk / b - std::lgamma(a) - a * std::log(b);
            const double lS0 = std::log(std::max(pfqn::detail::gammainc(a, tk / b).q, 1e-300));
            H0 = -lS0;
            logh0 = lf0 - lS0;
        } else if (f == "lognormal") {
            const double mu = th[0], sg = std::exp(th[1]);
            const double z = (std::log(tk) - mu) / sg;
            const double lf0 = -std::log(tk) - std::log(sg) - 0.5 * log2pi - z * z / 2.0;
            const double lS0 = std::log(std::max(0.5 * std::erfc(z / std::sqrt(2.0)), 1e-300));
            H0 = -lS0;
            logh0 = lf0 - lS0;
        } else {
            infer_smp_family_npar(f);
        }
        const double e = eta.empty() ? 0.0 : (eta.size() == 1 ? eta[0] : eta[k]);
        const double H = H0 * std::exp(e);
        logS[k] = -H;
        logf[k] = logh0 + e - H;
    }
}
 
/** `infer_smp_family(F, 'natural', THETA)`: rate | shape, scale | mu, sigma. */
inline std::vector<double> infer_smp_family_natural(const std::string& family,
                                                    const std::vector<double>& th) {
    const std::string f = smp_family_name(family);
    if (f == "exp") return std::vector<double>{std::exp(th[0])};
    if (f == "weibull" || f == "gamma") return std::vector<double>{std::exp(th[0]), std::exp(th[1])};
    if (f == "lognormal") return std::vector<double>{th[0], std::exp(th[1])};
    infer_smp_family_npar(f);
    return std::vector<double>();
}
 
/**
 * `infer_smp_family(F, 'theta0', TEV, TEXP)`: the starting point from the
 * event times `tev` and the total exposure `texp`.
 */
inline std::vector<double> infer_smp_family_theta0(const std::string& family,
                                                   const std::vector<double>& tev, double texp) {
    const std::string f = smp_family_name(family);
    const double d = static_cast<double>(tev.size());
    if (f == "exp") return std::vector<double>{std::log(d / texp)};
    if (f == "weibull" || f == "gamma") return std::vector<double>{0.0, std::log(texp / d)};
    if (f == "lognormal") {
        std::vector<double> lt(tev.size());
        double mean = 0.0;
        for (std::size_t k = 0; k < tev.size(); ++k) {
            lt[k] = std::log(tev[k]);
            mean += lt[k];
        }
        mean /= d;
        double sg = 1.0;
        if (tev.size() >= 2) {
            // MATLAB's std: the sample standard deviation, n - 1 in the denominator.
            double ss = 0.0;
            for (double v : lt) ss += (v - mean) * (v - mean);
            const double sd = std::sqrt(ss / (d - 1.0));
            if (sd > 0.0) sg = sd;
        }
        return std::vector<double>{mean, std::log(sg)};
    }
    infer_smp_family_npar(f);
    return std::vector<double>();
}
 
/**
 * `infer_smp_family(F, 'dist', THETA, ETA)`: the LINE distribution of the law
 * with hazard h0(t) exp(eta). Only 'exp' and 'weibull' are closed under a
 * proportional-hazards factor; the other two are refused for eta != 0.
 */
inline lang::Distrib<double> infer_smp_family_dist(const std::string& family,
                                                   const std::vector<double>& th, double eta) {
    const std::string f = smp_family_name(family);
    const std::vector<double> p = infer_smp_family_natural(f, th);
    if (f == "exp") return lang::Distrib<double>::exp_rate(p[0] * std::exp(eta));
    // MATLAB's Weibull(shape, scale) is the C++ Distrib::weibull(scale, shape).
    if (f == "weibull") return lang::Distrib<double>::weibull(p[1] * std::exp(-eta / p[0]), p[0]);
    if (eta != 0.0)
        throw InputError("infer_smp_family: A '" + f +
                         "' law scaled by a proportional-hazards factor is not a '" + f +
                         "' law, so it has no LINE distribution: build the environment at the "
                         "baseline (z = 0), or fit 'exp' or 'weibull'.");
    if (f == "gamma") return lang::Distrib<double>::gamma_dist(p[0], p[1]);
    return lang::Distrib<double>::lognormal(p[0], p[1]);
}
 
// ---- infer_smp -------------------------------------------------------------
 
/** Options of `infer_smp`, the reference's name/value pairs. */
struct SmpOptions {
    /** 'intensity' (Approach II, default) or 'sojourn' (Approach I). */
    std::string approach = "intensity";
    /** The family of every transition: 'weibull' (default), 'exp', 'gamma' or 'lognormal'. */
    std::string dist = "weibull";
    /**
     * Per transition, `dist_table[i][j]` for i -> j (0-based states), when not
     * empty; it must then be E x E, and an empty entry takes 'weibull'.
     */
    std::vector<std::vector<std::string>> dist_table;
    /** n x p covariates Z, one row per sojourn; an empty matrix means none. */
    Matrix<double> covariates;
    /**
     * `covmask[i][j][c]`: whether covariate c enters transition i -> j (0-based
     * states); empty means all, otherwise E x E x p.
     */
    std::vector<std::vector<std::vector<bool>>> covmask;
    /** E; 0 means max(from, to). */
    std::size_t nstates = 0;
    double tolx = 1e-9;
    double tolf = 1e-12;
    /** Nelder-Mead iteration cap per run; 0 means 2000 per parameter. */
    std::size_t maxiter = 0;
};
 
/** One fitted transition. States are 1-based. */
struct SmpArc {
    std::size_t from = 0, to = 0;
    std::string family;
    std::vector<double> theta;   ///< unconstrained parameters
    std::vector<double> params;  ///< natural parameters
    std::vector<double> beta;    ///< p covariate effects, zero where masked out
    std::vector<double> se;      ///< standard errors of [theta; beta]
    /** Approach I: p_ij; Approach II: the embedded probability at z = 0. */
    double prob = std::numeric_limits<double>::quiet_NaN();
    std::size_t nevents = 0;
};
 
/** What `infer_smp` returns. */
struct SmpFit {
    std::string approach;
    std::size_t nstates = 0;
    std::size_t ncov = 0;
    std::vector<SmpArc> arcs;
    double loglik = 0.0;
    std::size_t nparams = 0;
    double aic = 0.0;
    Matrix<double> P;  ///< E x E embedded jump chain at z = 0
};
 
namespace smp_infer_detail {
 
using Objective = std::function<double(const std::vector<double>&)>;
 
/** MATLAB's `max` over a vector: NaNs ignored, NaN when every entry is one. */
inline double nanmax(const std::vector<double>& v) {
    double m = std::numeric_limits<double>::quiet_NaN();
    for (double x : v) {
        if (std::isnan(x)) continue;
        if (std::isnan(m) || x > m) m = x;
    }
    return m;
}
 
/** One run of the reference's Nelder-Mead (`nmRun`). */
inline void nm_run(const Objective& f, const std::vector<double>& x0, double tolx, double tolf,
                   std::size_t maxit, std::vector<double>& xb, double& fb) {
    const std::size_t n = x0.size();
    std::vector<std::vector<double>> v(n + 1, x0);
    std::vector<double> fv(n + 1, 0.0);
    fv[0] = f(x0);
    for (std::size_t j = 0; j < n; ++j) {
        std::vector<double> y = x0;
        if (y[j] != 0.0)
            y[j] = 1.05 * y[j];
        else
            y[j] = 0.00025;
        v[j + 1] = y;
        fv[j + 1] = f(y);
    }
    // MATLAB's sort is stable, so ties keep their order.
    const auto sort_simplex = [&]() {
        std::vector<std::size_t> ix(n + 1);
        for (std::size_t k = 0; k <= n; ++k) ix[k] = k;
        std::stable_sort(ix.begin(), ix.end(), [&](std::size_t a, std::size_t b) {
            const double fa = fv[a], fb2 = fv[b];
            if (std::isnan(fb2)) return !std::isnan(fa);
            if (std::isnan(fa)) return false;
            return fa < fb2;
        });
        std::vector<std::vector<double>> v2(n + 1);
        std::vector<double> f2(n + 1);
        for (std::size_t k = 0; k <= n; ++k) {
            v2[k] = v[ix[k]];
            f2[k] = fv[ix[k]];
        }
        v.swap(v2);
        fv.swap(f2);
    };
    sort_simplex();
    std::size_t it = 0;
    while (it < maxit) {
        std::vector<double> df(n), dx;
        dx.reserve(n * n);
        for (std::size_t j = 1; j <= n; ++j) {
            df[j - 1] = std::fabs(fv[0] - fv[j]);
            for (std::size_t a = 0; a < n; ++a) dx.push_back(std::fabs(v[j][a] - v[0][a]));
        }
        if (nanmax(df) <= tolf && nanmax(dx) <= tolx) break;
        std::vector<double> xbar(n, 0.0);
        for (std::size_t a = 0; a < n; ++a) {
            double s = 0.0;
            for (std::size_t j = 0; j < n; ++j) s += v[j][a];
            xbar[a] = s / static_cast<double>(n);
        }
        std::vector<double>& vw = v[n];
        std::vector<double> xr(n);
        for (std::size_t a = 0; a < n; ++a) xr[a] = 2.0 * xbar[a] - vw[a];
        const double fr = f(xr);
        if (fr < fv[0]) {
            std::vector<double> xe(n);
            for (std::size_t a = 0; a < n; ++a) xe[a] = 3.0 * xbar[a] - 2.0 * vw[a];
            const double fe = f(xe);
            if (fe < fr) {
                vw = xe;
                fv[n] = fe;
            } else {
                vw = xr;
                fv[n] = fr;
            }
        } else if (fr < fv[n - 1]) {
            vw = xr;
            fv[n] = fr;
        } else {
            bool shrink = false;
            if (fr < fv[n]) {
                std::vector<double> xc(n);
                for (std::size_t a = 0; a < n; ++a) xc[a] = 1.5 * xbar[a] - 0.5 * vw[a];
                const double fc = f(xc);
                if (fc <= fr) {
                    vw = xc;
                    fv[n] = fc;
                } else {
                    shrink = true;
                }
            } else {
                std::vector<double> xcc(n);
                for (std::size_t a = 0; a < n; ++a) xcc[a] = 0.5 * xbar[a] + 0.5 * vw[a];
                const double fcc = f(xcc);
                if (fcc < fv[n]) {
                    vw = xcc;
                    fv[n] = fcc;
                } else {
                    shrink = true;
                }
            }
            if (shrink)
                for (std::size_t j = 1; j <= n; ++j) {
                    for (std::size_t a = 0; a < n; ++a)
                        v[j][a] = v[0][a] + 0.5 * (v[j][a] - v[0][a]);
                    fv[j] = f(v[j]);
                }
        }
        sort_simplex();
        ++it;
    }
    xb = v[0];
    fb = fv[0];
}
 
/** `nelderMead`: restart `nm_run` from its own optimum, at most four runs. */
inline std::vector<double> nelder_mead(const Objective& f, const std::vector<double>& x0,
                                       const SmpOptions& opt) {
    const std::size_t maxit = opt.maxiter > 0 ? opt.maxiter : 2000 * x0.size();
    std::vector<double> x = x0;
    double fx = f(x);
    for (int restart = 0; restart < 4; ++restart) {
        std::vector<double> xn;
        double fn = 0.0;
        nm_run(f, x, opt.tolx, opt.tolf, maxit, xn, fn);
        const bool improved = fx - fn > opt.tolf;
        if (fn <= fx) {
            x = xn;
            fx = fn;
        }
        if (!improved) break;
    }
    return x;
}
 
/**
 * `hessianSE`: central-difference Hessian with steps 1e-4 max(1, |x|), its
 * Cholesky factor H = R'R, and se = sqrt(rowsum((R \ I).^2)); NaN everywhere
 * when H is not finite or not positive definite.
 */
inline std::vector<double> hessian_se(const Objective& f, const std::vector<double>& x) {
    const std::size_t n = x.size();
    std::vector<double> h(n);
    for (std::size_t a = 0; a < n; ++a) h[a] = 1e-4 * std::max(1.0, std::fabs(x[a]));
    Matrix<double> H(n, n, 0.0);
    const double f0 = f(x);
    for (std::size_t a = 0; a < n; ++a) {
        std::vector<double> xp = x, xm = x;
        xp[a] += h[a];
        xm[a] -= h[a];
        H(a, a) = (f(xp) - 2.0 * f0 + f(xm)) / (h[a] * h[a]);
        for (std::size_t b = a + 1; b < n; ++b) {
            std::vector<double> pp = xp, pm = xp, mp = xm, mm = xm;
            pp[b] += h[b];
            pm[b] -= h[b];
            mp[b] += h[b];
            mm[b] -= h[b];
            H(a, b) = (f(pp) - f(pm) - f(mp) + f(mm)) / (4.0 * h[a] * h[b]);
            H(b, a) = H(a, b);
        }
    }
    const std::vector<double> nan(n, std::numeric_limits<double>::quiet_NaN());
    for (std::size_t a = 0; a < n; ++a)
        for (std::size_t b = 0; b < n; ++b)
            if (!std::isfinite(H(a, b))) return nan;
    // chol(H): the upper factor R, reading the upper triangle as MATLAB does.
    Matrix<double> R(n, n, 0.0);
    for (std::size_t j = 0; j < n; ++j) {
        double s = H(j, j);
        for (std::size_t k = 0; k < j; ++k) s -= R(k, j) * R(k, j);
        if (!(s > 0.0)) return nan;
        R(j, j) = std::sqrt(s);
        for (std::size_t i = j + 1; i < n; ++i) {
            double t = H(j, i);
            for (std::size_t k = 0; k < j; ++k) t -= R(k, j) * R(k, i);
            R(j, i) = t / R(j, j);
        }
    }
    // Ri = R \ I by back substitution, column by column.
    Matrix<double> Ri(n, n, 0.0);
    for (std::size_t c = 0; c < n; ++c)
        for (std::size_t ii = n; ii-- > 0;) {
            double s = (ii == c) ? 1.0 : 0.0;
            for (std::size_t k = ii + 1; k < n; ++k) s -= R(ii, k) * Ri(k, c);
            Ri(ii, c) = s / R(ii, ii);
        }
    std::vector<double> se(n, 0.0);
    for (std::size_t a = 0; a < n; ++a) {
        double s = 0.0;
        for (std::size_t c = 0; c < n; ++c) s += Ri(a, c) * Ri(a, c);
        se[a] = std::sqrt(s);
    }
    return se;
}
 
/** `Z(rows, cov) * beta`, one value per row. */
inline std::vector<double> linpred(const Matrix<double>& Z, const std::vector<std::size_t>& rows,
                                   const std::vector<std::size_t>& cov, const double* beta) {
    std::vector<double> eta(rows.size(), 0.0);
    for (std::size_t r = 0; r < rows.size(); ++r) {
        double s = 0.0;
        for (std::size_t c = 0; c < cov.size(); ++c) s += Z(rows[r], cov[c]) * beta[c];
        eta[r] = s;
    }
    return eta;
}
 
/** `fitIntensityArc`: one transition of Approach II, its events `ev` among the sojourns `t`. */
inline void fit_intensity_arc(const std::string& family, const std::vector<double>& t,
                              const std::vector<char>& ev, const Matrix<double>& Z,
                              const std::vector<std::size_t>& rows,
                              const std::vector<std::size_t>& cov, const SmpOptions& opt,
                              std::vector<double>& th, std::vector<double>& be,
                              std::vector<double>& se, double& ll) {
    const std::size_t q = cov.size();
    const std::size_t npar = infer_smp_family_npar(family);
    std::vector<double> tev;
    double tsum = 0.0;
    for (std::size_t k = 0; k < t.size(); ++k) {
        tsum += t[k];
        if (ev[k]) tev.push_back(t[k]);
    }
    const Objective nll = [&](const std::vector<double>& x) {
        std::vector<double> eta;
        if (q > 0) eta = linpred(Z, rows, cov, x.data() + npar);
        std::vector<double> lf, lS;
        infer_smp_family_logfs(family, std::vector<double>(x.begin(), x.begin() + npar), t, eta,
                               lf, lS);
        double sev = 0.0, scen = 0.0;
        for (std::size_t k = 0; k < t.size(); ++k) {
            if (ev[k])
                sev += lf[k];
            else
                scen += lS[k];
        }
        double v = -(sev + scen);
        if (!std::isfinite(v)) v = std::numeric_limits<double>::infinity();
        return v;
    };
    std::vector<double> x;
    if (smp_family_name(family) == "exp" && q == 0) {
        x.push_back(std::log(static_cast<double>(tev.size()) / tsum));
    } else {
        std::vector<double> x0 = infer_smp_family_theta0(family, tev, tsum);
        x0.resize(npar + q, 0.0);
        x = nelder_mead(nll, x0, opt);
    }
    th.assign(x.begin(), x.begin() + npar);
    be.assign(x.begin() + npar, x.end());
    ll = -nll(x);
    se = hessian_se(nll, x);
}
 
/** `expandBeta`: the effects of the selected covariates, zero elsewhere. */
inline std::vector<double> expand_beta(const std::vector<double>& be,
                                       const std::vector<std::size_t>& cov, std::size_t p) {
    std::vector<double> b(p, 0.0);
    for (std::size_t c = 0; c < cov.size(); ++c) b[cov[c]] = be[c];
    return b;
}
 
/**
 * `fitSojournState`: state i of Approach I. The parameter vector is laid out
 * [logits of destinations 2..m; theta_1; beta_1; theta_2; beta_2; ...].
 */
inline void fit_sojourn_state(std::size_t i, const std::vector<std::size_t>& dests,
                              const std::vector<std::vector<std::string>>& fam,
                              const std::vector<double>& t, const std::vector<std::size_t>& to,
                              const Matrix<double>& Z, const std::vector<std::size_t>& rows,
                              const std::vector<std::vector<std::vector<bool>>>& mask,
                              std::size_t p, const SmpOptions& opt, std::vector<SmpArc>& arcs,
                              double& ll, std::size_t& np) {
    const std::size_t m = dests.size(), n = t.size();
    std::vector<std::vector<std::size_t>> cov(m);
    std::vector<std::size_t> npar(m, 0);
    std::vector<std::string> famk(m);
    for (std::size_t k = 0; k < m; ++k) {
        for (std::size_t c = 0; c < p; ++c)
            if (mask[i - 1][dests[k] - 1][c]) cov[k].push_back(c);
        famk[k] = fam[i - 1][dests[k] - 1];
        npar[k] = infer_smp_family_npar(famk[k]);
    }
    std::vector<char> cens(n, 0);
    bool anycens = false;
    for (std::size_t r = 0; r < n; ++r)
        if (to[r] == 0) {
            cens[r] = 1;
            anycens = true;
        }
    std::vector<std::size_t> off(m, 0);
    std::size_t pos = m - 1;
    for (std::size_t k = 0; k < m; ++k) {
        off[k] = pos;
        pos += npar[k] + cov[k].size();
    }
    // Per destination, the rows that ended there.
    std::vector<std::vector<std::size_t>> evidx(m);
    for (std::size_t k = 0; k < m; ++k)
        for (std::size_t r = 0; r < n; ++r)
            if (to[r] == dests[k]) evidx[k].push_back(r);
 
    const auto nll_dest = [&](std::size_t k, const std::vector<double>& y) {
        std::vector<double> tev(evidx[k].size()), eta;
        std::vector<std::size_t> erows(evidx[k].size());
        for (std::size_t q = 0; q < evidx[k].size(); ++q) {
            tev[q] = t[evidx[k][q]];
            erows[q] = rows[evidx[k][q]];
        }
        if (!cov[k].empty()) eta = linpred(Z, erows, cov[k], y.data() + npar[k]);
        std::vector<double> lf, lS;
        infer_smp_family_logfs(famk[k], std::vector<double>(y.begin(), y.begin() + npar[k]), tev,
                               eta, lf, lS);
        double s = 0.0;
        for (double v : lf) s += v;
        double v = -s;
        if (!std::isfinite(v)) v = std::numeric_limits<double>::infinity();
        return v;
    };
    const auto softmax = [&](const std::vector<double>& x) {
        std::vector<double> l(m, 0.0);
        for (std::size_t k = 1; k < m; ++k) l[k] = x[k - 1];
        const double mx = *std::max_element(l.begin(), l.end());
        double tot = 0.0;
        for (std::size_t k = 0; k < m; ++k) {
            l[k] -= mx;
            tot += std::exp(l[k]);
        }
        std::vector<double> pk(m);
        for (std::size_t k = 0; k < m; ++k) pk[k] = std::exp(l[k]) / tot;
        return pk;
    };
    std::vector<std::size_t> crow;
    for (std::size_t r = 0; r < n; ++r)
        if (cens[r]) crow.push_back(r);
    const Objective nll_state = [&](const std::vector<double>& x) {
        const std::vector<double> pk = softmax(x);
        double v = 0.0;
        std::vector<double> Sc(crow.size(), 0.0);
        for (std::size_t kk = 0; kk < m; ++kk) {
            const std::vector<double> th(x.begin() + off[kk], x.begin() + off[kk] + npar[kk]);
            std::vector<double> eta(n, 0.0);
            if (!cov[kk].empty()) eta = linpred(Z, rows, cov[kk], x.data() + off[kk] + npar[kk]);
            std::vector<double> lf, lS;
            infer_smp_family_logfs(famk[kk], th, t, eta, lf, lS);
            double s = 0.0;
            const double lp = std::log(pk[kk]);
            for (std::size_t r : evidx[kk]) s += lp + lf[r];
            v = v - s;
            for (std::size_t q = 0; q < crow.size(); ++q) Sc[q] = Sc[q] + pk[kk] * std::exp(lS[crow[q]]);
        }
        double s = 0.0;
        for (double c : Sc) s += std::log(std::max(c, 1e-300));
        v = v - s;
        if (!std::isfinite(v)) v = std::numeric_limits<double>::infinity();
        return v;
    };
 
    std::vector<double> x(pos, 0.0);
    std::vector<std::size_t> cnt(m, 0);
    for (std::size_t k = 0; k < m; ++k) {
        cnt[k] = evidx[k].size();
        std::vector<double> tev(cnt[k]);
        double texp = 0.0;
        for (std::size_t q = 0; q < cnt[k]; ++q) {
            tev[q] = t[evidx[k][q]];
            texp += tev[q];
        }
        std::vector<double> xk = infer_smp_family_theta0(famk[k], tev, texp);
        xk.resize(npar[k] + cov[k].size(), 0.0);
        if (!(smp_family_name(famk[k]) == "exp" && cov[k].empty()))
            xk = nelder_mead([&](const std::vector<double>& y) { return nll_dest(k, y); }, xk, opt);
        for (std::size_t q = 0; q < xk.size(); ++q) x[off[k] + q] = xk[q];
    }
    for (std::size_t k = 1; k < m; ++k)
        x[k - 1] = std::log(static_cast<double>(cnt[k]) / static_cast<double>(cnt[0]));
    if (anycens) x = nelder_mead(nll_state, x, opt);
    ll = -nll_state(x);
    np = pos;
    const std::vector<double> Sfull = hessian_se(nll_state, x);
    const std::vector<double> pr = softmax(x);
    for (std::size_t k = 0; k < m; ++k) {
        const std::size_t b0 = off[k], b1 = off[k] + npar[k] + cov[k].size();
        SmpArc a;
        a.from = i;
        a.to = dests[k];
        a.family = famk[k];
        a.theta.assign(x.begin() + off[k], x.begin() + off[k] + npar[k]);
        a.params = infer_smp_family_natural(famk[k], a.theta);
        a.beta = expand_beta(std::vector<double>(x.begin() + off[k] + npar[k], x.begin() + b1),
                             cov[k], p);
        if (anycens) {
            a.se.assign(Sfull.begin() + b0, Sfull.begin() + b1);
        } else {
            a.se = hessian_se([&](const std::vector<double>& y) { return nll_dest(k, y); },
                              std::vector<double>(x.begin() + b0, x.begin() + b1));
        }
        a.prob = pr[k];
        a.nevents = cnt[k];
        arcs.push_back(a);
    }
}
 
}  // namespace smp_infer_detail
 
/**
 * `infer_smp(from, to, dur, ...)`: the fit, under `opt.approach`.
 *
 * `from` holds states 1..E and `to` states 1..E or 0 for a censored sojourn;
 * `dur` the sojourn durations. See the file comment for the two approaches.
 */
inline SmpFit infer_smp(const std::vector<std::size_t>& from, const std::vector<std::size_t>& to,
                        const std::vector<double>& dur, const SmpOptions& opt = SmpOptions()) {
    using namespace smp_infer_detail;
    const std::size_t n = from.size();
    if (to.size() != n || dur.size() != n || n == 0)
        throw InputError("infer_smp: FROM, TO and DUR must be nonempty and of equal length.");
    std::size_t smax = 0;
    for (std::size_t k = 0; k < n; ++k) {
        if (from[k] < 1)
            throw InputError("infer_smp: FROM must hold states 1..E and TO states 1..E, or 0 for "
                             "a censored sojourn.");
        if (to[k] == from[k])
            throw InputError("infer_smp: A sojourn cannot end in the state it started from.");
        if (!std::isfinite(dur[k]) || !(dur[k] > 0.0))
            throw InputError("infer_smp: Every sojourn duration must be finite and positive.");
        smax = std::max(smax, std::max(from[k], to[k]));
    }
    const std::size_t E = opt.nstates > 0 ? opt.nstates : smax;
    if (smax > E) throw InputError("infer_smp: A state exceeds 'nstates'.");
    Matrix<double> Z = opt.covariates;
    if (Z.rows() == 0 && Z.cols() == 0) Z = Matrix<double>(n, 0, 0.0);
    if (Z.rows() != n)
        throw InputError("infer_smp: The covariates must be a finite n x p matrix, one row per "
                         "sojourn.");
    for (std::size_t r = 0; r < Z.rows(); ++r)
        for (std::size_t c = 0; c < Z.cols(); ++c)
            if (!std::isfinite(Z(r, c)))
                throw InputError("infer_smp: The covariates must be a finite n x p matrix, one row "
                                 "per sojourn.");
    const std::size_t p = Z.cols();
    std::vector<std::vector<std::vector<bool>>> mask = opt.covmask;
    if (mask.empty()) {
        mask.assign(E, std::vector<std::vector<bool>>(E, std::vector<bool>(p, true)));
    } else if (p > 0) {
        bool ok = mask.size() == E;
        for (std::size_t i = 0; ok && i < E; ++i) {
            ok = mask[i].size() == E;
            for (std::size_t j = 0; ok && j < E; ++j) ok = mask[i][j].size() == p;
        }
        if (!ok) throw InputError("infer_smp: The covariate mask must be E x E x p.");
    }
    // familyTable
    std::vector<std::vector<std::string>> fam(E, std::vector<std::string>(E));
    if (opt.dist_table.empty()) {
        for (std::size_t i = 0; i < E; ++i)
            for (std::size_t j = 0; j < E; ++j) fam[i][j] = smp_family_name(opt.dist);
    } else {
        bool ok = opt.dist_table.size() == E;
        for (std::size_t i = 0; ok && i < E; ++i) ok = opt.dist_table[i].size() == E;
        if (!ok)
            throw InputError("infer_smp: The 'dist' option must be a family name or an E x E cell "
                             "of names.");
        for (std::size_t i = 0; i < E; ++i)
            for (std::size_t j = 0; j < E; ++j)
                fam[i][j] = opt.dist_table[i][j].empty() ? std::string("weibull")
                                                         : smp_family_name(opt.dist_table[i][j]);
    }
    for (std::size_t i = 0; i < E; ++i)
        for (std::size_t j = 0; j < E; ++j) infer_smp_family_npar(fam[i][j]);
    const std::string approach = smp_family_name(opt.approach);
    if (approach != "intensity" && approach != "sojourn")
        throw InputError("infer_smp: The approach must be 'intensity' or 'sojourn'.");
 
    SmpFit fit;
    fit.approach = approach;
    fit.nstates = E;
    fit.ncov = p;
    double loglik = 0.0;
    std::size_t nparams = 0;
    for (std::size_t i = 1; i <= E; ++i) {
        std::vector<std::size_t> rows;
        for (std::size_t k = 0; k < n; ++k)
            if (from[k] == i) rows.push_back(k);
        if (rows.empty()) continue;
        std::vector<std::size_t> dests;
        for (std::size_t r : rows)
            if (to[r] > 0) dests.push_back(to[r]);
        std::sort(dests.begin(), dests.end());
        dests.erase(std::unique(dests.begin(), dests.end()), dests.end());
        if (dests.empty()) continue;
        std::vector<double> t(rows.size());
        std::vector<std::size_t> tor(rows.size());
        for (std::size_t q = 0; q < rows.size(); ++q) {
            t[q] = dur[rows[q]];
            tor[q] = to[rows[q]];
        }
        if (approach == "intensity") {
            for (std::size_t j : dests) {
                std::vector<std::size_t> cov;
                for (std::size_t c = 0; c < p; ++c)
                    if (mask[i - 1][j - 1][c]) cov.push_back(c);
                std::vector<char> ev(rows.size(), 0);
                std::size_t nev = 0;
                for (std::size_t q = 0; q < rows.size(); ++q)
                    if (tor[q] == j) {
                        ev[q] = 1;
                        ++nev;
                    }
                SmpArc a;
                std::vector<double> be;
                double ll = 0.0;
                fit_intensity_arc(fam[i - 1][j - 1], t, ev, Z, rows, cov, opt, a.theta, be, a.se,
                                  ll);
                a.from = i;
                a.to = j;
                a.family = fam[i - 1][j - 1];
                a.params = infer_smp_family_natural(a.family, a.theta);
                a.beta = expand_beta(be, cov, p);
                a.nevents = nev;
                fit.arcs.push_back(a);
                loglik += ll;
                nparams += a.theta.size() + cov.size();
            }
        } else {
            double ll = 0.0;
            std::size_t np = 0;
            fit_sojourn_state(i, dests, fam, t, tor, Z, rows, mask, p, opt, fit.arcs, ll, np);
            loglik += ll;
            nparams += np;
        }
    }
    fit.loglik = loglik;
    fit.nparams = nparams;
    fit.aic = 2.0 * static_cast<double>(nparams) - 2.0 * loglik;
    fit.P = Matrix<double>(E, E, 0.0);
    if (approach == "sojourn") {
        for (const SmpArc& a : fit.arcs) fit.P(a.from - 1, a.to - 1) = a.prob;
    } else {
        for (std::size_t i = 1; i <= E; ++i) {
            std::vector<lang::Distrib<double>> stage(E);
            bool any = false;
            for (const SmpArc& a : fit.arcs)
                if (a.from == i) {
                    stage[a.to - 1] = infer_smp_family_dist(a.family, a.theta, 0.0);
                    any = true;
                }
            if (!any) continue;
            const mc::SmpStageQuad q = mc::smp_stage_quad(stage, std::vector<double>());
            for (std::size_t j = 0; j < E; ++j) fit.P(i - 1, j) = q.P[j];
        }
        for (SmpArc& a : fit.arcs) a.prob = fit.P(a.from - 1, a.to - 1);
    }
    return fit;
}
 
// ---- infer_smp_env ---------------------------------------------------------
 
/** Options of `infer_smp_env`. */
struct SmpEnvOptions {
    std::string name = "SMPFit";
    std::vector<std::string> names;  ///< default S1..SE
    std::vector<std::string> types;  ///< default 'stage'
    std::vector<double> z;           ///< covariate vector, default the baseline 0
};
 
/**
 * `infer_smp_env(fit, models, ...)`: the Environment of a fit, `models[e]` the
 * network of stage e. An 'intensity' fit becomes competing arcs
 * (`add_transition`), each the fitted latent clock; a 'sojourn' fit becomes
 * kernel arcs (`add_kernel_transition`) with the fitted embedded probability and
 * conditional sojourn law. Either is solved by SolverENV method 'smp'.
 */
inline env::Environment<double> infer_smp_env(const SmpFit& fit,
                                              const std::vector<qn::NetworkStruct<double>>& models,
                                              const SmpEnvOptions& opt = SmpEnvOptions()) {
    const std::size_t E = fit.nstates;
    if (models.size() != E)
        throw InputError("infer_smp_env: MODELS must be a cell of " + std::to_string(E) +
                         " networks, one per state.");
    std::vector<std::string> names = opt.names, types = opt.types;
    if (names.empty())
        for (std::size_t e = 1; e <= E; ++e) names.push_back("S" + std::to_string(e));
    if (types.empty()) types.assign(E, "stage");
    if (names.size() != E || types.size() != E)
        throw InputError("infer_smp_env: 'names' and 'types' must have " + std::to_string(E) +
                         " entries.");
    std::vector<double> z = opt.z;
    if (z.empty()) z.assign(fit.ncov, 0.0);
    if (z.size() != fit.ncov)
        throw InputError("infer_smp_env: 'z' must have " + std::to_string(fit.ncov) +
                         " entries, one per covariate.");
    env::Environment<double> envm(opt.name, E);
    for (std::size_t e = 0; e < E; ++e) envm.set_stage(e, names[e], types[e], models[e]);
    for (const SmpArc& a : fit.arcs) {
        double eta = 0.0;
        for (std::size_t c = 0; c < z.size(); ++c) eta += a.beta[c] * z[c];
        const lang::Distrib<double> d = infer_smp_family_dist(a.family, a.theta, eta);
        if (fit.approach == "sojourn")
            envm.add_kernel_transition(a.from - 1, a.to - 1, a.prob, d);
        else
            envm.add_transition(a.from - 1, a.to - 1, d);
    }
    return envm;
}
 
// ---- infer_smp_hazard ------------------------------------------------------
 
/** Both hazard functions of one transition, rows aligned with the instants asked for. */
struct SmpHazard {
    std::size_t from = 0, to = 0;
    std::vector<double> intensity;   ///< a~_ij(t), the transition intensity
    std::vector<double> sojourn;     ///< alpha_ij(t), the hazard of the conditional sojourn law
    std::vector<double> cumulative;  ///< A~_ij(t) = int_0^t a~_ij(s) ds
    std::vector<double> cif;         ///< CIF_ij(t) = P(sojourn <= t, next state j)
};
 
/**
 * `infer_smp_hazard(fit, t, 'z', z)`: both hazard functions of every arc of a
 * fit on the ascending instants t > 0 (Asanjarani et al., Eqs. 11-15), related
 * by a~_ij(t) = p_ij f_ij(t) / S_i(t). An 'intensity' fit has a~ and A~ in
 * closed form and the CIF by `mc::smp_stage_cif`, so alpha_ij =
 * a~_ij S_i / (p_ij - CIF_ij). A 'sojourn' fit has alpha and the CIF in closed
 * form, a~ from the relation, and A~ as a midpoint sum on the grid
 * [0, logspace(log10(tmax)-8, log10(tmax), 2001)] joined with
 * linspace(0, tmax, 20001) and t. The result is aligned with `fit.arcs`.
 */
inline std::vector<SmpHazard> infer_smp_hazard(const SmpFit& fit, const std::vector<double>& t,
                                               const std::vector<double>& zin = std::vector<double>()) {
    std::vector<double> z = zin;
    if (z.empty()) z.assign(fit.ncov, 0.0);
    if (z.size() != fit.ncov)
        throw InputError("infer_smp_hazard: 'z' must have " + std::to_string(fit.ncov) +
                         " entries, one per covariate.");
    bool ok = !t.empty();
    for (std::size_t k = 0; ok && k < t.size(); ++k) {
        if (!std::isfinite(t[k]) || !(t[k] > 0.0)) ok = false;
        if (k > 0 && !(t[k] - t[k - 1] > 0.0)) ok = false;
    }
    if (!ok) throw InputError("infer_smp_hazard: T must be ascending, finite and positive.");
    const std::size_t E = fit.nstates, nt = t.size();
    std::vector<SmpHazard> hz(fit.arcs.size());
    std::vector<std::size_t> froms;
    for (const SmpArc& a : fit.arcs) froms.push_back(a.from);
    std::sort(froms.begin(), froms.end());
    froms.erase(std::unique(froms.begin(), froms.end()), froms.end());
    const auto eta_of = [&](const SmpArc& a) {
        double s = 0.0;
        for (std::size_t c = 0; c < z.size(); ++c) s += a.beta[c] * z[c];
        return s;
    };
    for (std::size_t i : froms) {
        std::vector<std::size_t> idx;
        for (std::size_t a = 0; a < fit.arcs.size(); ++a)
            if (fit.arcs[a].from == i) idx.push_back(a);
        const std::size_t m = idx.size();
        std::vector<std::vector<double>> lf(m), lS(m);
        for (std::size_t k = 0; k < m; ++k) {
            const SmpArc& a = fit.arcs[idx[k]];
            infer_smp_family_logfs(a.family, a.theta, t, std::vector<double>{eta_of(a)}, lf[k],
                                   lS[k]);
        }
        if (fit.approach == "intensity") {
            std::vector<lang::Distrib<double>> stage(E);
            for (std::size_t k = 0; k < m; ++k) {
                const SmpArc& a = fit.arcs[idx[k]];
                stage[a.to - 1] = infer_smp_family_dist(a.family, a.theta, eta_of(a));
            }
            const mc::SmpStageQuad q = mc::smp_stage_quad(stage, std::vector<double>());
            const mc::SmpStageCif C = mc::smp_stage_cif(stage, std::vector<double>(), t);
            std::vector<double> Si(nt, 0.0);
            for (std::size_t j = 0; j < nt; ++j) {
                double s = 0.0;
                for (std::size_t k = 0; k < m; ++k) s += lS[k][j];
                Si[j] = std::exp(s);
            }
            for (std::size_t k = 0; k < m; ++k) {
                const SmpArc& a = fit.arcs[idx[k]];
                SmpHazard& h = hz[idx[k]];
                h.from = i;
                h.to = a.to;
                h.intensity.resize(nt);
                h.sojourn.resize(nt);
                h.cumulative.resize(nt);
                h.cif.resize(nt);
                for (std::size_t j = 0; j < nt; ++j) {
                    const double hh = std::exp(lf[k][j] - lS[k][j]);
                    const double cif = C.C(a.to - 1, j);
                    h.intensity[j] = hh;
                    h.sojourn[j] = hh * Si[j] / std::max(q.P[a.to - 1] - cif, 0.0);
                    h.cumulative[j] = -lS[k][j];
                    h.cif[j] = cif;
                }
            }
        } else {
            std::vector<double> pk(m);
            for (std::size_t k = 0; k < m; ++k) pk[k] = fit.arcs[idx[k]].prob;
            std::vector<double> Si(nt, 0.0);
            for (std::size_t j = 0; j < nt; ++j)
                for (std::size_t k = 0; k < m; ++k) Si[j] += std::exp(lS[k][j]) * pk[k];
            const double tmax = t.back();
            std::vector<double> g(1, 0.0);
            const std::vector<double> ls =
                mc::smp_detail::linspace(std::log10(tmax) - 8.0, std::log10(tmax), 2001);
            for (double v : ls) g.push_back(std::pow(10.0, v));
            const std::vector<double> lin = mc::smp_detail::linspace(0.0, tmax, 20001);
            g.insert(g.end(), lin.begin(), lin.end());
            g.insert(g.end(), t.begin(), t.end());
            std::sort(g.begin(), g.end());
            g.erase(std::unique(g.begin(), g.end()), g.end());
            const std::size_t ng = g.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>> lfm(m), lSm(m);
            for (std::size_t k = 0; k < m; ++k) {
                const SmpArc& a = fit.arcs[idx[k]];
                infer_smp_family_logfs(a.family, a.theta, gm, std::vector<double>{eta_of(a)},
                                       lfm[k], lSm[k]);
            }
            std::vector<double> Sim(ng - 1, 0.0);
            for (std::size_t j = 0; j + 1 < ng; ++j)
                for (std::size_t k = 0; k < m; ++k) Sim[j] += std::exp(lSm[k][j]) * pk[k];
            std::vector<std::size_t> at(nt);
            for (std::size_t j = 0; j < nt; ++j)
                at[j] = static_cast<std::size_t>(std::lower_bound(g.begin(), g.end(), t[j]) -
                                                 g.begin());
            for (std::size_t k = 0; k < m; ++k) {
                const SmpArc& a = fit.arcs[idx[k]];
                std::vector<double> cum(ng, 0.0);
                for (std::size_t j = 0; j + 1 < ng; ++j)
                    cum[j + 1] = cum[j] + pk[k] * std::exp(lfm[k][j]) / Sim[j] * (g[j + 1] - g[j]);
                SmpHazard& h = hz[idx[k]];
                h.from = i;
                h.to = a.to;
                h.intensity.resize(nt);
                h.sojourn.resize(nt);
                h.cumulative.resize(nt);
                h.cif.resize(nt);
                for (std::size_t j = 0; j < nt; ++j) {
                    h.intensity[j] = pk[k] * std::exp(lf[k][j]) / Si[j];
                    h.sojourn[j] = std::exp(lf[k][j] - lS[k][j]);
                    h.cumulative[j] = cum[at[j]];
                    h.cif[j] = pk[k] * (1.0 - std::exp(lS[k][j]));
                }
            }
        }
    }
    return hz;
}
 
// ---- infer_smp_nelson_aalen ------------------------------------------------
 
/** One Nelson-Aalen estimate, on the distinct event times of its transition. */
struct SmpNelsonAalen {
    std::size_t from = 0, to = 0;
    std::vector<double> time;  ///< the distinct event times of i -> j, ascending
    std::vector<double> A;     ///< the estimate there (right-continuous step function)
    std::vector<double> var;   ///< its Greenwood-type variance sum d/Y^2
};
 
/**
 * `infer_smp_nelson_aalen(from, to, dur, 'nstates', E)`: the Nelson-Aalen
 * estimate of the cumulative transition intensity of every observed
 * transition, on the clock reset at each jump. The risk set at t of i -> j is
 * the number of sojourns in i of duration at least t, so this is the
 * nonparametric benchmark of the 'intensity' approach (Asanjarani et al.,
 * Sec. 4). `nstates` 0 means max(from, to).
 */
inline std::vector<SmpNelsonAalen> infer_smp_nelson_aalen(const std::vector<std::size_t>& from,
                                                          const std::vector<std::size_t>& to,
                                                          const std::vector<double>& dur,
                                                          std::size_t nstates = 0) {
    const std::size_t n = from.size();
    if (to.size() != n || dur.size() != n || n == 0)
        throw InputError("infer_smp_nelson_aalen: FROM, TO and DUR must be nonempty and of equal "
                         "length.");
    std::size_t smax = 0;
    for (std::size_t k = 0; k < n; ++k) {
        if (!std::isfinite(dur[k]) || !(dur[k] > 0.0) || from[k] < 1 || to[k] == from[k])
            throw InputError("infer_smp_nelson_aalen: Invalid sojourn data: see INFER_SMP for the "
                             "coding.");
        smax = std::max(smax, std::max(from[k], to[k]));
    }
    const std::size_t E = nstates > 0 ? nstates : smax;
    std::vector<SmpNelsonAalen> na;
    for (std::size_t i = 1; i <= E; ++i) {
        std::vector<double> ti;
        std::vector<std::size_t> di;
        for (std::size_t k = 0; k < n; ++k)
            if (from[k] == i) {
                ti.push_back(dur[k]);
                di.push_back(to[k]);
            }
        std::vector<std::size_t> dests;
        for (std::size_t d : di)
            if (d > 0) dests.push_back(d);
        std::sort(dests.begin(), dests.end());
        dests.erase(std::unique(dests.begin(), dests.end()), dests.end());
        for (std::size_t j : dests) {
            std::vector<double> tev;
            for (std::size_t q = 0; q < ti.size(); ++q)
                if (di[q] == j) tev.push_back(ti[q]);
            std::sort(tev.begin(), tev.end());
            tev.erase(std::unique(tev.begin(), tev.end()), tev.end());
            SmpNelsonAalen r;
            r.from = i;
            r.to = j;
            r.time = tev;
            r.A.resize(tev.size());
            r.var.resize(tev.size());
            double a = 0.0, v = 0.0;
            for (std::size_t k = 0; k < tev.size(); ++k) {
                double d = 0.0, Y = 0.0;
                for (std::size_t q = 0; q < ti.size(); ++q) {
                    if (ti[q] == tev[k] && di[q] == j) d += 1.0;
                    if (ti[q] >= tev[k]) Y += 1.0;
                }
                a = a + d / Y;
                v = v + d / (Y * Y);
                r.A[k] = a;
                r.var[k] = v;
            }
            na.push_back(r);
        }
    }
    return na;
}
 
}  // namespace infer
}  // namespace line
 
#endif  // LINE_API_INFER_INFER_SMP_H