intersection math on DimArray

[Beforerr]

I tried to mimic xarray behavorior for doing basic math (result only contains intersected coordinates)

import Base: -
-(a::DimArray, b::DimArray) = begin
    shared_selectors = intersect(DimSelectors(a), DimSelectors(b))
    invoke(-, Tuple{AbstractArray,AbstractArray}, a[shared_selectors], -b[shared_selectors])
end

However the shared_selectors does not keep dimension as original one. Is there anyway to do the basic math?

[rafaqz]

Basic math on a DimArray works the same as math on a base julia Array because DimArray <: AbstractArray and we get most of it for free.

Most likely math is more consistently and broadly implemented here than with xarray. Even functions from a third party julia packages we don't know about will usually work.

So - is already defined:

julia> a = rand(X(3), Y(4))
┌ 3×4 DimArray{Float64, 2} ┐
├──────────────────────────┴─────────────────────────────── dims ┐
  ↓ X, → Y
└────────────────────────────────────────────────────────────────┘
 0.356113  0.848334  0.60161    0.401837
 0.722693  0.546814  0.668194   0.0447609
 0.502677  0.820793  0.0470841  0.657557

julia> b = rand(X(3), Y(4))
┌ 3×4 DimArray{Float64, 2} ┐
├──────────────────────────┴─────────────────────────────── dims ┐
  ↓ X, → Y
└────────────────────────────────────────────────────────────────┘
 0.371906  0.209425  0.971638  0.263029
 0.290974  0.489066  0.827357  0.997566
 0.704057  0.636763  0.480227  0.365513

julia> a - b
┌ 3×4 DimArray{Float64, 2} ┐
├──────────────────────────┴─────────────────────────────── dims ┐
  ↓ X, → Y
└────────────────────────────────────────────────────────────────┘
 -0.0157936  0.63891    -0.370027   0.138808
  0.431719   0.0577473  -0.159164  -0.952805
 -0.20138    0.184029   -0.433143   0.292044

Additionally, its usually a bad idea to change methods like these yourself.

[rafaqz]

What you are talking about is custom operations that don't really map to Julia AbstractArray interface or the usual meaning of array math.

In Julia we try not to change the semantics of base methods or cause any ambiguities with what would happen to a regular array (because this means DimArray will work everywhere a regular array will without bugs). What I would do instead, and we could even add to the package is add something like an intersecting(-, a, b) that does what you want for arbitrary functions, with the semantics of your intersection rather than breaking base math semantics.

[rafaqz]

You are losing dimensions because of intersect function returns a vector.

This does a bit better:

function intersecting(f, a::DimArray, b::DimArray)
    shared_selectors = DimSelectors(a)[DimSelectors(b)]
    data = f(a[shared_selectors], b[shared_selectors])
end

Then:

julia> intersecting(-, a, b)
3×4 Matrix{Float64}:
 -0.0157936  0.63891    -0.370027   0.138808
  0.431719   0.0577473  -0.159164  -0.952805
 -0.20138    0.184029   -0.433143   0.292044

But selecting DimSelectors from DimSelectors and then selecting with the result loses our dimensions because we don't know the selectors are ordered. We can make indexing DimSelectors with DimSelectors maintain dimensions - its just not implemented but isn't complicated to do.

[Beforerr]

Hi Rafael, thank you for your thoughtful response regarding the implementation of custom operations.
I agree with your idea of maintaining the semantics of base methods and yes I prefer the syntax intersecting here. It would be great to see related intersection function implemented

[rafaqz]

Great. I will first look at fixing DimSelectors/DimIndices etc so the output is a DimArray through multiple selections, then we can add something like intersecting.

Can you think of a better name? That's really the hardest thing here.

[Beforerr]

I would suggest align_apply as this name clearly conveys the two-step process: first aligning the arrays to their common coordinates (intersection) and then applying the specified operation.

[rafaqz]

Ah good point. map is a little more julian than apply for this kind of thing, so maybe mapaligned ?

[Beforerr]

Maybe amap like pmap or more verbosedalign_map? Put aligned after map is a little bit strange.

[rafaqz]

Ok good ideas. Will see what seems best in the PR

[felixcremer]

I am not entirely convinced that we need a special function for that. If this is only going to be two lines of code that might also be properly documented in a how to and not added as a special function. Especially since there are the keyword arguments to the two DimSelectors that we would have to deal with.

[rafaqz]

The biggest thing is the DimSelector etc indexing chains properly so the above code returns a DimArray

But I'm not sure everyone will be able to figure that out even when it works? I don't know if lazy
generators of selectors are a common concept... Having a function could be a bit more accessible to people.

[rafaqz]

The first problem was just a minor bug, fixed in #915

Then we can see if we need the function too.

But the output is now:

julia> intersecting(-, a, b)
┌ 3×4 DimArray{Float64, 2} ┐
├──────────────────────────┴─────────────────────────────── dims ┐
  ↓ X, → Y
└────────────────────────────────────────────────────────────────┘
 -0.537523   -0.768356   0.284044   0.269776
 -0.0890164  -0.535661  -0.246372  -0.647103
 -0.656387    0.370359  -0.681565  -0.266145