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 As i understand your comment, you are mostly ok with the convention, but have thoughts on nan data. Am I right? THis is what im currently not that opinionated about, so im happy to hear your suggestions.

[FBumann]

@brynpickering Thanks for the detailed feedback.

Fill values depend on context (efficiency=1 vs cost=0) — this is id like to make it a user choice via .fillna(value) rather than picking a default.

Misalignment doesn't have to raise. The escape hatches are lightweight:

total = cost_a.add(cost_b, join="outer")  # absent terms → zero at solve time
scaled = var * efficiency.fillna(1)        # you choose the fill

As far as I know, Memory overhead is the same either way — fillna() on a parameter is not much larger than the parameter itself, and join="outer" does alignment inside xarray's concat anyway. But maybe im wrong here.

For expressions with mismatched labels, join="outer" already works: absent terms carry NaN coefficients that are filtered out when flattening for the solver. No explicit fill value needed.

Maybe adding a fill_value= parameter on .mul() / .add() resolves your issue if the fillna + join= pattern proves too verbose in practice?

[FBumann]

@coroa @FabianHofmann Im quite happy with this. But i think it desperately needs discussion.

[brynpickering]

@FBumann your fillna suggestion is fine for parameters, but doesn't address nan filling later down the line of variables or expressions. As far as I understand linopy objects, it isn't possible to nan fill in the same way. This risks NaNs making their way into arithmetic when one would rather have some numeric value that has no effect to be there.

[FBumann]

@brynpickering

I added a (temporary?) notebook to document the behaviour more clearly. It would be great if you could give it a look and check if you find sth inconsistent.
The mental model of nan in linopy objects isnt trivial imo. The internal representation is nan, but it doesnt behave like nan!

As far as i understand, it boils down to as single design question
Take this setup:

import pandas as pd
import linopy

linopy.options["arithmetic_convention"] = "v1"
m = linopy.Model()
time = pd.RangeIndex(5, name="time")
x = m.add_variables(lower=0, coords=[time], name="x")

shiftintroduces an empty slot. Should arithmetics revive this (with a fill value defined by linopy (zero))?

expr_=  x.shift(time=1) + 2
print(expr)

Revived

LinearExpression [time: 5]:
---------------------------
[0]: +2
[1]: +5 x[0] + 2
[2]: +5 x[1] + 2
[3]: +5 x[2] + 2
[4]: +5 x[3] + 2

nan propagation

LinearExpression [time: 5]:
---------------------------
[0]: None
[1]: +5 x[0] + 2
[2]: +5 x[1] + 2
[3]: +5 x[2] + 2
[4]: +5 x[3] + 2

Did I understand you correctly?

Currently, Option A (no propagation) is implemented

[FBumann]

@brynpickering I think the current state (updated description) should satisfy your concerns.