New functions try_strides, is_ptr_loadable, and is_ptr_storable

[nhz2]

This is an alternative to #60894 and #61807

It adds three new functions to the strided array interface.

try_strides, is_ptr_loadable, and is_ptr_storable.

try_strides is similar to strides except it returns nothing when the array isn't strided. is_ptr_loadable and is_ptr_storable can be used to check if a pointer to an array element can be used to get or set the value.

XRef: #54715 #59435 #10889

[nhz2]

I've moved the test reorganization to #61250

[nhz2]

Question: is it a breaking change to replace strides with try_strides? I think that was the decision made in #30432 (comment)

[nhz2]

JuliaLang/LinearAlgebra.jl#1619 is the companion PR in LinearAlgebra.jl

[adienes]

little confused how this interacts with

a pointer to any element of the array can be obtained

so we might have an array where we can obtain a pointer to an element, but we can neither load nor store from that pointer?

[nhz2]

Yes, for example, due to padding alignment it is possible to have a write only reinterpret array. If this were wrapped around a readonly array the result would be a both unwriteable and unreadable array.

[gbaraldi]

This violates strict aliasing so is not legal

[nhz2]

By "wrapped" I don't mean unsafe_wrap I mean using a new wrapper array type, this avoids any TBAA issues.
Here is a more explicit example.

# A read only vector wrapper type
struct ReadOnlyWrapper{T, A<:AbstractVector{T}} <: AbstractVector{T}
    parent::A
end
Base.parent(A::ReadOnlyWrapper) = A.parent
Base.size(A::ReadOnlyWrapper) = size(A.parent)
Base.axes(A::ReadOnlyWrapper) = axes(A.parent)
Base.getindex(A::ReadOnlyWrapper, i::Int) = A.parent[i]
Base.cconvert(::Type{Ptr{T}}, A::ReadOnlyWrapper{T}) where {T} = Base.cconvert(Ptr{T}, A.parent)
Base.strides(A::ReadOnlyWrapper) = strides(A.parent)
Base.elsize(::Type{ReadOnlyWrapper{T,A}}) where {T,A} = Base.elsize(A)
Base.try_strides(A::ReadOnlyWrapper) = try_strides(A.parent)
Base.is_ptr_loadable(A::ReadOnlyWrapper) = is_ptr_loadable(A.parent)
# An array with padding
a = [(0x01, 0x0001)]
@show is_ptr_loadable(a)
@show is_ptr_storable(a)
b = reinterpret(UInt8, a);
@show is_ptr_loadable(b)
@show is_ptr_storable(b)
c = ReadOnlyWrapper(b)
@show is_ptr_loadable(c)
@show is_ptr_storable(c)

This results in:

is_ptr_loadable(a) = true
is_ptr_storable(a) = true
is_ptr_loadable(b) = false
is_ptr_storable(b) = true
is_ptr_loadable(c) = false
is_ptr_storable(c) = false

The ordering of the wrapping can be changed:

# An array with padding
a = [(0x01, 0x0001)]
@show is_ptr_loadable(a)
@show is_ptr_storable(a)
b = ReadOnlyWrapper(a)
@show is_ptr_loadable(b)
@show is_ptr_storable(b)
c = reinterpret(UInt8, b);
@show is_ptr_loadable(c)
@show is_ptr_storable(c)

Results in:

is_ptr_loadable(a) = true
is_ptr_storable(a) = true
is_ptr_loadable(b) = true
is_ptr_storable(b) = false
is_ptr_loadable(c) = false
is_ptr_storable(c) = false

[adienes]

maybe the example was just intended for the flavor and not as a full motivator, but even here I think some cracks are showing. namely that the is_ptr_*able definitions here are making assumptions about the padding / layout of Tuple{UInt8, UInt16}. but if Julia were to, say, choose to start representing that type without padding then the traits would become wrong.

there's a reason that array_subpadding in reinterpretarray.jl is doing a bunch of nontrivial work, because to query this correctly I think we have to be careful about reading the runtime values of datatype_alignment and sizeof

[nhz2]

The is_ptr_*able methods match the checks in getindex and setindex!, and this is tested. Please compare is_ptr_*able methods to getindex and setindex! methods for ReinterpretArray.

julia/base/reinterpretarray.jl

Lines 416 to 420 in c07cf8f

	@propagate_inbounds function getindex(a::ReinterpretArray{T,N,S}, inds::Vararg{Int, N}) where {T,N,S}
	check_readable(a)
	check_ptr_indexable(a) && return _getindex_ptr(a, inds...)
	_getindex_ra(a, inds[1], tail(inds))
	end

This calls check_readable(a) which is:

julia/base/reinterpretarray.jl

Lines 226 to 231 in c07cf8f

	function check_readable(a::ReinterpretArray{T, N, S} where N) where {T,S}
	# See comment in check_writable
	if !a.readable && !array_subpadding(T, S)
	throw(PaddingError(T, S))
	end
	end

Julia cannot arbitrarily change the padding of isbits because they are required to be C compatible.

[adienes]

do they match the checks though? because check_readable and check_writable call array_subpadding which is doing a bunch of work to read the type layout. or put another way:

function is_ptr_loadable(a::ReinterpretArray{T,N,S} where N) where {T,S}
    is_ptr_loadable(parent(a)) && (a.readable || array_subpadding(T, S))
end

function is_ptr_storable(a::ReinterpretArray{T,N,S} where N) where {T,S}
    is_ptr_storable(parent(a)) && (a.writable || array_subpadding(S, T))
end

why isn't is_ptr_*able(a) = is_ptr_*able(parent(a)) sufficient?

[nhz2]

Only checking the parent isn't enough to ensure there are no padding issues caused by the reinterpret. For that reason (a.readable || array_subpadding(T, S)) is also checked.

The tests in test/reinterpretarray.jl lines 392 and 395 fail with is_ptr_*able(a) = is_ptr_*able(parent(a)).

[nhz2]

The practical point is that I want zip_crc32(reinterpret_array_with_strange_padding_issues) to throw some error instead of silently returning something undefined.

[adienes]

is this Int necessary?

[nhz2]

Yes, try_strides should return Union{Dims, Nothing}, but step(first(I)) might not be an Int, for example:

julia> a = view([1,2,3], Int128(1):Int128(2):Int128(3))
2-element view(::Vector{Int64}, 1:2:3) with eltype Int64:
 1
 3

julia> step(first(a.indices)) |> typeof
Int128

[adienes]

maybe want _checkcontiguous fast path like strides has? same for ReshapedArray

[nhz2]

Yes, this could be done as a future optimization.

[adienes]

Suggested change

	try_strides(A::AbstractArray) = nothing
	function try_strides(A::AbstractArray)
	try
	strides(A)
	catch e
	(e isa MethodError) ? return nothing : rethrow()
	end
	end

should we "try" harder ? also this means the invariant you claim in the docstring

Return a tuple of the memory strides in each dimension if `A` is strided.
Otherwise return `nothing`.

will be more true. otherwise, any existing strided array type that doesn't yet define try_strides will violate this.

[nhz2]

I think it is better to go through and add these methods to existing strided array types. This is good chance to check they are actually following this new strided array interface. For example, JuliaLang/LinearAlgebra.jl#1618 #53996

[nhz2]

JuliaIO/InputBuffers.jl#16 is an example use case

[nhz2]

JuliaPy/PythonCall.jl#777
JuliaArrays/FixedSizeArrays.jl#191
and JuliaIO/ZipArchives.jl#106
are more examples. Here is a basic benchmark of the new version of zip_crc32 using try_strides.

using ZipArchives, FixedSizeArrays, PythonCall, BenchmarkTools, Random
data = rand(Xoshiro(1234), UInt8, 300000000);
f_data = FixedSizeArray(copy(data));
py_data = PyArray(copy(data));
@btime zip_crc32($data)
@btime zip_crc32($f_data)
@btime zip_crc32($py_data)

Julia 1.12.6:

70.287 ms (0 allocations: 0 bytes)
82.036 ms (3 allocations: 24.07 KiB)
147.411 ms (3 allocations: 24.07 KiB)

With try_strides:

70.127 ms (0 allocations: 0 bytes)
71.035 ms (0 allocations: 0 bytes)
71.421 ms (0 allocations: 0 bytes)

[gbaraldi]

Triage thinks this should behave more like a trait.

[nhz2]

Triage thinks this should behave more like a trait.

What would that look like?

[adienes]

(relaying what I remember from the discussion, not my own thoughts)

the argument in was that this function is being proposed for use to switch to strided algorithm. but this is what traits do, and the try_ family of functions is for cases where the operation might be unsuccessful for reasons completely independent of the algorithm design (e.g. data is malformed, or already locked somewhere, or some other exogenous factor)

so this should be closer to something that looks like isstrided or HasStride and then downstream users that want to implement fast paths for strided arrays dispatch on that.

furthermore, is_ptr_loadable and is_ptr_storable "should just be pointer" and either that works or it gives MethodError. although there is a pointer(::AbstractArray) fallback, that method is not generically correct and probably should be formally deprecated.

(back to my personal thought) my interpretation of the discussion was that the desired incarnation of this feature would look more like something that can close the referenced issues, not just update them.

XRef: #54715 #59435 #10889

[nhz2]

So basically implementing something like #10889 (comment)

abstract MemLayout
abstract Strided <: MemLayout
abstract UnitStrided <: Strided
abstract Contiguous <: UnitStrided

And then a function mapping from AbstractArray type to MemLayout type?

That sounds useful but still not complete. If the memory layout cannot be computed based on the type, I want to be able to get the strides if they exist at runtime. Also once the memory layout is known I want to know if loads and stores to a pointer to an element are equivalent to using setindex! or getindex! (this is two separate questions one for load and one for store).

Maybe these requirements are not needed by Base? A parallel strided array interface could live in a separate package with extensions for various array packages I am using.

[adienes]

something like that, although I'm not sure what the difference between UnitStrided and Contiguous would be. you are right that afaik, Base mostly doesn't need this, although it might be useful to replace some big unions.

I want to be able to get the strides if they exist at runtime

at the very least, the type should be able to give you enough information to know if a strides method will work or not, hopefully? then the caller can e.g. check strides(A) == (1,) as the example in ZipArchives does

it might be useful (at least for me, who remains a bit confused) to come up with a table of types/values that would occupy different spots in that proposed type hierarchy or different results from is_ptr_loadable and is_ptr_storable. I have to admit --- thinking about the proposal again --- it seems surprising that we'd want to recommend to users working with pointers to array locations directly given how easy it is to do incorrectly. even the proposed usage in this PR does not seem fully legal (#60964 (comment))

if getindex or setindex! are slow, maybe there is some simpler solution like making sure the user's methods are properly @propagate_inbounds ? Base.ReshapedArray and Base.ReinterpretArray are implemented very carefully to avoid UB; maybe users can have their array types wrap those directly rather than trying to fork / recreate them ?

[nhz2]

Yes using pointers is not recommended, but they are needed when calling C and Julia functions that operate on pointers.

I'm also not sure what UnitStrided means exactly.

I'm assuming Contiguous means that if the element type is isbits, the memory layout matches the corresponding Array.

Here are some examples of arrays that do not have Contiguous types but that are contiguous:

transpose(zeros(1,5))
view(zeros(5), 2:1:4)

I think these types of arrays are easy to run into.

Here are examples of arrays that do not have Strided types but that are strided:

reshape(view(zeros(4,4), 1:4, 1:2), 8)
reinterpret(UInt16, view(zeros(UInt8, 9, 8), 1:8, 1:2:8))

I'm not sure if these exist in the wild, but the existing strides methods handle these cases.

codeunits("abc") is the main example in Base where is_ptr_storable is false and is_ptr_loadable is true.

In general setindex! and getindex may include various side effects or checks that need to be run, and cannot be safely bypassed by directly using the pointer to an element. is_ptr_storable and is_ptr_loadable act as a way for new types to opt in to allowing setindex! and getindex to be bypassed.

[nhz2]

Looking at @stevengj 's function in Blosc.jl

maybe Contiguous shouldn't have to exactly match Array's memory layout?