Improve documentation on overriding coordinates feature

[apfelix]

Hi,

I'm currently stuck trying to archive the following rather specific problem in linopy:

import linopy
m = linopy.Model()

# 2 days, hourly (48 hours)
t = list(range(1,49))
n = ["a", "b", "c", "d"]

level = m.add_variables(coords=[t, n], name="level", dims=["time", "node"])

change = m.add_variables(coords=[t, n], name="change", dims=["time", "node"])

# previous_level should by cyclic per day, so have time indices
#  24, 1, 2, ..., 23, 48, 25, 26, ..., 47
# |---- 1st day ----|----  2nd day   ----|
# constraint should be valid for each node
level == previous_level + change

So problem is in the creation of previous_level: I don't know have to do that. I know there is Variable.roll() for creating a shift in specific dimensions, but I basically want this operation but restricted to a subset of the time index. My plan was to use the result of roll() and assign the correct variable value to the positions 0 and 24 myself, but assignment does not seem to work. Is there any way to create a self-defined variable vector that has the form of the result of roll() (single variables per coordinate tuple), but using an arbitrary assignment?

PS: My fallback solution would be to create single constraints per point in time, but I would like to avoid that if possible ^^

[apfelix]

After some more time spend on this, I found at least two solutions/work arounds for this. Both create a corresponding linopy.LinearExpression instead of a linopy.Variable, but that's fine for me.

The first one uses the Pyomo-style LinearExpression.from_rule constructor. However, it has to be called for each n and does not use any vectorised operations.

import numpy as np
# create array [24, 1, 2, ..., 23, 48, 25, 26, ..., 47]
prev_t_map = np.array(range(48))
prev_t_map[0] += 24
prev_t_map[24] += 24

def previous_level_rule(m, this_t, this_n):
    prev_t = prev_t_map[this_t - 1]
    return 1.0 * level[prev_t, this_n]

previous_level_1 = linopy.LinearExpression.from_rule(m, previous_level_rule, level.coords)

The second one is kind of complicated as it builds a copy of the time dimension (with different name) and creates a square coefficient matrix between both dimensions, multiplies this onto the level, sums over the old time dimension (both operation could be replaced by the newly supported @ operator), and then renames the remaining copied dimension. The problem is the explicit creation of the square matrix, which requires a large amount of memory for larger dimension sizes. Also previous_level_2 still contains a lot of variables with zero value coefficients.

import xarray as xr
def get_row(i):
    row = np.zeros(48)
    row[prev_t_map[i] - 1] = 1.0
    return row

coeffcients = xr.DataArray(
    [get_row(i) for i in range(48)], coords=[t, t], dims=["time2", "time"]
)
previous_level_2 = (coeffcients * level).sum(dims="time").rename({"time2": "time"})

After calling densify_terms() on previous_level_2 both give the same result:

>>> print(previous_level_1.equals(previous_level_2.densify_terms()))
True

Both solutions do not really satisfy me due to the mentioned drawbacks. Also the execution time for large dimensions compared to a simple level.roll(time=1) is really large.

So I would still appreciate any tips on how to improve the construction of the model 🙂

[FabianHofmann]

Hey @apfelix, have you played around with .groupby("time.day") (this requires to have a date time index) with a subsequent rolling operation?

Something like

previous_level = level.groupby("time.day").roll(time=1)

Or you can create the index first and add the reindexed level

ind = [24, 1, ...]
level == change + level.loc[ind]

(None of this is tested, but might help anyway)

[apfelix]

@FabianHofmann Thanks for the fast response! level.loc[ind] was exactly what I was looking for, thanks a lot! It's also reasonably fast compared to standard .roll(time=1).

However, I'm a little confused on why this works to be honest ^^ When I print level.loc[prev_t_map] it gives:

Variable (time: 48, node: 4)
----------------------------
[24, a]: level[24, a] ∈ [-inf, inf]
[24, b]: level[24, b] ∈ [-inf, inf]
[24, c]: level[24, c] ∈ [-inf, inf]
[24, d]: level[24, d] ∈ [-inf, inf]
[1, a]: level[1, a] ∈ [-inf, inf]
[1, b]: level[1, b] ∈ [-inf, inf]
[1, c]: level[1, c] ∈ [-inf, inf]
		...
[46, b]: level[46, b] ∈ [-inf, inf]
[46, c]: level[46, c] ∈ [-inf, inf]
[46, d]: level[46, d] ∈ [-inf, inf]
[47, a]: level[47, a] ∈ [-inf, inf]
[47, b]: level[47, b] ∈ [-inf, inf]
[47, c]: level[47, c] ∈ [-inf, inf]
[47, d]: level[47, d] ∈ [-inf, inf]

So, the array value positions are as I need them, but the coordinates in front is also adjusted. I guess the coordinates are also the reason that previous_level_1.equals(level.loc[prev_t_map]) gives False. However, when using it in an equation as
change + level.loc[prev_t_map] this works (even though the coordinates are different):

LinearExpression (time: 48, node: 4):
-------------------------------------
[1, a]: +1 change[1, a] + 1 level[24, a]
[1, b]: +1 change[1, b] + 1 level[24, b]
[1, c]: +1 change[1, c] + 1 level[24, c]
[1, d]: +1 change[1, d] + 1 level[24, d]
[2, a]: +1 change[2, a] + 1 level[1, a]
[2, b]: +1 change[2, b] + 1 level[1, b]
[2, c]: +1 change[2, c] + 1 level[1, c]
		...
[47, b]: +1 change[47, b] + 1 level[46, b]
[47, c]: +1 change[47, c] + 1 level[46, c]
[47, d]: +1 change[47, d] + 1 level[46, d]
[48, a]: +1 change[48, a] + 1 level[47, a]
[48, b]: +1 change[48, b] + 1 level[47, b]
[48, c]: +1 change[48, c] + 1 level[47, c]
[48, d]: +1 change[48, d] + 1 level[47, d]

I guess that the coordinates are taken from the first of the terms, as switching the order in the sum results in a LinearExpression with coordinates as in level.loc[prev_t_map]:

[24, a]: +1 level[24, a] + 1 change[1, a]
[24, b]: +1 level[24, b] + 1 change[1, b]
[24, c]: +1 level[24, c] + 1 change[1, c]
[24, d]: +1 level[24, d] + 1 change[1, d]
[1, a]: +1 level[1, a] + 1 change[2, a]
[1, b]: +1 level[1, b] + 1 change[2, b]
[1, c]: +1 level[1, c] + 1 change[2, c]
		...
[46, b]: +1 level[46, b] + 1 change[47, b]
[46, c]: +1 level[46, c] + 1 change[47, c]
[46, d]: +1 level[46, d] + 1 change[47, d]
[47, a]: +1 level[47, a] + 1 change[48, a]
[47, b]: +1 level[47, b] + 1 change[48, b]
[47, c]: +1 level[47, c] + 1 change[48, c]
[47, d]: +1 level[47, d] + 1 change[48, d]

When closely reading the documentation, it even mentions this overwrite behaviour in https://linopy.readthedocs.io/en/latest/creating-expressions.html#Using-.loc-to-select-a-subset
I just didn't expect this, as in the cases I experienced during my short time of using linopy/xarray, coordinates are usually aligned and broadcasted.
Edit: A little above I now even found a large red box labeled "Important" explaining this, so I should have found this myself ...

Anyway, thanks for the help and for the general work you do on linopy! I really like the library and am especially grateful that you respond and fix things so quickly! Nice job!

[FabianHofmann]

Thanks for the nice feedback @apfelix :) happy to hear that. You are right about making this overriding feature clearer to the user base.

From your direct experience now, where would you have expected this quite important info? We could think about a runtime logger info the first time executing such an operation. Or a more prominent spot in the doc?

[apfelix]

Good question. I'm not sure about the logger, as it would produce just noise for those who are aware. Also it wouldn't have helped me, as you need to accidentally use the feature to get the corresponding notification explaining it.

Wrt to documentation, I would have two proposals:

1. Make the part below the "Important"-box an separate section in between "Arithmetic Operations" and "Using .loc to select a subset", to increase its visibility.
2. Next to the .shift explanation, add an explanation to .roll (as its pretty similar), and then show how the same effect as .roll could be achieved with the .loc[ind] + coordinate override functionality. So if you want to do something similar as .shift or .roll, you already have an alternative recipe listed there to achieve this (even though that may be too specifically target at my personal use case ^^).

[FabianHofmann]

very good, thanks for the feedback! I will keep this issue then open (renamed it) and try to tackle that asap. If in the meanwhile you would like to make a PR, feel free to do so. We are happy about every contribution.

[apfelix]

I would be happy to adjust the documentation (and also incorporate the result of the discussion in #257 ) 🙂 However, I will probably not have time today and next week is fully blocked, so I can only do that the week after.

[FabianHofmann]

Sounds great! I will make a PR that fixes the remaining inconsistency discussed in #257

[FBumann]

#591

[FBumann]

I know this is very late, but this serves as documentation:

The inconsistencies you observed (alpha * other_time_level vs alpha * other_time_level.to_linexpr() giving different results) come from the legacy "override" behavior — when two objects share a dimension with the same size but different coordinate order, linopy silently aligns by position. Which object's coords "win" depends on the operation type and operand order, leading to the confusing differences you found. This is still the case on current master.

The recommended pattern for .sel() reindexing is to reassign the target coordinates afterward:

other_times = [2, 0, 3, 1]
other_time_level = level.sel(time=other_times).assign_coords(time=t)

Coords now match — plain arithmetic just works

import linopy                                                                                                                                                                   
m = linopy.Model()

t = list(range(1, 49))
level = m.add_variables(coords=[t], name="level", dims=["time"])
change = m.add_variables(coords=[t], name="change", dims=["time"])

# Per-day cyclic shift: previous hour, wrapping within each day
# Day 1: [1..24], Day 2: [25..48]
prev_t = [24 if h == 1 else 48 if h == 25 else h - 1 for h in t]

previous_level = level.sel(time=prev_t).assign_coords(time=t)
m.add_constraints(level == previous_level + change, name="balance")
print(m.constraints["balance"]) 

Constraint `balance` [time: 48]:
--------------------------------
[1]: +1 level[1] - 1 level[24] - 1 change[1]    = -0.0
[2]: +1 level[2] - 1 level[1] - 1 change[2]     = -0.0
[3]: +1 level[3] - 1 level[2] - 1 change[3]     = -0.0
[4]: +1 level[4] - 1 level[3] - 1 change[4]     = -0.0
[5]: +1 level[5] - 1 level[4] - 1 change[5]     = -0.0
[6]: +1 level[6] - 1 level[5] - 1 change[6]     = -0.0
[7]: +1 level[7] - 1 level[6] - 1 change[7]     = -0.0
		...
[42]: +1 level[42] - 1 level[41] - 1 change[42] = -0.0
[43]: +1 level[43] - 1 level[42] - 1 change[43] = -0.0
[44]: +1 level[44] - 1 level[43] - 1 change[44] = -0.0
[45]: +1 level[45] - 1 level[44] - 1 change[45] = -0.0
[46]: +1 level[46] - 1 level[45] - 1 change[46] = -0.0
[47]: +1 level[47] - 1 level[46] - 1 change[47] = -0.0
[48]: +1 level[48] - 1 level[47] - 1 change[48] = -0.0

.sel() picks the variables by label, .assign_coords() stamps the target coordinates back on so everything aligns naturally. This works on current master.

We're working on a v1 arithmetic convention (#591) that fixes the root cause: without .assign_coords(), v1 raises a ValueError when coordinates don't match on a shared dimension — instead of silently overriding positionally. This would have caught the inconsistencies in this issue at the point of operation, rather than producing silently wrong results