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 8251736..b630382 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)) @@ -190,6 +262,445 @@ 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``. + + The Jacobian of the substitution t = log ω is ω and is applied + internally, so plain kernels ``K(p, ω)`` can be passed directly. + """ + + 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 + self.jac = make_row_vector(self.w) # log-space Jacobian: ω_i, shape (1, n) + + def doubleIntegrationSymmetric( + self, + constraint: LinearEquality, + kernel: Callable, + derivative: bool = False # unused: this rule is linear in `kernel` + ) -> np.ndarray: + return self.doubleIntegration(constraint, kernel, constraint, derivative=derivative) + + def doubleIntegration( + self, + constraint1: LinearEquality, + kernel: Callable, + 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)) + @ kernel(self.w, self.w) + @ (self.jac * 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 * ( + self.jac * 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). + + The log-space Jacobian factor ω is absorbed into the quadrature weights, + so plain kernels ``K(p, ω)`` can be passed directly. + """ + + def __init__( + self, + w_min: float = 0.01, + w_max: float = 10., + 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) + + # 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 and the log-space Jacobian ω + self.weights = make_row_vector(wi * half_range) * self.w.T # shape (1, int_n) + + def doubleIntegrationSymmetric( + self, + constraint: LinearEquality, + kernel: Callable, + derivative: bool = False + ) -> np.ndarray: + return self.doubleIntegration(constraint, kernel, constraint, derivative=derivative) + + def doubleIntegration( + self, + constraint1: LinearEquality, + kernel: Callable, + constraint2: LinearEquality, + derivative: bool = False + ) -> np.ndarray: + 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, + 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) + ) + + @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) plus an analytic UV tail. + + **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: + + 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))`` 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 + applying (*) to whichever argument lies in the tail gives + + A_i[m] = ∫_bulk C_i(p_m, ω) θ_uv(ω) f_uv(ω) dω + T_i[m] = ∫_{w_uv}^∞ C_i(p_m, ω) f_uv(ω) dω + + Σ_12 = BB + A_1 T_2^T + T_1 A_2^T + T_1 T_2^T (**) + + 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 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``. 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__( + self, + w_min: float, + w_uv: float, + int_n: int, + uv_func: Callable, + tail_moment: Callable[[np.ndarray], np.ndarray] + ) -> None: + super().__init__(w_min, w_uv, int_n) + 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) + + def uv_tail_moment( + self, + constraint: LinearEquality + ) -> np.ndarray: + 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 + + def uv_anchor(self) -> tuple: + return (self.w_uv, self.f_uv_anchor) + + def singleIntegration( + self, + constraint: LinearEquality, + kernel: Callable, + w_pred: np.ndarray + ) -> np.ndarray: + numerical = super().singleIntegration(constraint, kernel, w_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 + + +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²) and the ω factor are both absorbed into + the quadrature weights, so plain kernels ``K(p, ω)`` can be passed + directly. + + 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)) + # 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, + constraint: LinearEquality, + kernel: Callable, + derivative: bool = False # unused: this rule is linear in `kernel` + ) -> np.ndarray: + return self.doubleIntegration(constraint, kernel, constraint, derivative=derivative) + + def doubleIntegration( + self, + constraint1: LinearEquality, + kernel: Callable, + 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)) + @ 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""" @@ -215,15 +726,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 * ( @@ -270,7 +783,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/kernels.py b/fredipy/kernels.py index 3a2c87c..33826c9 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,33 @@ 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) 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. """ 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 +307,49 @@ 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; + 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 + self._x, self._y = np.array([None]), np.array([None]) + self.kernel._empty_cache() + def add_asymptotics( self, region: str, @@ -319,13 +391,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) ) @@ -370,6 +446,208 @@ 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. + + 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 + ---------- + 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, 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) + 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 = 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 = [] + + # ---- 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 = 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) + + # ---- 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 = 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) + 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 + 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)) + + def _dK_dl_uv(x, y): + x1 = make_column_vector(x[:, 0]) + y1 = make_row_vector(y[:, 0]) + 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 + 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) + 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) + + 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]) + 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) + 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)) + + def _dK_dl_ir(x, y): + x1 = make_column_vector(x[:, 0]) + y1 = make_row_vector(y[:, 0]) + 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 + 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) + 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) + + grads.extend([_dK_dmu_ir, _dK_dl_ir]) + + return grads + class Matern12(Kernel): """Matern-1/2 kernel @@ -477,6 +755,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 @@ -539,6 +899,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/fredipy/models.py b/fredipy/models.py index 12c04b6..7d1f9aa 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.""" @@ -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/fredipy/util.py b/fredipy/util.py index c7c4d3a..69af53a 100644 --- a/fredipy/util.py +++ b/fredipy/util.py @@ -53,7 +53,68 @@ 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) + + +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_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) diff --git a/tests/test_kernels.py b/tests/test_kernels.py index a886dd8..0cd93eb 100644 --- a/tests/test_kernels.py +++ b/tests/test_kernels.py @@ -332,3 +332,310 @@ 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) + + +# --------------------------------------------------------------------------------------------- +# 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 (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_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 diff --git a/tests/test_util.py b/tests/test_util.py index b7af1c9..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(): @@ -27,6 +28,91 @@ 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_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])