LINE Solver (C++)
Templated C++ port of the LINE queueing solver
Loading...
Searching...
No Matches
map_block.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_API_MAM_MAP_BLOCK_H
6#define LINE_API_MAM_MAP_BLOCK_H
7
8/**
9 * @file
10 * @ingroup api_mam
11 * Fit a MAP(2) to three moments and an autocorrelation decay rate.
12 *
13 * Templated port of matlab/lib/kpctoolbox/map/map_block.m and map_feasblock.m.
14 * The four rates of the general (non-MMPP) MAP(2)
15 *
16 * D0 = diag(-mu00-mu01, -mu10-mu11), D1 = [mu00, mu01; mu10, mu11]
17 *
18 * are the closed-form inverse of (E1, E2, E3, G2). That inverse is Maple output,
19 * some 57 KB of algebra across the four expressions.
20 *
21 * THE ALGEBRA IS MACHINE-TRANSCRIBED, NOT RETYPED, by
22 * `cpp/tools/matlab_expr_to_cpp.py`, which parses the MATLAB expression, binds
23 * the repeated subexpressions to `cseN` temporaries, and re-emits it. Do not
24 * hand-edit the generated block; regenerate it.
25 *
26 * IT MUST BE EVALUATED IN COMPLEX ARITHMETIC, and this is the trap. The single
27 * radicand the four expressions share goes NEGATIVE on perfectly feasible
28 * moment sets -- (E1, SCV, E3/E3min, G2) = (1, 2, 2, 0.3) is one -- and the
29 * imaginary parts then CANCEL in the four rates. MATLAB evaluates in complex
30 * arithmetic throughout and only afterwards asks whether any entry retains an
31 * imaginary part above 1e-4; that residual test, not the sign of the radicand,
32 * is what rejects a fit. An earlier version of this port branched on the
33 * radicand's sign and sent three of three feasible test cases to the fallback
34 * while MATLAB fitted all three exactly. The generated block is therefore
35 * instantiated at `fitdetail::Cplx<T>` and the reference's residual test is
36 * applied to the result.
37 *
38 * THE FALLBACK IS THE REFERENCE'S, and it drops the third moment:
39 * - SCV >= 1: a hyperexponential-shaped MAP(2) matching E1 and E2 only, with
40 * the autocorrelation carried by the switching probability p = (1 - G2)/2 and
41 * D1 = -D0 P over P = [1-p, p; p, 1-p];
42 * - SCV < 1: the exponential of mean E1. The reference's commented-out general
43 * MAP(2) branch for SCV < 1 is NOT reinstated -- it is commented out there,
44 * so reinstating it would answer a different model than every other codebase.
45 *
46 * THE FALLBACK IS ITSELF BOUNDED ABOVE BY SCV 3. Its first branch rate is
47 * E1 (1 - sqrt((SCV-1)/2)), which vanishes at SCV = 3 and is negative above it.
48 * MATLAB returns the resulting non-generator regardless (measured in R2025a: an
49 * infinite diagonal at SCV = 3 with a NaN mean, a POSITIVE diagonal at SCV = 5,
50 * `map_isfeasible` 0 for both). This port refuses by name instead, since handing
51 * back a matrix that is not a MAP is worse than saying so.
52 *
53 * `map_feasblock` is the same fit behind a moment repair: an SCV at or below one
54 * is raised to 1 + tol and a third moment below (3/2) E2^2 / E1 is raised to it.
55 *
56 * ARITHMETIC: transcendental, for the radical.
57 */
58
59#include <cmath>
60#include <cstddef>
61
62#include "line/api/mam/map_dist.h"
63#include "line/api/mam/map_fit_detail.h"
64#include "line/api/mam/map_moment.h"
65#include "line/api/mam/map_transform.h"
66#include "line/num/number.h"
67#include "line/util/error.h"
68#include "line/util/matrix.h"
69
70namespace line {
71namespace mam {
72
73namespace blockdetail {
74
75/** The four Maple rates, evaluated in complex arithmetic as MATLAB does. */
76template <class T>
77struct BlockRates {
78 fitdetail::Cplx<T> mu00, mu10, mu01, mu11;
79};
80
81template <class T>
82BlockRates<T> block_rates(const fitdetail::Cplx<T>& E1, const fitdetail::Cplx<T>& E2,
83 const fitdetail::Cplx<T>& E3, const fitdetail::Cplx<T>& G2) {
84 typedef fitdetail::Cplx<T> C;
85 // Local spellings the generated block uses; see the header note.
86 auto CI = [](long v) { return C(num_traits<T>::from_int(v)); };
87 auto CD = [](double v) { return C(num_traits<T>::from_double(v)); };
88 auto cpw = [](const C& b, unsigned e) {
89 C r(num_traits<T>::from_int(1));
90 for (unsigned i = 0; i < e; ++i) r = r * b;
91 return r;
92 };
93
94 // ---- BEGIN GENERATED (matlab_expr_to_cpp.py, map_block.m) ----
95 const C cse30 = (CI(216) * cpw(E1, 7));
96 const C cse20 = (cse30 * cpw(E2, 2));
97 const C cse21 = ((CI(24) * cpw(E1, 5)) * cpw(E3, 3));
98 const C cse24 = ((CI(6) * E1) * cpw(E3, 3));
99 const C cse28 = (CI(162) * cpw(E2, 7));
100 const C cse29 = (CI(567) * cpw(E1, 2));
101 const C cse31 = ((CI(270) * E1) * E3);
102 const C cse34 = (CI(24) * cpw(E1, 4));
103 const C cse35 = (CI(18) * cpw(E1, 2));
104 const C cse47 = (CI(9) * cpw(E2, 4));
105 const C cse18 = ((((((cse34 * E3) - ((CI(27) * cpw(E1, 3)) * cpw(E2, 2))) - ((cse35 * E2) * E3)) + ((CI(18) * E1) * cpw(E2, 3))) + (cpw(E3, 2) * E1)) + fitdetail::cplx_sqrt(((((((((((((((((((-CI(243)) * cpw(E1, 6)) * cpw(E2, 4)) - (cse24 * cpw(E2, 2))) + cse28) + (cse47 * cpw(E3, 2))) + cse21) + ((CI(648) * cpw(E1, 4)) * cpw(E2, 5))) - (cse29 * cpw(E2, 6))) + (cse20 * E3)) + (((CI(144) * cpw(E1, 6)) * cpw(E3, 2)) * E2)) - (((CI(756) * cpw(E1, 5)) * cpw(E2, 3)) * E3)) - (((CI(270) * cpw(E1, 4)) * cpw(E3, 2)) * cpw(E2, 2))) - (((CI(12) * cpw(E1, 3)) * cpw(E3, 3)) * E2)) + (((CI(810) * cpw(E1, 3)) * E3) * cpw(E2, 4))) + (((CI(108) * cpw(E1, 2)) * cpw(E3, 2)) * cpw(E2, 3))) - (cse31 * cpw(E2, 5))) + (cpw(E3, 4) * cpw(E1, 2)))));
106 const C cse19 = ((((((((((-CI(3)) * cpw(E3, 2)) * cpw(E2, 2)) + ((CI(48) * cpw(E3, 2)) * cpw(E1, 4))) + ((CI(81) * cpw(E2, 4)) * cpw(E1, 2))) - (((CI(36) * cpw(E3, 2)) * cpw(E1, 2)) * E2)) + (((CI(90) * cpw(E2, 3)) * E3) * E1)) - (((CI(126) * cpw(E2, 2)) * cpw(E1, 3)) * E3)) - (CI(54) * cpw(E2, 5))) + ((CI(2) * E1) * cpw(E3, 3)));
107 const C cse0 = (((((CI(9072) * cpw(E1, 9)) * cpw(E2, 2)) * E3) * cse18) / cse19);
108 const C cse1 = (((((CI(1512) * cpw(E1, 7)) * cpw(E2, 2)) * E3) * cse18) / cse19);
109 const C cse2 = (((((CI(30) * cpw(E1, 2)) * cpw(E2, 4)) * E3) * cse18) / cse19);
110 const C cse3 = (((((CI(12) * E1) * cpw(E3, 3)) * cse18) / cse19) * cpw(E2, 2));
111 const C cse4 = (((((CI(2) * E1) * cse18) / cse19) * cpw(E3, 2)) * cpw(E2, 3));
112 const C cse5 = ((((CI(3456) * cpw(E1, 10)) * cpw(E3, 2)) * cse18) / cse19);
113 const C cse6 = (((((CI(378) * E1) * E3) * cse18) / cse19) * cpw(E2, 5));
114 const C cse7 = ((((CI(576) * cpw(E1, 8)) * cpw(E3, 2)) * cse18) / cse19);
115 const C cse50 = (CI(4) * cpw(E1, 2));
116 const C cse8 = (((cse50 * cse18) / cse19) * cpw(E3, 4));
117 const C cse9 = ((((CI(9) * cpw(E3, 2)) * cse18) / cse19) * cpw(E2, 4));
118 const C cse10 = ((((CI(9) * E1) * cpw(E2, 6)) * cse18) / cse19);
119 const C cse17 = ((CI(3) * cse18) / cse19);
120 const C cse16 = (cse17 * E3);
121 const C cse11 = (cse16 * cpw(E2, 5));
122 const C cse12 = (((CI(5832) * cpw(E1, 8)) * cse18) / cse19);
123 const C cse13 = (((CI(486) * cse18) / cse19) * cpw(E2, 8));
124 const C cse14 = (((CI(162) * cse18) / cse19) * cpw(E2, 7));
125 const C cse15 = (((CI(972) * cpw(E1, 6)) * cse18) / cse19);
126 const C cse22 = (((CI(16) * cpw(E1, 5)) * E3) * E2);
127 const C cse23 = ((CI(27) * E1) * cpw(E2, 5));
128 const C cse25 = (cse47 * E3);
129 const C cse26 = (CI(1728) * cpw(E1, 10));
130 const C cse27 = (CI(1296) * cpw(E1, 9));
131 const C cse32 = (CI(648) * cpw(E1, 3));
132 const C cse33 = (CI(288) * cpw(E1, 8));
133 const C cse36 = (CI(16) * cpw(E1, 7));
134 const C cse37 = (CI(12) * cpw(E2, 4));
135 const C cse38 = (CI(36) * cpw(E1, 3));
136 const C cse39 = (CI(36) * cpw(E1, 5));
137 const C cse40 = (CI(84) * cpw(E1, 4));
138 const C cse41 = (CI(40) * cpw(E1, 5));
139 const C cse42 = (CI(72) * cpw(E1, 6));
140 const C cse43 = (CI(16) * cpw(E1, 3));
141 const C cse44 = (CI(32) * cpw(E1, 7));
142 const C cse45 = (CI(96) * cpw(E1, 5));
143 const C cse46 = (CI(72) * cpw(E1, 2));
144 const C cse48 = (CI(3) * cpw(E2, 2));
145 const C cse49 = (CI(3) * cpw(E2, 5));
146
147 const C mu00 = ((((-((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((((CI(6) * cpw(E3, 4)) * cpw(E2, 2)) * cse18) / cse19) + ((CI(6) * cpw(E3, 3)) * cpw(E2, 2))) + cse13) - ((((((CI(18) * E1) * cpw(E3, 3)) * cse18) / cse19) * cpw(E2, 3)) * G2)) + ((((((CI(5994) * cpw(E1, 3)) * E3) * cse18) / cse19) * cpw(E2, 5)) * G2)) - ((((((CI(810) * E1) * E3) * cse18) / cse19) * cpw(E2, 6)) * G2)) - ((((((CI(216) * cpw(E1, 5)) * cpw(E3, 3)) * cse18) / cse19) * E2) * G2)) + ((((((CI(108) * cpw(E1, 3)) * cpw(E3, 3)) * cse18) / cse19) * cpw(E2, 2)) * G2)) + ((((((CI(162) * cpw(E1, 2)) * cpw(E3, 2)) * cse18) / cse19) * cpw(E2, 4)) * G2)) + ((((((CI(20088) * cpw(E1, 7)) * cpw(E2, 3)) * E3) * cse18) / cse19) * G2)) - ((((((CI(7776) * cpw(E1, 8)) * cpw(E3, 2)) * cse18) / cse19) * E2) * G2)) - (cse0 * G2)) + (((((CI(16848) * cpw(E2, 4)) * E3) * cse18) / cse19) * cpw(E1, 5))) + (cse5 * G2)) + (((((CI(9504) * cpw(E1, 8)) * cpw(E3, 2)) * cse18) / cse19) * E2)) + cse0) + (((((CI(3024) * cpw(E1, 2)) * cpw(E3, 2)) * cse18) / cse19) * cpw(E2, 4))) - (((((CI(288) * E1) * cpw(E3, 3)) * cse18) / cse19) * cpw(E2, 3))) + ((((cse32 * E3) * cse18) / cse19) * cpw(E2, 5))) - (((((CI(2430) * E1) * E3) * cse18) / cse19) * cpw(E2, 6))) + ((((((CI(6264) * cpw(E1, 6)) * cpw(E3, 2)) * cse18) / cse19) * cpw(E2, 2)) * G2)) - ((((((CI(2052) * cpw(E1, 4)) * cpw(E3, 2)) * cse18) / cse19) * cpw(E2, 3)) * G2)) - (((((CI(828) * cpw(E1, 3)) * cpw(E3, 3)) * cse18) / cse19) * cpw(E2, 2))) - (((((CI(3645) * cse18) / cse19) * cpw(E2, 7)) * cpw(E1, 2)) * G2)) - (((((CI(3024) * cpw(E1, 6)) * cpw(E3, 2)) * cse18) / cse19) * cpw(E2, 2))) - (((((CI(5940) * cpw(E1, 4)) * cpw(E3, 2)) * cse18) / cse19) * cpw(E2, 3))) + (((((CI(3024) * cpw(E1, 5)) * cpw(E3, 3)) * cse18) / cse19) * E2)) + (((((CI(10206) * cpw(E1, 4)) * cse18) / cse19) * cpw(E2, 6)) * G2)) + ((cse12 * G2) * cpw(E2, 4))) - (((((CI(12636) * cpw(E1, 6)) * cse18) / cse19) * G2) * cpw(E2, 5))) - ((((((CI(16524) * cpw(E2, 4)) * E3) * cse18) / cse19) * cpw(E1, 5)) * G2)) - (((((CI(24624) * cpw(E1, 7)) * cpw(E2, 3)) * E3) * cse18) / cse19)) + (cse13 * G2)) - ((((CI(168) * cpw(E1, 4)) * cse18) / cse19) * cpw(E3, 4))) + ((((CI(135) * cpw(E3, 2)) * cse18) / cse19) * cpw(E2, 5))) - ((((CI(1872) * cpw(E1, 7)) * cpw(E3, 3)) * cse18) / cse19)) - ((((CI(12150) * cpw(E1, 4)) * cse18) / cse19) * cpw(E2, 6))) - (cse12 * cpw(E2, 4))) - ((((CI(4) * E1) * cse18) / cse19) * cpw(E3, 5))) + ((((CI(2187) * cse18) / cse19) * cpw(E2, 7)) * cpw(E1, 2))) + ((((CI(15552) * cpw(E1, 6)) * cse18) / cse19) * cpw(E2, 5))) + (((((CI(132) * cpw(E1, 2)) * cse18) / cse19) * cpw(E3, 4)) * E2)) + (((((CI(27) * cpw(E3, 2)) * cse18) / cse19) * cpw(E2, 5)) * G2)) + (((((CI(144) * cpw(E1, 7)) * cpw(E3, 3)) * cse18) / cse19) * G2)) - cse5) + (((CI(27) * E3) * cpw(E2, 5)) * G2)) - ((CI(5184) * cpw(E1, 5)) * cpw(E2, 4))) + ((CI(4536) * cpw(E1, 7)) * cpw(E2, 3))) + ((CI(576) * cpw(E1, 7)) * cpw(E3, 2))) + ((CI(81) * E3) * cpw(E2, 5))) + (cse26 * E3)) - (cse27 * cpw(E2, 2))) - ((cse26 * G2) * E3)) + ((cse27 * G2) * cpw(E2, 2))) - (((CI(5616) * cpw(E1, 8)) * E2) * E3)) - (((CI(2592) * cpw(E1, 7)) * cpw(E2, 3)) * G2)) + (((CI(5832) * cpw(E1, 6)) * cpw(E2, 2)) * E3)) - (((CI(1080) * cpw(E1, 5)) * cpw(E3, 2)) * E2)) + (((CI(1944) * cpw(E1, 5)) * G2) * cpw(E2, 4))) - (((CI(2160) * cpw(E1, 4)) * cpw(E2, 3)) * E3)) + ((cse32 * cpw(E3, 2)) * cpw(E2, 2))) - ((cse32 * cpw(E2, 5)) * G2)) - (((CI(24) * cpw(E1, 2)) * cpw(E3, 3)) * E2)) + (((CI(54) * cpw(E1, 2)) * E3) * cpw(E2, 4))) - (((CI(126) * E1) * cpw(E3, 2)) * cpw(E2, 3))) + (((CI(81) * E1) * cpw(E2, 6)) * G2)) + (cse34 * cpw(E3, 3))) + ((CI(2430) * cpw(E1, 3)) * cpw(E2, 5))) - ((CI(405) * E1) * cpw(E2, 6))) + ((((CI(3888) * cpw(E1, 8)) * E2) * G2) * E3)) - ((((CI(3456) * cpw(E1, 6)) * E3) * G2) * cpw(E2, 2))) + ((((CI(1512) * cpw(E1, 4)) * E3) * cpw(E2, 3)) * G2)) - ((((CI(324) * cpw(E1, 2)) * E3) * G2) * cpw(E2, 4)))) / (((CI(2) * E3) * E1) - cse48)) / ((((((CI(24) * E3) * cpw(E1, 3)) - ((CI(27) * cpw(E2, 2)) * cpw(E1, 2))) - (((CI(18) * E3) * E2) * E1)) + cpw(E3, 2)) + (CI(18) * cpw(E2, 3)))) / ((((((((((-CI(12)) * cpw(E1, 4)) + ((((CI(24) * cse18) / cse19) * cpw(E1, 4)) * E3)) - ((((CI(36) * cse18) / cse19) * cpw(E1, 3)) * cpw(E2, 2))) - (((((CI(18) * cse18) / cse19) * cpw(E1, 2)) * E2) * E3)) + ((CI(12) * E2) * cpw(E1, 2))) + ((((CI(2) * cse18) / cse19) * cpw(E3, 2)) * E1)) + ((((CI(27) * cse18) / cse19) * E1) * cpw(E2, 3))) - (cse16 * cpw(E2, 2))) - cse48));
148
149 const C mu10 = ((cse17 * ((CI(2) * cpw(E1, 2)) - E2)) * (G2 - CI(1)));
150
151 const C mu01 = ((CI(9) * ((((((((((((((((((((((((((((((((((cse36 * E3) - ((cse37 * G2) * cpw(E1, 2))) - ((CI(12) * cpw(E2, 3)) * cpw(E1, 4))) + (cse49 * G2)) - cse49) + cse4) + (cse10 * G2)) + (cse11 * G2)) + (cse37 * cpw(E1, 2))) - ((((cse38 * cpw(E2, 5)) * cse18) / cse19) * G2)) + ((((cse39 * cpw(E2, 4)) * cse18) / cse19) * G2)) + cse2) - ((((cse40 * cse18) / cse19) * E3) * cpw(E2, 3))) + ((((cse41 * cse18) / cse19) * cpw(E3, 2)) * E2)) + (((cse38 * cse18) / cse19) * cpw(E2, 5))) - cse10) + ((((cse42 * cse18) / cse19) * E3) * cpw(E2, 2))) - ((((cse43 * cpw(E3, 2)) * cpw(E2, 2)) * cse18) / cse19)) + ((((cse44 * cpw(E3, 2)) * cse18) / cse19) * G2)) - (cse4 * G2)) - (cse2 * G2)) - (((cse44 * cse18) / cse19) * cpw(E3, 2))) + (((((cse43 * cpw(E2, 2)) * cpw(E3, 2)) * cse18) / cse19) * G2)) - (((((cse41 * E2) * cpw(E3, 2)) * cse18) / cse19) * G2)) - (((cse39 * cse18) / cse19) * cpw(E2, 4))) - cse11) + (((((cse40 * cpw(E2, 3)) * E3) * cse18) / cse19) * G2)) - (((((cse42 * cpw(E2, 2)) * E3) * cse18) / cse19) * G2)) + (((CI(4) * cpw(E2, 2)) * cpw(E1, 3)) * E3)) - ((cse36 * G2) * E3)) - cse22) + (((CI(12) * G2) * cpw(E2, 3)) * cpw(E1, 4))) + (cse22 * G2)) - ((((CI(4) * cpw(E1, 3)) * G2) * E3) * cpw(E2, 2)))) / (((((((((((((((((((((-CI(48)) * cpw(E1, 5)) * cpw(E3, 2)) + (((cse45 * cpw(E3, 3)) * cse18) / cse19)) + (((CI(108) * cpw(E1, 4)) * E3) * cpw(E2, 2))) - (((((CI(396) * cpw(E1, 4)) * cpw(E3, 2)) * cse18) / cse19) * cpw(E2, 2))) + ((cse38 * cpw(E3, 2)) * E2)) - (((((CI(72) * cpw(E1, 3)) * cpw(E3, 3)) * cse18) / cse19) * E2)) + (((((CI(540) * cpw(E2, 4)) * E3) * cse18) / cse19) * cpw(E1, 3))) - ((CI(54) * cpw(E1, 3)) * cpw(E2, 4))) - ((cse46 * cpw(E2, 3)) * E3)) - ((((CI(243) * cpw(E1, 2)) * cse18) / cse19) * cpw(E2, 6))) + (((((CI(288) * cpw(E1, 2)) * cpw(E3, 2)) * cse18) / cse19) * cpw(E2, 3))) + cse8) - (((CI(6) * E1) * cpw(E3, 2)) * cpw(E2, 2))) - cse6) - cse3) + cse23) + cse9) + cse14) + cse25));
152
153 const C mu11 = ((-((((((((((((((((((((((((((((((((((((((((((((((((((((((((-CI(60)) * cpw(E1, 3)) * cpw(E3, 2)) * E2) + (((CI(198) * cpw(E1, 2)) * cpw(E2, 3)) * E3)) - (((CI(288) * cpw(E1, 4)) * E3) * cpw(E2, 2))) + cse25) - ((cse47 * G2) * E3)) - ((((CI(324) * cpw(E1, 4)) * G2) * E3) * cpw(E2, 2))) + ((((CI(90) * cpw(E1, 2)) * G2) * cpw(E2, 3)) * E3)) - ((CI(189) * E1) * cpw(E2, 5))) - (((CI(12) * E1) * cpw(E3, 2)) * cpw(E2, 2))) - ((cse33 * G2) * E3)) + ((cse30 * G2) * cpw(E2, 2))) - (((CI(216) * cpw(E1, 6)) * E3) * E2)) - (((CI(324) * cpw(E1, 5)) * cpw(E2, 3)) * G2)) + (((CI(162) * cpw(E1, 3)) * G2) * cpw(E2, 4))) - (cse23 * G2)) - ((((CI(120) * cpw(E1, 5)) * cpw(E3, 3)) * cse18) / cse19)) - cse14) + (((((CI(84) * cpw(E1, 3)) * cpw(E3, 3)) * cse18) / cse19) * E2)) - ((((CI(81) * cpw(E1, 2)) * cse18) / cse19) * cpw(E2, 6))) - cse8) - cse9) - cse7) - (cse15 * cpw(E2, 4))) + (((((CI(216) * cpw(E1, 4)) * cpw(E3, 2)) * cse18) / cse19) * cpw(E2, 2))) - (((((CI(306) * cpw(E1, 2)) * cpw(E3, 2)) * cse18) / cse19) * cpw(E2, 3))) + cse6) + cse3) + ((((CI(1134) * cpw(E1, 4)) * cse18) / cse19) * cpw(E2, 5))) - (((cse28 * cse18) / cse19) * G2)) - ((((((CI(1458) * cpw(E1, 3)) * E3) * cse18) / cse19) * G2) * cpw(E2, 4))) - (((((cse46 * cpw(E3, 2)) * cse18) / cse19) * G2) * cpw(E2, 3))) + ((((cse24 * cse18) / cse19) * G2) * cpw(E2, 2))) + ((cse15 * G2) * cpw(E2, 4))) + (((cse21 * cse18) / cse19) * G2)) - (((((CI(1836) * cpw(E1, 5)) * cpw(E2, 3)) * E3) * cse18) / cse19)) - (((((CI(1620) * cpw(E1, 4)) * cse18) / cse19) * G2) * cpw(E2, 5))) + (((((CI(891) * cpw(E1, 2)) * cse18) / cse19) * cpw(E2, 6)) * G2)) - (((((CI(9) * cse18) / cse19) * cpw(E3, 2)) * cpw(E2, 4)) * G2)) + (cse7 * G2)) + ((((((CI(540) * cpw(E1, 4)) * cpw(E3, 2)) * cse18) / cse19) * G2) * cpw(E2, 2))) - ((((((CI(24) * cpw(E1, 3)) * cpw(E3, 3)) * cse18) / cse19) * G2) * E2)) + ((((cse31 * cse18) / cse19) * G2) * cpw(E2, 5))) + cse1) + (((((CI(720) * cpw(E1, 6)) * E2) * cpw(E3, 2)) * cse18) / cse19)) + ((((((CI(2592) * cpw(E1, 5)) * cpw(E2, 3)) * E3) * cse18) / cse19) * G2)) - (cse1 * G2)) - ((((((CI(1008) * cpw(E1, 6)) * E2) * cpw(E3, 2)) * cse18) / cse19) * G2)) + ((((CI(504) * cpw(E1, 6)) * E3) * E2) * G2)) + (cse45 * cpw(E3, 2))) + ((CI(378) * cpw(E1, 3)) * cpw(E2, 4))) + (cse33 * E3)) - cse20) + (cse50 * cpw(E3, 3)))) / ((((((((((((((((((((((((-CI(96)) * cpw(E3, 2)) * cpw(E1, 6)) + ((((CI(192) * cpw(E1, 6)) * cpw(E3, 3)) * cse18) / cse19)) - (((((CI(648) * cpw(E1, 5)) * cse18) / cse19) * cpw(E3, 2)) * cpw(E2, 2))) + (((CI(108) * E3) * cpw(E2, 2)) * cpw(E1, 5))) + (((((CI(702) * cpw(E1, 4)) * cse18) / cse19) * E3) * cpw(E2, 4))) - (((((CI(192) * cpw(E1, 4)) * cpw(E3, 3)) * cse18) / cse19) * E2)) + (((CI(72) * cpw(E3, 2)) * cpw(E1, 4)) * E2)) - ((CI(4) * cpw(E3, 3)) * cpw(E1, 3))) + ((((CI(8) * cpw(E1, 3)) * cse18) / cse19) * cpw(E3, 4))) - ((((CI(243) * cpw(E1, 3)) * cse18) / cse19) * cpw(E2, 6))) + (((((CI(594) * cpw(E1, 3)) * cse18) / cse19) * cpw(E3, 2)) * cpw(E2, 3))) + ((((cse35 * cpw(E3, 3)) * cse18) / cse19) * cpw(E2, 2))) - ((CI(81) * cpw(E2, 5)) * cpw(E1, 2))) - ((((cse29 * cpw(E2, 5)) * E3) * cse18) / cse19)) - (((((CI(2) * E1) * cpw(E3, 4)) * E2) * cse18) / cse19)) - (((((CI(81) * E1) * cse18) / cse19) * cpw(E3, 2)) * cpw(E2, 4))) + ((((CI(162) * E1) * cpw(E2, 7)) * cse18) / cse19)) - (((CI(54) * E3) * cpw(E2, 4)) * E1)) + ((((CI(54) * cse18) / cse19) * E3) * cpw(E2, 6))) + (CI(54) * cpw(E2, 6))) + ((CI(3) * cpw(E2, 3)) * cpw(E3, 2))) + ((((CI(3) * cpw(E3, 3)) * cse18) / cse19) * cpw(E2, 3))));
154 // ---- END GENERATED ----
155
156 BlockRates<T> r;
157 r.mu00 = mu00;
158 r.mu10 = mu10;
159 r.mu01 = mu01;
160 r.mu11 = mu11;
161 return r;
162}
163
164/**
165 * The reference's fallback when the exact fit is infeasible: match E1 and E2
166 * only, carrying G2 in the switching probability.
167 */
168template <class T>
169Map<T> block_fallback(const T& E1, const T& E2, const T& G2) {
170 using fitdetail::num_sqrt;
171 const T one = num_traits<T>::from_int(1), two = num_traits<T>::from_int(2);
172 const T zero = num_traits<T>::from_int(0), four = num_traits<T>::from_int(4);
173 const T SCV = T((E2 - E1 * E1) / (E1 * E1));
174
175 if (num_traits<T>::to_double(SCV) < 1.0) {
176 // SCV < 1: the reference returns the exponential; its general MAP(2)
177 // branch is commented out there and is not reinstated here.
178 Map<T> m;
179 m.D0 = Matrix<T>(1, 1, T(-one / E1));
180 m.D1 = Matrix<T>(1, 1, T(one / E1));
181 return m;
182 }
183
184 // disc = 2 E1^2 (SCV - 1), so mu1 = E1 (1 - sqrt((SCV-1)/2)).
185 const T disc = T(-four * E1 * E1 + two * E2);
186 const T r = num_sqrt(disc > zero ? disc : zero);
187 const T mu1 = T(E1 - r / two);
188 const T mu2 = T(E1 + r / two);
189 if (!(num_traits<T>::to_double(mu1) > 0.0))
190 throw InputError(
191 "map_block: the moment set is infeasible for the exact fit, and the reference's "
192 "hyperexponential fallback is defined only below SCV 3, where its first branch rate "
193 "E1 (1 - sqrt((SCV-1)/2)) is still positive. MATLAB returns the resulting "
194 "non-generator regardless; refusing here rather than handing back a matrix that is "
195 "not a MAP. Supply a feasible third moment, or fit with map_mmpp2");
196 const T p = T(one / two - G2 / two);
197
198 Matrix<T> D0(2, 2, zero);
199 D0(0, 0) = T(-one / mu1);
200 D0(1, 1) = T(-one / mu2);
201 Matrix<T> P(2, 2, zero);
202 P(0, 0) = T(one - p);
203 P(0, 1) = p;
204 P(1, 0) = p;
205 P(1, 1) = T(one - p);
206
207 Map<T> m;
208 m.D0 = D0;
209 m.D1 = Matrix<T>(2, 2, zero);
210 for (std::size_t i = 0; i < 2; ++i)
211 for (std::size_t j = 0; j < 2; ++j) {
212 T acc = zero;
213 for (std::size_t q = 0; q < 2; ++q) acc += -D0(i, q) * P(q, j);
214 m.D1(i, j) = acc;
215 }
216 return m;
217}
218
219} // namespace blockdetail
220
221/**
222 * @brief Fit a MAP(2) to three moments and an autocorrelation decay rate.
223 *
224 * @param E1 first moment
225 * @param E2 second moment
226 * @param E3 third moment
227 * @param G2 autocorrelation decay rate rho(i)/rho(i-1)
228 */
229template <class T>
230Map<T> map_block(const T& E1, const T& E2, const T& E3, const T& G2) {
231 static_assert(num_traits<T>::has_transcendental,
232 "map_block inverts the moment equations through a radical");
233 typedef fitdetail::Cplx<T> C;
234 const T zero = num_traits<T>::from_int(0);
235
236 const blockdetail::BlockRates<T> r =
237 blockdetail::block_rates<T>(C(E1), C(E2), C(E3), C(G2));
238
239 // The reference's own test: a residual imaginary part above 1e-4 anywhere
240 // means the fit is not real and is discarded.
241 const C* all[4] = {&r.mu00, &r.mu10, &r.mu01, &r.mu11};
242 for (int i = 0; i < 4; ++i) {
243 const double im = num_traits<T>::to_double(all[i]->im);
244 if (!(std::fabs(im) <= 1e-4) || std::isnan(im)) return blockdetail::block_fallback(E1, E2, G2);
245 }
246 const T mu00 = r.mu00.re, mu10 = r.mu10.re, mu01 = r.mu01.re, mu11 = r.mu11.re;
247 for (int i = 0; i < 4; ++i)
248 if (std::isnan(num_traits<T>::to_double(all[i]->re)))
249 return blockdetail::block_fallback(E1, E2, G2);
250 if (mu00 < zero || mu11 < zero || mu01 < zero || mu10 < zero)
251 return blockdetail::block_fallback(E1, E2, G2);
252
253 Map<T> m;
254 m.D0 = Matrix<T>(2, 2, zero);
255 m.D1 = Matrix<T>(2, 2, zero);
256 m.D0(0, 0) = T(-mu00 - mu01);
257 m.D0(1, 1) = T(-mu10 - mu11);
258 m.D1(0, 0) = mu00;
259 m.D1(0, 1) = mu01;
260 m.D1(1, 0) = mu10;
261 m.D1(1, 1) = mu11;
262 if (!map_isfeasible(m)) return blockdetail::block_fallback(E1, E2, G2);
263 return m;
264}
265
266/** `map_block` with the SCV spelling of the second argument. */
267template <class T>
268Map<T> map_block_scv(const T& E1, const T& SCV, const T& E3, const T& G2) {
269 return map_block(E1, T((num_traits<T>::from_int(1) + SCV) * E1 * E1), E3, G2);
270}
271
272/**
273 * `map_feasblock`: repair the moments into the feasible region, then fit.
274 *
275 * An E2 at or below the exponential value makes the SCV non-positive, and an E3
276 * below (3/2) E2^2 / E1 is outside what any MAP(2) admits; both are raised to
277 * their limits plus a tolerance, exactly as the reference does.
278 */
279template <class T>
280Map<T> map_feasblock(const T& E1, const T& E2_in, const T& E3_in, const T& G2) {
281 const T one = num_traits<T>::from_int(1), two = num_traits<T>::from_int(2);
282 const T tol = num_traits<T>::from_double(1e-10);
283 T E2 = E2_in, E3 = E3_in;
284
285 // The exponential boundary E2 == 2 E1^2 is returned as the Poisson process,
286 // rescaled to the requested mean.
287 if (E2 == two * E1 * E1) {
288 Map<T> m;
289 m.D0 = Matrix<T>(2, 2, num_traits<T>::from_int(0));
290 m.D1 = Matrix<T>(2, 2, num_traits<T>::from_int(0));
291 m.D0(0, 0) = -one;
292 m.D0(1, 1) = -one;
293 m.D1(0, 0) = num_traits<T>::from_rational(1, 2);
294 m.D1(0, 1) = num_traits<T>::from_rational(1, 2);
295 m.D1(1, 0) = num_traits<T>::from_rational(1, 2);
296 m.D1(1, 1) = num_traits<T>::from_rational(1, 2);
297 return map_scale(m, E1);
298 }
299 if (E2 <= two * E1 * E1) E2 = T((two + tol) * E1 * E1);
300 const T e3min = T(num_traits<T>::from_rational(3, 2) * E2 * E2 / E1);
301 if (E3 <= e3min) E3 = T((num_traits<T>::from_rational(3, 2) + tol) * E2 * E2 / E1);
302 return map_block(E1, E2, E3, G2);
303}
304
305} // namespace mam
306} // namespace line
307
308#endif // LINE_API_MAM_MAP_BLOCK_H
line::InputError::InputError
InputError(const std::string &what)
Definition error.h:39
error.h
The exception types the port throws.
map_dist.h
Analytic distances between continuous-time MAPs.
map_fit_detail.h
Scalar helpers shared by the MAP/PH moment-matching headers.
map_moment.h
Markovian arrival process descriptors: stationary vectors, rate, moments, autocorrelation and the ind...
map_transform.h
MAP constructors and structural transformations.
matrix.h
Dense matrix and non-owning view.
line::mam::map_block_scv
Map< T > map_block_scv(const T &E1, const T &SCV, const T &E3, const T &G2)
map_block with the SCV spelling of the second argument.
Definition map_block.h:268
line::mam::map_block
Map< T > map_block(const T &E1, const T &E2, const T &E3, const T &G2)
Fit a MAP(2) to three moments and an autocorrelation decay rate.
Definition map_block.h:230
line::mam::map_isfeasible
bool map_isfeasible(const Map< T > &m, const T &tol)
Structural feasibility of a MAP within a tolerance (map_isfeasible.m): off-diagonal D0 and all of D1 ...
Definition map_transform.h:498
line::mam::map_scale
Map< T > map_scale(const Map< T > &in, const T &new_mean)
Rescale time so that the mean inter-arrival time becomes new_mean.
Definition map_transform.h:75
line::mam::map_feasblock
Map< T > map_feasblock(const T &E1, const T &E2_in, const T &E3_in, const T &G2)
map_feasblock: repair the moments into the feasible region, then fit.
Definition map_block.h:280
number.h
Number-type abstraction for the templated API port.
line::mam::Map
A MAP as the pair of matrices (D0, D1).
Definition map_moment.h:53
line::mam::Map::D1
Matrix< T > D1
Definition map_moment.h:55
line::mam::Map::D0
Matrix< T > D0
Definition map_moment.h:54
line::num_traits
Definition number.h:111