From c18a090445f6b6f60a54250cc5622559f9ccd941 Mon Sep 17 00:00:00 2001 From: timholy Date: Thu, 4 Sep 2014 12:59:06 -0500 Subject: [PATCH] WIP: add views founded on cartesian indexing --- base/sysimg.jl | 3 + base/views.jl | 205 +++++++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 2 +- test/views.jl | 64 +++++++++++++++ 4 files changed, 273 insertions(+), 1 deletion(-) create mode 100644 base/views.jl create mode 100644 test/views.jl diff --git a/base/sysimg.jl b/base/sysimg.jl index c9b6694af79fc..2f1aa5b3d9ad4 100644 --- a/base/sysimg.jl +++ b/base/sysimg.jl @@ -278,6 +278,9 @@ function __init__() fdwatcher_init() end +# Views +include("views.jl") + include("precompile.jl") include = include_from_node1 diff --git a/base/views.jl b/base/views.jl new file mode 100644 index 0000000000000..775dcf7f12a54 --- /dev/null +++ b/base/views.jl @@ -0,0 +1,205 @@ +module Views + +import Base: copy, eltype, getindex, length, ndims, setindex!, similar, size + +export + # types + View, + # functions + sliceview, + subview + +DIMS = 1:6 + +typealias ViewIndex Union(Int, UnitRange{Int}, StepRange{Int,Int}, Vector{Int}) +typealias NonSliceIndex Union(UnitRange{Int}, StepRange{Int,Int}, Vector{Int}) + +# Since there are no multidimensional range objects, we only permit 1d indexes +immutable View{T,N,P<:AbstractArray,I<:(ViewIndex...)} <: AbstractArray{T,N} + parent::P + indexes::I + dims::NTuple{N,Int} +end + +# Simple utilities +eltype{T,N,P,I}(V::View{T,N,P,I}) = T +eltype{T,N,P,I}(::Type{View{T,N,P,I}}) = T +ndims{T,N,P,I}(V::View{T,N,P,I}) = N +ndims{T,N,P,I}(::Type{View{T,N,P,I}}) = N +size(V::View) = V.dims +size(V::View, d::Integer) = d <= ndims(V) ? (@inbounds ret = V.dims[d]; ret) : 1 +length(V::View) = prod(V.dims) +similar(V::View, T, dim::Dims) = similar(V.parent, T, dims) + +# View creation +# We have two modes of creation: subview and sliceview. subview drops trailing singleton +# dimensions, sliceview drops all singleton dimensions. +for N in DIMS + vars = Symbol[symbol(string("I_",k)) for k = 1:N] + # Input arguments to subview and sliceview + createargtypes = Array(Any, N) + createargs = Array(Any, N) + # Expressions for creating the indexes field + indexesslice = Array(Any, N) + indexessub = Array(Any, N) + # Size of each supplied argument, and info about whether to drop it + sizeexprs = Array(Any, N) + keepslice = Array(Bool, N) + keepsub = Array(Bool, N) + for i = 0:2^N-1 # treat scalar and nonscalar, leading to 2^N possibilities + Mslice = 0 # dimensionality of created slice object + Msub = 0 # dimensionality of created sub object + trailing = true # is the current dimension still a trailing singleton? + fill!(keepslice, false) + fill!(keepsub, false) + for k = N:-1:1 # starting from the back lets us update trailing + if (i>>(k-1))&0x01 > 0 + # Non-scalar argument + createargtypes[k] = NonSliceIndex + indexesslice[k] = indexessub[k] = vars[k] + keepslice[k] = keepsub[k] = true + trailing = false + sizeex = :(length($(vars[k]))) + sizeexprs[k] = sizeex + Mslice += 1 + Msub += 1 + else + # Scalar argument + createargtypes[k] = Real + sizeexprs[k] = 1 + keepsub[k] = !trailing + Msub += !trailing + indexesslice[k] = indexessub[k] = :(int($(vars[k]))) + if !trailing + indexessub[k] = :(int($(vars[k])):int($(vars[k]))) + end + end + createargs[k] = :($(vars[k])::$(createargtypes[k])) + end + indexesslice_t = Expr(:tuple, indexesslice...) + indexessub_t = Expr(:tuple, indexessub...) + dimsslice = Expr(:tuple, sizeexprs[keepslice]...) + dimssub = Expr(:tuple, sizeexprs[keepsub]...) + eval(quote + sliceview{T}(A::AbstractArray{T,$N}, $(createargs...)) = + Views.View{T,$Mslice,typeof(A),typeof($indexesslice_t)}(A, $indexesslice_t, $dimsslice) + subview{T}(A::AbstractArray{T,$N}, $(createargs...)) = + Views.View{T,$Msub, typeof(A),typeof($indexessub_t)}(A, $indexessub_t, $dimssub) + end) + end +end + +unsafe_getindex(v::Real, ind::Int) = v +unsafe_getindex(v::Range, ind::Int) = first(v) + (ind-1)*step(v) +unsafe_getindex(v::BitArray, ind::Int) = Base.unsafe_bitgetindex(v.chunks, ind) +unsafe_getindex(v::AbstractArray, ind::Int) = v[ind] +unsafe_getindex(v, ind::Real) = unsafe_getindex(v, to_index(ind)) + + +# TODO: Constructing from another View + + +function index_generate(Nd, Itypes, Vsym, Isyms) + if isempty(Isyms) + Isyms = Any[1] # this handles the syntax getindex(V) + end + exhead = :nothing + if length(Isyms) < Nd + # Linear indexing in the last index + n = Nd - length(Isyms) + m = length(Isyms) + strides = [gensym() for i = 1:n] + indexes = [gensym() for i = 1:n+1] + resid = gensym() + linblock = Array(Expr, 2n+2) + linblock[1] = :($(strides[1]) = size($Vsym, $m)) + for k = 2:n + m += 1 + linblock[k] = :($(strides[k]) = $(strides[k-1]) * size($Vsym, $m)) + end + k = n+1 + linblock[k] = :($resid = $(Isyms[end])-1) + for i = n:-1:1 + k += 1 + linblock[k] = quote + $(indexes[i+1]), $resid = divrem($resid, $(strides[i])) + $(indexes[i+1]) += 1 + end + end + linblock[end] = :($(indexes[1]) = $resid+1) + exhead = Expr(:block, linblock...) + pop!(Isyms) + append!(Isyms, indexes) + end + NP = length(Itypes) + indexexprs = Array(Any, NP) + j = 0 + for i = 1:NP + if Itypes[i] <: Real + indexexprs[i] = :($Vsym.indexes[$i]) + else + j += 1 + indexexprs[i] = :(unsafe_getindex($Vsym.indexes[$i], $(Isyms[j]))) # TODO: make Range bounds-checking respect @inbounds + end + end + # Append any extra indexes. Must be trailing 1s or it will cause a BoundsError. + for k = j+1:length(Isyms) + push!(indexexprs, :($(Isyms[k]))) + end + exhead, :($Vsym.parent[$(indexexprs...)]) +end + +# Scalar indexing +# It's possible to index a 3D object with only 2 indexes, in which case the last index is interpreted as a linear index. +# Likewise, it's possible to index a 2D object with 3 indexes, in which case the last had better be 1 or +# you'll get a bounds error. +# So we have to support DIMS x DIMS possibilities, specializing each on all possible typeof(V.indexes). +# That's a lot of methods. +for Nargs = 0:maximum(DIMS) + args = Symbol[symbol(string("I_",k)) for k = 1:Nargs] + typedargs = map(x->:($x::Real), args) + for Nd in DIMS + Vindextypes = Array(Any, Nd) + for i = 0:2^Nd-1 # treat scalar and nonscalar, leading to 2^Nd possibilities + np = 0 + for k = 1:Nd + if (i>>(k-1))&0x01 > 0 + Vindextypes[k] = NonSliceIndex + np += 1 + s = symbol(string("P", np)) + else + Vindextypes[k] = Int + end + end + body = index_generate(np, Vindextypes, :V, copy(args)) + tupleexpr = Expr(:tuple, Vindextypes...) + eval(quote + function getindex{T,A,TT<:$tupleexpr}(V::View{T,$np,A,TT}, $(typedargs...)) + $(body[1]) + $(body[2]) + end + function setindex!{T,A,TT<:$tupleexpr}(V::View{T,$np,A,TT}, val, $(typedargs...)) + $(body[1]) + $(body[2]) = val + end + end) + end + end +end + + +# # More utility functions +# stagedfunction copy(V::View) +# T, N = eltype(V), ndims(V) +# quote +# A = Array($T, V.dims) +# k = 1 +# Base.Cartesian.@nloops $N i A begin +# @inbounds A[k] = Base.Cartesian.@nref($N, V, i) +# k += 1 +# end +# A +# end +# end + +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 68ebca260097c..a7a8e114c0e9b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,7 +9,7 @@ testnames = [ "floatapprox", "readdlm", "regex", "float16", "combinatorics", "dates", "sysinfo", "rounding", "ranges", "mod2pi", "euler", "show", "lineedit", "replcompletions", "repl", "test", "examples", "goto", - "llvmcall", "grisu" + "llvmcall", "grisu", "views" ] @unix_only push!(testnames, "unicode") diff --git a/test/views.jl b/test/views.jl new file mode 100644 index 0000000000000..3874a1d4e314b --- /dev/null +++ b/test/views.jl @@ -0,0 +1,64 @@ +import Base.Views +using Iterators, Base.Test + +A = reshape(1:75, 5, 5, 3) +# Copy the values to a subarray. Do this manually since getindex(A, ...) will have a new meaning +indexes1 = (2:4, 3, 1:2) +Aslice1 = Array(eltype(A), length(indexes1[1]), length(indexes1[3])) +i = 1 +for p in product(indexes1...) + Aslice1[i] = A[p...] + i += 1 +end +Asub1 = reshape(Aslice1, map(length, indexes1)) +indexes2 = (2:4, 1:2, 3) +Aslice2 = Array(eltype(A), length(indexes2[1]), length(indexes2[2])) +i = 1 +for p in product(indexes2...) + Aslice2[i] = A[p...] + i += 1 +end +Asub2 = Aslice2 # drop the last singleton dimension + +# sliceview +for (indexes,Aslice) in ((indexes1,Aslice1), (indexes2,Aslice2)) + B = Views.sliceview(A, indexes...) + @test ndims(B) == 2 + @test B[1,1] == Aslice[1,1] + @test B[2,2] == Aslice[2,2] + @test B[3] == Aslice[3] + @test B[4] == Aslice[4] + @test Aslice == B + +# C = Views.sliceview(B, 1:3, 1:2) +# @test C == B +# C = Views.sliceview(B, 2:3, 2) +# @test ndims(C) == 1 +# @test C[1] == B[2,2] +# @test C[2] == B[3,2] +# C = Views.sliceview(B, 2, 1:2) +# @test ndims(C) == 1 +# @test C[1] == B[2,1] +# @test C[2] == B[2,2] +# +# @test C[1,1] == B[2,1] +# @test_throws BoundsError C[1,2] +end + +# subview +for (indexes,Asub) in ((indexes1,Asub1), (indexes2,Asub2)) + B = Views.subview(A, indexes...) + @test Asub == B + +# C = Views.subview(B, 1:3, 1, 1:2) +# @test C == B +# C = Views.subview(B, 2:3, 1, 2) +# @test ndims(C) == 1 +# @test C[1] == B[2,2] +# @test C[2] == B[3,2] +# C = Views.subview(B, 2, 1, 1:2) +# @test ndims(C) == 3 +# @test C[1] == B[2,1] +# @test C[2] == B[2,2] +end +