pfqn_mci.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_PFQN_MCI_H
#define LINE_API_PFQN_MCI_H
 
/**
 * @file
 * @ingroup api_pfqn
 * Monte Carlo Integration estimate of the normalizing constant of a closed
 * product-form network (Ross, Wang and Yao; MonteQueue 2.0).
 *
 * Templated port of matlab/src/api/pfqn/pfqn_mci.m, cross-checked against
 * jar/src/main/java/jline/api/pfqn/nc/Pfqn_mci.java and
 * python/line_solver/api/pfqn/asymptotic.py. Variants: mci, imci, amci,
 * lhsmci, rm.
 *
 * The McKenna-Mitra integral form of the constant,
 *
 *   G(N) = 1/prod_r N_r! int_{R_+^M} e^{-sum_i v_i} prod_r (sum_i v_i D(i,r) + Z_r)^{N_r} dv,
 *
 * is estimated by importance sampling with an independent exponential proposal
 * of rate gamma_i per station. Writing V ~ prod_i Exp(gamma_i), one draw
 * contributes
 *
 *   lZ = -sum_i (1 - gamma_i) V_i - sum_i log gamma_i - sum_r log N_r!
 *        + sum_r N_r log(sum_i V_i D(i,r) + Z_r),
 *
 * and lG = log mean exp(lZ). The proposal rates come from a Bard-Schweitzer
 * pre-solve: gamma_i = max(0.01, 1 - U_i) in the IMCI variant, and the
 * saturation-aware gamma_i = 1/sqrt(max N) when U_i > 0.9 in the plain MCI
 * variant. Overshooting rates make the estimator heavy-tailed, which is why
 * the reference clamps them.
 *
 * The deviates come from the caller's generator, one draw per station per
 * sample. (MATLAB used to cache its uniforms in a persistent across calls; that
 * cache was removed on 2026-07-28, so MATLAB now redraws per call as this port,
 * the JAR and python do.)
 *
 * Variant notes:
 *
 *  - The 'rm' variant reads the balanced-job-bound rate with the COLUMN maximum
 *    of the demands, tput_r = N_r / (sum_i D(i,r) + Z_r + max_i D(i,r) (sum N - 1)),
 *    for every M, as MATLAB (max(D,[],1)), the JAR and python do. Before
 *    2026-10-21 MATLAB wrote max(D,1), the elementwise maximum with the scalar 1,
 *    which for M > 1 built the rates from the FIRST station's clamped demands
 *    only, and the JAR copied those numbers. Only the proposal, hence the
 *    variance, was affected: every reading targets the same constant.
 *
 *  - 'amci' and 'lhsmci' use the 'imci' tilt and differ only in the sampler.
 *    'amci' draws antithetic pairs (u, 1-u) in rows s and ceil(I/2)+s, and
 *    'lhsmci' stratifies each station's uniforms by a Latin hypercube. Both fill
 *    the uniforms STATION-MAJOR (all draws for station 1, then station 2, ...),
 *    in the order of MATLAB's column-major rand(Ih,M) and of the JAR; the
 *    hypercube draws its Fisher-Yates permutation (as the JAR does, by
 *    mc_uniform_int) before the I jitters of that station. The plain samplers
 *    ('mci', 'imci', 'rm') keep the sample-major order this port has always
 *    used, so their seeded outputs are unchanged.
 *
 * Arithmetic: INEXACT BY CONSTRUCTION. The value is a random variable, the
 * integrand is evaluated through logs, and the proposal rates come from a
 * tolerance-stopped Bard-Schweitzer solve; all three require transcendental
 * arithmetic.
 *
 * RNG contract: see pfqn_mc_common.h. Comparable to MATLAB, the JAR and python
 * only in distribution, never stream for stream (MATLAB's rand/randperm,
 * java.util.Random and numpy's legacy generator are three different streams);
 * reproducible within this port only when the generator is passed in the same
 * state. The cross-codebase check is therefore statistical: the mean of
 * independent replications must sit within its own standard error of the
 * exact constant (line-test.git cpp/tests/test_pfqn_mci_variants.cpp).
 */
 
#include <cmath>
#include <cstddef>
#include <limits>
#include <vector>
 
#include "line/api/pfqn/pfqn_bs.h"
#include "line/api/pfqn/pfqn_ca.h"
#include "line/api/pfqn/pfqn_mc_common.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace pfqn {
 
/** The proposal-rate rules the reference selects between. */
enum class MciVariant {
    Imci,  ///< gamma = max(0.01, 1 - U), the MonteQueue 2.0 recommendation
    Mci,   ///< gamma = 1/sqrt(max N) where U > 0.9, else 1 - U
    Rm,    ///< repairman: rates from the balanced job bound, column-max reading
    Amci,  ///< the Imci tilt, antithetic uniform pairs (u, 1-u) (Ross and Wang)
    Lhsmci ///< the Imci tilt, Latin hypercube stratified uniforms per station
};
 
/**
 * @brief Monte Carlo Integration estimate of the normalizing constant of a
 *        closed product-form network (Ross, Wang and Yao; MonteQueue 2.0).
 *
 * @param D       (M x R) service demands
 * @param N       (R) population per class
 * @param Z       (R) think times; empty for none
 * @param samples number of Monte Carlo samples
 * @param variant proposal-rate rule
 * @param rng     explicit generator, advanced by the call
 */
template <class T>
NcResult<T> pfqn_mci(const Matrix<T>& D, const std::vector<int>& N, const std::vector<T>& Z,
                     std::size_t samples, MciVariant variant, McRng& rng) {
    static_assert(num_traits<T>::has_transcendental,
                  "pfqn_mci requires transcendental arithmetic: it is a Monte Carlo estimator over "
                  "an exponential proposal, its integrand is evaluated in the log domain, and its "
                  "proposal rates come from a tolerance-stopped Bard-Schweitzer solve");
 
    const std::size_t M = D.empty() ? 0 : D.rows();
    const std::size_t R = N.size();
    if (!D.empty() && D.cols() != R)
        throw InputError("pfqn_mci: D and N disagree on the class count");
    if (!Z.empty() && Z.size() != R) throw InputError("pfqn_mci: Z has the wrong length");
    for (int n : N)
        if (n < 0) throw InputError("pfqn_mci: negative population");
    if (samples == 0) throw InputError("pfqn_mci: at least one sample is required");
 
    const T zero = num_traits<T>::from_int(0);
    std::vector<T> Zv(R, zero);
    for (std::size_t r = 0; r < Z.size(); ++r) Zv[r] = Z[r];
 
    // ---- degenerate model: pure delay ------------------------------------
    T dsum = zero;
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t r = 0; r < R; ++r) dsum += D(i, r);
    if (M == 0 || dsum < num_traits<T>::from_double(1e-4)) {
        double lG = 0.0;
        for (std::size_t r = 0; r < R; ++r) {
            lG -= mc_log_factorial<T>(N[r]);
            if (N[r] > 0) lG += N[r] * num_traits<T>::log_as_double(Zv[r]);
        }
        return {mc_exp<T>(lG), lG};
    }
 
    // ---- proposal rates ---------------------------------------------------
    std::vector<double> gamma(M);
    if (variant == MciVariant::Rm) {
        long Ntot = 0;
        int Nmax = 0;
        for (int n : N) {
            Ntot += n;
            if (n > Nmax) Nmax = n;
        }
        std::vector<T> tput(R, zero);
        for (std::size_t r = 0; r < R; ++r) {
            T dsumr = zero;
            T dmax = D(0, r);
            for (std::size_t i = 0; i < M; ++i) {
                dsumr += D(i, r);
                if (D(i, r) > dmax) dmax = D(i, r);
            }
            const T den = dsumr + Zv[r] + dmax * num_traits<T>::from_int(Ntot - 1);
            if (den == zero) throw NumericError(ErrorKind{"Degenerate"}, "pfqn_mci: zero denominator in the 'rm' bound");
            tput[r] = num_traits<T>::from_int(N[r]) / den;
        }
        for (std::size_t i = 0; i < M; ++i) {
            T util = zero;
            for (std::size_t r = 0; r < R; ++r) util += D(i, r) * tput[r];
            const double u = num_traits<T>::to_double(util);
            gamma[i] = u > 0.9 ? 1.0 / std::sqrt(static_cast<double>(Nmax)) : 1.0 - u;
        }
    } else {
        std::vector<T> Nt(R), Zt(R);
        for (std::size_t r = 0; r < R; ++r) {
            Nt[r] = num_traits<T>::from_int(N[r]);
            Zt[r] = Zv[r];
        }
        const AmvaResult<T> bs = pfqn_bs(D, Nt, Zt);
        int Nmax = 0;
        for (int n : N)
            if (n > Nmax) Nmax = n;
        for (std::size_t i = 0; i < M; ++i) {
            T util = zero;
            for (std::size_t r = 0; r < R; ++r) util += D(i, r) * bs.XN[r];
            const double u = num_traits<T>::to_double(util);
            if (variant != MciVariant::Mci) {  // Imci, Amci and Lhsmci share the tilt
                gamma[i] = 1.0 - u > 0.01 ? 1.0 - u : 0.01;
            } else {
                gamma[i] = u > 0.9 ? 1.0 / std::sqrt(static_cast<double>(Nmax)) : 1.0 - u;
            }
        }
    }
    for (std::size_t i = 0; i < M; ++i)
        if (!(gamma[i] > 0.0))
            throw NumericError(ErrorKind{"Instability"},
                "pfqn_mci: a proposal rate is non-positive, the station is saturated beyond what "
                "the variant's clamp covers");
 
    // ---- constant part of every log-weight ---------------------------------
    double lconst = 0.0;
    for (std::size_t i = 0; i < M; ++i) lconst -= std::log(gamma[i]);
    for (std::size_t r = 0; r < R; ++r) lconst -= mc_log_factorial<T>(N[r]);
 
    // ---- exponential draws V(s,i) = -log(u)/gamma_i -------------------------
    const double tiny = 1.0 / 9007199254740992.0;  // keeps log(u) finite at u = 0
    std::vector<double> Vd(samples * M);            // row-major: sample s, station i
    if (variant == MciVariant::Amci) {
        const std::size_t Ih = (samples + 1) / 2;
        for (std::size_t i = 0; i < M; ++i) {
            for (std::size_t s = 0; s < Ih; ++s) {
                double u = mc_uniform01(rng);
                if (u <= 0.0) u = tiny;
                Vd[s * M + i] = -std::log(u) / gamma[i];
                if (Ih + s < samples) Vd[(Ih + s) * M + i] = -std::log(1.0 - u) / gamma[i];
            }
        }
    } else if (variant == MciVariant::Lhsmci) {
        std::vector<std::size_t> perm(samples);
        for (std::size_t i = 0; i < M; ++i) {
            for (std::size_t s = 0; s < samples; ++s) perm[s] = s;
            for (std::size_t s = samples - 1; s > 0; --s) {
                const std::size_t k = static_cast<std::size_t>(mc_uniform_int(rng, s + 1));
                const std::size_t t = perm[s];
                perm[s] = perm[k];
                perm[k] = t;
            }
            for (std::size_t s = 0; s < samples; ++s) {
                double u = (static_cast<double>(perm[s]) + mc_uniform01(rng)) /
                           static_cast<double>(samples);
                if (u <= 0.0) u = tiny;
                Vd[s * M + i] = -std::log(u) / gamma[i];
            }
        }
    } else {
        for (std::size_t s = 0; s < samples; ++s) {
            for (std::size_t i = 0; i < M; ++i) {
                double u = mc_uniform01(rng);
                if (u <= 0.0) u = tiny;
                Vd[s * M + i] = -std::log(u) / gamma[i];
            }
        }
    }
 
    std::vector<double> lZ(samples);
    std::vector<T> V(M);
    for (std::size_t s = 0; s < samples; ++s) {
        double lz = lconst;
        for (std::size_t i = 0; i < M; ++i) {
            const double v = Vd[s * M + i];
            V[i] = num_traits<T>::from_double(v);
            lz -= (1.0 - gamma[i]) * v;
        }
        for (std::size_t r = 0; r < R; ++r) {
            if (N[r] == 0) continue;
            T inner = Zv[r];
            for (std::size_t i = 0; i < M; ++i) inner += V[i] * D(i, r);
            lz += N[r] * num_traits<T>::log_as_double(inner);
        }
        lZ[s] = lz;
    }
 
    double lG = mc_logmeanexp(lZ);
    if (!std::isfinite(lG)) {
        // Floating-point range exception: the reference falls back to the
        // largest single log-weight, which is the Laplace-style lower bound.
        double m = -std::numeric_limits<double>::infinity();
        for (double v : lZ)
            if (v > m) m = v;
        lG = m;
    }
    return {mc_exp<T>(lG), lG};
}
 
/** Reference defaults: 1e5 samples, the IMCI proposal. */
template <class T>
NcResult<T> pfqn_mci(const Matrix<T>& D, const std::vector<int>& N, const std::vector<T>& Z,
                     McRng& rng) {
    return pfqn_mci(D, N, Z, static_cast<std::size_t>(100000), MciVariant::Imci, rng);
}
 
}  // namespace pfqn
}  // namespace line
 
#endif  // LINE_API_PFQN_MCI_H