Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

getindex with ColVecs #417

Open
willtebbutt opened this issue Dec 19, 2021 · 3 comments
Open

getindex with ColVecs #417

willtebbutt opened this issue Dec 19, 2021 · 3 comments

Comments

@willtebbutt
Copy link
Member

We currently return a view to the underlying data when getindex is called on a ColVecs. See

Base.getindex(D::ColVecs, i) = ColVecs(view(D.X, :, i))

@theogf was the last person to touch this bit of code, but I'm not sure who wrote it initially (I don't think it was me). Perhaps @devmotion ?

My concern is that if you do getindex, rather thanview, my intuition is that I wouldn't get a view of the underlying data back, but a copy in a new location in memory. Any thoughts on whether my intuition is off, or should we change this?

@devmotion
Copy link
Member

devmotion commented Dec 19, 2021

I think it wasn't me (actually I thought you implemented ColVecs and RowVecs initially @willtebbutt 😄 ) but in my opinion it is correct to return views here. I view ColVecs and RowVecs as an optimized implementation of eachcol and eachrow with the AbstractArray interface that allows for more performant operations by using the underlying matrix. This is also how (hopefully) eachcol and eachrow will be implemented at some point: JuliaLang/julia#32310

@willtebbutt
Copy link
Member Author

willtebbutt commented Dec 23, 2021

Ah, interesting. How would you propose that view(::ColVecs, idx...) and getindex(::ColVecs, idx...) differ in their output then (assuming that they should)?

@devmotion
Copy link
Member

I think getindex should satisfy the same properties as for Vectors of views:

julia> A = rand(5, 4);

julia> x = [view(A, :, i) for i in axes(A, 2)];

julia> getindex(x, 1)
5-element view(::Matrix{Float64}, :, 1) with eltype Float64:
 0.832858822923345
 0.1486458121118549
 0.20187209965755426
 0.17211775415001118
 0.8813222986233952

julia> getindex(x, 2:3)
2-element Vector{SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}}:
 [0.75958434210484, 0.28127270898481205, 0.11477999962460317, 0.6093522955819292, 0.30000872694236713]
 [0.0028572719840013194, 0.6307387952125958, 0.25339090743930304, 0.2795101302467987, 0.5211441981223575]

Mainly, it should satisfy

julia> getindex(x, 1) isa eltype(x)
true

julia> getindex(x, 1:2) isa typeof(x)
true

While the first property is satisfied for ColVecs and RowVecs (

Base.getindex(D::ColVecs, i::Int) = view(D.X, :, i)
Base.getindex(D::ColVecs, i::CartesianIndex{1}) = view(D.X, :, i)
and
Base.getindex(D::RowVecs, i::Int) = view(D.X, i, :)
Base.getindex(D::RowVecs, i::CartesianIndex{1}) = view(D.X, i, :)
), the second property is violated by
Base.getindex(D::ColVecs, i) = ColVecs(view(D.X, :, i))
and
Base.getindex(D::RowVecs, i) = RowVecs(view(D.X, i, :))
(the type of the wrapped matrix is changed).

On the other hand, view may return more optimized representations without additional allocations:

julia> view(x, 1)
0-dimensional view(::Vector{SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}}, 1) with eltype SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{
Base.OneTo{Int64}}, Int64}, true}:
[0.832858822923345, 0.1486458121118549, 0.20187209965755426, 0.17211775415001118, 0.8813222986233952]

julia> view(x, 2:3)
2-element view(::Vector{SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}}, 2:3) with eltype SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}:
 [0.75958434210484, 0.28127270898481205, 0.11477999962460317, 0.6093522955819292, 0.30000872694236713]
 [0.0028572719840013194, 0.6307387952125958, 0.25339090743930304, 0.2795101302467987, 0.5211441981223575]

I think it would be reasonable to expect

julia> all(y isa typeof(getindex(x, 1)) for y in view(x, 2:3))
true

However, the exact return type of view is not clearly defined as far as I know. I think a natural definition for non-integer index/indices i would be

Base.view(x::ColVecs, i) = ColVecs(view(x.X, :, i))
Base.view(x::RowVecs, i) = RowVecs(view(x.X, :, i))

I'm not completely sure what to expect from view(::ColVecs, ::Int) though, maybe just the same as getindex(::ColVecs, ::Int) - but then it would not be consistent with the example above in which a 0-dimensional view is returned.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants