Implement indexing by value with atvalue

README.md

@@ -13,13 +13,13 @@ Collaboration is welcome! This is still a work-in-progress. See [the roadmap](ht
 ## Example of currently-implemented behavior:
 
 ```julia
-julia> Pkg.clone("https://github.com/JuliaArrays/AxisArrays.jl")
+julia> Pkg.add("AxisArrays")
        using AxisArrays, Unitful
        import Unitful: s, ms, µs
 
 julia> fs = 40000 # Generate a 40kHz noisy signal, with spike-like stuff added for testing
        y = randn(60*fs+1)*3
-       for spk = (sin.(0.8:0.2:8.6) .* [0:0.01:.1; .15:.1:.95; 1:-.05:.05]   .* 50,
+       for spk = (sin.(0.8:0.2:8.6) .* [0:0.01:.1; .15:.1:.95; 1:-.05:.05] .* 50,
                   sin.(0.8:0.4:8.6) .* [0:0.02:.1; .15:.1:1; 1:-.2:.1] .* 50)
            i = rand(round(Int,.001fs):1fs)
            while i+length(spk)-1 < length(y)
@@ -62,7 +62,7 @@ And data, a 2-element Array{Float64,1}:
 julia> A[Axis{:chan}(:c2), Axis{:time}(1:5)]
 1-dimensional AxisArray{Float64,1,...} with axes:
     :time, 0.0 s:2.5e-5 s:0.0001 s
-A[Axis{:chan}(:c2), Axis{:time}(1:5)]:
+And data, a 5-element Array{Float64,1}:
  -6.12181
   0.304668
  15.7366
@@ -92,8 +92,23 @@ julia> axes(ans, 1)
 AxisArrays.Axis{:time,StepRangeLen{Quantity{Float64, Dimensions:{𝐓}, Units:{s}},Base.TwicePrecision{Quantity{Float64, Dimensions:{𝐓}, Units:{s}}},Base.TwicePrecision{Quantity{Float64, Dimensions:{𝐓}, Units:{s}}}}}(5.0e-5 s:2.5e-5 s:0.0002 s)
 ```
 
+You can also index by a single value on an axis using `atvalue`. This will drop
+a dimension. Indexing with an `Interval` type retains dimensions, even
+when the ends of the interval are equal:
+
+```jl
+julia> A[atvalue(2.5e-5s), :c1]
+0.152334
+
+julia> A[2.5e-5s..2.5e-5s, :c1]
+1-dimensional AxisArray{Float64,1,...} with axes:
+    :time, 2.5e-5 s:2.5e-5 s:2.5e-5 s
+And data, a 1-element Array{Float64,1}:
+0.152334
+```
+
 Sometimes, though, what we're really interested in is a window of time about a
-specific index. The operation above (looking for values in the window from 40µs
+specific index. One of the operations above (looking for values in the window from 40µs
 to 220µs) might be more clearly expressed as a symmetrical window about a
 specific index where we know something interesting happened. To represent this,
 we use the `atindex` function:

docs/src/index.md

@@ -1,6 +1,6 @@
 # AxisArrays
 
-[![Build Status](https://travis-ci.org/JuliaArrays/AxisArrays.jl.svg?branch=master)](https://travis-ci.org/JuliaArrays/AxisArrays.jl) [![Coverage Status](https://coveralls.io/repos/JuliaArrays/AxisArrays.jl/badge.svg?branch=master)](https://coveralls.io/r/JuliaArrays/AxisArrays.jl?branch=master)
+[![Build Status](https://travis-ci.org/JuliaArrays/AxisArrays.jl.svg?branch=master)](https://travis-ci.org/JuliaArrays/AxisArrays.jl) [![Coverage Status](https://coveralls.io/repos/github/JuliaArrays/AxisArrays.jl/badge.svg?branch=master)](https://coveralls.io/github/JuliaArrays/AxisArrays.jl?branch=master)
 
 This package for the Julia language provides an array type (the `AxisArray`) that knows about its dimension names and axis values.
 This allows for indexing with the axis name without incurring any runtime overhead.
@@ -13,14 +13,14 @@ Collaboration is welcome! This is still a work-in-progress. See [the roadmap](ht
 ## Example of currently-implemented behavior:
 
 ```julia
-julia> Pkg.clone("https://github.com/JuliaArrays/AxisArrays.jl")
-       using AxisArrays, SIUnits
-       import SIUnits.ShortUnits: s, ms, µs
+julia> Pkg.add("AxisArrays")
+       using AxisArrays, Unitful
+       import Unitful: s, ms, µs
 
 julia> fs = 40000 # Generate a 40kHz noisy signal, with spike-like stuff added for testing
        y = randn(60*fs+1)*3
-       for spk = (sin(0.8:0.2:8.6) .* [0:0.01:.1; .15:.1:.95; 1:-.05:.05]   .* 50,
-                  sin(0.8:0.4:8.6) .* [0:0.02:.1; .15:.1:1; 1:-.2:.1] .* 50)
+       for spk = (sin.(0.8:0.2:8.6) .* [0:0.01:.1; .15:.1:.95; 1:-.05:.05] .* 50,
+                  sin.(0.8:0.4:8.6) .* [0:0.02:.1; .15:.1:1; 1:-.2:.1] .* 50)
            i = rand(round(Int,.001fs):1fs)
            while i+length(spk)-1 < length(y)
                y[i:i+length(spk)-1] += spk
@@ -54,16 +54,15 @@ indices in *any* order, just so long as we annotate them with the axis name:
 
 ```jl
 julia> A[Axis{:time}(4)]
-2-dimensional AxisArray{Float64,2,...} with axes:
-    :time, 7.5e-5 s:2.5e-5 s:7.5e-5 s
-    :chan, [:c1,:c2]
-And data, a 1x2 SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},Colon},2}:
+2-dimensional AxisArray{Float64,1,...} with axes:
+    :chan, Symbol[:c1,:c2]
+And data, a 2-element Array{Float64,1}:
  -1.4144  -2.82879
 
 julia> A[Axis{:chan}(:c2), Axis{:time}(1:5)]
 1-dimensional AxisArray{Float64,1,...} with axes:
     :time, 0.0 s:2.5e-5 s:0.0001 s
-And data, a 5-element SubArray{Float64,1,Array{Float64,2},Tuple{UnitRange{Int64},Int64},2}:
+And data, a 5-element Array{Float64,1}:
  -6.12181
   0.304668
  15.7366
@@ -80,7 +79,7 @@ still has the correct time information for those datapoints!
 julia> A[40µs .. 220µs, :c1]
 1-dimensional AxisArray{Float64,1,...} with axes:
     :time, 5.0e-5 s:2.5e-5 s:0.0002 s
-And data, a 7-element SubArray{Float64,1,Array{Float64,2},Tuple{UnitRange{Int64},Int64},2}:
+And data, a 7-element Array{Float64,1}:
   7.86831
  -1.4144
  -2.02881
@@ -90,11 +89,26 @@ And data, a 7-element SubArray{Float64,1,Array{Float64,2},Tuple{UnitRange{Int64}
  -1.97716
 
 julia> axes(ans, 1)
-AxisArrays.Axis{:time,SIUnits.SIRange{FloatRange{Float64},Float64,0,0,1,0,0,0,0,0,0}}(5.0e-5 s:2.5e-5 s:0.0002 s)
+AxisArrays.Axis{:time,StepRangeLen{Quantity{Float64, Dimensions:{𝐓}, Units:{s}},Base.TwicePrecision{Quantity{Float64, Dimensions:{𝐓}, Units:{s}}},Base.TwicePrecision{Quantity{Float64, Dimensions:{𝐓}, Units:{s}}}}}(5.0e-5 s:2.5e-5 s:0.0002 s)
+```
+
+You can also index by a single value on an axis using `atvalue`. This will drop
+a dimension. Indexing with an `Interval` type retains dimensions, even
+when the ends of the interval are equal:
+
+```jl
+julia> A[atvalue(2.5e-5s), :c1]
+0.152334
+
+julia> A[2.5e-5s..2.5e-5s, :c1]
+1-dimensional AxisArray{Float64,1,...} with axes:
+    :time, 2.5e-5 s:2.5e-5 s:2.5e-5 s
+And data, a 1-element Array{Float64,1}:
+0.152334
 ```
 
 Sometimes, though, what we're really interested in is a window of time about a
-specific index. The operation above (looking for values in the window from 40µs
+specific index. One of the operations above (looking for values in the window from 40µs
 to 220µs) might be more clearly expressed as a symmetrical window about a
 specific index where we know something interesting happened. To represent this,
 we use the `atindex` function:
@@ -125,7 +139,7 @@ julia> idxs = find(diff(A[:,:c1] .< -15) .> 0)
 julia> spks = A[atindex(-200µs .. 800µs, idxs), :c1]
 2-dimensional AxisArray{Float64,2,...} with axes:
     :time_sub, -0.000175 s:2.5e-5 s:0.000775 s
-    :time_rep, SIUnits.SIQuantity{Float64,0,0,1,0,0,0,0,0,0}[0.178725 s,0.806825 s,0.88305 s,1.47485 s,1.50465 s,1.53805 s,1.541025 s,2.16365 s,2.368425 s,2.739 s  …  57.797925 s,57.924075 s,58.06075 s,58.215125 s,58.6403 s,58.96215 s,58.990225 s,59.001325 s,59.48395 s,59.611525 s]
+    :time_rep, Quantity{Float64, Dimensions:{𝐓}, Units:{s}}[0.178725 s,0.806825 s,0.88305 s,1.47485 s,1.50465 s,1.53805 s,1.541025 s,2.16365 s,2.368425 s,2.739 s  …  57.797925 s,57.924075 s,58.06075 s,58.215125 s,58.6403 s,58.96215 s,58.990225 s,59.001325 s,59.48395 s,59.611525 s]
 And data, a 39x242 Array{Float64,2}:
  -1.53038     4.72882     5.8706    …  -0.231564      0.624714   3.44076
  -2.24961     2.12414     5.69936       7.00179       2.30993    5.20432

src/AxisArrays.jl

@@ -6,7 +6,7 @@ using Base: tail
 using RangeArrays, IntervalSets
 using Compat
 
-export AxisArray, Axis, axisnames, axisvalues, axisdim, axes, atindex
+export AxisArray, Axis, axisnames, axisvalues, axisdim, axes, atindex, atvalue
 
 # From IntervalSets:
 export ClosedInterval, ..

src/indexing.jl

@@ -2,6 +2,23 @@ const Idx = Union{Real,Colon,AbstractArray{Int}}
 
 using Base: ViewIndex, @propagate_inbounds, tail
 
+immutable Value{T}
+    val::T
+    tol::T
+end
+Value(x, tol=Base.rtoldefault(typeof(x))*abs(x)) = Value(promote(x,tol)...)
+atvalue(x; rtol=Base.rtoldefault(typeof(x)), atol=zero(x)) = Value(x, atol+rtol*abs(x))
+
+# For throwing a BoundsError with a Value index, we need to define the following
+# (note that we could inherit them for free, were Value <: Number)
+Base.start(::Value) = false
+Base.next(x::Value, state) = (x, true)
+Base.done(x::Value, state) = state
+
+# How to show Value objects (e.g. in a BoundsError)
+Base.show(io::IO, v::Value) =
+    print(io, string("Value(", v.val, ", tol=", v.tol, ")"))
+
 # Defer IndexStyle to the wrapped array
 @compat Base.IndexStyle{T,N,D,Ax}(::Type{AxisArray{T,N,D,Ax}}) = IndexStyle(D)
 
@@ -170,14 +187,27 @@ axisindexes(ax, idx) = axisindexes(axistrait(ax.val), ax.val, idx)
 axisindexes(::Type{Unsupported}, ax, idx) = error("elementwise indexing is not supported for axes of type $(typeof(ax))")
 axisindexes(t, ax, idx) = error("cannot index $(typeof(ax)) with $(typeof(idx)); expected $(eltype(ax)) axis value or Integer index")
 
-# Dimensional axes may be indexed directy by their elements if Non-Real and unique
+# Dimensional axes may be indexed directly by their elements if Non-Real and unique
 # Maybe extend error message to all <: Numbers if Base allows it?
-axisindexes{T<:Real}(::Type{Dimensional}, ax::AbstractVector{T}, idx::T) = error("indexing by axis value is not supported for axes with $(eltype(ax)) elements; use an ClosedInterval instead")
+axisindexes(::Type{Dimensional}, ax::AbstractVector, idx::Real) =
+    throw(ArgumentError("invalid index: $idx. Use `atvalue` when indexing by value."))
 function axisindexes(::Type{Dimensional}, ax::AbstractVector, idx)
     idxs = searchsorted(ax, ClosedInterval(idx,idx))
     length(idxs) > 1 && error("more than one datapoint lies on axis value $idx; use an interval to return all values")
     idxs[1]
 end
+# Dimensional axes may always be indexed by value if in a Value type wrapper.
+function axisindexes(::Type{Dimensional}, ax::AbstractVector, idx::Value)
+    idxs = searchsorted(ax, ClosedInterval(idx.val,idx.val))
+    length(idxs) > 1 && error("more than one datapoint lies on axis value $idx; use an interval to return all values")
+    if length(idxs) == 1
+        idxs[1]
+    else # it's zero
+        last(idxs) > 0 && abs(ax[last(idxs)] - idx.val) < idx.tol && return last(idxs)
+        first(idxs) <= length(ax) && abs(ax[first(idxs)] - idx.val) < idx.tol && return first(idxs)
+        throw(BoundsError(ax, idx))
+    end
+end
 
 # Dimensional axes may be indexed by intervals to select a range
 axisindexes{T}(::Type{Dimensional}, ax::AbstractVector{T}, idx::ClosedInterval) = searchsorted(ax, idx)

test/REQUIRE

@@ -1,2 +1,2 @@
 OffsetArrays
-Unitful 0.1.0
+Unitful 0.2.6

test/indexing.jl

@@ -165,3 +165,38 @@ A = AxisArray(rand(2,2), :x, :y)
 acc = zeros(Int, 4, 1, 2)
 Base.mapreducedim!(x->x>5, +, acc, A3)
 @test acc == reshape([1 3; 2 3; 2 3; 2 3], 4, 1, 2)
+
+# Indexing by value using `atvalue`
+A = AxisArray([1 2; 3 4], Axis{:x}([1.0,4.0]), Axis{:y}([2.0,6.1]))
+@test @inferred(A[atvalue(1.0)]) == @inferred(A[atvalue(1.0), :]) == [1,2]
+# `atvalue` doesn't require same type:
+@test @inferred(A[atvalue(1)]) == @inferred(A[atvalue(1), :]) ==[1,2]
+@test A[atvalue(4.0)] == A[atvalue(4.0),:] == [3,4]
+@test A[atvalue(4)] == A[atvalue(4),:] == [3,4]
+@test_throws BoundsError A[atvalue(5.0)]
+@test @inferred(A[atvalue(1.0), atvalue(2.0)]) == 1
+@test @inferred(A[:, atvalue(2.0)]) == [1,3]
+@test @inferred(A[Axis{:x}(atvalue(4.0))]) == [3,4]
+@test @inferred(A[Axis{:y}(atvalue(6.1))]) == [2,4]
+@test @inferred(A[Axis{:x}(atvalue(4.00000001))]) == [3,4]
+@test @inferred(A[Axis{:x}(atvalue(2.0, atol=5))]) == [1,2]
+@test_throws BoundsError A[Axis{:x}(atvalue(4.00000001, rtol=0))]
+
+# Indexing by value into an OffsetArray
+A = AxisArray(OffsetArrays.OffsetArray([1 2; 3 4], 0:1, 1:2),
+    Axis{:x}([1.0,4.0]), Axis{:y}([2.0,6.1]))
+@test_broken @inferred(A[atvalue(4.0)]) == [3,4]
+@test @inferred(A[:, atvalue(2.0)]) == OffsetArrays.OffsetArray([1,3], 0:1)
+@test_throws BoundsError A[atvalue(5.0)]
+
+# Indexing by value directly is forbidden for indexes that are Real
+@test_throws ArgumentError A[4.0]
+@test_throws ArgumentError A[BigFloat(1.0)]
+@test_throws ArgumentError A[1.0f0]
+if VERSION == v"0.5.2"
+    # Cannot compose @test_broken with @test_throws (Julia #21098)
+    # A[:,6.1] incorrectly throws a BoundsError instead of an ArgumentError on Julia 0.5.2
+    @test_broken A[:,6.1]
+else
+    @test_throws ArgumentError A[:,6.1]
+end

test/readme.jl

@@ -22,7 +22,9 @@ axes(ax, 1)
 A[atindex(-90µs .. 90µs, 5), :c2]
 idxs = find(diff(A[:,:c1] .< -15) .> 0)
 spks = A[atindex(-200µs .. 800µs, idxs), :c1]
-
+A[atvalue(2.5e-5s), :c1]
+A[2.5e-5s..2.5e-5s, :c1]
+A[atvalue(25.0µs)]
 
 # # A possible "dynamic verification" strategy
 # const readmefile = joinpath(dirname(dirname(@__FILE__)), "README.md")