Classes
struct	WeeksParams
	A Laguerre expansion: the damping, the scaling and the coefficients. More...

Typedefs
using	Cplx = std::complex<double>
using	LaplaceFn = std::function<Cplx(Cplx)>
	The transform, evaluated at complex argument.
using	RealLaplaceFn = std::function<double(double)>
	A transform that can only be evaluated on the real axis.

Enumerations
enum class	LaplaceMethod { Euler = 0 , Talbot , GaverStehfest , Cme , Weeks }
	The methods laplace_invert accepts. More...

Functions
std::vector< Cplx >	euler_get_alpha (std::size_t n)
	Euler nodes: a vertical line at Re = (n-1) log(10) / 6.
std::vector< double >	euler_get_eta (std::size_t n)
	Euler weights before the alternating sign and the scale.
std::vector< Cplx >	euler_get_omega (std::size_t n)
	Euler weights: eta, alternating in sign, scaled by 10^((n-1)/6).
std::vector< Cplx >	talbot_get_alpha (std::size_t n)
	Talbot nodes: the cotangent contour, bending into the left half plane.
std::vector< Cplx >	talbot_get_omega (std::size_t n, const std::vector< Cplx > &alpha)
	Talbot weights, which carry the contour's own derivative.
std::vector< double >	gaver_stehfest_get_alpha (std::size_t n)
	Gaver-Stehfest nodes: k log 2, on the REAL axis.
std::vector< double >	gaver_stehfest_get_omega (std::size_t n)
	Gaver-Stehfest weights.
double	laplace_invert_euler (const LaplaceFn &F, double t, std::size_t n=41)
	Euler inversion of F at t.
double	laplace_invert_talbot (const LaplaceFn &F, double t, std::size_t n=32)
	Talbot inversion of F at t.
double	laplace_invert_gaver_stehfest (const RealLaplaceFn &F, double t, std::size_t n=12)
	Gaver-Stehfest inversion of F at t.
std::vector< double >	laplace_weeks_coeffs (const LaplaceFn &F, double sigma=0.0, double b=1.0, std::size_t p0=200)
	Laguerre coefficients q_n, n = 0..2*p0-1, of f_{sigma,b}(t) = exp(-sigma t) f(t/b), whose generating function is.
WeeksParams	laplace_weeks_scaling (const LaplaceFn &F, std::size_t p0=200, double tol=1e-10)
	The automatic (sigma, b) search of Fig.
double	laplace_invert_weeks (const WeeksParams &w, double t)
	Invert by the Laguerre series f(t) = sum_n q_n l_n(t), recovered as exp(sigmabt) f_{sigma,b}(b*t).
double	laplace_invert_weeks (const LaplaceFn &F, double t, std::size_t p0=200)
	Convenience overload: build the parameters, then invert at one point.
LaplaceMethod	laplace_method (const std::string &s)
	Parse the reference's method names, including its two Gaver spellings.
double	laplace_invert (const LaplaceFn &F, double t, LaplaceMethod method=LaplaceMethod::Euler, std::size_t n=0)
	Invert F at t by the named method.
std::vector< double >	laplace_invert_pdf (const LaplaceFn &F, const std::vector< double > &t, LaplaceMethod method=LaplaceMethod::Euler, std::size_t n=0)
	The DENSITY on a grid: the inversion clamped at zero.
std::vector< double >	laplace_invert_cdf (const LaplaceFn &F, const std::vector< double > &t, LaplaceMethod method=LaplaceMethod::Euler, std::size_t n=0)
	The DISTRIBUTION on a grid, from the transform of the DENSITY.

Typedef Documentation

◆ Cplx

using line::lti::Cplx = std::complex<double>

Definition at line 62 of file laplace_invert.h.

◆ LaplaceFn

using line::lti::LaplaceFn = std::function<Cplx(Cplx)>

The transform, evaluated at complex argument.

Definition at line 65 of file laplace_invert.h.

◆ RealLaplaceFn

using line::lti::RealLaplaceFn = std::function<double(double)>

A transform that can only be evaluated on the real axis.

Definition at line 67 of file laplace_invert.h.

Enumeration Type Documentation

◆ LaplaceMethod

enum class line::lti::LaplaceMethod

strong

The methods laplace_invert accepts.

Enumerator
Euler
Talbot
GaverStehfest
Cme
Weeks

Definition at line 387 of file laplace_invert.h.

Function Documentation

◆ euler_get_alpha()

std::vector< Cplx > line::lti::euler_get_alpha ( std::size_t n )

inline

Euler nodes: a vertical line at Re = (n-1) log(10) / 6.

Definition at line 83 of file laplace_invert.h.

References euler_get_alpha().

Referenced by euler_get_alpha(), and laplace_invert_euler().

◆ euler_get_eta()

std::vector< double > line::lti::euler_get_eta ( std::size_t n )

inline

Euler weights before the alternating sign and the scale.

The tail is the binomial partial sums that give the Euler acceleration; it is filled BACKWARDS from the last entry, which is what makes the running sum correct.

Definition at line 98 of file laplace_invert.h.

References euler_get_eta(), and line::InputError::InputError().

Referenced by euler_get_eta(), and euler_get_omega().

◆ euler_get_omega()

std::vector< Cplx > line::lti::euler_get_omega ( std::size_t n )

inline

Euler weights: eta, alternating in sign, scaled by 10^((n-1)/6).

Definition at line 111 of file laplace_invert.h.

References euler_get_eta(), and euler_get_omega().

Referenced by euler_get_omega(), and laplace_invert_euler().

◆ gaver_stehfest_get_alpha()

std::vector< double > line::lti::gaver_stehfest_get_alpha ( std::size_t n )

inline

Gaver-Stehfest nodes: k log 2, on the REAL axis.

Definition at line 148 of file laplace_invert.h.

References gaver_stehfest_get_alpha().

Referenced by gaver_stehfest_get_alpha(), and laplace_invert_gaver_stehfest().

◆ gaver_stehfest_get_omega()

std::vector< double > line::lti::gaver_stehfest_get_omega ( std::size_t n )

inline

Gaver-Stehfest weights.

They alternate in sign and grow rapidly with n, which is why the method needs more precision than the others rather than more terms.

Definition at line 161 of file laplace_invert.h.

References gaver_stehfest_get_omega(), and line::InputError::InputError().

Referenced by gaver_stehfest_get_omega(), and laplace_invert_gaver_stehfest().

◆ laplace_invert()

double line::lti::laplace_invert	(	const LaplaceFn &	F,
		double	t,
		LaplaceMethod	method = LaplaceMethod::Euler,
		std::size_t	n = 0 )

inline

Invert F at t by the named method.

Parameters

n	0 takes the method's own default: 41 Euler (see laplace_invert_euler), 32 Talbot, 12 Gaver-Stehfest, 25 CME

Definition at line 407 of file laplace_invert.h.

References Cme, line::mam::Cme, Euler, GaverStehfest, line::InputError::InputError(), laplace_invert(), laplace_invert_euler(), laplace_invert_gaver_stehfest(), laplace_invert_talbot(), laplace_invert_weeks(), line::mam::matlab_ilt(), Talbot, and Weeks.

Referenced by laplace_invert(), laplace_invert_cdf(), laplace_invert_pdf(), and line::qsys::qsys_mgisrgi_whitt().

◆ laplace_invert_cdf()

std::vector< double > line::lti::laplace_invert_cdf	(	const LaplaceFn &	F,
		const std::vector< double > &	t,
		LaplaceMethod	method = LaplaceMethod::Euler,
		std::size_t	n = 0 )

inline

The DISTRIBUTION on a grid, from the transform of the DENSITY.

F(s)/s is the transform of the CDF, so that is what is inverted – passing the CDF's own transform here would invert it twice. The result is clamped into [0,1] and made monotone by a running maximum, because a numerical inversion is pointwise and nothing in it enforces either property; a non-monotone "CDF" then yields negative probabilities downstream.

Definition at line 464 of file laplace_invert.h.

References Euler, laplace_invert(), laplace_invert_cdf(), laplace_invert_weeks(), laplace_weeks_scaling(), and Weeks.

Referenced by line::mc::ctmc_passage_time(), laplace_invert_cdf(), and line::mc::smp_passage_time().

◆ laplace_invert_euler()

double line::lti::laplace_invert_euler	(	const LaplaceFn &	F,
		double	t,
		std::size_t	n = 41 )

inline

Euler inversion of F at t.

n is rounded UP to odd, as the reference does.

THE DEFAULT IS 41, NOT THE OLD 99, AND THAT IS A CORRECTION RATHER THAN A PREFERENCE. The weights carry a factor 10^((n-1)/6), and the sum they multiply alternates in sign, so the method's accuracy is a race between the series converging and the cancellation eating the mantissa. Measured on F(s) = 2/(s+2), whose inverse is 2 exp(-2t), as the worst relative error over t in {0.1, 0.5, 1, 2}:

n = 11 21 31 41 51 71 99 err = 4.4e-3 2.1e-6 1.6e-9 1.6e-10 4.7e-8 1.7e-4 1.4e+0

At the reference's 99 the scale is 2.2e16, past what a double resolves, and the answer is 140 per cent wrong – negative at some t. Native Python has the same default and the same behaviour (measured: 1.997 against an exact 1.637 at t = 0.1, and -0.033 at t = 0.5), and api/lti has no MATLAB twin, so nothing else in the tree catches it. Shipping a default that returns noise is not a convention worth preserving; the reference's own defect is recorded for its maintainer rather than reproduced here.

Definition at line 201 of file laplace_invert.h.

References euler_get_alpha(), euler_get_omega(), line::InputError::InputError(), and laplace_invert_euler().

Referenced by laplace_invert(), and laplace_invert_euler().

◆ laplace_invert_gaver_stehfest()

double line::lti::laplace_invert_gaver_stehfest	(	const RealLaplaceFn &	F,
		double	t,
		std::size_t	n = 12 )

inline

Gaver-Stehfest inversion of F at t.

The transform is sampled on the REAL axis only, which is the whole reason to choose this method. n is rounded DOWN to even.

Definition at line 227 of file laplace_invert.h.

References gaver_stehfest_get_alpha(), gaver_stehfest_get_omega(), line::InputError::InputError(), and laplace_invert_gaver_stehfest().

Referenced by laplace_invert(), and laplace_invert_gaver_stehfest().

◆ laplace_invert_pdf()

std::vector< double > line::lti::laplace_invert_pdf	(	const LaplaceFn &	F,
		const std::vector< double > &	t,
		LaplaceMethod	method = LaplaceMethod::Euler,
		std::size_t	n = 0 )

inline

The DENSITY on a grid: the inversion clamped at zero.

A density cannot be negative, and a numerical inversion can undershoot near the origin or in a tail; the reference clamps, and so does this.

Definition at line 437 of file laplace_invert.h.

References Euler, laplace_invert(), laplace_invert_pdf(), laplace_invert_weeks(), laplace_weeks_scaling(), and Weeks.

Referenced by line::mc::ctmc_passage_time(), laplace_invert_pdf(), and line::mc::smp_passage_time().

◆ laplace_invert_talbot()

double line::lti::laplace_invert_talbot	(	const LaplaceFn &	F,
		double	t,
		std::size_t	n = 32 )

inline

Talbot inversion of F at t.

Definition at line 211 of file laplace_invert.h.

References line::InputError::InputError(), laplace_invert_talbot(), talbot_get_alpha(), and talbot_get_omega().

Referenced by laplace_invert(), and laplace_invert_talbot().

◆ laplace_invert_weeks() [1/2]

double line::lti::laplace_invert_weeks	(	const LaplaceFn &	F,
		double	t,
		std::size_t	p0 = 200 )

inline

Convenience overload: build the parameters, then invert at one point.

Definition at line 382 of file laplace_invert.h.

References laplace_invert_weeks(), and laplace_weeks_scaling().

◆ laplace_invert_weeks() [2/2]

double line::lti::laplace_invert_weeks	(	const WeeksParams &	w,
		double	t )

inline

Invert by the Laguerre series f(t) = sum_n q_n l_n(t), recovered as exp(sigma*b*t) f_{sigma,b}(b*t).

Unlike Euler and Talbot the coefficients do not depend on t, so ONE parameter set serves an arbitrary number of time points: the transform is evaluated 2*p0 times in total, not 2*p0 times per t. That is the property this method is here for, so build the WeeksParams once and reuse it on a grid.

Definition at line 372 of file laplace_invert.h.

References line::lti::WeeksParams::b, laplace_invert_weeks(), line::lti::WeeksParams::q, and line::lti::WeeksParams::sigma.

Referenced by laplace_invert(), laplace_invert_cdf(), laplace_invert_pdf(), laplace_invert_weeks(), and laplace_invert_weeks().

◆ laplace_method()

LaplaceMethod line::lti::laplace_method ( const std::string & s )

inline

Parse the reference's method names, including its two Gaver spellings.

Definition at line 390 of file laplace_invert.h.

References Cme, Euler, GaverStehfest, line::InputError::InputError(), laplace_method(), Talbot, and Weeks.

Referenced by line::mc::ctmc_passage_time(), laplace_method(), line::qsys::qsys_mgisrgi_whitt(), and line::mc::smp_passage_time().

◆ laplace_weeks_coeffs()

std::vector< double > line::lti::laplace_weeks_coeffs	(	const LaplaceFn &	F,
		double	sigma = 0.0,
		double	b = 1.0,
		std::size_t	p0 = 200 )

inline

Laguerre coefficients q_n, n = 0..2*p0-1, of f_{sigma,b}(t) = exp(-sigma t) f(t/b), whose generating function is.

Q_{sigma,b}(z) = b/(1-z) * L( b(1+z)/(2(1-z)) + b*sigma ).

NOTE ON THE PAPER. Eq. 10 as printed carries the factor (1-z) rather than 1/(1-z). The scaled form above, printed later in the same section, carries 1/(1-z) and is the correct one: with l_n(t) = exp(-t/2) L_n(t) the transform of l_n is (s-1/2)^n/(s+1/2)^{n+1}, so L(s) = Q(z)/(s+1/2) with z = (s-1/2)/(s+1/2) and s+1/2 = 1/(1-z). Implementing the printed (1-z) is wrong at every t (163 per cent at t = 0.1 on Exp(2)).

Sec. 4.3 fixes the trapezoid count at 2*p0 and the radius at r = 0.1^(4/p0) for every n, so the quadrature is one discrete Fourier transform of Q sampled on the circle and the transform is evaluated 2*p0 times IN TOTAL rather than per coefficient. The DFT is evaluated directly: at 2*p0 = 400 points that is 160k complex multiplies, which is not worth a dependency.

Definition at line 269 of file laplace_invert.h.

References line::InputError::InputError(), and laplace_weeks_coeffs().

Referenced by laplace_weeks_coeffs(), and laplace_weeks_scaling().

◆ laplace_weeks_scaling()

WeeksParams line::lti::laplace_weeks_scaling	(	const LaplaceFn &	F,
		std::size_t	p0 = 200,
		double	tol = 1e-10 )

inline

The automatic (sigma, b) search of Fig.

1: accept the first pair at which the coefficients have decayed by term p0, doubling sigma from 0.001 and stepping b by 4 whenever sigma passes 0.2.

REFUSES BY NAME when the box is exhausted. Raising b further is counterproductive and excessive damping is unstable in finite precision, and a density with a discontinuity in itself or its derivatives has no usable Laguerre representation at all (Sec. 4.2). Returning the last iterate would report noise as an answer; Euler handles those cases instead.

Definition at line 311 of file laplace_invert.h.

References line::lti::WeeksParams::b, laplace_weeks_coeffs(), laplace_weeks_scaling(), line::NumericError::NumericError(), line::lti::WeeksParams::q, and line::lti::WeeksParams::sigma.

Referenced by laplace_invert_cdf(), laplace_invert_pdf(), laplace_invert_weeks(), and laplace_weeks_scaling().

◆ talbot_get_alpha()

std::vector< Cplx > line::lti::talbot_get_alpha ( std::size_t n )

inline

Talbot nodes: the cotangent contour, bending into the left half plane.

Definition at line 121 of file laplace_invert.h.

References line::InputError::InputError(), and talbot_get_alpha().

Referenced by laplace_invert_talbot(), and talbot_get_alpha().

◆ talbot_get_omega()

std::vector< Cplx > line::lti::talbot_get_omega	(	std::size_t	n,
		const std::vector< Cplx > &	alpha )

inline

Talbot weights, which carry the contour's own derivative.

Definition at line 134 of file laplace_invert.h.

References line::InputError::InputError(), and talbot_get_omega().

Referenced by laplace_invert_talbot(), and talbot_get_omega().

Classes

Typedefs

Enumerations

Functions

Typedef Documentation

◆ Cplx

◆ LaplaceFn

◆ RealLaplaceFn

Enumeration Type Documentation

◆ LaplaceMethod

Function Documentation

◆ euler_get_alpha()

◆ euler_get_eta()

◆ euler_get_omega()

◆ gaver_stehfest_get_alpha()

◆ gaver_stehfest_get_omega()

◆ laplace_invert()

◆ laplace_invert_cdf()

◆ laplace_invert_euler()

◆ laplace_invert_gaver_stehfest()

◆ laplace_invert_pdf()

◆ laplace_invert_talbot()

◆ laplace_invert_weeks() [1/2]

◆ laplace_invert_weeks() [2/2]

◆ laplace_method()

◆ laplace_weeks_coeffs()

◆ laplace_weeks_scaling()

◆ talbot_get_alpha()

◆ talbot_get_omega()