abstractarray.jl: hcat/vcat/cat are not sufficiently abstract?

[mcg1969]

I just ran into some of these issues and posted about it over on julia-users.

Given the intent to allow AbstractArray to be extremely flexible, some of the default implementations seem unnecessarily specialized and potentially very inefficient. For instance, in theory, hcat and vcat are just cat(2,...) and cat(1,...), perhaps with some exceptions to handle scalars thrown in. So an AbstractArray implementation could be nothing but a direct reference to cat. Instead, hcat and vcat (for example) have the following specializations:

function hcat{T}(V::AbstractVector{T}...)
function vcat{T}(V::AbstractVector{T}...)
function hcat{T}(A::AbstractVecOrMat{T}...)
function vcat{T}(A::AbstractMatrix{T}...)

I can totally understand why these specializations are worthwhile for DenseArray objects. For instance, in the first three of these specalizations, you can take advantage of the fact that with Fortran ordering the input matrices will be laid out sequentially with no interleaving. And indeed, if you look at the implementation of these methods, the implementations are taking advantage of this fact.

But what happens if the concrete matrix uses C ordering, or perhaps some sort of novel/compressed/distributed storage method? These implementations be suboptimal---and not just because of the layout assumption, but also because they use straight for-loops or list comprehensions to move elements one element at a time. Again, for Array objects, those loops are going to be nice and tight when compiled down. But consider how inefficient that is for sparse matrices, for instance.

I would suggest, therefore, that

• Methods that are tuned to be fast for StridedArray and/or DenseArray be moved up to those classes.
• Methosd that remain in AbstractArray should use block accesses as much as possible, so that subclasses with non-contiguous storage schemes are more likely to see decent performance.
• Find the solution to the heterogeneous array problem.

[JeffBezanson]

We have gone back and forth on this a bit. At one point we changed most of the AbstractArray definitions to apply only to Array for some of the reasons you cite. Then some of those were changed back. But you're right that *cat for AbstractArray should use block assignments.

I only see two that do not use block access:

function hcat{T}(V::AbstractVector{T}...)
function vcat{T}(V::AbstractVector{T}...)

[mcg1969]

It seems to me that the AbstractMatrix/AbstractVecOrMat implementations are unnecessarily fine-grained as well. The general case for cat does look fine.

[mcg1969]

As for cat_t, this gets to the heterogeneous array problem I was referring to earlier. Ideally, there should be an elegant way to decide dynamically whether or not to use an Array or some other non-standard storage type. One could conceive of a hack to similar that returns a different, non-contiguous array type for certain values of eltype. In that case, this code might actually not work as well:

    C = similar(full(A[1]), typeC, dimsC)
    range = 1
    for k=1:nargs
        nextrange = range+cat_ranges[k]
        cat_one = [ i != catdim ? (1:dimsC[i]) : (range:nextrange-1) for i=1:ndimsC ]
        C[cat_one...] = A[k]
        range = nextrange
    end
    return C

Consider what would happen, for instance, if C were a sparse matrix, and catdim=1. Now I believe you've already overloaded things to catch most of the sparse matrix cases, but I'm just offering that as an example of how this code might not in fact be the best way to go about it in the general abstract case.

So this is probably an appropriate default implementation for cat_t. Indeed, I actually think it's the best way to handle all of the special cases above, too. But there needs to be a way to override it cleanly. Currently I'm overloading Base.cat_t to detect arrays of my particular eltype. That's working pretty well for me.

[JeffBezanson]

The AbstractMatrix/AbstractVecOrMat implementations do a single assignment operation per array, as far as I can tell.

https://github.com/JuliaLang/julia/blob/master/base/abstractarray.jl#L582
https://github.com/JuliaLang/julia/blob/master/base/abstractarray.jl#L601

[mcg1969]

Ah yes, you're right. I think I just misread. For any Fortran ordering or compressed column storage format, this should work fine. Not so much for compressed row storage or other format that's not optimized for column wise access, but then we're basically back to my more esoteric issue similar to cat_t above.

[mcg1969]

Nice. Thanks!