Julia: SubArray Supertype

Created on 19 Aug 2016  Â·  11Comments  Â·  Source: JuliaLang/julia

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.

linear algebra

Most helpful comment

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.

All 11 comments

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.

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

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?

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.

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.

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

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?

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 [email protected] 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
https://github.com/JuliaLang/julia/issues/18131#issuecomment-241095192,
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAe0qZVxAUnI9CpDiwJkLcn5M85YID66ks5qhfPSgaJpZM4JoDbO
.

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.

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.

Was this page helpful?
0 / 5 - 0 ratings

Related issues

tkoolen picture tkoolen  Â·  3Comments

omus picture omus  Â·  3Comments

wilburtownsend picture wilburtownsend  Â·  3Comments

TotalVerb picture TotalVerb  Â·  3Comments

arshpreetsingh picture arshpreetsingh  Â·  3Comments