Problem with nested array types in similar_arr_type

[RainerHeintzmann]

julia> similar_arr_type(typeof(reshape(view(rand(200), 1:100), 1, 100)))
ERROR: MethodError: no method matching Base.ReshapedArray{Float64, 2, SubArray{…}, Tuple{}}(::UndefInitializer, ::Tuple{Int64, Int64})

Closest candidates are:
  (::Type{Base.ReshapedArray{T, N, P, MI}} where {T, N, P<:AbstractArray, MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64}}}})(::Any, ::Any, ::Any)
   @ Base reshapedarray.jl:6

Stacktrace:
 [1] similar_arr_type(::Type{Base.ReshapedArray{Float64, 2, SubArray{…}, Tuple{}}}; dims::Int64, dtype::Type)
   @ NDTools C:\Users\pi96doc\.julia\packages\NDTools\Wk0Cg\src\type_tools.jl:85
 [2] similar_arr_type(::Type{Base.ReshapedArray{Float64, 2, SubArray{Float64, 1, Vector{…}, Tuple{…}, true}, Tuple{}}})
   @ NDTools C:\Users\pi96doc\.julia\packages\NDTools\Wk0Cg\src\type_tools.jl:84ulia> similar_arr_type(typeof(reshape(view(rand(200), 1:100), 1, 100)))
ERROR: MethodError: no method matching Base.ReshapedArray{Float64, 2, SubArray{…}, Tuple{}}(::UndefInitializer, ::Tuple{Int64, Int64})

A similar problem exists for the SubArray type.

julia> similar_arr_type(typeof(view(rand(200), 1:100)))
ERROR: MethodError: no method matching SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}(::UndefInitializer, ::Tuple{Int64})

Closest candidates are:
  SubArray{T, N, P, I, L}(::Any, ::Any, ::Any, ::Any) where {T, N, P, I, L}
   @ Base subarray.jl:19

Stacktrace:
 [1] similar_arr_type(::Type{SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}; dims::Int64, dtype::Type)
   @ NDTools C:\Users\pi96doc\.julia\packages\NDTools\Wk0Cg\src\type_tools.jl:85
 [2] similar_arr_type(::Type{SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}})
   @ NDTools C:\Users\pi96doc\.julia\packages\NDTools\Wk0Cg\src\type_tools.jl:84

... will try to fix it by introducing type specializations for these cases.

[roflmaostc]

We should delete those methods.

They are by design not type stable since the dimensionality depends on a value of the arguments:

julia> @code_warntype similar_arr_type(typeof(x); dtype=Float32)
MethodInstance for Core.kwcall(::@NamedTuple{dtype::DataType}, ::typeof(similar_arr_type), ::Type{Vector{Int64}})
  from kwcall(::NamedTuple, ::typeof(similar_arr_type), ::Type{TA}) where {T, N, TA<:AbstractArray{T, N}} @ NDTools ~/.julia/packages/NDTools/Wk0Cg/src/type_tools.jl:84
Static Parameters
  T = Int64
  N = 1
  TA = Vector{Int64}
Arguments
  _::Core.Const(Core.kwcall)
  @_2::@NamedTuple{dtype::DataType}
  @_3::Core.Const(NDTools.similar_arr_type)
  @_4::Core.Const(Vector{Int64})
Locals
  dims::Union{}
  dtype::Union{}
  @_7::Union{Int64, DataType}
Body::DataType
1 ─       Core.NewvarNode(:(dims))
│         Core.NewvarNode(:(dtype))
│         Core.NewvarNode(:(@_7))
│   %4  = Core.isdefined(@_2, :dims)::Core.Const(false)
└──       goto #3 if not %4
2 ─       Core.Const(:(@_7 = Core.getfield(@_2, :dims)))
└──       Core.Const(:(goto %9))
3 ┄       (@_7 = $(Expr(:static_parameter, 2)))
│   %9  = @_7::Core.Const(1)
│   %10 = Core.isdefined(@_2, :dtype)::Core.Const(true)
└──       goto #5 if not %10
4 ─       (@_7 = Core.getfield(@_2, :dtype))
└──       goto #6
5 ─       Core.Const(:(@_7 = $(Expr(:static_parameter, 1))))
6 ┄ %15 = @_7::DataType
│   %16 = (:dims, :dtype)::Core.Const((:dims, :dtype))
│   %17 = Core.apply_type(Core.NamedTuple, %16)::Core.Const(NamedTuple{(:dims, :dtype)})
│   %18 = Base.structdiff(@_2, %17)::Core.Const(NamedTuple())
│   %19 = Base.pairs(%18)::Core.Const(Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}())
│   %20 = Base.isempty(%19)::Core.Const(true)
└──       goto #8 if not %20
7 ─       goto #9
8 ─       Core.Const(:(Base.kwerr(@_2, @_3, @_4)))
9 ┄ %24 = NDTools.:(var"#similar_arr_type#5")(%9, %15, @_3, @_4)::DataType
└──       return %24


julia> @code_warntype similar_arr_type(typeof(x); dims=5, dtype=Float32)
MethodInstance for Core.kwcall(::@NamedTuple{dims::Int64, dtype::DataType}, ::typeof(similar_arr_type), ::Type{Vector{Int64}})
  from kwcall(::NamedTuple, ::typeof(similar_arr_type), ::Type{TA}) where {T, N, TA<:AbstractArray{T, N}} @ NDTools ~/.julia/packages/NDTools/Wk0Cg/src/type_tools.jl:84
Static Parameters
  T = Int64
  N = 1
  TA = Vector{Int64}
Arguments
  _::Core.Const(Core.kwcall)
  @_2::@NamedTuple{dims::Int64, dtype::DataType}
  @_3::Core.Const(NDTools.similar_arr_type)
  @_4::Core.Const(Vector{Int64})
Locals
  dims::Union{}
  dtype::Union{}
  @_7::Union{Int64, DataType}
Body::DataType
1 ──       Core.NewvarNode(:(dims))
│          Core.NewvarNode(:(dtype))
│          Core.NewvarNode(:(@_7))
│    %4  = Core.isdefined(@_2, :dims)::Core.Const(true)
└───       goto #3 if not %4
2 ──       (@_7 = Core.getfield(@_2, :dims))
└───       goto #4
3 ──       Core.Const(:(@_7 = $(Expr(:static_parameter, 2))))
4 ┄─ %9  = @_7::Int64
│    %10 = Core.isdefined(@_2, :dtype)::Core.Const(true)
└───       goto #6 if not %10
5 ──       (@_7 = Core.getfield(@_2, :dtype))
└───       goto #7
6 ──       Core.Const(:(@_7 = $(Expr(:static_parameter, 1))))
7 ┄─ %15 = @_7::DataType
│    %16 = (:dims, :dtype)::Core.Const((:dims, :dtype))
│    %17 = Core.apply_type(Core.NamedTuple, %16)::Core.Const(NamedTuple{(:dims, :dtype)})
│    %18 = Base.structdiff(@_2, %17)::Core.Const(NamedTuple())
│    %19 = Base.pairs(%18)::Core.Const(Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}())
│    %20 = Base.isempty(%19)::Core.Const(true)
└───       goto #9 if not %20
8 ──       goto #10
9 ──       Core.Const(:(Base.kwerr(@_2, @_3, @_4)))
10 ┄ %24 = NDTools.:(var"#similar_arr_type#5")(%9, %15, @_3, @_4)::DataType
└───       return %24

[RainerHeintzmann]

I see. We should introduce the Val type for dims and than it is fine. I can add this.

[RainerHeintzmann]

... it actually is already working. No need to change. But we should mention this in the documentation.

@code_warntype similar_arr_type(typeof(rand(100)), dims=Val(2))
MethodInstance for Core.kwcall(::@NamedTuple{dims::Val{2}}, ::typeof(similar_arr_type), ::Type{Vector{Float64}})
  from kwcall(::NamedTuple, ::typeof(similar_arr_type), ::Type{TA}) where {T, N, TA<:AbstractArray{T, N}} @ NDTools C:\Users\pi96doc\Documents\Programming\Julia\NDTools.jl\src\type_tools.jl:84
Static Parameters
  T = Float64
  N = 1
  TA = Vector{Float64}
Arguments
  _::Core.Const(Core.kwcall)
  @_2::Core.Const((dims = Val{2}(),))
  @_3::Core.Const(NDTools.similar_arr_type)
  @_4::Core.Const(Vector{Float64})
Locals
  dims::Union{}
  dtype::Union{}
  @_7::Union{Val{2}, Type{Float64}}
Body::Type{Matrix{Float64}}
1 ─       Core.NewvarNode(:(dims))
│         Core.NewvarNode(:(dtype))
│         Core.NewvarNode(:(@_7))
│   %4  = Core.isdefined(@_2, :dims)::Core.Const(true)
└──       goto #3 if not %4
2 ─       (@_7 = Core.getfield(@_2, :dims))
└──       goto #4
3 ─       Core.Const(:(@_7 = $(Expr(:static_parameter, 2))))
4 ┄ %9  = @_7::Core.Const(Val{2}())
│   %10 = Core.isdefined(@_2, :dtype)::Core.Const(false)
└──       goto #6 if not %10
5 ─       Core.Const(:(@_7 = Core.getfield(@_2, :dtype)))
└──       Core.Const(:(goto %15))
6 ┄       (@_7 = $(Expr(:static_parameter, 1)))
│   %15 = @_7::Core.Const(Float64)
│   %16 = (:dims, :dtype)::Core.Const((:dims, :dtype))
│   %17 = Core.apply_type(Core.NamedTuple, %16)::Core.Const(NamedTuple{(:dims, :dtype)})
│   %18 = Base.structdiff(@_2, %17)::Core.Const(NamedTuple())
│   %19 = Base.pairs(%18)::Core.Const(Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}())
│   %20 = Base.isempty(%19)::Core.Const(true)
└──       goto #8 if not %20
7 ─       goto #9
8 ─       Core.Const(:(Base.kwerr(@_2, @_3, @_4)))
9 ┄ %24 = NDTools.:(var"#similar_arr_type#5")(%9, %15, @_3, @_4)::Core.Const(Matrix{Float64})
└──       return %24

[RainerHeintzmann]

here is the new code:

function similar_arr_type(::Type{TA}; dims=N, dtype=T) where {T, N, TA<:AbstractArray{T,N}}
    typeof(similar(TA(undef, ntuple(x->0, N)), dtype, ntuple(x->0, dims)))
end

function similar_arr_type(::Type{TA}; dims=N, dtype=T) where {T, N, P, I, L, TA<:SubArray{T,N,P,I,L}}
    similar_arr_type(P, dims=dims, dtype=dtype)
end

function similar_arr_type(::Type{TA}; dims=N, dtype=T) where {T, N, P, MI, TA<:Base.ReshapedArray{T,N,P,MI}}
    similar_arr_type(P, dims=dims, dtype=dtype)
end

# note that T refers to the new type (if not explicitely specified) and therefore replaces the eltype of the array as defined by P
function similar_arr_type(::Type{TA}; dims=N, dtype=T) where {T, N, O, P, B, TA<:Base.ReinterpretArray{T,N,O,P,B}}
    similar_arr_type(P, dims=dims, dtype=dtype)
end

function similar_arr_type(::Type{TA}; dims=Val(1), dtype=eltype(TA)) where {TA<:AbstractArray}
    typeof(similar(TA(undef), dtype, ntuple(x->0, dims)))
end

can I push on the main branch?

[RainerHeintzmann]

    @test similar_arr_type(Array{ComplexF64,1}, dims=Val(2), dtype=Int) == Matrix{Int}
    @test similar_arr_type(typeof(view(ones(10,10),2:5,2:5)), dims=Val(1)) == Vector{Float64}
    @test similar_arr_type(typeof(reinterpret(Int, ones(10))), dims=Val(2), dtype=Float32) == Matrix{Float32}
    @test similar_arr_type(typeof(reshape(view(ones(25),1:25), 5,5)), dims=Val(1), dtype=Int) == Vector{Int}

are running fine.

[RainerHeintzmann]

Here is an argument why this functionality is needed:
Some toolboxes (e.g. SeparableFunctions.jl) create objects. They are often parameterized by types (not objects).
Similar has only limited type changing capabilities, and they recently added one type-based functionality, based on the code in NDTools.jl I would guess.

[roflmaostc]

Which code has been recently added?
I highly doubt many people are using this 😅

Can you send a concrete example?

[roflmaostc]

Make a PR and then I review parts of the code.

[RainerHeintzmann]

CuArrays and probably many other types are no problem, since it has a type-based constructor of the form TA(undef, Tuple)
However, ReshapedArrays and SubArray and Base.ReinterpretArray do not.
I make the PR now.

[RainerHeintzmann]

Which code has been recently added? I highly doubt many people are using this 😅

Can you send a concrete example?

Not sure, if recently added, but this is the documentation:

 similar(storagetype, axes)


  Create an uninitialized mutable array analogous to that specified by storagetype, but with axes specified by the last argument.

  Examples:

  similar(Array{Int}, axes(A))


  creates an array that "acts like" an Array{Int} (and might indeed be backed by one), but which is indexed identically to A. If A has conventional indexing, this will be identical to  
  Array{Int}(undef, size(A)), but if A has unconventional indexing then the indices of the result will match A.

but it has the same bug:

julia> similar(typeof(A), axes(A));

julia> A = ones(10,10)[1:20];

julia> similar(typeof(A), axes(A)); # works fine

julia> A = @view ones(10,10)[1:20];

julia> similar(typeof(A), axes(A));
ERROR: MethodError: no method matching SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}(::UndefInitializer, ::Tuple{Int64})

Closest candidates are:
  SubArray{T, N, P, I, L}(::Any, ::Any, ::Any, ::Any) where {T, N, P, I, L}
   @ Base subarray.jl:19

Stacktrace:
 [1] similar(::Type{SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}, dims::Tuple{Int64})
   @ Base .\abstractarray.jl:877
 [2] similar(::Type{SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}}, shape::Tuple{Base.OneTo{Int64}})
   @ Base .\abstractarray.jl:876
 [3] top-level scope
   @ REPL[57]:1

[roflmaostc]

Why not similar(A, axes(A))?

[RainerHeintzmann]

See example above. It also does not work. In addition it is nice to be able to change the dimensions of the array and the type.
similar needs to axes to agree to the type.

julia> A = ones(10,10);

julia> similar(typeof(A), (Base.OneTo(30), ))
ERROR: MethodError: no method matching Matrix{Float64}(::UndefInitializer, ::Tuple{Int64})

[roflmaostc]

Use similar(A, (Base.OneTo(30),)), you get what you need.

[RainerHeintzmann]

As mentioned above, there are generator functions, which do not have a concrete object to apply similar to. This is exactly what the similar_arr_type code does and it does allow to change dimensions as well as base type. Yes, you can do it by copying the code inside the functions, but since it is not trivial to do this, its nice to have it encapsulated in such a function.