Re: [julia-users] irregular indexing/slicing for subarrays?

[Tim Holy]

Subarrays must be evenly-spaced. It's part of the requirements for them to be
"StridedArrays" which can be passed to Lapack. The following works:

julia> x=reshape([1:100],10,10);

julia> y = sub(x, 1:3:9, 1:2:8)
3x4 SubArray of 10x10 Int64 Array:
1 21 41 61
4 24 44 64
7 27 47 67



If you wanted a "view" using arbitrary AbstractVector slicing, that would need
to be a different type---the ability to pass to Lapack is far too important to
sacrifice.

--Tim

[dr]

Ah, I understand. What I really need is the ability to arbitrarily slice subarrays to minimize memory usage.

I am not currently concerned with Lapack compatibility, although I agree it should not be sacrificed.

A "view" type would be nice ... and looking at the subarray code, doesn't seem like it would be too difficult.

[Tom Short]

In the DataFrames package, we also have a need for a more general SubArray (or at least SubVectors). We've been waiting to see if it would be implemented in Base before rolling our own. There's a Julia issue on this:

https://github.com/JuliaLang/julia/issues/2193

[Viral Shah]

The original subArray implementation did support this, but later on was constrained for the ability to pass to LAPACK.

I believe the underlying implementation is still general, and we should expose the more general "View" interface. This goes back to the redesign of the Array type hierarchy.

-viral

[Tim Holy]

A bit of archaeology here...if someone wants to add a "view" type that
supports irregular indexing, it seems that

git show 49462addba9555810fce8938ccbedad8a558db05:j/tensor.j

from the main julia directory will get you a version of SubArray that did
this.

Probably the main challenge will be coming up with a decent name for the type.

--Tim

[dr]

Thanks for that pointer. I tried for a bit to get this to work but could not quite get there. My julia knowhow is a bit too basic to understand the problem.

I have attached my attempt at an "ArrayView".

julia> require("ArrayView.jl")

julia> x=reshape([1:100],10,10);

julia> xx=view(x,[1,3,5],[2,5,7])

ERROR: Error showing value of type AbstractArrayView{Int64,2,Array{Int64,2},(Array{Int64,1},Array{Int64,1})}:

MethodError(size,(SYSTEM: show(lasterr) caused an error

ERROR: no method size(AbstractArrayView{Int64,2,Array{Int64,2},(Array{Int64,1},Array{Int64,1})},)

 in length at abstractarray.jl:15

 in show_delim_array at show.jl:61

 in show at show.jl:85

 in show at show.jl:6

 in show at base.jl:83

 in error_show at repl.jl:22

 in error_show at repl.jl:67

 in anonymous at no file:57

 in with_output_color at error.jl:36

 in display_error at client.jl:55

julia> xx[1,1]

ERROR: no method ref(AbstractArrayView{Int64,2,Array{Int64,2},(Array{Int64,1},Array{Int64,1})},Int64,Int64)

[dr]

I hate this google groups interface; it dropped my attachment. Here it is, embedded.

===

typealias Indices{T<:Int} Union(Int, AbstractVector{T})

type AbstractArrayView{T,N,A<:AbstractArray,I<:(AbstractVector...)} <: AbstractArray{T,N}

    parent::A

    indexes::I

    dims::Dims

    AbstractArrayView(p::A, i::I) = new(p, i, map(length, i))

end

view{T,N}(A::AbstractArray{T,N}, i::NTuple{N,AbstractVector}) =

    AbstractArrayView{T,N,typeof(A),typeof(i)}(A, i)

#change integer indexes into Range1 objects                                                                                                                                                              

view(A::AbstractArray, i::Indices...) =

    view(A, ntuple(length(i), j -> isa(i[j],AbstractVector) ? i[j] :

                                                            (i[j]:i[j])))

view(A::AbstractArrayView, i::Indices...) =

    view(A.parent, ntuple(length(i), j -> isa(i[j],AbstractVector) ? A.indexes[j][i[j]] :

                                                                    (A.indexes[j][i[j]]):(A.indexes[j][i[j]])))

size(s::AbstractArrayView) = s.dims

ndims{T,N}(s::AbstractArrayView{T,N}) = N

copy(s::AbstractArrayView) = copy_to(similar(s.parent, size(s)), s)

similar(s::AbstractArrayView, T::Type, dims::Dims) = similar(s.parent, T, dims)

ref(s::AbstractArrayView) = s

ref{T}(s::AbstractArrayView{T,1}, i::Int) = s.parent[s.indexes[1][i]]

ref{T}(s::AbstractArrayView{T,2}, i::Int, j::Int) =

    s.parent[s.indexes[1][i], s.indexes[2][j]]

ref(s::AbstractArrayView, is::Int...) = s.parent[map(ref, s.indexes, is)...]

ref(s::AbstractArrayView, i::Int) = s[ind2sub(size(s), i)...]

[Tim Holy]

I'm guessing it's because you've (properly) made it a type of AbstractArray,
but you haven't provided a size(A, d) method. The generic show() method for
AbstractArray relies on that being available.

I'm guessing you should be able to at least see whether it's been properly
created using dump().

--Tim

[René Hiemstra]

Hi Tim,

I don't really understand the problem with passing the more flexible SubMatrix (or SubArray) to Lapack. Suppose you have something as follows,

typealias Index Union(Int, Vector{Int}, RangeIndex, ...)

type SubMatrix{T,I<:Index, J<:Index}

     data  :: Matrix{T}

     rows  :: I

     cols  :: J

end

A = rand(10,10)

B = SubMatrix(A,[1,3,5], [1,6,7,10])

C = B.data[B.rows, B.cols] 

Then we can always pass C to Lapack, right?

Rene

[Tim Holy]

Old thread alert :-).

SubArrays have changed a lot in 0.4, and now you can create a SubArray with
irregular indexes. You just can't pass them directly to Lapack---only
SubArrays created with Range-type indexes have the uniform strides that Lapack
needs.

Your example exploits the fact that you can always copy the data to a new
array (which has uniform strides) and pass that to Lapack. You could imagine
doing that "automatically" behind the scenes, but then the scary situation
would be that users might inadvertently make such a copy 10^6 times when, had
the user realized it were necessary, it could be done just once. So SubArrays
with irregular indexing will simply throw an error if you try to pass them to
Lapack. To circumvent the error, just call copy on your SubArray and then call
matrix algebra routines.

Best,
--Tim