Breaking: Implement AbstractArray interface for AbstractDimStack

[tiemvanderdeure]

See #870

[rafaqz]

Lots of method ambiguities to work through.

Possibly we can get even more out of this by making it an AbstractBasicDimArray ? They are like AbstractDimArray but don't necessarily have an AbstractArray parent (Like DimIndices).

[tiemvanderdeure]

Yes I was already working on that. It mean I had to add the DimTuple to the type signature of AbstractDimStack, but I don't see why that should be a problem. This is breaking anyway.

[tiemvanderdeure]

except for all the ambiguities the only tests that are broken are in broadcast_dims. Honestly it's much less bad than I expected

[tiemvanderdeure]

Okay I actually managed to get rid of all ambiguities, I think. Most of it I could just delete, but I also explicitly specified that the eltype of a AbstractDimArray is a NamedTuple, because Aqua was complaining about ambiguous dispatches on AbstractVector{<:AbstractArray}. I don't know if there are every cases where this is not the case, though?

I ran some of the tests locally but let's see what they say.

[rafaqz]

Yeah its always a NamedTuple so fine to limit it.

One day we can make a MathyNamedTuple that + etc work on with a scalar. Like a mixed-type static array.

[tiemvanderdeure]

I think this should work for the most part now, but would like to suggest one more breaking change.

Currently the way DimArrays are treated in broadcast_dims is that it broadcasts for each of the layers, or throws an error if multiple stacks are broadcasted but have different layers. This is inconsistent with the behaviour of a normal broadcast, and not at all what I would expect when I broadcast an AbstractArray{<:NamedTuple}. I think we should just change it to it broadcasting over NamedTuples. I also think the current behaviour is undocumented.

[rafaqz]

The problem is broadcasting over NamedTuple is a very inefficient and annoying way to e.g. multiply all layers of a stack by a scalar, especially if you have mixed size layers.

That change will probably break mask or something similar in Rasters.

(There are reasons it's not AbstractArray ;)

Trying the Rasters tests may help clarify this

[tiemvanderdeure]

The problem is broadcasting over NamedTuple is a very inefficient and annoying way to e.g. multiply all layers of a stack by a scalar,

Wouldn't that already require something like maplayers(l -> l .* 2, s). Isn't that already preferred over broadcast_dims(x -> x * 2, s)?

(There are reasons it's not AbstractArray ;)

I know, but I'm positively surprised that deleting a bunch of methods and some small tweaks made tests pass.

Trying the Rasters tests may help clarify this

I'll give it a try and see how much of is broken on this branch

[rafaqz]

Ok maybe a little clearer is multiplying by another raster that is a subset of the dimensions, often a bitmask.

broadcast_dims means the same code works for a Raster and a RasterStack but that will break with your approach.

We will at least need to provide a function that does that.

Maybe or a bylayer keyword is a solution?

[tiemvanderdeure]

Trying the Rasters tests may help clarify this

I'll give it a try and see how much of is broken on this branch

For whatever it's worth: I got methods and stack tests to pass just by changing a few maps to maplayers (which I think we had to do anyway?)

[tiemvanderdeure]

Okay with that example I see the problem more clearly. A keyword is an option but also makes a pretty basic and intuitive function more complicated.

On this branch I think the easiest way to solve a case like this is

st_ma = maplayers(l -> (@d ma .* l), st)

where st is a stack and ma a mask. I think this is not awful, but also not great.

[rafaqz]

Except maplayers won't work on a Raster (so we lose the current generality). Maybe we make it so that just returns a Raster?

bylayer will also be a keyword elsewhere, so not the worst addition either

[tiemvanderdeure]

Both those changes work well together, I think? So maplayers would take AbstractDimStacks and/or AbstractDimArrays. If there no stacks at all then that would just be maplayers(f, A...) = f(A...) and maybe throw a warning. And if it's mixed then we wrap the arrays.

Then we can define broadcast_dims with bylayer = true to be just maplayers(As -> broadcast_dims(f, As...), st).

[rafaqz]

Why a warning?

[tiemvanderdeure]

We don't have to, I was thinking it could be worth warning because the function doesn't do anything in that case. Anyway, I gave it a shot and you can tell me what you think about that implementation.

Multiple dispatch got complicated because AbstractDimStack is a subtype of AbstractBasicDimArray, so I just used if-else logic, which still should compile away.

[rafaqz]

Right. We just don't want warnings because we'll actually use this on Rasters.

[tiemvanderdeure]

Aargh what are these Plots segfaults/compilation errors. So annoying

[tiemvanderdeure]

What do you think here @rafaqz? The only thing here that I ended up actually breaking at the user end is the behaviour of broadcast_dims, where you can get the old behaviour by setting bylayer=true. And I think this has some nice upsides like vec(stack) just works and stack .= f(stack) just works.

[rafaqz]

There are still some conceptual problems, like mean(stack), sum(stack) is even worse for mixed size layers.

[tiemvanderdeure]

There are still some conceptual problems

But we should we try to figure those out? I'm sure you have a much better understanding than I do.

like mean(stack), sum(stack) is even worse for mixed size layers.

How so worse? The implementation didn't change at all so I think those should be the exact same?

[rafaqz]

sum on an AbstractDimStack is currently calculated by layer. But as an AbstractArray some layers are repeated over missing dimensions so sum over the indexed values would be a much larger number (mean will still have the same answer in contrast). It's the same now but as our own object that's fine, it's just what it is.

But if we switch to AbstractArray we have to lose the mixed dimensionality in a lot of contexts to be consistent. And I'm not sure it's worth it. Instead we can implement broadcast on AbstractDimStack to be over all layers and always return a mixed size stack like the original. And other things like that, because we define the behaviour of the object rather than following an existing definition.