refac: Use consistent "outer" default alignment for all arithmetic operations

[FBumann]

Breaking Change

This PR changes the default coordinate alignment for all arithmetic operations from a shape-dependent heuristic to a consistent "inner" join (intersection of coordinates). This matches xarray's own arithmetic_join default. Users who relied on the old positional matching for same-shape operands with different coordinates will see reduced or empty results instead. The old behavior is available via explicit join="override".

Convention

All arithmetic (+, -, *, /) defaults to "inner" (intersection of coordinates).

This is the only choice that:

1. Is safe: only shared coordinates appear in results — no silent fill values
2. Makes mistakes visible: disjoint coordinates produce empty results, immediately flagging errors
3. Is consistent: the same default applies regardless of whether operands happen to have the same shape
4. Preserves associativity: (a + b) + c == a + (b + c) regardless of operand shapes
5. Matches xarray: xarray's own arithmetic_join default is "inner"

No fill values are needed with inner join — every position in the result has data from both operands.

Constraints with DataArray RHS default to "left" — the expression defines where variables exist; missing RHS values become NaN (masked out).

Context	Default	Rationale
All arithmetic (+, -, *, /)	"inner"	Safe, consistent, matches xarray
Constraint with DataArray RHS	"left"	Expression defines variable positions; extra RHS coords → trivial constraints
Constraint with expression RHS	"inner"	Goes through arithmetic sub, inherits the arithmetic default

What changes from the status quo

The primary behavioral change is when operands have the same shape but different coordinates:

Scenario	Before (heuristic)	After
x(i=[0,1,2]) + z(i=[5,6,7])	"override": 3 entries, positional	"inner": 0 entries, empty (visible error)
x(i=[0,1,2]) + c(i=[5,6,7])	"override": 3 entries, positional	"inner": 0 entries, empty (visible error)
x(i=[0,1,2]) * c(i=[5,6,7])	"override": 3 entries, positional	"inner": 0 entries, empty (visible error)

When operands have subset/superset relationships (e.g., x has 20 coords, subset has 2):

• Before: result had 20 entries (filled with zeros/ones)
• After: result has 2 entries (intersection only)

When operands have matching coordinates: no change.

Available join modes

For explicit control, use .add(), .sub(), .mul(), .div(), .le(), .ge(), .eq() with a join parameter:

join	Coordinates	Use case
"inner" (default)	Intersection only	Safe default, no fill needed
"outer"	Union with fill values	When you want all positions
"left"	Left operand's	Keep expression's positions
"right"	Right operand's	Keep other operand's positions
"override"	Left operand's (positional)	Ignore labels, match by position
"exact"	Must match exactly	Raises error if different

Why this is necessary

The previous heuristic (check_common_keys_values) chose alignment based on operand shapes:

• Same shape → "override" (positional matching, ignoring labels)
• Different shape → "outer" (for expr+expr) or "left" (for expr+constant)

This caused two bugs:

1. Silent positional matching: x(i=[0,1,2]) + z(i=[5,6,7]) matched by position even though coordinates were entirely different, producing wrong results under x's labels
2. Non-associative addition: (y + factor) + x != y + (x + factor) because "left" for expr+constant dropped the constant's extra coordinates before they could be recovered by a subsequent expr addition

Changes

linopy/expressions.py

• merge(): Remove check_common_keys_values heuristic; when join is None, always use "inner". Pre-pad helper dimensions (_term) when concatenating along non-helper dimensions to ensure inner join only affects coordinate dimensions.
• _align_constant(): Default to "inner"; use xr.align for non-override/non-left joins
• _apply_constant_op(): Fill NaN coefficients with 0 before applying op (prevents NaN propagation from reindex fill at new positions)
• _multiply_by_constant(): Use fill_value=1 (identity element) for explicit outer/left joins
• _divide_by_constant(): Use fill_value=1 (identity element) for explicit outer/left joins
• to_constraint(): Handle DataArray RHS directly (without sub) to preserve intentional NaN masking for constraint positions; supports join parameter for all join modes

linopy/variables.py

• Variable.__mul__(): Use to_linexpr()._multiply_by_constant(other) instead of to_linexpr(other) so a * factor and (1 * a) * factor produce identical results

test/test_linear_expression.py

• Update TestSubsetCoordinateAlignment for inner defaults (subset operations give intersection)
• Update test_linear_expression_sum for inner default on disjoint slices (empty result)
• Rename TestAssociativity tests for inner semantics; add test_outer_gives_union
• Update TestJoinParameter quadratic tests for inner default on cross-product expressions

test/test_optimization.py

• Update test_non_aligned_variables for inner default (uncovered variables have no constraint)

examples/coordinate-alignment.ipynb

• Default is now "inner" for all arithmetic, "left" for constraint DataArray RHS
• Explain disjoint coordinates produce empty results (visible error)
• Show join="outer" as explicit opt-in for union behavior
• Update summary table

Checklist

• Code changes are sufficiently documented
• Unit tests for new features were added
• A note for the release notes doc/release_notes.rst of the upcoming release is included
• I consent to the release of this PR's code under the MIT license

[FabianHofmann]

thanks @FBumann, just FYI, the shape dependent overwrite was introduced as a convention and the documentation should state is clearly, but I also understand that this is not very intuitive and user might be confused. Since changing the behavior is very breaking (in pypsa for example we rely in the overwrite mechanism) I would like to make a clear transition path to the new behavior. So, we could pull in #572 and then add a global options tag in linopy to allow for new arithmetic alignments which users to opt in. default then in next major version, likely v2 ?

[FBumann]

thanks @FBumann, just FYI, the shape dependent overwrite was introduced as a convention and the documentation should state is clearly, but I also understand that this is not very intuitive and user might be confused. Since changing the behavior is very breaking (in pypsa for example we rely in the overwrite mechanism) I would like to make a clear transition path to the new behavior. So, we could pull in #572 and then add a global options tag in linopy to allow for new arithmetic alignments which users to opt in. default then in next major version, likely v2 ?

First of all, this PR is only a sketch up.
Im planing on opening another one with a better convention, inspired by @pyoframe.
About the positional alignment: I agree, and it IS convenient.
My issue is mainly that it could happen anytime by accident and is no opt in for this operation. This feels brittle.
I also rely on such, so we definitely need to keep it, but id like it as an opt in.
Ignore this PR for now please