"""Core full-history and Gauss-Jacobi time-memory representations."""
from __future__ import annotations
from collections.abc import Mapping
from dataclasses import dataclass
from math import pi, sin
from types import MappingProxyType
from typing import Any
import numpy as np
from yonderdrake.time.coefficients import FloatArray, gauss_jacobi
[docs]
@dataclass(frozen=True, slots=True)
class FullHistory:
"""Direct time history with piecewise-linear interpolation."""
interpolant: str = "linear"
def __post_init__(self) -> None:
if self.interpolant != "linear":
raise ValueError("only interpolant='linear' is supported")
def describe(self) -> dict[str, str]:
return {
"representation": "FullHistory",
"interpolant": self.interpolant,
}
def _immutable_float_array(values: Any) -> FloatArray:
result = np.array(values, dtype=np.float64, copy=True)
result.setflags(write=False)
return result
@dataclass(frozen=True)
class DiffusiveSpectrum:
"""Normalized positive-rate approximation of a Caputo kernel."""
rates: FloatArray
weights: FloatArray
metadata: MappingProxyType[str, Any]
def __post_init__(self) -> None:
rates = _immutable_float_array(self.rates)
weights = _immutable_float_array(self.weights)
if rates.ndim != 1 or weights.ndim != 1 or rates.shape != weights.shape:
raise ValueError(
"rates and weights must be one-dimensional and equal-sized"
)
if rates.size == 0:
raise ValueError("a spectrum must contain at least one mode")
if not np.all(np.isfinite(rates)) or not np.all(rates > 0.0):
raise ValueError("all rates must be finite and positive")
if not np.all(np.isfinite(weights)) or not np.all(weights > 0.0):
raise ValueError("all weights must be finite and positive")
object.__setattr__(self, "rates", rates)
object.__setattr__(self, "weights", weights)
object.__setattr__(self, "metadata", MappingProxyType(dict(self.metadata)))
class _SingleExponential:
"""Exact one-mode representation of an exponential memory kernel."""
num_modes = 1
@staticmethod
def spectrum(decay_rate: float) -> DiffusiveSpectrum:
try:
rate = float(decay_rate)
except (TypeError, ValueError) as error:
raise TypeError("decay_rate must be a real scalar") from error
if not np.isfinite(rate) or rate <= 0.0:
raise ValueError("decay_rate must be finite and positive")
metadata = MappingProxyType(
{
"representation": "SingleExponential",
"reference": "doi:10.12785/pfda/010201",
"decay_rate": rate,
"num_modes": 1,
"quadrature": "none-exact",
"quadrature_nodes": (rate,),
"ordering": "increasing_rate",
}
)
return DiffusiveSpectrum(
np.asarray([rate], dtype=np.float64),
np.asarray([1.0], dtype=np.float64),
metadata,
)
def describe(self, decay_rate: float | None = None) -> dict[str, Any]:
if decay_rate is None:
return {
"representation": "SingleExponential",
"num_modes": 1,
"reference": "doi:10.12785/pfda/010201",
"status": "exact-exponential-memory",
}
return dict(self.spectrum(decay_rate).metadata)
def _validate_alpha(alpha: float) -> float:
try:
value = float(alpha)
except (TypeError, ValueError) as error:
raise TypeError("alpha must be a real scalar") from error
if not np.isfinite(value) or not 0.0 < value < 1.0:
raise ValueError("alpha must satisfy 0 < alpha < 1")
return value
def validate_checkpoint_representation(
metadata: Any,
representation: Any,
alpha: float,
) -> None:
"""Validate that checkpoint modes use the stepper's exact spectrum."""
if not isinstance(metadata, Mapping):
raise ValueError("checkpoint representation metadata is missing")
expected = representation.describe(alpha)
for field, expected_value in expected.items():
if field == "quadrature_nodes":
continue
if field not in metadata:
raise ValueError(
"checkpoint representation metadata is incomplete"
)
if metadata[field] != expected_value:
raise ValueError(
"checkpoint representation does not match the stepper"
)
try:
observed_nodes = tuple(float(value) for value in metadata["quadrature_nodes"])
except (KeyError, TypeError, ValueError) as error:
raise ValueError(
"checkpoint representation metadata is incomplete"
) from error
if observed_nodes != tuple(expected["quadrature_nodes"]):
raise ValueError("checkpoint representation does not match the stepper")
class _JacobiRepresentation:
_name: str
_reference: str
def __init__(self, num_modes: int, **method_parameters: Any) -> 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 <= 256:
raise ValueError("num_modes must lie between 1 and 256")
unknown = set(method_parameters) - {"rate_scale"}
if unknown:
names = ", ".join(sorted(unknown))
raise TypeError(f"unsupported representation parameter(s): {names}")
rate_scale = float(method_parameters.get("rate_scale", 1.0))
if not np.isfinite(rate_scale) or rate_scale <= 0.0:
raise ValueError("rate_scale must be finite and positive")
self.num_modes = num_modes
self._rate_scale = rate_scale
def _jacobi_exponents(self, alpha: float) -> tuple[float, float]:
raise NotImplementedError
def _native_coefficients(
self,
alpha: float,
nodes: FloatArray,
quadrature_weights: FloatArray,
) -> tuple[FloatArray, FloatArray]:
raise NotImplementedError
def spectrum(self, alpha: float) -> DiffusiveSpectrum:
"""Generate a normalized spectrum for one immutable Caputo order."""
order = _validate_alpha(alpha)
jacobi_alpha, jacobi_beta = self._jacobi_exponents(order)
nodes, quadrature_weights = gauss_jacobi(
self.num_modes,
jacobi_alpha,
jacobi_beta,
)
rates, weights = self._native_coefficients(
order,
nodes,
quadrature_weights,
)
rates = self._rate_scale * rates
weights = self._rate_scale**order * weights
permutation = np.argsort(rates, kind="stable")
rates = rates[permutation]
weights = weights[permutation]
ordered_nodes = nodes[permutation]
metadata = MappingProxyType(
{
"representation": self._name,
"reference": self._reference,
"alpha": order,
"num_modes": self.num_modes,
"rate_scale": self._rate_scale,
"quadrature": "Gauss-Jacobi",
"jacobi_alpha": jacobi_alpha,
"jacobi_beta": jacobi_beta,
"quadrature_nodes": tuple(float(node) for node in ordered_nodes),
"ordering": "increasing_rate",
}
)
return DiffusiveSpectrum(rates, weights, metadata)
def describe(self, alpha: float | None = None) -> dict[str, Any]:
"""Describe configuration, or full generated metadata when ordered."""
if alpha is None:
return {
"representation": self._name,
"num_modes": self.num_modes,
"rate_scale": self._rate_scale,
"reference": self._reference,
"configurable_parameters": ("rate_scale",),
}
return dict(self.spectrum(alpha).metadata)
[docs]
class Diethelm(_JacobiRepresentation):
"""Diethelm's Gauss-Jacobi improvement of the diffusive representation."""
_name = "Diethelm"
_reference = "doi:10.1007/s11075-008-9193-8"
def _jacobi_exponents(self, alpha: float) -> tuple[float, float]:
transformed_order = 2.0 * alpha - 1.0
return transformed_order, -transformed_order
def _native_coefficients(
self,
alpha: float,
nodes: FloatArray,
quadrature_weights: FloatArray,
) -> tuple[FloatArray, FloatArray]:
denominator = 1.0 + nodes
ratio = (1.0 - nodes) / denominator
rates = np.square(ratio)
weights = (
(4.0 * sin(pi * alpha) / pi)
* quadrature_weights
/ np.square(denominator)
)
return rates, weights
[docs]
class BirkSong(_JacobiRepresentation):
"""Birk-Song's squared Cayley-transform Gauss-Jacobi spectrum."""
_name = "BirkSong"
_reference = "doi:10.1007/s00466-010-0510-4"
def _jacobi_exponents(self, alpha: float) -> tuple[float, float]:
transformed_order = 2.0 * alpha - 1.0
return 2.0 * transformed_order + 1.0, 1.0 - 2.0 * transformed_order
def _native_coefficients(
self,
alpha: float,
nodes: FloatArray,
quadrature_weights: FloatArray,
) -> tuple[FloatArray, FloatArray]:
denominator = 1.0 + nodes
ratio = (1.0 - nodes) / denominator
rates = np.power(ratio, 4)
weights = (
(8.0 * sin(pi * alpha) / pi)
* quadrature_weights
/ np.power(denominator, 4)
)
return rates, weights