SubArray Supertype

[JaredCrean2]

SubArrays are not subtypes of DenseArray, even when they could be:

julia> arr = rand(3,3)
3x3 Array{Float64,2}:
 0.838499  0.649594  0.810356
 0.865377  0.204888  0.224049
 0.842879  0.747802  0.67063

julia> s_arr = sub(arr, :, 1)
3-element SubArray{Float64,1,Array{Float64,2},Tuple{Colon,Int64},2}:
 0.838499
 0.865377
 0.842879

julia> super(typeof(s_arr))
AbstractArray{Float64,1}

This is a problem when trying to pass the SubArray to a function that only accepts DenseArrays.

[timholy]

This can't be solved elegantly until we have traits; once we do, you could write

foo(A::AbstractArray::IsDense) = ...

and it will Just Work.

But you can do this now if you're willing to go to a bit of effort in your dispatch hierarchy:

abstract DenseKind
immutable IsDense <: DenseKind end
immutable NotDense <: DenseKind end

DenseKind(::DenseArray) = IsDense()
DenseKind(::AbstractArray) = NotDense()
DenseKind{T,N,P<:DenseArray}(A::Base.FastContiguousSubArray{T,N,P}) = IsDense()

foo(A::AbstractArray) = _foo(DenseKind(A), A)
_foo(::IsDense, A::AbstractArray) = "I'm really dense"
_foo(::NotDense, A::AbstractArray) = "I'm a bit fluffy"

Testing:

julia> A = rand(3,3)
3×3 Array{Float64,2}:
 0.8386    0.203996  0.542984
 0.78251   0.573883  0.377574
 0.633186  0.622487  0.253908

julia> S = view(A, :, 1:2);

julia> S2 = view(A, [1,3], 1:2);

julia> foo(A)
"I'm really dense"

julia> foo(S)
"I'm really dense"

julia> foo(S2)
"I'm a bit fluffy"

So it's a bit unwieldy, but it works.

[andreasnoack]

Could the method in question be changed to work for StridedArrays instead of DenseArrays?

[JaredCrean2]

Could the method in question be changed to work for StridedArrays instead of DenseArrays?

No, it has to be DenseArray

So it's a bit unwieldy, but it works.

That will work for now. Leaving the issue open until traits are implemented?

[andreasnoack]

No, it has to be DenseArray

Just curious. Why is that? StridedArray is the union of DenseArray and SubArray of a DenseArray so StridedArray seems like the obvious candidate for your signature.

[JaredCrean2]

At one point the SubArray gets passed into C like this:

function (arr::DenseArray)
  ccall( ...stuff.., arr, length(arr))
end

If the underlying array isn't contiguous in memory, the C code can't use it.

[andreasnoack]

Wouldn't you need to check the strides of the DenseArray anyway before passing it to ccall?

[JaredCrean2]

I just checked the docs, and it looks like you are right. Now I am a bit confused though. I thought DenseArray implied all strides = 1, and StridedArray implied non-unit strides. According to the docs, DenseArray implies a general strided memory layout, and StridedArray is defined as you described. But then a StridedArray does not imply a strided memory layout because this is valid:

julia> a = rand(4,4)
4x4 Array{Float64,2}:
 0.706608   0.0505246  0.951182  0.678708
 0.0962837  0.321637   0.224594  0.887127
 0.406915   0.061083   0.554764  0.170252
 0.130242   0.943475   0.438224  0.647376

julia> v = [1, 3, 4]
3-element Array{Int64,1}:
 1
 3
 4

julia> sv = sub(a, v, 2:3)
3x2 SubArray{Float64,2,Array{Float64,2},Tuple{Array{Int64,1},UnitRange{Int64}},0}:
 0.0505246  0.951182
 0.061083   0.554764
 0.943475   0.438224

Interesting, further experimentation shows:

julia> stride(sv, 1)
ERROR: ArgumentError: strides valid only for RangeIndex indexing
 in strides at subarray.jl:359
 in stride at subarray.jl:375

julia> pointer(sv)
ERROR: MethodError: `convert` has no method matching convert(::Type{Ptr{Float64}}, ::SubArray{Float64,2,Array{Float64,2},Tuple{Array{Int64,1},UnitRange{Int64}},0})
This may have arisen from a call to the constructor Ptr{Float64}(...),
since type constructors fall back to convert methods.

julia> sv2 = sub(a, :, 1)
4-element SubArray{Float64,1,Array{Float64,2},Tuple{Colon,Int64},2}:
 0.706608 
 0.0962837
 0.406915 
 0.130242 

julia> pointer(sv2)
Ptr{Float64} @0x00007f905759b280

However

julia> v2 = 1:2:4
1:2:3

julia> sub(a, v2, 2:3)
2x2 SubArray{Float64,2,Array{Float64,2},Tuple{StepRange{Int64,Int64},UnitRange{Int64}},1}:
 0.0505246  0.951182
 0.061083   0.554764

julia> pointer(ans)
Ptr{Float64} @0x00007f905759b2a0

So it seems there is no abstract type that implies contiguous memory layout, so the next best thing is to use StridedArray and check the stride values at runtime?

[andreasnoack]

Notice that sv is not a StridedArray. I think the best way for now is to
use StridedArray and check the strides. Maybe there is a better way and
maybe we'll need to revisit this when we get a buffer type .

On Friday, August 19, 2016, Jared Crean notifications@github.com wrote:

I just checked the docs, and it looks like you are right. Now I am a bit
confused though. I thought DenseArray implied all strides = 1, and
StridedArray implied non-unit strides. According to the docs
http://docs.julialang.org/en/release-0.4/manual/arrays/, DenseArray
implies a general strided memory layout, and StridedArray is defined as you
described. But then a StridedArray does not imply a strided memory layout
because this is valid:

julia> a = rand(4,4)4x4 Array{Float64,2}:
0.706608 0.0505246 0.951182 0.678708
0.0962837 0.321637 0.224594 0.887127
0.406915 0.061083 0.554764 0.170252
0.130242 0.943475 0.438224 0.647376

julia> v = [1, 3, 4]3-element Array{Int64,1}:
1
3
4

julia> sv = sub(a, v, 2:3)3x2 SubArray{Float64,2,Array{Float64,2},Tuple{Array{Int64,1},UnitRange{Int64}},0}:
0.0505246 0.951182
0.061083 0.554764
0.943475 0.438224

Interesting, further experimentation shows:

julia> stride(sv, 1)
ERROR: ArgumentError: strides valid only for RangeIndex indexing
in strides at subarray.jl:359
in stride at subarray.jl:375

julia> pointer(sv)
ERROR: MethodError: convert has no method matching convert(::Type{Ptr{Float64}}, ::SubArray{Float64,2,Array{Float64,2},Tuple{Array{Int64,1},UnitRange{Int64}},0})
This may have arisen from a call to the constructor Ptr{Float64}(...),
since type constructors fall back to convert methods.

julia> sv2 = sub(a, :, 1)4-element SubArray{Float64,1,Array{Float64,2},Tuple{Colon,Int64},2}:
0.706608
0.0962837
0.406915
0.130242

julia> pointer(sv2)
Ptr{Float64} @0x00007f905759b280

However

julia> v2 = 1:2:41:2:3

julia> sub(a, v2, 2:3)2x2 SubArray{Float64,2,Array{Float64,2},Tuple{StepRange{Int64,Int64},UnitRange{Int64}},1}:
0.0505246 0.951182
0.061083 0.554764

julia> pointer(ans)
Ptr{Float64} @0x00007f905759b2a0

So it seems there is no abstract type that implies contiguous memory
layout, so the next best thing is to use StridedArray and check the stride
values at runtime?

—
You are receiving this because you commented.
Reply to this email directly, view it on GitHub
#361,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAe0qZVxAUnI9CpDiwJkLcn5M85YID66ks5qhfPSgaJpZM4JoDbO
.

[timholy]

The traits I supplied in my previous post should also do the right thing. Fundamentally StridedArray is a broken concept, because it's not extensible: it's defined as a Union in base, so no one can ever add more types to it that aren't defined in Base. In contrast, the traits solution is robust and extensible.

[JaredCrean2]

Notice that sv is not a StridedArray

Ah, that's clever and makes the hierarchy a bit more precise

Fundamentally StridedArray is a broken concept

I was getting that feeling. Definitely looking forward to seeing how traits can revamp the array and linear algebra code.

[mbauman]

Closing as a duplicate of JuliaLang/julia#10889