Do not require DenseArray

[orenbenkiki]

Right not the constructor of StridedView requires the base array to be DenseArray.
However there are types which support strides and are not DenseArray,
so would work in a StridedView - if it were not for the requirements.

How about removing the <:DenseArray from the signature? If the argument has strides, it "should" just work; if it does not, then it would fail with a more informative error message (missing strides) than today (not having a StridedView constructor).

Right now StridedView knows about some special cases (e.g. Transpose) which can be made strided w/o being DenseArray, but this doesn't work with non-Base types. I encountered this trying to pass PyArray from PythonCall to @strided, but this is just one example.

julia> using BenchmarkTools

julia> using Strided

julia> using PythonCall

julia> a = rand(4000, 4000);

julia> pa = PyArray(a);

julia> b = similar(a);

julia> pb = PyArray(b);

julia> @btime b .= (a .+ a) ./ 2;
  18.413 ms (4 allocations: 128 bytes)

julia> @btime pb .= (a .+ a) ./ 2;
  18.694 ms (4 allocations: 128 bytes)

julia> @btime @strided b .= (a .+ a) ./ 2;
  8.684 ms (17 allocations: 1.06 KiB)

julia> @btime @strided pb .= (pa .+ pa) ./ 2;
ERROR: MethodError: no method matching StridedView(::PyArray{Float64, 2, true, true, Float64})

[Jutho]

The problem is that, the way all the methods are implemented, is by assuming that in the end there is a dense memory region underlying the data, as all indices are in the end transformed into a linear index into this dense parent array. I guess I could add a generic fallback for AbstractArray that then tries to apply parent to see if this yields (possibly after several levels of recursion) a DenseArray.

I should check that would work with the PyArray implementation.

[orenbenkiki]

Why not simply call strides on the array? If it works, great, you can use your implementation. If it doesn't, the user gets a clear error that the array isn't strided, so it makes no sense to use @strided on it.

If you want to fall back to an implementation for any AbstractArray, that's also possible, but that's a separate feature.

[Jutho]

My apologies for not responding sooner. Maybe this is a good strategy nowadays. I think there used to be a time that strides did return something also for arrays that were not strided.

Still, in the end the current implementation needs to have an object in which it can index linearly using the strides to compute the linear index. Maybe that feels for some cases which are not DenseArray?

[orenbenkiki]

I don't see what strides can (meaningfully) return for arrays that are not strided. Being paranoid one can test that not only it returns (and not missing), that it also returns a tuple of the right amount of integers - otherwise, consider the array non-strided. If it does return a tuple of the right amount of integers, then it is pretty safe to say the array elements can be indexed linearly using the strides, right?

[Jutho]

My apologies for the late response. I am afraid that what you say is not true. Consider the following:

julia> A = randn(10,10);

julia> B = view(A, 1:2:10, 1:2:10); strides(B)
(2, 20)

julia> B[3,2], B[1 + (3-1)*2 + (2-1)*20], B.parent[1 + (3-1)*2+(2-1)*20]
(-1.3798464342130103, -1.7417661865291487, -1.3798464342130103)

julia> B[3,4], B[1 + (3-1)*2 + (4-1)*20]
ERROR: BoundsError: attempt to access 5×5 view(::Matrix{Float64}, 1:2:9, 1:2:9) with eltype Float64 at index [65]
Stacktrace:
 [1] throw_boundserror(A::SubArray{Float64, 2, Matrix{Float64}, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}}, false}, I::Tuple{Int64})
   @ Base ./abstractarray.jl:744
 [2] checkbounds
   @ ./abstractarray.jl:709 [inlined]
 [3] _getindex
   @ ./abstractarray.jl:1340 [inlined]
 [4] getindex(A::SubArray{Float64, 2, Matrix{Float64}, Tuple{StepRange{Int64, Int64}, StepRange{Int64, Int64}}, false}, I::Int64)
   @ Base ./abstractarray.jl:1296
 [5] top-level scope
   @ REPL[10]:1

As you can see, it is only because I can do linearly indexing in the parent array associated with the view, or more generally, in the memory associated with that array and strides, that I can use the strides to compute all the index information. Hence, I take DenseArray as a substitute for a continuous block of memory in which I can compute index information using strides. I should check the PyArray implementation, but I assume it has something similar, despite maybe not being declared as DenseArray.

[Jutho]

To add to that, we are loosening the restriction of DenseArray in the type definition, but still there will need to be something like a continuous block of memory in which linear indexing using stride calculations makes sense, in order for a custom array to be compatible with StridedViews.jl

[orenbenkiki]

I don't follow the code in the issue above (what's "B"?)

If the point is that an array can be strided, but still not contiguous, well, of course.
To be contiguous, one stride must be 1 element, another must be the dimension size, another the product of this dimension size by another dimension size, etc.

That is, if I have a 5x6 array, then strides (1, 5) and (6, 1) indicate a contiguous layout in memory, and any other strides will indicate a "Swiss cheese" array - e.g., (2, 10).

It is reasonable that some operations will only work for strided arrays that are also contiguous. Still, not all such arrays happen to be DenseArray. It still seems to me that the most robust way to check would be to call strides and examine the results.

Am I missing something?

[Jutho]

My example is indeed missing a line, my apologies for that. Not sure how that happend. I have now corrected it.

The point is that linear indexing into an array does not mean that you are accessing the underlying memory. Julia's indexing behaviour makes it such that you can do linear indexing into any AbstractArray as if it were an Array, i.e. as if it were a contiguous block of memory with the strides given by the cumulative product of the previous sizes.

There is no general way to directly index into the memory that backs this AbstractArray, for which one would compute a "linear memory index" using the actual stride information. For the Base wrapper types, this is only possible by recursively calling parent, until the underlying DenseArray is obtained.

[orenbenkiki]

I thought the whole point of the strides function was that it only existed in case one can compute the memory address of each element using the linear indexing using the strides. For any array where this isn't the case (e.g., sparse arrays or whatever), then the strides function is simply not implemented. All this has nothing to do with the fact that accessing the array using a single index still works.

Consider this session:

julia> using SparseArrays

julia> a = randn(3, 4)
3x4 Matrix{Float64}:
 -0.811094    0.205791   0.649853   -0.393014
  0.0505602   1.39507   -0.23117     0.158512
  0.838679   -0.235745  -0.0231283   0.623796

julia> a[a .< 0] .= 0
5-element view(::Vector{Float64}, ith eltype Float64:
 0.0
 0.0
 0.0
 0.0
 0.0

julia> a
3x4 Matrix{Float64}:
 0.0        0.205791  0.649853  0.0
 0.0505602  1.39507   0.0       0.158512
 0.838679   0.0       0.0       0.623796

julia> s = sparse(a)
3x4 SparseMatrixCSC{Float64, Int64} with 7 stored entries:
 .          0.205791  0.649853   .
 0.0505602  1.39507   .          0.158512
 0.838679   .         .          0.623796

julia> a[1]
0.0

julia> s[1]
0.0

julia> a[2]
0.05056023785407835

julia> s[2]
0.05056023785407835

julia> a[3]
0.8386790273652533

julia> s[3]
0.8386790273652533

julia> a[4]
0.2057912550575144

julia> s[4]
0.2057912550575144

julia> a[5]
1.39507

julia> s[5]
1.39507

julia> a[6]
0.0

julia> s[6]
0.0

julia> strides(a)
(1, 3)

julia> strides(s)
ERROR: MethodError: no method matching strides(::SparseMatrixCSC{Float64, Int64})
...

Yes, one can write both a[4] and s[4] and get the same (correct) result, but s is not strided while a is. That is, for a the implementation of a[4] would be trivial and efficient - the address of the element in memory is exactly the 4th consecutive element starting from the base array pointer; while for s the implementation of s[4] would be complex and inefficient as it needs to account for the unstored 0 values (the offset in memory would be 3, not 4).

Looking at the parent doesn't tell us anything because anyone is free to define a non-dense array type that has a DenseArray parent. For example, a hypothetical zero-copy view of a triangular matrix out of a full matrix would have a DenseArray parent, but would most definitely not be a strided array.

It seems to me that the only safe way to distinguish between the cases is by checking whether strides is implemented. If it is, we know that the memory address of elements by using a linear combination of the indices using the strides; if it is not implemented, then the array isn't strided. If the array is strided, we can further check whether the strides are such that the array is contiguous in memory or not (e.g., a zero-copy view of every other row of a DenseArray would be strided, but not contiguous).

[akriegman]

I have a different use case with a different solution. I was hoping to compose multiple StridedViews, similar to what they do in tinygrad. So to index into a composition of strided views, you would do outer multiindex -> linear index -> inner multiindex -> linear index, and then finally use that to index into memory. So in other words we could pretend that the parent is a dense block of memory, and in the cases where it actually is hopefully the compiler would recognize that. Sounds like maybe that's not the point of this package though.

[Jutho]

The point is that I still need to know, if I compute a linear index using the strides, what to index into. Consider the following:

julia> A = randn(5,5)
5×5 Matrix{Float64}:
 -1.83299     0.373954  -0.0752648   0.504695    1.99652
 -0.442616    0.411342   1.24364    -0.0993078  -0.670033
  0.394759    0.672543   1.53753    -0.569871   -1.55445
 -0.714782    0.910227   0.0708801  -0.317265    0.563627
 -0.0564112  -1.16758   -1.93046     2.44655     1.28567

julia> B = view(A, 1:3, 1:5);

julia> strides(B)
(1, 5)

julia> (x, y) = (3, 3)
(3, 3)

julia> B[x, y]
1.537533054119266

julia> B[1 + (x-1) * stride(B, 1) + (y-1) * stride(B, 2)]
1.9965179361688334

julia> parent(B)[1 + (x-1) * stride(B, 1) + (y-1) * stride(B, 2)]
1.537533054119266

The only reason I can get the correct element B[x,y] via linear indexing using the strides is because I know that B has an underlying DenseArray parent, for which the "memory index" computed via the strides is equivalent to linear indexing as provided by Julia. If I wouldn't know and use this, I would only have B to index into, and I could not use the "memory index" computed via the strides. Hence, even though B is a strided array, Julia does not provide an established high-level interface to use those strides to do indexing. One could do this using pointers and unsafe_ operations, but ultimately that amounts to assuming that there is a DenseArray underlying the given input array.

[akriegman]

I don't fully understand the considerations here, but I was imagining we could define linear indexing into sparse arrays to agree with their multiindexing behavior. So we could define B[idx] by

(s1, s2, ...) := size(B)
i1 := idx % s1
i2 := idx ÷ s1 % s2
i3 := idx ÷ (s1*s2) % s3
...
B[idx] := B[i1, i2, ...]

instead of just treating idx as a memory offset.

[Jutho]

Yes, that is what Julia's linear indexing already allows you to do. The point is that Strided.jl is specifically about about data that is layed out according to a strided pattern in memory, and that it will optimize the data access pattern utilizing that information.