ldes_warm_start.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_SOLVERS_WRAPPERS_LDES_LDES_WARM_START_H
#define LINE_SOLVERS_WRAPPERS_LDES_LDES_WARM_START_H
 
/**
 * @file
 * @ingroup line_solvers
 * The LDES warm start: `SolverLDES(model, initSolver)` / `initFromSolver`.
 *
 * SOURCE. `@@NetworkSolver/warmStartPlacement.m` and
 * `@@SolverLDES/initFromSolver.m`. An auxiliary solver's steady state decides
 * an integer job placement (nstations x nclasses), which reaches the engine as
 * `LdesOptions.init_sol` (`--initsol`, STATION-MAJOR) with the warmup filter
 * off, since a run that starts near steady state has no transient to remove.
 *
 * TWO ARMS, as in the reference, and the caller picks one by which function it
 * calls rather than by the dynamic type of a solver object:
 *
 *  - `warm_start_placement_from_qlen`: any solver's mean queue lengths, floored
 *    per station and completed by LARGEST REMAINDER so each closed population is
 *    conserved; an open class is rounded and not completed.
 *  - `warm_start_placement_from_ctmc`: the MODE of the exact stationary law
 *    aggregated over per-(station, class) job counts. Ties go to the
 *    lexicographically smallest aggregate state, which is the row `unique` puts
 *    first and `max` then returns in the reference.
 *
 * Only a Queue or a Delay receives jobs; every other station keeps a zero row.
 */
 
#include <cmath>
#include <cstddef>
#include <map>
#include <vector>
 
#include "line/lang/qn/network_struct.h"
#include "line/solvers/ctmc/solver_ctmc.h"
#include "line/solvers/ctmc/solver_ctmc_analyzer.h"
#include "line/solvers/wrappers/ldes/ldes_options.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace ldes {
 
namespace warm_detail {
 
template <class T>
bool is_service_station(const qn::NetworkStruct<T>& sn, std::size_t i) {
    const lang::NodeType nt = sn.stations[i].nodetype;
    return nt == lang::NodeType::Queue || nt == lang::NodeType::Delay;
}
 
}  // namespace warm_detail
 
/**
 * Port of `placementFromMeanQLen`: an integer placement from mean queue lengths.
 *
 * @param sn the refreshed network struct the means were computed on
 * @param QN (nstations x nclasses) steady-state mean queue lengths
 */
template <class T>
Matrix<double> warm_start_placement_from_qlen(const qn::NetworkStruct<T>& sn,
                                              const Matrix<T>& QN) {
    const std::size_t M = sn.nstations, K = sn.nclasses;
    if (QN.rows() != M || QN.cols() != K)
        throw InputError("warm_start_placement_from_qlen: QN must be (nstations x nclasses)");
    Matrix<double> P(M, K, 0.0);
    for (std::size_t r = 0; r < K; ++r) {
        const double njobs = sn.classes[r].population;
        const bool closed = std::isfinite(njobs);
        std::vector<double> floors(M, 0.0), fracs(M, -1.0);
        std::vector<bool> eligible(M, false);
        for (std::size_t i = 0; i < M; ++i) {
            eligible[i] = warm_detail::is_service_station(sn, i);
            if (!eligible[i]) continue;
            const double m = std::max(0.0, num_traits<T>::to_double(QN(i, r)));
            if (closed) {
                floors[i] = std::floor(m);
                fracs[i] = m - floors[i];
            } else {
                floors[i] = std::round(m);
            }
            P(i, r) = floors[i];
        }
        if (!closed) continue;
        double placed = 0.0;
        for (std::size_t i = 0; i < M; ++i) placed += floors[i];
        long residual = std::lround(njobs) - std::lround(placed);
        while (residual > 0) {
            // `max` returns the FIRST maximum, so ties go to the lower station.
            std::size_t bi = 0;
            for (std::size_t i = 1; i < M; ++i)
                if (fracs[i] > fracs[bi]) bi = i;
            if (fracs[bi] < 0.0) {
                // Every remainder is spent: the leftover jobs go to the reference
                // station (or the first eligible one) so the population holds.
                std::size_t ref = sn.classes[r].refstat;  // 1-based
                std::size_t at = M;
                if (ref >= 1 && ref <= M && eligible[ref - 1]) at = ref - 1;
                for (std::size_t i = 0; at == M && i < M; ++i)
                    if (eligible[i]) at = i;
                if (at == M)
                    throw InputError("warm_start_placement_from_qlen: closed class '" +
                                     sn.classes[r].name + "' has no Queue or Delay to start at");
                P(at, r) += static_cast<double>(residual);
                break;
            }
            P(bi, r) += 1.0;
            fracs[bi] = -1.0;
            --residual;
        }
    }
    return P;
}
 
/**
 * Port of `placementFromCtmcSteadyState`: the mode of the aggregate stationary law.
 *
 * @param sn  the refreshed network struct
 * @param opt the CTMC options of the auxiliary solve (cutoff, method)
 */
template <class T>
Matrix<double> warm_start_placement_from_ctmc(const qn::NetworkStruct<T>& sn,
                                              const ctmc::CtmcOptions& opt) {
    const ctmc::CtmcSolution<T> sol = ctmc::solver_ctmc_analyzer(sn, opt);
    // A fork-join chain lives on the augmented struct, whose classes are not the
    // model's; the aggregate state would be read against the wrong columns.
    if (!sol.fjclassmap.empty())
        throw UnsupportedError(
            "warm_start_placement_from_ctmc: a fork-join model is solved on its tag-augmented "
            "chain; warm-start it from a mean-value solver instead");
    const std::size_t M = sn.nstations, K = sn.nclasses;
    const Matrix<T> SSq = ctmc::ctmc_state_space_aggr(sn, sol.chain.space);
    std::map<std::vector<long>, double> aggr;
    for (std::size_t s = 0; s < SSq.rows(); ++s) {
        double p = num_traits<T>::to_double(sol.pi[s]);
        if (p < 1e-14) p = 0.0;  // GlobalConstants.Zero, as the reference clips it
        std::vector<long> key(M * K);
        for (std::size_t c = 0; c < M * K; ++c)
            key[c] = std::lround(num_traits<T>::to_double(SSq(s, c)));
        aggr[key] += p;
    }
    if (aggr.empty())
        throw InputError("warm_start_placement_from_ctmc: the chain has no states");
    // Ordered keys, strict comparison: the first (smallest) of several modes wins.
    std::map<std::vector<long>, double>::const_iterator best = aggr.begin();
    for (std::map<std::vector<long>, double>::const_iterator it = aggr.begin(); it != aggr.end();
         ++it)
        if (it->second > best->second) best = it;
    Matrix<double> P(M, K, 0.0);
    for (std::size_t i = 0; i < M; ++i) {
        if (!warm_detail::is_service_station(sn, i)) continue;
        for (std::size_t r = 0; r < K; ++r) P(i, r) = static_cast<double>(best->first[i * K + r]);
    }
    return P;
}
 
/**
 * Port of `initFromSolver`'s last step: the placement becomes `init_sol`,
 * STATION-MAJOR (`reshape(placement', 1, M*K)`), and the transient filter is
 * disabled with `tranfilter = fixed`, `warmupfrac = 0`.
 */
inline void init_from_placement(LdesOptions& o, const Matrix<double>& P) {
    o.init_sol.clear();
    for (std::size_t i = 0; i < P.rows(); ++i)
        for (std::size_t r = 0; r < P.cols(); ++r) o.init_sol.push_back(P(i, r));
    o.tranfilter = "fixed";
    o.warmupfrac = 0.0;
}
 
}  // namespace ldes
}  // namespace line
 
#endif  // LINE_SOLVERS_WRAPPERS_LDES_LDES_WARM_START_H