RFC: Drop dimensions indexed by scalars

[mbauman]

This is the first step towards APL indexing semantics. It does not yet allow indexing with multiple multi-dimensional arrays (the new feature with the catchy "rank of the result is the sum of the rank of the indices" semantic); it simply drops all scalar dimensions (a major breaking change).

There is an interesting interaction between these array semantics and transposes that I didn't fully appreciate until looking at this more in-depth. The tricky part about APL indexing is that it loses the orientation of row vectors, and so we can no longer define complex matrix multiplication in terms of slices directly:

• This makes A[1,:] * B[:,1] an illegal (n,) × (n,) operation.
• Similarly, dot(A[1,:], B[:,1]) and A[1,:]' * B[:,1] erroneously conjugate the row.

Working with complex row slices becomes a little dicier. When do you conjugate? Or (c)transpose?

Fortunately, the behavior here is a strict superset of our current indexing rules. You can explicitly request that the shapes of the slices be maintained in both dimensions: A[1:1,:] * B[:,1:1] if you wish to use the result as a linear algebra construct.

[mbauman]

AppVeyor failure looks like #9572.

[IainNZ]

Very cute

[malmaud]

🐨

[andreasnoack]

@mbauman Would you rebase this one. I think we should consider merging this soon such that people will start using it.

[mbauman]

Done.

[mschauer]

Thanks for keeping it rolling here

[andreasnoack]

Okay to merge this?

[mbauman]

Sure, let's give it a shot.

[timholy]

Sorry I've not had time to review this. Looks fine to me. As you noted, it depends on a lot more splatting (esp in checksize), but all that is O(1) and not O(N) in the array elements. So I'm fine with seeing how it goes (and it might be further incentive to fix the splatting penalty), though if you've run any benchmarks it might be interesting to learn more about any performance hits.

[yuyichao]

Should we also drop the scalar dimension for SubArray (especially if we'd like to return SubArray as the result of indexing).

This breaks the equivalence between getindex and sub

julia> a = rand(2, 2)
2x2 Array{Float64,2}:
 0.377862   0.85667 
 0.0177678  0.695366

julia> a[2, 1:2] == sub(a, 2, 1:2)
false

[timholy]

slice already does that---essentially, now getindex is equivalent to slice rather than sub.

[mschauer]

Are there plans to adapt mapslices to slice rather than to sub?

[andreasnoack]

I don't think it has been discussed. Can you explain further?

[mschauer]

Currently, mapslices and special cased reductions over dimensions keep the reduced dimension as singleton dimension. As the indexing behavior now changed I thought mapslices and friends should follow soonish. This would mean changing the standard behavior of mapslices to dimension dropping as hinted at by

julia/base/abstractarray.jl

Line 1217 in dc5c974

# TODO: maybe support removing dimensions

[timholy]

I'm not sure I agree that those two design decisions are coupled. Meaning, I don't think it's inconsistent for getindex to drop dimensions but reductions to retain them: they are different function calls, and reductions are obviously different from indexing.

It's so common to say

Xnorm = X ./ maximum(X, 2)

that it would be a shame to make that any harder. (And dropping dimensions definitely does make that harder.) When you want to drop the dimensions, it's easier than ever to do so.

[StefanKarpinski]

I agree – I would prefer that reductions keep reduced dimensions for exactly this reason.

[mschauer]

I kind of agree for the monadic case, for a dyadic mapslices function the current scheme does not translate well. I now recall why I noted this, I have even RFC #13317 but I forgot about it. Maybe someone could comment there? Some feedback would be appreciated.

[stevengj]

I think this should be moved to the "Breaking Changes" part of the NEWS file. This just caused silently incorrect results in one of my packages (Hadamard.jl)

[davidavdav]

Hello,

For NamedArrays I need to copy the exact behaviour of the dropping of dimensions (for julia-0.4 and 0.5). The AbstractArray contained in a n::NamedArray is sliced by calling getindex(n.array, index...)---so there I get the same behaviour automatically. But for the names of the dimensions in n, I need to know exactly which types of index... I need to keep.

Is there a Base function that will tell me this? I was looking at Base.index_shape but that isn't defined for all index types (e.g., CartesianIndex). Currently I use a helper function

dimkeepingtype(x) = false
dimkeepingtype(x::Vector) = true
dimkeepingtype(x::Range) = true
dimkeepingtype(x::BitVector) = true

that tells me if an index type x (∈ index) is non-scalar or not. I am not sure this covers all index types in the same way that AbstractArray keeps dimensions, so I'd rather refer similar functionality in Base.

[timholy]

Can't you add the missing index_shape_dim method(s) to NamedArrays? For julia-0.5 we could add such a method to Base, though---probably a good idea.

[davidavdav]

Thanks. But I do not really understand the logic in index_shape_dim---it seems a recursive helper function for index_shape. It doesn't look like the way towards telling me which dimensions are kept in slices. I need these to know which dimension names I need to select, and from which dimensions I have to slice the names of the indices in the NamedArray.

Somewhat related, I saw that in Julia-0.5 the dimension of a sliced array is the sum of the dimensions of the indices. That is great, but I don't really know yet how to come up with dimension and index names, with cases like NamedArray[Array] and NamedArray[NamedArray]. Does anybody know of a precedence for this case? I only know of R named arrays, but R doesn't have this fancy indexing scheme.

[mbauman]

I don't know about precedence in other languages, but perhaps you can take some inspiration from AxisArrays:

• AxisArray[Array1, Array2] returns an AxisArray with axis names row_1, row_2, …, row_$(ndims(Array1)), and col_1, col_2, …, col_$(ndims(Array2)). The left-hand side of the _ comes from the names of the parent AxisArray.
• AxisArray[AxisArray] is similar, except instead of using numbers on the right-hand side of the _, it uses the axis names of the index. This is where the axis names time_sub and time_rep come from on the README example.

[davidavdav]

Thanks, that is indeed a useful way of naming the dimensions.

For the names of the indexes along the dimensions, though, things are more complicated. I think I will opt for copying the names of the indexes of the indexing array (if this has ndims() > 1), or defaulting to 1..size(index) if the indexing array itself doesn't have names.

Further challenges will occur when indexing a NamedArray with an AxisArray or vice versa. I don't even want to think about this at this moment...