Should every AbstractArray be iterated over in cloumn-major order only? In other words, what relies on iterate preceding in a set order for AbstractArrays?
To clarify, I mean iteration order in a very general sense. Could I make an array type that iterates over its values in an inward going spiral for example? (Without breaking assumptions somewhere)
rand(2,2)' is row major.
julia> a = rand(2,2)'
2×2 LinearAlgebra.Adjoint{Float64,Array{Float64,2}}:
0.258786 0.342334
0.381858 0.474121
julia> for i in eachindex(a) @show a[i] end
a[i] = 0.2587855392761107
a[i] = 0.3818578636243877
a[i] = 0.34233445327094136
a[i] = 0.47412107821130967
eachindex doesn’t seem to specialize on that, though I would think it should.
No, AbstractArray does not tell you anything about data layout. For example, transpose of Matrix, index-based views, or sparse matrices all claim conformance to AbstractArray, but have nonstandard (i.e. not column-major) layout.
Afaiu we don’t have a consensus holy-trait for querying the layout of types yet.
@Elrod did you ever add a column major trait to ArrayInterface.jl? I think that’s an interesting omission.
I ask because for arrays stored on disk, I’d like to iterate chunk by chunk, but I’m not sure if I can do that without implicitly breaking someones code in a way that’s quite hard to diagnose.
If we can also have some more general interface for iteration order, that would be awesome.
julia> using ArrayInterface, StaticArrays, Test
julia> using ArrayInterface: stride_rank
julia> @test @inferred(stride_rank(@SArray(zeros(2,2,2)))) === ArrayInterface.StrideRank((1, 2, 3))
Test Passed
julia> @test @inferred(stride_rank(A)) === ArrayInterface.StrideRank((1,2,3))
Test Passed
julia> @test @inferred(stride_rank(PermutedDimsArray(A,(3,1,2)))) === ArrayInterface.StrideRank((3, 1, 2))
Test Passed
julia> @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:]))) === ArrayInterface.StrideRank((1, 2))
Test Passed
julia> @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[2,1:2,:])')) === ArrayInterface.StrideRank((2, 1))
Test Passed
julia> @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[2:3,1:2,:]))) === ArrayInterface.StrideRank((3, 1, 2))
Test Passed
julia> @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:]))) === ArrayInterface.StrideRank((3, 2))
Test Passed
julia> @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[2:3,2,:])')) === ArrayInterface.StrideRank((2, 3))
Test Passed
julia> @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[:,1:2,1])')) === ArrayInterface.StrideRank((1, 3))
Test Passed
julia> @test @inferred(stride_rank(@view(PermutedDimsArray(A,(3,1,2))[:,2,1])')) === ArrayInterface.StrideRank((2, 1))
Test Passed
If it returns StrideRank((1,2,3,4,5)) for a 5-d array, that array is column major.
Should we add a convenience function for checking if the StrideRank is ordered?
To the initial question, I don’t know how much this behavior is relied upon. It is the behavior, though, and changing it will be tough/breaking — largely in user code. It’d take a survey to see how much this breaks — my hunch is that it’ll largely be in loops that iterate over i in eachindex(A) but use i as in index into other arrays. It’ll also apply to functions that manually iterate through (potentially multiple) arrays.
Ref https://github.com/JuliaLang/julia/issues/38284 for some more context.
I’d say yes.
That isn’t quite what I meant here though. You could have a matrix, that, when iterated iterates in an inward spiral, but still have eachindex return indices column by column.
As has become increasingly clear to me, iteration order and indexing are actually quite separate issues.
And if in fact iterate ought to iterate in some specific order, should this be added to the manual / the docs, to make this a clear part of the interface?
Is returning a Val{true}() for any increasing StrideRank and a Val{false}() otherwise fine?
For the op’s question, the StrideRank actually yields the correct answer: it gives you the order of the inner most to outer most loops if you want to iterate as close to contiguously as possible.
The example:
A = @view(PermutedDimsArray(A,(3,1,2))[2:3,1:2,:])
@test @inferred(stride_rank(A)) === ArrayInterface.StrideRank((3, 1, 2))
Means that we should have
for i in axes(A,1) # dim 1 has rank 3
for j in axes(A,3) # dim 3 has rank 2
for k in axes(A,2) # dim 2 has rank 1
end
end
end
Perhaps we should add something to help iterate in this order? But given that CartesianIndices are often slow, this would be as well.
(Note that the next major release of LoopVectorization will use StrideRanks of all accessed arrays to help choose orderings; currently it used assumes column major for everything except a few hard coded exceptions like Transpose, Adjoint, or PermutedDimsArray.)
I think so
So, an AbstractArray can specify how it ought to be iterated, at least in terms of specifying whether cloumns, rows etc. should be iterated in first, last, etc. using StrideRank but this still says nothing about what iterate is expected to do. @ChrisRackauckas original answer seems to suggest that, arrays are free to switch column / row ordering. But does this generalise to more complex iteration? What if I made an array type that iterated over its elements in random order (assuming every element was iterated over exactly once)? Is that something I could do without breaking assumptions? I realise this is a bit of a squishy question. Of course, some user somewhere could have written code that relies on iteration orrder. My question is really whether this is part of the interface. Should users expect a certain iteration order?
@Elrod If you’ll allow me going off topic for a sec, would it be possible to extend this concept of StrideRank to a more general interface specifying how arrays ought to be iterated? Say, to allow arrays to specify that they ought to be iterated in blocks without having to subtype an AbstractBlockArray, but using a trait instead.
When you have operations with more than one array, you won’t necessarily be able to access elements of all of them in order. As an obvious example:
A = rand(64,64);
A .+ A'
What is the best order here?
The best order there is probably something like 4x4 or 8x8 blocks of both arrays, which corresponds to the memory layout of neither array.
For any given array, there is no one how it ought to be iterated. How we ought to iterate across it depends on what we’re doing with it, and what we’re doing with other arrays at the same time.
This is why ArrayInterface doesn’t currently tell people what to do with the arrays. It just describes their internal memory layout in a way that LoopVectorization.jl can use (and will in the next major release) to pick what it actually ought to do in that particular context.
The behavior of the particular code is also what will determine if the iteration order matters.
But if someone writes eachindex(A), presumably they don’t care much about the order; the docs promise to [c]reate an iterable object for visiting each index of an AbstractArray A in an efficient manner.
But on specifying that an Array’s memory is stored in blocks…perhaps I should add that. Currently, it only has support for saying that it is stored in vectors:
julia> A = rand(4,4);
julia> ArrayInterface.contiguous_axis(A)
ArrayInterface.Contiguous{1}()
If this returns a dimension other than the one for which StrideRank is 1, then that means there are contiguous vectors embedded there, of length ArrayInterface.contiguous_batch_size(A).
For example, if we have, hypothetically, an array A such that
julia> ArrayInterface.stride_rank(A)
ArrayInterface.StrideRank{(1, 2)}()
julia> ArrayInterface.contiguous_axis(A)
ArrayInterface.Contiguous{2}()
julia> ArrayInterface.contiguous_batch_size(A)
ArrayInterface.ContiguousBatch{8}()
Then if A were a 2 x 16 matrix, the order of the indices in memory would be:
A[1,1]
A[1,2]
A[1,3]
A[1,4]
A[1,5]
A[1,6]
A[1,7]
A[1,8]
A[2,1]
A[2,2]
A[2,3]
A[2,4]
A[2,5]
A[2,6]
A[2,7]
A[2,8]
A[1,9]
A[1,10]
A[1,11]
A[1,12]
A[1,13]
A[1,14]
A[1,15]
A[1,16]
A[2,9]
A[2,10]
A[2,11]
A[2,12]
A[2,13]
A[2,14]
A[2,15]
A[2,16]
That is, we have sets of 8 elements from the second axis mixed in between.
Perhaps this should be generalized to blocks.
Thanks for that answer. I learned a few things there.
So, as to my specific question: since eachindex doesn’t even need to adhere to any specifc order, iterate doesn’t either. In other words, users should not assume that eachindex or iterate iterate in the same order for any AbstractArray. Is that summary correct?
+1 for adding a blocked storage layout trait. I’d be happy to help implement that, I just have no clue where to start 
. I want blocked memory layout as a trait so I can more easily make array stored on disk in chunks work. BlockArrays.jl requires subtyping, which isn’t ideal for me.
I would replace the current contiguous_axis and contiguous_batch_size with a version generalized to an arbitrary number of indices.
I’d be inclined to require the block sizes to be known at compile time, but if you’d prefer this not to be the case I would define them to return tuples of ArrayInterface.StaticInt for those known at compile time, and Int otherwise.
Let me know if you’d like to take a stab at it soon. Otherwise, I’ll get to it shortly.
I’d be happy to take a stab at it weekend after next if you haven’t gotten to it by then.