LINE Solver (C++)
Templated C++ port of the LINE queueing solver
Loading...
Searching...
No Matches
lang_types.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
3 * All rights reserved.
4 */
5#ifndef LINE_LANG_LANG_TYPES_H
6#define LINE_LANG_LANG_TYPES_H
7
8/**
9 * @file
10 * @ingroup line_lang
11 * Enumerations and the minimal distribution descriptor shared by the model
12 * layer of the C++ port.
13 *
14 * The numeric values are the MATLAB ones (matlab/src/lang/constant), not a
15 * fresh numbering, because they cross the JSON boundary and appear in the
16 * dumps used as parity oracles. _kb/11 records that MATLAB and Python disagree
17 * on ProcessType numbering (MATLAB starts at EXP=0, Python at EXP=1); this port
18 * follows MATLAB, the reference implementation, so a numeric comparison against
19 * a MATLAB dump is meaningful and one against a Python dump is not -- compare
20 * by name there.
21 *
22 * SCOPE: the model layer added here exists to run SolverLN over a layered
23 * queueing network whose layers are solved by SolverMVA. It carries the
24 * scheduling disciplines, node kinds and precedence types that path reaches
25 * and REFUSES the rest by name rather than silently mapping them onto a
26 * neighbour, because a discipline that is silently treated as FCFS returns a
27 * plausible number that is wrong.
28 */
29
30#include <cmath>
31#include <cstddef>
32#include <functional>
33#include <limits>
34#include <memory>
35#include <string>
36#include <vector>
37
38#include "line/num/number.h"
39#include "line/util/error.h"
40#include "line/util/matrix.h"
41
42namespace line {
43namespace lang {
44
45/** Solver output metrics, with the numeric values of MATLAB `MetricType`. */
46enum class MetricType {
47 ResidT = 0,
48 RespT = 1,
49 DropRate = 2,
50 QLen = 3,
51 QueueT = 4,
52 FCRWeight = 5,
53 FCRMemOcc = 6,
54 FJQLen = 7,
55 FJRespT = 8,
56 RespTSink = 9,
57 SysQLen = 10,
58 SysRespT = 11,
59 SysTput = 12,
60 Tput = 13,
61 ArvR = 14,
62 TputSink = 15,
63 Util = 16,
64 TranQLen = 17,
65 TranUtil = 18,
66 TranTput = 19,
67 TranRespT = 20,
68 Tard = 21,
69 SysTard = 22
70};
71
72/** Port of `MetricType.toText`. */
73inline const char* metric_to_text(MetricType metric) {
74 switch (metric) {
75 case MetricType::ResidT: return "Residence Time";
76 case MetricType::RespT: return "Response Time";
77 case MetricType::DropRate: return "Drop Rate";
78 case MetricType::QLen: return "Number of Customers";
79 case MetricType::QueueT: return "Queue Time";
80 case MetricType::FCRWeight: return "FCR Total Weight";
81 case MetricType::FCRMemOcc: return "FCR Memory Occupation";
82 case MetricType::FJQLen: return "Fork Join Number of Customers";
83 case MetricType::FJRespT: return "Fork Join Response Time";
84 case MetricType::RespTSink: return "Response Time per Sink";
85 case MetricType::SysQLen: return "System Number of Customers";
86 case MetricType::SysRespT: return "System Response Time";
87 case MetricType::SysTput: return "System Throughput";
88 case MetricType::Tput: return "Throughput";
89 case MetricType::ArvR: return "Arrival Rate";
90 case MetricType::TputSink: return "Throughput per Sink";
91 case MetricType::Util: return "Utilization";
92 case MetricType::TranQLen: return "Tran Number of Customers";
93 case MetricType::TranUtil: return "Tran Utilization";
94 case MetricType::TranTput: return "Tran Throughput";
95 case MetricType::TranRespT: return "Tran Response Time";
96 case MetricType::Tard: return "Tardiness";
97 case MetricType::SysTard: return "System Tardiness";
98 default: return "Unknown Metric";
99 }
100}
101
102/**
103 * The events a state can undergo, with the values of MATLAB EventType.
104 *
105 * An event is ACTIVE at the node that schedules it and PASSIVE at the node
106 * that receives it: a DEP at one station is the ARV at the next, and only the
107 * active half carries a rate. The passive half is marked with a rate of -1,
108 * which the generator assembly replaces with the active rate -- a convention
109 * that only reads as a sentinel because a rate can never be negative.
110 */
111enum class EventType {
112 INIT = -1, ///< the model is initialized, t = 0
113 LOCAL = 0, ///< dummy event, no state change outside the node
114 ARV = 1, ///< a job arrives
115 DEP = 2, ///< a job departs
116 PHASE = 3, ///< service advances a phase WITHOUT departing
117 READ = 4, ///< a cache item is read
118 STAGE = 5, ///< a random environment changes stage
119 ENABLE = 6, ///< an SPN mode becomes enabled
120 FIRE = 7, ///< an SPN mode fires
121 PRE = 8, ///< consume from a place or queue buffer, no server effect
122 POST = 9, ///< produce to a place or queue buffer
123 RENEGE = 10, ///< a waiting job abandons the queue (impatience)
124 RETRY = 11, ///< an orbiting job retries entry at a retrial station
125 SWITCH = 12, ///< a polling server advances its switchover timer
126 FAILURE = 13, ///< the server breaks down, going from up to down
127 REPAIR = 14, ///< the server is repaired, going from down to up, resuming
128 ///< the held job, which is why it emits no START
129 START = 15, ///< a job begins or resumes holding a server
130 PREEMPT = 16 ///< a job holding a server is pushed back into the buffer
131};
132// START and PREEMPT are instantaneous tags on the arc of the ARV or DEP that
133// causes them, never the active half of an sn.sync entry: no clock, no state,
134// no change to any numerical result. PREEMPT is spelled in full because PRE
135// already names the Petri-net pre-arc.
136
137inline const char* event_to_text(EventType e) {
138 switch (e) {
139 case EventType::ARV: return "ARV";
140 case EventType::DEP: return "DEP";
141 case EventType::PHASE: return "PHASE";
142 case EventType::READ: return "READ";
143 case EventType::LOCAL: return "LOCAL";
144 case EventType::STAGE: return "STAGE";
145 case EventType::ENABLE: return "ENABLE";
146 case EventType::FIRE: return "FIRE";
147 case EventType::PRE: return "PRE";
148 case EventType::POST: return "POST";
149 case EventType::RENEGE: return "RENEGE";
150 case EventType::RETRY: return "RETRY";
151 case EventType::SWITCH: return "SWITCH";
152 case EventType::FAILURE: return "FAILURE";
153 case EventType::REPAIR: return "REPAIR";
154 case EventType::START: return "START";
155 case EventType::PREEMPT: return "PREEMPT";
156 default: return "INIT";
157 }
158}
159
160/**
161 * G-network signal classes, with the values of MATLAB SignalType.
162 *
163 * A signal is not a job: it never joins a station, it acts on the jobs already
164 * there and is annihilated. REPLY is the odd one out -- it completes a
165 * synchronous call and then joins as an ordinary job.
166 */
167enum class SignalType {
168 REPLY = 0, ///< completes a synchronous call, releasing a held server
169 NEGATIVE = 1, ///< removes a batch of jobs (Gelenbe's negative customer)
170 CATASTROPHE = 2 ///< removes EVERY job at the station
171};
172
173/** Which job a negative signal removes, with the values of MATLAB RemovalPolicy. */
174enum class RemovalPolicy {
175 RANDOM = 0, ///< uniform over waiting AND in-service jobs
176 FCFS = 1, ///< the oldest waiting job; servers only once nobody waits
177 LCFS = 2 ///< the newest waiting job; servers only once nobody waits
178};
179
180/** Scheduling disciplines, with the values of MATLAB SchedStrategy. */
181enum class SchedStrategy {
182 INF = 0,
183 FCFS = 1,
184 LCFS = 2,
185 SIRO = 3,
186 SJF = 4,
187 LJF = 5,
188 PS = 6,
189 DPS = 7,
190 GPS = 8,
191 SEPT = 9,
192 LEPT = 10,
193 HOL = 11,
194 FORK = 12,
195 EXT = 13,
196 REF = 14,
197 LCFSPR = 15,
198 POLLING = 16,
199 // The preemptive and priority families. PR resumes an interrupted job in
200 // the phase it held; PI restarts it from the entry phase, which is why the
201 // two cannot share an encoding: PR must carry the phase of EVERY preempted
202 // job, PI need not. FCFSPRIO is MATLAB's alias for HOL, so it is not a
203 // distinct enumerator here.
204 PSPRIO = 17,
205 DPSPRIO = 18,
206 GPSPRIO = 19,
207 LCFSPI = 20,
208 LCFSPRIO = 21,
209 LCFSPRPRIO = 22,
210 LCFSPIPRIO = 23,
211 FCFSPR = 24,
212 FCFSPI = 25,
213 FCFSPRPRIO = 26,
214 FCFSPIPRIO = 27,
215 SRPT = 28,
216 SRPTPRIO = 29,
217 EDD = 30,
218 EDF = 31,
219 LPS = 32,
220 PSJF = 33,
221 FB = 34,
222 LRPT = 35,
223 SETF = 36,
224 FSP = 37,
225 PAS = 38,
226 OI = 39,
227 NONE = -1
228};
229
230inline const char* sched_to_text(SchedStrategy s) {
231 switch (s) {
232 case SchedStrategy::INF: return "inf";
233 case SchedStrategy::FCFS: return "fcfs";
234 case SchedStrategy::LCFS: return "lcfs";
235 case SchedStrategy::SIRO: return "siro";
236 case SchedStrategy::SJF: return "sjf";
237 case SchedStrategy::LJF: return "ljf";
238 case SchedStrategy::PS: return "ps";
239 case SchedStrategy::DPS: return "dps";
240 case SchedStrategy::GPS: return "gps";
241 case SchedStrategy::SEPT: return "sept";
242 case SchedStrategy::LEPT: return "lept";
243 case SchedStrategy::HOL: return "hol";
244 case SchedStrategy::FORK: return "fork";
245 case SchedStrategy::EXT: return "ext";
246 case SchedStrategy::REF: return "ref";
247 case SchedStrategy::LCFSPR: return "lcfspr";
248 case SchedStrategy::POLLING: return "polling";
249 case SchedStrategy::SRPT: return "srpt";
250 case SchedStrategy::LPS: return "lps";
251 case SchedStrategy::PSJF: return "psjf";
252 case SchedStrategy::FB: return "fb";
253 case SchedStrategy::LRPT: return "lrpt";
254 case SchedStrategy::SETF: return "setf";
255 case SchedStrategy::FSP: return "fsp";
256 case SchedStrategy::PSPRIO: return "psprio";
257 case SchedStrategy::DPSPRIO: return "dpsprio";
258 case SchedStrategy::GPSPRIO: return "gpsprio";
259 case SchedStrategy::LCFSPI: return "lcfspi";
260 case SchedStrategy::LCFSPRIO: return "lcfsprio";
261 case SchedStrategy::LCFSPRPRIO: return "lcfsprprio";
262 case SchedStrategy::LCFSPIPRIO: return "lcfspiprio";
263 case SchedStrategy::FCFSPR: return "fcfspr";
264 case SchedStrategy::FCFSPI: return "fcfspi";
265 case SchedStrategy::FCFSPRPRIO: return "fcfsprprio";
266 case SchedStrategy::FCFSPIPRIO: return "fcfspiprio";
267 case SchedStrategy::SRPTPRIO: return "srptprio";
268 case SchedStrategy::EDD: return "edd";
269 case SchedStrategy::EDF: return "edf";
270 case SchedStrategy::PAS: return "pas";
271 case SchedStrategy::OI: return "oi";
272 default: return "none";
273 }
274}
275
276/** Parse the `scheduling` attribute of an .lqnx processor or task. */
277inline SchedStrategy sched_from_lqnx(const std::string& s) {
278 if (s == "inf" || s == "INF") return SchedStrategy::INF;
279 if (s == "fcfs" || s == "FCFS") return SchedStrategy::FCFS;
280 if (s == "ps" || s == "PS") return SchedStrategy::PS;
281 if (s == "ref" || s == "REF") return SchedStrategy::REF;
282 if (s == "hol" || s == "HOL") return SchedStrategy::HOL;
283 // lqns spells preemptive priority resume (SCHEDULE_PPR) "pri"; "pp" is the
284 // stale lqn-core.xsd spelling, absent from the lqns 6.2.31 sources. It is
285 // preemptive, so it is FCFSPRPRIO and not the non-preemptive HOL
286 if (s == "pri" || s == "PRI" || s == "pp") return SchedStrategy::FCFSPRPRIO;
287 if (s == "rand" || s == "siro") return SchedStrategy::SIRO;
288 if (s == "sjf") return SchedStrategy::SJF;
289 if (s == "ljf") return SchedStrategy::LJF;
290 if (s == "lcfs") return SchedStrategy::LCFS;
291 if (s == "burst" || s == "poll") return SchedStrategy::FCFS;
292 throw UnsupportedError("lqn reader: unsupported scheduling discipline '" + s + "'");
293}
294
295/**
296 * The `scheduling` attribute an .lqnx processor or task carries for a strategy.
297 *
298 * The inverse of sched_from_lqnx over the disciplines the schema spells, and a
299 * refusal by name for every other one. It refuses rather than falling back on
300 * fcfs because the file is handed to an external solver: a task written as fcfs
301 * when the model says lcfspr is answered, not rejected, and the discipline
302 * would be lost inside a number that looks ordinary.
303 */
304inline std::string sched_to_lqnx(SchedStrategy s) {
305 switch (s) {
306 case SchedStrategy::INF: return "inf";
307 case SchedStrategy::FCFS: return "fcfs";
308 case SchedStrategy::PS: return "ps";
309 case SchedStrategy::REF: return "ref";
310 case SchedStrategy::HOL: return "hol";
311 case SchedStrategy::FCFSPRPRIO: return "pri";
312 case SchedStrategy::SIRO: return "rand";
313 case SchedStrategy::SJF: return "sjf";
314 case SchedStrategy::LJF: return "ljf";
315 case SchedStrategy::LCFS: return "lcfs";
316 default:
317 throw UnsupportedError(
318 std::string("the LQN XML schema has no spelling for scheduling discipline '") +
319 sched_to_text(s) + "'; it carries inf, fcfs, ps, ref, hol, pri, rand, sjf, ljf and lcfs");
320 }
321}
322
323/** Node kinds, with the values of MATLAB NodeType. */
324enum class NodeType {
325 Queue = 0,
326 Source = 1,
327 Delay = 2,
328 ClassSwitch = 3,
329 Logger = 4,
330 Cache = 5,
331 Router = 6,
332 Fork = 7,
333 Place = 8,
334 Transition = 9,
335 Region = 10,
336 Join = 11,
337 Sink = 12
338};
339
340/** Name of a node kind, for diagnostics. The JSON spelling lives in the writer. */
341inline const char* node_type_to_text(NodeType t) {
342 switch (t) {
343 case NodeType::Queue: return "Queue";
344 case NodeType::Source: return "Source";
345 case NodeType::Delay: return "Delay";
346 case NodeType::ClassSwitch: return "ClassSwitch";
347 case NodeType::Logger: return "Logger";
348 case NodeType::Cache: return "Cache";
349 case NodeType::Router: return "Router";
350 case NodeType::Fork: return "Fork";
351 case NodeType::Place: return "Place";
352 case NodeType::Transition: return "Transition";
353 case NodeType::Region: return "Region";
354 case NodeType::Join: return "Join";
355 case NodeType::Sink: return "Sink";
356 }
357 return "Unknown";
358}
359
360/** SPN transition timing, with the values of MATLAB TimingStrategy. */
361enum class TimingStrategy {
362 TIMED = 0, ///< fires after its firing distribution elapses
363 IMMEDIATE = 1 ///< fires with zero delay, resolved by weight and priority
364};
365
366/** Job class kinds, with the values of MATLAB JobClassType. */
367enum class JobClassType { OPEN = 0, CLOSED = 1 };
368
369/** Polling service disciplines, with the values of MATLAB PollingType. */
370enum class PollingType {
371 GATED = 0, ///< serve exactly the jobs present at the polling instant
372 EXHAUSTIVE = 1, ///< serve until the queue empties
373 KLIMITED = 2, ///< serve at most K per visit (K in pollingPar)
374 DECREMENTING = 3 ///< serve until the queue is one shorter than at arrival
375};
376
377/** Cache replacement policies, with the values of MATLAB ReplacementStrategy. */
378enum class ReplacementStrategy {
379 RR = 0, ///< random replacement
380 FIFO = 1, ///< first in, first out
381 SFIFO = 2, ///< strict FIFO
382 LRU = 3, ///< least recently used
383 HLRU = 4, ///< h-LRU / LRU(m): h lists, promote i -> i+1 on a hit
384 CLIMB = 5, ///< move up one position on a hit (transposition rule)
385 QLRU = 6 ///< q-LRU: LRU with probabilistic admission on a miss
386};
387
388/** Routing strategies, with the values of MATLAB RoutingStrategy. */
389enum class RoutingStrategy {
390 RAND = 0,
391 PROB = 1,
392 RROBIN = 2,
393 WRROBIN = 3,
394 JSQ = 4,
395 FIRING = 5,
396 SQ = 6,
397 /** Krzesinski (1987) product-form state-dependent routing. */
398 SDR = 7,
399 DISABLED = -1
400};
401
402inline const char* routing_to_text(RoutingStrategy r) {
403 switch (r) {
404 case RoutingStrategy::RAND: return "rand";
405 case RoutingStrategy::PROB: return "prob";
406 case RoutingStrategy::RROBIN: return "rrobin";
407 case RoutingStrategy::WRROBIN: return "wrrobin";
408 case RoutingStrategy::JSQ: return "jsq";
409 case RoutingStrategy::FIRING: return "firing";
410 case RoutingStrategy::SQ: return "sq";
411 case RoutingStrategy::SDR: return "sdr";
412 default: return "disabled";
413 }
414}
415
416/**
417 * Blocking and loss rules, with the values of MATLAB DropStrategy.
418 *
419 * WAITQ is -1 and is also the marker `refreshCapacity` writes where the rule is
420 * never consulted (an unbounded station, or a closed class), so it means two
421 * different things depending on the station's capacity; see the comment in
422 * refresh_capacity().
423 */
424enum class DropStrategy {
425 WAITQ = -1,
426 DROP = 1,
427 BAS = 2,
428 BBS = 3,
429 RSRD = 4,
430 RETRIAL = 5,
431 RETRIAL_WITH_LIMIT = 6
432};
433
434/**
435 * Impatience kinds, with the values of MATLAB ImpatienceType.
436 *
437 * RENEGING is a timer a job started at a queue; BALKING is a decision taken
438 * BEFORE joining, on the state of the queue, and is parameterized by
439 * `Station::balking` rather than by a distribution; RETRIAL sends the job to an
440 * orbit and is parameterized by `RetrialParam`.
441 */
442enum class ImpatienceType { NONE = 0, RENEGING = 1, BALKING = 2, RETRIAL = 3 };
443
444/** Balking rules, with the values of MATLAB BalkingStrategy. */
445enum class BalkingStrategy { NONE = 0, QUEUE_LENGTH = 1, EXPECTED_WAIT = 2, COMBINED = 3 };
446
447/**
448 * How a heterogeneous station picks among its server types, MATLAB
449 * HeteroSchedPolicy. ORDER is the default: the declared order of the types.
450 */
451enum class HeteroSchedPolicy { ORDER = 0, ALIS = 1, ALFS = 2, FAIRNESS = 3, FSF = 4, RAIS = 5 };
452
453/**
454 * When a Place releases a served token, MATLAB DepartureDiscipline. NORMAL is
455 * the standard queueing-Petri-net rule (available on completion); FIFO holds
456 * it until every earlier arrival to the depository has been released.
457 */
458enum class DepartureDiscipline { NORMAL = 0, FIFO = 1 };
459
460/** Join rules, with the values of MATLAB JoinStrategy. */
461enum class JoinStrategy { STD = 1, PARTIAL = 2 };
462
463/** LQN element kinds, with the values of MATLAB LayeredNetworkElement. */
464enum class LqnElement { HOST = 0, TASK = 1, ENTRY = 2, ACTIVITY = 3, CALL = 4 };
465
466/** Call kinds, with the values of MATLAB CallType. */
467enum class CallType { NONE = 0, SYNC = 1, ASYNC = 2, FWD = 3 };
468
469/** Activity precedence kinds, with the values of MATLAB ActivityPrecedenceType. */
481
482/** Distribution kinds, with the values of MATLAB ProcessType. */
483enum class ProcessType {
484 EXP = 0,
485 ERLANG = 1,
486 HYPEREXP = 2,
487 PH = 3,
488 APH = 4,
489 MAP = 5,
490 UNIFORM = 6,
491 DET = 7,
492 COXIAN = 8,
493 GAMMA = 9,
494 PARETO = 10,
495 MMPP2 = 11,
496 REPLAYER = 12,
497 IMMEDIATE = 13,
498 DISABLED = 14,
499 COX2 = 15,
500 WEIBULL = 16,
501 LOGNORMAL = 17,
502 DUNIFORM = 18,
503 BERNOULLI = 19,
504 /**
505 * A `Prior`: a weighted set of ALTERNATIVE distributions, or a density over
506 * a scalar parameter plus a factory from it. It is not a mixture -- each
507 * alternative is a separate model realization -- and only SolverUQ consumes
508 * it; every other solver refuses it through Feature::Prior.
509 */
510 PRIOR = 20,
511 BINOMIAL = 21,
512 POISSON = 22,
513 GEOMETRIC = 23,
514 BMAP = 24,
515 ME = 25,
516 RAP = 26,
517 DISCRETESAMPLER = 27,
518 ZIPF = 28,
519 DMAP = 29,
520 MMAP = 31,
521 EMPIRICALCDF = 32,
522 /**
523 * The time-INHOMOGENEOUS families of Ko and Pender (ORL 45, 2017): an
524 * NHPP is a rate schedule lambda(t), a MAPt a (D0(t), D1(t)) schedule and a
525 * PHt an (alpha(t), S(t)) one, all piecewise constant on one breakpoint
526 * vector and optionally cyclic. The numeric values are MATLAB's
527 * (`ProcessType.m:41-43`).
528 *
529 * THEY CARRY A NOMINAL PAIR TOO. `Distrib::D0`/`D1` hold the width-weighted
530 * time average of the schedule, which is what `sn_schedule_nominal` returns
531 * as its first two outputs and what every consumer that has no notion of
532 * time -- the phase count, the rate, the fluid layout -- reads. The schedule
533 * itself lives in `sched_bp`/`sched_D0`/`sched_D1` beside it, and only a
534 * solver that integrates in time looks at it.
535 */
536 NHPP = 33,
537 MAPT = 34,
538 PHT = 35,
539 /**
540 * A Gaussian, and the ONE family whose value is not MATLAB's, because
541 * MATLAB has none to copy: `ProcessType.m` stops at 35 and `Normal.m` is a
542 * `ContinuousDistribution` with no id at all, exactly as `Normal.java` and
543 * the python `Normal` have none. That is not an oversight in the reference
544 * -- a Gaussian has mass below zero, so it can never be a service or
545 * interarrival process and can never appear in `sn.proc`. It reaches this
546 * port only as the PARAMETER density of a continuous `Prior`, where it is
547 * read through `dist_cdf` and `dist_quantile` and never through
548 * `dist_to_map`.
549 *
550 * The value is deliberately far outside 0..35 so that it can never collide
551 * with an id MATLAB assigns later; `sn.procid` must never carry it, and
552 * `dist_to_map` refuses it by name rather than handing back the Erlang fit
553 * its default arm would otherwise produce for a law with negative support.
554 */
555 NORMAL = 100,
556 NONE = -1
557};
558
559/** The MATLAB ProcessType name, as `sn.procid` prints it. */
560inline const char* process_to_text(ProcessType p) {
561 switch (p) {
562 case ProcessType::EXP: return "Exp";
563 case ProcessType::ERLANG: return "Erlang";
564 case ProcessType::HYPEREXP: return "HyperExp";
565 case ProcessType::PH: return "PH";
566 case ProcessType::APH: return "APH";
567 case ProcessType::MAP: return "MAP";
568 case ProcessType::UNIFORM: return "Uniform";
569 case ProcessType::DET: return "Det";
570 case ProcessType::COXIAN: return "Coxian";
571 case ProcessType::GAMMA: return "Gamma";
572 case ProcessType::PARETO: return "Pareto";
573 case ProcessType::MMPP2: return "MMPP2";
574 case ProcessType::REPLAYER: return "Replayer";
575 case ProcessType::IMMEDIATE: return "Immediate";
576 case ProcessType::DISABLED: return "Disabled";
577 case ProcessType::COX2: return "Cox2";
578 case ProcessType::WEIBULL: return "Weibull";
579 case ProcessType::LOGNORMAL: return "Lognormal";
580 case ProcessType::DUNIFORM: return "DiscreteUniform";
581 case ProcessType::BERNOULLI: return "Bernoulli";
582 case ProcessType::BINOMIAL: return "Binomial";
583 case ProcessType::POISSON: return "Poisson";
584 case ProcessType::GEOMETRIC: return "Geometric";
585 case ProcessType::BMAP: return "BMAP";
586 case ProcessType::ME: return "ME";
587 case ProcessType::RAP: return "RAP";
588 case ProcessType::DISCRETESAMPLER: return "DiscreteSampler";
589 case ProcessType::ZIPF: return "Zipf";
590 case ProcessType::DMAP: return "DMAP";
591 case ProcessType::MMAP: return "MMAP";
592 case ProcessType::EMPIRICALCDF: return "EmpiricalCDF";
593 case ProcessType::PRIOR: return "Prior";
594 case ProcessType::NHPP: return "NHPP";
595 case ProcessType::MAPT: return "MAPt";
596 case ProcessType::PHT: return "PHt";
597 case ProcessType::NORMAL: return "Normal";
598 default: return "none";
599 }
600}
601
602/**
603 * `ProcessType.isMarkovian`: true when `sn.proc` carries an exact matrix
604 * representation of the law, rather than the Erlang fit `convertToMAP` leaves
605 * there for the parameter-only families. ME and RAP count -- their
606 * representation is the matrix-exponential analogue, not a generator -- and the
607 * discrete families do not, so a solver reading `sn.proc` as the law must gate
608 * on this exactly as MATLAB `ProcessType.m` does.
609 */
610inline bool process_is_markovian(ProcessType p) {
611 switch (p) {
612 case ProcessType::EXP:
613 case ProcessType::ERLANG:
614 case ProcessType::HYPEREXP:
615 case ProcessType::PH:
616 case ProcessType::APH:
617 case ProcessType::MAP:
618 case ProcessType::COXIAN:
619 case ProcessType::COX2:
620 case ProcessType::MMPP2:
621 case ProcessType::ME:
622 case ProcessType::RAP:
623 case ProcessType::DMAP:
624 case ProcessType::BMAP:
625 case ProcessType::MMAP: return true;
626 default: return false;
627 }
628}
629
630/**
631 * A class-dependent scaling map, `sn.cdscaling`.
632 *
633 * It takes the per-class population vector at one station and returns the
634 * per-class rate multipliers, which is the signature `pfqn_cdfun` consumes; the
635 * alias resolves to the same std::function type as `pfqn::CdScaling`, so a map
636 * built here is passed straight through to the api layer.
637 */
638template <class T>
639using CdScaling = std::function<std::vector<T>(const std::vector<T>&)>;
640
641/**
642 * A globally state-dependent scaling, `sn.gdscaling`.
643 *
644 * Unlike CdScaling it is declared on the NETWORK, not on a station: the argument
645 * is the FULL population matrix, given row-major as nstations rows of nclasses
646 * entries, and the result is either one scalar, one entry per station, or one
647 * entry per (station, class) in the same row-major order. This is the Whittle
648 * primitive -- a rate that reads the whole state -- and no per-station scaling
649 * can express it when one route holds several resources at once.
650 */
651template <class T>
652using GdScaling = std::function<std::vector<T>(const std::vector<T>&)>;
653
654// ---------------------------------------------------------------------------
655// Global constants
656// ---------------------------------------------------------------------------
657
658/**
659 * The MATLAB GlobalConstants, as reported by lineStart at its defaults.
660 *
661 * These are doubles on purpose even in the exact instantiation: they are
662 * tolerances and sentinels of the ALGORITHM, not quantities of the model, and
663 * `lineStart` prints exactly these values. Converting them through
664 * num_traits<T>::from_double keeps the exact backend reproducing the same
665 * branch decisions as the reference rather than a mathematically cleaner set.
666 */
667struct GlobalConstants {
668 static constexpr double FineTol = 1e-8;
669 static constexpr double CoarseTol = 1e-3;
670 static constexpr double Zero = 1e-14;
671 /** Below this an off-diagonal entry is NO ARC of the phase / state graph. */
672 static constexpr double ArcTol = 1e-12;
673 /** Rate of an Immediate distribution; its mean is 1/Immediate = 1e-8. */
674 static constexpr double Immediate = 1e8;
675 /**
676 * Stand-in for an unbounded COUNT, MATLAB `GlobalConstants.MaxInt`. Used
677 * where a state row must hold a server count and Inf is not a count.
678 */
679 static constexpr double MaxInt = 2147483647.0;
680};
681
682// ---------------------------------------------------------------------------
683// Distribution descriptor
684// ---------------------------------------------------------------------------
685
686/**
687 * A LINE Distribution, as the model layer and `sn` carry it.
688 *
689 * WHAT EACH CONSUMER READS, which is why all of it is here:
690 * sn.rates, sn.scv the first two moments -- every AMVA path
691 * sn.procid the type tag -- the qsys and QNA dispatch
692 * sn.proc, sn.pie the (D0,D1) pair -- QNA, RQNA, cache, polling
693 * sn.mu, sn.phi the phase rates and completion probabilities
694 * sn.phases the order of that representation
695 * sn.lst the Laplace-Stieltjes transform -- M/G/1 analyzers
696 *
697 * `params` holds the constructor arguments in MATLAB's getParam order, so a
698 * dump can be compared parameter by parameter rather than through the moments,
699 * which two different distributions can share.
700 *
701 * Two values are special and must not be confused, because they enter the
702 * struct as opposite extremes:
703 * Immediate mean = 1/GlobalConstants.Immediate = 1e-8, rate = 1e8
704 * Disabled rate = NaN, which marks a (station, class) pair the class never
705 * visits; the chain and visit machinery keys on it.
706 *
707 * ARITHMETIC. The Markovian families are rational in their parameters and are
708 * built exactly. Gamma, Weibull and Lognormal are not -- their moments call
709 * tgamma or exp -- so their factories refuse by name under exact arithmetic
710 * rather than returning a rounded rational that would look exact.
711 */
712template <class T>
713struct PriorSpec;
714
715template <class T>
716struct Distrib {
717 ProcessType type = ProcessType::DISABLED;
718 /**
719 * The law as DECLARED, when this one is a surrogate fitted over it.
720 *
721 * `sn_nonmarkov_toph` installs a fitted (D0,D1) over a Gamma or a Lognormal
722 * and retags `type` PH or ME, after which nothing names or evaluates the law
723 * the user wrote. MMAP[K]/G[K]/1 reads that law's TRANSFORM rather than the
724 * fit, so it needs the original; every other consumer wants the surrogate
725 * and reads this struct as before. Null when no substitution has happened.
726 * A POINTER, and shared, for the reason `prior` below is one: an inline
727 * member would make the type self-embedding.
728 */
729 std::shared_ptr<Distrib<T>> declared;
730 T mean = num_traits<T>::from_int(0);
731 T scv = num_traits<T>::from_int(1);
732 bool disabled = true;
733 /** Constructor arguments, in MATLAB getParam order. */
734 std::vector<T> params;
735 /** Replayer / Trace samples; empty for every other type. */
736 std::vector<T> trace;
737 /**
738 * The trace FILE a Replayer was read from, when there was one.
739 *
740 * The samples above are what every solver in this port uses, so the path is
741 * carried only for the exporters: `saveServiceStrategy` hands JMT a
742 * `ReplayerPar` naming a file, and a Replayer exported without it is a JMT
743 * model that reads nothing. Empty when the samples were supplied directly.
744 * `network_writer.h` emits it for the same reason: the samples have no wire
745 * form, so without the path a Replayer is written back as the moments and
746 * reloads as a different law.
747 */
748 std::string trace_file;
749 /**
750 * The (D0,D1) pair when the type carries one directly.
751 *
752 * EMPTY for Det, Uniform, Pareto, Gamma, Weibull, Lognormal and Replayer:
753 * MATLAB's getProcess returns their raw PARAMETERS there, and
754 * refreshProcessRepresentations replaces them with an Erlang approximation
755 * (`convertToMAP`) on the way into sn.proc. That conversion is a property
756 * of the refresh, not of the distribution, so it is not done here; see
757 * dist_to_map() in lang/distribution.h.
758 */
759 Matrix<T> D0, D1;
760 /** MMAP per-class D1 blocks / BMAP per-batch-size blocks; empty otherwise. */
761 std::vector<Matrix<T>> Dmark;
762 /**
763 * The alternatives of a `Prior`, set only when `type == PRIOR`.
764 *
765 * A POINTER, and shared: `PriorSpec` holds `Distrib<T>` values, so an
766 * inline member would make the type self-embedding, and the spec is
767 * immutable once built, so copying a service table copies a pointer rather
768 * than a design. `mean` and `scv` beside it are the MIXTURE moments, as
769 * MATLAB's `Prior.getMean`/`getSCV` return: a Prior that reaches
770 * `refresh_rates` therefore lowers to a rate rather than to a NaN. That is
771 * for honesty of the struct dump only -- the featset gate refuses the model
772 * before any solver reads those rates, and `dist_to_map`, `dist_lst` and
773 * `dist_moment` refuse a Prior by name.
774 */
775 std::shared_ptr<PriorSpec<T>> prior;
776 bool is_prior() const { return type == ProcessType::PRIOR; }
777
778 /**
779 * `sn.proc{i}{r} = {breakpoints, A, B, cyclic}` of a MAPt / PHt / NHPP.
780 *
781 * `sched_bp` has one more entry than there are segments -- it is the
782 * BOUNDARY vector, so segment k is in force on [bp(k), bp(k+1)) -- and
783 * `sched_D0[k]`, `sched_D1[k]` are the pair of segment k, ALREADY in MAP
784 * form. A PHt is stored converted, D0 = S and D1 = (-S e) alpha, because
785 * `sn_schedule_nominal` converts it on every read and keeping the raw
786 * (alpha, S) here would make every consumer repeat that conversion and one
787 * of them eventually forget. The raw form is not needed: the conversion is
788 * lossless and nothing downstream asks for alpha again.
789 *
790 * EMPTY FOR EVERY OTHER TYPE. `has_schedule()` is the test, and a solver
791 * with no notion of time simply never calls it -- the nominal pair in
792 * `D0`/`D1` is a complete, time-averaged answer for such a solver.
793 */
794 std::vector<T> sched_bp;
795 std::vector<Matrix<T>> sched_D0, sched_D1;
796 bool sched_cyclic = false;
797 bool has_schedule() const { return !sched_D0.empty(); }
798
799 static Distrib exp_mean(const T& m) {
800 const T one = num_traits<T>::from_int(1);
801 Distrib d;
802 d.type = ProcessType::EXP;
803 d.mean = m;
804 d.scv = one;
805 d.disabled = false;
806 const T lambda = m > num_traits<T>::from_int(0) ? T(one / m)
807 : num_traits<T>::from_double(
808 GlobalConstants::Immediate);
809 d.params.push_back(lambda);
810 d.D0 = Matrix<T>(1, 1, T(-lambda));
811 d.D1 = Matrix<T>(1, 1, lambda);
812 return d;
813 }
814 static Distrib exp_rate(const T& r) {
815 const T zero = num_traits<T>::from_int(0);
816 if (r <= zero) {
817 // Exp.fitRate(0) rationale: see _kb/04-networkstruct.md (cpp port notes)
818 return exp_mean(num_traits<T>::from_double(1.0 / GlobalConstants::Zero));
819 }
820 return exp_mean(T(num_traits<T>::from_int(1) / r));
821 }
822 /**
823 * The Immediate singleton. Its MEAN is zero and its RATE is 1e8, and the
824 * two are deliberately not reciprocal: MATLAB's Immediate.getMean() returns
825 * 0 while Immediate.getRate() returns GlobalConstants.Immediate, and both
826 * are read, by different callers. SolverLN reads the mean (a task with an
827 * Immediate think time contributes no think time); Network.refreshRates
828 * reads the rate (the station serves the class in 1e-8 time units, not in
829 * zero, which would be an infinite service rate the MVA recursion cannot
830 * carry). Collapsing them onto 1/mean or 1/rate breaks one caller or the
831 * other, so the type tag decides.
832 */
833 /**
834 * The point mass at zero.
835 *
836 * Its SCV is 1, NOT the 0 of a degenerate distribution. The reference makes
837 * this explicit -- `Immediate.getSCV` returns 1 in MATLAB, in the JAR and in
838 * Python -- because Immediate is realised downstream as an exponential of
839 * rate GlobalConstants.Immediate rather than as a Dirac: `rate()` below
840 * returns 1e8, and the SCV has to be the one that goes with it. Declaring 0
841 * here is invisible on a layer of PS or infinite-server stations, where the
842 * AMVA correction does not read the SCV at all, and shows up only once a
843 * layer holds an FCFS or multiserver station -- so it survives a model like
844 * lqn_ofbiz and breaks a model like lqn_basic.
845 */
846 static Distrib immediate() {
847 Distrib d;
848 d.type = ProcessType::IMMEDIATE;
849 d.mean = num_traits<T>::from_int(0);
850 d.scv = num_traits<T>::from_int(1);
851 d.disabled = false;
852 const T imm = num_traits<T>::from_double(GlobalConstants::Immediate);
853 d.D0 = Matrix<T>(1, 1, T(-imm));
854 d.D1 = Matrix<T>(1, 1, imm);
855 return d;
856 }
857 static Distrib disabled_dist() { return Distrib(); }
858 static Distrib det(const T& m) {
859 Distrib d;
860 d.type = ProcessType::DET;
861 d.mean = m;
862 d.scv = num_traits<T>::from_int(0);
863 d.disabled = false;
864 d.params.push_back(m);
865 return d;
866 }
867
868 // -----------------------------------------------------------------------
869 // The Markovian families: (D0,D1) is built here, exactly
870 // -----------------------------------------------------------------------
871
872 /** Erlang(alpha, r): r phases of rate alpha, as MATLAB's Erlang(phaseRate, nphases). */
873 static Distrib erlang(const T& phase_rate, std::size_t r) {
874 const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
875 if (r == 0) throw InputError("Erlang: the number of phases must be positive");
876 if (!(phase_rate > zero)) throw InputError("Erlang: the phase rate must be positive");
877 Distrib d;
878 d.type = ProcessType::ERLANG;
879 d.disabled = false;
880 d.params.push_back(phase_rate);
881 d.params.push_back(num_traits<T>::from_int(static_cast<long>(r)));
882 d.mean = T(num_traits<T>::from_int(static_cast<long>(r)) / phase_rate);
883 d.scv = T(one / num_traits<T>::from_int(static_cast<long>(r)));
884 d.D0 = Matrix<T>(r, r, zero);
885 d.D1 = Matrix<T>(r, r, zero);
886 for (std::size_t i = 0; i < r; ++i) {
887 d.D0(i, i) = T(-phase_rate);
888 if (i + 1 < r) d.D0(i, i + 1) = phase_rate;
889 }
890 d.D1(r - 1, 0) = phase_rate;
891 return d;
892 }
893
894 /**
895 * Erlang fitted to a mean and an SCV, as MATLAB's Erlang.fitMeanAndSCV.
896 *
897 * AN SCV ABOVE ONE IS REFUSED, not answered. An Erlang of order r has
898 * SCV = 1/r, so the family reaches 1 and no higher; ceil(1/c2) is 1 for
899 * every c2 > 1, and returning that means handing back an EXPONENTIAL under
900 * the name of the distribution the caller asked for. MATLAB errors here and
901 * this port now does too, so a mis-specified SCV is a diagnostic rather
902 * than a silently different service process.
903 */
904 static Distrib erlang_fit(const T& m, const T& c2) {
905 const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
906 if (!(c2 > zero)) throw InputError("Erlang.fitMeanAndSCV: the SCV must be positive");
907 if (c2 > one)
908 throw InputError(
909 "Erlang.fitMeanAndSCV: the Erlang distribution requires a squared coefficient "
910 "of variation <= 1; use HyperExp, Coxian or APH above 1");
911 // MATLAB: r = ceil(1/scv); alpha = r/mean
912 const double r_d = std::ceil(1.0 / num_traits<T>::to_double(c2));
913 const std::size_t r = static_cast<std::size_t>(r_d < 1.0 ? 1.0 : r_d);
914 return erlang(T(num_traits<T>::from_int(static_cast<long>(r)) / m), r);
915 }
916
917 /**
918 * HyperExp(p, lambda1, lambda2): phase i chosen with probability p_i.
919 *
920 * D1(i,j) = mu(i) p(j) -- the OUTER product. Associating the other way
921 * gives D1(i,j) = mu(i) p(i) replicated across the row, whose rows no
922 * longer sum with D0 to zero; _kb records that trap in the MATLAB class.
923 */
924 /**
925 * The same for any number of branches, which is what MATLAB `HyperExp`
926 * accepts and what the writers emit as vector `p` and `lambda`. The
927 * two-branch entry point stays because it carries MATLAB's getParam order
928 * (p, lambda1, lambda2), which a parameter-by-parameter dump compares
929 * against.
930 */
931 static Distrib hyperexp_n(const std::vector<T>& p, const std::vector<T>& lambda) {
932 const T zero = num_traits<T>::from_int(0), two = num_traits<T>::from_int(2);
933 const std::size_t n = p.size();
934 if (n == 0 || lambda.size() != n)
935 throw InputError("HyperExp: p and lambda must be non-empty and of equal length");
936 Distrib d;
937 d.type = ProcessType::HYPEREXP;
938 d.disabled = false;
939 for (const T& v : p) d.params.push_back(v);
940 for (const T& v : lambda) d.params.push_back(v);
941 d.D0 = Matrix<T>(n, n, zero);
942 d.D1 = Matrix<T>(n, n, zero);
943 T m1 = zero, m2 = zero;
944 for (std::size_t i = 0; i < n; ++i) {
945 if (!(lambda[i] > zero)) throw InputError("HyperExp: the phase rates must be positive");
946 d.D0(i, i) = T(-lambda[i]);
947 for (std::size_t j = 0; j < n; ++j) d.D1(i, j) = T(lambda[i] * p[j]);
948 m1 += T(p[i] / lambda[i]);
949 m2 += T(two * p[i] / (lambda[i] * lambda[i]));
950 }
951 d.mean = m1;
952 d.scv = T((m2 - m1 * m1) / (m1 * m1));
953 return d;
954 }
955
956 static Distrib hyperexp(const T& p, const T& lambda1, const T& lambda2) {
957 const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
958 if (!(lambda1 > zero) || !(lambda2 > zero))
959 throw InputError("HyperExp: the phase rates must be positive");
960 if (p < zero || p > one) throw InputError("HyperExp: p is not a probability");
961 Distrib d;
962 d.type = ProcessType::HYPEREXP;
963 d.disabled = false;
964 d.params.push_back(p);
965 d.params.push_back(lambda1);
966 d.params.push_back(lambda2);
967 const T q = T(one - p);
968 d.mean = T(p / lambda1 + q / lambda2);
969 const T m2 = T(num_traits<T>::from_int(2) *
970 (p / (lambda1 * lambda1) + q / (lambda2 * lambda2)));
971 d.scv = T((m2 - d.mean * d.mean) / (d.mean * d.mean));
972 d.D0 = Matrix<T>(2, 2, zero);
973 d.D1 = Matrix<T>(2, 2, zero);
974 d.D0(0, 0) = T(-lambda1);
975 d.D0(1, 1) = T(-lambda2);
976 d.D1(0, 0) = T(lambda1 * p);
977 d.D1(0, 1) = T(lambda1 * q);
978 d.D1(1, 0) = T(lambda2 * p);
979 d.D1(1, 1) = T(lambda2 * q);
980 return d;
981 }
982
983 /**
984 * Coxian(mu, phi): phase i completes with probability phi(i) and otherwise
985 * moves to phase i+1. The last phase always completes.
986 */
987 static Distrib coxian(const std::vector<T>& mu, const std::vector<T>& phi) {
988 const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
989 const std::size_t n = mu.size();
990 if (n == 0 || phi.size() != n)
991 throw InputError("Coxian: mu and phi must be non-empty and of equal length");
992 Distrib d;
993 d.type = n == 2 ? ProcessType::COX2 : ProcessType::COXIAN;
994 d.disabled = false;
995 for (const T& v : mu) d.params.push_back(v);
996 for (const T& v : phi) d.params.push_back(v);
997 d.D0 = Matrix<T>(n, n, zero);
998 d.D1 = Matrix<T>(n, n, zero);
999 for (std::size_t i = 0; i < n; ++i) {
1000 if (!(mu[i] > zero)) throw InputError("Coxian: the phase rates must be positive");
1001 d.D0(i, i) = T(-mu[i]);
1002 if (i + 1 < n) d.D0(i, i + 1) = T(mu[i] * (one - phi[i]));
1003 d.D1(i, 0) = T(mu[i] * phi[i]);
1004 }
1005 // moments of the absorbing chain started in phase 1: m_k = k! e_1 (-D0)^-k e
1006 d.mean = ph_moment_from(d.D0, 1, 1);
1007 const T m2 = ph_moment_from(d.D0, 1, 2);
1008 d.scv = T((m2 - d.mean * d.mean) / (d.mean * d.mean));
1009 return d;
1010 }
1011
1012 /** Cox2(mu1, mu2, phi1), MATLAB's two-phase Coxian constructor. */
1013 static Distrib cox2(const T& mu1, const T& mu2, const T& phi1) {
1014 std::vector<T> mu, phi;
1015 mu.push_back(mu1);
1016 mu.push_back(mu2);
1017 phi.push_back(phi1);
1018 phi.push_back(num_traits<T>::from_int(1));
1019 return coxian(mu, phi);
1020 }
1021
1022 /**
1023 * PH / APH given by (alpha, A): D0 = A and D1 = (-A e) alpha.
1024 *
1025 * `acyclic` selects the type tag only; the representation is the same, and
1026 * no consumer of sn.proc distinguishes them.
1027 */
1028 static Distrib phase_type(const std::vector<T>& alpha, const Matrix<T>& A, bool acyclic) {
1029 const T zero = num_traits<T>::from_int(0);
1030 const std::size_t n = alpha.size();
1031 if (n == 0 || A.rows() != n || A.cols() != n)
1032 throw InputError("PH: alpha and the subgenerator have inconsistent sizes");
1033 Distrib d;
1034 d.type = acyclic ? ProcessType::APH : ProcessType::PH;
1035 d.disabled = false;
1036 for (const T& v : alpha) d.params.push_back(v);
1037 d.D0 = A;
1038 d.D1 = Matrix<T>(n, n, zero);
1039 for (std::size_t i = 0; i < n; ++i) {
1040 T out = zero;
1041 for (std::size_t j = 0; j < n; ++j) out += A(i, j);
1042 for (std::size_t j = 0; j < n; ++j) d.D1(i, j) = T(-out * alpha[j]);
1043 }
1044 d.mean = ph_moment(alpha, A, 1);
1045 const T m2 = ph_moment(alpha, A, 2);
1046 d.scv = T((m2 - d.mean * d.mean) / (d.mean * d.mean));
1047 return d;
1048 }
1049
1050 /** A MAP given by its two matrices; the moments are those of its stationary phase. */
1051 static Distrib map_dist(const Matrix<T>& D0, const Matrix<T>& D1, ProcessType tag) {
1052 if (D0.rows() != D0.cols() || D1.rows() != D1.cols() || D0.rows() != D1.rows())
1053 throw InputError("MAP: D0 and D1 must be square and of the same order");
1054 Distrib d;
1055 d.type = tag;
1056 d.disabled = false;
1057 d.D0 = D0;
1058 d.D1 = D1;
1059 // MAP moment rationale: see _kb/04-networkstruct.md (cpp port notes)
1060 d.mean = num_traits<T>::from_int(0);
1061 d.scv = num_traits<T>::from_int(1);
1062 return d;
1063 }
1064
1065 // -----------------------------------------------------------------------
1066 // The time-inhomogeneous families (Ko-Pender)
1067 // -----------------------------------------------------------------------
1068
1069 /** The shared constructor of the three schedule families. */
1070 static Distrib sched_dist(const std::vector<T>& breakpoints,
1071 const std::vector<Matrix<T>>& D0segs,
1072 const std::vector<Matrix<T>>& D1segs, bool cyclic, ProcessType tag) {
1073 const T zero = num_traits<T>::from_int(0);
1074 const std::string who = process_to_text(tag);
1075 const std::size_t n = D0segs.size();
1076 if (n == 0 || D1segs.size() != n)
1077 throw InputError(who + ": the two segment lists must be non-empty and equally long");
1078 if (breakpoints.size() != n + 1)
1079 throw InputError(who +
1080 ": breakpoints is the BOUNDARY vector and must hold one more entry "
1081 "than there are segments");
1082 for (std::size_t k = 0; k + 1 < breakpoints.size(); ++k)
1083 if (!(breakpoints[k + 1] > breakpoints[k]))
1084 throw InputError(who + ": breakpoints must be strictly increasing");
1085 const std::size_t order = D0segs[0].rows();
1086 for (std::size_t k = 0; k < n; ++k)
1087 if (D0segs[k].rows() != order || D0segs[k].cols() != order ||
1088 D1segs[k].rows() != order || D1segs[k].cols() != order)
1089 throw InputError(who +
1090 ": every segment must have the same order; the schedule "
1091 "modulates one phase structure and does not switch between them");
1092
1093 Distrib d;
1094 d.type = tag;
1095 d.disabled = false;
1096 d.sched_bp = breakpoints;
1097 d.sched_D0 = D0segs;
1098 d.sched_D1 = D1segs;
1099 d.sched_cyclic = cyclic;
1100 // The nominal pair: the width-weighted time average, i.e. the first two
1101 // outputs of `sn_schedule_nominal`.
1102 T total = zero;
1103 for (std::size_t k = 0; k < n; ++k) total += T(breakpoints[k + 1] - breakpoints[k]);
1104 d.D0 = Matrix<T>(order, order, zero);
1105 d.D1 = Matrix<T>(order, order, zero);
1106 for (std::size_t k = 0; k < n; ++k) {
1107 const T w = T(T(breakpoints[k + 1] - breakpoints[k]) / total);
1108 for (std::size_t a = 0; a < order; ++a)
1109 for (std::size_t b = 0; b < order; ++b) {
1110 d.D0(a, b) += w * D0segs[k](a, b);
1111 d.D1(a, b) += w * D1segs[k](a, b);
1112 }
1113 }
1114 // Same convention as `map_dist`: the moments of a MAP are computed by
1115 // `dist_moment`, not stored.
1116 d.mean = zero;
1117 d.scv = num_traits<T>::from_int(1);
1118 return d;
1119 }
1120
1121 /**
1122 * MAPt(breakpoints, {D0_k}, {D1_k}, cyclic): a piecewise-constant
1123 * (D0(t), D1(t)).
1124 *
1125 * `breakpoints` is the boundary vector, so it holds one more entry than
1126 * there are segments and must be strictly increasing. Every segment must
1127 * have the SAME ORDER: the schedule modulates one phase structure, it does
1128 * not switch between structures, and a solver that integrated across a
1129 * change of order would have no way to map the phase occupancy across the
1130 * boundary. That is the reference's own constructor requirement.
1131 *
1132 * `D0` / `D1` are set to the WIDTH-WEIGHTED TIME AVERAGE of the segments,
1133 * which is `sn_schedule_nominal`'s nominal pair: it is the stationary
1134 * carrier of the phase structure the schedule modulates, so the phase count
1135 * and the mean rate a time-blind consumer reads are the ones the model
1136 * actually has.
1137 */
1138 static Distrib mapt(const std::vector<T>& breakpoints, const std::vector<Matrix<T>>& D0segs,
1139 const std::vector<Matrix<T>>& D1segs, bool cyclic) {
1140 return sched_dist(breakpoints, D0segs, D1segs, cyclic, ProcessType::MAPT);
1141 }
1142
1143 /**
1144 * PHt(breakpoints, {alpha_k}, {S_k}, cyclic), stored as its equivalent MAP
1145 * schedule: D0 = S and D1 = (-S e) alpha, the pair `sn_schedule_nominal`
1146 * builds from a PHt slot.
1147 */
1148 static Distrib pht(const std::vector<T>& breakpoints, const std::vector<std::vector<T>>& alphas,
1149 const std::vector<Matrix<T>>& Ssegs, bool cyclic) {
1150 const T zero = num_traits<T>::from_int(0);
1151 if (alphas.size() != Ssegs.size() || alphas.empty())
1152 throw InputError("PHt: alpha and S must have the same, non-zero number of segments");
1153 std::vector<Matrix<T>> D0segs, D1segs;
1154 for (std::size_t k = 0; k < Ssegs.size(); ++k) {
1155 const Matrix<T>& S = Ssegs[k];
1156 const std::vector<T>& a = alphas[k];
1157 if (S.rows() != S.cols() || S.rows() != a.size())
1158 throw InputError("PHt: alpha and the sub-generator disagree in order");
1159 Matrix<T> D1(S.rows(), S.cols(), zero);
1160 for (std::size_t i = 0; i < S.rows(); ++i) {
1161 T s = zero;
1162 for (std::size_t j = 0; j < S.cols(); ++j) s += S(i, j);
1163 for (std::size_t j = 0; j < S.cols(); ++j) D1(i, j) = T(-s * a[j]);
1164 }
1165 D0segs.push_back(S);
1166 D1segs.push_back(D1);
1167 }
1168 Distrib d = sched_dist(breakpoints, D0segs, D1segs, cyclic, ProcessType::PHT);
1169 return d;
1170 }
1171
1172 /**
1173 * NHPP(breakpoints, rates, cyclic): a MAPt of ORDER ONE, which is what an
1174 * inhomogeneous Poisson process is. Building it through the same path is
1175 * what makes every schedule consumer see one representation.
1176 */
1177 static Distrib nhpp(const std::vector<T>& breakpoints, const std::vector<T>& rates,
1178 bool cyclic) {
1179 std::vector<Matrix<T>> D0segs, D1segs;
1180 for (std::size_t k = 0; k < rates.size(); ++k) {
1181 Matrix<T> a(1, 1, T(-rates[k])), b(1, 1, rates[k]);
1182 D0segs.push_back(a);
1183 D1segs.push_back(b);
1184 }
1185 return sched_dist(breakpoints, D0segs, D1segs, cyclic, ProcessType::NHPP);
1186 }
1187
1188 // -----------------------------------------------------------------------
1189 // The non-Markovian families: parameters only, as MATLAB's getProcess
1190 // -----------------------------------------------------------------------
1191
1192 /** Uniform(a, b). */
1193 static Distrib uniform(const T& a, const T& b) {
1194 if (!(b > a)) throw InputError("Uniform: the upper bound must exceed the lower bound");
1195 Distrib d;
1196 d.type = ProcessType::UNIFORM;
1197 d.disabled = false;
1198 d.params.push_back(a);
1199 d.params.push_back(b);
1200 d.mean = T((a + b) / num_traits<T>::from_int(2));
1201 const T w = T(b - a);
1202 d.scv = T((w * w / num_traits<T>::from_int(12)) / (d.mean * d.mean));
1203 return d;
1204 }
1205
1206 /** Pareto(shape, scale), with the MATLAB parameter order (alpha, k). */
1207 static Distrib pareto(const T& shape, const T& scale) {
1208 const T one = num_traits<T>::from_int(1), two = num_traits<T>::from_int(2);
1209 if (!(shape > two))
1210 throw InputError("Pareto: the shape must exceed 2 for a finite variance");
1211 Distrib d;
1212 d.type = ProcessType::PARETO;
1213 d.disabled = false;
1214 d.params.push_back(shape);
1215 d.params.push_back(scale);
1216 d.mean = T(shape * scale / (shape - one));
1217 const T var = T(scale * scale * shape / ((shape - one) * (shape - one)) / (shape - two));
1218 d.scv = T(var / (d.mean * d.mean));
1219 return d;
1220 }
1221
1222 /**
1223 * Gamma(shape, scale), Weibull(scale, shape) and Lognormal(mu, sigma).
1224 *
1225 * Their moments are values of the gamma function or of exp, so they exist
1226 * only where the arithmetic has transcendentals. Under exact arithmetic the
1227 * factory REFUSES rather than storing a rounded rational: a rational that
1228 * came out of tgamma is not the exact moment of the distribution, and every
1229 * downstream claim of exactness would be false.
1230 */
1231 static Distrib gamma_dist(const T& shape, const T& scale) {
1232 const T one = num_traits<T>::from_int(1);
1233 Distrib d;
1234 d.type = ProcessType::GAMMA;
1235 d.disabled = false;
1236 d.params.push_back(shape);
1237 d.params.push_back(scale);
1238 d.mean = T(shape * scale);
1239 d.scv = T(one / shape);
1240 return d;
1241 }
1242
1243 static Distrib weibull(const T& scale, const T& shape) {
1244 if constexpr (!num_traits<T>::has_transcendental) {
1245 throw UnsupportedError(
1246 "Weibull: its moments are values of the gamma function, which exact arithmetic "
1247 "has no representation for; use the double or real backend");
1248 } else {
1249 const double a = num_traits<T>::to_double(scale);
1250 const double r = num_traits<T>::to_double(shape);
1251 if (!(a > 0.0) || !(r > 0.0))
1252 throw InputError("Weibull: the scale and the shape must be positive");
1253 const double g1 = std::tgamma(1.0 + 1.0 / r);
1254 const double g2 = std::tgamma(1.0 + 2.0 / r);
1255 Distrib d;
1256 d.type = ProcessType::WEIBULL;
1257 d.disabled = false;
1258 d.params.push_back(scale);
1259 d.params.push_back(shape);
1260 d.mean = num_traits<T>::from_double(a * g1);
1261 d.scv = num_traits<T>::from_double((g2 - g1 * g1) / (g1 * g1));
1262 return d;
1263 }
1264 }
1265
1266 static Distrib lognormal(const T& logmean, const T& logsigma) {
1267 if constexpr (!num_traits<T>::has_transcendental) {
1268 throw UnsupportedError(
1269 "Lognormal: its moments are values of exp, which exact arithmetic has no "
1270 "representation for; use the double or real backend");
1271 } else {
1272 const double mu = num_traits<T>::to_double(logmean);
1273 const double sg = num_traits<T>::to_double(logsigma);
1274 if (!(sg > 0.0)) throw InputError("Lognormal: sigma must be positive");
1275 Distrib d;
1276 d.type = ProcessType::LOGNORMAL;
1277 d.disabled = false;
1278 d.params.push_back(logmean);
1279 d.params.push_back(logsigma);
1280 d.mean = num_traits<T>::from_double(std::exp(mu + sg * sg / 2.0));
1281 d.scv = num_traits<T>::from_double(std::exp(sg * sg) - 1.0);
1282 return d;
1283 }
1284 }
1285
1286 /**
1287 * `Normal(mu, sigma)`: the Gaussian, for use as a continuous `Prior`'s
1288 * parameter density.
1289 *
1290 * `scv` is Inf at mu = 0, which is `Normal.m:52-63`'s own answer and not a
1291 * degradation: the SCV of a zero-mean law is not defined, and the reference
1292 * says so rather than dividing.
1293 */
1294 static Distrib normal(const T& mu, const T& sigma) {
1295 if (!(num_traits<T>::to_double(sigma) > 0.0))
1296 throw InputError("Normal: sigma must be positive");
1297 Distrib d;
1298 d.type = ProcessType::NORMAL;
1299 d.disabled = false;
1300 d.params.push_back(mu);
1301 d.params.push_back(sigma);
1302 d.mean = mu;
1303 const double m = num_traits<T>::to_double(mu), s = num_traits<T>::to_double(sigma);
1304 d.scv = std::abs(m) < GlobalConstants::FineTol
1305 ? num_traits<T>::from_double(std::numeric_limits<double>::infinity())
1306 : num_traits<T>::from_double((s * s) / (m * m));
1307 return d;
1308 }
1309
1310 /**
1311 * Replayer / Trace read FROM A FILE, which keeps the path beside the samples.
1312 *
1313 * Every solver in this port replays `trace`; the path is what the EXPORTERS
1314 * need. JMT is handed a `ReplayerPar` naming a file and has no way to take
1315 * samples inline, so a Replayer exported without it is a JMT model that
1316 * reads nothing -- `oqn_trace_driven` was refused outright for exactly that.
1317 */
1318 static Distrib replayer_from(const std::vector<T>& samples, const std::string& path) {
1319 Distrib d = replayer(samples);
1320 d.trace_file = path;
1321 return d;
1322 }
1323
1324 /** Replayer / Trace: the samples, with their empirical first two moments. */
1325 static Distrib replayer(const std::vector<T>& samples) {
1326 if (samples.empty()) throw InputError("Replayer: the trace is empty");
1327 Distrib d;
1328 d.type = ProcessType::REPLAYER;
1329 d.disabled = false;
1330 d.trace = samples;
1331 const T n = num_traits<T>::from_int(static_cast<long>(samples.size()));
1332 T s1 = num_traits<T>::from_int(0), s2 = num_traits<T>::from_int(0);
1333 for (const T& x : samples) {
1334 s1 += x;
1335 s2 += T(x * x);
1336 }
1337 d.mean = T(s1 / n);
1338 d.scv = T((s2 / n - d.mean * d.mean) / (d.mean * d.mean));
1339 return d;
1340 }
1341
1342 // -----------------------------------------------------------------------
1343 // The discrete families
1344 //
1345 // MATLAB's getProcess for each of these returns [mean, scv] and nothing
1346 // else, so they carry NO (D0,D1): refreshProcessRepresentations replaces
1347 // them with the Erlang fit convertToMAP, which `dist_to_map` reproduces
1348 // from the two moments alone. Storing an invented representation here
1349 // would make sn.proc disagree with the reference for the same model.
1350 // -----------------------------------------------------------------------
1351
1352 /** DiscreteUniform(a, b) over the integers a..b inclusive. */
1353 static Distrib discrete_uniform(const T& a, const T& b) {
1354 const T two = num_traits<T>::from_int(2), one = num_traits<T>::from_int(1);
1355 if (!(b >= a)) throw InputError("DiscreteUniform: the upper bound must not be below the lower");
1356 Distrib d;
1357 d.type = ProcessType::DUNIFORM;
1358 d.disabled = false;
1359 d.params.push_back(a);
1360 d.params.push_back(b);
1361 d.mean = T((a + b) / two);
1362 const T w = T(b - a + one);
1363 const T var = T((w * w - one) / num_traits<T>::from_int(12));
1364 d.scv = T(var / (d.mean * d.mean));
1365 return d;
1366 }
1367
1368 /** Bernoulli(p): one trial, mean p and variance p(1-p). */
1369 static Distrib bernoulli(const T& p) {
1370 const T one = num_traits<T>::from_int(1);
1371 Distrib d;
1372 d.type = ProcessType::BERNOULLI;
1373 d.disabled = false;
1374 d.params.push_back(p);
1375 d.mean = p;
1376 d.scv = T((one - p) / p);
1377 return d;
1378 }
1379
1380 /** Binomial(n, p). */
1381 static Distrib binomial(const T& n, const T& p) {
1382 const T one = num_traits<T>::from_int(1);
1383 Distrib d;
1384 d.type = ProcessType::BINOMIAL;
1385 d.disabled = false;
1386 d.params.push_back(n);
1387 d.params.push_back(p);
1388 d.mean = T(n * p);
1389 d.scv = T((one - p) / (n * p));
1390 return d;
1391 }
1392
1393 /**
1394 * Poisson(lambda), whose SCV is 1/lambda -- the count's variance is lambda
1395 * and its mean is lambda, so this is NOT the exponential's SCV of 1.
1396 */
1397 static Distrib poisson(const T& lambda) {
1398 const T one = num_traits<T>::from_int(1);
1399 Distrib d;
1400 d.type = ProcessType::POISSON;
1401 d.disabled = false;
1402 d.params.push_back(lambda);
1403 d.mean = lambda;
1404 d.scv = T(one / lambda);
1405 return d;
1406 }
1407
1408 /**
1409 * Geometric(p) on the MATLAB convention: the NUMBER OF TRIALS to the first
1410 * success, support {1, 2, ...}, so the mean is 1/p and the SCV is 1-p.
1411 */
1412 static Distrib geometric(const T& p) {
1413 const T one = num_traits<T>::from_int(1);
1414 Distrib d;
1415 d.type = ProcessType::GEOMETRIC;
1416 d.disabled = false;
1417 d.params.push_back(p);
1418 d.mean = T(one / p);
1419 d.scv = T(one - p);
1420 return d;
1421 }
1422
1423 /**
1424 * Zipf(s, n) over the ranks 1..n, with the generalized harmonic moments
1425 * H(s-1,n)/H(s,n) and H(s-2,n)/H(s,n) MATLAB `Zipf.m` uses.
1426 *
1427 * The harmonic sums call pow for a non-integer shape, so the factory is
1428 * gated on transcendental arithmetic exactly as Weibull and Lognormal are.
1429 */
1430 static Distrib zipf(const T& s, std::size_t n) {
1431 if constexpr (!num_traits<T>::has_transcendental) {
1432 (void)s;
1433 (void)n;
1434 throw UnsupportedError(
1435 "Zipf: its moments are generalized harmonic sums of a real exponent, which exact "
1436 "arithmetic has no representation for; use the double or real backend");
1437 } else {
1438 if (n == 0) throw InputError("Zipf: the item count must be positive");
1439 const double sv = num_traits<T>::to_double(s);
1440 auto harmonic = [n](double e) {
1441 double acc = 0.0;
1442 for (std::size_t k = 1; k <= n; ++k) acc += std::pow(double(k), -e);
1443 return acc;
1444 };
1445 const double h0 = harmonic(sv), h1 = harmonic(sv - 1.0), h2 = harmonic(sv - 2.0);
1446 Distrib d;
1447 d.type = ProcessType::ZIPF;
1448 d.disabled = false;
1449 d.params.push_back(s);
1450 d.params.push_back(num_traits<T>::from_int(static_cast<long>(n)));
1451 const double m1 = h1 / h0;
1452 d.mean = num_traits<T>::from_double(m1);
1453 d.scv = num_traits<T>::from_double((h2 / h0 - m1 * m1) / (m1 * m1));
1454 return d;
1455 }
1456 }
1457
1458 /**
1459 * DiscreteSampler(p, x): the pmf p over the points x.
1460 *
1461 * THE MOMENTS ARE TAKEN OVER x. This port, and MATLAB
1462 * `DiscreteSampler.getMean` with it, used to weight by the RANKS 1..n
1463 * instead; the two agree on the default x = 1:n, which is the form the
1464 * cache popularity vectors are written in, so the rank form survived
1465 * unnoticed until a fork's jobs-per-link distribution arrived on a shifted
1466 * support. The JAR and native Python already weighted by x.
1467 */
1468 static Distrib discrete_sampler(const std::vector<T>& p, const std::vector<T>& x) {
1469 if (p.empty()) throw InputError("DiscreteSampler: the probability vector is empty");
1470 if (!x.empty() && x.size() != p.size())
1471 throw InputError("DiscreteSampler: p and x must have the same length");
1472 Distrib d;
1473 d.type = ProcessType::DISCRETESAMPLER;
1474 d.disabled = false;
1475 d.params = p;
1476 d.trace = x;
1477 T m1 = num_traits<T>::from_int(0), m2 = num_traits<T>::from_int(0);
1478 T tot = num_traits<T>::from_int(0);
1479 for (std::size_t k = 0; k < p.size(); ++k) {
1480 // an absent x is the default support 1..n
1481 const T pt = x.empty() ? num_traits<T>::from_int(static_cast<long>(k + 1)) : x[k];
1482 m1 += T(p[k] * pt);
1483 m2 += T(p[k] * pt * pt);
1484 tot += p[k];
1485 }
1486 m1 = T(m1 / tot);
1487 m2 = T(m2 / tot);
1488 d.mean = m1;
1489 d.scv = T((m2 - m1 * m1) / (m1 * m1));
1490 return d;
1491 }
1492
1493 /**
1494 * EmpiricalCDF(x, F): the moments of the MIDPOINT rule over the CDF bins,
1495 * which is what MATLAB `EmpiricalCDF.getMoments` integrates -- each bin
1496 * contributes its midpoint raised to the moment order, weighted by the CDF
1497 * increment. The rows are the (F, x) pairs in the order they arrive.
1498 */
1499 static Distrib empirical_cdf(const std::vector<T>& x, const std::vector<T>& F) {
1500 if (x.size() != F.size() || x.size() < 2)
1501 throw InputError("EmpiricalCDF: x and F must be equally long and hold at least two points");
1502 const T two = num_traits<T>::from_int(2);
1503 Distrib d;
1504 d.type = ProcessType::EMPIRICALCDF;
1505 d.disabled = false;
1506 d.trace = x;
1507 d.params = F;
1508 T m1 = num_traits<T>::from_int(0), m2 = num_traits<T>::from_int(0);
1509 for (std::size_t i = 0; i + 1 < x.size(); ++i) {
1510 const T mid = T((x[i + 1] - x[i]) / two + x[i]);
1511 const T w = T(F[i + 1] - F[i]);
1512 m1 += T(mid * w);
1513 m2 += T(mid * mid * w);
1514 }
1515 d.mean = m1;
1516 d.scv = T(m2 / (m1 * m1) - num_traits<T>::from_int(1));
1517 return d;
1518 }
1519
1520 // -----------------------------------------------------------------------
1521 // The matrix-exponential and discrete-time Markovian families
1522 // -----------------------------------------------------------------------
1523
1524 /**
1525 * ME(alpha, A): the matrix-exponential distribution, whose moments are the
1526 * phase-type ones -- k! alpha (-A)^-k e -- evaluated by DEFINITION rather
1527 * than through the stationary vector of A + (-Ae)alpha. MATLAB `ME.getMean`
1528 * makes the same choice and says why: the stationary solve is a
1529 * probabilistic object that a non-Markovian A degrades badly (a CME of
1530 * order 101 lost 2.6e-4 in the SCV that way, against 1e-13 by definition).
1531 */
1532 static Distrib me(const std::vector<T>& alpha, const Matrix<T>& A) {
1533 Distrib d = phase_type(alpha, A, false);
1534 d.type = ProcessType::ME;
1535 return d;
1536 }
1537
1538 /** RAP(H0, H1): a rational arrival process, whose moments are the MAP ones. */
1539 static Distrib rap(const Matrix<T>& H0, const Matrix<T>& H1) {
1540 return map_dist(H0, H1, ProcessType::RAP);
1541 }
1542
1543 /**
1544 * DMAP(D0, D1): a DISCRETE-time MAP, where D0 + D1 is stochastic rather
1545 * than a generator. Its moments cannot come from `dist_refresh_moments`'s
1546 * continuous formulas; `dmap_moments` in lang/distribution.h fills them.
1547 */
1548 static Distrib dmap(const Matrix<T>& D0, const Matrix<T>& D1) {
1549 return map_dist(D0, D1, ProcessType::DMAP);
1550 }
1551
1552 /**
1553 * MMAP: D0 plus one D1 block per mark. `D1` is their sum, the aggregate
1554 * arrival matrix every unmarked consumer reads, and the blocks stay in
1555 * `Dmark` for the ones that distinguish marks.
1556 */
1557 static Distrib mmap(const Matrix<T>& D0, const std::vector<Matrix<T>>& D1k) {
1558 if (D1k.empty()) throw InputError("MMAP: no marked arrival block was given");
1559 Matrix<T> agg(D0.rows(), D0.cols(), num_traits<T>::from_int(0));
1560 for (const Matrix<T>& Dk : D1k) {
1561 if (Dk.rows() != D0.rows() || Dk.cols() != D0.cols())
1562 throw InputError("MMAP: every marked block must have the order of D0");
1563 for (std::size_t i = 0; i < agg.rows(); ++i)
1564 for (std::size_t j = 0; j < agg.cols(); ++j) agg(i, j) += Dk(i, j);
1565 }
1566 Distrib d = map_dist(D0, agg, ProcessType::MMAP);
1567 d.Dmark = D1k;
1568 return d;
1569 }
1570
1571 /**
1572 * BMAP: the batch-size blocks D0, D1, ..., Dk, where Dj carries an arrival
1573 * of batch size j. The wire form is the whole list including D0, so the
1574 * head is split off here.
1575 */
1576 static Distrib bmap(const std::vector<Matrix<T>>& D) {
1577 if (D.size() < 2) throw InputError("BMAP: the block list must carry D0 and at least one batch block");
1578 const std::vector<Matrix<T>> batches(D.begin() + 1, D.end());
1579 Distrib d = mmap(D[0], batches);
1580 d.type = ProcessType::BMAP;
1581 return d;
1582 }
1583
1584 bool is_immediate() const { return type == ProcessType::IMMEDIATE; }
1585
1586 /** True when the type carries a (D0,D1) pair of its own. */
1587 bool has_map() const { return D0.rows() > 0; }
1588
1589 /**
1590 * Phase rates, MATLAB's getMu: the total outgoing rate of each phase.
1591 *
1592 * Empty when the type carries no representation, which is what
1593 * refreshProcessPhases writes as NaN for a Fork or a Join.
1594 */
1595 std::vector<T> mu_vec() const {
1596 std::vector<T> v;
1597 if (!has_map()) {
1598 if (disabled) return v;
1599 v.push_back(rate()); // one phase at the mean rate, as MATLAB does
1600 return v;
1601 }
1602 for (std::size_t i = 0; i < D0.rows(); ++i) v.push_back(T(-D0(i, i)));
1603 return v;
1604 }
1605
1606 /** Completion probabilities, MATLAB's getPhi: (D1 e) ./ (-diag(D0)). */
1607 std::vector<T> phi_vec() const {
1608 std::vector<T> v;
1609 if (!has_map()) {
1610 if (disabled) return v;
1611 v.push_back(num_traits<T>::from_int(1));
1612 return v;
1613 }
1614 const T zero = num_traits<T>::from_int(0);
1615 for (std::size_t i = 0; i < D1.rows(); ++i) {
1616 T s = zero;
1617 for (std::size_t j = 0; j < D1.cols(); ++j) s += D1(i, j);
1618 const T out = T(-D0(i, i));
1619 v.push_back(out == zero ? num_traits<T>::from_int(1) : T(s / out));
1620 }
1621 return v;
1622 }
1623
1624 /** The order of the representation, MATLAB's sn.phases. */
1625 std::size_t phases() const {
1626 if (disabled) return 0;
1627 return has_map() ? D0.rows() : 1;
1628 }
1629
1630 /**
1631 * The k-th raw moment of a phase-type (alpha, A): k! alpha (-A)^-k e.
1632 *
1633 * The inverse is never formed: the powers are accumulated by repeated
1634 * solves of (-A) x = b, which is exact in rational arithmetic and stable
1635 * in floating point.
1636 */
1637 static T ph_moment(const std::vector<T>& alpha, const Matrix<T>& A, unsigned k) {
1638 const std::size_t n = alpha.size();
1639 std::vector<T> x(n, num_traits<T>::from_int(1));
1640 for (unsigned i = 0; i < k; ++i) x = solve_neg(A, x);
1641 T acc = num_traits<T>::from_int(0);
1642 for (std::size_t i = 0; i < n; ++i) acc += alpha[i] * x[i];
1643 T fact = num_traits<T>::from_int(1);
1644 for (unsigned i = 2; i <= k; ++i) fact *= num_traits<T>::from_int(static_cast<long>(i));
1645 return T(fact * acc);
1646 }
1647
1648 /** The same, for a representation entered in a single phase (1-based). */
1649 static T ph_moment_from(const Matrix<T>& A, std::size_t start, unsigned k) {
1650 std::vector<T> alpha(A.rows(), num_traits<T>::from_int(0));
1651 alpha[start - 1] = num_traits<T>::from_int(1);
1652 return ph_moment(alpha, A, k);
1653 }
1654
1655 private:
1656 /** Solve (-A) x = b by Gaussian elimination with partial pivoting. */
1657 static std::vector<T> solve_neg(const Matrix<T>& A, const std::vector<T>& b) {
1658 const T zero = num_traits<T>::from_int(0);
1659 const std::size_t n = A.rows();
1660 if (A.cols() != n || b.size() != n)
1661 throw InputError("phase-type moment: the subgenerator is not square");
1662 Matrix<T> M(n, n, zero);
1663 std::vector<T> x = b;
1664 for (std::size_t i = 0; i < n; ++i)
1665 for (std::size_t j = 0; j < n; ++j) M(i, j) = T(-A(i, j));
1666 for (std::size_t col = 0; col < n; ++col) {
1667 std::size_t best = col;
1668 double bv = std::fabs(num_traits<T>::to_double(M(col, col)));
1669 for (std::size_t r = col + 1; r < n; ++r) {
1670 const double v = std::fabs(num_traits<T>::to_double(M(r, col)));
1671 if (v > bv) {
1672 bv = v;
1673 best = r;
1674 }
1675 }
1676 if (best != col) {
1677 for (std::size_t j = 0; j < n; ++j) std::swap(M(col, j), M(best, j));
1678 std::swap(x[col], x[best]);
1679 }
1680 if (M(col, col) == zero)
1681 throw NumericError("phase-type moment: the subgenerator is singular");
1682 for (std::size_t r = 0; r < n; ++r) {
1683 if (r == col) continue;
1684 const T f = T(M(r, col) / M(col, col));
1685 if (f == zero) continue;
1686 for (std::size_t j = 0; j < n; ++j) M(r, j) = T(M(r, j) - f * M(col, j));
1687 x[r] = T(x[r] - f * x[col]);
1688 }
1689 }
1690 for (std::size_t i = 0; i < n; ++i) x[i] = T(x[i] / M(i, i));
1691 return x;
1692 }
1693
1694 public:
1695
1696
1697 /**
1698 * The rate MATLAB's refreshRates would store: 1/mean, with the Immediate
1699 * singleton short-circuited to its declared rate so that the reciprocal of
1700 * 1e-8 is exactly 1e8 in every arithmetic rather than 1e8 plus rounding.
1701 */
1708};
1709
1710/**
1711 * What a `Prior` carries, in either of its two forms.
1712 *
1713 * DISCRETE: an explicit set of alternative distributions and their prior
1714 * weights, which must sum to one. This is the form the model.json wire carries
1715 * (`{"type":"Prior","distributions":[...],"probabilities":[...]}`), because a
1716 * factory cannot cross JSON.
1717 *
1718 * CONTINUOUS: a density over a scalar parameter theta plus a map theta ->
1719 * Distribution, the form the epistemic propagation of Trivedi and Bobbio
1720 * (2017), Sec. 3.4 needs. It is reduced to the discrete form by
1721 * `prior_discretize` (lang/prior.h) before anything downstream sees it, so both
1722 * forms are consumed identically. It can only be BUILT programmatically.
1723 *
1724 * IT IS NOT A MIXTURE. Each alternative is a separate model realization whose
1725 * weight is a prior probability over models, not a branching probability inside
1726 * one model. The mixture moments are still computed (see `Distrib::prior`)
1727 * because MATLAB's `Prior.getMean`/`getSCV` do, but they are a summary of the
1728 * epistemic uncertainty and not the law any station serves.
1729 */
1730template <class T>
1731struct PriorSpec {
1732 /** True for the parameter-density form, false for the alternative-set form. */
1733 bool continuous = false;
1734 /** The alternatives and their weights; the discrete form only. */
1735 std::vector<Distrib<T>> alternatives;
1736 std::vector<T> probabilities;
1737 /** The law of the scalar parameter; the continuous form only. */
1738 Distrib<T> param_dist;
1739 /**
1740 * theta -> Distribution; the continuous form only.
1741 *
1742 * A `std::function` and not a serializable description, exactly as MATLAB's
1743 * `distFactory` is a function handle: the map is arbitrary code (a rate
1744 * becomes an Exp, a scale becomes an Erlang of fixed order), and no wire
1745 * format in this codebase encodes it. That is why the JSON reader builds
1746 * the discrete form only.
1747 */
1748 std::function<Distrib<T>(const T&)> factory;
1749};
1750
1751} // namespace lang
1752} // namespace line
1753
1754#endif // LINE_LANG_LANG_TYPES_H
line::Cache
Cache(model, name, params).
Definition nodes.h:208
line::ClassSwitch
ClassSwitch(model, name, C).
Definition nodes.h:177
line::Delay
Delay(model, name): the infinite-server station.
Definition nodes.h:147
line::Fork
Fork(model, name).
Definition nodes.h:190
line::InputError::InputError
InputError(const std::string &what)
Definition error.h:39
line::Join
Join(model, name, fork).
Definition nodes.h:197
line::Matrix::cols
std::size_t cols() const
Definition matrix.h:90
line::Matrix::rows
std::size_t rows() const
Definition matrix.h:89
line::NumericError::NumericError
NumericError(const std::string &what)
Definition error.h:45
line::Place
Place(model, name): a Petri-net place.
Definition nodes.h:220
line::Queue
Queue(model, name, strategy).
Definition nodes.h:139
line::Router
Router(model, name): a stateless routing node.
Definition nodes.h:171
line::Sink
Sink(model, name): the external departure node, which holds no jobs.
Definition nodes.h:165
line::Source
Source(model, name): the external arrival station.
Definition nodes.h:153
line::Transition
Transition(model, name, params): a Petri-net transition.
Definition nodes.h:230
line::UnsupportedError::UnsupportedError
UnsupportedError(const std::string &what)
Definition error.h:51
error.h
The exception types the port throws.
matrix.h
Dense matrix and non-owning view.
line::lang::SchedStrategy
SchedStrategy
Scheduling disciplines, with the values of MATLAB SchedStrategy.
Definition lang_types.h:181
line::lang::DropStrategy
DropStrategy
Blocking and loss rules, with the values of MATLAB DropStrategy.
Definition lang_types.h:424
line::lang::TimingStrategy
TimingStrategy
SPN transition timing, with the values of MATLAB TimingStrategy.
Definition lang_types.h:361
line::lang::TimingStrategy::IMMEDIATE
@ IMMEDIATE
fires with zero delay, resolved by weight and priority
Definition lang_types.h:363
line::lang::TimingStrategy::TIMED
@ TIMED
fires after its firing distribution elapses
Definition lang_types.h:362
line::lang::BalkingStrategy
BalkingStrategy
Balking rules, with the values of MATLAB BalkingStrategy.
Definition lang_types.h:445
line::lang::SignalType
SignalType
G-network signal classes, with the values of MATLAB SignalType.
Definition lang_types.h:167
line::lang::SignalType::REPLY
@ REPLY
completes a synchronous call, releasing a held server
Definition lang_types.h:168
line::lang::SignalType::NEGATIVE
@ NEGATIVE
removes a batch of jobs (Gelenbe's negative customer)
Definition lang_types.h:169
line::lang::SignalType::CATASTROPHE
@ CATASTROPHE
removes EVERY job at the station
Definition lang_types.h:170
line::lang::JoinStrategy
JoinStrategy
Join rules, with the values of MATLAB JoinStrategy.
Definition lang_types.h:461
line::lang::PrecedenceType
PrecedenceType
Activity precedence kinds, with the values of MATLAB ActivityPrecedenceType.
Definition lang_types.h:470
line::lang::LqnElement
LqnElement
LQN element kinds, with the values of MATLAB LayeredNetworkElement.
Definition lang_types.h:464
line::lang::RoutingStrategy
RoutingStrategy
Routing strategies, with the values of MATLAB RoutingStrategy.
Definition lang_types.h:389
line::lang::RoutingStrategy::SDR
@ SDR
Krzesinski (1987) product-form state-dependent routing.
Definition lang_types.h:398
line::lang::sched_from_lqnx
SchedStrategy sched_from_lqnx(const std::string &s)
Parse the scheduling attribute of an .lqnx processor or task.
Definition lang_types.h:277
line::lang::DepartureDiscipline
DepartureDiscipline
When a Place releases a served token, MATLAB DepartureDiscipline.
Definition lang_types.h:458
line::lang::CallType
CallType
Call kinds, with the values of MATLAB CallType.
Definition lang_types.h:467
line::lang::metric_to_text
const char * metric_to_text(MetricType metric)
Port of MetricType.toText.
Definition lang_types.h:73
line::lang::RemovalPolicy
RemovalPolicy
Which job a negative signal removes, with the values of MATLAB RemovalPolicy.
Definition lang_types.h:174
line::lang::RemovalPolicy::LCFS
@ LCFS
the newest waiting job; servers only once nobody waits
Definition lang_types.h:177
line::lang::RemovalPolicy::RANDOM
@ RANDOM
uniform over waiting AND in-service jobs
Definition lang_types.h:175
line::lang::EventType
EventType
The events a state can undergo, with the values of MATLAB EventType.
Definition lang_types.h:111
line::lang::EventType::PHASE
@ PHASE
service advances a phase WITHOUT departing
Definition lang_types.h:116
line::lang::EventType::READ
@ READ
a cache item is read
Definition lang_types.h:117
line::lang::EventType::FAILURE
@ FAILURE
the server breaks down, going from up to down
Definition lang_types.h:126
line::lang::EventType::STAGE
@ STAGE
a random environment changes stage
Definition lang_types.h:118
line::lang::EventType::SWITCH
@ SWITCH
a polling server advances its switchover timer
Definition lang_types.h:125
line::lang::EventType::LOCAL
@ LOCAL
dummy event, no state change outside the node
Definition lang_types.h:113
line::lang::EventType::RENEGE
@ RENEGE
a waiting job abandons the queue (impatience)
Definition lang_types.h:123
line::lang::EventType::POST
@ POST
produce to a place or queue buffer
Definition lang_types.h:122
line::lang::EventType::DEP
@ DEP
a job departs
Definition lang_types.h:115
line::lang::EventType::START
@ START
a job begins or resumes holding a server
Definition lang_types.h:129
line::lang::EventType::ENABLE
@ ENABLE
an SPN mode becomes enabled
Definition lang_types.h:119
line::lang::EventType::FIRE
@ FIRE
an SPN mode fires
Definition lang_types.h:120
line::lang::EventType::RETRY
@ RETRY
an orbiting job retries entry at a retrial station
Definition lang_types.h:124
line::lang::EventType::REPAIR
@ REPAIR
the server is repaired, going from down to up, resuming the held job, which is why it emits no START
Definition lang_types.h:127
line::lang::EventType::ARV
@ ARV
a job arrives
Definition lang_types.h:114
line::lang::EventType::PRE
@ PRE
consume from a place or queue buffer, no server effect
Definition lang_types.h:121
line::lang::EventType::INIT
@ INIT
the model is initialized, t = 0
Definition lang_types.h:112
line::lang::EventType::PREEMPT
@ PREEMPT
a job holding a server is pushed back into the buffer
Definition lang_types.h:130
line::lang::MetricType
MetricType
Solver output metrics, with the numeric values of MATLAB MetricType.
Definition lang_types.h:46
line::lang::HeteroSchedPolicy
HeteroSchedPolicy
How a heterogeneous station picks among its server types, MATLAB HeteroSchedPolicy.
Definition lang_types.h:451
line::lang::PollingType
PollingType
Polling service disciplines, with the values of MATLAB PollingType.
Definition lang_types.h:370
line::lang::PollingType::KLIMITED
@ KLIMITED
serve at most K per visit (K in pollingPar)
Definition lang_types.h:373
line::lang::PollingType::EXHAUSTIVE
@ EXHAUSTIVE
serve until the queue empties
Definition lang_types.h:372
line::lang::PollingType::GATED
@ GATED
serve exactly the jobs present at the polling instant
Definition lang_types.h:371
line::lang::PollingType::DECREMENTING
@ DECREMENTING
serve until the queue is one shorter than at arrival
Definition lang_types.h:374
line::lang::JobClassType
JobClassType
Job class kinds, with the values of MATLAB JobClassType.
Definition lang_types.h:367
line::lang::sched_to_lqnx
std::string sched_to_lqnx(SchedStrategy s)
The scheduling attribute an .lqnx processor or task carries for a strategy.
Definition lang_types.h:304
line::lang::node_type_to_text
const char * node_type_to_text(NodeType t)
Name of a node kind, for diagnostics.
Definition lang_types.h:341
line::lang::ProcessType
ProcessType
Distribution kinds, with the values of MATLAB ProcessType.
Definition lang_types.h:483
line::lang::ProcessType::NORMAL
@ NORMAL
A Gaussian, and the ONE family whose value is not MATLAB's, because MATLAB has none to copy: ProcessT...
Definition lang_types.h:555
line::lang::ProcessType::NHPP
@ NHPP
The time-INHOMOGENEOUS families of Ko and Pender (ORL 45, 2017): an NHPP is a rate schedule lambda(t)...
Definition lang_types.h:536
line::lang::ProcessType::PRIOR
@ PRIOR
A Prior: a weighted set of ALTERNATIVE distributions, or a density over a scalar parameter plus a fac...
Definition lang_types.h:510
line::lang::GdScaling
std::function< std::vector< T >(const std::vector< T > &)> GdScaling
A globally state-dependent scaling, sn.gdscaling.
Definition lang_types.h:652
line::lang::sched_to_text
const char * sched_to_text(SchedStrategy s)
Definition lang_types.h:230
line::lang::process_to_text
const char * process_to_text(ProcessType p)
The MATLAB ProcessType name, as sn.procid prints it.
Definition lang_types.h:560
line::lang::process_is_markovian
bool process_is_markovian(ProcessType p)
ProcessType.isMarkovian: true when sn.proc carries an exact matrix representation of the law,...
Definition lang_types.h:610
line::lang::CdScaling
std::function< std::vector< T >(const std::vector< T > &)> CdScaling
A class-dependent scaling map, sn.cdscaling.
Definition lang_types.h:639
line::lang::routing_to_text
const char * routing_to_text(RoutingStrategy r)
Definition lang_types.h:402
line::lang::NodeType
NodeType
Node kinds, with the values of MATLAB NodeType.
Definition lang_types.h:324
line::lang::event_to_text
const char * event_to_text(EventType e)
Definition lang_types.h:137
line::lang::ReplacementStrategy
ReplacementStrategy
Cache replacement policies, with the values of MATLAB ReplacementStrategy.
Definition lang_types.h:378
line::lang::ReplacementStrategy::HLRU
@ HLRU
h-LRU / LRU(m): h lists, promote i -> i+1 on a hit
Definition lang_types.h:383
line::lang::ReplacementStrategy::CLIMB
@ CLIMB
move up one position on a hit (transposition rule)
Definition lang_types.h:384
line::lang::ReplacementStrategy::QLRU
@ QLRU
q-LRU: LRU with probabilistic admission on a miss
Definition lang_types.h:385
line::lang::ReplacementStrategy::LRU
@ LRU
least recently used
Definition lang_types.h:382
line::lang::ReplacementStrategy::FIFO
@ FIFO
first in, first out
Definition lang_types.h:380
line::lang::ImpatienceType
ImpatienceType
Impatience kinds, with the values of MATLAB ImpatienceType.
Definition lang_types.h:442
number.h
Number-type abstraction for the templated API port.
line::lang::Distrib::is_immediate
bool is_immediate() const
Definition lang_types.h:1584
line::lang::Distrib::sched_dist
static Distrib sched_dist(const std::vector< T > &breakpoints, const std::vector< Matrix< T > > &D0segs, const std::vector< Matrix< T > > &D1segs, bool cyclic, ProcessType tag)
The shared constructor of the three schedule families.
Definition lang_types.h:1070
line::lang::Distrib::empirical_cdf
static Distrib empirical_cdf(const std::vector< T > &x, const std::vector< T > &F)
EmpiricalCDF(x, F): the moments of the MIDPOINT rule over the CDF bins, which is what MATLAB Empirica...
Definition lang_types.h:1499
line::lang::Distrib::nhpp
static Distrib nhpp(const std::vector< T > &breakpoints, const std::vector< T > &rates, bool cyclic)
NHPP(breakpoints, rates, cyclic): a MAPt of ORDER ONE, which is what an inhomogeneous Poisson process...
Definition lang_types.h:1177
line::lang::Distrib::mu_vec
std::vector< T > mu_vec() const
Phase rates, MATLAB's getMu: the total outgoing rate of each phase.
Definition lang_types.h:1595
line::lang::Distrib::has_map
bool has_map() const
True when the type carries a (D0,D1) pair of its own.
Definition lang_types.h:1587
line::lang::Distrib< double >::replayer
static Distrib replayer(const std::vector< double > &samples)
Definition lang_types.h:1325
line::lang::Distrib< double >::sched_bp
std::vector< double > sched_bp
Definition lang_types.h:794
line::lang::Distrib::normal
static Distrib normal(const T &mu, const T &sigma)
Normal(mu, sigma): the Gaussian, for use as a continuous Prior's parameter density.
Definition lang_types.h:1294
line::lang::Distrib::dmap
static Distrib dmap(const Matrix< T > &D0, const Matrix< T > &D1)
DMAP(D0, D1): a DISCRETE-time MAP, where D0 + D1 is stochastic rather than a generator.
Definition lang_types.h:1548
line::lang::Distrib::exp_rate
static Distrib exp_rate(const T &r)
Definition lang_types.h:814
line::lang::Distrib::phase_type
static Distrib phase_type(const std::vector< T > &alpha, const Matrix< T > &A, bool acyclic)
PH / APH given by (alpha, A): D0 = A and D1 = (-A e) alpha.
Definition lang_types.h:1028
line::lang::Distrib< double >::Dmark
std::vector< Matrix< double > > Dmark
Definition lang_types.h:761
line::lang::Distrib::mapt
static Distrib mapt(const std::vector< T > &breakpoints, const std::vector< Matrix< T > > &D0segs, const std::vector< Matrix< T > > &D1segs, bool cyclic)
MAPt(breakpoints, {D0_k}, {D1_k}, cyclic): a piecewise-constant (D0(t), D1(t)).
Definition lang_types.h:1138
line::lang::Distrib< double >::declared
std::shared_ptr< Distrib< double > > declared
Definition lang_types.h:729
line::lang::Distrib::weibull
static Distrib weibull(const T &scale, const T &shape)
Definition lang_types.h:1243
line::lang::Distrib::replayer_from
static Distrib replayer_from(const std::vector< T > &samples, const std::string &path)
Replayer / Trace read FROM A FILE, which keeps the path beside the samples.
Definition lang_types.h:1318
line::lang::Distrib::has_schedule
bool has_schedule() const
Definition lang_types.h:797
line::lang::Distrib::bmap
static Distrib bmap(const std::vector< Matrix< T > > &D)
BMAP: the batch-size blocks D0, D1, ..., Dk, where Dj carries an arrival of batch size j.
Definition lang_types.h:1576
line::lang::Distrib::mmap
static Distrib mmap(const Matrix< T > &D0, const std::vector< Matrix< T > > &D1k)
MMAP: D0 plus one D1 block per mark.
Definition lang_types.h:1557
line::lang::Distrib::disabled_dist
static Distrib disabled_dist()
Definition lang_types.h:857
line::lang::Distrib::pht
static Distrib pht(const std::vector< T > &breakpoints, const std::vector< std::vector< T > > &alphas, const std::vector< Matrix< T > > &Ssegs, bool cyclic)
PHt(breakpoints, {alpha_k}, {S_k}, cyclic), stored as its equivalent MAP schedule: D0 = S and D1 = (-...
Definition lang_types.h:1148
line::lang::Distrib< double >::ph_moment
static double ph_moment(const std::vector< double > &alpha, const Matrix< double > &A, unsigned k)
Definition lang_types.h:1637
line::lang::Distrib::erlang_fit
static Distrib erlang_fit(const T &m, const T &c2)
Erlang fitted to a mean and an SCV, as MATLAB's Erlang.fitMeanAndSCV.
Definition lang_types.h:904
line::lang::Distrib::pareto
static Distrib pareto(const T &shape, const T &scale)
Pareto(shape, scale), with the MATLAB parameter order (alpha, k).
Definition lang_types.h:1207
line::lang::Distrib< double >::params
std::vector< double > params
Definition lang_types.h:734
line::lang::Distrib::gamma_dist
static Distrib gamma_dist(const T &shape, const T &scale)
Gamma(shape, scale), Weibull(scale, shape) and Lognormal(mu, sigma).
Definition lang_types.h:1231
line::lang::Distrib::cox2
static Distrib cox2(const T &mu1, const T &mu2, const T &phi1)
Cox2(mu1, mu2, phi1), MATLAB's two-phase Coxian constructor.
Definition lang_types.h:1013
line::lang::Distrib::map_dist
static Distrib map_dist(const Matrix< T > &D0, const Matrix< T > &D1, ProcessType tag)
A MAP given by its two matrices; the moments are those of its stationary phase.
Definition lang_types.h:1051
line::lang::Distrib::poisson
static Distrib poisson(const T &lambda)
Poisson(lambda), whose SCV is 1/lambda – the count's variance is lambda and its mean is lambda,...
Definition lang_types.h:1397
line::lang::Distrib< double >::trace
std::vector< double > trace
Definition lang_types.h:736
line::lang::Distrib::bernoulli
static Distrib bernoulli(const T &p)
Bernoulli(p): one trial, mean p and variance p(1-p).
Definition lang_types.h:1369
line::lang::Distrib< double >::ph_moment_from
static double ph_moment_from(const Matrix< double > &A, std::size_t start, unsigned k)
Definition lang_types.h:1649
line::lang::Distrib::hyperexp_n
static Distrib hyperexp_n(const std::vector< T > &p, const std::vector< T > &lambda)
HyperExp(p, lambda1, lambda2): phase i chosen with probability p_i.
Definition lang_types.h:931
line::lang::Distrib::discrete_uniform
static Distrib discrete_uniform(const T &a, const T &b)
DiscreteUniform(a, b) over the integers a..b inclusive.
Definition lang_types.h:1353
line::lang::Distrib< double >::prior
std::shared_ptr< PriorSpec< double > > prior
Definition lang_types.h:775
line::lang::Distrib::geometric
static Distrib geometric(const T &p)
Geometric(p) on the MATLAB convention: the NUMBER OF TRIALS to the first success, support {1,...
Definition lang_types.h:1412
line::lang::Distrib::det
static Distrib det(const T &m)
Definition lang_types.h:858
line::lang::Distrib::rap
static Distrib rap(const Matrix< T > &H0, const Matrix< T > &H1)
RAP(H0, H1): a rational arrival process, whose moments are the MAP ones.
Definition lang_types.h:1539
line::lang::Distrib::uniform
static Distrib uniform(const T &a, const T &b)
Uniform(a, b).
Definition lang_types.h:1193
line::lang::Distrib< double >::sched_D1
std::vector< Matrix< double > > sched_D1
Definition lang_types.h:795
line::lang::Distrib::lognormal
static Distrib lognormal(const T &logmean, const T &logsigma)
Definition lang_types.h:1266
line::lang::Distrib::hyperexp
static Distrib hyperexp(const T &p, const T &lambda1, const T &lambda2)
Definition lang_types.h:956
line::lang::Distrib< double >::sched_D0
std::vector< Matrix< double > > sched_D0
Definition lang_types.h:795
line::lang::Distrib::me
static Distrib me(const std::vector< T > &alpha, const Matrix< T > &A)
ME(alpha, A): the matrix-exponential distribution, whose moments are the phase-type ones – k!
Definition lang_types.h:1532
line::lang::Distrib::immediate
static Distrib immediate()
The Immediate singleton.
Definition lang_types.h:846
line::lang::Distrib::phi_vec
std::vector< T > phi_vec() const
Completion probabilities, MATLAB's getPhi: (D1 e) .
Definition lang_types.h:1607
line::lang::Distrib::erlang
static Distrib erlang(const T &phase_rate, std::size_t r)
Erlang(alpha, r): r phases of rate alpha, as MATLAB's Erlang(phaseRate, nphases).
Definition lang_types.h:873
line::lang::Distrib::exp_mean
static Distrib exp_mean(const T &m)
Definition lang_types.h:799
line::lang::Distrib::phases
std::size_t phases() const
The order of the representation, MATLAB's sn.phases.
Definition lang_types.h:1625
line::lang::Distrib::coxian
static Distrib coxian(const std::vector< T > &mu, const std::vector< T > &phi)
Coxian(mu, phi): phase i completes with probability phi(i) and otherwise moves to phase i+1.
Definition lang_types.h:987
line::lang::Distrib::is_prior
bool is_prior() const
Definition lang_types.h:776
line::lang::Distrib::discrete_sampler
static Distrib discrete_sampler(const std::vector< T > &p, const std::vector< T > &x)
DiscreteSampler(p, x): the pmf p over the points x.
Definition lang_types.h:1468
line::lang::Distrib::binomial
static Distrib binomial(const T &n, const T &p)
Binomial(n, p).
Definition lang_types.h:1381
line::lang::Distrib::zipf
static Distrib zipf(const T &s, std::size_t n)
Zipf(s, n) over the ranks 1..n, with the generalized harmonic moments H(s-1,n)/H(s,...
Definition lang_types.h:1430
line::lang::GlobalConstants
The MATLAB GlobalConstants, as reported by lineStart at its defaults.
Definition lang_types.h:667
line::lang::GlobalConstants::Immediate
static constexpr double Immediate
Rate of an Immediate distribution; its mean is 1/Immediate = 1e-8.
Definition lang_types.h:674
line::lang::GlobalConstants::FineTol
static constexpr double FineTol
Definition lang_types.h:668
line::lang::GlobalConstants::ArcTol
static constexpr double ArcTol
Below this an off-diagonal entry is NO ARC of the phase / state graph.
Definition lang_types.h:672
line::lang::GlobalConstants::Zero
static constexpr double Zero
Definition lang_types.h:670
line::lang::GlobalConstants::CoarseTol
static constexpr double CoarseTol
Definition lang_types.h:669
line::lang::GlobalConstants::MaxInt
static constexpr double MaxInt
Stand-in for an unbounded COUNT, MATLAB GlobalConstants.MaxInt.
Definition lang_types.h:679
line::lang::PriorSpec
A LINE Distribution, as the model layer and sn carry it.
Definition lang_types.h:1731
line::lang::PriorSpec::factory
std::function< Distrib< T >(const T &)> factory
theta -> Distribution; the continuous form only.
Definition lang_types.h:1748
line::lang::PriorSpec::continuous
bool continuous
True for the parameter-density form, false for the alternative-set form.
Definition lang_types.h:1733
line::lang::PriorSpec::param_dist
Distrib< T > param_dist
The law of the scalar parameter; the continuous form only.
Definition lang_types.h:1738
line::lang::PriorSpec::probabilities
std::vector< T > probabilities
Definition lang_types.h:1736
line::lang::PriorSpec::alternatives
std::vector< Distrib< T > > alternatives
The alternatives and their weights; the discrete form only.
Definition lang_types.h:1735
line::num_traits
Definition number.h:111