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.
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.647376julia> v = [1, 3, 4]3-element Array{Int64,1}:
1
3
4julia> 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.438224Interesting, 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:375julia> pointer(sv)
ERROR: MethodError:converthas 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.130242julia> pointer(sv2)
Ptr{Float64} @0x00007f905759b280However
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.554764julia> pointer(ans)
Ptr{Float64} @0x00007f905759b2a0So 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.
Closing as a duplicate of https://github.com/JuliaLang/julia/issues/10889
Most helpful comment
The traits I supplied in my previous post should also do the right thing. Fundamentally
StridedArrayis a broken concept, because it's not extensible: it's defined as aUnionin 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.