Fit a MAP(2) to three moments and an autocorrelation decay rate. More...

#include <cmath>
#include <cstddef>
#include "line/api/mam/map_dist.h"
#include "line/api/mam/map_fit_detail.h"
#include "line/api/mam/map_moment.h"
#include "line/api/mam/map_transform.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for map_block.h:

Namespaces
namespace	line
namespace	line::mam

Functions
template<class T>
Map< T >	line::mam::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.
template<class T>
Map< T >	line::mam::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.
template<class T>
Map< T >	line::mam::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.

Detailed Description

Fit a MAP(2) to three moments and an autocorrelation decay rate.

Templated port of matlab/lib/kpctoolbox/map/map_block.m and map_feasblock.m. The four rates of the general (non-MMPP) MAP(2)

D0 = diag(-mu00-mu01, -mu10-mu11), D1 = [mu00, mu01; mu10, mu11]

are the closed-form inverse of (E1, E2, E3, G2). That inverse is Maple output, some 57 KB of algebra across the four expressions.

THE ALGEBRA IS MACHINE-TRANSCRIBED, NOT RETYPED, by cpp/tools/matlab_expr_to_cpp.py, which parses the MATLAB expression, binds the repeated subexpressions to cseN temporaries, and re-emits it. Do not hand-edit the generated block; regenerate it.

IT MUST BE EVALUATED IN COMPLEX ARITHMETIC, and this is the trap. The single radicand the four expressions share goes NEGATIVE on perfectly feasible moment sets – (E1, SCV, E3/E3min, G2) = (1, 2, 2, 0.3) is one – and the imaginary parts then CANCEL in the four rates. MATLAB evaluates in complex arithmetic throughout and only afterwards asks whether any entry retains an imaginary part above 1e-4; that residual test, not the sign of the radicand, is what rejects a fit. An earlier version of this port branched on the radicand's sign and sent three of three feasible test cases to the fallback while MATLAB fitted all three exactly. The generated block is therefore instantiated at fitdetail::Cplx<T> and the reference's residual test is applied to the result.

THE FALLBACK IS THE REFERENCE'S, and it drops the third moment:

SCV >= 1: a hyperexponential-shaped MAP(2) matching E1 and E2 only, with the autocorrelation carried by the switching probability p = (1 - G2)/2 and D1 = -D0 P over P = [1-p, p; p, 1-p];
SCV < 1: the exponential of mean E1. The reference's commented-out general MAP(2) branch for SCV < 1 is NOT reinstated – it is commented out there, so reinstating it would answer a different model than every other codebase.

THE FALLBACK IS ITSELF BOUNDED ABOVE BY SCV 3. Its first branch rate is E1 (1 - sqrt((SCV-1)/2)), which vanishes at SCV = 3 and is negative above it. MATLAB returns the resulting non-generator regardless (measured in R2025a: an infinite diagonal at SCV = 3 with a NaN mean, a POSITIVE diagonal at SCV = 5, map_isfeasible 0 for both). This port refuses by name instead, since handing back a matrix that is not a MAP is worse than saying so.

map_feasblock is the same fit behind a moment repair: an SCV at or below one is raised to 1 + tol and a third moment below (3/2) E2^2 / E1 is raised to it.

ARITHMETIC: transcendental, for the radical.

Definition in file map_block.h.

Namespaces

Functions

Detailed Description