NEGM stacked envelope: complete right-germ tie ownership

The first right derivative alone does not identify the branch that wins
immediately right of an exact value tie: two limited Hermite pieces can
agree to first order and differ in curvature, and a piece whose right
bracket ends in -inf got derivative zero — tying a genuinely finite
branch — so the lowest-index fallback could publish the marginal of a
candidate that loses for every small step right.

The value read's local pieces are cubics or constant clamps, so the
right germ is finite-dimensional and the comparison can be exhaustive:

- interp_right_germ_on_{padded,prepared}_grid replaces the
  derivative-only read, returning (right_finite, first, second, third)
  — the same bracket location, limiter, and clamp rays; secant fallback
  carries zero curvature; a non-finite bracket endpoint clears the
  right-finiteness flag.
- right_germ_winner (shared by outer_envelope_at_query and the
  production _collapse_stacked_candidates) ranks the tie set by staged
  exact comparisons: right-finiteness, then the three derivatives, then
  the lowest index — now reached only by locally identical pieces.

The scalar NumPy host reference is extended to the full germ; the
equal-slope/different-curvature and dying-right-bracket ties are
red-first regressions in both readers.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_0123hk3uzBgBeZown4eBBQ2C

Author: hmgaudecker

src/_lcm/egm/continuation.py

@@ -26,10 +26,11 @@
 from _lcm.egm.carry import EGMCarry
 from _lcm.egm.interp import (
     interp_on_prepared_grid,
-    interp_right_derivative_on_prepared_grid,
+    interp_right_germ_on_prepared_grid,
     locate_on_grid,
     prepare_padded_grid,
 )
+from _lcm.egm.outer_envelope import right_germ_winner
 from _lcm.egm.regime_introspection import (
     _get_child_discrete_actions,
     _get_child_resources_arg_names,
@@ -56,9 +57,11 @@
 )
 from lcm.regime import Regime as UserRegime
 from lcm.typing import (
+    BoolND,
     Float1D,
     FloatND,
     IntND,
+    ScalarBool,
     ScalarFloat,
     ScalarInt,
 )
@@ -920,18 +923,21 @@ def interp_row(
         # Below a candidate's own first finite coh node its support has not
         # started: mask the read to `-inf` so the edge clamp cannot hand an
         # infeasible lifted candidate a boundary value that wins the max. The
-        # `-inf` also pins the marginal to zero below. The right derivative of
-        # each candidate's value read feeds the tie rule at the candidate max.
-        def right_slope_row(
+        # `-inf` also pins the marginal to zero below. The right germ of each
+        # candidate's value read feeds the tie rule at the candidate max; it
+        # needs no support mask of its own — a below-support candidate enters
+        # the tie set only when every candidate is below support, where all
+        # published marginals are exactly zero regardless of the winner.
+        def right_germ_row(
             search_grid: Float1D,
             valid_length: ScalarInt,
             xp: Float1D,
             fp: Float1D,
             fp_slopes: Float1D,
             x_query: ScalarFloat,
-        ) -> ScalarFloat:
-            """Right-derive one value row at its query; positional per `jax.vmap`."""
-            return interp_right_derivative_on_prepared_grid(
+        ) -> tuple[ScalarBool, ScalarFloat, ScalarFloat, ScalarFloat]:
+            """Right germ of one value row at its query; positional per `jax.vmap`."""
+            return interp_right_germ_on_prepared_grid(
                 x_query=x_query,
                 search_grid=search_grid,
                 valid_length=valid_length,
@@ -940,7 +946,7 @@ def right_slope_row(
                 fp_slopes=fp_slopes,
             )
 
-        right_slope_at_child = jax.vmap(right_slope_row)(
+        right_germ_at_child = jax.vmap(right_germ_row)(
             search_rows, valid_rows, grid_rows, value_rows, marginal_rows, queries_flat
         )
         row_lower = jnp.min(
@@ -952,7 +958,6 @@ def right_slope_row(
         # ordinary infeasible `(-inf, 0)` pair.
         below_row_support = (queries_flat < row_lower) & jnp.isfinite(row_lower)
         value_at_child = jnp.where(below_row_support, -jnp.inf, value_at_child)
-        right_slope_at_child = jnp.where(below_row_support, 0.0, right_slope_at_child)
     # `-inf` entries interpolate pointwise to `-inf` (never NaN) and carry
     # exactly-zero marginal utility, so an infeasible-everywhere row reads as
     # the `-inf` / zero pair while a row with isolated `-inf` nodes (e.g. a
@@ -968,7 +973,9 @@ def right_slope_row(
         value_at_child, marginal_at_child = _collapse_stacked_candidates(
             value_at_child=value_at_child,
             marginal_at_child=marginal_at_child,
-            right_slope_at_child=right_slope_at_child.reshape(block_shape),
+            right_germ_at_child=tuple(
+                component.reshape(block_shape) for component in right_germ_at_child
+            ),
         )
 
     for passive_value, passive_grid in zip(
@@ -1016,47 +1023,41 @@ def _collapse_stacked_candidates(
     *,
     value_at_child: FloatND,
     marginal_at_child: FloatND,
-    right_slope_at_child: FloatND,
+    right_germ_at_child: tuple[BoolND, FloatND, FloatND, FloatND],
 ) -> tuple[FloatND, FloatND]:
     """Collapse the candidate axis by the exact hard max at the query.
 
     Publishes the winner's marginal (Danskin). This runs *before* the passive
     blend so the blend interpolates the nodewise outer maximum. An exact value
-    tie resolves right-continuously in the value read itself: among tied
-    candidates the largest right derivative of the value interpolant wins — the
-    Fritsch-Carlson-limited slope of the bracket right of the query, exactly
-    zero on the clamp ray at or above a candidate's last node — so the branch
-    whose read actually wins immediately to the right owns the (economic)
-    marginal the parent's Euler inversion consumes; `argmax` resolves exact
-    rank ties to the lowest index. (The rank uses the unscaled interpolant
-    slope: the composed gradient is shared by all candidates of a cell and
-    positive, so scaling could never reorder them.) A cell whose candidates are
-    all `-inf` (no live support) keeps the `(-inf, 0)` infeasible contract:
-    every masked marginal is exactly zero.
+    tie resolves right-continuously in the value read itself: the tied
+    candidates are compared by the complete right germ of their own value
+    interpolants (`right_germ_winner` — right-finiteness, then the first three
+    one-sided derivatives, then the lowest index among locally identical
+    pieces), so the branch whose read actually wins immediately to the right
+    owns the (economic) marginal the parent's Euler inversion consumes. (The
+    germ uses the unscaled interpolant derivatives: the composed gradient is
+    shared by all candidates of a cell and positive, so scaling could never
+    reorder them.) A cell whose candidates are all `-inf` (no live support)
+    keeps the `(-inf, 0)` infeasible contract: every masked marginal is
+    exactly zero.
 
     Args:
         value_at_child: Candidate value reads; the candidate axis is last.
         marginal_at_child: Gradient-scaled candidate marginal reads, same shape.
-        right_slope_at_child: Right derivatives of the candidate value reads,
-            same shape.
+        right_germ_at_child: Tuple of the right-finiteness flag and the first
+            three right derivatives of the candidate value reads, same shape.
 
     Returns:
         Tuple of the winner's value and marginal, with the candidate axis
         collapsed.
 
     """
-    best_value = jnp.max(value_at_child, axis=-1, keepdims=True)
-    rank = jnp.where(
-        value_at_child >= best_value,
-        right_slope_at_child,
-        -jnp.inf,
-    )
-    winner = jnp.argmax(rank, axis=-1, keepdims=True)
+    winner = right_germ_winner(value=value_at_child, right_germ=right_germ_at_child)
     winner_marginal = jnp.take_along_axis(marginal_at_child, winner, axis=-1)[..., 0]
     # The published value is the maximum itself: identical to the winner's read
     # at any tie, and NaN-propagating when a poisoned candidate row (whose NaN
     # empties the tie set) must surface fail-loud.
-    return best_value[..., 0], winner_marginal
+    return jnp.max(value_at_child, axis=-1), winner_marginal
 
 
 def _build_child_reads(

src/_lcm/egm/interp.py

@@ -12,7 +12,7 @@
 
 import jax.numpy as jnp
 
-from lcm.typing import Float1D, FloatND, ScalarFloat, ScalarInt
+from lcm.typing import BoolND, Float1D, FloatND, ScalarFloat, ScalarInt
 
 
 def interp_on_padded_grid(
@@ -170,31 +170,32 @@ def interp_on_prepared_grid(
     return jnp.where(valid_length == 0, jnp.nan, result)
 
 
-def interp_right_derivative_on_padded_grid(
+def interp_right_germ_on_padded_grid(
     *,
     x_query: FloatND,
     xp: Float1D,
     fp: Float1D,
     fp_slopes: Float1D,
-) -> FloatND:
-    """Compute the right derivative of the Hermite value read on a padded row.
+) -> tuple[BoolND, FloatND, FloatND, FloatND]:
+    """Compute the right germ of the Hermite value read on a padded row.
 
     Same row contract as `interp_on_padded_grid`; see
-    `interp_right_derivative_on_prepared_grid` for the derivative semantics.
+    `interp_right_germ_on_prepared_grid` for the germ semantics.
 
     Args:
-        x_query: Points at which to evaluate the right derivative; any shape.
+        x_query: Points at which to evaluate the right germ; any shape.
         xp: Weakly ascending grid row with NaNs only in the tail.
         fp: Function values on `xp`, NaN-padded in lockstep with `xp`.
         fp_slopes: Derivatives of `fp` with respect to `xp` at the `xp` nodes,
             NaN-padded in lockstep.
 
     Returns:
-        Right derivatives with the shape of `x_query`.
+        Tuple of the right-finiteness flag and the first, second, and third
+        right derivatives, each with the shape of `x_query`.
 
     """
     search_grid, valid_length = prepare_padded_grid(xp)
-    return interp_right_derivative_on_prepared_grid(
+    return interp_right_germ_on_prepared_grid(
         x_query=x_query,
         search_grid=search_grid,
         valid_length=valid_length,
@@ -204,37 +205,40 @@ def interp_right_derivative_on_padded_grid(
     )
 
 
-def interp_right_derivative_on_prepared_grid(
+def interp_right_germ_on_prepared_grid(
     *,
     x_query: FloatND,
     search_grid: Float1D,
     valid_length: ScalarInt,
     xp: Float1D,
     fp: Float1D,
     fp_slopes: Float1D,
-) -> FloatND:
-    """Compute the right derivative of the Hermite value read at each query.
-
-    This is the one-sided derivative of the *value interpolant* that
-    `interp_on_prepared_grid` evaluates with `fp_slopes` — the local slope
-    governing the read immediately to the right of the query — not a read of
-    the slope row itself. The two differ exactly where a tie-owner decision
-    needs the truth: the Fritsch-Carlson limiter may cap a node's raw slope,
-    and the edge clamps flatten the read outside the valid range. Semantics:
-
-    - Strictly inside a bracket: the derivative of that bracket's limited
-      cubic Hermite (the secant where the correction is inapplicable — the
-      linear fallback).
-    - Exactly on a node: the derivative at the left edge of the node's *right*
-      bracket (the bracket search is right-continuous), i.e. the node's
-      limited slope with respect to that bracket's secant.
-    - Strictly below the first node, and at or above the last valid node: zero
-      — the read clamps to a constant there.
-    - A bracket with a non-finite value difference (e.g. a `-inf` endpoint):
-      zero, matching the linear rule's flat, clamp-like behavior there.
+) -> tuple[BoolND, FloatND, FloatND, FloatND]:
+    """Compute the right germ of the Hermite value read at each query.
+
+    The germ is the complete local description of the *value interpolant* that
+    `interp_on_prepared_grid` evaluates with `fp_slopes` immediately to the
+    right of the query — not a read of the slope row itself. Each local piece
+    is a cubic (or a constant clamp), so the germ is finite-dimensional: a
+    right-finiteness flag plus the first, second, and third one-sided
+    derivatives determine the read on a right neighborhood exactly. The germ
+    differs from the raw slope row exactly where a tie-owner decision needs the
+    truth: the Fritsch-Carlson limiter may cap a node's raw slope, the edge
+    clamps flatten the read outside the valid range, and a `-inf` bracket
+    endpoint kills the read immediately right of a finite node. Semantics:
+
+    - Strictly inside a bracket: the derivatives of that bracket's limited
+      cubic Hermite (the secant and zero curvature where the correction is
+      inapplicable — the linear fallback).
+    - Exactly on a node: the derivatives at the left edge of the node's *right*
+      bracket (the bracket search is right-continuous).
+    - Strictly below the first node, and at or above the last valid node: the
+      read clamps to a constant — right-finite with all derivatives zero.
+    - A bracket with a non-finite endpoint value: not right-finite (the read
+      is `-inf` on the bracket's interior), derivatives zero.
 
     Args:
-        x_query: Points at which to evaluate the right derivative; any shape.
+        x_query: Points at which to evaluate the right germ; any shape.
         search_grid: The row's `+inf`-padded search key from
             `prepare_padded_grid`.
         valid_length: The row's non-NaN prefix length from
@@ -244,12 +248,13 @@ def interp_right_derivative_on_prepared_grid(
         fp_slopes: Node derivatives, NaN-padded in lockstep.
 
     Returns:
-        Right derivatives with the shape of `x_query`.
+        Tuple of the right-finiteness flag and the first, second, and third
+        right derivatives, each with the shape of `x_query`.
 
     """
     # Identical bracket location to `interp_on_prepared_grid`: `side="right"`
     # puts an on-node query into the node's right bracket, which is exactly the
-    # bracket whose derivative the right-continuous read needs.
+    # bracket whose germ the right-continuous read needs.
     upper = jnp.clip(
         jnp.searchsorted(search_grid, x_query, side="right"),
         1,
@@ -274,37 +279,52 @@ def interp_right_derivative_on_prepared_grid(
     safe_df = jnp.where(jnp.isfinite(df), df, 0.0)
     secant = safe_df / safe_width
 
-    # The same limiter as `_hermite_correction`, so this derivative is the
-    # derivative of exactly the polynomial the value read evaluates.
+    # The same limiter as `_hermite_correction`, so these derivatives are the
+    # derivatives of exactly the polynomial the value read evaluates. In the
+    # bracket's local coordinate `t` the read is
+    # `p(t) = f_l + Δf t + c_l t + (c_u - 2 c_l) t² + (c_l - c_u) t³`.
     def limit(slope: FloatND) -> FloatND:
         same_sign = slope * secant > 0.0
         limited = jnp.sign(secant) * jnp.minimum(jnp.abs(slope), 3.0 * jnp.abs(secant))
         return jnp.where(same_sign, limited, 0.0)
 
     coeff_lower = safe_width * limit(slope_lower) - safe_df
     coeff_upper = safe_df - safe_width * limit(slope_upper)
-    hermite_slope = (
+    hermite_first = (
         safe_df
         + (1.0 - 2.0 * relative_position)
         * ((1.0 - relative_position) * coeff_lower + relative_position * coeff_upper)
         + relative_position * (1.0 - relative_position) * (coeff_upper - coeff_lower)
     ) / safe_width
+    hermite_second = (
+        2.0 * (coeff_upper - 2.0 * coeff_lower)
+        + 6.0 * (coeff_lower - coeff_upper) * relative_position
+    ) / safe_width**2
+    hermite_third = 6.0 * (coeff_lower - coeff_upper) / safe_width**3
     applicable = (
         (bracket_width > 0.0)
         & jnp.isfinite(fp_lower)
         & jnp.isfinite(fp_upper)
         & jnp.isfinite(slope_lower)
         & jnp.isfinite(slope_upper)
     )
-    derivative = jnp.where(applicable, hermite_slope, secant)
+    first = jnp.where(applicable, hermite_first, secant)
+    second = jnp.where(applicable, hermite_second, 0.0)
+    third = jnp.where(applicable, hermite_third, 0.0)
     # The read clamps to a constant strictly below the first node and at or
-    # above the last valid one; its right derivative there is exactly zero.
-    # (`search_grid` is `+inf` on the pad, so a poisoned all-NaN row lands on
-    # the lower clamp and stays zero.)
+    # above the last valid one; the germ there is right-finite with all
+    # derivatives exactly zero. (`search_grid` is `+inf` on the pad, so a
+    # poisoned all-NaN row lands on the lower clamp.)
     first_node = search_grid[0]
     last_node = search_grid[jnp.maximum(valid_length - 1, 0)]
     on_clamp_ray = (x_query < first_node) | (x_query >= last_node)
-    return jnp.where(on_clamp_ray, 0.0, derivative)
+    right_finite = on_clamp_ray | (jnp.isfinite(fp_lower) & jnp.isfinite(fp_upper))
+    return (
+        right_finite,
+        jnp.where(on_clamp_ray, 0.0, first),
+        jnp.where(on_clamp_ray, 0.0, second),
+        jnp.where(on_clamp_ray, 0.0, third),
+    )
 
 
 def _interp_between_nodes(