From 8ad67479040e5494adf47ea20874f2e6515f04e0 Mon Sep 17 00:00:00 2001 From: Jonas Wessely Date: Tue, 26 May 2026 16:49:11 +0200 Subject: [PATCH 1/8] new integrators, logscaled and Gauss-legendre. Minor bugfux in kernely.py regarding the cache, added analytic parrameter gradients for kernels with asymptotics --- fredipy/integrators.py | 272 +++++++++++++++++++++++++++++++++++++++++ fredipy/kernels.py | 136 +++++++++++++++++++++ 2 files changed, 408 insertions(+) diff --git a/fredipy/integrators.py b/fredipy/integrators.py index 8251736..3d44456 100644 --- a/fredipy/integrators.py +++ b/fredipy/integrators.py @@ -190,6 +190,278 @@ def singleIntegration( ) +class Riemann_1D_log(Integrator): + """Riemann integration in 1D on a logarithmic (geomspace) grid. + + Uses midpoints equally spaced in log-space and a scalar step + ``dw = d(log ω) = log(w_max/w_min) / int_n``. + + Because the Jacobian of the substitution t = log ω is ω, every kernel + passed to this integrator must include that factor, i.e. the caller is + responsible for passing ``ω · K(p, ω)`` instead of ``K(p, ω)``. + """ + + def __init__( + self, + w_min: float = 0.01, + w_max: float = 10., + int_n: int = 1000 + ) -> None: + + w_edges = make_column_vector(np.geomspace(w_min, w_max, int_n + 1)) + self.dw = np.log(w_max / w_min) / int_n # uniform step in log-space + self.w = np.sqrt(w_edges[:-1] * w_edges[1:]) # geometric midpoints + + def doubleIntegrationSymmetric( + self, + constraint: LinearEquality, + kernel: Callable + ) -> np.ndarray: + return self.doubleIntegration(constraint, kernel, constraint) + + def doubleIntegration( + self, + constraint1: LinearEquality, + kernel: Callable, + constraint2: LinearEquality + ) -> np.ndarray: + return self.dw**2 * ( + constraint1(make_row_vector(self.w), x=make_column_vector(constraint1.x)) + @ kernel(self.w, self.w) + @ (constraint2(make_row_vector(self.w), x=make_column_vector(constraint2.x))).T + ) + + def singleIntegration( + self, + constraint: LinearEquality, + kernel: Callable, + w_pred: np.ndarray + ) -> np.ndarray: + return self.dw * ( + constraint(make_row_vector(self.w), x=make_column_vector(constraint.x)) + @ kernel(self.w, w_pred) + ) + + +class GaussLegendre_1D_log(Integrator): + """Gauss-Legendre quadrature in 1D on a logarithmic grid. + + Maps the standard Gauss-Legendre nodes from [-1, 1] to log-space + [log(w_min), log(w_max)] and then to ω-space via exponentiation. + + Achieves exponential convergence in the number of quadrature points for + smooth integrands, requiring far fewer points than Riemann rules for the + same accuracy (typically 50–200 instead of 1000). + + As with ``Riemann_1D_log``, every kernel passed to this integrator must + include the log-space Jacobian factor ω, i.e. the caller is responsible + for passing ``ω · K(p, ω)`` instead of ``K(p, ω)``. + """ + + def __init__( + self, + w_min: float = 0.01, + w_max: float = 10., + int_n: int = 100 + ) -> None: + + # GL nodes (xi in [-1,1]) and weights + xi, wi = np.polynomial.legendre.leggauss(int_n) + + # Map nodes from [-1,1] to log-space [log(w_min), log(w_max)] + log_min = np.log(w_min) + log_max = np.log(w_max) + half_range = 0.5 * (log_max - log_min) + + t = half_range * xi + 0.5 * (log_max + log_min) # nodes in log-space + self.w = make_column_vector(np.exp(t)) # nodes in ω-space + # Weights absorb the half-range Jacobian of the linear change of variables + self.weights = make_row_vector(wi * half_range) # shape (1, int_n) + + def doubleIntegrationSymmetric( + self, + constraint: LinearEquality, + kernel: Callable + ) -> np.ndarray: + return self.doubleIntegration(constraint, kernel, constraint) + + def doubleIntegration( + self, + constraint1: LinearEquality, + kernel: Callable, + constraint2: LinearEquality + ) -> np.ndarray: + return ( + self.weights * constraint1(make_row_vector(self.w), x=make_column_vector(constraint1.x)) + @ kernel(self.w, self.w) + @ (self.weights * constraint2(make_row_vector(self.w), x=make_column_vector(constraint2.x))).T + ) + + def singleIntegration( + self, + constraint: LinearEquality, + kernel: Callable, + w_pred: np.ndarray + ) -> np.ndarray: + return ( + self.weights * constraint(make_row_vector(self.w), x=make_column_vector(constraint.x)) + @ kernel(self.w, w_pred) + ) + + +class GaussLegendre_1D_log_UVtail(GaussLegendre_1D_log): + """GL log-space quadrature on (w_min, w_uv) with an analytic UV tail. + + For ω > w_uv the AsymptoticKernel factorises as + K(ω, ω') ≈ f_uv(ω) f_uv(ω'), + so the analytic tail T = ∫_{log w_uv}^∞ (2t)^{-35/22} dt can be + folded in exactly via corrections to the integration methods. + + The substitution u = (2 log ω)^{-13/22} maps the UV integrand to a + constant, giving the closed form + + T = (2 log w_uv)^{-13/22} / (2 × 13/22). + + **Correction formulas** (A_θ = UV-component of the bulk integral): + + * ``doubleIntegrationSymmetric``: + full = bulk + 2 A_θ T + T² + where A_θ[m] ≈ (1/f_uv(w_uv)) × Σ_i W_i C_m(ω_i) K(ω_i, w_uv) + + * ``singleIntegration``: + full = numerical + T × K(w_uv, ω_pred) / f_uv(w_uv) + (K(w_uv, ω_pred)/f_uv(w_uv) ≈ θ_uv(ω_pred) f_uv(ω_pred) for + prediction points far from w_uv, where the RBF part of K vanishes) + + Parameters + ---------- + w_min : float + w_uv : float + UV split where UV asymptotics fully apply. Recommended: ≥ 1000. + int_n : int + uv_func : callable + The UV asymptotic function f_uv(ω), e.g. ``uv_asymptotics``. + """ + + def __init__( + self, + w_min: float, + w_uv: float, + int_n: int, + uv_func: Callable + ) -> None: + super().__init__(w_min, w_uv, int_n) + self.uv_func = uv_func + # Anchor point at the UV split boundary + self.w_uv = make_column_vector(np.array([w_uv])) # (1, 1) + self.f_uv_anchor = float(uv_func(self.w_uv).flat[0]) # scalar: f_uv(w_uv) + # Analytic tail: ∫_{log w_uv}^∞ (2t)^{-35/22} dt + self.T_uv = (2.0 * np.log(w_uv)) ** (-13.0 / 22.0) / (2.0 * 13.0 / 22.0) + + def doubleIntegrationSymmetric( + self, + constraint: LinearEquality, + kernel: Callable + ) -> np.ndarray: + bulk = super().doubleIntegrationSymmetric(constraint, kernel) # (M, M) + + # A_θ ≈ (W C @ K_col) / f_uv(w_uv), shape (M, 1) + # K_col[i] = K(ω_i, w_uv) ≈ θ_uv(ω_i) f_uv(ω_i) f_uv(w_uv) (for large w_uv) + K_col = kernel(self.w, self.w_uv) # (N, 1) + c_row = constraint(make_row_vector(self.w), x=make_column_vector(constraint.x)) + A = (self.weights * c_row @ K_col) / self.f_uv_anchor # (M, 1) + + # correction_{mn} = T (A_m + A_n) + T² + ones_col = np.ones((A.shape[0], 1)) + correction = self.T_uv * (A @ ones_col.T + ones_col @ A.T) + self.T_uv ** 2 + return bulk + correction + + def singleIntegration( + self, + constraint: LinearEquality, + kernel: Callable, + w_pred: np.ndarray + ) -> np.ndarray: + numerical = super().singleIntegration(constraint, kernel, w_pred) + # tail = T × K(w_uv, ω_pred) / f_uv(w_uv), shape (1, N_pred) + # K(w_uv, ω_pred) / f_uv(w_uv) ≈ θ_uv(ω_pred) f_uv(ω_pred) for w_uv >> ω_pred + tail = (self.T_uv / self.f_uv_anchor) * kernel(self.w_uv, w_pred) # (1, N_pred) + return numerical + tail + + +class GaussLegendre_1D_semiinf(Integrator): + """Gauss-Legendre quadrature on the full semi-infinite interval (0, ∞). + + Uses the rational change of variables + + ω = w_scale · (1 + u) / (1 − u), u ∈ (−1, 1) + + which maps (−1, 1) exactly onto (0, ∞). The log-space Jacobian + d(log ω)/du = 2/(1 − u²) is absorbed into the quadrature weights, + so every kernel passed here must include the log-space factor ω, + i.e. pass ``ω · K(p, ω)`` instead of ``K(p, ω)``. + + Parameters + ---------- + w_scale : float + The scale point: u = 0 maps to ω = w_scale. Choose it near the + centre of the integrand in log-space for the best node distribution. + int_n : int + Number of Gauss-Legendre nodes. + + Notes + ----- + For n = 100 nodes and w_scale = 1, the outermost nodes lie at + roughly ω ≈ 7×10⁻⁵ and ω ≈ 1.4×10⁴. The GL weights at those + extreme nodes automatically vanish proportionally to (1 − u²), + exactly cancelling the Jacobian divergence, so the product weight + W_i = w_i · 2/(1 − u_i²) is bounded O(1/n²) at all nodes. + """ + + def __init__( + self, + w_scale: float = 1.0, + int_n: int = 100 + ) -> None: + + xi, wi = np.polynomial.legendre.leggauss(int_n) + + # rational map (-1,1) → (0,∞) + self.w = make_column_vector(w_scale * (1.0 + xi) / (1.0 - xi)) + # log-space weights: dlog(ω)/du = 2/(1-u²) + self.weights = make_row_vector(wi * 2.0 / (1.0 - xi**2)) + + def doubleIntegrationSymmetric( + self, + constraint: LinearEquality, + kernel: Callable + ) -> np.ndarray: + return self.doubleIntegration(constraint, kernel, constraint) + + def doubleIntegration( + self, + constraint1: LinearEquality, + kernel: Callable, + constraint2: LinearEquality + ) -> np.ndarray: + return ( + self.weights * constraint1(make_row_vector(self.w), x=make_column_vector(constraint1.x)) + @ kernel(self.w, self.w) + @ (self.weights * constraint2(make_row_vector(self.w), x=make_column_vector(constraint2.x))).T + ) + + def singleIntegration( + self, + constraint: LinearEquality, + kernel: Callable, + w_pred: np.ndarray + ) -> np.ndarray: + return ( + self.weights * constraint(make_row_vector(self.w), x=make_column_vector(constraint.x)) + @ kernel(self.w, w_pred) + ) + + class Simpson_1D(Integrator): """Implementation of Simpson's rule for 1D integration""" diff --git a/fredipy/kernels.py b/fredipy/kernels.py index 3a2c87c..638fc98 100644 --- a/fredipy/kernels.py +++ b/fredipy/kernels.py @@ -278,6 +278,12 @@ def __init__( self._K_asymp = None self.dim = kernel.dim + def _empty_cache(self) -> None: + """Clear this kernel's cache and propagate to the inner kernel.""" + self._K_asymp = None + self._x, self._y = np.array([None]), np.array([None]) + self.kernel._empty_cache() + def add_asymptotics( self, region: str, @@ -370,6 +376,136 @@ def set_params( self.mu_uv = asymp_params[0] self.l_uv = asymp_params[1] + def params_gradient(self) -> List[Callable]: + """Analytic gradient of the asymptotic kernel w.r.t. all hyperparameters. + + Returns a list of functions ``f_i(x, y)`` such that + ``dK_asymp / dtheta_i = f_i(x, y)``. + + Parameter ordering (identical to ``set_params`` / flat-vector convention): + [base_kernel_params..., mu_uv, l_uv (if UV), mu_ir, l_ir (if IR)] + + Derivation + ---------- + With the shorthands + σ_+(x) = softtheta(x, mu, l, +1) (sigmoid, 0 → 1) + σ_-(x) = softtheta(x, mu, l, −1) (complement, 1 → 0) + and the three kernel regions + A = σ_-(x;ir)·σ_-(y;ir)·f_ir(x)·f_ir(y) + B = σ_+(x;ir)·σ_+(y;ir)·K_rbf·σ_-(x;uv)·σ_-(y;uv) + C = σ_+(x;uv)·σ_+(y;uv)·f_uv(x)·f_uv(y) + the sigmoid derivative identity dσ_±/dmu = ∓(1/l)·σ_+·σ_− gives: + + dK/dmu_uv = (1/l_uv)·[B·(σ_+(x;uv)+σ_+(y;uv)) − C·(σ_-(x;uv)+σ_-(y;uv))] + dK/dl_uv = (1/l_uv²)·[B·((x−mu_uv)·σ_+(x;uv)+(y−mu_uv)·σ_+(y;uv)) + −C·((x−mu_uv)·σ_-(x;uv)+(y−mu_uv)·σ_-(y;uv))] + dK/dmu_ir = (1/l_ir)·[A·(σ_+(x;ir)+σ_+(y;ir)) − B·(σ_-(x;ir)+σ_-(y;ir))] + dK/dl_ir = (1/l_ir²)·[A·((x−mu_ir)·σ_+(x;ir)+(y−mu_ir)·σ_+(y;ir)) + −B·((x−mu_ir)·σ_-(x;ir)+(y−mu_ir)·σ_-(y;ir))] + + The base-kernel gradients are windowed by the RBF transition region: + dK/d(rbf_i) = σ_+(x;ir)·σ_+(y;ir)·(dK_rbf/dtheta_i)·σ_-(x;uv)·σ_-(y;uv) + + When only one asymptotic region is active the absent softthetas reduce to + 1 (sign=0 path in softtheta), so the formulas remain valid throughout. + """ + grads = [] + + # ---- base-kernel parameter gradients -------------------------------- + base_grads = self.kernel.params_gradient() + for bg in base_grads: + def _wrap(x, y, _bg=bg): + x1 = make_column_vector(x[:, 0]) + y1 = make_row_vector(y[:, 0]) + tau_x = softtheta(x1, self.mu_ir, self.l_ir, -self.ir) # σ_+(x;ir) + tau_y = softtheta(y1, self.mu_ir, self.l_ir, -self.ir) + sig_x = softtheta(x1, self.mu_uv, self.l_uv, -self.uv) # σ_-(x;uv) + sig_y = softtheta(y1, self.mu_uv, self.l_uv, -self.uv) + return tau_x * tau_y * _bg(x, y) * sig_x * sig_y + grads.append(_wrap) + + # ---- UV asymptotic parameter gradients: mu_uv, l_uv ---------------- + if self.uv: + def _dK_dmu_uv(x, y): + x1 = make_column_vector(x[:, 0]) + y1 = make_row_vector(y[:, 0]) + tau_x = softtheta(x1, self.mu_ir, self.l_ir, -self.ir) + tau_y = softtheta(y1, self.mu_ir, self.l_ir, -self.ir) + sp_x = softtheta(x1, self.mu_uv, self.l_uv, self.uv) # σ_+(x;uv) + sp_y = softtheta(y1, self.mu_uv, self.l_uv, self.uv) + sm_x = softtheta(x1, self.mu_uv, self.l_uv, -self.uv) # σ_-(x;uv) + sm_y = softtheta(y1, self.mu_uv, self.l_uv, -self.uv) + K_rbf = self.kernel(x, y) + f_x = self.uv_asymptotics(x1) + f_y = self.uv_asymptotics(y1) + return (1.0 / self.l_uv) * ( + tau_x * tau_y * K_rbf * sm_x * sm_y * (sp_x + sp_y) + - sp_x * sp_y * f_x * f_y * (sm_x + sm_y) + ) + + def _dK_dl_uv(x, y): + x1 = make_column_vector(x[:, 0]) + y1 = make_row_vector(y[:, 0]) + tau_x = softtheta(x1, self.mu_ir, self.l_ir, -self.ir) + tau_y = softtheta(y1, self.mu_ir, self.l_ir, -self.ir) + sp_x = softtheta(x1, self.mu_uv, self.l_uv, self.uv) + sp_y = softtheta(y1, self.mu_uv, self.l_uv, self.uv) + sm_x = softtheta(x1, self.mu_uv, self.l_uv, -self.uv) + sm_y = softtheta(y1, self.mu_uv, self.l_uv, -self.uv) + K_rbf = self.kernel(x, y) + f_x = self.uv_asymptotics(x1) + f_y = self.uv_asymptotics(y1) + dx = x1 - self.mu_uv + dy = y1 - self.mu_uv + return (1.0 / self.l_uv**2) * ( + tau_x * tau_y * K_rbf * sm_x * sm_y * (dx * sp_x + dy * sp_y) + - sp_x * sp_y * f_x * f_y * (dx * sm_x + dy * sm_y) + ) + + grads.extend([_dK_dmu_uv, _dK_dl_uv]) + + # ---- IR asymptotic parameter gradients: mu_ir, l_ir ---------------- + if self.ir: + def _dK_dmu_ir(x, y): + x1 = make_column_vector(x[:, 0]) + y1 = make_row_vector(y[:, 0]) + tp_x = softtheta(x1, self.mu_ir, self.l_ir, -self.ir) # σ_+(x;ir) + tp_y = softtheta(y1, self.mu_ir, self.l_ir, -self.ir) + tm_x = softtheta(x1, self.mu_ir, self.l_ir, self.ir) # σ_-(x;ir) + tm_y = softtheta(y1, self.mu_ir, self.l_ir, self.ir) + sm_x = softtheta(x1, self.mu_uv, self.l_uv, -self.uv) + sm_y = softtheta(y1, self.mu_uv, self.l_uv, -self.uv) + K_rbf = self.kernel(x, y) + f_x = self.ir_asymptotics(x1) + f_y = self.ir_asymptotics(y1) + return (1.0 / self.l_ir) * ( + tm_x * tm_y * f_x * f_y * (tp_x + tp_y) + - tp_x * tp_y * K_rbf * sm_x * sm_y * (tm_x + tm_y) + ) + + def _dK_dl_ir(x, y): + x1 = make_column_vector(x[:, 0]) + y1 = make_row_vector(y[:, 0]) + tp_x = softtheta(x1, self.mu_ir, self.l_ir, -self.ir) + tp_y = softtheta(y1, self.mu_ir, self.l_ir, -self.ir) + tm_x = softtheta(x1, self.mu_ir, self.l_ir, self.ir) + tm_y = softtheta(y1, self.mu_ir, self.l_ir, self.ir) + sm_x = softtheta(x1, self.mu_uv, self.l_uv, -self.uv) + sm_y = softtheta(y1, self.mu_uv, self.l_uv, -self.uv) + K_rbf = self.kernel(x, y) + f_x = self.ir_asymptotics(x1) + f_y = self.ir_asymptotics(y1) + dx = x1 - self.mu_ir + dy = y1 - self.mu_ir + return (1.0 / self.l_ir**2) * ( + tm_x * tm_y * f_x * f_y * (dx * tp_x + dy * tp_y) + - tp_x * tp_y * K_rbf * sm_x * sm_y * (dx * tm_x + dy * tm_y) + ) + + grads.extend([_dK_dmu_ir, _dK_dl_ir]) + + return grads + class Matern12(Kernel): """Matern-1/2 kernel From d644bdbf74e7f073c71d33b5ab67399365ad870a Mon Sep 17 00:00:00 2001 From: Jonas Wessely <104987725+Jpwessely@users.noreply.github.com> Date: Tue, 2 Jun 2026 16:00:47 +0200 Subject: [PATCH 2/8] jacobian of log-trafo into internal weights --- fredipy/integrators.py | 34 ++++++++++++++++++---------------- 1 file changed, 18 insertions(+), 16 deletions(-) diff --git a/fredipy/integrators.py b/fredipy/integrators.py index 3d44456..d847f51 100644 --- a/fredipy/integrators.py +++ b/fredipy/integrators.py @@ -196,9 +196,8 @@ class Riemann_1D_log(Integrator): Uses midpoints equally spaced in log-space and a scalar step ``dw = d(log ω) = log(w_max/w_min) / int_n``. - Because the Jacobian of the substitution t = log ω is ω, every kernel - passed to this integrator must include that factor, i.e. the caller is - responsible for passing ``ω · K(p, ω)`` instead of ``K(p, ω)``. + The Jacobian of the substitution t = log ω is ω and is applied + internally, so plain kernels ``K(p, ω)`` can be passed directly. """ def __init__( @@ -211,6 +210,7 @@ def __init__( w_edges = make_column_vector(np.geomspace(w_min, w_max, int_n + 1)) self.dw = np.log(w_max / w_min) / int_n # uniform step in log-space self.w = np.sqrt(w_edges[:-1] * w_edges[1:]) # geometric midpoints + self.jac = make_row_vector(self.w) # log-space Jacobian: ω_i, shape (1, n) def doubleIntegrationSymmetric( self, @@ -226,9 +226,9 @@ def doubleIntegration( constraint2: LinearEquality ) -> np.ndarray: return self.dw**2 * ( - constraint1(make_row_vector(self.w), x=make_column_vector(constraint1.x)) + self.jac * constraint1(make_row_vector(self.w), x=make_column_vector(constraint1.x)) @ kernel(self.w, self.w) - @ (constraint2(make_row_vector(self.w), x=make_column_vector(constraint2.x))).T + @ (self.jac * constraint2(make_row_vector(self.w), x=make_column_vector(constraint2.x))).T ) def singleIntegration( @@ -238,7 +238,7 @@ def singleIntegration( w_pred: np.ndarray ) -> np.ndarray: return self.dw * ( - constraint(make_row_vector(self.w), x=make_column_vector(constraint.x)) + self.jac * constraint(make_row_vector(self.w), x=make_column_vector(constraint.x)) @ kernel(self.w, w_pred) ) @@ -253,9 +253,8 @@ class GaussLegendre_1D_log(Integrator): smooth integrands, requiring far fewer points than Riemann rules for the same accuracy (typically 50–200 instead of 1000). - As with ``Riemann_1D_log``, every kernel passed to this integrator must - include the log-space Jacobian factor ω, i.e. the caller is responsible - for passing ``ω · K(p, ω)`` instead of ``K(p, ω)``. + The log-space Jacobian factor ω is absorbed into the quadrature weights, + so plain kernels ``K(p, ω)`` can be passed directly. """ def __init__( @@ -275,8 +274,8 @@ def __init__( t = half_range * xi + 0.5 * (log_max + log_min) # nodes in log-space self.w = make_column_vector(np.exp(t)) # nodes in ω-space - # Weights absorb the half-range Jacobian of the linear change of variables - self.weights = make_row_vector(wi * half_range) # shape (1, int_n) + # Weights absorb the half-range Jacobian and the log-space Jacobian ω + self.weights = make_row_vector(wi * half_range) * self.w.T # shape (1, int_n) def doubleIntegrationSymmetric( self, @@ -322,6 +321,9 @@ class GaussLegendre_1D_log_UVtail(GaussLegendre_1D_log): T = (2 log w_uv)^{-13/22} / (2 × 13/22). + As with the parent class, plain kernels ``K(p, ω)`` are passed directly; + the log-space Jacobian ω is absorbed into the quadrature weights. + **Correction formulas** (A_θ = UV-component of the bulk integral): * ``doubleIntegrationSymmetric``: @@ -397,9 +399,9 @@ class GaussLegendre_1D_semiinf(Integrator): ω = w_scale · (1 + u) / (1 − u), u ∈ (−1, 1) which maps (−1, 1) exactly onto (0, ∞). The log-space Jacobian - d(log ω)/du = 2/(1 − u²) is absorbed into the quadrature weights, - so every kernel passed here must include the log-space factor ω, - i.e. pass ``ω · K(p, ω)`` instead of ``K(p, ω)``. + d(log ω)/du = 2/(1 − u²) and the ω factor are both absorbed into + the quadrature weights, so plain kernels ``K(p, ω)`` can be passed + directly. Parameters ---------- @@ -428,8 +430,8 @@ def __init__( # rational map (-1,1) → (0,∞) self.w = make_column_vector(w_scale * (1.0 + xi) / (1.0 - xi)) - # log-space weights: dlog(ω)/du = 2/(1-u²) - self.weights = make_row_vector(wi * 2.0 / (1.0 - xi**2)) + # weights: log-space d(log ω)/du = 2/(1-u²), plus ω Jacobian absorbed + self.weights = make_row_vector(wi * 2.0 / (1.0 - xi**2)) * self.w.T def doubleIntegrationSymmetric( self, From 27972c22959cf7a2b65b3cbc961514db0aa139d4 Mon Sep 17 00:00:00 2001 From: Jonas Wessely <104987725+Jpwessely@users.noreply.github.com> Date: Thu, 13 Aug 2026 11:56:53 +0200 Subject: [PATCH 3/8] predict/predict_data: clip variance before sqrt to avoid NaN np.sqrt(np.diag(C)) can receive a spuriously negative diagonal entry from floating-point cancellation in the posterior covariance (e.g. at constrained data points or deep in the UV/IR asymptotic tails where variance -> 0), returning nan instead of ~0. Clip to [0, inf) before the sqrt. See fredipy_diff_doc.md #6 in the reconstructions repo for full rationale. --- fredipy/models.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fredipy/models.py b/fredipy/models.py index 12c04b6..8b8fa26 100644 --- a/fredipy/models.py +++ b/fredipy/models.py @@ -106,7 +106,7 @@ def predict( if full_cov: return mu, C else: - return mu, np.sqrt(np.diag(C)) + return mu, np.sqrt(np.clip(np.diag(C), a_min=0.0, a_max=None)) def predict_data( self, @@ -142,7 +142,7 @@ def predict_data( if full_cov: return mu, C else: - return mu, np.sqrt(np.diag(C)) + return mu, np.sqrt(np.clip(np.diag(C), a_min=0.0, a_max=None)) def log_likelihood(self) -> float: """Returns the log-likelihood of the posterior GP.""" From cea49320defa79178207e622269f24bbd1df3b7a Mon Sep 17 00:00:00 2001 From: Jonas Wessely <104987725+Jpwessely@users.noreply.github.com> Date: Thu, 13 Aug 2026 11:56:56 +0200 Subject: [PATCH 4/8] softtheta: fix sign<0 branch overflowing to nan for x << mu0 The sign<0 branch computed exp(-(x-mu0)/l0)/(exp(-(x-mu0)/l0)+1) directly. For x far below mu0 relative to l0 (e.g. deep-IR kernel evaluations with a narrow l_ir), exp(...) overflows to inf in both numerator and denominator, giving inf/inf = nan even though the correct limit is exactly 1.0. Rewritten as the algebraic complement of the overflow-safe sign>0 branch instead. Adds a regression test. See fredipy_diff_doc.md #7 in the reconstructions repo for full rationale. --- fredipy/util.py | 6 +++++- tests/test_util.py | 13 +++++++++++++ 2 files changed, 18 insertions(+), 1 deletion(-) diff --git a/fredipy/util.py b/fredipy/util.py index c7c4d3a..addce83 100644 --- a/fredipy/util.py +++ b/fredipy/util.py @@ -53,7 +53,11 @@ def softtheta( return 1/(np.exp(-1*(x - mu0)/l0) + 1) elif sign < 0: - return np.exp(-1*(x - mu0)/l0)/(np.exp(-1*(x - mu0)/l0) + 1) + # Written as the exact algebraic complement of the sign>0 branch above, rather than as + # exp(...)/(exp(...)+1) directly: for x far below mu0 (deep cutoff region), the direct + # form divides inf by inf and returns nan, even though the correct limit is exactly 1. + # The complement form only ever evaluates the already-overflow-safe sign>0 branch. + return 1 - 1/(np.exp(-1*(x - mu0)/l0) + 1) else: return np.ones_like(x) diff --git a/tests/test_util.py b/tests/test_util.py index b7af1c9..c2d377c 100644 --- a/tests/test_util.py +++ b/tests/test_util.py @@ -27,6 +27,19 @@ def test_softtheta(): assert np.allclose(result, expected) +def test_softtheta_sign_negative_does_not_overflow_to_nan(): + # Regression test: x deep below mu0 relative to l0 (e.g. deep-IR kernel evaluations with a + # narrow l_ir) previously overflowed exp() in the sign<0 branch's direct exp/(exp+1) form, + # producing inf/inf = nan instead of the correct limit of 1.0. + x = np.array([0.0, 1e-3]) + mu0, l0 = 1.0, 1e-3 # (mu0 - x)/l0 >> 709 for both x values: exp() overflows + + result = softtheta(x, mu0, l0, -1) + + assert np.all(np.isfinite(result)) + assert np.allclose(result, 1.0) + + def test_allclose(): a = np.array([1, 2, 3]) b = np.array([1, 2, 3]) From 95aa10b6f27e1c616fa71865df41390cd257a834 Mon Sep 17 00:00:00 2001 From: Jonas Wessely <104987725+Jpwessely@users.noreply.github.com> Date: Thu, 20 Aug 2026 00:25:32 +0200 Subject: [PATCH 5/8] UVtail: apply the analytic tail correction to cross-covariance blocks GaussLegendre_1D_log_UVtail overrode singleIntegration and doubleIntegrationSymmetric but not doubleIntegration, so the cross-covariance block between the sum-rule constraint and any other integral constraint silently dropped the analytic tail correction that both the diagonal block and the prediction operator carried. The GP covariance matrix was therefore inconsistent with the operator predicted through, making predict() and predict_data() describe different models (44% relative error on the reintegration identity: reintegrating the posterior mean over the integrator's own nodes must reproduce predict_data(), and did not). Derivation: for w' > w_uv >> mu_uv the kernel is rank-1 in the tail, so splitting both integrals at w_uv gives Sigma_12 = BB + A_1 T_2^T + T_1 A_2^T + T_1 T_2^T with A the bulk-UV overlap and T the tail moment, both per-constraint. The existing symmetric formula bulk + 2AT + T^2 is the C_1 = C_2 special case, so the override is deleted rather than kept in parallel. T is a per-constraint, per-row vector -- the sum-rule and KL kernels have different UV falloffs, so no single scalar serves both sides. The correction lives in GaussLegendre_1D_log driven by per-constraint uv_tail_moment() / uv_anchor() hooks, so blocks (i,j) and (j,i) are the same expression with roles swapped and symmetry is structural. This matters because models.py's np.linalg.cholesky reads only the lower triangle and would silently accept an asymmetric matrix. A node guard raises NotImplementedError for genuinely different grids, and an all-zero-T fast path keeps every non-UVtail model bit-identical. Also fixed here: - T's exponent 13/22 was hardcoded and gluon-specific, but T is the sum-rule-weighted tail and depends on the observable's anomalous dimension; the ghost project (9/44) was off by a factor 10.4. tail_moment is now a required constructor argument -- deliberately no default, since the silent default is what caused this. - singleIntegration's tail was a broadcast where the general form is an outer product (latent: the sum-rule constraint has a single row). - The analytic NLL gradient added T1 @ T2.T unchanged, but that term is hyperparameter-independent so its derivative is zero. Pre-existing (the deleted symmetric override did the same to the sr-sr block); corrected via an explicit gradient mode. Verified: reintegration identity 4.02e-16 on a well-conditioned model (4.4e-2 uncorrected, so the test is not vacuous); covariance asymmetry 6.66e-16 against matrix scale 6.2; analytic gradient agrees with finite differences to 5.2e-07 for all six hyperparameters. New regression suite in tests/test_integrators_uvtail.py; 83 passed. Co-Authored-By: Claude Sonnet 5 --- fredipy/covariance.py | 47 ++- fredipy/integrators.py | 356 ++++++++++++++++++---- fredipy/models.py | 5 +- tests/test_integrators_uvtail.py | 495 +++++++++++++++++++++++++++++++ 4 files changed, 829 insertions(+), 74 deletions(-) create mode 100644 tests/test_integrators_uvtail.py diff --git a/fredipy/covariance.py b/fredipy/covariance.py index a11a91f..a70af8b 100644 --- a/fredipy/covariance.py +++ b/fredipy/covariance.py @@ -20,6 +20,17 @@ class TwoSided: New rules for user-defined operators can be defined by creating a class that inherits from this class, adding the appropriate functions following the naming scheme, and passing an instance to the model constructor. + + Every combiner takes a ``derivative`` flag, which is threaded through + unchanged from :meth:`__call__`. It is ``False`` when the covariance matrix + itself is being assembled and ``True`` when ``kernel`` is a *derivative* of + the kernel with respect to one hyperparameter (see + ``GaussianProcess.log_likelihood_grad``). Combiners that build a block by + plain linear application of ``kernel`` can ignore it; only pieces that add a + term which is *not* linear in ``kernel`` -- e.g. the constant analytic UV + tail moment in ``GaussLegendre_1D_log._uv_tail_correction``, whose + hyperparameter derivative is zero rather than itself -- need to act on it. + Custom combiners must accept the flag. """ def __init__(self): pass @@ -27,7 +38,8 @@ def __init__(self): def __call__( self, kernel: Kernel, - constraints: List[LinearEquality] + constraints: List[LinearEquality], + derivative: bool = False ) -> np.ndarray: rows = [] for c1 in constraints: @@ -36,9 +48,9 @@ def __call__( combiner12 = getattr(self, f"_{type(c1.op).__name__}_{type(c2.op).__name__}", None) combiner21 = getattr(self, f"_{type(c2.op).__name__}_{type(c1.op).__name__}", None) if combiner12: - entry = combiner12(c1, kernel, c2) + entry = combiner12(c1, kernel, c2, derivative=derivative) elif combiner21: - entry = combiner21(c2, kernel, c1).T + entry = combiner21(c2, kernel, c1, derivative=derivative).T else: raise NotImplementedError( f"No rule found to combine operators of types \ @@ -47,28 +59,39 @@ def __call__( rows.append(np.concatenate(columns, axis=1)) return np.concatenate(rows) - def _Integral_Integral(self, c1, k, c2): + def _Integral_Integral(self, c1, k, c2, derivative: bool = False): + # NOTE: c1's integrator drives both sides of the cross block, so the two constraints + # must discretise omega compatibly. That is not checked here: integrators that apply an + # analytic UV tail (GaussLegendre_1D_log._uv_tail_correction) enforce matching nodes and + # weights downstream and raise NotImplementedError otherwise. Blocks (i, j) and (j, i) + # are built independently, each through its own left constraint's integrator, so any + # correction added downstream must be symmetric by construction. + # `derivative` is forwarded because that UV tail correction contains one term (T1 T2^T) + # that is independent of the kernel and must therefore vanish, not be reproduced, when + # k is dK/dtheta. if c1 == c2: - return c1.op.integrator.doubleIntegrationSymmetric(c1, k) + return c1.op.integrator.doubleIntegrationSymmetric(c1, k, derivative=derivative) else: - return c1.op.integrator.doubleIntegration(c1, k, c2) + return c1.op.integrator.doubleIntegration(c1, k, c2, derivative=derivative) - def _Identity_Identity(self, c1, k, c2): + def _Identity_Identity(self, c1, k, c2, derivative: bool = False): return k(c1.x, c2.x) - def _Derivative_Derivative(self, c1, k, c2): + def _Derivative_Derivative(self, c1, k, c2, derivative: bool = False): return k.d2K_dxdy(c1.x, c2.x) - def _Integral_Identity(self, c1, k, c2): + def _Integral_Identity(self, c1, k, c2, derivative: bool = False): + # No `derivative` handling needed: singleIntegration's analytic tail, T @ K(w_uv, x) / + # f_uv(w_uv), is linear in the kernel, so substituting dK/dtheta already differentiates it. return c1.op.integrator.singleIntegration(c1, k, c2.x) - def _Integral_Derivative(self, c1, k, c2): + def _Integral_Derivative(self, c1, k, c2, derivative: bool = False): return c1.op.integrator.singleIntegration(c1, k.dK_dy, c2.x) - def _Identity_Derivative(self, c1, k, c2): + def _Identity_Derivative(self, c1, k, c2, derivative: bool = False): return k.dK_dy(c1.x, c2.x) - def _ConstMatrix_ConstMatrix(self, c1, k, c2): + def _ConstMatrix_ConstMatrix(self, c1, k, c2, derivative: bool = False): return c1.op() @ k(c1.x, c2.x) @ c2.op().T diff --git a/fredipy/integrators.py b/fredipy/integrators.py index d847f51..382201d 100644 --- a/fredipy/integrators.py +++ b/fredipy/integrators.py @@ -21,7 +21,8 @@ def __init__( def doubleIntegrationSymmetric( self, constraint: LinearEquality, - kernel: Callable + kernel: Callable, + derivative: bool = False ) -> np.ndarray: r"""Symmetric double integration of a constraint C with a GP kernel K, i.e. @@ -31,6 +32,8 @@ def doubleIntegrationSymmetric( ---------- constraint : Some integral constraint kernel : Some Gaussian process kernel + derivative : True when ``kernel`` is dK/dtheta rather than K, see + :meth:`doubleIntegration`. Returns ------- @@ -42,7 +45,8 @@ def doubleIntegration( self, constraint1: LinearEquality, kernel: Callable, - constraint2: LinearEquality + constraint2: LinearEquality, + derivative: bool = False ) -> np.ndarray: r"""Double integration of two (different) constraints C1 and C2 with a GP kernel K, i.e. @@ -52,6 +56,16 @@ def doubleIntegration( ---------- constraint1, constraint2 : Some integral constraints kernel : Some Gaussian process kernel + derivative : Flag set by ``GaussianProcess.log_likelihood_grad`` (via + ``TwoSided``) to say that ``kernel`` is the derivative + dK/dtheta of the kernel with respect to one hyperparameter + rather than the kernel itself. Purely quadrature-based + implementations are linear in ``kernel`` and can ignore it; + any term that is *not* linear in ``kernel`` (the constant + analytic tail moment of + :meth:`GaussLegendre_1D_log._uv_tail_correction`) must be + dropped when it is set, since the derivative of a constant + is zero. Returns ------- @@ -81,6 +95,61 @@ def singleIntegration( """ raise NotImplementedError + def uv_tail_moment( + self, + constraint: LinearEquality + ) -> np.ndarray: + r"""Analytic UV tail moment ``T`` of this integrator for one constraint. + + .. math:: + T_i[m] = \int_{w_{uv}}^{\infty} C_i(p_m, w) f_{uv}(w)\, dw + + i.e. the constraint kernel weighted by the *bare* UV asymptotic shape + ``f_uv`` over the part of the ω axis that the quadrature grid does not + cover. It is the only ingredient (besides the bulk-UV overlap ``A``, + see :meth:`uv_anchor`) of the general two-constraint tail correction + + .. math:: + \Sigma_{12} = BB + A_1 T_2^T + T_1 A_2^T + T_1 T_2^T \qquad (**) + + implemented in :meth:`GaussLegendre_1D_log._uv_tail_correction`. + + ``T`` is per-constraint **and** per-row: two constraints whose kernels + have different UV falloffs have different tail moments, so a single + shared scalar is *not* valid for the cross block (see + ``docs/fredipy_uvtail_fix_plan.md`` in the reconstructions repo). + + Returns + ------- + 2D array of shape ``(len(constraint.x), 1)``. The default is all + zeros — integrators that carry no analytic tail, which lets + ``_uv_tail_correction`` short-circuit and stay bit-identical to the + uncorrected quadrature. + """ + return np.zeros((constraint.x.shape[0], 1)) + + def uv_anchor(self) -> tuple | None: + r"""Anchor point used to evaluate the bulk-UV overlap ``A`` of ``(**)``. + + .. math:: + A_i[m] = \frac{1}{f_{uv}(w_{uv})} + \sum_j W_j\, C_i(p_m, w_j)\, K(w_j, w_{uv}) + + which is *exact*, not approximate, whenever ``w_uv >> mu_uv``: the UV + branch of ``AsymptoticKernel`` factorises there as + ``K(w, w_uv) = θ_uv(w) f_uv(w) f_uv(w_uv)``, so dividing by the anchor + value ``f_uv(w_uv)`` recovers ``Σ_j W_j C_i(w_j) θ_uv(w_j) f_uv(w_j)`` + to all digits. + + Returns + ------- + ``None`` for tail-free integrators (the default). Tail-carrying + integrators return ``(w_uv, f_uv_anchor)`` with ``w_uv`` a ``(1, 1)`` + column vector holding the UV split point and ``f_uv_anchor`` the scalar + ``f_uv(w_uv)``. + """ + return None + class Riemann(Integrator): """Implementation of the Riemann integration in arbitrary dimensions and kernels. @@ -103,7 +172,8 @@ def __init__( def doubleIntegrationSymmetric( self, constraint: LinearEquality, - kernel: Callable + kernel: Callable, + derivative: bool = False # unused: this rule is linear in `kernel` ) -> np.ndarray: p = constraint.x @@ -161,15 +231,17 @@ def __init__( def doubleIntegrationSymmetric( self, constraint: LinearEquality, - kernel: Callable + kernel: Callable, + derivative: bool = False # unused: this rule is linear in `kernel` ) -> np.ndarray: - return self.doubleIntegration(constraint, kernel, constraint) + return self.doubleIntegration(constraint, kernel, constraint, derivative=derivative) def doubleIntegration( self, constraint1: LinearEquality, kernel: Callable, - constraint2: LinearEquality + constraint2: LinearEquality, + derivative: bool = False # unused: this rule is linear in `kernel` ) -> np.ndarray: return self.dw**2 * ( constraint1(make_row_vector(self.w), x=make_column_vector(constraint1.x)) @@ -215,15 +287,17 @@ def __init__( def doubleIntegrationSymmetric( self, constraint: LinearEquality, - kernel: Callable + kernel: Callable, + derivative: bool = False # unused: this rule is linear in `kernel` ) -> np.ndarray: - return self.doubleIntegration(constraint, kernel, constraint) + return self.doubleIntegration(constraint, kernel, constraint, derivative=derivative) def doubleIntegration( self, constraint1: LinearEquality, kernel: Callable, - constraint2: LinearEquality + constraint2: LinearEquality, + derivative: bool = False # unused: this rule is linear in `kernel` ) -> np.ndarray: return self.dw**2 * ( self.jac * constraint1(make_row_vector(self.w), x=make_column_vector(constraint1.x)) @@ -264,6 +338,12 @@ def __init__( int_n: int = 100 ) -> None: + # Grid definition kept as attributes purely for diagnostics: the node-compatibility + # guard in _uv_tail_correction names these in its error message. + self.w_min = w_min # lower end of the quadrature range in ω + self.w_max = w_max # upper end of the quadrature range in ω + self.int_n = int_n # number of Gauss-Legendre nodes + # GL nodes (xi in [-1,1]) and weights xi, wi = np.polynomial.legendre.leggauss(int_n) @@ -280,21 +360,29 @@ def __init__( def doubleIntegrationSymmetric( self, constraint: LinearEquality, - kernel: Callable + kernel: Callable, + derivative: bool = False ) -> np.ndarray: - return self.doubleIntegration(constraint, kernel, constraint) + return self.doubleIntegration(constraint, kernel, constraint, derivative=derivative) def doubleIntegration( self, constraint1: LinearEquality, kernel: Callable, - constraint2: LinearEquality + constraint2: LinearEquality, + derivative: bool = False ) -> np.ndarray: - return ( + bulk = ( self.weights * constraint1(make_row_vector(self.w), x=make_column_vector(constraint1.x)) @ kernel(self.w, self.w) @ (self.weights * constraint2(make_row_vector(self.w), x=make_column_vector(constraint2.x))).T ) + # The analytic UV tail lives here (not as an override on the UVtail subclass) so that + # both orderings of a constraint pair go through the same code and the assembled + # covariance matrix is symmetric *by construction* — np.linalg.cholesky reads only the + # lower triangle and would silently accept an asymmetric matrix otherwise. + return bulk + self._uv_tail_correction( + constraint1, kernel, constraint2, derivative=derivative) def singleIntegration( self, @@ -307,42 +395,182 @@ def singleIntegration( @ kernel(self.w, w_pred) ) + @staticmethod + def _same_quadrature( + integrator1: Integrator, + integrator2: Integrator + ) -> bool: + """True if two integrators share (numerically) the same nodes and weights.""" + w1, w2 = getattr(integrator1, 'w', None), getattr(integrator2, 'w', None) + v1, v2 = getattr(integrator1, 'weights', None), getattr(integrator2, 'weights', None) + if w1 is None or w2 is None or v1 is None or v2 is None: + return False + if w1.shape != w2.shape or v1.shape != v2.shape: + return False + return bool(np.allclose(w1, w2) and np.allclose(v1, v2)) + + @staticmethod + def _describe_grid(integrator: Integrator) -> str: + """Human-readable grid summary used in the node-mismatch error message.""" + return (f"{type(integrator).__name__}(w_min={getattr(integrator, 'w_min', '?')}, " + f"w_max={getattr(integrator, 'w_max', '?')}, " + f"int_n={getattr(integrator, 'int_n', '?')})") + + def _uv_tail_correction( + self, + constraint1: LinearEquality, + kernel: Callable, + constraint2: LinearEquality, + derivative: bool = False + ) -> np.ndarray | float: + r"""Analytic UV tail correction to the (constraint1, constraint2) covariance block. + + Splitting *each* of the two integrals at the UV point ``w_uv`` into a bulk part + ``(w_min, w_uv)`` — which the quadrature grid covers — and a tail part + ``(w_uv, ∞)`` — which it does not — and using that ``AsymptoticKernel`` + factorises as ``K(w, w') = [θ_uv(w) f_uv(w)] f_uv(w')`` whenever ``w' > w_uv`` + with ``w_uv >> mu_uv`` (rank-1 in the tail, exact up to + ``O(exp(-(w_uv - mu_uv)/l_uv))``), the four pieces give + + .. math:: + \Sigma_{12} = BB + A_1 T_2^T + T_1 A_2^T + T_1 T_2^T \qquad (**) + + with ``BB`` the plain double quadrature computed by + :meth:`doubleIntegration`, ``A_i`` the bulk-UV overlap (see + :meth:`Integrator.uv_anchor`) and ``T_i`` the per-constraint, per-row tail + moment (see :meth:`Integrator.uv_tail_moment`). + + Setting ``C_1 = C_2`` and ``T`` row-independent collapses ``(**)`` to the + legacy symmetric closed form ``bulk + 2 A T + T²``, so the symmetric case is + a special case of this method and needs no separate override. + + Gradient mode (``derivative=True``) + ----------------------------------- + ``GaussianProcess.log_likelihood_grad`` reassembles the same block with + ``dK/dtheta`` substituted for ``K``. Under that substitution the two ``A`` + terms differentiate themselves correctly — ``A_i`` is linear in the kernel + (``A_i = W C_i K(ω, w_uv) / f_uv(w_uv)`` with a hyperparameter-independent + anchor value), so ``A_i[dK/dtheta] = dA_i/dtheta`` exactly. The term + ``T_1 T_2^T`` however contains no kernel at all: ``T`` is the *analytic* + tail moment, a function of ``w_uv`` and the anomalous dimension only, so + ``d(T_1 T_2^T)/dtheta = 0``. Adding it unchanged in gradient mode would + treat the derivative of a constant as the constant itself, which is why it + is dropped here. The likelihood *value* is unaffected either way; only the + gradient is. + """ + integrator1 = constraint1.op.integrator + integrator2 = constraint2.op.integrator + tail1 = integrator1.uv_tail_moment(constraint1) + tail2 = integrator2.uv_tail_moment(constraint2) + + # Fast path: no constraint carries an analytic tail, so (**) reduces to BB. + # This short circuit is load-bearing — it keeps every plain-quadrature model + # (three-gluon-vert, ghost-gluon-vert, fgvert-sd, all pre-existing fredipy + # tests) bit-identical to the uncorrected result. Do not drop it. + if not np.any(tail1) and not np.any(tail2): + return 0.0 + + # (**) mixes the two constraints on a *shared* set of quadrature nodes; it is + # not well defined if the two integrators discretise ω differently. + if integrator1 is not integrator2 and not self._same_quadrature(integrator1, integrator2): + raise NotImplementedError( + "UV tail correction requires both integral constraints to share the same " + "quadrature nodes and weights, but got " + f"{self._describe_grid(integrator1)} and {self._describe_grid(integrator2)}." + ) + + anchors = [a for a in (integrator1.uv_anchor(), integrator2.uv_anchor()) if a is not None] + if not anchors: + raise NotImplementedError( + "A non-zero UV tail moment was reported but no integrator supplied a UV " + f"anchor via uv_anchor(): {self._describe_grid(integrator1)} and " + f"{self._describe_grid(integrator2)}." + ) + w_uv, f_uv_anchor = anchors[0] + for other_w_uv, other_f_uv in anchors[1:]: + assert np.allclose(w_uv, other_w_uv) and np.isclose(f_uv_anchor, other_f_uv), \ + "Integral constraints report incompatible UV anchors." + + # A_i = (W C_i) @ K(ω, w_uv) / f_uv(w_uv), both evaluated on *this* integrator's + # nodes (guaranteed above to coincide with the other integrator's). + c1_row = constraint1(make_row_vector(self.w), x=make_column_vector(constraint1.x)) + c2_row = constraint2(make_row_vector(self.w), x=make_column_vector(constraint2.x)) + K_col = kernel(self.w, w_uv) # (N, 1) + A1 = (self.weights * c1_row @ K_col) / f_uv_anchor # (M1, 1) + A2 = (self.weights * c2_row @ K_col) / f_uv_anchor # (M2, 1) + + correction = A1 @ tail2.T + tail1 @ A2.T + if derivative: + # d(T_1 T_2^T)/dtheta = 0: T carries no hyperparameter dependence. See the + # "Gradient mode" section of the docstring. + return correction + return correction + tail1 @ tail2.T + class GaussLegendre_1D_log_UVtail(GaussLegendre_1D_log): - """GL log-space quadrature on (w_min, w_uv) with an analytic UV tail. + """GL log-space quadrature on (w_min, w_uv) plus an analytic UV tail. - For ω > w_uv the AsymptoticKernel factorises as - K(ω, ω') ≈ f_uv(ω) f_uv(ω'), - so the analytic tail T = ∫_{log w_uv}^∞ (2t)^{-35/22} dt can be - folded in exactly via corrections to the integration methods. + **Rank-1 tail assumption.** For ω' > w_uv with w_uv >> mu_uv the soft + thetas of ``AsymptoticKernel`` saturate (θ_uv(ω') = 1 to machine + precision), so the kernel becomes rank-1 in the tail: - The substitution u = (2 log ω)^{-13/22} maps the UV integrand to a - constant, giving the closed form + K(ω, ω') = [θ_uv(ω) f_uv(ω)] · f_uv(ω') (*) - T = (2 log w_uv)^{-13/22} / (2 × 13/22). + This is the *only* assumption made here; it holds up to + ``O(exp(-(w_uv - mu_uv)/l_uv))``. - As with the parent class, plain kernels ``K(p, ω)`` are passed directly; - the log-space Jacobian ω is absorbed into the quadrature weights. + **General two-constraint correction.** Splitting each of the two integrals + of a covariance block at w_uv into bulk (w_min, w_uv) and tail (w_uv, ∞) and + applying (*) to whichever argument lies in the tail gives - **Correction formulas** (A_θ = UV-component of the bulk integral): + A_i[m] = ∫_bulk C_i(p_m, ω) θ_uv(ω) f_uv(ω) dω + T_i[m] = ∫_{w_uv}^∞ C_i(p_m, ω) f_uv(ω) dω - * ``doubleIntegrationSymmetric``: - full = bulk + 2 A_θ T + T² - where A_θ[m] ≈ (1/f_uv(w_uv)) × Σ_i W_i C_m(ω_i) K(ω_i, w_uv) + Σ_12 = BB + A_1 T_2^T + T_1 A_2^T + T_1 T_2^T (**) - * ``singleIntegration``: - full = numerical + T × K(w_uv, ω_pred) / f_uv(w_uv) - (K(w_uv, ω_pred)/f_uv(w_uv) ≈ θ_uv(ω_pred) f_uv(ω_pred) for - prediction points far from w_uv, where the RBF part of K vanishes) + where ``BB`` is the plain double quadrature of the parent class. Two points + that are easy to get wrong: + + * ``T`` is **per constraint and per row**, not a single shared scalar: two + constraints whose kernels have different UV falloffs (e.g. a sum-rule + kernel ``C(p, ω) = ω`` versus a Källén-Lehmann data kernel) have entirely + different tail moments. This is why ``tail_moment`` is a required + constructor argument rather than a hardcoded closed form. + * ``A`` obtained via the anchor, ``A = (W C @ K(ω, w_uv)) / f_uv(w_uv)``, is + **exact**, not approximate: by (*) the anchor column is + ``K(ω_i, w_uv) = θ_uv(ω_i) f_uv(ω_i) f_uv(w_uv)``, so the division + recovers the bulk-UV overlap to all digits. + + The symmetric case ``C_1 = C_2`` with row-independent ``T`` collapses (**) + to the legacy closed form ``bulk + 2 A T + T²``, so this class deliberately + does **not** override ``doubleIntegration``/``doubleIntegrationSymmetric``: + the correction is applied by the parent class for *both* orderings of a + constraint pair, which is what keeps the assembled covariance matrix + symmetric by construction (``np.linalg.cholesky`` reads only the lower + triangle and would silently accept an asymmetric matrix). + + As with the parent class, plain kernels ``K(p, ω)`` are passed directly; + the log-space Jacobian ω is absorbed into the quadrature weights. Parameters ---------- w_min : float + Lower end of the quadrature range in ω. w_uv : float - UV split where UV asymptotics fully apply. Recommended: ≥ 1000. + UV split where the UV asymptotics fully apply, i.e. the upper end of the + quadrature range; everything above it is covered by ``tail_moment``. + Recommended: ≥ 1000, and in any case >> mu_uv of the kernel. int_n : int + Number of Gauss-Legendre nodes on (w_min, w_uv). uv_func : callable - The UV asymptotic function f_uv(ω), e.g. ``uv_asymptotics``. + The UV asymptotic function f_uv(ω), e.g. ``uv_asymptotics``. Used only + to evaluate the anchor value f_uv(w_uv). + tail_moment : callable + Maps ``constraint.x`` of shape (M, d) to the tail moments T of shape + (M, 1), i.e. ``T[m] = ∫_{w_uv}^∞ C(p_m, ω) f_uv(ω) dω`` for the + constraint this integrator is attached to. Required (no default): a + silent default is exactly what let a gluon-specific closed form be + applied to the ghost sum rule, off by a factor 10.4. """ def __init__( @@ -350,33 +578,33 @@ def __init__( w_min: float, w_uv: float, int_n: int, - uv_func: Callable + uv_func: Callable, + tail_moment: Callable[[np.ndarray], np.ndarray] ) -> None: super().__init__(w_min, w_uv, int_n) - self.uv_func = uv_func + self.uv_func = uv_func # UV asymptotic shape f_uv(ω), used for the anchor value + self.tail_moment = tail_moment # constraint.x -> (M, 1) analytic tail moments T # Anchor point at the UV split boundary self.w_uv = make_column_vector(np.array([w_uv])) # (1, 1) self.f_uv_anchor = float(uv_func(self.w_uv).flat[0]) # scalar: f_uv(w_uv) - # Analytic tail: ∫_{log w_uv}^∞ (2t)^{-35/22} dt - self.T_uv = (2.0 * np.log(w_uv)) ** (-13.0 / 22.0) / (2.0 * 13.0 / 22.0) - def doubleIntegrationSymmetric( + def uv_tail_moment( self, - constraint: LinearEquality, - kernel: Callable + constraint: LinearEquality ) -> np.ndarray: - bulk = super().doubleIntegrationSymmetric(constraint, kernel) # (M, M) - - # A_θ ≈ (W C @ K_col) / f_uv(w_uv), shape (M, 1) - # K_col[i] = K(ω_i, w_uv) ≈ θ_uv(ω_i) f_uv(ω_i) f_uv(w_uv) (for large w_uv) - K_col = kernel(self.w, self.w_uv) # (N, 1) - c_row = constraint(make_row_vector(self.w), x=make_column_vector(constraint.x)) - A = (self.weights * c_row @ K_col) / self.f_uv_anchor # (M, 1) + tail = np.asarray(self.tail_moment(constraint.x), dtype=float) + if tail.ndim == 1: + tail = tail.reshape(-1, 1) + expected = (constraint.x.shape[0], 1) + if tail.shape != expected: + raise ValueError( + f"tail_moment must return an array of shape {expected} for this constraint, " + f"got {tail.shape}." + ) + return tail - # correction_{mn} = T (A_m + A_n) + T² - ones_col = np.ones((A.shape[0], 1)) - correction = self.T_uv * (A @ ones_col.T + ones_col @ A.T) + self.T_uv ** 2 - return bulk + correction + def uv_anchor(self) -> tuple: + return (self.w_uv, self.f_uv_anchor) def singleIntegration( self, @@ -385,9 +613,10 @@ def singleIntegration( w_pred: np.ndarray ) -> np.ndarray: numerical = super().singleIntegration(constraint, kernel, w_pred) - # tail = T × K(w_uv, ω_pred) / f_uv(w_uv), shape (1, N_pred) - # K(w_uv, ω_pred) / f_uv(w_uv) ≈ θ_uv(ω_pred) f_uv(ω_pred) for w_uv >> ω_pred - tail = (self.T_uv / self.f_uv_anchor) * kernel(self.w_uv, w_pred) # (1, N_pred) + # tail[m, n] = T[m] × K(w_uv, ω_pred[n]) / f_uv(w_uv) — an outer product, since T is + # per-row: broadcasting a single scalar T is only correct for a one-row constraint. + # K(w_uv, ω_pred) / f_uv(w_uv) ≈ θ_uv(ω_pred) f_uv(ω_pred) for w_uv >> ω_pred. + tail = self.uv_tail_moment(constraint) @ (kernel(self.w_uv, w_pred) / self.f_uv_anchor) return numerical + tail @@ -436,15 +665,17 @@ def __init__( def doubleIntegrationSymmetric( self, constraint: LinearEquality, - kernel: Callable + kernel: Callable, + derivative: bool = False # unused: this rule is linear in `kernel` ) -> np.ndarray: - return self.doubleIntegration(constraint, kernel, constraint) + return self.doubleIntegration(constraint, kernel, constraint, derivative=derivative) def doubleIntegration( self, constraint1: LinearEquality, kernel: Callable, - constraint2: LinearEquality + constraint2: LinearEquality, + derivative: bool = False # unused: this rule is linear in `kernel` ) -> np.ndarray: return ( self.weights * constraint1(make_row_vector(self.w), x=make_column_vector(constraint1.x)) @@ -489,15 +720,17 @@ def __init__( def doubleIntegrationSymmetric( self, constraint: LinearEquality, - kernel: Callable + kernel: Callable, + derivative: bool = False # unused: this rule is linear in `kernel` ) -> np.ndarray: - return self.doubleIntegration(constraint, kernel, constraint) + return self.doubleIntegration(constraint, kernel, constraint, derivative=derivative) def doubleIntegration( self, constraint1: LinearEquality, kernel: Callable, - constraint2: LinearEquality + constraint2: LinearEquality, + derivative: bool = False # unused: this rule is linear in `kernel` ) -> np.ndarray: return self.dw**2 / 9 * ( @@ -544,7 +777,8 @@ def __init__( def doubleIntegrationSymmetric( self, constraint: LinearEquality, - kernel: Callable + kernel: Callable, + derivative: bool = False # unused: this rule is linear in `kernel` ) -> np.ndarray: p = constraint.x diff --git a/fredipy/models.py b/fredipy/models.py index 8b8fa26..7d1f9aa 100644 --- a/fredipy/models.py +++ b/fredipy/models.py @@ -175,7 +175,10 @@ def log_likelihood_grad(self) -> List[float]: kernel_grad = self.kernel.params_gradient() for i in range(self.kernel.dim): - OpKer_grad = self.OpKerOp(kernel_grad[i], self.constraints) + # derivative=True tells the covariance machinery that kernel_grad[i] is dK/dtheta_i, + # not K: blocks assembled with terms that do not depend on the kernel (the constant + # analytic UV tail moment T) must drop those terms rather than reproduce them. + OpKer_grad = self.OpKerOp(kernel_grad[i], self.constraints, derivative=True) loglik_grad[i] = 0.5 * np.trace((alpha @ alpha.T - K_inv) @ OpKer_grad) return loglik_grad diff --git a/tests/test_integrators_uvtail.py b/tests/test_integrators_uvtail.py new file mode 100644 index 0000000..74b9bb2 --- /dev/null +++ b/tests/test_integrators_uvtail.py @@ -0,0 +1,495 @@ +"""Coverage for ``GaussLegendre_1D_log``'s analytic UV tail correction. + +The correction implements + + Sigma_12 = BB + A_1 T_2^T + T_1 A_2^T + T_1 T_2^T (**) + +for a pair of integral constraints, where BB is the plain double quadrature, T_i are the +per-constraint, per-row analytic tail moments above the UV split point w_uv, and A_i are the +bulk-UV overlaps evaluated through the anchor K(w, w_uv) / f_uv(w_uv). + +The tests below pin down, in order: that (**) makes the cross block consistent with the +prediction operator (T1), that the assembled covariance stays symmetric (T2), that the legacy +symmetric closed form is unchanged (T3), that the sum-rule tail moment is the right analytic +value (T4), that the end-to-end posterior is self-consistent (T5), that plain quadrature models +are bit-identical to their pre-fix values (T6), and the two error/shape guards (T7, T8). +""" + +import numpy as np +import pytest + +from fredipy import constraints, integrators, kernels, models, operators +from fredipy.covariance import OneSided, TwoSided +from fredipy.util import make_column_vector, make_row_vector + + +# -------------------------------------------------------------------------------------------- +# Synthetic, deliberately well-conditioned model. +# +# Tolerances in an end-to-end test are limited by cond(OpKerOp + cov_y); the production gluon and +# ghost models run at data_cov = 1e-10 where cond ~ 7e12 and no identity can be checked below +# ~1e-5. These settings keep the model far from that regime so the tests measure the correction, +# not the conditioning. +# -------------------------------------------------------------------------------------------- + +W_MIN = 1e-3 # Lower end of the omega quadrature range. +W_UV = 1e5 # Upper end of the quadrature range == UV split point w_uv. +INT_N = 120 # Number of Gauss-Legendre nodes on (W_MIN, W_UV). +DATA_COV = 1e-6 # Diagonal data covariance; large enough to keep the model well conditioned. +SR_COV = 1e-8 # Sum-rule target covariance. +N_DATA = 15 # Number of synthetic momentum points. +P_MIN, P_MAX = 0.5, 20.0 # Momentum window of the synthetic data, in the same units as omega. + +RBF_VARIANCE = 2.0 # Base RBF variance of the synthetic kernel. +RBF_LENGTHSCALE = 1.5 # Base RBF lengthscale of the synthetic kernel. +MU_UV = 1.1 # Centre of the UV blend-in; must satisfy W_UV >> MU_UV for (*) to hold. +L_UV = 0.15 # Width of the UV blend-in. + +GAMMA_GLUON = 13.0 / 22.0 # Gluon anomalous dimension (the exponent formerly hardcoded). +GAMMA_GHOST = 9.0 / 44.0 # Ghost anomalous dimension; differs, which is why it is a parameter. + + +def propagator_uv_shape(w, gamma=GAMMA_GLUON): + """Propagator-type UV asymptotic shape f_uv(w) = w^-2 (2 log w)^-(1+gamma). + + Uses the same smooth-maximum regularisation as the reconstructions repo's ``uv_asymptotics`` + so that the log stays finite below w = 1; the regularisation is irrelevant at w >= w_uv. + """ + delta = 0.05 # Smoothing width of the regularised maximum. + w1 = 0.5 * (w + 1.15 + np.sqrt((w - 1.15) ** 2 + delta ** 2)) + return 1.0 / (w1 ** 2 * np.log(w1 ** 2) ** (1.0 + gamma)) + + +def ghost_uv_shape(w): + """Propagator-type UV shape at the ghost anomalous dimension.""" + return propagator_uv_shape(w, gamma=GAMMA_GHOST) + + +def kl_p2(p, w): + """Kallen-Lehmann data kernel multiplied by p^2.""" + return p ** 2 * w / (w ** 2 + p ** 2) / np.pi + + +def sr_kernel(p, w): + """Sum-rule kernel C(p, w) = w, broadcast over the (ignored) momentum axis.""" + return w * np.ones_like(p) + + +def analytic_sum_rule_tail(w_uv, gamma): + """Closed form (2 log w_uv)^-gamma / (2 gamma) of the sum-rule tail moment.""" + return (2.0 * np.log(w_uv)) ** (-gamma) / (2.0 * gamma) + + +def sum_rule_tail_moment(w_uv, gamma): + """Row-independent ``tail_moment`` callable for the sum-rule kernel C(p, w) = w.""" + value = analytic_sum_rule_tail(w_uv, gamma) + + def _tail_moment(x): + return np.full((np.atleast_2d(x).shape[0], 1), value) + + return _tail_moment + + +def make_kernel(uv_shape=propagator_uv_shape): + """AsymptoticKernel with a UV branch only (no IR branch) at the synthetic hyperparameters.""" + rbf = kernels.RadialBasisFunction(RBF_VARIANCE, RBF_LENGTHSCALE) + kernel = kernels.AsymptoticKernel(rbf) + kernel.add_asymptotics(region="UV", asymptotics=uv_shape) + kernel.set_params(asymp_params=[MU_UV, L_UV]) + return kernel + + +def synthetic_data(): + """Momenta and correlator values of the synthetic data set.""" + p = np.geomspace(P_MIN, P_MAX, N_DATA) + return p, p ** 2 / (p ** 2 + 1.0) + + +def build_model(with_tail=True, gamma=GAMMA_GLUON, uv_shape=propagator_uv_shape): + """Two-constraint model: a KL data constraint plus a one-row sum rule. + + ``with_tail=False`` gives both constraints a plain ``GaussLegendre_1D_log`` on the same grid, + which is what the tgvert/ghgvert/fgvertsd projects do and what T6 pins down. + """ + p, G = synthetic_data() + kernel = make_kernel(uv_shape) + + data_integrator = integrators.GaussLegendre_1D_log(W_MIN, W_UV, INT_N) + c_data = constraints.LinearEquality( + operators.Integral(kl_p2, data_integrator), + {"x": p, "y": G, "cov_y": DATA_COV * np.ones_like(G)}, + ) + + if with_tail: + sr_integrator = integrators.GaussLegendre_1D_log_UVtail( + W_MIN, W_UV, INT_N, uv_shape, tail_moment=sum_rule_tail_moment(W_UV, gamma) + ) + else: + sr_integrator = integrators.GaussLegendre_1D_log(W_MIN, W_UV, INT_N) + c_sr = constraints.LinearEquality( + operators.Integral(sr_kernel, sr_integrator), + {"x": np.array([0.0]), "y": np.array([0.0]), "cov_y": np.array([SR_COV])}, + ) + + model = models.GaussianProcess(kernel, [c_data, c_sr]) + return model, c_data, c_sr, data_integrator, sr_integrator + + +def bulk_double_integration(integrator, c1, kernel, c2): + """Uncorrected double quadrature BB, written out so the tests do not lean on the fix.""" + row1 = integrator.weights * c1(make_row_vector(integrator.w), x=make_column_vector(c1.x)) + row2 = integrator.weights * c2(make_row_vector(integrator.w), x=make_column_vector(c2.x)) + return row1 @ kernel(integrator.w, integrator.w) @ row2.T + + +# -------------------------------------------------------------------------------------------- +# T1 -- the sharpest test: pure linear algebra, no linear solve, so it must hold to machine +# precision. The cross block of OpKerOp is, by definition, the data-side reintegration of the +# sum-rule row of the prediction operator OpKer. That identity is exactly what the missing +# cross-block correction used to break. +# -------------------------------------------------------------------------------------------- + +def test_cross_block_equals_reintegrated_opker_row(): + model, c_data, c_sr, data_integrator, _ = build_model() + kernel = model.kernel + n_data = c_data.x.shape[0] + + opkerop = TwoSided()(kernel, [c_data, c_sr]) + cross_block = opkerop[:n_data, n_data:] # (M, 1) + + nodes = data_integrator.w + opker = OneSided()(kernel, [c_data, c_sr], nodes) # (M + 1, N) + sr_row = opker[n_data:, :] # (1, N) + + data_weighted = data_integrator.weights * kl_p2( + make_column_vector(c_data.x), make_row_vector(nodes) + ) + reintegrated = data_weighted @ sr_row.T # (M, 1) + + np.testing.assert_allclose(cross_block, reintegrated, rtol=1e-12) + + +# -------------------------------------------------------------------------------------------- +# T2 -- np.linalg.cholesky in models.py reads only the lower triangle, so an asymmetric OpKerOp +# would never raise. This is the only guard against that. +# -------------------------------------------------------------------------------------------- + +def test_covariance_matrix_is_symmetric(): + model, c_data, c_sr, _, _ = build_model() + kernel = model.kernel + n_data = c_data.x.shape[0] + + opkerop = TwoSided()(kernel, [c_data, c_sr]) + scale = np.abs(opkerop).max() + np.testing.assert_allclose(opkerop, opkerop.T, rtol=1e-14, atol=1e-14 * scale) + + # Building the same model with the constraint order reversed must give the permuted matrix, + # i.e. block (sr, data) must carry exactly the same correction as block (data, sr). + reversed_opkerop = TwoSided()(kernel, [c_sr, c_data]) + perm = np.concatenate([np.arange(n_data, n_data + 1), np.arange(n_data)]) + np.testing.assert_allclose( + reversed_opkerop, opkerop[np.ix_(perm, perm)], rtol=1e-14, atol=1e-14 * scale + ) + + +# -------------------------------------------------------------------------------------------- +# T3 -- the UVtail subclass no longer overrides doubleIntegrationSymmetric. This locks in that +# the inherited path reproduces the deleted override's closed form bulk + 2 A T + T^2 for a +# multi-row constraint with a row-independent tail moment. +# -------------------------------------------------------------------------------------------- + +def test_symmetric_block_matches_legacy_closed_form(): + kernel = make_kernel() + tail_value = analytic_sum_rule_tail(W_UV, GAMMA_GLUON) + integrator = integrators.GaussLegendre_1D_log_UVtail( + W_MIN, W_UV, INT_N, propagator_uv_shape, + tail_moment=sum_rule_tail_moment(W_UV, GAMMA_GLUON), + ) + p = np.array([0.7, 3.0, 11.0]) # Three rows, so a broadcast-vs-outer-product bug would show. + constraint = constraints.LinearEquality( + operators.Integral(kl_p2, integrator), + {"x": p, "y": np.zeros_like(p), "cov_y": DATA_COV * np.ones_like(p)}, + ) + + symmetric = integrator.doubleIntegrationSymmetric(constraint, kernel) + general = integrator.doubleIntegration(constraint, kernel, constraint) + np.testing.assert_allclose(symmetric, general, rtol=1e-12) + + bulk = bulk_double_integration(integrator, constraint, kernel, constraint) + c_row = constraint(make_row_vector(integrator.w), x=make_column_vector(constraint.x)) + A = (integrator.weights * c_row @ kernel(integrator.w, integrator.w_uv)) / integrator.f_uv_anchor + ones_col = np.ones((A.shape[0], 1)) + legacy = bulk + tail_value * (A @ ones_col.T + ones_col @ A.T) + tail_value ** 2 + + np.testing.assert_allclose(symmetric, legacy, rtol=1e-12) + + +# -------------------------------------------------------------------------------------------- +# T4 -- the analytic value itself, for both anomalous dimensions, plus a quadrature cross-check +# that extending the numerical range towards infinity really does converge onto +# bulk(w_uv) + T_sr. Convergence in log(w) is slow, hence the loose few-percent tolerance. +# -------------------------------------------------------------------------------------------- + +@pytest.mark.parametrize("gamma, uv_shape, max_truncation_error", [ + # max_truncation_error: how far the integral truncated at 10^20 still is from bulk + T. The + # tail decays only like (2 log w)^-gamma, so the smaller gamma the slower the convergence -- + # ~1% for the gluon, ~18% for the ghost. That slowness is exactly why the analytic tail is + # mandatory and extending the quadrature range is not a viable alternative. + (GAMMA_GLUON, propagator_uv_shape, 0.05), + (GAMMA_GHOST, ghost_uv_shape, 0.20), +]) +def test_tail_moment_matches_analytic_sum_rule_value(gamma, uv_shape, max_truncation_error): + integrator = integrators.GaussLegendre_1D_log_UVtail( + W_MIN, W_UV, INT_N, uv_shape, tail_moment=sum_rule_tail_moment(W_UV, gamma) + ) + constraint = constraints.LinearEquality( + operators.Integral(sr_kernel, integrator), + {"x": np.array([0.0]), "y": np.array([0.0]), "cov_y": np.array([SR_COV])}, + ) + tail = integrator.uv_tail_moment(constraint) + assert tail.shape == (1, 1) + np.testing.assert_allclose(tail.item(), analytic_sum_rule_tail(W_UV, gamma), rtol=1e-14) + + # Cross-check: sum_rule integral over (W_MIN, 10^k) -> integral over (W_MIN, w_uv) + T. + # n_nodes=2000: the two truncations use different grids, so their shared bulk parts only + # cancel to quadrature accuracy; 2000 nodes brings that below 4e-6 relative (800 leaves 2e-2). + def numeric_sum_rule_integral(w_max, n_nodes=2000): + integ = integrators.GaussLegendre_1D_log(W_MIN, w_max, n_nodes) + return float((integ.weights @ (integ.w * uv_shape(integ.w))).item()) + + target = numeric_sum_rule_integral(W_UV) + analytic_sum_rule_tail(W_UV, gamma) + shortfalls = [target - numeric_sum_rule_integral(10.0 ** k) for k in (6, 8, 12, 20)] + + # Whatever the truncated quadrature misses must be precisely the closed form's own remainder + # above the truncation point -- the direct check of the closed form against quadrature of the + # actual (regularised) f_uv, independent of how slowly the tail converges. + for k, shortfall in zip((6, 8, 12, 20), shortfalls): + np.testing.assert_allclose(shortfall, analytic_sum_rule_tail(10.0 ** k, gamma), rtol=1e-5) + + errors = [abs(s) / abs(target) for s in shortfalls] + assert errors == sorted(errors, reverse=True), \ + f"truncated integral should approach bulk + T monotonically, got {errors}" + assert errors[-1] < max_truncation_error, \ + f"k=20 still {errors[-1]:.3%} away from bulk + T" + + +# -------------------------------------------------------------------------------------------- +# T5 -- end-to-end self-consistency of the posterior. predict_data() and predict() must describe +# the same model: reintegrating the posterior mean over fredipy's own nodes and weights has to +# reproduce predict_data(). This is the identity the bug broke by ~1e0. +# -------------------------------------------------------------------------------------------- + +def test_predict_data_matches_reintegrated_posterior_mean(): + model, c_data, c_sr, data_integrator, sr_integrator = build_model() + n_data = c_data.x.shape[0] + + nodes = data_integrator.w + rho, _ = model.predict(nodes) # (N, 1) + mu, _ = model.predict_data() # (M + 1, 1) + + data_weighted = data_integrator.weights * kl_p2( + make_column_vector(c_data.x), make_row_vector(nodes) + ) + np.testing.assert_allclose(data_weighted @ rho, mu[:n_data], rtol=1e-9) + + # Sum-rule row: the reintegration only covers the bulk, so the analytic tail T * A_uv has to + # be added back, with A_uv = rho(w_uv) / f_uv(w_uv) the UV amplitude of the posterior mean. + sr_weighted = sr_integrator.weights * sr_kernel( + make_column_vector(c_sr.x), make_row_vector(nodes) + ) + rho_at_w_uv, _ = model.predict(sr_integrator.w_uv) + tail_value = sr_integrator.uv_tail_moment(c_sr).item() + reintegrated_sr = (sr_weighted @ rho).item() + reintegrated_sr += tail_value * rho_at_w_uv.item() / sr_integrator.f_uv_anchor + + np.testing.assert_allclose(reintegrated_sr, mu[n_data:].item(), rtol=1e-4) + + +# -------------------------------------------------------------------------------------------- +# T6 -- reference values captured with the pre-fix code, before the tail hook existed. A model +# whose constraints all carry plain GaussLegendre_1D_log integrators must stay bit-identical: +# this is what protects three-gluon-vert, ghost-gluon-vert and fgvert-sd from the change. +# -------------------------------------------------------------------------------------------- + +REFERENCE_LOG_LIKELIHOOD = -278560.62089715304 +REFERENCE_PREDICT_DATA = np.array([ + 0.11838105684728362, 0.2086624375442625, 0.35474911867640913, 0.5465432399942074, + 0.75084480180521496, 0.92431427465635352, 1.034895190037787, 1.0745688577007968, + 1.0557558063883334, 0.99917723311227746, 0.92389010410988703, 0.8432306828617584, + 0.76490387672674842, 0.69258282921509817, 0.62751788803143427, 0.0017991249333135784, +]) +REFERENCE_PREDICT_MEAN = np.array([ + 32.848601571109612, -11.554450084397104, -0.77481442881980911, + -0.00017162581324248194, -1.0987647334001352e-09, +]) +REFERENCE_PREDICT_GRID = np.array([0.01, 0.5, 2.0, 50.0, 1e4]) # Prediction points of the reference. +REFERENCE_OPKEROP_TRACE = 10.170139984656482 + + +def test_plain_integrators_are_unchanged_by_the_tail_hook(): + model, _, _, _, _ = build_model(with_tail=False) + + np.testing.assert_allclose(model.log_likelihood(), REFERENCE_LOG_LIKELIHOOD, rtol=1e-15) + mu, _ = model.predict_data() + np.testing.assert_allclose(mu.ravel(), REFERENCE_PREDICT_DATA, rtol=1e-15) + mean, _ = model.predict(REFERENCE_PREDICT_GRID) + np.testing.assert_allclose(mean.ravel(), REFERENCE_PREDICT_MEAN, rtol=1e-15) + trace = np.trace(model._posterior_cache["OpKerOp"]) + np.testing.assert_allclose(trace, REFERENCE_OPKEROP_TRACE, rtol=1e-15) + + +# -------------------------------------------------------------------------------------------- +# T7 -- (**) mixes both constraints on a shared quadrature grid, so it is undefined when the two +# integrators discretise omega differently. Fail loudly rather than silently mixing grids. +# -------------------------------------------------------------------------------------------- + +def test_mismatched_integrator_grids_raise(): + kernel = make_kernel() + p, G = synthetic_data() + + data_integrator = integrators.GaussLegendre_1D_log(W_MIN, 1e4, INT_N + 30) # Different grid. + c_data = constraints.LinearEquality( + operators.Integral(kl_p2, data_integrator), + {"x": p, "y": G, "cov_y": DATA_COV * np.ones_like(G)}, + ) + sr_integrator = integrators.GaussLegendre_1D_log_UVtail( + W_MIN, W_UV, INT_N, propagator_uv_shape, + tail_moment=sum_rule_tail_moment(W_UV, GAMMA_GLUON), + ) + c_sr = constraints.LinearEquality( + operators.Integral(sr_kernel, sr_integrator), + {"x": np.array([0.0]), "y": np.array([0.0]), "cov_y": np.array([SR_COV])}, + ) + + with pytest.raises(NotImplementedError) as excinfo: + TwoSided()(kernel, [c_data, c_sr]) + message = str(excinfo.value) + assert "10000.0" in message and "100000.0" in message, message + assert str(INT_N) in message and str(INT_N + 30) in message, message + + +# -------------------------------------------------------------------------------------------- +# T8 -- singleIntegration's tail is T[m] * K(w_uv, w_pred) / f_uv(w_uv), an outer product. The +# pre-fix code broadcast a single scalar, which happened to be right only because the sum-rule +# constraint has exactly one row. +# -------------------------------------------------------------------------------------------- + +def test_single_integration_tail_is_row_wise_outer_product(): + kernel = make_kernel() + row_moments = np.array([[0.25], [0.5], [1.75]]) # Deliberately row-dependent tail moments. + integrator = integrators.GaussLegendre_1D_log_UVtail( + W_MIN, W_UV, INT_N, propagator_uv_shape, + tail_moment=lambda x: row_moments.copy(), + ) + p = np.array([0.7, 3.0, 11.0]) + constraint = constraints.LinearEquality( + operators.Integral(kl_p2, integrator), + {"x": p, "y": np.zeros_like(p), "cov_y": DATA_COV * np.ones_like(p)}, + ) + w_pred = make_column_vector(np.geomspace(0.05, 500.0, 7)) + + full = integrator.singleIntegration(constraint, kernel, w_pred) + assert full.shape == (3, w_pred.shape[0]) + + bulk = ( + integrator.weights * constraint(make_row_vector(integrator.w), x=make_column_vector(p)) + @ kernel(integrator.w, w_pred) + ) + anchor_row = kernel(integrator.w_uv, w_pred) / integrator.f_uv_anchor # (1, N_pred) + for m in range(3): + expected = bulk[m] + row_moments[m, 0] * anchor_row.ravel() + np.testing.assert_allclose(full[m], expected, rtol=1e-12) + + +# -------------------------------------------------------------------------------------------- +# T9 -- the analytic NLL gradient against finite differences of the NLL itself. +# +# log_likelihood_grad() reassembles OpKerOp with dK/dtheta substituted for K. The A terms of +# (**) are linear in the kernel and differentiate themselves correctly, but T_1 T_2^T contains +# no kernel at all -- T is the analytic tail moment, a function of w_uv and gamma only -- so its +# hyperparameter derivative is zero. Adding it unchanged (which is what the code did before the +# derivative flag existed) put the gradient off by factors of 34 to 1675 on this very model, with +# the wrong sign on the RBF lengthscale. This test is the acceptance criterion for that fix: no +# other test in either suite can see a wrong gradient, since the likelihood *value* is correct +# either way. +# -------------------------------------------------------------------------------------------- + +FD_REL_STEP = 1e-4 # Relative central-difference step; near the optimum of truncation vs roundoff. +FD_RTOL = 1e-5 # Agreement demanded per component; the FD truncation floor here is ~5e-7. + +# Hyperparameter vector [rbf_variance, rbf_lengthscale, mu_uv, l_uv] at which the check is run. +FD_PARAMS = [RBF_VARIANCE, RBF_LENGTHSCALE, MU_UV, L_UV] + + +def _set_flat_params(model, params): + """Push a flat [rbf_variance, rbf_lengthscale, mu_uv, l_uv] vector into the model.""" + model.reset() + model.kernel.set_params(kernel_params=list(params[:2]), asymp_params=list(params[2:])) + + +def _finite_difference_gradient(model, params, rel_step=FD_REL_STEP): + """Central-difference gradient of log_likelihood() w.r.t. every entry of ``params``.""" + grad = [] + for i in range(len(params)): + h = rel_step * abs(params[i]) + plus = list(params) + plus[i] += h + _set_flat_params(model, plus) + f_plus = model.log_likelihood() + minus = list(params) + minus[i] -= h + _set_flat_params(model, minus) + f_minus = model.log_likelihood() + grad.append((f_plus - f_minus) / (2.0 * h)) + _set_flat_params(model, params) + return np.array(grad) + + +@pytest.mark.parametrize("with_tail", [True, False]) +def test_log_likelihood_grad_matches_finite_differences(with_tail): + """Every active hyperparameter, with and without the UVtail sum-rule constraint.""" + model, _, _, _, _ = build_model(with_tail=with_tail) + assert model.kernel.dim == len(FD_PARAMS) + + _set_flat_params(model, FD_PARAMS) + analytic = np.array(model.log_likelihood_grad(), dtype=float) + numeric = _finite_difference_gradient(model, FD_PARAMS) + + np.testing.assert_allclose(analytic, numeric, rtol=FD_RTOL) + + +def test_gradient_flag_is_a_noop_for_plain_integrators(): + """The derivative flag must not perturb the plain-quadrature gradient path at all. + + Plain ``GaussLegendre_1D_log`` blocks are linear in the kernel, so the assembled gradient + covariance has to be bit-identical with and without the flag. This is what guarantees that + three-gluon-vert, ghost-gluon-vert and fgvert-sd optimize against exactly the gradients they + did before the flag was introduced. + """ + model, _, _, _, _ = build_model(with_tail=False) + _set_flat_params(model, FD_PARAMS) + for kernel_grad in model.kernel.params_gradient(): + with_flag = TwoSided()(kernel_grad, model.constraints, derivative=True) + without_flag = TwoSided()(kernel_grad, model.constraints, derivative=False) + assert np.array_equal(with_flag, without_flag) + + +def test_gradient_flag_drops_only_the_constant_tail_term(): + """With a UVtail constraint the flag must remove exactly T_1 T_2^T and nothing else.""" + model, c_data, c_sr, _, sr_integrator = build_model(with_tail=True) + _set_flat_params(model, FD_PARAMS) + n_data = c_data.x.shape[0] + tail = sr_integrator.uv_tail_moment(c_sr) # (1, 1) + + for kernel_grad in model.kernel.params_gradient(): + with_flag = TwoSided()(kernel_grad, model.constraints, derivative=True) + without_flag = TwoSided()(kernel_grad, model.constraints, derivative=False) + difference = without_flag - with_flag + + # Only the sum-rule diagonal block differs, by exactly T^2: the data constraint carries a + # plain integrator, hence a zero tail moment, hence no constant term in any block it enters. + expected = np.zeros_like(difference) + expected[n_data:, n_data:] = tail @ tail.T + # rtol 1e-12: the sum-rule block is a difference of two O(1) sums, so it carries a few + # ulp of cancellation. atol 0: every other block must be bitwise unchanged by the flag. + np.testing.assert_allclose(difference, expected, rtol=1e-12, atol=0.0) From b1c8fd2e8bfd3431f9bf496b76128290e858d11c Mon Sep 17 00:00:00 2001 From: Jonas Wessely Date: Fri, 18 Sep 2026 10:00:40 +0200 Subject: [PATCH 6/8] Add analytic gradients and derivatives to Matern32/Matern52 Adds params_gradient() (variance and per-lengthscale partials) plus dK_dx/dK_dy/d2K_dxdy to both Matern32 and Matern52, giving them full derivative parity with RadialBasisFunction. This is load-bearing for optimizer support: local.py's L-BFGS-B step calls with jac=True unconditionally and objective.py raises if params_gradient is missing, so neither Matern kernel could be used in any optimizer run before this. Formulas per docs/matern_kernel_plan.md section 4 (Phase 1), derived and cross-checked analytically against central finite differences: - dK/dvariance = 2K/v for both kernels. - dK/dlengthscale_i = 3 v^2 e^{-z} D_i/l_i^3 (Matern32) and (5/3) v^2 (1+z) e^{-z} D_i/l_i^3 (Matern52). - dK_dx/dK_dy/d2K_dxdy use the same 1-D convention (lengthscale[0], broadcast) as RadialBasisFunction's existing implementations. New tests in tests/test_kernels.py, modeled on test_rbf_params_gradient: test_matern{32,52}_params_gradient (central finite difference, eps=1e-5, rtol=1e-4, parametrized 1-D/2-D) and test_matern{32,52}_derivatives (dK_dx/dK_dy vs. finite difference of make(), d2K_dxdy vs. second-order finite difference, eps=1e-4 per the plan's noted second-derivative noise floor). Matern-1/2 is intentionally excluded (plan section 2): its covariance has a kink at r=0, so d2K_dxdy does not exist as a finite function there. fredipy/tests: 89 passed (baseline 83 + 6 new test cases). Co-Authored-By: Claude Sonnet 5 --- fredipy/kernels.py | 167 ++++++++++++++++++++++++++++++++++++++++++ tests/test_kernels.py | 151 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 318 insertions(+) diff --git a/fredipy/kernels.py b/fredipy/kernels.py index 638fc98..aa49cde 100644 --- a/fredipy/kernels.py +++ b/fredipy/kernels.py @@ -613,6 +613,88 @@ def make( self._K = _K return _K + def params_gradient(self) -> List[Callable]: + """Returns the gradient of the kernel in the direction of the parameters. + Currently no implementation for the gradient when having derivative constraints. + + Returns + ------- + grad : list of the gradient in the directions [variance, lengthscale1, lengthscale2, ...] + """ + def empy(x, y): None + grad = [empy for _ in range(self.dim)] + grad[0] = lambda x, y: 2 / self.variance * self.make(x, y) + for i in range(self.dim - 1): + grad[i + 1] = lambda x, y, i=i: 3 * self.variance**2 * np.exp( + -np.sqrt(3 * self._sqdist(x / self.lengthscale, y / self.lengthscale)) + ) * ( + np.sum(make_column_vector(x[:, i])**2, 1)[:, None] + + np.sum(make_column_vector(y[:, i])**2, 1) + - 2 * make_column_vector(x[:, i]) @ make_row_vector(y[:, i]) + ) / self.lengthscale[i]**3 + return grad + + def dK_dx( + self, + x: np.ndarray, + y: np.ndarray + ) -> np.ndarray: + """Derivative of the kernel with respect to the FIRST argument + + Parameters + ---------- + x : first argument of the kernel, numpy array + y : second argument of the kernel, numpy array + + Returns + ------- + Matrix of the kernel derivative, numpy array + """ + _x = x / self.lengthscale + _y = y / self.lengthscale + dist = np.sqrt(3 * self._sqdist(_x, _y)) + return -3 * self.variance**2 * np.exp(-dist) * \ + (_x / self.lengthscale - (_y / self.lengthscale).T) + + def dK_dy( + self, + x: np.ndarray, + y: np.ndarray + ) -> np.ndarray: + """Derivative of the kernel with respect to the SECOND argument + + Parameters + ---------- + x : first argument of the kernel, numpy array + y : second argument of the kernel, numpy array + + Returns + ------- + Matrix of the kernel derivative, numpy array + """ + return -1 * self.dK_dx(x, y) + + def d2K_dxdy( + self, + x: np.ndarray, + y: np.ndarray + ) -> np.ndarray: + """ + Derivative of the kernel with respect to both arguments + Parameters + ---------- + x : first argument of the kernel, numpy array + y : second argument of the kernel, numpy array + + Returns + ------- + Matrix of the kernel derivative, numpy array + """ + _x = x / self.lengthscale + _y = y / self.lengthscale + dist = np.sqrt(3 * self._sqdist(_x, _y)) + return 3 * self.variance**2 * (1 - dist) * np.exp(-dist) / self.lengthscale[0]**2 + def set_params( self, params: List @@ -675,6 +757,91 @@ def make( self._K = _K return _K + def params_gradient(self) -> List[Callable]: + """Returns the gradient of the kernel in the direction of the parameters. + Currently no implementation for the gradient when having derivative constraints. + + Returns + ------- + grad : list of the gradient in the directions [variance, lengthscale1, lengthscale2, ...] + """ + def empy(x, y): None + grad = [empy for _ in range(self.dim)] + grad[0] = lambda x, y: 2 / self.variance * self.make(x, y) + for i in range(self.dim - 1): + grad[i + 1] = lambda x, y, i=i: (5. / 3.) * self.variance**2 * ( + 1 + np.sqrt(5 * self._sqdist(x / self.lengthscale, y / self.lengthscale)) + ) * np.exp( + -np.sqrt(5 * self._sqdist(x / self.lengthscale, y / self.lengthscale)) + ) * ( + np.sum(make_column_vector(x[:, i])**2, 1)[:, None] + + np.sum(make_column_vector(y[:, i])**2, 1) + - 2 * make_column_vector(x[:, i]) @ make_row_vector(y[:, i]) + ) / self.lengthscale[i]**3 + return grad + + def dK_dx( + self, + x: np.ndarray, + y: np.ndarray + ) -> np.ndarray: + """Derivative of the kernel with respect to the FIRST argument + + Parameters + ---------- + x : first argument of the kernel, numpy array + y : second argument of the kernel, numpy array + + Returns + ------- + Matrix of the kernel derivative, numpy array + """ + _x = x / self.lengthscale + _y = y / self.lengthscale + dist = np.sqrt(5 * self._sqdist(_x, _y)) + return -(5. / 3.) * self.variance**2 * (1 + dist) * np.exp(-dist) * \ + (_x / self.lengthscale - (_y / self.lengthscale).T) + + def dK_dy( + self, + x: np.ndarray, + y: np.ndarray + ) -> np.ndarray: + """Derivative of the kernel with respect to the SECOND argument + + Parameters + ---------- + x : first argument of the kernel, numpy array + y : second argument of the kernel, numpy array + + Returns + ------- + Matrix of the kernel derivative, numpy array + """ + return -1 * self.dK_dx(x, y) + + def d2K_dxdy( + self, + x: np.ndarray, + y: np.ndarray + ) -> np.ndarray: + """ + Derivative of the kernel with respect to both arguments + Parameters + ---------- + x : first argument of the kernel, numpy array + y : second argument of the kernel, numpy array + + Returns + ------- + Matrix of the kernel derivative, numpy array + """ + _x = x / self.lengthscale + _y = y / self.lengthscale + dist = np.sqrt(5 * self._sqdist(_x, _y)) + return (5. / 3.) * self.variance**2 * (1 + dist - dist**2) * np.exp(-dist) \ + / self.lengthscale[0]**2 + def set_params( self, params: List diff --git a/tests/test_kernels.py b/tests/test_kernels.py index a886dd8..b8786b9 100644 --- a/tests/test_kernels.py +++ b/tests/test_kernels.py @@ -332,3 +332,154 @@ def test_rbf_params_gradient( # restore kernel.set_params(params) + + +@pytest.mark.parametrize("variance, lengthscale", [ + [2.0, 1.5], # 1D + [2.0, [0.5, 1.5]], # 2D -- catches late-binding closure bugs + ]) +def test_matern32_params_gradient( + variance: float, + lengthscale: float | list[float], + ) -> None: + """Check analytical gradient against finite differences for each hyperparameter. + + Cloned from test_rbf_params_gradient's structure (plan docs/matern_kernel_plan.md section 5). + """ + kernel = kernels.Matern32(variance=variance, lengthscale=lengthscale) + n_dim = len(kernel.lengthscale) + n_test = 5 + X = rng.randn(n_test, n_dim) + Y = rng.randn(n_test, n_dim) + + params = [kernel.variance] + list(kernel.lengthscale) + + # collect analytical gradients before mutating kernel state + grad = kernel.params_gradient() + analytical = [grad[idx](X, Y) for idx in range(kernel.dim)] + + eps = 1e-5 # Central finite-difference step for the hyperparameter gradient check (plan section 5). + for idx in range(kernel.dim): + # central finite difference + params_plus = list(params) + params_plus[idx] += eps + params_minus = list(params) + params_minus[idx] -= eps + + kernel.set_params(params_plus) + kernel._empty_cache() + K_plus = kernel.make(X, Y, cache=False) + kernel.set_params(params_minus) + kernel._empty_cache() + K_minus = kernel.make(X, Y, cache=False) + + numerical = (K_plus - K_minus) / (2 * eps) + # rtol per plan section 5's Phase 1 bullet. + assert_allclose(analytical[idx], numerical, rtol=1e-4) + + # restore + kernel.set_params(params) + + +@pytest.mark.parametrize("variance, lengthscale", [ + [2.0, 1.5], # 1D + [2.0, [0.5, 1.5]], # 2D -- catches late-binding closure bugs + ]) +def test_matern52_params_gradient( + variance: float, + lengthscale: float | list[float], + ) -> None: + """Check analytical gradient against finite differences for each hyperparameter. + + Cloned from test_rbf_params_gradient's structure (plan docs/matern_kernel_plan.md section 5). + """ + kernel = kernels.Matern52(variance=variance, lengthscale=lengthscale) + n_dim = len(kernel.lengthscale) + n_test = 5 + X = rng.randn(n_test, n_dim) + Y = rng.randn(n_test, n_dim) + + params = [kernel.variance] + list(kernel.lengthscale) + + # collect analytical gradients before mutating kernel state + grad = kernel.params_gradient() + analytical = [grad[idx](X, Y) for idx in range(kernel.dim)] + + eps = 1e-5 # Central finite-difference step for the hyperparameter gradient check (plan section 5). + for idx in range(kernel.dim): + # central finite difference + params_plus = list(params) + params_plus[idx] += eps + params_minus = list(params) + params_minus[idx] -= eps + + kernel.set_params(params_plus) + kernel._empty_cache() + K_plus = kernel.make(X, Y, cache=False) + kernel.set_params(params_minus) + kernel._empty_cache() + K_minus = kernel.make(X, Y, cache=False) + + numerical = (K_plus - K_minus) / (2 * eps) + # rtol per plan section 5's Phase 1 bullet. + assert_allclose(analytical[idx], numerical, rtol=1e-4) + + # restore + kernel.set_params(params) + + +def _check_matern_derivatives( + kernel_cls, + variance: float, + lengthscale: float, + ) -> None: + """Shared body for test_matern{32,52}_derivatives below (not collected by pytest itself, + since its name does not start with 'test_'; called directly by each named test). + + Checks dK_dx/dK_dy against finite differences of make(), and d2K_dxdy against a + second-order finite difference, per plan docs/matern_kernel_plan.md section 5. + + dK_dx/dK_dy/d2K_dxdy are only defined in the 1-D convention (matching + RadialBasisFunction's own dK_dx/dK_dy, which likewise use lengthscale[0] and broadcast; + see plan section 4), so this check is 1-D only. + """ + kernel = kernel_cls(variance=variance, lengthscale=lengthscale) + n_test = 6 + X = rng.randn(n_test, 1) + Y = rng.randn(n_test, 1) + + eps = 1e-4 # Finite-difference step; 1e-6 gives ~1e-3 noise on the 2nd derivative (plan section 5). + + analytical_dK_dx = kernel.dK_dx(X, Y) + analytical_dK_dy = kernel.dK_dy(X, Y) + analytical_d2K_dxdy = kernel.d2K_dxdy(X, Y) + + # dK_dx: central finite difference of make() in the first argument. + K_x_plus = kernel.make(X + eps, Y, cache=False) + K_x_minus = kernel.make(X - eps, Y, cache=False) + numerical_dK_dx = (K_x_plus - K_x_minus) / (2 * eps) + assert_allclose(analytical_dK_dx, numerical_dK_dx, rtol=1e-4, atol=1e-6) + + # dK_dy: central finite difference of make() in the second argument. + K_y_plus = kernel.make(X, Y + eps, cache=False) + K_y_minus = kernel.make(X, Y - eps, cache=False) + numerical_dK_dy = (K_y_plus - K_y_minus) / (2 * eps) + assert_allclose(analytical_dK_dy, numerical_dK_dy, rtol=1e-4, atol=1e-6) + + # d2K_dxdy: second-order (mixed) finite difference of make() in both arguments. + K_pp = kernel.make(X + eps, Y + eps, cache=False) + K_pm = kernel.make(X + eps, Y - eps, cache=False) + K_mp = kernel.make(X - eps, Y + eps, cache=False) + K_mm = kernel.make(X - eps, Y - eps, cache=False) + numerical_d2K_dxdy = (K_pp - K_pm - K_mp + K_mm) / (4 * eps**2) + assert_allclose(analytical_d2K_dxdy, numerical_d2K_dxdy, rtol=1e-3, atol=1e-3) + + +def test_matern32_derivatives() -> None: + """dK_dx/dK_dy/d2K_dxdy vs. finite differences for Matern-3/2 (plan section 5).""" + _check_matern_derivatives(kernels.Matern32, 2.0, 1.5) + + +def test_matern52_derivatives() -> None: + """dK_dx/dK_dy/d2K_dxdy vs. finite differences for Matern-5/2 (plan section 5).""" + _check_matern_derivatives(kernels.Matern52, 2.0, 1.5) From bc986fb56db5c8118d12df7c7701400a807f7d88 Mon Sep 17 00:00:00 2001 From: Jonas Wessely <104987725+Jpwessely@users.noreply.github.com> Date: Tue, 22 Sep 2026 13:44:45 +0200 Subject: [PATCH 7/8] AsymptoticKernel: selectable blend ('logistic' | 'hill'), Hill rising weights In blend='hill' the rising weights (sigma_+^IR on the base kernel, sigma_+^UV on the UV branch) become a logistic in ln(w), 1/(1+(mu/w)^(1/l)), vanishing like (w/mu)^(1/l) at w -> 0; falling weights stay the linear-w logistic. 'logistic' is unchanged bit for bit. Adds hill_softtheta, gradient formulas and finite-difference tests. Co-Authored-By: Claude Opus 5 --- fredipy/integrators.py | 8 +- fredipy/kernels.py | 241 +++++++++++++++++++++++++++++------------ fredipy/util.py | 57 ++++++++++ tests/test_kernels.py | 154 ++++++++++++++++++++++++++ tests/test_util.py | 75 ++++++++++++- 5 files changed, 464 insertions(+), 71 deletions(-) diff --git a/fredipy/integrators.py b/fredipy/integrators.py index 382201d..b630382 100644 --- a/fredipy/integrators.py +++ b/fredipy/integrators.py @@ -517,7 +517,13 @@ class GaussLegendre_1D_log_UVtail(GaussLegendre_1D_log): K(ω, ω') = [θ_uv(ω) f_uv(ω)] · f_uv(ω') (*) This is the *only* assumption made here; it holds up to - ``O(exp(-(w_uv - mu_uv)/l_uv))``. + ``O(exp(-(w_uv - mu_uv)/l_uv))`` under ``AsymptoticKernel``'s default ``blend="logistic"``. + Under ``blend="hill"`` (docs/rewrite_plan.md section 4, "selectable Hill blend") theta_uv + itself is the one *rising* weight this bound depends on, and it saturates as a power law + instead, ``O((mu_uv/w_uv)**(1/l_uv))`` -- a slower asymptotic order, but still utterly + negligible at the ``w_uv >> mu_uv`` separations this repo actually uses (e.g. gluon: + ``mu_uv ~ 1``, ``w_uv ~ 1e5``, ``1/l_uv ~ 5`` gives ``~1e-25``), so no change to this + integrator or its callers was needed for the Hill blend. **General two-constraint correction.** Splitting each of the two integrals of a covariance block at w_uv into bulk (w_min, w_uv) and tail (w_uv, ∞) and diff --git a/fredipy/kernels.py b/fredipy/kernels.py index aa49cde..f39d501 100644 --- a/fredipy/kernels.py +++ b/fredipy/kernels.py @@ -1,7 +1,17 @@ from typing import List, Callable import numpy as np -from .util import allclose, make_column_vector, make_row_vector, softtheta +from .util import allclose, make_column_vector, make_row_vector, softtheta, hill_softtheta + + +# The complete set of legal AsymptoticKernel ``blend`` values (docs/rewrite_plan.md section 4, +# "Taken 2026-09-22 (gluon coworking: data_cov scan, selectable Hill blend)"). "logistic" is the +# kernel's original, always-linear-in-x softtheta blend, kept bit-identical. "hill" switches only +# the two *rising* transition weights (the base-kernel weight above mu_ir, the UV-branch weight +# above mu_uv) to a logistic in ln(x) (hill_softtheta); the two *falling* weights (the IR-branch +# weight, the base-kernel weight below mu_uv) stay the linear-x softtheta in both modes, so every +# large-x cutoff stays exponential regardless of blend. +BLEND_SHAPES: tuple[str, ...] = ("logistic", "hill") class Kernel: @@ -257,14 +267,32 @@ class AsymptoticKernel(Kernel): """Kernel with optional UV and IR asymptotics The asymptotics are only applied in the first (omega) dimension. + + Parameters + ---------- + kernel : the 'non-asymptotic' base kernel. + blend : one of ``BLEND_SHAPES`` ("logistic", "hill"), selecting the functional form of the + transition weights between regions. "logistic" (default) is the original behaviour: every + transition weight is ``softtheta``, a logistic in x itself. "hill" switches only the two + *rising* weights -- the base-kernel weight above mu_ir, the UV-branch weight above mu_uv -- + to ``hill_softtheta``, a logistic in ln(x); the two *falling* weights (the IR-branch weight, + the base-kernel weight below mu_uv) are left as softtheta in both modes, so every large-x + cutoff stays exponential regardless of blend (docs/rewrite_plan.md section 4, "Taken + 2026-09-22 (gluon coworking ... selectable Hill blend)"). "logistic" is bit-identical to + the kernel's pre-blend-selector behaviour. """ def __init__( self, kernel: Kernel, + blend: str = "logistic", ) -> None: self.kernel = kernel # the 'non-asymptotic' kernel part + if blend not in BLEND_SHAPES: + raise ValueError(f"blend must be one of {BLEND_SHAPES}, got {blend!r}") + self.blend = blend # transition-weight functional form, see class docstring + self.uv_asymptotics = lambda x: np.zeros_like(x) # some dummy function so we avoid errors self.mu_uv = 1 self.l_uv = 1 @@ -278,6 +306,25 @@ def __init__( self._K_asymp = None self.dim = kernel.dim + def _theta( + self, + x: np.ndarray, + mu0: float, + l0: float, + sign: int, + ) -> np.ndarray: + """Dispatches a single transition-weight evaluation to hill_softtheta or softtheta. + + Only sign > 0 ever switches under ``blend == "hill"``: every call site below that passes a + strictly-positive sign is exactly one of the two *rising* weights (self.uv when the UV + branch is on is always +1; -self.ir when the IR branch is on is always +1 since self.ir is + always -1) -- see the class docstring and the params_gradient docstring below for the + region bookkeeping. sign <= 0 (falling weight, or the branch absent) always stays softtheta. + """ + if self.blend == "hill" and sign > 0: + return hill_softtheta(x, mu0, l0, sign) + return softtheta(x, mu0, l0, sign) + def _empty_cache(self) -> None: """Clear this kernel's cache and propagate to the inner kernel.""" self._K_asymp = None @@ -325,13 +372,17 @@ def make( else: x1 = make_column_vector(x[:, 0]) y1 = make_row_vector(y[:, 0]) + # self.ir (falling IR-branch weight, sign=self.ir) and -self.uv (falling base-kernel + # UV-side weight, sign=-self.uv) are always <= 0 -> always softtheta via self._theta. + # -self.ir (rising base-kernel IR-side weight) and self.uv (rising UV-branch weight) + # are always >= 0 -> hill_softtheta under blend="hill" via self._theta. _K_asymp = ( - softtheta(x1, self.mu_ir, self.l_ir, self.ir) * softtheta(y1, self.mu_ir, self.l_ir, self.ir) + self._theta(x1, self.mu_ir, self.l_ir, self.ir) * self._theta(y1, self.mu_ir, self.l_ir, self.ir) * self.ir_asymptotics(x1) * self.ir_asymptotics(y1) - + softtheta(x1, self.mu_ir, self.l_ir, -self.ir) * softtheta(y1, self.mu_ir, self.l_ir, -self.ir) + + self._theta(x1, self.mu_ir, self.l_ir, -self.ir) * self._theta(y1, self.mu_ir, self.l_ir, -self.ir) * self.kernel(x, y, cache=cache) - * softtheta(x1, self.mu_uv, self.l_uv, -self.uv) * softtheta(y1, self.mu_uv, self.l_uv, -self.uv) - + softtheta(x1, self.mu_uv, self.l_uv, self.uv) * softtheta(y1, self.mu_uv, self.l_uv, self.uv) + * self._theta(x1, self.mu_uv, self.l_uv, -self.uv) * self._theta(y1, self.mu_uv, self.l_uv, -self.uv) + + self._theta(x1, self.mu_uv, self.l_uv, self.uv) * self._theta(y1, self.mu_uv, self.l_uv, self.uv) * self.uv_asymptotics(x1) * self.uv_asymptotics(y1) ) @@ -376,6 +427,34 @@ def set_params( self.mu_uv = asymp_params[0] self.l_uv = asymp_params[1] + def _rising_l_term( + self, + x1: np.ndarray, + mu0: float, + l0: float, + sign: int, + ) -> np.ndarray: + """``t_R * (1 - R)`` where ``R = self._theta(x1, mu0, l0, sign)`` (``sign > 0``, i.e. R is + one of the two *rising* weights) and ``t_R`` is R's own internal logit variable, so that + ``dR/dl0 = -(t_R/l0) * R * (1-R)``: ``t_R = (x1-mu0)/l0`` under ``blend="logistic"`` + (R = softtheta), or ``t_R = (ln(x1)-ln(mu0))/l0`` under ``blend="hill"`` (R = hill_softtheta). + + The hill case is guarded at ``x1 = 0``: the analytic limit of ``t_R*(1-R)`` there is exactly + 0 (R -> 0 exponentially, beating t_R's logarithmic divergence), but a direct evaluation would + compute ``t_R = -inf`` (from ``ln(0)``) times ``(1-R) = 1``, giving ``-inf`` there, which a + caller then multiplies by the region's own (correctly zero) prefactor -- ``0 * (-inf) = nan`` + -- unless resolved to 0 first, here. + """ + R = self._theta(x1, mu0, l0, sign) + if self.blend == "hill" and sign > 0: + with np.errstate(divide="ignore"): + t = (np.log(x1) - np.log(mu0)) / l0 + term = t * (1.0 - R) + return np.where(x1 == 0, 0.0, term) + else: + t = (x1 - mu0) / l0 + return t * (1.0 - R) + def params_gradient(self) -> List[Callable]: """Analytic gradient of the asymptotic kernel w.r.t. all hyperparameters. @@ -387,27 +466,44 @@ def params_gradient(self) -> List[Callable]: Derivation ---------- - With the shorthands - σ_+(x) = softtheta(x, mu, l, +1) (sigmoid, 0 → 1) - σ_-(x) = softtheta(x, mu, l, −1) (complement, 1 → 0) + Every transition weight below is evaluated through ``self._theta``, which dispatches to + ``hill_softtheta`` only for the two *rising* weights (sign > 0) when ``blend == "hill"``; + every *falling* weight (sign <= 0) is always ``softtheta``, regardless of blend + (docs/rewrite_plan.md section 4, "Taken 2026-09-22 ... selectable Hill blend"). With + R_uv(x) = self._theta(x, mu_uv, l_uv, uv) (rising UV-branch weight; hill-eligible) + F_uv(x) = self._theta(x, mu_uv, l_uv, -uv) (falling base-kernel weight; always softtheta) + R_ir(x) = self._theta(x, mu_ir, l_ir, -ir) (rising base-kernel weight; hill-eligible) + F_ir(x) = self._theta(x, mu_ir, l_ir, ir) (falling IR-branch weight; always softtheta) and the three kernel regions - A = σ_-(x;ir)·σ_-(y;ir)·f_ir(x)·f_ir(y) - B = σ_+(x;ir)·σ_+(y;ir)·K_rbf·σ_-(x;uv)·σ_-(y;uv) - C = σ_+(x;uv)·σ_+(y;uv)·f_uv(x)·f_uv(y) - the sigmoid derivative identity dσ_±/dmu = ∓(1/l)·σ_+·σ_− gives: - - dK/dmu_uv = (1/l_uv)·[B·(σ_+(x;uv)+σ_+(y;uv)) − C·(σ_-(x;uv)+σ_-(y;uv))] - dK/dl_uv = (1/l_uv²)·[B·((x−mu_uv)·σ_+(x;uv)+(y−mu_uv)·σ_+(y;uv)) - −C·((x−mu_uv)·σ_-(x;uv)+(y−mu_uv)·σ_-(y;uv))] - dK/dmu_ir = (1/l_ir)·[A·(σ_+(x;ir)+σ_+(y;ir)) − B·(σ_-(x;ir)+σ_-(y;ir))] - dK/dl_ir = (1/l_ir²)·[A·((x−mu_ir)·σ_+(x;ir)+(y−mu_ir)·σ_+(y;ir)) - −B·((x−mu_ir)·σ_-(x;ir)+(y−mu_ir)·σ_-(y;ir))] - - The base-kernel gradients are windowed by the RBF transition region: - dK/d(rbf_i) = σ_+(x;ir)·σ_+(y;ir)·(dK_rbf/dtheta_i)·σ_-(x;uv)·σ_-(y;uv) - - When only one asymptotic region is active the absent softthetas reduce to - 1 (sign=0 path in softtheta), so the formulas remain valid throughout. + A = F_ir(x)·F_ir(y)·f_ir(x)·f_ir(y) + B = R_ir(x)·R_ir(y)·K_rbf·F_uv(x)·F_uv(y) + C = R_uv(x)·R_uv(y)·f_uv(x)·f_uv(y) + every weight obeys a sigmoid-derivative identity in its own logit variable t (``t = (x-mu)/l`` + for softtheta; ``t = (ln x - ln mu)/l`` for hill_softtheta, so ``dt/dmu = -1/(l*mu)`` instead + of ``-1/l``): ``dW/dmu = -(1/l_eff)*W*(1-W)`` (rising) / ``+(1/l)*W*(1-W)`` (falling), and + ``dW/dl = -(t/l)*W*(1-W)`` (rising) / ``+(t/l)*W*(1-W)`` (falling), with ``l_eff = l`` under + "logistic" or ``l*mu`` under "hill" for the rising weights only (falling weights are always + softtheta, so their own ``l_eff`` is always plain ``l``). Because F_uv and F_ir are always + softtheta, ``1 - F_uv`` and ``1 - F_ir`` used below are exactly the "shadow" logistic rising + value (softtheta's own algebraic complement) even under ``blend="hill"`` -- a different + quantity from R_uv/R_ir themselves whenever hill is active, which is exactly why the two + blends' formulas below only coincide under "logistic" (their bit-identity gate). + + dK/dmu_uv = (1/l_uv)·B·((1-F_uv(x))+(1-F_uv(y))) − (1/l_eff_uv)·C·((1-R_uv(x))+(1-R_uv(y))) + dK/dl_uv = (1/l_uv²)·B·((x−mu_uv)(1-F_uv(x))+(y−mu_uv)(1-F_uv(y))) + − (1/l_uv)·C·(t_R_uv(x)+t_R_uv(y)) [t_R_uv = self._rising_l_term(..., uv)] + dK/dmu_ir = (1/l_ir)·A·((1-F_ir(x))+(1-F_ir(y))) − (1/l_eff_ir)·B·((1-R_ir(x))+(1-R_ir(y))) + dK/dl_ir = (1/l_ir²)·A·((x−mu_ir)(1-F_ir(x))+(y−mu_ir)(1-F_ir(y))) + − (1/l_ir)·B·(t_R_ir(x)+t_R_ir(y)) [t_R_ir = self._rising_l_term(..., -ir)] + + The base-kernel gradients are windowed the same way the kernel itself is: + dK/d(rbf_i) = R_ir(x)·R_ir(y)·(dK_rbf/dtheta_i)·F_uv(x)·F_uv(y) + + Under ``blend="logistic"`` every R equals its region's own ``1 - F`` (softtheta's algebraic + complement) and every ``l_eff`` reduces to plain ``l``, reproducing the original single-blend + formulas (and thus the original kernel matrices) exactly. When only one asymptotic region is + active the absent weights reduce to 1 (sign=0 path, unaffected by blend), so the formulas + remain valid throughout. """ grads = [] @@ -417,10 +513,10 @@ def params_gradient(self) -> List[Callable]: def _wrap(x, y, _bg=bg): x1 = make_column_vector(x[:, 0]) y1 = make_row_vector(y[:, 0]) - tau_x = softtheta(x1, self.mu_ir, self.l_ir, -self.ir) # σ_+(x;ir) - tau_y = softtheta(y1, self.mu_ir, self.l_ir, -self.ir) - sig_x = softtheta(x1, self.mu_uv, self.l_uv, -self.uv) # σ_-(x;uv) - sig_y = softtheta(y1, self.mu_uv, self.l_uv, -self.uv) + tau_x = self._theta(x1, self.mu_ir, self.l_ir, -self.ir) # R_ir(x) + tau_y = self._theta(y1, self.mu_ir, self.l_ir, -self.ir) + sig_x = self._theta(x1, self.mu_uv, self.l_uv, -self.uv) # F_uv(x) + sig_y = self._theta(y1, self.mu_uv, self.l_uv, -self.uv) return tau_x * tau_y * _bg(x, y) * sig_x * sig_y grads.append(_wrap) @@ -429,38 +525,42 @@ def _wrap(x, y, _bg=bg): def _dK_dmu_uv(x, y): x1 = make_column_vector(x[:, 0]) y1 = make_row_vector(y[:, 0]) - tau_x = softtheta(x1, self.mu_ir, self.l_ir, -self.ir) - tau_y = softtheta(y1, self.mu_ir, self.l_ir, -self.ir) - sp_x = softtheta(x1, self.mu_uv, self.l_uv, self.uv) # σ_+(x;uv) - sp_y = softtheta(y1, self.mu_uv, self.l_uv, self.uv) - sm_x = softtheta(x1, self.mu_uv, self.l_uv, -self.uv) # σ_-(x;uv) - sm_y = softtheta(y1, self.mu_uv, self.l_uv, -self.uv) + tau_x = self._theta(x1, self.mu_ir, self.l_ir, -self.ir) + tau_y = self._theta(y1, self.mu_ir, self.l_ir, -self.ir) + R_x = self._theta(x1, self.mu_uv, self.l_uv, self.uv) # rising UV-branch weight + R_y = self._theta(y1, self.mu_uv, self.l_uv, self.uv) + F_x = self._theta(x1, self.mu_uv, self.l_uv, -self.uv) # falling base-kernel weight + F_y = self._theta(y1, self.mu_uv, self.l_uv, -self.uv) K_rbf = self.kernel(x, y) f_x = self.uv_asymptotics(x1) f_y = self.uv_asymptotics(y1) - return (1.0 / self.l_uv) * ( - tau_x * tau_y * K_rbf * sm_x * sm_y * (sp_x + sp_y) - - sp_x * sp_y * f_x * f_y * (sm_x + sm_y) - ) + B = tau_x * tau_y * K_rbf * F_x * F_y + C = R_x * R_y * f_x * f_y + # Rising weight's own dt/dmu normalization: -1/l (logistic) or -1/(l*mu) (hill). + l_eff = self.l_uv * self.mu_uv if self.blend == "hill" else self.l_uv + return (1.0 / self.l_uv) * B * ((1.0 - F_x) + (1.0 - F_y)) \ + - (1.0 / l_eff) * C * ((1.0 - R_x) + (1.0 - R_y)) def _dK_dl_uv(x, y): x1 = make_column_vector(x[:, 0]) y1 = make_row_vector(y[:, 0]) - tau_x = softtheta(x1, self.mu_ir, self.l_ir, -self.ir) - tau_y = softtheta(y1, self.mu_ir, self.l_ir, -self.ir) - sp_x = softtheta(x1, self.mu_uv, self.l_uv, self.uv) - sp_y = softtheta(y1, self.mu_uv, self.l_uv, self.uv) - sm_x = softtheta(x1, self.mu_uv, self.l_uv, -self.uv) - sm_y = softtheta(y1, self.mu_uv, self.l_uv, -self.uv) + tau_x = self._theta(x1, self.mu_ir, self.l_ir, -self.ir) + tau_y = self._theta(y1, self.mu_ir, self.l_ir, -self.ir) + F_x = self._theta(x1, self.mu_uv, self.l_uv, -self.uv) + F_y = self._theta(y1, self.mu_uv, self.l_uv, -self.uv) + R_x = self._theta(x1, self.mu_uv, self.l_uv, self.uv) + R_y = self._theta(y1, self.mu_uv, self.l_uv, self.uv) K_rbf = self.kernel(x, y) f_x = self.uv_asymptotics(x1) f_y = self.uv_asymptotics(y1) dx = x1 - self.mu_uv dy = y1 - self.mu_uv - return (1.0 / self.l_uv**2) * ( - tau_x * tau_y * K_rbf * sm_x * sm_y * (dx * sp_x + dy * sp_y) - - sp_x * sp_y * f_x * f_y * (dx * sm_x + dy * sm_y) - ) + B = tau_x * tau_y * K_rbf * F_x * F_y + C = R_x * R_y * f_x * f_y + tR_x = self._rising_l_term(x1, self.mu_uv, self.l_uv, self.uv) # t_R*(1-R), x=0-safe + tR_y = self._rising_l_term(y1, self.mu_uv, self.l_uv, self.uv) + return (1.0 / self.l_uv**2) * B * (dx * (1.0 - F_x) + dy * (1.0 - F_y)) \ + - (1.0 / self.l_uv) * C * (tR_x + tR_y) grads.extend([_dK_dmu_uv, _dK_dl_uv]) @@ -469,38 +569,41 @@ def _dK_dl_uv(x, y): def _dK_dmu_ir(x, y): x1 = make_column_vector(x[:, 0]) y1 = make_row_vector(y[:, 0]) - tp_x = softtheta(x1, self.mu_ir, self.l_ir, -self.ir) # σ_+(x;ir) - tp_y = softtheta(y1, self.mu_ir, self.l_ir, -self.ir) - tm_x = softtheta(x1, self.mu_ir, self.l_ir, self.ir) # σ_-(x;ir) - tm_y = softtheta(y1, self.mu_ir, self.l_ir, self.ir) - sm_x = softtheta(x1, self.mu_uv, self.l_uv, -self.uv) - sm_y = softtheta(y1, self.mu_uv, self.l_uv, -self.uv) + R_x = self._theta(x1, self.mu_ir, self.l_ir, -self.ir) # rising base-kernel weight + R_y = self._theta(y1, self.mu_ir, self.l_ir, -self.ir) + F_x = self._theta(x1, self.mu_ir, self.l_ir, self.ir) # falling IR-branch weight + F_y = self._theta(y1, self.mu_ir, self.l_ir, self.ir) + sm_x = self._theta(x1, self.mu_uv, self.l_uv, -self.uv) + sm_y = self._theta(y1, self.mu_uv, self.l_uv, -self.uv) K_rbf = self.kernel(x, y) f_x = self.ir_asymptotics(x1) f_y = self.ir_asymptotics(y1) - return (1.0 / self.l_ir) * ( - tm_x * tm_y * f_x * f_y * (tp_x + tp_y) - - tp_x * tp_y * K_rbf * sm_x * sm_y * (tm_x + tm_y) - ) + A = F_x * F_y * f_x * f_y + B = R_x * R_y * K_rbf * sm_x * sm_y + l_eff = self.l_ir * self.mu_ir if self.blend == "hill" else self.l_ir + return (1.0 / self.l_ir) * A * ((1.0 - F_x) + (1.0 - F_y)) \ + - (1.0 / l_eff) * B * ((1.0 - R_x) + (1.0 - R_y)) def _dK_dl_ir(x, y): x1 = make_column_vector(x[:, 0]) y1 = make_row_vector(y[:, 0]) - tp_x = softtheta(x1, self.mu_ir, self.l_ir, -self.ir) - tp_y = softtheta(y1, self.mu_ir, self.l_ir, -self.ir) - tm_x = softtheta(x1, self.mu_ir, self.l_ir, self.ir) - tm_y = softtheta(y1, self.mu_ir, self.l_ir, self.ir) - sm_x = softtheta(x1, self.mu_uv, self.l_uv, -self.uv) - sm_y = softtheta(y1, self.mu_uv, self.l_uv, -self.uv) + R_x = self._theta(x1, self.mu_ir, self.l_ir, -self.ir) + R_y = self._theta(y1, self.mu_ir, self.l_ir, -self.ir) + F_x = self._theta(x1, self.mu_ir, self.l_ir, self.ir) + F_y = self._theta(y1, self.mu_ir, self.l_ir, self.ir) + sm_x = self._theta(x1, self.mu_uv, self.l_uv, -self.uv) + sm_y = self._theta(y1, self.mu_uv, self.l_uv, -self.uv) K_rbf = self.kernel(x, y) f_x = self.ir_asymptotics(x1) f_y = self.ir_asymptotics(y1) dx = x1 - self.mu_ir dy = y1 - self.mu_ir - return (1.0 / self.l_ir**2) * ( - tm_x * tm_y * f_x * f_y * (dx * tp_x + dy * tp_y) - - tp_x * tp_y * K_rbf * sm_x * sm_y * (dx * tm_x + dy * tm_y) - ) + A = F_x * F_y * f_x * f_y + B = R_x * R_y * K_rbf * sm_x * sm_y + tR_x = self._rising_l_term(x1, self.mu_ir, self.l_ir, -self.ir) # t_R*(1-R), x=0-safe + tR_y = self._rising_l_term(y1, self.mu_ir, self.l_ir, -self.ir) + return (1.0 / self.l_ir**2) * A * (dx * (1.0 - F_x) + dy * (1.0 - F_y)) \ + - (1.0 / self.l_ir) * B * (tR_x + tR_y) grads.extend([_dK_dmu_ir, _dK_dl_ir]) diff --git a/fredipy/util.py b/fredipy/util.py index addce83..69af53a 100644 --- a/fredipy/util.py +++ b/fredipy/util.py @@ -61,3 +61,60 @@ def softtheta( else: return np.ones_like(x) + + +def hill_softtheta( + x: np.ndarray, + mu0: float, + l0: float, + sign: int + ) -> np.ndarray: + """Hill-type smooth step function: a logistic in ``t = (ln(x) - ln(mu0)) / l0`` rather than in + ``x`` itself (``softtheta`` above is the linear-``x`` logistic). Selected per kernel by + ``AsymptoticKernel(kernel, blend="hill")`` (docs/rewrite_plan.md section 4, "Taken 2026-09-22 + (gluon coworking ... selectable Hill blend)"): only the *rising* weights switch to this shape, + so it only ever needs to handle sign > 0 (and sign == 0) in practice, but sign < 0 is + implemented too, as the exact algebraic complement, mirroring ``softtheta``'s own structure. + + Closed form (sign > 0): ``1 / (1 + (mu0/x)^(1/l0)) == x^n / (x^n + mu0^n)`` with ``n = 1/l0`` -- + vanishes like ``(x/mu0)^n`` as ``x -> 0`` and saturates to 1 as ``x -> infinity``. + + Parameters + ---------- + x : array of position values, must be >= 0 (this function is only ever evaluated on the + omega axis, which is non-negative by construction; a negative value would make ``ln(x)`` + complex rather than raising a clear error). + mu0 : position of the 'cutoff', in the same (linear) units as x; must be > 0 (``ln(mu0)`` must + be finite). + l0 : strength of the 'cutoff' in ln(x)-space, must be > 0 (mirrors softtheta's l0, but here it + sets a width in ln(x), not in x). + sign: if positive, from 0 to 1 as x increases; if negative, the other way around (the exact + complement, matching softtheta's sign<0 convention). If sign == 0, no cutoff, returns ones. + """ + if mu0 <= 0 or not np.isfinite(mu0): + raise ValueError(f"hill_softtheta requires mu0 > 0 (finite), got {mu0}") + if l0 <= 0 or not np.isfinite(l0): + raise ValueError(f"hill_softtheta requires l0 > 0 (finite), got {l0}") + if np.any(np.asarray(x) < 0): + raise ValueError("hill_softtheta requires x >= 0 (the omega axis is non-negative)") + + if sign == 0: + return np.ones_like(x) + + # sign > 0 branch, computed once and reused (via the complement) for sign < 0 -- same pattern + # as softtheta's own sign<0-as-complement-of-sign>0 comment above. + # t = (ln x - ln mu0) / l0 diverges to -inf exactly at x = 0; np.log(0) itself raises a + # 'divide by zero' RuntimeWarning that is not a real error here (the -inf limit is the + # mathematically correct input to the sigmoid below), so it is suppressed locally. + with np.errstate(divide="ignore"): + t = (np.log(x) - np.log(mu0)) / l0 + # Overflow-safe logistic: e = exp(-|t|) is always in [0, 1] (never overflows), and the two + # equivalent forms 1/(1+e) (t>=0) / e/(1+e) (t<0) each stay well-conditioned on their own side. + # At x = 0 (t = -inf): e = exp(-inf) = 0.0 exactly (no warning), giving e/(1+e) = 0 exactly. + e = np.exp(-np.abs(t)) + pos = np.where(t >= 0, 1.0 / (1.0 + e), e / (1.0 + e)) + + if sign > 0: + return pos + else: + return 1.0 - pos diff --git a/tests/test_kernels.py b/tests/test_kernels.py index b8786b9..af60aca 100644 --- a/tests/test_kernels.py +++ b/tests/test_kernels.py @@ -483,3 +483,157 @@ def test_matern32_derivatives() -> None: def test_matern52_derivatives() -> None: """dK_dx/dK_dy/d2K_dxdy vs. finite differences for Matern-5/2 (plan section 5).""" _check_matern_derivatives(kernels.Matern52, 2.0, 1.5) + + +# --------------------------------------------------------------------------------------------- +# Selectable Hill blend (docs/rewrite_plan.md section 4, "Taken 2026-09-22 (gluon coworking: +# data_cov scan, selectable Hill blend)"). "logistic" is the kernel's original, always-linear-in-x +# blend; "hill" switches only the two *rising* transition weights to a logistic in ln(x) +# (hill_softtheta), leaving the two *falling* weights (and hence every large-w cutoff) unchanged. +# --------------------------------------------------------------------------------------------- + +def _make_asymp_kernel( + blend: str, + uv: bool = True, + ir: bool = True, + mu_uv: float = 3.0, + l_uv: float = 0.3, + mu_ir: float = 0.5, + l_ir: float = 0.1, + variance: float = 1.3, + lengthscale: float = 1.1, + ) -> "kernels.AsymptoticKernel": + """AsymptoticKernel with a plain RBF base and the requested blend/branches, for the gradient + and validation tests below. UV shape (exp(-w), bounded everywhere including w=0) and IR shape + (w**2, the repo's actual propagator-type IR form) are deliberately simple synthetic stand-ins, + not physics kernels -- this file tests the generic AsymptoticKernel machinery only.""" + base = kernels.RadialBasisFunction(variance=variance, lengthscale=lengthscale) + kernel = kernels.AsymptoticKernel(kernel=base, blend=blend) + if uv: + kernel.add_asymptotics(region="UV", asymptotics=lambda w: np.exp(-w)) + if ir: + kernel.add_asymptotics(region="IR", asymptotics=lambda w: w**2) + if uv and ir: + asymp_params = [mu_uv, l_uv, mu_ir, l_ir] + elif ir: + asymp_params = [mu_ir, l_ir] + else: + asymp_params = [mu_uv, l_uv] + kernel.set_params(asymp_params=asymp_params, kernel_params=[variance, lengthscale]) + return kernel + + +def test_asymptotic_kernel_blend_defaults_to_logistic() -> None: + rbf = kernels.RadialBasisFunction(variance=1.0, lengthscale=1.0) + kernel = kernels.AsymptoticKernel(kernel=rbf) + assert kernel.blend == "logistic" + assert kernel.blend in kernels.BLEND_SHAPES + + +@pytest.mark.parametrize("blend", kernels.BLEND_SHAPES) +def test_asymptotic_kernel_blend_selectable(blend) -> None: + rbf = kernels.RadialBasisFunction(variance=1.0, lengthscale=1.0) + kernel = kernels.AsymptoticKernel(kernel=rbf, blend=blend) + assert kernel.blend == blend + + +def test_asymptotic_kernel_invalid_blend_raises() -> None: + rbf = kernels.RadialBasisFunction(variance=1.0, lengthscale=1.0) + with pytest.raises(ValueError): + kernels.AsymptoticKernel(kernel=rbf, blend="not-a-blend") + + +@pytest.mark.parametrize("blend, uv, ir", [ + ("logistic", True, True), + ("logistic", True, False), + ("logistic", False, True), + ("hill", True, True), + ("hill", True, False), + ("hill", False, True), + ]) +def test_asymptotic_kernel_params_gradient_matches_finite_differences(blend, uv, ir) -> None: + """AsymptoticKernel.params_gradient vs. central finite differences, for both blends and every + branch combination, at nodes spanning w = 0 (the deep-IR limit the hill x=0 guard exists for) + to far UV (docs/rewrite_plan.md section 4 "selectable Hill blend"; gates the per-factor + hill/logistic derivative split documented in params_gradient's own docstring).""" + kernel = _make_asymp_kernel(blend, uv=uv, ir=ir) + + # Nodes deliberately include x = 0 and a far-UV node, plus mid-range points straddling both + # transition windows (mu_uv=3.0/l_uv=0.3, mu_ir=0.5/l_ir=0.1 from _make_asymp_kernel). + X = make_column_vector(np.array([0.0, 0.01, 0.2, 0.5, 1.0, 3.0, 10.0])) + Y = make_column_vector(np.array([0.0, 0.05, 0.3, 0.8, 2.0, 5.0, 20.0])) + + if uv and ir: + params = [1.3, 1.1, 3.0, 0.3, 0.5, 0.1] + names = ["variance", "lengthscale", "mu_uv", "l_uv", "mu_ir", "l_ir"] + elif ir: + params = [1.3, 1.1, 0.5, 0.1] + names = ["variance", "lengthscale", "mu_ir", "l_ir"] + else: + params = [1.3, 1.1, 3.0, 0.3] + names = ["variance", "lengthscale", "mu_uv", "l_uv"] + + def _set(p): + kernel.set_params(asymp_params=p[2:], kernel_params=p[:2]) + + _set(params) + grads = kernel.params_gradient() + analytical = [g(X, Y) for g in grads] + for a in analytical: + assert np.all(np.isfinite(a)) + + eps = 1e-6 + for idx, name in enumerate(names): + p_plus = list(params) + p_plus[idx] += eps + p_minus = list(params) + p_minus[idx] -= eps + + _set(p_plus) + kernel._empty_cache() + K_plus = kernel.make(X, Y, cache=False) + + _set(p_minus) + kernel._empty_cache() + K_minus = kernel.make(X, Y, cache=False) + + numerical = (K_plus - K_minus) / (2 * eps) + assert_allclose( + analytical[idx], numerical, rtol=1e-4, atol=1e-8, + err_msg=f"blend={blend} uv={uv} ir={ir} param={name}", + ) + + _set(params) + + +def test_asymptotic_kernel_hill_mode_base_kernel_term_vanishes_at_zero() -> None: + """Regression guard for the deep-IR plateau bug (docs/rewrite_plan.md section 4, "Deep-IR + plateau root cause": the base kernel has a floor on the IR side too, sigma_+^IR(0), the + ``_ir_leak_note`` rule). Isolates the base-kernel (B) term by turning the UV branch off + (self.uv=0 makes its softtheta(sign=0)=1 windowing a no-op) and using an IR shape (constant 1) + that does not itself vanish at w=0, so the IR-asymptotics (A) term does not mask B's own + behaviour there. Under blend="hill" the base-kernel weight (rising, hill-eligible) is exactly + zero at w=0, so K(0,0) is exactly the A term alone; under blend="logistic" it floors at the + nonzero constant predicted by fredipy.util.softtheta directly.""" + from fredipy.util import softtheta + + mu_ir, l_ir, variance, lengthscale = 0.3, 0.15, 1.3, 1.1 # mu_ir/l_ir=2: a sizeable, non-tiny floor + + def build(blend): + base = kernels.RadialBasisFunction(variance=variance, lengthscale=lengthscale) + kernel = kernels.AsymptoticKernel(kernel=base, blend=blend) + kernel.add_asymptotics(region="IR", asymptotics=lambda w: np.ones_like(w)) + kernel.set_params(asymp_params=[mu_ir, l_ir], kernel_params=[variance, lengthscale]) + return kernel + + x0 = make_column_vector(np.array([0.0])) + K_hill = build("hill").make(x0, x0, cache=False).item() + K_logistic = build("logistic").make(x0, x0, cache=False).item() + + F_ir_0 = softtheta(np.array([0.0]), mu_ir, l_ir, -1).item() # A term's own weight at w=0 + A_0 = F_ir_0 ** 2 # ir_asymptotics(0)**2 == 1 + R_ir_logistic_0 = softtheta(np.array([0.0]), mu_ir, l_ir, 1).item() + B_logistic_0 = R_ir_logistic_0 ** 2 * variance ** 2 # base kernel value at (0,0) is variance**2 + + assert K_hill == pytest.approx(A_0, abs=1e-12) # B term exactly zero under hill + assert K_logistic == pytest.approx(A_0 + B_logistic_0, rel=1e-10) # B term the nonzero floor diff --git a/tests/test_util.py b/tests/test_util.py index c2d377c..3d23942 100644 --- a/tests/test_util.py +++ b/tests/test_util.py @@ -1,6 +1,7 @@ import numpy as np +import pytest -from fredipy.util import softtheta, make_column_vector, make_row_vector, allclose +from fredipy.util import softtheta, hill_softtheta, make_column_vector, make_row_vector, allclose def test_softtheta(): @@ -40,6 +41,78 @@ def test_softtheta_sign_negative_does_not_overflow_to_nan(): assert np.allclose(result, 1.0) +def test_hill_softtheta_closed_form(): + # Closed form for sign > 0: 1/(1+(mu0/x)^(1/l0)) == x^n/(x^n+mu0^n) with n = 1/l0. + x = np.array([0.5, 1.0, 2.0, 10.0]) + mu0, l0 = 2.0, 0.5 # n = 1/l0 = 2 + n = 1.0 / l0 + + result = hill_softtheta(x, mu0, l0, 1) + expected = x**n / (x**n + mu0**n) + assert np.allclose(result, expected) + + # sign < 0 is the exact algebraic complement. + result_neg = hill_softtheta(x, mu0, l0, -1) + assert np.allclose(result_neg, 1.0 - expected) + + # sign == 0: no cutoff, matches softtheta's own sign==0 convention. + result_zero = hill_softtheta(x, mu0, l0, 0) + assert np.allclose(result_zero, np.ones_like(x)) + + +def test_hill_softtheta_zero_is_exact_without_warnings(): + # x = 0 is the deep-IR/deep-UV limit this function exists for: sign > 0 must vanish exactly, + # sign < 0 must saturate to exactly 1, and neither may emit a RuntimeWarning (log(0), 0*inf). + x = np.array([0.0, 1.0, 5.0]) + mu0, l0 = 1.0, 0.2 + + with np.errstate(all="raise"): + pos = hill_softtheta(x, mu0, l0, 1) + neg = hill_softtheta(x, mu0, l0, -1) + + assert pos[0] == 0.0 + assert neg[0] == 1.0 + assert np.all(np.isfinite(pos)) + assert np.all(np.isfinite(neg)) + + +def test_hill_softtheta_extremes_finite(): + # A very narrow transition (small l0) and points far on either side of mu0 must stay finite -- + # the overflow-safe e = exp(-|t|) construction is exactly what this guards. + x = np.array([1e-12, 1e-6, 1e6, 1e12]) + mu0, l0 = 1.0, 1e-3 + + for sign in (1, -1): + result = hill_softtheta(x, mu0, l0, sign) + assert np.all(np.isfinite(result)) + assert np.all(result >= 0.0) and np.all(result <= 1.0) + + +def test_hill_softtheta_negative_x_raises(): + with pytest.raises(ValueError): + hill_softtheta(np.array([-1.0, 1.0]), 1.0, 0.5, 1) + + +def test_hill_softtheta_nonpositive_mu0_or_l0_raises(): + x = np.array([1.0]) + with pytest.raises(ValueError): + hill_softtheta(x, 0.0, 0.5, 1) + with pytest.raises(ValueError): + hill_softtheta(x, -1.0, 0.5, 1) + with pytest.raises(ValueError): + hill_softtheta(x, 1.0, 0.0, 1) + with pytest.raises(ValueError): + hill_softtheta(x, 1.0, -0.5, 1) + + +def test_softtheta_unchanged_by_hill_addition(): + # Regression guard: adding hill_softtheta must not perturb softtheta's own behaviour. + x = np.array([-1e2, 0, 1e2]) + assert np.allclose(softtheta(x, 0, 1, 1), np.array([0, 0.5, 1])) + assert np.allclose(softtheta(x, 0, 1, -1), np.array([1, 0.5, 0])) + assert np.allclose(softtheta(x, 0, 1, 0), np.ones_like(x)) + + def test_allclose(): a = np.array([1, 2, 3]) b = np.array([1, 2, 3]) From 11b85bb497ee1ce1f828594473e40411577e4a2d Mon Sep 17 00:00:00 2001 From: Jonas Wessely <104987725+Jpwessely@users.noreply.github.com> Date: Tue, 22 Sep 2026 14:24:01 +0200 Subject: [PATCH 8/8] AsymptoticKernel hill mode: slope-matched falling width l*mu Under blend='hill' the falling (logistic) weights use the linear width l*mu, the same slope at mu as the Hill rise, so l is the width in ln(w) on both sides; gradients include the width's mu-dependence. 'logistic' unchanged bit for bit. Co-Authored-By: Claude Opus 5 --- fredipy/kernels.py | 47 +++++++++++++++++++++++++++++++++++++++---- tests/test_kernels.py | 6 ++++-- 2 files changed, 47 insertions(+), 6 deletions(-) diff --git a/fredipy/kernels.py b/fredipy/kernels.py index f39d501..33826c9 100644 --- a/fredipy/kernels.py +++ b/fredipy/kernels.py @@ -276,8 +276,9 @@ class AsymptoticKernel(Kernel): transition weight is ``softtheta``, a logistic in x itself. "hill" switches only the two *rising* weights -- the base-kernel weight above mu_ir, the UV-branch weight above mu_uv -- to ``hill_softtheta``, a logistic in ln(x); the two *falling* weights (the IR-branch weight, - the base-kernel weight below mu_uv) are left as softtheta in both modes, so every large-x - cutoff stays exponential regardless of blend (docs/rewrite_plan.md section 4, "Taken + the base-kernel weight below mu_uv) stay softtheta in both modes, so every large-x cutoff + stays exponential regardless of blend -- under "hill" with the slope-matched linear width + l*mu, so l is the width in ln(x) on both sides of a transition (docs/rewrite_plan.md section 4, "Taken 2026-09-22 (gluon coworking ... selectable Hill blend)"). "logistic" is bit-identical to the kernel's pre-blend-selector behaviour. """ @@ -319,12 +320,30 @@ def _theta( strictly-positive sign is exactly one of the two *rising* weights (self.uv when the UV branch is on is always +1; -self.ir when the IR branch is on is always +1 since self.ir is always -1) -- see the class docstring and the params_gradient docstring below for the - region bookkeeping. sign <= 0 (falling weight, or the branch absent) always stays softtheta. + region bookkeeping. sign <= 0 (falling weight, or the branch absent) always stays softtheta; + under ``blend == "hill"`` a falling weight (sign < 0) uses the linear width ``l0*mu0``, which + gives it the same slope at mu0 as the Hill rise, so l0 is the width in ln(x) on both sides + (docs/rewrite_plan.md section 4, "Correction decided ... slope-matched falling width"). """ if self.blend == "hill" and sign > 0: return hill_softtheta(x, mu0, l0, sign) + if self.blend == "hill" and sign < 0: + return softtheta(x, mu0, l0 * mu0, sign) return softtheta(x, mu0, l0, sign) + def _falling_mu_factor(self, x1: np.ndarray, mu0: float, l0: float) -> np.ndarray | float: + """``c`` with ``dF/dmu0 = c * F * (1-F)`` for a falling weight F: ``1/l0`` under "logistic" + (t = (x-mu0)/l0); ``x/(l0*mu0**2)`` under "hill", where F uses width l0*mu0 so + t = (x-mu0)/(l0*mu0) and mu0 also enters through the width.""" + if self.blend == "hill": + return x1 / (l0 * mu0**2) + return 1.0 / l0 + + def _falling_l_width(self, mu0: float, l0: float) -> float: + """``w`` with ``dF/dl0 = ((x-mu0)/(l0*w)) * F * (1-F)`` for a falling weight F: ``l0`` under + "logistic", ``l0*mu0`` (the falling weight's own linear width) under "hill".""" + return l0 * mu0 if self.blend == "hill" else l0 + def _empty_cache(self) -> None: """Clear this kernel's cache and propagate to the inner kernel.""" self._K_asymp = None @@ -468,7 +487,11 @@ def params_gradient(self) -> List[Callable]: ---------- Every transition weight below is evaluated through ``self._theta``, which dispatches to ``hill_softtheta`` only for the two *rising* weights (sign > 0) when ``blend == "hill"``; - every *falling* weight (sign <= 0) is always ``softtheta``, regardless of blend + every *falling* weight (sign <= 0) is always ``softtheta``, regardless of blend, but under + "hill" with the linear width ``l*mu`` (slope-matched to the Hill rise), so its logit is + ``t = (x-mu)/(l*mu)``: ``dF/dmu = (x/(l*mu^2))*F*(1-F)`` and + ``dF/dl = ((x-mu)/(l^2*mu))*F*(1-F)`` (``_falling_mu_factor`` / ``_falling_l_width``); the + formulas below are written for "logistic" and carry these two replacements under "hill" (docs/rewrite_plan.md section 4, "Taken 2026-09-22 ... selectable Hill blend"). With R_uv(x) = self._theta(x, mu_uv, l_uv, uv) (rising UV-branch weight; hill-eligible) F_uv(x) = self._theta(x, mu_uv, l_uv, -uv) (falling base-kernel weight; always softtheta) @@ -538,6 +561,10 @@ def _dK_dmu_uv(x, y): C = R_x * R_y * f_x * f_y # Rising weight's own dt/dmu normalization: -1/l (logistic) or -1/(l*mu) (hill). l_eff = self.l_uv * self.mu_uv if self.blend == "hill" else self.l_uv + if self.blend == "hill": # Falling weight's width l*mu also depends on mu_uv. + return B * (self._falling_mu_factor(x1, self.mu_uv, self.l_uv) * (1.0 - F_x) + + self._falling_mu_factor(y1, self.mu_uv, self.l_uv) * (1.0 - F_y)) \ + - (1.0 / l_eff) * C * ((1.0 - R_x) + (1.0 - R_y)) return (1.0 / self.l_uv) * B * ((1.0 - F_x) + (1.0 - F_y)) \ - (1.0 / l_eff) * C * ((1.0 - R_x) + (1.0 - R_y)) @@ -559,6 +586,10 @@ def _dK_dl_uv(x, y): C = R_x * R_y * f_x * f_y tR_x = self._rising_l_term(x1, self.mu_uv, self.l_uv, self.uv) # t_R*(1-R), x=0-safe tR_y = self._rising_l_term(y1, self.mu_uv, self.l_uv, self.uv) + if self.blend == "hill": # Falling weight's linear width is l_uv*mu_uv. + w_fall = self._falling_l_width(self.mu_uv, self.l_uv) + return (1.0 / (self.l_uv * w_fall)) * B * (dx * (1.0 - F_x) + dy * (1.0 - F_y)) \ + - (1.0 / self.l_uv) * C * (tR_x + tR_y) return (1.0 / self.l_uv**2) * B * (dx * (1.0 - F_x) + dy * (1.0 - F_y)) \ - (1.0 / self.l_uv) * C * (tR_x + tR_y) @@ -581,6 +612,10 @@ def _dK_dmu_ir(x, y): A = F_x * F_y * f_x * f_y B = R_x * R_y * K_rbf * sm_x * sm_y l_eff = self.l_ir * self.mu_ir if self.blend == "hill" else self.l_ir + if self.blend == "hill": # Falling weight's width l*mu also depends on mu_ir. + return A * (self._falling_mu_factor(x1, self.mu_ir, self.l_ir) * (1.0 - F_x) + + self._falling_mu_factor(y1, self.mu_ir, self.l_ir) * (1.0 - F_y)) \ + - (1.0 / l_eff) * B * ((1.0 - R_x) + (1.0 - R_y)) return (1.0 / self.l_ir) * A * ((1.0 - F_x) + (1.0 - F_y)) \ - (1.0 / l_eff) * B * ((1.0 - R_x) + (1.0 - R_y)) @@ -602,6 +637,10 @@ def _dK_dl_ir(x, y): B = R_x * R_y * K_rbf * sm_x * sm_y tR_x = self._rising_l_term(x1, self.mu_ir, self.l_ir, -self.ir) # t_R*(1-R), x=0-safe tR_y = self._rising_l_term(y1, self.mu_ir, self.l_ir, -self.ir) + if self.blend == "hill": # Falling weight's linear width is l_ir*mu_ir. + w_fall = self._falling_l_width(self.mu_ir, self.l_ir) + return (1.0 / (self.l_ir * w_fall)) * A * (dx * (1.0 - F_x) + dy * (1.0 - F_y)) \ + - (1.0 / self.l_ir) * B * (tR_x + tR_y) return (1.0 / self.l_ir**2) * A * (dx * (1.0 - F_x) + dy * (1.0 - F_y)) \ - (1.0 / self.l_ir) * B * (tR_x + tR_y) diff --git a/tests/test_kernels.py b/tests/test_kernels.py index af60aca..0cd93eb 100644 --- a/tests/test_kernels.py +++ b/tests/test_kernels.py @@ -630,10 +630,12 @@ def build(blend): K_hill = build("hill").make(x0, x0, cache=False).item() K_logistic = build("logistic").make(x0, x0, cache=False).item() - F_ir_0 = softtheta(np.array([0.0]), mu_ir, l_ir, -1).item() # A term's own weight at w=0 + F_ir_0 = softtheta(np.array([0.0]), mu_ir, l_ir, -1).item() # A term's own weight at w=0 (logistic) A_0 = F_ir_0 ** 2 # ir_asymptotics(0)**2 == 1 + # Under hill the falling IR-branch weight is softtheta with the slope-matched width l*mu. + A_0_hill = softtheta(np.array([0.0]), mu_ir, l_ir * mu_ir, -1).item() ** 2 R_ir_logistic_0 = softtheta(np.array([0.0]), mu_ir, l_ir, 1).item() B_logistic_0 = R_ir_logistic_0 ** 2 * variance ** 2 # base kernel value at (0,0) is variance**2 - assert K_hill == pytest.approx(A_0, abs=1e-12) # B term exactly zero under hill + assert K_hill == pytest.approx(A_0_hill, abs=1e-12) # B term exactly zero under hill assert K_logistic == pytest.approx(A_0 + B_logistic_0, rel=1e-10) # B term the nonzero floor