"""Log-rate and Gauss-Laguerre time-memory representations."""
from __future__ import annotations
from functools import cache
from math import isfinite, log, pi, sin, sqrt
from types import MappingProxyType
from typing import Any
import numpy as np
from scipy.linalg import eigh_tridiagonal
from scipy.optimize import brentq, minimize_scalar
from scipy.special import roots_legendre
from yonderdrake.time.representations.core import (
DiffusiveSpectrum,
_validate_alpha,
)
_MAX_NUM_MODES = 16_385
_MIN_LOG_RATE = log(float(np.nextafter(0.0, 1.0)))
_MAX_LOG_RATE = log(float(np.finfo(np.float64).max))
def _positive_finite(value: Any, name: str) -> float:
try:
result = float(value)
except (TypeError, ValueError) as error:
raise TypeError(f"{name} must be a real scalar") from error
if not np.isfinite(result) or result <= 0.0:
raise ValueError(f"{name} must be finite and positive")
return result
@cache
def _laguerre_nodes_and_unweighted_weights(
num_modes: int,
) -> tuple[np.ndarray, np.ndarray]:
indices = np.arange(num_modes, dtype=np.float64)
nodes = eigh_tridiagonal(
2.0 * indices + 1.0,
indices[1:],
eigvals_only=True,
)
previous = np.ones(num_modes, dtype=np.float64)
current = 1.0 - nodes
log_scale = np.zeros(num_modes, dtype=np.float64)
scale = np.maximum(np.abs(previous), np.abs(current))
previous /= scale
current /= scale
log_scale += np.log(scale)
for degree in range(1, num_modes + 1):
following = (
(2.0 * degree + 1.0 - nodes) * current
- degree * previous
) / (degree + 1.0)
scale = np.maximum(np.abs(current), np.abs(following))
previous = current / scale
current = following / scale
log_scale += np.log(scale)
log_effective_weights = (
np.log(nodes)
- 2.0 * log(num_modes + 1.0)
- 2.0 * (np.log(np.abs(current)) + log_scale)
+ nodes
)
effective_weights = np.exp(log_effective_weights)
nodes.setflags(write=False)
effective_weights.setflags(write=False)
return nodes, effective_weights
[docs]
class YuanAgrawal:
"""Original Yuan-Agrawal Gauss-Laguerre diffusive spectrum.
Warning:
This representation is provided for expert use and comparison. It is
generally less robust per mode than :class:`BirkSong` or
:class:`Diethelm`.
"""
_name = "YuanAgrawal"
_reference = "doi:10.1115/1.1448322"
def __init__(self, num_modes: int, *, rate_scale: float = 1.0) -> None:
if isinstance(num_modes, bool) or not isinstance(num_modes, int):
raise TypeError("num_modes must be an integer")
if not 1 <= num_modes <= 2048:
raise ValueError("num_modes must lie between 1 and 2048")
self.num_modes = num_modes
self._rate_scale = _positive_finite(
rate_scale,
"rate_scale",
)
[docs]
def spectrum(self, alpha: float) -> DiffusiveSpectrum:
"""Generate the original squared-node Gauss-Laguerre spectrum."""
order = _validate_alpha(alpha)
nodes, effective_weights = _laguerre_nodes_and_unweighted_weights(
self.num_modes
)
rates = self._rate_scale * np.square(nodes)
weights = (
(2.0 * sin(pi * order) / pi)
* effective_weights
* np.power(nodes, 2.0 * order - 1.0)
* self._rate_scale**order
)
metadata = MappingProxyType(
{
"representation": self._name,
"reference": self._reference,
"alpha": order,
"num_modes": self.num_modes,
"rate_scale": self._rate_scale,
"quadrature": "Gauss-Laguerre",
"quadrature_nodes": tuple(float(node) for node in nodes),
"rate_map": "squared",
"ordering": "increasing_rate",
}
)
return DiffusiveSpectrum(rates, weights, metadata)
[docs]
def describe(self, alpha: float | None = None) -> dict[str, Any]:
"""Describe the configuration or one generated spectrum."""
if alpha is not None:
return dict(self.spectrum(alpha).metadata)
return {
"representation": self._name,
"num_modes": self.num_modes,
"rate_scale": self._rate_scale,
"reference": self._reference,
"configurable_parameters": ("rate_scale",),
"status": "supported-not-recommended",
}
[docs]
class Diethelm2022:
"""Truncate Diethelm's 2022 real-line representation for quadrature.
``quadrature`` selects composite trapezoidal, composite Simpson, or
Gauss-Legendre weights. Simpson quadrature requires an odd ``num_modes``.
Automatic truncation balances tail and grid estimates. An explicit
``truncation_radius`` bypasses that balance.
Warning:
This representation is provided for expert use and comparison. It is
more sensitive to mode count, time range, and scaling than BirkSong or
Diethelm.
"""
_name = "Diethelm2022"
_reference = "doi:10.1007/978-981-19-7716-9_1"
def __init__(
self,
num_modes: int,
*,
quadrature: str = "trapezoidal",
target_error: float = 1.0e-8,
decay_scale: float = 1.0,
truncation_radius: float | None = None,
rate_scale: float = 1.0,
) -> None:
if isinstance(num_modes, bool) or not isinstance(num_modes, int):
raise TypeError("num_modes must be an integer")
if not 3 <= num_modes <= _MAX_NUM_MODES:
raise ValueError(
f"num_modes must lie between 3 and {_MAX_NUM_MODES}"
)
if quadrature not in {
"trapezoidal",
"simpson",
"gauss-legendre",
}:
raise ValueError(
"quadrature must be 'trapezoidal', 'simpson', "
"or 'gauss-legendre'"
)
if quadrature == "simpson" and num_modes % 2 == 0:
raise ValueError(
"Simpson quadrature requires an odd num_modes"
)
self.num_modes = num_modes
self._quadrature = quadrature
self._target_error = _positive_finite(
target_error,
"target_error",
)
self._decay_scale = _positive_finite(
decay_scale,
"decay_scale",
)
self._rate_scale = _positive_finite(
rate_scale,
"rate_scale",
)
if truncation_radius is None:
self._truncation_radius = None
else:
self._truncation_radius = _positive_finite(
truncation_radius,
"truncation_radius",
)
def _automatic_radius(self, alpha: float) -> float:
slow_decay = min(alpha, 1.0 - alpha)
initial_radius = pi * sqrt(
(self.num_modes - 1) / (2.0 * slow_decay)
)
if not isfinite(initial_radius):
raise ValueError(
"balanced truncation radius is nonfinite for this alpha"
)
log_rate_scale = log(self._rate_scale)
float64_radius = min(
log_rate_scale - _MIN_LOG_RATE,
_MAX_LOG_RATE - log_rate_scale,
)
upper_radius = min(
float64_radius,
max(8.0, 4.0 * initial_radius),
)
if upper_radius <= 1.0:
raise ValueError(
"rate_scale leaves no finite interval for automatic truncation"
)
result = minimize_scalar(
lambda radius: self._estimated_total_bound(alpha, radius),
bounds=(1.0, upper_radius),
method="bounded",
options={"xatol": 1.0e-10},
)
if not result.success or not isfinite(float(result.x)):
raise ValueError("could not determine a balanced truncation radius")
balanced_radius = float(result.x)
if (
self._estimated_total_bound(alpha, balanced_radius)
> self._target_error
):
return balanced_radius
if self._estimated_total_bound(alpha, 1.0) <= self._target_error:
return 1.0
return float(
brentq(
lambda radius: (
self._estimated_total_bound(alpha, radius)
- self._target_error
),
1.0,
balanced_radius,
xtol=1.0e-12,
)
)
def _tail_bound(self, alpha: float, radius: float) -> float:
coefficient = (
self._decay_scale * abs(sin(pi * alpha)) / pi
)
return float(
coefficient
* (
np.exp(-alpha * radius) / alpha
+ np.exp(-(1.0 - alpha) * radius) / (1.0 - alpha)
)
)
def _discretization_bound(self, radius: float) -> float:
spacing = 2.0 * radius / (self.num_modes - 1)
trapezoidal_bound = np.exp(-(pi * pi) / spacing)
if self._quadrature == "trapezoidal":
return float(trapezoidal_bound)
if self._quadrature == "gauss-legendre":
# Radius-selection surrogate; unlike the trapezoidal estimate,
# this is not treated as an accuracy certificate.
return float(trapezoidal_bound)
# Composite Simpson is (4 T_h - T_{2h}) / 3. Bound both terms;
# the embedded coarser trapezoid controls the asymptotic exponent.
coarse_trapezoidal_bound = np.exp(
-(pi * pi) / (2.0 * spacing)
)
return float(
(
4.0 * trapezoidal_bound
+ coarse_trapezoidal_bound
)
/ 3.0
)
def _estimated_total_bound(self, alpha: float, radius: float) -> float:
return self._tail_bound(alpha, radius) + self._discretization_bound(
radius
)
[docs]
def spectrum(self, alpha: float) -> DiffusiveSpectrum:
"""Generate positive modes for the truncated real-line integral."""
order = _validate_alpha(alpha)
radius = (
self._automatic_radius(order)
if self._truncation_radius is None
else self._truncation_radius
)
log_rate_scale = log(self._rate_scale)
if (
log_rate_scale - radius < _MIN_LOG_RATE
or log_rate_scale + radius > _MAX_LOG_RATE
):
raise ValueError(
"truncation interval produces nonfinite float64 rates; "
"relax target_error or provide a smaller "
"truncation_radius"
)
nominal_spacing = 2.0 * radius / (self.num_modes - 1)
if self._quadrature == "gauss-legendre":
reference_nodes, quadrature_factors = roots_legendre(
self.num_modes
)
nodes = radius * reference_nodes
quadrature_scale = radius
quadrature_name = "Gauss-Legendre"
log_rate_spacing = None
maximum_log_rate_spacing = float(np.max(np.diff(nodes)))
else:
nodes = np.linspace(
-radius,
radius,
self.num_modes,
dtype=np.float64,
)
quadrature_factors = np.ones(
self.num_modes,
dtype=np.float64,
)
log_rate_spacing = nominal_spacing
maximum_log_rate_spacing = nominal_spacing
if self._quadrature == "trapezoidal":
quadrature_factors[[0, -1]] = 0.5
quadrature_scale = nominal_spacing
quadrature_name = "composite-trapezoidal"
else:
quadrature_factors[1:-1:2] = 4.0
quadrature_factors[2:-1:2] = 2.0
quadrature_scale = nominal_spacing / 3.0
quadrature_name = "composite-simpson"
estimated_tail_bound = self._tail_bound(order, radius)
estimated_discretization_bound = self._discretization_bound(radius)
estimated_total_bound = (
estimated_tail_bound + estimated_discretization_bound
)
target_achievable = bool(
estimated_total_bound
<= self._target_error
* (1.0 + 64.0 * float(np.finfo(np.float64).eps))
)
rates = np.exp(log_rate_scale + nodes)
weights = (
quadrature_scale
* quadrature_factors
* (sin(pi * order) / pi)
* np.exp(order * (log_rate_scale + nodes))
)
metadata = MappingProxyType(
{
"representation": self._name,
"reference": self._reference,
"alpha": order,
"num_modes": self.num_modes,
"rate_scale": self._rate_scale,
"quadrature": quadrature_name,
"quadrature_nodes": tuple(float(node) for node in nodes),
"ordering": "increasing_rate",
"target_error": self._target_error,
"decay_scale": self._decay_scale,
"truncation_radius": radius,
"truncation_radius_source": (
(
"target_error"
if target_achievable
else "balanced_tail_and_grid"
)
if self._truncation_radius is None
else "user"
),
"estimated_tail_bound": estimated_tail_bound,
"estimated_discretization_bound": (
estimated_discretization_bound
),
"estimated_total_bound": estimated_total_bound,
"target_achievable": (
None
if self._quadrature == "gauss-legendre"
else target_achievable
),
"discretization_error_model": (
"analytic_strip_radius_surrogate"
if self._quadrature == "gauss-legendre"
else "analytic_strip_bound"
),
"log_rate_spacing": log_rate_spacing,
"maximum_log_rate_spacing": maximum_log_rate_spacing,
}
)
return DiffusiveSpectrum(rates, weights, metadata)
[docs]
def describe(self, alpha: float | None = None) -> dict[str, Any]:
"""Describe the configuration or one ordered spectrum."""
if alpha is not None:
return dict(self.spectrum(alpha).metadata)
return {
"representation": self._name,
"num_modes": self.num_modes,
"rate_scale": self._rate_scale,
"quadrature": self._quadrature,
"target_error": self._target_error,
"decay_scale": self._decay_scale,
"truncation_radius": self._truncation_radius,
"reference": self._reference,
"configurable_parameters": (
"quadrature",
"target_error",
"decay_scale",
"truncation_radius",
"rate_scale",
),
"status": "supported-not-recommended",
}