refac: Harmonize linopy operations and introduce a new predictable and strict convention

[FBumann]

Harmonize linopy arithmetic with legacy/v1 convention transition

Why: silent bugs in the current (legacy) arithmetic

linopy's current coordinate alignment has several classes of silent correctness bugs. None of them raise errors — they produce results that look reasonable but are quietly wrong.

1. Positional alignment ignores labels (#586, #550, #257)

When two operands have the same shape on a shared dimension, linopy matches them by position, ignoring coordinate labels entirely:

x = m.add_variables(coords=[["costs", "penalty"]], name="x")
factors = xr.DataArray([2.0, 1.0], dims=["effect"],
                       coords={"effect": ["penalty", "costs"]})
x * factors
# Result: x["costs"] * 2.0, x["penalty"] * 1.0  — SWAPPED!
# Expected: x["costs"] * 1.0, x["penalty"] * 2.0

This affects all arithmetic operators and constraint creation. The optimization solves without error, but with the wrong constraints.

2. Subset constants break associativity (#572)

When a constant has different-sized coordinates, a left-join drops coordinates that might be needed by a later operation:

# a has 3 time steps, b has 5, factor has 5
a + factor + b    # factor at time=3,4 → 40, 50  ✓
a + b + factor    # factor at time=3,4 → LOST     ✗

(a + factor) + b ≠ a + (b + factor) — the result depends on operand order.

3. User NaN silently swallowed (#620)

NaN in user-supplied data (often indicating missing or erroneous data) is silently filled with inconsistent neutral elements:

data = xr.DataArray([1.0, np.nan, 3.0], ...)
x + data    # NaN → 0 (silent)
x * data    # NaN → 0 (kills the variable!)
x / data    # NaN → 1 (leaves variable unchanged!)

The fill values differ across operations with no consistent principle.

4. Absent slots indistinguishable from zero (#620)

Variable.to_linexpr() does not mark absent variable slots (from .shift(), .where()) as NaN, so multiplication and fillna() cannot distinguish "absent" from "zero":

xs = x.shift(time=1)   # time=0 is absent
xs * 3                  # absent slot looks like 0 — not NaN
xs.fillna(42)           # no-op — nothing to fill

5. Inconsistent Variable vs Expression paths (#569, #571)

x * subset_constant and (1 * x) * subset_constant previously gave different results (one crashed, one didn't).

---

What: the v1 convention

This PR introduces linopy.options["arithmetic_convention"] with two modes and a non-breaking transition:

• "legacy" (default) — reproduces current master behavior exactly. Emits LinopyDeprecationWarning on every legacy codepath.
• "v1" — strict coordinate matching, explicit NaN handling. Silent bugs become loud errors.

linopy.options["arithmetic_convention"] = "v1"

Strict coordinate matching

Arithmetic operators (+, -, *, /) require matching coordinates on shared dimensions. Mismatched coordinates raise ValueError with suggestions:

x[i=0,1,2] + y[i=1,2,3]  # ValueError: use .add(y, join="inner")

Named methods with explicit join — .add(), .sub(), .mul(), .div(), .le(), .ge(), .eq() accept join= parameter:

x.add(y, join="inner")      # intersection
x.add(y, join="outer")      # union with fill
x.add(y, join="left")       # keep x's coordinates
x.add(y, join="override")   # positional alignment (opt-in only)

Named methods with fill_value — .add(), .sub(), .mul(), .div() accept fill_value= for explicit NaN filling before the operation (see PR #620):

expr.add(5, fill_value=0)      # fill NaN const with 0 before adding
expr.mul(factor, fill_value=1) # fill NaN const with 1 before multiplying

Free broadcasting — Constants can introduce new dimensions without restriction.

Algebraic laws

All standard algebraic laws hold under v1 and under legacy for same-coordinate operands:

• Commutativity: a + b == b + a, a * c == c * a
• Associativity: (a + b) + c == a + (b + c)
• Distributivity: c * (a + b) == c*a + c*b, (a + b) / c == a/c + b/c
• Identity, negation, zero

Legacy breaks associativity only when operands have mismatched coordinate ranges (the subset-constant case above). The v1 convention prevents this by requiring explicit coordinate handling.

---

v1 NaN convention

NaN means "absent term" — never a numeric value.

NaN enters only from mask= at construction or structural operations (.shift(), .where(), .reindex(), .reindex_like(), .unstack()). Operations like .roll(), .sel(), .isel() do not produce NaN.

Arithmetic with absent slots:

Operation	Absent slot	Why
shifted + 5	+5	Fills const with 0 (additive identity), required for associativity
shifted - 5	-5	Same — subtraction is addition of negation
shifted * 3	absent	NaN propagates — no correct implicit fill (0 kills, 1 preserves)
shifted / 2	absent	Same as multiplication

When merging expressions (e.g., x + y.shift(time=1)), NaN marks individual terms, not entire coordinates. A coordinate is only fully absent when all terms are absent — isnull() checks this.

User-supplied NaN raises ValueError — users must handle NaN explicitly before arithmetic:

x + data.fillna(0)           # NaN = "no offset"
x * factor.fillna(1)         # NaN = "no scaling"
expr.mul(3, fill_value=0)    # fill_value= shorthand on named methods

fillna() — Variable.fillna(numeric) returns LinearExpression, Variable.fillna(Variable) returns Variable, Expression.fillna(value) fills const at absent slots.

---

Source changes

File	Change
config.py	LinopyDeprecationWarning, arithmetic_convention setting
expressions.py	All arithmetic paths branch on convention; merge() pre-validates user-dim coords under v1; merge() preserves NaN const when all inputs NaN; _add_constant fills const with additive identity; to_constraint has separate legacy/v1 paths; NaN validation at API boundaries; fill_value= on .add()/.sub()/.mul()/.div()
common.py	align() reads convention (legacy→inner, v1→exact)
variables.py	Scalar fast path in __mul__, explicit TypeError in __div__, .reindex() methods, Variable.fillna(numeric) returns LinearExpression, to_linexpr() sets const=NaN at absent slots in v1
piecewise.py	.fillna(0) on breakpoint data for v1 compatibility
monkey_patch_xarray.py	DataArray/Dataset arithmetic with linopy types
model.py	Convention-aware model methods

Documentation

Notebook	Content
arithmetic-convention.ipynb	Coordinate alignment rules, join parameter, migration guide
missing-data.ipynb	NaN convention principles, fillna patterns, masking with .sel() and mask=, legacy comparison
_nan-edge-cases.ipynb	Dev notebook: shift, roll, where, reindex, isnull, absent slot arithmetic, fillna/fill_value API, FILL_VALUE internals

Test structure

• Marker-based separation: @pytest.mark.v1_only and @pytest.mark.legacy_only for convention-specific tests
• Shared fixtures in conftest.py with auto-convention switching
• test_algebraic_properties.py (92 tests): formal specification of algebraic laws — commutativity, associativity, distributivity, identity, negation, zero, division/subtraction distributivity, multi-step constant folding, mixed-type commutativity, expression-expression laws. All pass under both conventions except 8 NaN-propagation tests (v1 only).
• test_legacy_violations.py (23 tests): catalog of concrete legacy bugs with paired legacy/v1 tests — positional alignment, subset associativity, user NaN handling, variable/expression inconsistency, absent-slot propagation. Each test traces to a specific issue number.
• Additional convention-specific tests in test_linear_expression.py, test_constraints.py, test_convention.py

Rollout plan

1. This PR: Default "legacy" — nothing breaks
2. Downstream: Users opt in with linopy.options["arithmetic_convention"] = "v1"
3. linopy v1: Flip default to "v1" and drop legacy mode

Open questions

• from_tuples / linexpr() — Currently follows the global convention. In practice always called with same-coord variables, so convention doesn't matter. Low-priority.
• Pipe operator — Only linopy objects, or also constants? (follow-up PR)

Sub-PRs

• Fix NaN propagation consistency and add fill_value API #620 — Fix NaN propagation consistency and add fill_value API

Test plan

• All tests pass under both conventions
• Legacy tests validate backward compatibility
• v1 tests validate strict coordinate matching and NaN raises
• Algebraic properties verified (92 tests): all laws hold under both conventions
• Legacy violations documented (23 tests): concrete bugs with issue references
• Piecewise and SOS tests run under both conventions
• Documentation notebooks execute cleanly

🤖 Generated with Claude Code

[FBumann]

@FabianHofmann Im quite happy with the notebook now. It showcases the convention and its consequences.
Tests need some work though. And migration as well.
Looking forward to your opinion on the convention

[FBumann]

The convention should be "exact" for all of +, -, *, /, with an additional check that neither side may introduce dimensions the other doesn't have — also for all operations.

Why "exact" instead of "inner" for * and /

"exact" still broadcasts freely over dimensions that only exist on one side — it only enforces strict matching on shared dimensions. So the common scaling pattern works fine:

cost = xr.DataArray([10, 20], coords=[("tech", ["wind", "solar"])])
capacity  # dims: (tech=["wind", "solar"], region=["A", "B"])

cost * capacity  # ✓ tech matches exactly, region broadcasts freely

"inner" is dangerous: if coords on a shared dimension don't match due to a typo or upstream change, it silently drops values. The explicit and safe way to subset before multiplying is:

capacity.sel(tech=["wind", "solar"]) * renewable_cost

No operation should introduce new dimensions

Neither side of any arithmetic operation should be allowed to introduce dimensions the other doesn't have. The same problem applies to + and - as to * and / — new dimensions silently expand the optimization problem in unintended ways:

cost_expr      # dims: (tech, time)
regional_expr  # dims: (tech, time, region)

cost_expr + regional_expr  # ✗ silently expands to (tech, time, region)

capacity  # dims: (tech, region, time)
risk      # dims: (tech, scenario)
risk * capacity  # ✗ silently expands to (tech, region, time, scenario)

An explicit pre-check on all operations:

asymmetric_dims = set(other.dims).symmetric_difference(set(self.dims))
if asymmetric_dims:
    raise ValueError(f"Operation introduces new dimensions: {asymmetric_dims}")

Summary

Operation	Convention
+, -, *, /	"exact" on shared dims; neither side may introduce dims the other doesn't have

[coroa]

The convention should be "exact" for all of +, -, *, /, with an additional check that neither side may introduce dimensions the other doesn't have — also for all operations.

Let's clearly differentiate between dimensions and labels.

labels

I agree with "exact" for labels by default, but we need an easy way to have inner or outer joining characteristics. I found the pyoframe conventions
strange at the beginning, but they grew on me:

x + y.keep_extras() to say that an outer join is in order and mismatches should fill with 0.

x + y.drop_extras() to say that you want an outer inner join.
x.drop_extras() + y does the same, though.

I have in a different project used | 0 to indicate keep_extras ie (x + y | 0).

dimensions

i am actually fond of the ability to auto broadcast over different dimensions. and would want to keep that (actually my main problem with pyoframe).

your first example actually implicitly assumes broadcasting.

[FBumann]

Dimensions and broadcasting

I agree that auto broadcasting is helpful in some cases.
I'm happy with allowing broadcasting of constants. We could allow this always...?
But I would enforce that the constant never has more dims than the variable/expression.
Or is there a use case for this?

So the full convention requires two separate things:
1. "exact" join — shared dims must have matching coords (xarray handles this)
2. Subset dim check — the constant side’s dims must be a subset of the variable/expression (custom pre-check needed)

labels

I'm not sure if I like this approach, as it's needs careful state management of the flags on expressions. The flag (keep or drop extras) needs to be handled.
I would rather enforce to reindex or fill data to the correct index.
I think aligning is the correct approach:

import linopy

# outer join — fill gaps with 0 before adding
x_aligned, y_aligned = linopy.align(x, y, join="outer", fill_value=0)
x_aligned + y_aligned

# inner join — drop non-matching coords before adding
x_aligned, y_aligned = linopy.align(x, y, join="inner")
x_aligned + y_aligned

Combining disjoint expressions would then still need the explicit methods though.
I'm interested about your take on this

[FBumann]

The proposed convention for all arithmetic operations in linopy:
1. "exact" join by default — shared coords must match exactly, raises on mismatch
2. Subset dim check — constants may introduce dimensions the variable/expression doesn’t have
3. No implicit inner join — use .sel() explicitly instead
4. Outer join with fill — use x + (y | 0) or .add(join="outer", fill_value=0)
The escape hatches in order of preference: .sel() for subsetting, | 0 for inline fill, named method .add(join=...) for everything else. No context manager needed.​​​​​​​​​​​​​​​​

I'm not sure how to implement the | operator yet. Might need some sort of flag/state for defered indexing

[FBumann]

I thought about the pipe operator:
I think it should only work with linopy internal types (Variables/expression), not constants (scalar, numpy, pandas, dataarray), as this would need monkey patching a lot and hard to get stable.

Would this be an issue for you?

[FBumann]

@brynpickering I worked on it a bit and treating nan as absent terms works quite well.
Regarding your examples, it removes the need for fillna(), except for the efficiency, where not calling fillna(1) silently masks out the constraint, which is intended, but a dangerous caveat.
Further, this would be a silent change in behaviour which would not raise.
See #627 for details.

I tried to address this, but its not that easy and makes the convention less clear. My resulting code in #627 pretty much only reduces the amount of fillna(0) calls in user code. But i thik its not work the complexity in the mental model and api complexity in linopy.

Im quite happy with this branch #591

[brynpickering]

@FBumann agree that pre-alignment is cleanest. Thanks for the example. It's effectively what we do in the calliope optimisation backend, which is to align all relevant arrays to consistent sets for each constraint and expression. It could cause memory spikes in cases where alignment creates a very large (usually very sparse) array, which might soon after be collapsed with an aggregation. However, that's more of an issue for the user to handle in how they define their math.

RE that second example. You're right, I can't reproduce it locally now either. Not sure what happened there.

[FBumann]

@brynpickering Thanks for the back and forth.
This really helps

I also ensure the alignment and indexing before working with linopy. This is a pain point which gets resolved a lot by the strictness of this convention. Linopy raises instead of guessing what might be correct

[MaykThewessen]

Downstream PyPSA-Eur-style consumer perspective. We run an 8760h × NL+neighbors LP with several custom constraint modules that all sit in the v1-arithmetic blast radius described in the PR body. Sharing concrete patterns in case they're useful as regression cases.

Patterns we'd test against v1:

1. NaN-mostly per-snapshot equality pins. We pin BESS terminal SoC via n.storage_units_t.state_of_charge_set as a DataFrame(NaN, columns=batteries) with exactly one non-NaN entry per chunk. PyPSA forwards this to linopy as a constraint RHS. Under the legacy "NaN silently swallowed → 0" rule (case 3 in the PR), the wrong semantics would force SoC=0 every hour. Today this works because PyPSA pre-filters NaN before constraint creation; under v1 we'd want NaN-as-absent to be a first-class signal so the pre-filter could be dropped. A regression test: build m.add_constraints(x == rhs) with rhs a DataArray that is NaN-everywhere-except-one-snapshot, confirm only one constraint row is materialized.

2. Per-snapshot time-varying marginal_cost_t with sparse non-NaN coverage. We override foreign offshore wind mc_t to equal landing-bus LMP only on snapshots where the gen would otherwise dispatch (NaN elsewhere → static base cost). Same NaN-handling concern as (1), but on the objective side rather than the constraint side. Currently works; v1 needs to keep working.

3. Label-indexed pandas Series multiplied by per-country/per-carrier DataFrames. Our scarcity-rent and neighbor-markup-propagation modules construct factors like pd.Series({'NL': 0.15, 'BE': 0.12, ...}) and multiply against per-snapshot generation expressions. This is exactly case 1 (positional alignment ignores labels) territory. Mixing pandas (label-aligned) and xarray (position-aligned) operands is where we feel least safe. Strict label alignment in v1 would be a clear win — we'd happily accept a ValueError on mismatch over silent swap.

4. phase_shift_extendable PST custom constraint (PyPSA fork PR #1661, currently in our pin). Writes Transformer phase-shift variables and flow constraints with mixed label/positional coords. Would benefit from a clear v0→v1 migration story — we'd rather port once, with strict errors guiding us, than chase silent dispatch shifts later.

5. Sweep-orchestrator multi-process reproducibility. Adjacent to fix: keep coords dimension order for DataArray bounds (#706) #710's _broadcast_points set-iteration fix: we run many variants in parallel processes. Any hash-randomized dim ordering produces apparent run-to-run drift that masks real model deltas. Strict v1 ordering is a correctness win even beyond the silent-bugs angle.

Migration concern: the breakage surface is large for downstream PyPSA-Eur derivatives. A deprecation cycle that emits a LinopyAlignmentWarning on any operation whose v0 vs v1 result differs — rather than blanket on all alignment — would let us find and port only the affected call sites. Otherwise the warning noise floor will be too high to act on.

Offer: happy to run our 8760h LP against a v1-enabled branch and report any silently shifted constraints (compare LP file hash + objective). Let us know which SHA you want exercised.

[FBumann]

@MaykThewessen Thanks for the perspective. We are working on this currently, but it will take a few weeks probably until we have it ready for you to run.

Would you be open for discussion about the nan stuff?

[MaykThewessen]

@FBumann yes, happy to. NaN semantics are the part that touches our code most
directly — constraint-RHS side (SoC pin via state_of_charge_set) and objective
side (marginal_cost_t and p_max_pu with sparse non-NaN coverage). Async on a
dedicated issue is probably most efficient; sync call also fine if it helps
unblock design.

Concrete starting point: we ran a controlled test of #627's draft against our
SoC-pin pattern this morning.

• Under arithmetic_convention="legacy" (default), behaviour is identical to
  current master — SoC pin holds (SoC[-1]=100.0 in a 24h LP).
• Under arithmetic_convention="v1", the same pattern raises
  ValueError: Constraint RHS contains NaN values at expressions.py:1325,
  before PyPSA's existing NaN-RHS auto-mask runs. Loud, not silent — which is
  great. But it means every PyPSA-Eur downstream that relies on the auto-mask
  (storage with state_of_charge_set, sparse marginal_cost_t, gapped
  p_max_pu) will fail at constraint-build time the moment they flip the
  option.

The design question is where the pre-mask should live. Two options:

1. Keep linopy strict, push mask to consumers. Every downstream
   pre-filters NaN before add_constraints. One rule, but high migration cost
   and duplicated masking logic across PyPSA, PyPSA-Eur, atlite-derived
   builders, custom forks.
2. add_constraints accepts NaN-as-absent natively, per-call. E.g.
   add_constraints(lhs == rhs, on_nan="raise" | "skip_row" | "propagate")
   with raise as the v1 default. One implementation, fused with the existing
   xarray alignment, no global mode, explicit at each call site.

I lean toward (2) because the NaN in question is created by linopy's own
xr.align(join="outer"), not by user intent — the consumer can't distinguish
"user meant skip" from "alignment fill", only linopy sees the align op. So
handling it in the alignment layer feels structurally correct. But I see the
argument for keeping linopy policy-free.

Happy to open a discussion issue with the repro + sketch of both options.
File on linopy, or prefer to scope on your side first?

[FBumann]

@FBumann yes, happy to. NaN semantics are the part that touches our code most directly — constraint-RHS side (SoC pin via state_of_charge_set) and objective side (marginal_cost_t and p_max_pu with sparse non-NaN coverage). Async on a dedicated issue is probably most efficient; sync call also fine if it helps unblock design.

Concrete starting point: we ran a controlled test of #627's draft against our SoC-pin pattern this morning.

• Under arithmetic_convention="legacy" (default), behaviour is identical to
  current master — SoC pin holds (SoC[-1]=100.0 in a 24h LP).
• Under arithmetic_convention="v1", the same pattern raises
  ValueError: Constraint RHS contains NaN values at expressions.py:1325,
  before PyPSA's existing NaN-RHS auto-mask runs. Loud, not silent — which is
  great. But it means every PyPSA-Eur downstream that relies on the auto-mask
  (storage with state_of_charge_set, sparse marginal_cost_t, gapped
  p_max_pu) will fail at constraint-build time the moment they flip the
  option.

The design question is where the pre-mask should live. Two options:

1. Keep linopy strict, push mask to consumers. Every downstream
   pre-filters NaN before add_constraints. One rule, but high migration cost
   and duplicated masking logic across PyPSA, PyPSA-Eur, atlite-derived
   builders, custom forks.
2. add_constraints accepts NaN-as-absent natively, per-call. E.g.
   add_constraints(lhs == rhs, on_nan="raise" | "skip_row" | "propagate")
   with raise as the v1 default. One implementation, fused with the existing
   xarray alignment, no global mode, explicit at each call site.

I lean toward (2) because the NaN in question is created by linopy's own xr.align(join="outer"), not by user intent — the consumer can't distinguish "user meant skip" from "alignment fill", only linopy sees the align op. So handling it in the alignment layer feels structurally correct. But I see the argument for keeping linopy policy-free.

Happy to open a discussion issue with the repro + sketch of both options. File on linopy, or prefer to scope on your side first?

Feel free to open one. I like your take on nan. Im leaning on defaulting to rais/strict, but allowing it via opt in seems reasonable. Im not sure if it should be in add_constraints though, because a user using pypsa cant reach that, and it would not work for expressions, as they are created bedore add_constraints() does its job.

[FBumann]

Superseeded by #717