Implement indexing by value with atvalue

README.md

@@ -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]
+-1.6820949242530765
+
+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}:
+ -1.68209
+```
+
 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

@@ -93,8 +93,23 @@ 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)
 ```
 
+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]
+-1.6820949242530765
+
+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}:
+ -1.68209
+```
+
 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:

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≈", (@sprintf "%0.2e" 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/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