PermutedDimsArray etc.

[mcabbott]

This treats PermutedDimsArray exactly like permutedims.

And it adds a method unname(A::NamedDimsArray, names) which is always a view of the parent, with axes ordered to match the given names.

[codecov]

Codecov Report

Merging #64 into master will not change coverage.
The diff coverage is 100%.

@@           Coverage Diff           @@
##           master      #64   +/-   ##
=======================================
  Coverage   86.42%   86.42%           
=======================================
  Files           8        8           
  Lines         221      221           
=======================================
  Hits          191      191           
  Misses         30       30

Impacted Files	Coverage Δ
src/functions_dims.jl	93.33% <100%> (ø)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 41d55be...decf606. Read the comment docs.

[oxinabox]

Can we discuss unname in a seperate issue?
I'm much more hestitant about introducing new public API functions,
than about fixing things that already work for AbstractArrays.

And there is more to discuss on that one,

[oxinabox]

thanks for this.

I think adding PermutedDimsArray is probably a good idea.

though maybe we do want to think a little bit more.
Since we are now overloading a constructor to not construct its type.
Which cases issues.
(I will link to an unstream bug report in a moment.)

@iamed2 what do you think?

[nickrobinson251]

Thanks for opening this! I have some concerns about the proposed behaviour:

[nickrobinson251]

Shouldn't this return a PermutedDimsArray (wrapping a NamedDimsArray)?

Wouldn't we just want to add a constructor that allows perm to contain Symbols?

Since PermuteDimsArray (with Ints) already works here, don't we just want:

julia> nda = NamedDimsArray{(:w, :x, :y, :z)}(ones(10, 20, 30, 40));

julia> pda = PermutedDimsArray(nda, (1, 2, 4, 3));  # already works

julia> typeof(pda)
PermutedDimsArray{Float64,4,(1, 2, 4, 3),(1, 2, 4, 3),NamedDimsArray{(:w, :x, :y, :z),Float64,4,Array{Float64,4}}}

julia> size(pda)
(10, 20, 40, 30)

julia> pda == PermutedDimsArray(nda, (:w, :x, :z, :y));  # this is the method that should be added?
true

[mcabbott]

I was hoping PermutedDimsArray could return something with names you can refer to. Although @oxinabox raises the issue that it's weird not to return the constructor's type. Perhaps another objection is that if you care about the order of indices, then you are moving to positional not named code anyway, and should unwrap.

[iamed2]

What if we keep PermutedDimsArray returning a PermutedDimsArray, but add some methods on PermutedDimsArray{<:NamedDimsArray} to return/use names? That seems like the alternative option to me. The issue with that approach is that there might be a bunch of methods necessary for dealing with NamedDimsArray × PermutedDimsArray interactions. Could be prohibitively troublesome.

[mcabbott]

Couldn't you argue the same for Transpose? The only difference is that the function usually called is named transpose.

[iamed2]

Ah right, that might mean that we can just include PermutedDimsArray with existing transpose-handling code (with some augmentation).

[nickrobinson251]

I think Eric's suggestion is nice, but is too much work and maintainence for the gain

[nickrobinson251]

To be clear, I think this is an interesting case... but I find PermutedDimsArray(::NamedDimsArray, perm) -> NamedDimsArray{PermuteDimsArray} to be really surprising behaviour. I expect to get a PermuteDimsArray{NamedDimsArrays} which is what we get on master. I do think allowing perm to be names i.e Symbols make sense.

I'd be happy with special behaviour for PermuteDimsArray{NamedDimsArrays} if the maintainence cost is proportionate .

Just a personal opinion - open to alternative views :)

[mcabbott]

It does seem odd that Julia has transpose which makes Transpose but no lazypermutedims which makes PermutedDimsArray, so the name of the type does double duty. Maybe that's the core here? This PR wants to extend lazypermutedims, although it happens not to have that spelling.

Calling lazypermutedims is supposed to do the same as permutedims but without moving the data; the same functions should work on the output. Although (in Base) this relies on these working on AbstractArray, not Array.

It would seem very strange to me that functions like sum(A, dims=:k) would be expected to work on things not of the form NamedDimsArray{Storage}. Anyone writing a function should then have to know to look inside, and not just dispatch on ::NamedDimsArray. They don't have to do this with NamedDimsArray{Transpose} etc.

[oxinabox]

It would seem very strange to me that functions like sum(A, dims=:k) would be expected to work on things not of the form NamedDimsArray{Storage}

For the case of wrapper arrays, then unless they are also messing with dims, I expect them to delegate that kwarg straight to the inner array type. And if they are, I kind of expect in many cases (Inc PermuteDimsArray) they can do their messing about with it and passes them modified version to the wrapped array.
And so it just works,in theory.

Orthogonal features added via wrapper arrays should work no matter the order of the wrapping

[mcabbott]

(I seem to have deleted my comment by accident, sorry)

Orthogonal features added via wrapper arrays should work no matter the order of the wrapping

I like this vision, but for now it seems many things do care. Some examples:

julia> xy = NamedDimsArray{(:x, :y)}(rand(2,2)); 

julia> yx = NamedDimsArray{(:y, :x)}(rand(2,2));

julia> xy * transpose(yx) # correctly errors, as NamedDimsArray is outermost
ERROR: DimensionMismatch("Cannot take matrix product of arrays with different inner dimension names. (:x, :y) vs (:x, :y)")

julia> xy * PermutedDimsArray(yx, (2,1)) # silently falls back to AbstractArray
2×2 NamedDimsArray{(:x, :_),Float64,2,Array{Float64,2}}:
 0.571736  0.94319 
 0.379215  0.585958

julia> sum(PermutedDimsArray(yx, (2,1)), dims=:y) # doesn't know how to pass this along
ERROR: MethodError: no method matching iterate(::Symbol)

julia> NamedDims.names(Transpose(yx)) # doesn't look inside
(:_, :_)

Perhaps one could avoid these by always dispatching on some big Union{NDA, Adjoint{NDA}, Transpose{NDA}, ...} type, I don't know how messy that would be.

But without commutativity of wrapping, I guess deciding who goes outermost might be a hard problem. I suppose my implicit rule here was that PermutedDimsArray is just a storage detail, like Transpose. User code shouldn't care, so it doesn't need to be exposed.

[mcabbott]

OK, I trimmed this to just PermutedDimsArray, with your good suggestions.

That's interesting about default constructors etc, I didn't realise this was anything beyond a convention. I can't in 5 minutes make Tracker produce such errors with tracked PermutedDimsArrays, should that be possible?

[oxinabox]

I can't in 5 minutes make Tracker produce such errors with tracked PermutedDimsArrays, should that be possible?

julia> x = [10a + b for a in 1:3, b in 1:4];

julia> trx = TrackedArray(x);

julia> Set(PermutedDimsArray(data) for data in [trx])
ERROR: MethodError: no method matching PermutedDimsArray(::TrackedArray{…,Array{Int64,2}})

[mcabbott]

But that's just complaining about the lack of a permutation, there's no one-arg method:

julia> PermutedDimsArray(x)
ERROR: MethodError: no method matching PermutedDimsArray(::Array{Int64,2})

julia> Set(PermutedDimsArray(data, (2,1)) for data in [trx])
Set(TrackedArray{…,PermutedDimsArray{Int64,2,(2, 1),(2, 1),Array{Int64,2}}}[[11 21 31; 12 22 32; 13 23 33; 14 24 34] (tracked)])

[oxinabox]

Ah, right, yes. This probably can't happen for things that don't have 1arg constructors. Just because of the nature that bug.

Nvm

[iamed2]

@nickrobinson251's position sounds reasonable to me

[mcabbott]

Closing in favour of someday making wrappers outside NamedDimsArray work smoothly, instead.

[mcabbott]

I had a go at making wrappers outside NamedDimsArray work, over here: https://github.com/mcabbott/NamedPlus.jl . Some things work, but things like sum(A; dims=:a) don't yet look inside the outermost type.