Error calling `similar(Array{ComplexF32}, Float32)`

Internals & Design Internals
, ,
1

writing code that is agnostic about the array type but can still change the element_type, one wants to keep the array type but change the element-type for example to real(eltype(ArrayType)) for efficiency reasons.
The function similar() claims to provide just this according to the documentation. However, this feature seems to work only for instances but not for types. I.e.

A = Array{ComplexF32}
similar(A(undef, 0), real(eltype(A)))

is fine, whereas

A = Array{ComplexF32}
similar(A, real(eltype(A)))

causes an error (tested on Julia Versions 1.7.1 and 1.8.0).
Maybe this could be in the documentation or changed in the next Julia Version?

2

We could do something like the following.

julia> Base.similar(::Type{A}, ::Type{S}) where {T, A <: Array{T}, S} = Array{S}

julia> similar(Array{ComplexF64},Float64)
Array{Float64}

Perhaps someone could make this even more generic.

I’m trying to decide is this should still be similar.

1 Like
3

Creating an “similar” array from the original type and a new eltype or size would be a very useful feature. Without it, for example, it’s not possible to implement reduce(vcat) or stack that would work for empty arrays.
However, this feature isn’t implemented in julia for now.

4

Great way of doing this! Much nicer than my hacky implementation:

function get_real_arr_type(::Type{TA}) where {TA}
    typeof(similar(TA(undef, ntuple(x->0, ndims(TA))), real(eltype(TA))))
end
5

I seem to recall that @StefanKarpinski once suggested generalising our current idiom for creating uninitialized arrays. Something like ArrayType{ElementType}(value, dims...) which cleanly separates out the four different concerns. The same pattern could be extended for zero initialisation and such.

If that general pattern was implemented, you could to typeof(MyArray){Float64}(undef, size(MyArray)). In fact, that might already work, depending on your array type.

6

This seems like strange syntax, as it’s going to have two type parameters, the second of which becomes the eltype. At the moment that’s an error:

julia> x = StructArray(MtlArray(rand(Float32, 3)).+im);

julia> similar(x, Int32, (1, 4))
1Ă—4 MtlArray{Int32, 2}:
 0  0  0  0

julia> typeof(x)
StructVector{ComplexF32, NamedTuple{(:re, :im), Tuple{MtlArray{Float32, 1}, MtlArray{Float32, 1}}}, Int64} (alias for StructArray{Complex{Float32}, 1, NamedTuple{(:re, :im), Tuple{MtlArray{Float32, 1}, MtlArray{Float32, 1}}}, Int64})

julia> typeof(x){Int32}
ERROR: TypeError: in Type{...} expression, expected UnionAll, got Type{StructVector{ComplexF32

It’s for complicated wrapped types that similar is so useful, and (IMO) the extension to accept the type not an instance seems pretty natural.

7

Ah yes, you are right about the extra type parameter. Yeah, perhaps this is the role of similar. Although in this case it sounds like it should be a constructor? I mean “How do I make an instance given a type and its parameter” sounds like the definition of a constructor for me.

8

Yes as a generalisation of constructors like Array{Int}(undef, 3) to fill with other values, this sounds fine.

I guess similar(::Type, size...) is very close to being a constructor. The types must match (e.g. you can’t even change ndims, similar(typeof([1,2,3]), (4, 5)) is an error) and it doesn’t exist for many wrapped types like similar(typeof([1,2,3]'), (1, 5)). I’m not sure why it exists, in addition to Array{T}(undef, ...).

But similar(instance, eltype, size...) often makes something which is not of the type of the first argument, and usually recurses into complicated types. I wonder if there is any obstruction to generalising it to accept the type… does the type always give you enough information? Maybe it doesn’t indicate (say) which of several GPUs an object lives on?

9

Here’s a slight improvement. I’m not very convinced that this should be called similar, so I’ll call it seteltype for now.

seteltype(::Type{<: Array{T,N} where T}, ::Type{S}) where {S,N} = Array{S,N}
seteltype(::Type{<: Array}, ::Type{S}) where {S} = Array{S}

Here’s a demonstration.

julia> seteltype(Array, Float64)
Array{Float64}

julia> seteltype(Array{Int}, Float64)
Array{Float64}

julia> seteltype(Array{Int,5}, Float64)
Array{Float64, 5}
10

I see that the implementations of @mkitti contain the final type Array, which is not really what I intended it to be, yet it is debatable. E.g. similar(CuArray(1:10), Float32, 5) yields a 5-element CuArray{Float32} and not Array{Float32}.
I would prefer that similar(CuArray{Int64}, Float32, 5) yields also a 5-element CuArray{Float32}.
This could be implemented as follows:

function replace_eltype(::Type{TA}, ::Type{S}, num_dims=1) where {TA<:AbstractArray, S}
    typeof(similar(TA(undef, ntuple(x->0, num_dims)), S))
end

Base.similar(::Type{A}, ::Type{S}, dims) where {T, A <: AbstractArray{T}, S} = (replace_eltype(A, S, length(dims)))(undef, dims)

but its not a very nice implementation. If there was a Metaprogramming way to access the wrapping array-type, this may be implemented in a much nicer way.
Another bad part of this implementation is, that it does not allow the user the change the number of dimensions. I.e. similar(CuArray{Int, 2, CUDA.Mem.DeviceBuffer}, Float32, 6) yiels an error, whereas similar(CuArray{Int, 1, CUDA.Mem.DeviceBuffer}, Float32, 6) works.

But I guess the main question is, whether it is a good idea to assume by default that arrays are mutable or whether the default should indeed (as in @mkitti 's implementation) go to the type Array but then implementers of mutable arrays need to always care for this case individually.

11

The issue is that not all AbstractArrays necessarily have two parameters in the order we expect. I think it might be good to do this explicitly for each concrete AbstractArray implementation.

We could make some code to generically replace the first parameter of a type without any guarantees that it is actually the element type.

Another approach might be a parameter substitution approach e.g. swap Int64 for Float64.

12

I see. I guess, you mean the parameters of the type (written in curly brackets), which makes a lot of sense. But then the original way of implementing it via an empty array, is maybe not such a bad solution. It exploits that every implemented array should have a similar implemented and it makes sense that the type is then propagating the same way. Here is a slightly more complete solution which also works for AbstractArrays that do not have a dimension specified. See the implementation here:

These functions return a similar array type but using as TA, but eltype and ndims can be changed.
Does it make sense to have something like this in Base?