ldes_busyperiod.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_SOLVERS_LDES_LDES_BUSYPERIOD_H
#define LINE_SOLVERS_LDES_LDES_BUSYPERIOD_H
 
/**
 * @file
 * @ingroup line_solvers
 * Busy period measurement for the native LDES engine.
 *
 * A busy period of ORDER n for a set of stations runs from the instant an
 * arrival raises the jobs the set holds to n, up to the instant it falls back
 * below n (Daduna, J. ACM 35(3), 1988). Order 1 is the ordinary busy period;
 * higher orders measure how long the set stays at least n deep, which is what
 * the analytical reference `pfqn_busyp` returns.
 *
 * ONE TARGET PER (set, class filter). Every station is a set of its own, every
 * declared subnetwork is a set, and each is measured both aggregated over
 * classes and per class, exactly as `initializeBusyPeriods` builds them.
 *
 * THE EXIT IS DEFERRED TO THE END OF THE INSTANT, and that is the whole
 * difficulty. A job moving from one member of a target to another DROPS the
 * count and RESTORES it at the same simulated instant; treating the drop as an
 * exit ends a busy period that never ended and starts a new one, which roughly
 * HALVES the reported mean on any multi-station target. So a drop only records
 * a PENDING exit, an arrival at the same instant cancels it, and the pending
 * exits are committed once the clock actually advances. A period of zero
 * duration is an artifact of a simultaneous pair of events, not an
 * observation, and is dropped.
 */
 
#include <cmath>
#include <cstddef>
#include <limits>
#include <string>
#include <vector>
 
namespace line {
namespace ldes {
namespace engine {
 
/** One measured target: a station set, optionally filtered to one class. */
struct BpTarget {
    std::vector<std::size_t> stations;  ///< 0-based station indexes
    int job_class = -1;                 ///< -1 aggregates every class
    std::string name;
};
 
/**
 * The busy-period tracker of one run.
 *
 * `orders` is the highest order measured; 0 disables the measurement and every
 * entry point below returns immediately, so the engine pays nothing for it.
 */
class BusyPeriods {
public:
    BusyPeriods() {}
 
    void init(const std::vector<BpTarget>& targets, int orders, std::size_t nstations) {
        orders_ = orders;
        if (orders_ <= 0) return;
        targets_ = targets;
        const std::size_t nt = targets_.size();
        const std::size_t no = static_cast<std::size_t>(orders_);
        jobs_.assign(nt, 0);
        entry_.assign(nt, std::vector<double>(no, nan()));
        pending_.assign(nt, std::vector<double>(no, nan()));
        sum_.assign(nt, std::vector<double>(no, 0.0));
        count_.assign(nt, std::vector<double>(no, 0.0));
        by_station_.assign(nstations, std::vector<std::size_t>());
        for (std::size_t ti = 0; ti < nt; ++ti)
            for (std::size_t s : targets_[ti].stations)
                if (s < nstations) by_station_[s].push_back(ti);
    }
 
    bool enabled() const { return orders_ > 0; }
 
    /**
     * Record a population change of `delta` (+1 or -1) for `cls` at `station`.
     *
     * `now` must be non-decreasing across calls: the first call at a strictly
     * later instant is what commits the exits pending from the previous one.
     */
    void track(std::size_t station, std::size_t cls, int delta, double now) {
        if (orders_ <= 0 || station >= by_station_.size()) return;
        if (now > last_time_) {
            commit();
            last_time_ = now;
        }
        for (std::size_t ti : by_station_[station]) {
            const BpTarget& t = targets_[ti];
            if (t.job_class >= 0 && static_cast<std::size_t>(t.job_class) != cls) continue;
            jobs_[ti] += delta;
            // A rise to n opens the period of order n; a fall FROM n closes it,
            // and after the decrement the count is n-1, hence the +1.
            const int order = (delta > 0) ? jobs_[ti] : jobs_[ti] + 1;
            if (order < 1 || order > orders_) continue;
            const std::size_t oi = static_cast<std::size_t>(order - 1);
            if (delta > 0) {
                if (std::isnan(pending_[ti][oi]))
                    entry_[ti][oi] = now;
                else
                    // The level came back within the same instant: the drop was
                    // a move inside the target and the period never ended.
                    pending_[ti][oi] = nan();
            } else {
                pending_[ti][oi] = now;
                pending_list_.push_back(ti);
                pending_list_.push_back(oi);
            }
        }
    }
 
    /** Close every exit no longer contradicted by a same-instant arrival. */
    void commit() {
        for (std::size_t p = 0; p + 1 < pending_list_.size(); p += 2) {
            const std::size_t ti = pending_list_[p], oi = pending_list_[p + 1];
            const double exit = pending_[ti][oi];
            if (std::isnan(exit)) continue;
            const double entry = entry_[ti][oi];
            // A NaN entry belongs to a period straddling the warmup reset: its
            // start is unknown, so the period is dropped rather than guessed.
            if (!std::isnan(entry) && exit > entry) {
                sum_[ti][oi] += exit - entry;
                count_[ti][oi] += 1.0;
            }
            entry_[ti][oi] = nan();
            pending_[ti][oi] = nan();
        }
        pending_list_.clear();
    }
 
    /**
     * Discard everything observed during the warmup.
     *
     * A period IN PROGRESS at the reset has an unknown start, so its entry is
     * invalidated and the period is dropped when it ends -- counting it from
     * the reset instant would report a truncated period as a whole one and
     * bias every mean downwards.
     */
    void reset() {
        if (orders_ <= 0) return;
        for (std::size_t ti = 0; ti < targets_.size(); ++ti)
            for (std::size_t oi = 0; oi < static_cast<std::size_t>(orders_); ++oi) {
                sum_[ti][oi] = 0.0;
                count_[ti][oi] = 0.0;
                entry_[ti][oi] = nan();
                pending_[ti][oi] = nan();
            }
        pending_list_.clear();
    }
 
    const std::vector<BpTarget>& targets() const { return targets_; }
    int orders() const { return orders_; }
    /** Mean duration of order `oi`+1 at target `ti`; 0 when nothing was observed. */
    double mean(std::size_t ti, std::size_t oi) const {
        return (count_[ti][oi] > 0.0) ? sum_[ti][oi] / count_[ti][oi] : 0.0;
    }
    double count(std::size_t ti, std::size_t oi) const { return count_[ti][oi]; }
 
private:
    static double nan() { return std::numeric_limits<double>::quiet_NaN(); }
 
    int orders_ = 0;
    std::vector<BpTarget> targets_;
    std::vector<int> jobs_;
    std::vector<std::vector<double>> entry_, pending_, sum_, count_;
    std::vector<std::vector<std::size_t>> by_station_;
    std::vector<std::size_t> pending_list_;
    double last_time_ = 0.0;
};
 
}  // namespace engine
}  // namespace ldes
}  // namespace line
 
#endif  // LINE_SOLVERS_LDES_LDES_BUSYPERIOD_H