sn_schedule_nominal.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_SN_SN_SCHEDULE_NOMINAL_H
#define LINE_API_SN_SN_SCHEDULE_NOMINAL_H
 
/**
 * @file
 * @ingroup api_sn
 * Port of `matlab/src/api/sn/sn_schedule_nominal.m`: unpack the MAPt or PHt slot
 * of `sn.proc` at one (station, class) into the per-segment (D0, D1) pairs and
 * their width-weighted time average.
 *
 * IT IS AN ACCESSOR HERE AND A DECODER THERE, and the difference is where the
 * conversion happens. MATLAB stores the slot as `{breakpoints, A, B, cyclic}`
 * with A and B holding (D0, D1) for a MAPt and (alpha, S) for a PHt, so every
 * reader must know which family it is looking at and rebuild the MAP pair
 * (S, (-S e) alpha) itself. This port converts a PHt ONCE, in
 * `Distrib::pht`, so what is stored is always the MAP schedule and every reader
 * sees one representation. The nominal pair is likewise computed once, in the
 * constructor, and lives in `Distrib::D0` / `D1` where a time-blind consumer --
 * the phase count, `sn.rates`, the fluid layout -- already looks for it.
 *
 * What is left for this function is the CONTRACT: return the same eight outputs
 * the reference returns (the width-weighted nominal pair, recomputed from the
 * segments as the reference computes it, breakpoints, the per-segment MAP
 * pairs, cyclic, and the marked / batch segment blocks), and refuse a (station, class) that carries no schedule
 * rather than inventing a one-segment one. A caller that wants the nominal of an
 * ordinary process reads `sn.service[i][r].D0` directly, which is what the
 * reference's own `nomD0`/`nomD1` fallback does.
 */
 
#include <cstddef>
#include <vector>
 
#include "line/lang/lang_types.h"
#include "line/lang/qn/network_struct.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace sn {
 
/** What `sn_schedule_nominal` returns, in the reference's own output order. */
template <class T>
struct ScheduleNominal {
    Matrix<T> D0bar, D1bar;             ///< the width-weighted time average
    std::vector<T> breakpoints;         ///< the boundary vector, nseg + 1 long
    std::vector<Matrix<T> > segD0;      ///< per-segment D0, already in MAP form
    std::vector<Matrix<T> > segD1;      ///< per-segment D1
    bool cyclic = false;
    /**
     * The reference's 7th and 8th outputs: the per-mark blocks of a marked
     * schedule, `segMark[c][k]` being mark c+1 in segment k, and the per-batch
     * blocks of a batch schedule, `segBatch[c][b][k]`. Empty on an unmarked
     * schedule, as `segMark = {}` / `segBatch = {}` are there.
     */
    std::vector<std::vector<Matrix<T> > > segMark;
    std::vector<std::vector<std::vector<Matrix<T> > > > segBatch;
};
 
/** True when (ist, r) carries a MAPt / PHt / NHPP schedule. Both 0-based. */
template <class T>
bool sn_has_schedule(const qn::NetworkStruct<T>& sn, std::size_t ist, std::size_t r) {
    if (ist >= sn.service.size() || r >= sn.service[ist].size()) return false;
    return sn.service[ist][r].has_schedule();
}
 
/** Port of `sn_schedule_nominal(sn, ist, r)`; `ist` and `r` are 0-based here. */
template <class T>
ScheduleNominal<T> sn_schedule_nominal(const qn::NetworkStruct<T>& sn, std::size_t ist,
                                       std::size_t r) {
    if (!sn_has_schedule(sn, ist, r))
        throw InputError("sn_schedule_nominal: station " + std::to_string(ist + 1) + ", class " +
                         std::to_string(r + 1) +
                         " carries no MAPt/PHt schedule; read sn.service[i][r].D0 for an "
                         "ordinary process");
    const lang::Distrib<T>& d = sn.service[ist][r];
    ScheduleNominal<T> out;
    // The reference's width-weighted time average, D0bar = sum_k
    // (widths(k)/total) * segD0{k} with widths = diff(breakpoints); the
    // constructor stored the same pair in D0 / D1, which is the fallback when
    // the breakpoint vector does not frame the segments.
    const std::size_t nseg = d.sched_D0.size();
    if (d.sched_bp.size() == nseg + 1 && d.sched_D1.size() == nseg && nseg > 0) {
        const T zero = num_traits<T>::from_int(0);
        T total = zero;
        for (std::size_t k = 0; k < nseg; ++k) total += d.sched_bp[k + 1] - d.sched_bp[k];
        out.D0bar = Matrix<T>(d.sched_D0[0].rows(), d.sched_D0[0].cols(), zero);
        out.D1bar = Matrix<T>(d.sched_D1[0].rows(), d.sched_D1[0].cols(), zero);
        for (std::size_t k = 0; k < nseg; ++k) {
            const T w = T((d.sched_bp[k + 1] - d.sched_bp[k]) / total);
            for (std::size_t a = 0; a < out.D0bar.rows(); ++a)
                for (std::size_t b = 0; b < out.D0bar.cols(); ++b)
                    out.D0bar(a, b) = T(out.D0bar(a, b) + w * d.sched_D0[k](a, b));
            for (std::size_t a = 0; a < out.D1bar.rows(); ++a)
                for (std::size_t b = 0; b < out.D1bar.cols(); ++b)
                    out.D1bar(a, b) = T(out.D1bar(a, b) + w * d.sched_D1[k](a, b));
        }
    } else {
        out.D0bar = d.D0;
        out.D1bar = d.D1;
    }
    out.breakpoints = d.sched_bp;
    out.segD0 = d.sched_D0;
    out.segD1 = d.sched_D1;
    out.cyclic = d.sched_cyclic;
    out.segMark = d.sched_Dmark;
    out.segBatch = d.sched_Dbatch;
    return out;
}
 
}  // namespace sn
}  // namespace line
 
#endif  // LINE_API_SN_SN_SCHEDULE_NOMINAL_H