Ray (WKB) asymptotic expansion of the list-based cache normalizing constant. More...

#include <cmath>
#include <cstddef>
#include <functional>
#include <string>
#include <vector>
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for retrieval_rayint.h:

Classes
struct	line::retrieval::RetrievalRayintResult< T >
	Outcome of the expansion. More...

Namespaces
namespace	line
namespace	line::retrieval

Functions
template<class T>
RetrievalRayintResult< T >	line::retrieval::retrieval_rayint (const Matrix< T > &gamma, const std::vector< int > &m)
	Discrete (saddle) form.
template<class T>
RetrievalRayintResult< T >	line::retrieval::retrieval_rayint (const std::function< Matrix< T >(const std::vector< T > &)> &gfun, const std::vector< int > &m, int n, std::size_t nquad=4097)
	Continuum (ray-integral) form.

Detailed Description

Ray (WKB) asymptotic expansion of the list-based cache normalizing constant.

Templated port of matlab/src/api/retrieval/retrieval_rayint.m, cross-checked against jar/src/main/java/jline/api/retrieval/Retrieval_rayint.java.

Approximates the constant cache_erec computes exactly, in the SAME normalization, so the two are interchangeable:

E(m,n) = E(m,n-1) + sum_j gamma_{n,j} m_j E(m-1_j,n-1), E(0,0)=1

Writing E = prod_j m_j! * Et, the relaxation Et ~ H exp(phi/eps) with n = y/eps and m_j = x_j/eps gives the eikonal e^{phi_y} = 1 + sum_j gamma_j(y) e^{-phi_j}, whose rays carry the constants xi_j = e^{-phi_j}. With S(v) = 1 + sum_j gamma_j(v) xi_j,

x_j = int_0^y gamma_j(v) xi_j / S(v) dv (the saddle conditions) phi = int_0^y log S(v) dv - sum_j x_j log xi_j H = (2 pi)^{-h/2} sqrt(S(y)/S(0)) / sqrt(prod_j xi_j * det A) A_{ik} = d x_i / d xi_k

and E ~ prod_j m_j! * eps^{h/2} H exp(phi/eps).

DISCRETE (gamma an n x h matrix). The ray integrals are the sums they discretize and the expansion collapses to the Laplace form

Et ~ (2 pi)^{-h/2} exp(sum_k log D_k - sum_j m_j log xi_j) / sqrt(det Sigma)

with D_k = 1 + sum_j gamma_{k,j} xi_j, sum_k gamma_{k,j} xi_j / D_k = m_j and Sigma = A * diag(xi) the Hessian in log xi. This is the more accurate of the two forms; the sqrt(S(y)/S(0)) factor is exactly the Euler-Maclaurin term relating sum_k to int dv and is already accounted for.

CONTINUUM (gamma a callable profile on v in [0,1]). Composite Simpson quadrature on the profile itself, the form written in the note. Costs roughly a factor two in accuracy but does not need the n rows.

ACCURACY. The relative error is O(1/n) at fixed occupancy but is governed by the smallest occupancy rather than by n, tracking 0.14 * (1/min_j m_j + 1/(n - sum_j m_j)), so a per cent needs every m_j and n - sum_j m_j above about 15 and a part in a thousand needs them above about

Returned in the result's relerr_est. Lists with m_j = 0 contribute nothing and are dropped before the saddle is solved.

ARITHMETIC: the expansion is a Laplace approximation built out of logs, exps and a square root, so it is meaningless at exact arithmetic and is gated on has_transcendental. It is also an APPROXIMATION whatever the arithmetic – widening the type sharpens the saddle solve, never the O(1/n) model error. Use cache_erec or retrieval_nc when the exact constant is wanted.

This is the no-fetch (q=0) case, i.e. the same quantity as cache_erec. The delayed-hit extension carrying the fetch coordinates is NOT implemented: its eikonal is known but its amplitude has not been derived.

Definition in file retrieval_rayint.h.

Classes

Namespaces

Functions

Detailed Description