Memory layout of arrays

[wkearn]

numpy arrays can have either row- or column-major layouts. We currently handle these layouts inconsistently, and this is an ongoing source of confusion.

Here is how things should behave:

• libtopotoolbox functions should be agnostic to the memory layout of their arguments.
• libtopotoolbox functions should iterate in "column-major" order. Really this means iterating over the fastest-changing dimension first, and otherwise the dimensions shouldn't matter too much. You do have to swap the order of the dimensions if you pass a row-major array.
• dem.shape should return (rows, columns) for every array, regardless of its memory order, just as numpy arrays.
• dem.dims should return (rows, columns) for column-major ('F') arrays and (columns, rows) for row-major ('C') arrays.
• pytopotoolbox should pass dem.dims rather than dem.shape to libtopotoolbox.
• StreamObject and FlowObject refer to their nodes with 0-based linear indices in column-major order.

The one major exception here is for graphflood, because I don't think we were very careful about matching it to the other libtopotoolbox functions. Correct me if I am wrong, @bgailleton, but I think it works kind of the same way as the other libtopotoolbox functions, but the dim arguments are reversed relative to the libtopotoolbox -- (rows, columns) for row-major and (columns, rows) for column-major.

I am noticing that the addition of the dims argument did not propagate to every libtopotoolbox function in pytopotoolbox.

I have added subissues for tasks that we need to fix these problems.

[bgailleton]

@wkearn Yes,, you're right about graphflood. It took a lot of trial and error to get it working, so I probably mixed dims up at some points.

Making the changes isn't too hard, just quite time-consuming. I'll try to set aside some time to update the C code. The main thing is to make sure the iterations in C are optimized. For example, if we're iterating over a column-major array in a row-major way, it can slow things down. C is typically row-major (right?), so we need to keep that in mind. But if we set everything to column-major, it won't be an issue.

[wkearn]

The main thing is to make sure the iterations in C are optimized. For example, if we're iterating over a column-major array in a row-major way, it can slow things down. C is typically row-major (right?), so we need to keep that in mind. But if we set everything to column-major, it won't be an issue.

If the array is laid out in column-major order, then the libtopotoolbox dims argument should be {rows, columns}, and we iterate like

ptrdiff_t dims[2] = {rows, columns};

for (ptrdiff_t j = 0; j < dims[1]; j++) {
    // j is the column index
    for (ptrdiff_t i = 0; i < dims[0]; i++) {
        // i is the row index
        float z = dem[j * dims[0] + i];
    }
}

So we iterate over the rows first, then the columns, which is the faster iteration order for column-major arrays.

On the other hand, if the array is laid out in row-major order, we pass dims = {columns, rows}. We use the same loop,

ptrdiff_t dims[2] = {columns, rows};

for (ptrdiff_t j = 0; j < dims[1]; j++) {
    // j is now the row index
    for (ptrdiff_t i = 0; i < dims[0]; i++) {
        // i is now the column index
        float z = dem[j * dims[0] + i];
    }
}

but we are iterating over the columns first then the rows, which is again the faster iteration order for row-major arrays. The meaning of the i and j indices has changed -- mathematically, you might write $Z_{ji}$, but performance-wise everything should be the same.

We generally don't use the C multidimensional array indexing notation dem[i][j], so the fact that C arrays are typically stored as row-major shouldn't matter. To us, arrays are linear sequences of elements, and we interpret them multidimensionally by indexing with j * dims[0] + i.

The main difficulties with doing things this way are:

1. It is very confusing.
2. All the arrays we pass to a single invocation of a libtopotoolbox function have the same ordering. Passing a row-major input array and a column-major output array will cause problems when we try to access the output in Python.
3. libtopotoolbox does not usually have geographic referencing available to it. It can only reference data in a coordinate system given by the dimensions that we pass in. This can be a problem when doing something like hillshade (Add hillshade to GridObject #164) or aspect, which needs to be referred to a geographic coordinate system.

Making the changes isn't too hard, just quite time-consuming. I'll try to set aside some time to update the C code.

I can also take a look. Some test cases would be helpful to make sure I'm not breaking anything. I can set up some snapshot tests with the DEMs in our TopoToolbox/snapshot_data repository, but other test cases that you know exercise various corners of graphflood would also be helpful.

[bgailleton]

OK thanks, that's clearer.

About the difficulties:

1. It is confusing, but I guess a clear exemple/documentation page dedicated to the latter should suffice? Contributor to the C side of the code should be realtively aware of that kind of issues.
2. and 3. if the input is row major and linear, the output calculated from that should also be row major, or am I getting confused? Apart of course for geographic data, a flag and a set of utility function could help (e.g. get_left_neighbour). Maybe a global flag telling libttb if the code is called from a row major or a column major language?

I will add the test cases to my todo list, I can produce an analytical testing framework and an array calculated from my other library to reproduce - but no idea when I can dedicate time to do that.