juliaarrays / arrayinterface.jl

Designs for new Base array interface primitives, used widely through scientific machine learning (SciML) and other organizations

ArrayInterface.jl

The AbstractArray interface in Base Julia is still relatively young. The purpose of this library is to solidify extensions to the current AbstractArray interface, which are put to use in package ecosystems like DifferentialEquations.jl. Since these libraries are live, this package will serve as a staging ground for ideas before they are merged into Base Julia. For this reason, no functionality is exported so that if such functions are added and exported in a future Base Julia, there will be no issues with the upgrade.

List of things to add

• JuliaLang/julia#22216
• JuliaLang/julia#22218
• JuliaLang/julia#22622
• JuliaLang/julia#25760
• JuliaLang/julia#25107

Array Types to Handle

The following common array types are being understood and tested as part of this development.

• Array
• Various versions of sparse arrays
• SArray
• MArray
• FieldVector
• ArrayPartition
• VectorOfArray
• DistributedArrays
• GPUArrays (CLArrays and CuArrays)
• AFArrays
• MultiScaleArrays
• LabelledArrays

StaticArrayInterface.jl

If one is looking for an interface which includes functionality for statically-computed values, see StaticArrayInterface.jl. This was separated from ArrayInterface.jl because it includes a lot of functionality that does not give substantive improvements to the interface, and is likely to be deprecated in the near future as the compiler matures to automate a lot of its optimizations.

Breaking Release Notes

7.0: Setup to use the new v1.9 package extensions. All static interface functions were moved to StaticArrayInterface.jl. All packages using one of the subpackages (ArrayInterfaceCore, ArrayInterfaceStaticArrays, etc.) should update to v7 by simply depending on ArrayInterface.jl, unless static interface functions were used in which case a dependency on StaticArrayInterface.jl is also necessary.

6.0: ArrayInterface.jl completely removed all usage of Requires.jl and conditional dependencies due to compile time impact. All of the Requires.jl support changed to subpackages within the repository which are registered in the General registry. These subpackages are required by any packages which seek to use the additional functionality.

2.0: Changed the default of ismutable(array::AbstractArray) = true. We previously defaulted to Base.@pure ismutable(array::AbstractArray) = typeof(array).mutable, but there are a lot of cases where this tends to not work out in a way one would expect. For example, if you put a normal array into an immutable struct that adds more information to it, this is considered immutable, even if all of the setindex! methods work (by forwarding to the mutable array). Thus, it seems safer to just always assume mutability is standard for an array, and allow arrays to opt-out.

arrayinterface.jl's Issues

Minor typo

In the repo description, "interface primates" should probably be "interface primatives" (unless the plan is to hand the task over to some teams of lemurs) 😛

Bug in `device` for reshaped views

Reported by @chriselrod in #156 (comment):

 julia> using ArrayInterface
 
 julia> u_base = randn(10, 10); u_view = view(u_base, 2:3, :); u_reshaped_view = reshape(u_view, 1, length(u_view));
 
 julia> ArrayInterface.device(u_reshaped_view)
 ArrayInterface.CPUPointer()

This u_reshaped_view is not strided, so it should not be a CPUPointer(), but a CPUIndex() (which currently is not supported by LoopVectorization).

Add an `is_resizable` trait?

It'd be nice if there was an abstract array interface specifying whether an array is resizable in-place.

dimensional traits/predicates?

It would be useful to dispatch_on/query isvector ismatrix (without directly examining the type to find N, which may not exist).

requires.jl static arrays warning

┌ Warning: Error requiring `StaticArrays` from `ArrayInterface`
--
  | │   exception =
  | │    too many parameters for type
  | │    Stacktrace:
  | │      [1] top-level scope
  | │        @ ~/.cache/julia-buildkite-plugin/depots/609dd3f6-6d58-4fb2-a8ce-f94307987da5/packages/ArrayInterface/GV582/src/ArrayInterface.jl:930
  | │      [2] eval
  | │        @ ./boot.jl:360 [inlined]
  | │      [3] eval
  | │        @ ~/.cache/julia-buildkite-plugin/depots/609dd3f6-6d58-4fb2-a8ce-f94307987da5/packages/ArrayInterface/GV582/src/ArrayInterface.jl:1 [inlined]
  | │      [4] (::ArrayInterface.var"#53#80")()
  | │        @ ArrayInterface ~/.cache/julia-buildkite-plugin/depots/609dd3f6-6d58-4fb2-a8ce-f94307987da5/packages/Requires/7Ncym/src/require.jl:99
  | │      [5] err(f::Any, listener::Module, modname::String)
  | │        @ Requires ~/.cache/julia-buildkite-plugin/depots/609dd3f6-6d58-4fb2-a8ce-f94307987da5/packages/Requires/7Ncym/src/require.jl:47
  | │      [6] (::ArrayInterface.var"#52#79")()
  | │        @ ArrayInterface ~/.cache/julia-buildkite-plugin/depots/609dd3f6-6d58-4fb2-a8ce-f94307987da5/packages/Requires/7Ncym/src/require.jl:98
  | │      [7] withpath(f::Any, path::String)
  | │        @ Requires ~/.cache/julia-buildkite-plugin/depots/609dd3f6-6d58-4fb2-a8ce-f94307987da5/packages/Requires/7Ncym/src/require.jl:37
  | │      [8] (::ArrayInterface.var"#51#78")()
  | │        @ ArrayInterface ~/.cache/julia-buildkite-plugin/depots/609dd3f6-6d58-4fb2-a8ce-f94307987da5/packages/Requires/7Ncym/src/require.jl:97
  | │      [9] listenpkg(f::Any, pkg::Base.PkgId)
  | │        @ Requires ~/.cache/julia-buildkite-plugin/depots/609dd3f6-6d58-4fb2-a8ce-f94307987da5/packages/Requires/7Ncym/src/require.jl:20
  | │     [10] macro expansion
  | │        @ ~/.cache/julia-buildkite-plugin/depots/609dd3f6-6d58-4fb2-a8ce-f94307987da5/packages/Requires/7Ncym/src/require.jl:95 [inlined]
  | │     [11] __init__()
  | │        @ ArrayInterface ~/.cache/julia-buildkite-plugin/depots/609dd3f6-6d58-4fb2-a8ce-f94307987da5/packages/ArrayInterface/GV582/src/ArrayInterface.jl:869
  | │     [12] _include_from_serialized(path::String, depmods::Vector{Any})
  | │        @ Base ./loading.jl:674
  | │     [13] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String)
  | │        @ Base ./loading.jl:760
  | │     [14] _tryrequire_from_serialized(modkey::Base.PkgId, build_id::UInt64, modpath::String)
  | │        @ Base ./loading.jl:689
  | │     [15] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String)
  | │        @ Base ./loading.jl:749
  | │     [16] _tryrequire_from_serialized(modkey::Base.PkgId, build_id::UInt64, modpath::String)
  | │        @ Base ./loading.jl:689
  | │     [17] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String)
  | │        @ Base ./loading.jl:749
  | │     [18] _tryrequire_from_serialized(modkey::Base.PkgId, build_id::UInt64, modpath::String)
  | │        @ Base ./loading.jl:689
  | │     [19] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String)
  | │        @ Base ./loading.jl:749
  | │     [20] _tryrequire_from_serialized(modkey::Base.PkgId, build_id::UInt64, modpath::String)
  | │        @ Base ./loading.jl:689
  | │     [21] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String)
  | │        @ Base ./loading.jl:749
  | │     [22] _require(pkg::Base.PkgId)
  | │        @ Base ./loading.jl:998
  | │     [23] require(uuidkey::Base.PkgId)
  | │        @ Base ./loading.jl:914
  | │     [24] require(into::Module, mod::Symbol)
  | │        @ Base ./loading.jl:901
  | │     [25] include(mod::Module, _path::String)
  | │        @ Base ./Base.jl:386
  | │     [26] exec_options(opts::Base.JLOptions)
  | │        @ Base ./client.jl:285
  | │     [27] _start()
  | │        @ Base ./client.jl:485
  | └ @ Requires ~/.cache/julia-buildkite-plugin/depots/609dd3f6-6d58-4fb2-a8ce-f94307987da5/packages/Requires/7Ncym/src/require.jl:49

Chris R told me to ping you @chriselrod

`size` fails for view with discontiguous slices

For array views, size only seems to work when the slices are contiguous (of the form : or a:b). When there is a step (as in 2:2:10), the associated dimension is lost, which I guess is not intended.

Example:

u = rand(4, 6)
v = @view u[1:3, :]
w = @view u[2:2:4, :]

ArrayInterface.size(u)  # returns (4, 6) as expected
ArrayInterface.size(v)  # returns (3, 6) as expected
ArrayInterface.size(w)  # returns (6, ) -- the first dimension is lost

The issue seems to come from the _size function in stridelayout.jl, where possible range types are limited to Base.Slice and AbstractUnitRange:

@generated function _size(A::Tuple{Vararg{Any,N}}, inds::I, l::L) where {N, I<:Tuple, L}
    t = Expr(:tuple)
    for n in 1:N
        if (I.parameters[n] <: Base.Slice)
            push!(t.args, :(@inbounds(_try_static(A[$n], l[$n]))))
        elseif I.parameters[n] <: AbstractUnitRange
            push!(t.args, Expr(:ref, :l, n))
        end
    end
    Expr(:block, Expr(:meta, :inline), t)
end

One possible fix is to replace AbstractUnitRange by the more general OrdinalRange. Is this an acceptable solution?

Should setindex(x::AbstractArray, v, i) use mutate-or-widen strategy?

Current implementation of setindex(x::AbstractArray, v, i) does not support widening of the element type:

https://github.com/JuliaDiffEq/ArrayInterface.jl/blob/c8513cbc0b9dc24e6674206255838abec553b854/src/ArrayInterface.jl#L22-L26

I think it would be useful to support eltype-widening so that setindex([0], 1.2, 1) works. I implemented this in BangBang.jl internally (BangBang.jl provides a uniform API for mutable and immutable containers). It is useful to have type-changing setter because you can use it to, e.g., implement map based on mutate-or-widen strategy.

Actual code is pretty simple (from https://github.com/tkf/BangBang.jl/blob/94eab8d728194fac0fa41c840ca398600a50f903/src/NoBang/base.jl#L99-L107):

function _setindex(xs::AbstractArray, v, I...)
    T = promote_type(eltype(xs), typeof(v))
    ys = similar(xs, T)
    if eltype(xs) !== Union{}
        copy!(ys, xs)
    end
    ys[I...] = v
    return ys
end

One disadvantage of this approach is that it fails when there is an undef element. That's why Union{} was special-cased above. To support more general cases, I think we need something like copy_if_defined! (or maybe Base.copy! should do it automatically?).

Trouble using `indices` output for indexing.

Turns out similar doesn't like OptionallyStaticUnitRange.

julia> ones(10,10)[ArrayInterface.indices(1:10)]
ERROR: MethodError: no method matching similar(::Type{Array{Float64,1}}, ::Tuple{Base.IdentityUnitRange{ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1},Int64}}})
Closest candidates are:
  similar(::AbstractArray{T,N} where N, ::Tuple) where T at abstractarray.jl:630
  similar(::Type{T}, ::Union{Integer, AbstractUnitRange}...) where T<:AbstractArray at abstractarray.jl:673
  similar(::Type{T}, ::Tuple{Vararg{Int64,N}} where N) where T<:AbstractArray at abstractarray.jl:675
  ...
Stacktrace:
 [1] _array_for(::Type{Float64}, ::Base.Slice{ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1},Int64}}, ::Base.HasShape{1}) at ./array.jl:678
 [2] getindex(::Array{Float64,2}, ::Base.Slice{ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1},Int64}}) at ./array.jl:834
 [3] top-level scope at REPL[4]:1

Relying on constant prop for inference is unreliable

Most of the @inferred unit tests are on the edge of infer-ability. Add just a little more complexity -- such as by calling them from real code -- and inference fails.

│   │┌ @ /home/chriselrod/.julia/dev/ArrayInterface/src/stridelayout.jl:339 within `offsets'
│   ││┌ @ /home/chriselrod/.julia/dev/ArrayInterface/src/stridelayout.jl:342 within `macro expansion'
│   │││┌ @ /home/chriselrod/.julia/dev/ArrayInterface/src/stridelayout.jl:333 within `offsets'
│   ││││┌ @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:419 within `indices'
│   │││││┌ @ /home/chriselrod/.julia/dev/ArrayInterface/src/dimensions.jl:164 within `axes' @ abstractarray.jl:70
│   ││││││┌ @ int.jl:442 within `<='
│   │││││││ %3   = Base.sle_int(1, 2)::Bool
│   ││││││└
└───││││││        goto #3 if not %3
2 ──││││││        goto #4
3 ──││││││        goto #4
    ││││││ @ /home/chriselrod/.julia/dev/ArrayInterface/src/dimensions.jl:164 within `axes'
4 ┄─││││││ %7   = φ (#3 => 1)::Int64
│   ││││││ %8   = φ (#2 => $(QuoteNode(SOneTo(16))), #3 => $(QuoteNode(Base.OneTo(1))))::Union{SOneTo{16}, Base.OneTo{Int64}}
└───││││││        goto #5
    │││││└
5 ──│││││ %10  = (isa)(%8, SOneTo{16})::Bool
└───│││││        goto #7 if not %10
6 ──│││││        goto #10
7 ──│││││ %13  = (isa)(%8, Base.OneTo{Int64})::Bool
└───│││││        goto #9 if not %13
    │││││ @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:419 within `indices' @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:406
    │││││┌ @ abstractarray.jl:256 within `eachindex'
    ││││││┌ @ abstractarray.jl:109 within `axes1'
    │││││││┌ @ abstractarray.jl:89 within `axes'
    ││││││││┌ @ tuple.jl:213 within `map'
    │││││││││┌ @ range.jl:330 within `OneTo' @ range.jl:321
    ││││││││││┌ @ promotion.jl:421 within `max'
    │││││││││││┌ @ int.jl:83 within `<'
8 ──││││││││││││ %15  = Base.slt_int(%7, 0)::Bool
│   │││││││││││└
│   │││││││││││ %16  = Base.ifelse(%15, 0, %7)::Int64
│   │││││└└└└└└
│   │││││ @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:419 within `indices' @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:408
│   │││││┌ @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:94 within `OptionallyStaticUnitRange' @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:75
│   ││││││ %17  = %new(ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, Int64}, ArrayInterface.StaticInt{1}(), %16)::ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, Int64}
│   │││││└
│   │││││┌ @ indices.jl:348 within `Slice' @ indices.jl:348
│   ││││││ %18  = %new(Base.Slice{ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, Int64}}, %17)::Base.Slice{ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, Int64}}
│   │││││└
│   │││││ @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:419 within `indices'
└───│││││        goto #10
9 ──│││││        Core.throw(ErrorException("fatal error in type inference (type bound)"))
└───│││││        unreachable
10 ┄│││││ %22  = φ (#6 => $(QuoteNode(Base.Slice(ArrayInterface.StaticInt{1}():ArrayInterface.StaticInt{1}():ArrayInterface.StaticInt{16}()))), #8 => %18)::Union{Base.Slice{ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{16}}}, Base.Slice{ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, Int64}}}
└───│││││        goto #11
    │││└└
    │││┌ @ /home/chriselrod/.julia/dev/ArrayInterface/src/stridelayout.jl:334 within `offsets'
11 ─││││ %24  = (isa)(%22, Base.Slice{ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{16}}})::Bool
└───││││        goto #13 if not %24
12 ─││││        goto #16
13 ─││││ %27  = (isa)(%22, Base.Slice{ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, Int64}})::Bool
└───││││        goto #15 if not %27
14 ─││││        goto #16
15 ─││││        Core.throw(ErrorException("fatal error in type inference (type bound)"))
└───││││        unreachable
    │││└
    │││┌ @ /home/chriselrod/.julia/dev/ArrayInterface/src/stridelayout.jl:335 within `offsets'
16 ┄││││        goto #17
    │││└
    │││┌ @ /home/chriselrod/.julia/dev/ArrayInterface/src/stridelayout.jl:333 within `offsets'
    ││││┌ @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:419 within `indices'
    │││││┌ @ /home/chriselrod/.julia/dev/ArrayInterface/src/dimensions.jl:164 within `axes' @ abstractarray.jl:70
    ││││││┌ @ int.jl:442 within `<='
17 ─│││││││ %33  = Base.sle_int(2, 2)::Bool
│   ││││││└
└───││││││        goto #19 if not %33
18 ─││││││        goto #20
19 ─││││││        goto #20
    ││││││ @ /home/chriselrod/.julia/dev/ArrayInterface/src/dimensions.jl:164 within `axes'
20 ┄││││││ %37  = φ (#19 => 1)::Int64
│   ││││││ %38  = φ (#18 => $(QuoteNode(SOneTo(16))), #19 => $(QuoteNode(Base.OneTo(1))))::Union{SOneTo{16}, Base.OneTo{Int64}}
└───││││││        goto #21
    │││││└
21 ─│││││ %40  = (isa)(%38, SOneTo{16})::Bool
└───│││││        goto #23 if not %40
22 ─│││││        goto #26
23 ─│││││ %43  = (isa)(%38, Base.OneTo{Int64})::Bool
└───│││││        goto #25 if not %43
    │││││ @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:419 within `indices' @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:406
    │││││┌ @ abstractarray.jl:256 within `eachindex'
    ││││││┌ @ abstractarray.jl:109 within `axes1'
    │││││││┌ @ abstractarray.jl:89 within `axes'
    ││││││││┌ @ tuple.jl:213 within `map'
    │││││││││┌ @ range.jl:330 within `OneTo' @ range.jl:321
    ││││││││││┌ @ promotion.jl:421 within `max'
    │││││││││││┌ @ int.jl:83 within `<'
24 ─││││││││││││ %45  = Base.slt_int(%37, 0)::Bool
│   │││││││││││└
│   │││││││││││ %46  = Base.ifelse(%45, 0, %37)::Int64
│   │││││└└└└└└
│   │││││ @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:419 within `indices' @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:408
│   │││││┌ @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:94 within `OptionallyStaticUnitRange' @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:75
│   ││││││ %47  = %new(ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, Int64}, ArrayInterface.StaticInt{1}(), %46)::ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, Int64}
│   │││││└
│   │││││┌ @ indices.jl:348 within `Slice' @ indices.jl:348
│   ││││││ %48  = %new(Base.Slice{ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, Int64}}, %47)::Base.Slice{ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, Int64}}
│   │││││└
│   │││││ @ /home/chriselrod/.julia/dev/ArrayInterface/src/ranges.jl:419 within `indices'
└───│││││        goto #26
25 ─│││││        Core.throw(ErrorException("fatal error in type inference (type bound)"))
└───│││││        unreachable
26 ┄│││││ %52  = φ (#22 => $(QuoteNode(Base.Slice(ArrayInterface.StaticInt{1}():ArrayInterface.StaticInt{1}():ArrayInterface.StaticInt{16}()))), #24 => %48)::Union{Base.Slice{ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, ArrayInterface.StaticInt{16}}}, Base.Slice{ArrayInterface.OptionallyStaticUnitRange{ArrayInterface.StaticInt{1}, Int64}}}
└───│││││        goto #27
    │││└└

E.g., note the 4 failures above, along with Core.throw(ErrorException("fatal error in type inference (type bound)")).

TODO: Finish setting up Buildkite

Make accessing and mutating `parent` part of the API

It'd be nice if there was a common API for accessing and mutating parent data for wrapper array types.

Should we have pop, push and friends?

StaticArrays has all of these non-mutating versions of typically mutating methods implemented, which might be a good fit for this package since we have setindex. It would be nice if these were available outside of a dependency on StaticArrays and part of the general array interface.

Add is_strided

see JuliaLang/julia#2345 (comment)

use StructArrays.jl

Your currently just have

aos_to_soa(x) = x

I feel like it might be worth it to just add a dependancy on StructArrays.jl, which imo should really be a baked in part of Julia's array interface.

Register

@JuliaRegistrator register()

Safe co-iteration along axes

This came up with the following comment in from #61 .

Users will often check for equal sizes in front of a loop, and eachindex will guarantee this as well

for I in eachindex(A,C) # checks for equivalent size
   
end

Because eachindex isn't always appropriate (e.g., when you don't have exact correspondence between axes of the arrays), I think it'd be a good to have something like an eachindexalongaxis:

function mymatmulkernel!(C, A, B)
    @avx for m in eachindexalongaxis((C,A),(1,1)), n in eachindexalongaxis((C,B),(2,2)), k in eachindexalongaxis((A,B),(2,1))
        C[m,n] += A[m,k] * B[k,n]
    end
end

where it checks sizes. eachindexalongaxis or each_index_along_axis doesn't really roll off > the tongue, and I don't know where such a function should live.
I'd also have to think a little more about how to use it and actually pass the information of > dimension equality (to be combined with density guarantees) to deduce stride-equality. The hacky way would be to just pattern match eachindex and eachindexalongindex in the macro, handling dimension equality through the expression, while handling density info based on types.

I'd like to propose the use of indices for specifically retrieving the iterator that corresponds to each axis of an array. It would function like axes but if there's any sort of other data associated with an axis it throws it away so there's a clear relationship between memory layout and and each dimension.

In practice it would look like this:

julia> indices(ones(2,3))
(Base.OneTo(2), Base.OneTo(3))

julia> indices(ones(2,3), 1)
Base.OneTo(2)

julia> indices(ones(2,3), 2)
Base.OneTo(3)

julia> indices(ones(2,3), ones(2,3), 2)
Base.OneTo(3)

julia> indices(ones(2,3), ones(2,2), 2)
ERROR: DimensionMismatch("all inputs to indices must have the same indices, got Base.OneTo(3) and Base.OneTo(2)")

In this case the first several examples are the same as if using axes, but it helps for cases like those in AxisIndices.jl where each axis can have additional information (e.g., keys) associated with the indices. The last two examples would ensure that the iterator for each collection is the same along the second dimension.

WARNING: Wrapping Vararg directly in UnionAll is deprecated (wrap the tuple instead)

See also: JuliaSIMD/LoopVectorization.jl#175

cc: @chriselrod

interface for rewrapping Array wrappers

Some wrapper types like AxisArray, NamedDImsArray and DimArray (and probably LabelledArray?) need to be the "outside" wrapper for, e.g. keeping the axis index accurate and dispatch for getindex methods.

This issue just came up on discourse:
https://discourse.julialang.org/t/multiple-inheritance-carry-over-members/58355/8

People kind of expect a DimArray to keep working after wrapping with other types. It could be useful to have a trait and a shared method here that allows rewrapping:

rewrap(f, A) and rewrap(f, A, I...) could be used in other array constructors to allow rewrapping where necessary.

isouterwrapper could indicate if an object is an outer wrapper type.

Collaboration, transfer, and adding iterators

JuliaLang/julia#31563 (comment), posted by @AriMKatz, drew my attention to this package. I agree that there appear to be a subset of shared goals between this package and https://github.com/timholy/ArrayIteration.jl. This is a proposal to collaborate so we don't duplicate effort.

To summarize ArrayIteration (because the README is ancient and talks about some stuff that long ago made it into Base): if you had to boil it down to a single goal, it's to provide an efficient fallback for A*B for matrices A and B with arbitrary sparsity patterns and representation. The general problem is that if you have to specialize on both typeof(A) and typeof(B), you have to write O(N^2) methods where N is the number of array types. (And wrappers like Adjoint make the problem worse because it's essentially O(4N^2) or O(9N^2).) ArrayIteration takes the stance that you can solve the problem with O(N) methods if you instead create a fancy iterator for each matrix that allows synchronization with another iterator. You sync the row iterator of A to the column iterator of B and thereby mostly just visit entries that are nonzero in both.

ArrayIteration is basically waiting for JuliaLang/julia#34126 (so we don't pay the allocation penalty for wrappers) and probably JuliaLang/julia#34847 before I trying bringing it back to life. The state of the art is the teh/stored2 branch which was almost functional, before the performance problems due to wrapper creation put it on hold. (The README, however, was not updated to reflect the state of that branch.)

I'd propose we all get together in a package hosted at JuliaArrays, which seems like the natural home for generic array stuff. We could start by transferring this one, since I'm guessing it's working whereas ArrayIteration is stalled out. If anyone wants to just steal the code in ArrayIteration that's fine too, but I don't think keeping it under the JuliaDiffEq umbrella makes sense if it's supposed to be generic array handling. We could also keep them two packages, but in that case we'd want to make sure they're independent but interoperable.

Also CCing @vchuravy whose https://github.com/JuliaGPU/KernelAbstractions.jl was also mentioned. I am less certain there are true shared goals there, but just in case.

Register

@JuliaRegistrator register()

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Using `Base.@pure` in ArrayInterface causes 265 problems

JuliaLang/julia#265 (comment)

I think we should err on the side of finding any other way.

ismutable wrong for FillArrays

FillArrays are not mutable.

julia> ArrayInterface.ismutable(typeof(Fill(4, (2,2))))
true

I think that
https://github.com/SciML/ArrayInterface.jl/blob/730592e19c5e3effb086a5cc0b61ecc1a1936a76/src/ArrayInterface.jl#L75-L81

Should be

 function ismutable(::Type{T}) where {T<:AbstractArray} 
     if parent_type(T) <: T 
         return false  # no parent, and have not hit a type we know about 
     else 
         return ismutable(parent_type(T)) 
     end 
 end

ismutable{::Type{<:Array}} = true
ismutable{::Type{<:SparseVector}} = true
ismutable{::Type{<:SparseMatrixCSC}} = true
ismutable{::Type{<:BitArray}} = true

Alternatively
the fallback case could be: return T.mutable rather than false.
On the assumption that if you are a mutable struct that has subtypes AbstractArray you probably have defined setindex

using HybridArrays.jl

I'm trying to simulate interacting particle systems. The state vector has dimension '2N' or '3N'.
It seems beneficial to use an Array of StaticArrays for the 2D/3D points.

HybridArrays.jl would be super nice here, since it avoids a Vector of SVectors
and is still very fast. (timings are in the code below; notice that access is exactly as for a 2xN Array!).

Issues:

• HybridArrays does not work with all methods (for example not with Rosenbrock23(), Rosenbrock23(autodiff=false)).
• Missing plot recipe (should be easy to implement)

I can try to work on the integration if wanted. Of course, hints are very welcome :D

I think HybridArrays could be an alternative to an Array of StaticArrays and without problems like #47 and SciML/LabelledArrays.jl#68.

using HybridArrays, StaticArrays
using DifferentialEquations, BenchmarkTools, Random

N = 50
Dim = 2

# just some random dynamics to show speed benefit of static arrays
function f_ode(du, u, p, t)
    N = size(u,2)
    du[:] .= 0.
    for i = 1:N
        for j = 1:i-1
            v = u[:,i] - u[:,j]
            dist_squared = sum( x -> x^2, v)
            du[:,i] -= v * dist_squared
            du[:,j] += v * dist_squared
        end
    end
end

# same as f_ode but for vector of svectors
function f_ode_vec_svec(du, u, p, t)
    N = size(u,1)
    du[:] .-= du[:]  #(sorry, I don't know how to zero a vector of svectors fast.)
    for i = 1:N
        for j = 1:i-1
            v = u[i] - u[j]
            dist_squared = sum( x -> x^2, v)
            du[i] -= v * dist_squared
            du[j] += v * dist_squared
        end
    end
end


Random.seed!(1)
z0 = rand(Dim,N)
z0_hybrid = HybridArray{Tuple{Dim,StaticArrays.Dynamic()}}(z0)

Random.seed!(1)
z0_vec_svec = rand(SVector{Dim,Float64},N)


prob          = ODEProblem(f_ode,          z0,          (0., 10.))
prob_hybrid   = ODEProblem(f_ode,          z0_hybrid,   (0., 10.))
prob_vec_svec = ODEProblem(f_ode_vec_svec, z0_vec_svec, (0., 10.))

sol          = @btime solve(prob,          Euler(), dt=0.01) seconds=2
# 682.401 ms (33153208 allocations: 1.65 GiB)
sol_hybrid   = @btime solve(prob_hybrid,   Euler(), dt=0.01) seconds=2
# 7.337 ms (14070 allocations: 3.04 MiB)
sol_vec_svec = @btime solve(prob_vec_svec, Euler(), dt=0.01) seconds=2
# 5.019 ms (20074 allocations: 4.78 MiB)


sol2          = solve(prob,          Rosenbrock23())
# works
sol2_hybrid   = solve(prob_hybrid,   Rosenbrock23())
# MethodError: copyto!(::DiffEqBase.DiffEqBC{Array{Float64,2}}, ::Base.Broadcast.Broadcasted{HybridArrays.HybridArrayStyle{2},Tuple{Base.OneTo{Int64},Base.OneTo{Int64}},typeof(muladd),Tuple{Base.Broadcast.Broadcasted{HybridArrays.HybridArrayStyle{2},Nothing,typeof(DiffEqBase.ODE_DEFAULT_NORM),Tuple{HybridArray{Tuple{2,StaticArrays.Dynamic()},Float64,2,2,Array{Float64,2}},Float64}},Float64,Float64}})
# is ambiguous.
sol2_vec_svec = solve(prob_vec_svec, Rosenbrock23())
# ArgumentError: Cannot create a dual over scalar type SArray{Tuple{2},Float64,1,2}.
#  If the type behaves as a scalar, define FowardDiff.can_dual.

using Plots
plot( sol[1,1,:], sol[2,1,:])
plot( sol_hybrid[1,1,:], sol_hybrid[2,1,:])
plot( getindex.(sol_vec_svec[1,:],1) , getindex.(sol_vec_svec[1,:],2))
# plots look all the same

plot( sol )
# works
plot( sol_hybrid )
# MethodError: no method matching u_n(...)  
# Note: Here a ::CartesianIndex{2}  does not match the expected ::Int64
plot( sol_vec_svec )
# MethodError: no method matching one(::Type{SArray{Tuple{2},Float64,1,2}})

I can provide full error messages if needed.

SparseArrays?

Is the absence of SparseArrays intentional?

Invalidations

I checked this to see how ready this might be to include in ImageCore. Here are the invalidations from this package on Julia 1.6-rc1:

julia> trees
6-element Vector{SnoopCompile.MethodInvalidations}:
 inserting promote_rule(::Type{Bool}, ::Type{var"#s12"} where var"#s12"<:ArrayInterface.StaticBool) in ArrayInterface at /home/tim/.julia/dev/ArrayInterface/src/static.jl:251 invalidated:
   backedges: 1: superseding promote_rule(::Type{Bool}, ::Type{T}) where T<:Number in Base at bool.jl:4 with MethodInstance for promote_rule(::Type{Bool}, ::Type{var"#s444"} where var"#s444"<:Integer) (1 children)
   18 mt_cache

 inserting step(::ArrayInterface.OptionallyStaticUnitRange) in ArrayInterface at /home/tim/.julia/dev/ArrayInterface/src/ranges.jl:102 invalidated:
   backedges: 1: superseding step(r::AbstractUnitRange{T}) where T in Base at range.jl:544 with MethodInstance for step(::AbstractUnitRange{Int64}) (1 children)

 inserting eachindex(r::Union{var"#s12", var"#s9"} where {var"#s12"<:ArrayInterface.OptionallyStaticUnitRange, var"#s9"<:ArrayInterface.OptionallyStaticStepRange}) in ArrayInterface at /home/tim/.julia/dev/ArrayInterface/src/ranges.jl:401 invalidated:
   backedges: 1: superseding eachindex(A::AbstractVector{T} where T) in Base at abstractarray.jl:256 with MethodInstance for eachindex(::AbstractVector{T} where T) (1 children)

 inserting convert(::Type{T}, ::ArrayInterface.StaticInt{N}) where {T<:Number, N} in ArrayInterface at /home/tim/.julia/dev/ArrayInterface/src/static.jl:22 invalidated:
   backedges: 1: superseding convert(::Type{T}, x::Number) where T<:Number in Base at number.jl:7 with MethodInstance for convert(::Type{UInt64}, ::Integer) (2 children)

 inserting step(r::ArrayInterface.OptionallyStaticStepRange) in ArrayInterface at /home/tim/.julia/dev/ArrayInterface/src/ranges.jl:190 invalidated:
   mt_backedges: 1: signature Tuple{typeof(step), OrdinalRange{Int64, Int64}} triggered MethodInstance for (::Base.IteratorsMD.var"#21#23")(::OrdinalRange{Int64, Int64}) (3 children)

 inserting (::Type{T})(x::ArrayInterface.StaticInt{N}) where {T<:Integer, N} in ArrayInterface at /home/tim/.julia/dev/ArrayInterface/src/static.jl:26 invalidated:
   mt_backedges:  1: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for _array_for(::Type{Any}, ::Any, ::Base.HasLength) (0 children)
                  2: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for to_shape(::Integer) (0 children)
                  3: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for findnext(::Test.var"#3#5", ::AbstractString, ::Integer) (0 children)
                  4: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for findnext(::Test.var"#4#6", ::AbstractString, ::Integer) (0 children)
                  5: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for _array_for(::Type{T}, ::Any, ::Base.HasLength) where T (0 children)
                  6: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for _array_for(::Type{String}, ::Any, ::Base.HasLength) (0 children)
                  7: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for _array_for(::Type{Expr}, ::Any, ::Base.HasLength) (0 children)
                  8: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for _similar_for(::Vector{T} where T, ::Type, ::Base.Generator{_A, Test.var"#24#25"{Int64}} where _A, ::Base.HasLength) (0 children)
                  9: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for _similar_for(::Vector{T} where T, ::DataType, ::Base.Generator{_A, Test.var"#24#25"{Int64}} where _A, ::Base.HasLength) (0 children)
                 10: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for findnext(::Pkg.Types.var"#21#22"{String}, ::AbstractString, ::Integer) (0 children)
                 11: signature Tuple{Type{UInt64}, Integer} triggered MethodInstance for convert(::Type{UInt64}, ::Integer) (0 children)
                 12: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for var"#handle_message#2"(::Base.Iterators.Pairs{Symbol, _A, Tuple{Symbol}, _B} where {_A, _B}, ::typeof(Base.CoreLogging.handle_message), ::Base.CoreLogging.SimpleLogger, ::Base.CoreLogging.LogLevel, ::String, ::Nothing, ::Symbol, ::Symbol, ::Nothing, ::Int64) (1 children)
                 13: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for var"#handle_message#2"(::Base.Iterators.Pairs{Symbol, _A, Tuple{Symbol}, _B} where {_A, _B}, ::typeof(Base.CoreLogging.handle_message), ::Logging.ConsoleLogger, ::Base.CoreLogging.LogLevel, ::String, ::Module, ::Symbol, ::Symbol, ::String, ::Int64) (1 children)
                 14: signature Tuple{Type{UInt32}, Integer} triggered MethodInstance for VersionNumber(::Integer, ::Int64, ::Int64, ::Tuple{}, ::Tuple{}) (1 children)
                 15: signature Tuple{Type{Int32}, Integer} triggered MethodInstance for Int32(::Enum{T2} where T2<:Integer) (1 children)
                 16: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for var"#handle_message#2"(::Base.Iterators.Pairs{Symbol, _A, Tuple{Symbol}, _B} where {_A, _B}, ::typeof(Base.CoreLogging.handle_message), ::Base.CoreLogging.SimpleLogger, ::Base.CoreLogging.LogLevel, ::String, ::Module, ::Symbol, ::Symbol, ::String, ::Int64) (2 children)
                 17: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for rpad(::String, ::Integer, ::String) (5 children)
                 18: signature Tuple{Type{Int64}, Integer} triggered MethodInstance for var"#handle_message#2"(::Any, ::typeof(Base.CoreLogging.handle_message), ::Logging.ConsoleLogger, ::Base.CoreLogging.LogLevel, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any) (772 children)

Only the last one is bad, but it's pretty bad. (Among other things, it invalidates Base.require which means the next package to load has to recompile it.) The vulnerability arises in https://github.com/JuliaLang/julia/blob/b1fbe7f135f1b7b0f0c80a6484ea090eef1f1e72/stdlib/Logging/src/ConsoleLogger.jl#L104-L105; since all it knows is that maxlog isa Integer, it's possible that the (::Type{T})(x::ArrayInterface.StaticInt{N}) where {T<:Integer, N} method is applicable and more specific.

The only really good fix I can see is to restrict maxlog support to specific integer types, like the ones Julia ships with. This would be a fix for Base but I am filing this here as a reminder. When working on a package it's worth checking this periodically as it could be a deal-breaker for inclusion in other things.

Bug in `strides` of reinterpreted array with non-unit first stride

julia> using Revise, ArrayInterface, StaticArrays, BenchmarkHistograms; using ArrayInterface: StaticInt

julia> begin
           u_base = randn(1, 4, 4, 5)
           u_vectors = reshape(reinterpret(SVector{1, eltype(u_base)}, u_base),
                               Base.tail(size(u_base))...)
           u_view = view(u_vectors, 2, :, 3)
           u_view_reinterpreted = reinterpret(eltype(u_base), u_view)
           u_view_reshaped = reshape(u_view_reinterpreted, 1, length(u_view));
       end

julia> ArrayInterface.strides(u_view_reinterpreted)
(static(1),)

julia> ArrayInterface.strides(u_view)
(4,)

ArrayInterface.strides(u_view_reinterpreted) is wrong, it should have the same strides as u_view.

Promoting collaboration

It would be nice if more people could participate in development here, even if it was just reviewing PRs. Obviously we need documentation and examples to help people understand the scope of package, but we also have some things we've gravitated towards that aren't covered in ColProc.

I've come to think of these sorts of things before adding it here.

• Additions should typically only be dependent on trait types. This forces us to create an interface that is generic and robust to all the different array types we test.
• Would it make sense to eventually include this in base Julia or the standard library? Or would it be more suited as a package?

Thoughts?

Traits, functions, defaults

I'm excited someone is looking into this. :)

Apart from the design of a formal array interface (definition of traits such as mutability) I was wondering what the "meta" design decisions about how the traits are defined and discovered are.

Base has used "Tim Holy" traits in a number of places to help with e.g. iteration performation. I note the initial version here uses an ismutable function that returns a constant value.

I was wondering what the trade-offs are between type traits vs functions (and constant propagation / branch pruning) might be and what would be preferable. Are traits properties of types or of instances?

Also whether one wants "defaults" like ismutable(::AbstractArray) = true or false are desirable, or something to avoid. (One thing I like about the current system is one can create a simple abstract array with about 2 or 3 method definitions).

rename `has_sparsestruct` to `hassparsestruct`

I was reminded by JuliaDiff/BlueStyle#36

If this is to go into Base then the base convention is not to put _ after has
see:
hasfield haskey hasmethod hasproperty

Julia 1.3

Why does not this package accept Julia 1.3?

I'm a bit confused, but it looks like this causes Julia 1.3 to install the latest version of the package in Julia, without respecting any upper bound limit.
For instance, SparseDiffTools has an upper bound for ArrayInterface 2.0 but installing SparseDiffTools on my computer downloaded ArrayInterface 2.3.

Prettier printing for `OptionallyStaticUnitRange`?

julia> using PaddedMatrices

julia> @StrideArray rand(4,5)
4×5 StrideMatrix{Tuple{StaticInt{4}, StaticInt{5}}, (true, true), Float64, 1, 0, (1, 2), Tuple{StaticInt{8}, StaticInt{32}}, Tuple{StaticInt{1}, StaticInt{1}}, PaddedMatrices.MemoryBuffer{20, Float64}} with indices StaticInt{1}():StaticInt{1}():StaticInt{4}()×StaticInt{1}():StaticInt{1}():StaticInt{5}():
 0.301859  0.54085   0.222296  0.599135   0.767651
 0.336525  0.521614  0.120895  0.891117   0.58328
 0.469067  0.580745  0.745368  0.0832662  0.0853315
 0.489557  0.187212  0.228955  0.261394   0.632375

julia> StrideArray(rand(4,5))
4×5 StrideMatrix{Tuple{Int64, Int64}, (true, true), Float64, 1, 0, (1, 2), Tuple{StaticInt{8}, Int64}, Tuple{StaticInt{1}, StaticInt{1}}, Matrix{Float64}} with indices StaticInt{1}():StaticInt{1}():4×StaticInt{1}():StaticInt{1}():5:
 0.892622  0.862319   0.130142  0.308712   0.196662
 0.600617  0.0845287  0.752346  0.612079   0.534358
 0.545485  0.329272   0.305997  0.922645   0.650952
 0.618598  0.991797   0.60527   0.0707488  0.247227

The

with indices StaticInt{1}():StaticInt{1}():StaticInt{4}()×StaticInt{1}():StaticInt{1}():StaticInt{5}()
with indices StaticInt{1}():StaticInt{1}():4×StaticInt{1}():StaticInt{1}():5

Could be a lot prettier. Of course, I should probably make the rest of the StrideArray type print prettier, too.

AbstractRange traits

Would it be within the scope of this package to address range traits? I have several ideas that I figured would be more useful here than in a separate package.

Handling wrapped arrays

Today I gave up on something...

    function DiffEqBase.get_tmp(dc::DiffEqBase.DiffCache, u::LabelledArrays.LArray{T,N,D,Syms}) where {T,N,D,Syms}
      x = reinterpret(T, dc.dual_du.__x)
      LArray{T,N,D,Syms}(x)
    end

we need a better way to handle things like reinterpret and reshape on array types that are actually just wrapping another array type.

Move to JuliaArrays?

This has been brought up several times but most recently here https://discourse.julialang.org/t/converting-between-axisarray-like-packages/. I think it makes sense to move this package to JuliaArrays. I know that @timholy and @rafaqz are in favor.

@chriselrod , @ChrisRackauckas, @DilumAluthge. What do you think?

Should setindex(x::AbstractArray, v, i) support arbitrary indices?

https://github.com/JuliaDiffEq/ArrayInterface.jl/blob/c8513cbc0b9dc24e6674206255838abec553b854/src/ArrayInterface.jl#L22-L26

Why not remove type restriction from Base.setindex like this?

function Base.setindex(x::AbstractArray, v, I...)
  _x = copy(x)
  _x[I...] = v
  _x
end

Managing array memory

The formal interface for working with an array's memory is currently very limited. We mostly just have ways of accessing preallocated memory (e.g., unsafe_wrap and unsafe_pointer_to_objref). The rest of the time we need to rely on methods like push!/pop!/append!/resize to mutate a Vector down the line. Likewise, if we want to allocate memory for a new array we typically need to wrap a new instance of Array.

A couple reasons we this would be good to have here

1. Performance. The necessity to ensure that memory is safely managed also means that each call to something like resize! has to ensure that we don't create a situation where we are out of bounds. This probably only has minimal overhead in most cases but if we are doing something complicated like merging and or sorting two vectors then we may be calling resize! a lot. I recall this sort of thing coming up a lot when I was trying to work with graph algorithms last year that do a lot of insertion and deletion.
2. Necessary for creating some interfaces. One of the things that makes the indexing interface so nice is that there are formally defined places for new types to insert and propagate information (e.g., bounds checking). It also makes it easier to optimize indexing. In contrast, Vector directly uses methods that allocate and write to memory (e.g., Base._growend! and Base.arrayset!) while abstract types typically hope that resize! works and then use setindex!.

From where I'm standing it seams like this would be great to have here. The reason this is an issue instead of a PR is mainly because I'm unsure if the implementation for this would be too general or involved for this package. We could always define several methods here like unsafe_grow_end! and unsafe_shrink_end! but these aren't super helpful without implementations that can interact with pointers and references.

@chriselrod, do you think growing/shrinking/allocating memory in an efficient way would require a bunch of LLVM magic or could we do it pretty simply here?

Array Constructors - `similar`

I think we are in a good place to start working on array constructors (would also make documenting examples a whole lot easier). This brings up stuff related to similar:

• JuliaLang/julia#25107
• JuliaLang/julia#18161
• JuliaLang/julia#10064

I'm not convinced there's a silver bullet for array constructors, but I think we could at least find a solution to what similar is often trying to do, allow generic method definition without worrying about array specific constructors. I think we actually need to break up what similar does into several pieces though

1. allocate memory - a buffer/Ref that is safer than a pointer
2. do something to the memory - probably fastest on pointers
3. create user facing instance - e.g., Array

I haven't worked out all the details but here's some of the cleaner code I have so far that might support this.

""" allocate_memory(x, args...) """  # step 1 allocates memory with corresponding device to `x` and axes
function allocate_memory(x)  end

""" preserve(f, x, y) """  # step 2 operate on pointers and permit novel garbage collection approaches
preserve(f, x::X) where {X} = preserve(f, x, device(X))
preserve(f, x::X, ::CPUIndex) where {X} = f(x)
function preserve(f, x::X, ::CPUPointer) where {X}
    GC.@preserve x out = f(pointer(x))
    return out
end
preserve(f, x::X, y::Y) where {X,Y} = preserve(p_x -> preserve(Base.Fix1(op, p_x), y), x)


""" as_immutable """
function as_immutable(x) end

""" initialize(original, new_data) """ # step 3 turn processed buffer into user facing array
function initialize(x::X, data) where {X}
    if !ismutable(X)
        return as_immutable(data)
    else
        return data
    end
end

""" instantiate(f, x, allocator, initiializer) """
instantiate(f, x, allocator, initiializer) = initiializer(x, preserve(f, allocator(x),  x))
instantiate(f, x, allocator) = (f, x, allocator, initiialize)
instantiate(f, x) = (f, x, allocate_memory)

The reason I think separating things out like this is helpful is because it turns this function from base like this

function rotr90(A::AbstractMatrix)
    ind1, ind2 = axes(A)
    B = similar(A, (ind2,ind1))
    m = first(ind1)+last(ind1)
    for i=ind1, j=axes(A,2)
        B[j,m-i] = A[i,j]
    end
    return B
end

into this

function rotr90(A)
    ind1, ind2 = axes(A)
    return instantiate(A, x -> allocate_memory(x, (ind2, ind1))) do a, b
        m = first(ind1)+last(ind1)
        for i=ind1, j=axes(A,2)
            b[j,m-i] = a[i,j]
        end
    end
end

This means that new array types typically wouldn't need to change rotr90 but would just change their allocators and initializers. I'm not super familiar with how jagged arrays work but we could have allocate_memory take in device and memory layout info for this.

Error requiring `StaticArrays` from `ArrayInterface`: too many parameters for type

Seems to have just started happening with v3.1.8

julia> using StaticArrays

julia> using ArrayInterface
┌ Warning: Error requiring `StaticArrays` from `ArrayInterface`
│   exception =
│    too many parameters for type
│    Stacktrace:
│      [1] top-level scope
│        @ ~/.julia/packages/ArrayInterface/GV582/src/ArrayInterface.jl:930
│      [2] eval
│        @ ./boot.jl:360 [inlined]
│      [3] eval
│        @ ~/.julia/packages/ArrayInterface/GV582/src/ArrayInterface.jl:1 [inlined]
│      [4] (::ArrayInterface.var"#53#80")()
│        @ ArrayInterface ~/.julia/packages/Requires/7Ncym/src/require.jl:99
│      [5] err(f::Any, listener::Module, modname::String)
│        @ Requires ~/.julia/packages/Requires/7Ncym/src/require.jl:47
│      [6] (::ArrayInterface.var"#52#79")()
│        @ ArrayInterface ~/.julia/packages/Requires/7Ncym/src/require.jl:98
│      [7] withpath(f::Any, path::String)
│        @ Requires ~/.julia/packages/Requires/7Ncym/src/require.jl:37
│      [8] (::ArrayInterface.var"#51#78")()
│        @ ArrayInterface ~/.julia/packages/Requires/7Ncym/src/require.jl:97
│      [9] listenpkg(f::Any, pkg::Base.PkgId)
│        @ Requires ~/.julia/packages/Requires/7Ncym/src/require.jl:20
│     [10] macro expansion
│        @ ~/.julia/packages/Requires/7Ncym/src/require.jl:95 [inlined]
│     [11] __init__()
│        @ ArrayInterface ~/.julia/packages/ArrayInterface/GV582/src/ArrayInterface.jl:869
│     [12] _include_from_serialized(path::String, depmods::Vector{Any})
│        @ Base ./loading.jl:674
│     [13] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String)
│        @ Base ./loading.jl:760
│     [14] _require(pkg::Base.PkgId)
│        @ Base ./loading.jl:998
│     [15] require(uuidkey::Base.PkgId)
│        @ Base ./loading.jl:914
│     [16] require(into::Module, mod::Symbol)
│        @ Base ./loading.jl:901
│     [17] eval
│        @ ./boot.jl:360 [inlined]
│     [18] eval_user_input(ast::Any, backend::REPL.REPLBackend)
│        @ REPL /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:139
│     [19] repl_backend_loop(backend::REPL.REPLBackend)
│        @ REPL /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:200
│     [20] start_repl_backend(backend::REPL.REPLBackend, consumer::Any)
│        @ REPL /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:185
│     [21] run_repl(repl::REPL.AbstractREPL, consumer::Any; backend_on_current_task::Bool)
│        @ REPL /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:317
│     [22] run_repl(repl::REPL.AbstractREPL, consumer::Any)
│        @ REPL /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:305
│     [23] (::Base.var"#874#876"{Bool, Bool, Bool})(REPL::Module)
│        @ Base ./client.jl:387
│     [24] #invokelatest#2
│        @ ./essentials.jl:708 [inlined]
│     [25] invokelatest
│        @ ./essentials.jl:706 [inlined]
│     [26] run_main_repl(interactive::Bool, quiet::Bool, banner::Bool, history_file::Bool, color_set::Bool)
│        @ Base ./client.jl:372
│     [27] exec_options(opts::Base.JLOptions)
│        @ Base ./client.jl:302
│     [28] _start()
│        @ Base ./client.jl:485
└ @ Requires ~/.julia/packages/Requires/7Ncym/src/require.jl:49

(@v1.6) pkg> st
      Status `~/.julia/environments/v1.6/Project.toml`
  [4fba245c] ArrayInterface v3.1.8
  [90137ffa] StaticArrays v0.12.5

Goal: get to 100% code coverage for ArrayInterface.jl

For context, here is the current code coverage:

cc: @chriselrod

findstructralnz for BandedMatrix, BlockBandedMatrix, BandedBlockBandedMatrix

I think these complete the sparse differentiation story.

contiguous_axis_indicator for arrays with unknown contiguity

Currently, contiguous_axis returns nothing for an array type with unknown contiguity properties. However, for the same array type, contiguous_axis_indicator is undefined, as there is no defined contiguous_axis_indicator(::Nothing, ::Val) method.

Is this intentional? If not, would it make sense here to return nothing, just like for contiguous_axis?

Minimal example:

import ArrayInterface: contiguous_axis, contiguous_axis_indicator

# Define dummy array type with unknown contiguity properties
struct DummyArray{T,N} <: AbstractArray{T,N}
    dims :: Dims{N}
    DummyArray{T}(::UndefInitializer, dims...) where {T} = new{T,length(dims)}(dims)
end
Base.size(x::DummyArray) = x.dims
Base.getindex(::DummyArray{T}, inds...) where {T} = zero(T)

A = DummyArray{Float64}(undef, 3, 5)
contiguous_axis(A)            # returns nothing
contiguous_axis_indicator(A)  # MethodError: no method matching contiguous_axis_indicator(::Nothing, ::Val{2})

Register

@JuliaRegistrator register()

@JuliaRegistrator register subdir=lib/ArrayInterfaceBandedMatrices
@JuliaRegistrator register subdir=lib/ArrayInterfaceBlockBandedMatrices
@JuliaRegistrator register subdir=lib/ArrayInterfaceCUDA
@JuliaRegistrator register subdir=lib/ArrayInterfaceCore
@JuliaRegistrator register subdir=lib/ArrayInterfaceGPUArrays
@JuliaRegistrator register subdir=lib/ArrayInterfaceOffsetArrays
@JuliaRegistrator register subdir=lib/ArrayInterfaceStaticArrays
@JuliaRegistrator register subdir=lib/ArrayInterfaceTracker

Bug in `dense_dims` for reshaped views

Reported in #156 (comment):

julia> using ArrayInterface

julia> u_base = randn(10, 10); u_view = view(u_base, 3, :); u_reshaped_view = reshape(u_view, 1, size(u_base, 2));

julia> ArrayInterface.dense_dims(u_reshaped_view)

As confirmed by @chriselrod in #156 (comment), this is a bug:

The correct answer for u_reshaped_view would be (False(),False()).

StackOverflowError in DifferentialEquations

#144 did fix #142 here as well,
but did not help with SciML/DifferentialEquations.jl#740.
Is this relevant to this package ?

(try_array_interface) pkg> status ArrayInterface
      Status `~/share/prog/julia/dev/try_array_interface/Project.toml`
  [4fba245c] ArrayInterface v3.1.8 `https://github.com/JuliaArrays/ArrayInterface.jl.git#18aa87e`

Here is the new error message (still with julia-1.6.0)
(small change: "repeats 20127 times" instead of the previous 20123):

ERROR: StackOverflowError:
Stacktrace:
 [1] show(io::IOBuffer, x::Type)
   @ Base ./show.jl:814
 [2] print(io::IOBuffer, x::Type)
   @ Base ./strings/io.jl:35
 [3] print_to_string(xs::Type)
   @ Base ./strings/io.jl:135
 [4] string
   @ ./strings/io.jl:174 [inlined]
 [5] Symbol
   @ ./strings/basic.jl:229 [inlined]
 [6] issymbollike(x::Vector{Union{Nothing, Int64}})
   @ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/solutions/solution_interface.jl:294
 [7] getindex(A::ODESolution{Float64, 1, Vector{Float64}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Float64}}, ODEProblem{Float64, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, DiffEqBase.DEStats}, sym::Vector{Union{Nothing, Int64}}, args::Function)
   @ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/solutions/solution_interface.jl:56
 [8] getindex(A::ODESolution{Float64, 1, Vector{Float64}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Float64}}, ODEProblem{Float64, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, DiffEqBase.DEStats}, sym::Vector{Union{Nothing, Int64}}, args::Function) (repeats 20127 times)
   @ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/solutions/solution_interface.jl:69
 [9] top-level scope
   @ REPL[10]:1

---

Edit
With Julia-1.6.1 and the same project

(try_array_interface) pkg> status ArrayInterface
      Status `~/share/prog/julia/dev/try_array_interface/Project.toml`
  [4fba245c] ArrayInterface v3.1.8 `https://github.com/JuliaArrays/ArrayInterface.jl.git#18aa87e`

almost the same error, except back to "repeats 20123 times":

ERROR: StackOverflowError:
Stacktrace:
 [1] show(io::IOBuffer, x::Type)
   @ Base ./show.jl:814
 [2] print(io::IOBuffer, x::Type)
   @ Base ./strings/io.jl:35
 [3] print_to_string(xs::Type)
   @ Base ./strings/io.jl:135
 [4] string
   @ ./strings/io.jl:174 [inlined]
 [5] Symbol
   @ ./strings/basic.jl:229 [inlined]
 [6] issymbollike(x::Vector{Union{Nothing, Int64}})
   @ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/solutions/solution_interface.jl:294
 [7] getindex(A::ODESolution{Float64, 1, Vector{Float64}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Float64}}, ODEProblem{Float64, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, DiffEqBase.DEStats}, sym::Vector{Union{Nothing, Int64}}, args::Function)
   @ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/solutions/solution_interface.jl:56
 [8] getindex(A::ODESolution{Float64, 1, Vector{Float64}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Float64}}, ODEProblem{Float64, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5, OrdinaryDiffEq.InterpolationData{ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Float64}, Vector{Float64}, Vector{Vector{Float64}}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}}, DiffEqBase.DEStats}, sym::Vector{Union{Nothing, Int64}}, args::Function) (repeats 20123 times)
   @ SciMLBase ~/.julia/packages/SciMLBase/9EjAY/src/solutions/solution_interface.jl:69
 [9] top-level scope
   @ REPL[10]:1

iterate interface missing from structrualnz returns

@huanglangwen I assume the iterator interface may just be able to use the indexing?

UniqueVectors?

Does it make sense to include UniqueVectors here? It may be difficult to categorize; for instance, it is "mutable" but only in certain respects (e.g. it does not implement setindex!). Related, the traits defined here really ought to be defined very precisely.

Create checks to avoid a StackOverflowError

Hi!

I was using a structure inside a workspace for OrdinaryDiffEq and I was getting this error:

ERROR: StackOverflowError:
Stacktrace:
 [1] ismutable(::Type{T} where T) at /Users/ronisbr/.julia/packages/ArrayInterface/YFV07/src/ArrayInterface.jl:19 (repeats 79978 times)

It seems that something is not defined inside ArrayInterface. I can reproduce it with the following MWE:

julia> mutable struct A
       a
       end

julia> t = A(1)

julia> ArrayInterface.ismutable(t)
ERROR: StackOverflowError:
Stacktrace:
 [1] ismutable(::Type{T} where T) at /Users/ronisbr/.julia/packages/ArrayInterface/YFV07/src/ArrayInterface.jl:19 (repeats 79978 times)

Question: can we define a fallback function that prints an error when we do not know if the type is or isn't mutable?

ArrayInterface 3.1.7 causes `sum` to freeze on array with missings

On Julia 1.5.4 (and 1.5.2)

(@v1.5) pkg> add Feather, ArrayInterface@3.1.7

julia> using ArrayInterface, Feather
[ Info: Precompiling Feather [becb17da-46f6-5d3c-ad1b-1c5fe96bc73c]

julia> U2 = Matrix{Union{Missing, Float32}}(undef, 10, 10);

julia> my_sum(values::AbstractMatrix) = dropdims(sum(values, dims=2), dims=2)
my_sum (generic function with 1 method)

julia> my_sum(U2)
^C^C^C^C^C^C^C^CWARNING: Force throwing a SIGINT
^C^C^C^C^C^C^C^CWARNING: Force throwing a SIGINT
^C^C^C^C^CSegmentation fault (core dumped)

freezes forever. @etpinard (who is coauthor of this MWE, 👍) says that #137 is to blame. Not sure why Feather is necessary. @ExpandingMan, does it ring a bell? Importing Feather's dependencies doesn't cause the issue, so it should be an interaction with Feather's code...

Bug in `strides` of reshaped adjoint/transposed arrays

Reported by @chriselrod in #161 (comment):

julia> using ArrayInterface

julia> A = rand(5,12);

julia> B = reshape(A', (3,4,5));

julia> ArrayInterface.strides(B) # @chriselrod: should be `(5,15,StaticInt{1}())`
(static(1), 3, 12)