RFC/WIP: Support array type as input for functions returning AbstractArray instance

[martinholters]

This is a shot at #11557.

At the moment, only zeros, ones, and eye are supported. For example

eye(Bidiagonal, 4, 4)        # Bidiagonal{Float64}
eye(Bidiagonal{Int64}, 4, 4) # Bidiagonal{Int64}

will return 4×4 identity matrices of the indicated types. (EDIT: Changed example to use Bidiagonal.) There is still some work ahead, but I'd appreciate feedback on the approach taken before committing more time to this. To avoid code duplication, I have added Type(dims...) like constructors for all array types where it makes sense and then implemented zeros, ones, and eye for {T<:AbstractArray}(::Type{T}, dims...) using these.

To exercise things, spzeros now simply falls back to zeros while speye forwards to a custom eye(::Type{SparseMatrixCSC}, dims...) (which is more efficient than the generic implementation). So far, no syntax or functionality is invalidated/removed, with one exception: spzeros(Any, ...) now fails with ERROR: MethodError: no method matching zero(::Type{Any}), while previously it succeeded in constructing a matrix, which however was of rather limited use.

In addition to extending this to more functions (suggestions welcome) and adding tests and docs, there is a more conceptual aspect that might need some work: the way missing type parameters are handled. Presently, all new constructors just use Float64 as default eltype which just happens to be the right thing to do for these functions, but seems a bit arbitrary and requires special-casing Array, which defaults to Array{Any}. I'm thinking about something like fixate_eltype{T<:AbstractArray}(::Type{T}, default::Type) that would be implemented by array types like e.g.

fixate_eltype(::Type{SparseMatrixCSC}, default::Type) = SparseMatrixCSC{default,Int}
fixate_eltype{Tv}(::Type{SparseMatrixCSC{Tv}}, default::Type) = SparseMatrixCSC{Tv,Int}
fixate_eltype{Tv,Ti}(::Type{SparseMatrixCSC{Tv,Ti}}, default::Type) = SparseMatrixCSC{Tv,Ti}

Thoughts?

[andreasnoack]

We might need to distinguish between subtypes of AbstractArray that defines storage schemes, e.g. Array, SparseMatrixCSC, DArray etc. and view like types like SubArray, XTriangular, and Symmetric. The latter are parametric on an array type. Of course, they could also have a default array type like in your example where eye(LowerTriangular, 4, 4) defaults to Array{Float64,2} but, right now, I think it makes more sense to use the constructors for the view types, i.e. LowerTriangular(eye(Matrix{Float64,2},4,4)) instead which probably would have the shorter version LowerTriangular(eye(4)).

[martinholters]

@andreasnoack You're right, of course. I meant to exclude these view types (and did for SubArray, ReshapedArray, ...), but somehow was too focused on other things to recognize XTriangular as such.

[tkelman]

What is fixate supposed to mean?

Sparse matrices with Any eltype are useful once in a while, it would be nice not to lose the concise way of initializing them. Why does spzeros have to call zero on the eltype?

[martinholters]

The term fixate is surely questionable. I'm trying to say "if there is no eltype yet, use the given default".

The generic zeros fills with zero(Tv), but one could surely provide a specialized zeros(::Type{SparseMatrixCSC}, dims...). It's a bit of a misnomer, though, as the unset elements are not zero - accessing them gives an error. Is SparseMatrixCSC{Any}(m, n) so much worse?

[tkelman]

SparseMatrixCSC{Any}(m, n) wouldn't be so bad, but there's no such method currently. I see that you're adding one here.

[tkelman]

this should probably check for negative n instead of allowing n not to match length(colptr)-1

[martinholters]

That duplicates the check from the constructor above which it delegates to, but is probably better than the current code, I agree. Maybe even better to have a dedicated function to check all sizes are non-negative, as this check should occur in all the array constructors at some point.

[martinholters]

Anyone got an idea about the segfault on travis?

[martinholters]

Updated with the feedback I got so far (except for #16740 (comment), du2 of Tridiagonal).

[tkelman]

this seems wrong - it's same size, not same nonzero structure, right?

[martinholters]

Right. Same nonzero structure wouldn't make much sense anyway. Good opportunity to correct that.

[tkelman]

with explicit zeros I guess it could, but I think Y=copy(S); fill!(Y.nzval, 0); Y is fine for that at the moment

[martinholters]

You are thinking of spzeros now, aren't you? That doesn't support the spzeros(S) signature. (zeros(S) works, but doesn't preserve to nonzero structure, either).

[tkelman]

yes of course, don't review without coffee

[martinholters]

Updated again, incorporating remaining feedback.

New functionality is now having a fill that takes the return type:

fill(SparseMatrixCSC, 2, (3,3))          # gives a SparseMatrixCSC{Int64,Int64}
fill(SparseMatrixCSC{Float64}, 2, (3,3)) # gives a SparseMatrixCSC{Float64,Int64}

For this to work I need the fixate_eltype function (suggestions for better names welcome!), which I have implemented in a generic fashion, but it has become quite a monstrosity. Is there an easier way to achieve this?

I've also discovered an inconsistency of fill!, not introduced by this PR, but easier to demonstrate:

julia> fill(Diagonal, 1, (3,3))
3×3 Diagonal{Int64}:
 1  ⋅  ⋅
 ⋅  1  ⋅
 ⋅  ⋅  1

julia> fill(Bidiagonal, 1, (3,3))
ERROR: ArgumentError: Array A of type Bidiagonal{Int64} and size (3,3) can
    not be filled with x=1, since some of its entries are constrained.
 in fill!(::Bidiagonal{Int64}, ::Int64) at ./linalg/bidiag.jl:515
 in fill(::Type{Bidiagonal}, ::Int64, ::Tuple{Int64,Int64}) at ./abstractarray.jl:382
 in eval(::Module, ::Any) at ./boot.jl:225
 in macro expansion at ./REPL.jl:92 [inlined]
 in (::Base.REPL.##1#2{Base.REPL.REPLBackend})() at ./event.jl:46

julia> fill(Tridiagonal, 1, (3,3))
ERROR: ArgumentError: Array A of type Tridiagonal{Int64} and size (3,3) can
    not be filled with x=1, since some of its entries are constrained.
 in fill!(::Tridiagonal{Int64}, ::Int64) at ./linalg/bidiag.jl:515
 in fill(::Type{Tridiagonal}, ::Int64, ::Tuple{Int64,Int64}) at ./abstractarray.jl:382
 in eval(::Module, ::Any) at ./boot.jl:225
 in macro expansion at ./REPL.jl:92 [inlined]
 in (::Base.REPL.##1#2{Base.REPL.REPLBackend})() at ./event.jl:46

I'd say both behaviors may be desirable at times, I have no clear favorite, but I'd like this to be consistent. Are the "historical reasons" still relevant? Maybe export fillslots! (renamed to fillnz!?) and support it also for sparse matrices? CC @jw3126

[tkelman]

I've been meaning to report and/or fix that exact issue, essentially that fill! does not respect structural constraints - think about fill! on a UnitTriangular for example. Let's move that discussion out of this PR.

[martinholters]

Hmmm, segfault in sig_match_simple on travis again. I will rebase to see if that does any good and open an issue otherwise, unless someone can share any helpful insights here.

[martinholters]

Next bit is in place: typed_*cat and cat_t. This makes e.g.

Array[speye(10,10) ones(10)]
SparseMatrixCSC[speye(10,10) ones(10)]

return a 10×11 matrix of type Matrix{Float64} and SparseMatrixCSC{Float64,Int}, respectively.
However, I have not touched getindex yet, so

SparseVector[spzeros(10), ones(10)]

returns a Vector of two SparseVectors. The reason for this is that the getindex syntax is so often used to construct Arrays of AbstractArray subtypes and I tried not to break anything for now. The inconsistency thus introduced is, however, not nice either (especially for the zero- and one-argument case, where the syntactic difference is not that obvious...) Thoughts?

EDIT: Looks like I messed up the commits I have pushed somewhere and trying to fix it in a hurry didn't work (as usual), so it's a bit early for code review, but please share your thoughts about the general idea.
EDIT2: I think I got it now, hoping CI agrees....

[martinholters]

Except for everything *cat, tests and docs have been added. For the *cat changes, I'd appreciate feedback whether we even want to go in that direction, see remarks above. Apart from that, the code could probably use some review. (I'm sure all the rebasing did leave some unwanted traces...)

[pabloferz]

Love the idea. Maybe we can do it even more general? E.g.

function set_eltype{A<:AbstractArray}(::Type{A}, T::Type)
    Base.@_pure_meta
    AA = A
    AT = AA.name.primary
    while AT.name.primary != AbstractArray
        AA = supertype(AA)
        AT = supertype(AT)
    end
    if AA.parameters[1] !== T
        tvarname = AT.parameters[1].name
        params = A.parameters
        new_params = collect(A.name.primary.parameters)
        for i in eachindex(new_params)
            p = new_params[i]
            if isa(p, TypeVar) && p.name == tvarname
                new_params[i] = T
            else
                new_params[i] = params[i]
            end
        end
        return A.name.primary{new_params...}
    else
        return A
    end
end

so it also gives set_type(Diagonal{Int}, Float64) === Diagonal{Float64})

[martinholters]

I deliberately wanted fixate_eltype(Diagonal{Int}, Float64) === Diagonal{Int} to be able to specify a default element type. E.g. zeros(Diagonal, 2, 2) should give a Diagonal{Float64}, but zeros(Diagonal{Int}, 2, 2) should give a Diagonal{Int}. Maybe a third, boolean argument could be used to choose between "set always" and "set only if unset".

[pabloferz]

I see. I had in mind another use case, using something like this combined with something like deepcopy_internal to write a generic similar(::AbstractArray, ::Type). Handling the dimensions in a generic fashion might not be possible and we would still need specialization (as you have here).

With this in mind and what is needed here, a third boolean argument would seem like the way to cover both scenarios.

[vtjnash]

This function is fundamentally invalid. T.name.primary isn't a meaningful type and its type parameters don't have any fixed relation to the type parameters of AbtractArray. If you can't express a type relationship via dispatch, that is highly indicative that the relation you want is not valid.

I think you want to call the similar function, whose definition is similar to what it appears the function above was attempting to do.

[pabloferz]

I guess it could be something like

function set_eltype{A<:AbstractArray}(::Type{A}, T::Type, force=false)
    Base.@_pure_meta
    SA = A
    AA = A.name.primary
    while SA.name.primary != AbstractArray
        SA = supertype(SA)
        force && (AA = supertype(AA))
    end
    if (!force && !isa(SA.parameters[1], TypeVar)) ||
       (force && SA.parameters[1] === T)
        return A
    end
    varname = (force ? AA : SA).parameters[1].name
    params = A.parameters
    new_params = collect(force ? A.name.primary.parameters : params)
    for i in eachindex(new_params)
        p = new_params[i]
        if isa(p, TypeVar) && p.name == varname
            new_params[i] = T
        elseif force
            new_params[i] = params[i]
        end
    end
    return A.name.primary{new_params...}
end

This would give set_eltype(Diagonal, Float64) === Diagonal{Float64}, set_eltype(Diagonal{Int}, Float64) === Diagonal{Int}, set_eltype(Diagonal{Int}, Float64, true) === Diagonal{Float64}

[nalimilan]

I've often wished that similar(::Type, ...) returned a type rather than an instance, and that AbstractArray subtypes implemented only it. Then a generic fallback would automatically return instances when passed instances. Wouldn't that be all we need here?

[martinholters]

Thanks for sharing these insights, @vtjnash, it certainly is a valuable contribution to the manual. But I still don't see why the fixate_eltype from this PR would be invalid. As said, it seems to work as I want it:

julia> fixate_eltype(Array, Float16)
Array{Float16,N}

julia> fixate_eltype(Array{Int}, Float16)
Array{Int64,N}

julia> fixate_eltype(Vector, Float16)
Array{Float16,1}

julia> fixate_eltype(Vector{Int}, Float16)
Array{Int64,1}

julia> fixate_eltype(BitArray, Float16)
BitArray{N}

julia> fixate_eltype(BitVector, Float16)
BitArray{1}

julia> type FunkyArray{NumDims, Foo, ElType, ElNiño} <: AbstractArray{ElType, NumDims}; end

julia> fixate_eltype(FunkyArray, Float16)
FunkyArray{NumDims,Foo,Float16,ElNiño}

julia> fixate_eltype(FunkyArray{42}, Float16)
FunkyArray{42,Foo,Float16,ElNiño}

julia> fixate_eltype(FunkyArray{42, :bar}, Float16)
FunkyArray{42,:bar,Float16,ElNiño}

julia> fixate_eltype(FunkyArray{42, :bar, String}, Float16)
FunkyArray{42,:bar,String,ElNiño}

julia> fixate_eltype(FunkyArray{42, :bar, String, -1}, Float16)
FunkyArray{42,:bar,String,-1}

There is one corner case I could find:

julia> type UntypedArray <: AbstractArray; end

julia> fixate_eltype(UntypedArray, Float16)
UntypedArray

I'm not sure what would be desired in this case, though. Maybe throwing instead of returning a type with still unspecified element type would be better. Anyway, this kind of type does not look exactly useful, so I'm not worried too much.

T.name.primary isn't a meaningful type and its type parameters don't have any fixed relation to the type parameters of AbtractArray.

That sounds indeed worrying for the proposed fixate_eltype. However, T.name.primary is used in a couple of places around base, so there has to be something of value in it. Also, the relationship of the type parameters between, say, Array and AbstractArray has be to recorded somewhere. Why shouldn't we be able to exploit it here?

[nalimilan]

@vtjnash Comments?

[andreasnoack]

I think the conclusion is that this should be implemented with dispatch cf. the new paragraphs in the documentation and that T.name.primary should generally be avoided. The places in Base where T.name.primary are few and special.

[martinholters]

OTOH https://github.com/JuliaLang/julia/blob/7f230adcf13d7a2bcb76e55327f211b7b4f25046/doc/devdocs/types.rst#typenames reads as though the use of T.name.primary here is perfectly valid. (Maybe the relationship of the TypeVars is a real problem?) It feels a bit unsatisfying to require a method to be implemented for every array type that should support e.g. zeros(MyArrayType, (5, 3)) when there is this generic alternative available.

[tkelman]

bump, any plan to revisit?

[martinholters]

Not any time soon (for lack of time, I'm afraid), so if someone wants to take over, go ahead! Otherwise, I might come back to this eventually.

[martinholters]

This is mostly superseded. This main thing in this PR that might still be relevant is the change to the concatenation functions, but I'm not sure how much they were liked anyway.