fix(piecewise): declare ragged breakpoints with mask= instead of inferring NaN

doc/release_notes.rst

@@ -35,6 +35,7 @@ Upcoming Version
 
 * Default internal integer labels to ``int32``, cutting memory ~25% and speeding up model build 10-35%. Models exceeding the int32 maximum (~2.1 billion labels) widen to ``int64`` automatically with a ``UserWarning``; pass ``Model(dtypes={"labels": np.int64})`` upfront to avoid the mid-build upcast (exposed read-only via ``Model.dtypes``).
 * ``add_variables(binary=True, ...)`` now accepts ``lower``/``upper`` bounds, as long as they are 0 or 1. Previously binary bounds could only be set via the ``.lower``/``.upper`` setters after creation. (https://github.com/PyPSA/linopy/issues/776)
+* ``add_piecewise_formulation`` gained a ``mask`` parameter declaring which breakpoint slots hold a real breakpoint. It is needed for **ragged** curves — entities with different numbers of breakpoints — which are stored densely with the surplus slots left absent. Under v1 that absence must be declared (``mask=x_pts.notnull()``) rather than read off the NaN padding. (https://github.com/PyPSA/linopy/issues/884)
 * ``add_piecewise_formulation`` gained an ``active_fill`` parameter that gates a partial ``active`` (defined over a subset of the indexed dimension, or masked) as always-active (``1``) or always-off (``0``); without it, a partial ``active`` — which was previously zeroed silently — now raises. Useful when one formulation mixes gated and ungated entities (e.g. committable and non-committable units sharing a ``status``). ``active_fill`` is transitional and will be removed once v1 semantics make ``active.reindex(coords).fillna(value)`` sufficient. (https://github.com/PyPSA/linopy/issues/796)
 
 *Documentation*

linopy/piecewise.py

@@ -50,12 +50,13 @@
     EvolvingAPIWarning,
     sign_replace_dict,
 )
+from linopy.semantics import check_user_nan_breakpoints
 
 if TYPE_CHECKING:
     from linopy.constraints import Constraint, Constraints
     from linopy.expressions import LinearExpression
     from linopy.model import Model
-    from linopy.types import LinExprLike
+    from linopy.types import LinExprLike, MaskLike
     from linopy.variables import Variables
 
 logger = logging.getLogger(__name__)
@@ -736,12 +737,17 @@ def _tangent_lines_impl(
     x: LinExprLike,
     x_points: BreaksLike,
     y_points: BreaksLike,
+    piece_mask: DataArray | None = None,
 ) -> LinearExpression:
     """
     Chord-expression math — the body of ``tangent_lines`` without the
     :class:`EvolvingAPIWarning`.  Called internally by ``_add_lp`` so a
     single ``add_piecewise_formulation((y, y_pts, "<="), (x, x_pts))``
     emits exactly one warning, not two.
+
+    ``piece_mask`` marks the pieces that exist.  A padded piece of a ragged
+    curve gets a neutral zero chord, which the caller is expected to mask
+    out of the resulting constraint.
     """
     from linopy.expressions import LinearExpression as LinExpr
     from linopy.variables import Variable
@@ -753,15 +759,15 @@ def _tangent_lines_impl(
     dy = y_points.diff(BREAKPOINT_DIM)
     piece_index = np.arange(dx.sizes[BREAKPOINT_DIM])
 
-    slopes = _rename_to_pieces(dy / dx, piece_index)
+    slopes = _drop_absent(_rename_to_pieces(dy / dx, piece_index), piece_mask)
     x_base = _rename_to_pieces(
         x_points.isel({BREAKPOINT_DIM: slice(None, -1)}), piece_index
     )
     y_base = _rename_to_pieces(
         y_points.isel({BREAKPOINT_DIM: slice(None, -1)}), piece_index
     )
 
-    intercepts = y_base - slopes * x_base
+    intercepts = _drop_absent(y_base - slopes * x_base, piece_mask)
 
     if not isinstance(x, Variable | LinExpr):
         raise TypeError(f"x must be a Variable or LinearExpression, got {type(x)}")
@@ -830,6 +836,10 @@ def tangent_lines(
         "entirely with "
         '`warnings.filterwarnings("ignore", category=linopy.EvolvingAPIWarning)`.',
     )
+    if bool(_coerce_breaks(x_points).isnull().any()) or bool(
+        _coerce_breaks(y_points).isnull().any()
+    ):
+        check_user_nan_breakpoints()
     return _tangent_lines_impl(x, x_points, y_points)
 
 
@@ -983,6 +993,58 @@ def _paired_valid_points(*points: DataArray) -> DataArray:
     return points[0].where(~invalid)
 
 
+def _drop_absent(values: DataArray, mask: DataArray | None) -> DataArray:
+    """Zero out masked-absent slots; the variable's own mask already excludes them."""
+    if mask is None:
+        return values
+    return values.where(mask, 0.0)
+
+
+def _resolve_breakpoint_mask(
+    bp_list: list[DataArray], mask: MaskLike | None
+) -> tuple[list[DataArray], DataArray | None]:
+    """
+    Settle which slots hold a real breakpoint, and align the values to it.
+
+    A ragged curve is stored densely along ``BREAKPOINT_DIM`` with the
+    surplus slots absent. A shorter curve and a data error are
+    indistinguishable from NaN alone, so the caller declares the absence
+    with ``mask=``; without a declaration v1 raises and legacy keeps
+    inferring from NaN placement. The returned arrays carry ``NaN``
+    exactly where the mask says absent, so ``isnull()`` stays the single
+    absence predicate downstream.
+    """
+    combined_null = bp_list[0].isnull()
+    for bp in bp_list[1:]:
+        combined_null = combined_null | bp.isnull()
+
+    if mask is None:
+        if not bool(combined_null.any()):
+            return bp_list, None
+        check_user_nan_breakpoints()
+        return bp_list, ~combined_null
+
+    declared = (mask if isinstance(mask, DataArray) else DataArray(mask)).astype(bool)
+    bp_dims = set().union(*(bp.dims for bp in bp_list))
+    extra = set(declared.dims) - bp_dims
+    if extra:
+        raise ValueError(
+            f"`mask` is not broadcastable against the breakpoints: it has "
+            f"dimension(s) {sorted(map(str, extra))} which the breakpoints "
+            f"({sorted(map(str, bp_dims))}) do not have."
+        )
+    _, declared = xr.broadcast(bp_list[0], declared)
+
+    if bool((combined_null & declared).any()):
+        raise ValueError(
+            "NaN at a breakpoint slot that `mask` marks as present. Either "
+            "widen the mask to cover those slots or fix the values with "
+            "`.fillna(...)`."
+        )
+
+    return [bp.where(declared) for bp in bp_list], declared
+
+
 def _validate_shared_coords(points: Sequence[DataArray]) -> None:
     skip = {BREAKPOINT_DIM, SEGMENT_DIM} | set(HELPER_DIMS)
     for i, left in enumerate(points):
@@ -1074,6 +1136,7 @@ def add_piecewise_formulation(
     method: PWL_METHOD = "auto",
     active: LinExprLike | None = None,
     active_fill: int | None = None,
+    mask: MaskLike | None = None,
     name: str | None = None,
 ) -> PiecewiseFormulation:
     r"""
@@ -1183,6 +1246,16 @@ def add_piecewise_formulation(
         Transitional convenience: under v1 semantics, pad ``active``
         explicitly with ``active.reindex(coords).fillna(value)`` instead —
         this parameter is slated for removal then.
+    mask : MaskLike, optional
+        Which breakpoint slots hold a real breakpoint — ``True`` where one
+        exists, ``False`` where it is absent.  Shaped like the breakpoint
+        arrays, or anything that broadcasts against them.  Needed for
+        **ragged** curves (entities with different numbers of breakpoints),
+        stored densely with the surplus slots left absent: under v1 that
+        absence must be declared rather than read off the NaN padding, e.g.
+        ``mask=x_pts.notnull()``.  Absent slots get no auxiliary variable
+        and no constraint.  A mask hiding a present value drops that
+        breakpoint; a NaN at a slot the mask calls present raises.
     name : str, optional
         Base name for generated variables/constraints.
 
@@ -1319,10 +1392,7 @@ def add_piecewise_formulation(
     _validate_shared_coords(bp_list)
     _validate_expr_coords(bp_list, lin_exprs)
 
-    combined_null = bp_list[0].isnull()
-    for bp in bp_list[1:]:
-        combined_null = combined_null | bp.isnull()
-    bp_mask = ~combined_null if bool(combined_null.any()) else None
+    bp_list, bp_mask = _resolve_breakpoint_mask(bp_list, mask)
 
     if name is None:
         name = f"pwl{model._pwlCounter}"
@@ -1705,13 +1775,15 @@ def _add_sos2(
     )
 
     if links.eq_expr is not None and links.eq_bp is not None:
-        input_weighted = (lambda_var * links.eq_bp).sum(dim=dim)
+        eq_bp = _drop_absent(links.eq_bp, links.bp_mask)
+        input_weighted = (lambda_var * eq_bp).sum(dim=dim)
         model.add_constraints(
             links.eq_expr == input_weighted, name=f"{name}{PWL_LINK_SUFFIX}"
         )
 
     if links.signed_expr is not None and links.signed_bp is not None:
-        output_weighted = (lambda_var * links.signed_bp).sum(dim=dim)
+        signed_bp = _drop_absent(links.signed_bp, links.bp_mask)
+        output_weighted = (lambda_var * signed_bp).sum(dim=dim)
         _add_signed_link(
             model,
             links.signed_expr,
@@ -1796,6 +1868,7 @@ def _add_incremental(
     def _incremental_weighted(bp: DataArray) -> LinearExpression:
         steps = bp.diff(dim).rename({dim: piece_dim})
         steps[piece_dim] = piece_index
+        steps = _drop_absent(steps, delta_mask)
         # ``drop=True`` keeps the breakpoint coord from sticking around as a
         # scalar on ``bp0_term`` — otherwise §11 rejects it as an aux-coord
         # conflict against the constraint LHS.
@@ -1877,13 +1950,15 @@ def _add_disjunctive(
     )
 
     if links.eq_expr is not None and links.eq_bp is not None:
-        input_weighted = (lambda_var * links.eq_bp).sum(dim=[SEGMENT_DIM, dim])
+        eq_bp = _drop_absent(links.eq_bp, bp_mask)
+        input_weighted = (lambda_var * eq_bp).sum(dim=[SEGMENT_DIM, dim])
         model.add_constraints(
             links.eq_expr == input_weighted, name=f"{name}{PWL_LINK_SUFFIX}"
         )
 
     if links.signed_expr is not None and links.signed_bp is not None:
-        output_weighted = (lambda_var * links.signed_bp).sum(dim=[SEGMENT_DIM, dim])
+        signed_bp = _drop_absent(links.signed_bp, bp_mask)
+        output_weighted = (lambda_var * signed_bp).sum(dim=[SEGMENT_DIM, dim])
         _add_signed_link(
             model,
             links.signed_expr,
@@ -1939,7 +2014,7 @@ def _add_lp(
 
     # Use the internal impl so we don't fire a second EvolvingAPIWarning —
     # ``add_piecewise_formulation`` already warned on entry.
-    tangents = _tangent_lines_impl(x_expr, x_points, y_points)
+    tangents = _tangent_lines_impl(x_expr, x_points, y_points, piece_mask)
     _add_signed_link(
         model,
         y_expr,

linopy/semantics.py

@@ -122,6 +122,52 @@ def _legacy_nan_constant_message(op_kind: str) -> str:
     )
 
 
+_RAGGED_RESOLVE_HINT = (
+    "\n  Ragged:    declare which slots hold a breakpoint —"
+    "\n             `add_piecewise_formulation(..., mask=x_pts.notnull())`"
+    "\n  Data error: `.fillna(value)`"
+)
+
+
+def _user_nan_breakpoint_message(context: str) -> str:
+    """
+    User-NaN error text for breakpoint tables.
+
+    The generic §5 message points at `mask=` on a *variable*, which is the
+    wrong remedy here — breakpoints are plain data, and raggedness is the
+    reason a NaN usually shows up. Name that case instead.
+    """
+    return (
+        f"NaN found in user-supplied {context} (§5). linopy trusts NaN only "
+        "from its own structural operations (§4): here it cannot tell a "
+        "shorter curve (a padded slot) from a data error, and guessing would "
+        "silently change the model." + _RAGGED_RESOLVE_HINT
+    )
+
+
+def _legacy_nan_breakpoint_message(context: str) -> str:
+    """Legacy inferred padding from NaN placement; v1 wants it declared."""
+    return (
+        f"NaN in user-supplied {context} was silently read as a padded slot "
+        "(an absent breakpoint) by legacy. Under v1 this raises ValueError — "
+        "the raggedness has to be declared rather than inferred."
+        + _RAGGED_RESOLVE_HINT
+        + _OPT_IN_HINT
+    )
+
+
+def check_user_nan_breakpoints(*, context: str = "breakpoints") -> None:
+    """
+    Enforce §5 for a breakpoint table whose padding was not declared.
+
+    v1 raises; legacy warns and lets the caller infer the padding from NaN
+    placement, preserving pre-v1 behaviour.
+    """
+    if is_v1():
+        raise ValueError(_user_nan_breakpoint_message(context))
+    warn_legacy(_legacy_nan_breakpoint_message(context))
+
+
 def _legacy_coord_mismatch_message(
     context: str,
     dim: str | None = None,