StridedArray should be an abstract class, not a union

[stevengj]

Currently, Julia defines

typealias StridedArray{T,N,A<:Array}  Union(Array{T,N}, SubArray{T,N,A})

It seems like it would be better to define

abstract StridedArray{T,N} <: AbstractArray{T,N}

and have Array and SubArray both inherit from StridedArray.

The reason is that I would like to define a new StridedArray type (wrapping a NumPy ndarray). There is no reason, as far as I can tell, that I can't provide all the functionality of a SubArray and be used everywhere the latter is usable, but currently it does not seem possible to do this. Changing to abstract StridedArray would fix this.

[JeffBezanson]

See also #987 .

[stevengj]

It would be great to have a whole hierarchy, but we don't need to have it all at once. We can start with StridedArray since this is used everywhere now (just as a Union) and add more subtypes later. This will also make a good laboratory to see if any warts appear when the change to an abstract type is made.

[ViralBShah]

See also #1276

[ViralBShah]

One difficulty here to even experiment is that AbstractArray and Array are defined in C.

Once jb/immutable is merged, perhaps we can have these definitions in abstractarray.jl and array.jl. Also, if AbstractArray can contain fields, we may be able to do away with StridedArray altogether and just have Array and SubArray.

[JeffBezanson]

No, there might be other StridedArray types like NumPyArray. We might also need a StridedSubArray that only views StridedArrays, and only with index types that preserve stridedness. There can be a more general SubArray that allows any indexes.

[ViralBShah]

Ok. How about getting rid of Array in favour of StridedArray as the default storage type? We could then have SubArray <: StridedArray <: AbstractArray.

[pao]

Renaming Array will break a lot of code and seems likely to confuse new users; also in the hierarchy you show that would make StridedArray abstract, wouldn't it?

[ViralBShah]

Just thinking if we can reduce one of these things. Array could be made to include the functionality of StridedArray, for example - in order to avoid changing the name.

[stevengj]

You need at least three types (regardless of what you call them):

• An AbstractArray type for any data structure that can be accessed as a[i1,i2,...], regardless of how it is stored (or whether it is stored at all, in the case of Range).
• An abstract StridedArray <: AbstractArray type, representing any data structure where a[i1,i2,...] is accessed as unsafe_ref(ptr, 1 + (i1-1)*stride1+(i2-1)*stride2+...) for some pointer and strides. This needs to be abstract so that different concrete types (e.g. wrappers for NumPy arrays) can be bolted on underneath.
• A concrete Array <: StridedArray type, which is an implementation of StridedArray with contiguous column-major data.

As much as possible, our routines should work with the most abstract possible type. For StridedArray types, we should have special-case code for many operations when the strides correspond to contiguous column-major order or contiguous row-major order (e.g. things like copying, scaling, etcetera, can be more efficient).

The LAPACK routines should work with StridedArray types where possible. Note that many of the LAPACK routines require the stride of the first dimension to be 1, but the stride of the second dimension can be arbitrary. Alternatively, they can also work if the stride of the second dimension is 1, as long as you interpret the data as transposed, possibly transposing the output if necessary. The routines could make an Array copy if the input is not in one of these two cases.

In order to avoid the need for explicit transposes in many cases, it might be nice to have a concrete TransposedArray <: StridedArray type that represents an Array with transposed (i.e. row-major/C-order, i.e. reversed dimensions) storage. But this is a second-order priority compared to making StridedArray abstract. (Once StridedArray is abstract and most operations work on StridedArray arguments, adding new StridedArray subclasses will be a lot less work.)

[timholy]

Thoughts triggered by #8235 (comment) reminded me that because StridedArray is a Union, it's not extensible. Perhaps we can fix this.

Example:

• an Image is a wrapper for an AbstractArray,

type Image{T,N,A<:AbstractArray} <: AbstractArrary{T,N}
    data::A
    properties::Dict
end

• I have real-world cases (most of the data collected in my lab) where there is reason not to restrict A to be a DenseArray, hence an Image can't be a DenseArray. For this reason, an Image is not a StridedArray.
• However, for most people, most of the time, an Image will in fact wrap a DenseArray. It seems to be a shame not to get "credit" for that, especially if one thinks about passing image data to C functions.

One can fix it this way: rewrite

function somefunction(A::StridedArray, ...)
    ...
end

as

somefunction(A::AbstractArray, ...) = _somefunction(A, ..., isstrided(A))
function _somefunction(A, ..., ::Type{TypeConst{true}})
    ...
end

using the definitions

isstrided(A::DenseArray) = TypeConst{true}
isstrided{T,N,A<:DenseArray}(A::View{T,N,A}) = TypeConst{true}
isstrided{T,N,A<:DenseArray}(img::Image{T,N,A}) = TypeConst{true}
isstrided(A) = TypeConst{false}

It's a little ugly, but perhaps with a @stridefunction macro it wouldn't be too bad, at least for functions that have all-strided arguments. (An alternative is #8154.)

Thoughts?

[stevengj]

@timholy, could you expand upon why your images cannot be DenseArrays? Note that DenseArray does not mean contiguous, it means "arrays that are laid out at regular offsets in memory" (i.e. strided).

[timholy]

Briefly, our raw data files represent 4d images (3d over time), but they require some computation to be performed upon them before you get out what you'd call an "image". So there's no physical memory address corresponding to a single voxel at a particular time.