sn_set_service_coc.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_INFER_SN_SET_SERVICE_COC_H
#define LINE_API_INFER_SN_SET_SERVICE_COC_H
 
/**
 * @file
 * @ingroup api_infer
 * Port of matlab/src/api/infer/sn_set_service_coc.m: overwrite the service
 * process of one (station, class) pair of a NetworkStruct with a Markovian law
 * fitted to a rate and an SCV, the way the inference routines (infer_fmlps,
 * infer_mlps and the Gibbs samplers) vary a demand without rebuilding the
 * model.
 *
 * The reference writes `sn.rates`, `sn.scv` and then, for a finite positive
 * rate, the representation `sn.proc`, `sn.procid`, `sn.phases`, `sn.phasessz`,
 * `sn.phaseshift`, `sn.mu`, `sn.phi` and `sn.pie`. The C++ NetworkStruct
 * derives every one of those from the stored `Distrib` on demand (`procid`,
 * `phases_of`, `phasessz_of`, `dist_to_map`), so the port stores the fitted
 * `Distrib` through `set_service` and writes `rates`/`scv` directly, which is
 * the same state.
 *
 * The fit is the reference's: SCV within 1e-10 of one (or NaN) gives an
 * exponential; SCV below one an Erlang with k = max(1, ceil(1/scv)) phases
 * (map_erlang); SCV above one the two-phase hyperexponential of kpctoolbox's
 * map_hyperexp(mean, scv), whose balanced construction starts from p = 0.99 and
 * divides p by ten until a feasible pair of rates appears, falling back to an
 * exponential when none does.
 */
 
#include <cmath>
#include <cstddef>
#include <limits>
 
#include "line/lang/lang_types.h"
#include "line/lang/qn/network_struct.h"
#include "line/num/number.h"
#include "line/util/error.h"
 
namespace line {
namespace api {
 
namespace detail {
 
/**
 * kpctoolbox map_hyperexp(MEAN, SCV, p): the two-phase hyperexponential with
 * mixing probability p matching the mean and SCV, trying the two roots of the
 * second-moment equation and lowering p by a factor ten when neither is
 * feasible. Returns false where the reference returns an empty MAP.
 */
inline bool map_hyperexp_fit(double mean, double scv, double p, double* mu1_out, double* mu2_out,
                             double* p_out) {
    const double E2 = (1.0 + scv) * mean * mean;
    const double Delta = -4.0 * p * mean * mean + 4.0 * p * p * mean * mean + 2.0 * E2 * p -
                         2.0 * E2 * p * p;
    const double den = E2 * p - 2.0 * mean * mean;
    auto feasible = [](double mu1, double mu2, double pp) {
        return std::isfinite(mu1) && std::isfinite(mu2) && mu1 > 0.0 && mu2 > 0.0 &&
               pp >= 0.0 && pp <= 1.0;
    };
    if (Delta >= 0.0 && den != 0.0) {
        double mu2 = (-2.0 * mean + 2.0 * p * mean + std::sqrt(Delta)) / den;
        double mu1 = mu2 * p / (p - 1.0 + mean * mu2);
        if (feasible(mu1, mu2, p)) {
            *mu1_out = mu1;
            *mu2_out = mu2;
            *p_out = p;
            return true;
        }
        mu2 = (-2.0 * mean + 2.0 * p * mean - std::sqrt(Delta)) / den;
        mu1 = mu2 * p / (p - 1.0 + mean * mu2);
        if (feasible(mu1, mu2, p)) {
            *mu1_out = mu1;
            *mu2_out = mu2;
            *p_out = p;
            return true;
        }
    }
    if (p > 1e-6) return map_hyperexp_fit(mean, scv, p / 10.0, mu1_out, mu2_out, p_out);
    return false;
}
 
}  // namespace detail
 
/**
 * @brief Set the service process of (stationIdx, classIdx) from a rate and an
 *        SCV, the reference's sn_set_service_coc(sn, stationIdx, classIdx, rate, scv).
 *
 * @param sn         the network struct, modified in place
 * @param stationIdx 1-based station
 * @param classIdx   1-based class
 * @param rate       service rate; NaN, non-positive or infinite leaves the
 *                   representation untouched after writing rates/scv
 * @param scv        squared coefficient of variation, default 1
 */
template <class T>
void sn_set_service_coc(qn::NetworkStruct<T>& sn, std::size_t stationIdx, std::size_t classIdx,
                        const T& rate, const T& scv) {
    if (stationIdx == 0 || stationIdx > sn.nstations)
        throw InputError("sn_set_service_coc: station index out of range");
    if (classIdx == 0 || classIdx > sn.nclasses)
        throw InputError("sn_set_service_coc: class index out of range");
    const std::size_t i = stationIdx - 1, r = classIdx - 1;
    sn.rates(i, r) = rate;
    sn.scv(i, r) = scv;
    const double rd = num_traits<T>::to_double(rate);
    const double sd = num_traits<T>::to_double(scv);
    if (std::isnan(rd) || rd <= 0.0 || std::isinf(rd)) return;
 
    const T one = num_traits<T>::from_int(1);
    const T meanVal = T(one / rate);
    lang::Distrib<T> d;
    if (std::isnan(sd) || std::fabs(sd - 1.0) < 1e-10) {
        d = lang::Distrib<T>::exp_mean(meanVal);
    } else if (sd < 1.0) {
        const double kd = std::ceil(1.0 / sd);
        const std::size_t k = static_cast<std::size_t>(kd < 1.0 ? 1.0 : kd);
        // map_erlang(mean, k): k phases of rate k/mean
        d = lang::Distrib<T>::erlang(T(num_traits<T>::from_int(static_cast<long>(k)) / meanVal), k);
    } else {
        double mu1 = 0.0, mu2 = 0.0, p = 0.99;
        if (detail::map_hyperexp_fit(num_traits<T>::to_double(meanVal), sd, 0.99, &mu1, &mu2, &p)) {
            d = lang::Distrib<T>::hyperexp(num_traits<T>::from_double(p),
                                           num_traits<T>::from_double(mu1),
                                           num_traits<T>::from_double(mu2));
        } else {
            d = lang::Distrib<T>::exp_mean(meanVal);
        }
    }
    sn.set_service(stationIdx, classIdx, d);
    if (i < sn.disabled.size() && r < sn.disabled[i].size()) sn.disabled[i][r] = false;
}
 
/** sn_set_service_coc with the reference's default scv = 1. */
template <class T>
void sn_set_service_coc(qn::NetworkStruct<T>& sn, std::size_t stationIdx, std::size_t classIdx,
                        const T& rate) {
    sn_set_service_coc(sn, stationIdx, classIdx, rate, num_traits<T>::from_int(1));
}
 
}  // namespace api
}  // namespace line
 
#endif  // LINE_API_INFER_SN_SET_SERVICE_COC_H