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

Implement indexing by value with atvalue #52

Merged
merged 10 commits into from
Jul 6, 2017
17 changes: 16 additions & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
17 changes: 16 additions & 1 deletion docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
2 changes: 1 addition & 1 deletion src/AxisArrays.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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, ..
Expand Down
24 changes: 22 additions & 2 deletions src/indexing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,13 @@ const Idx = Union{Real,Colon,AbstractArray{Int}}

using Base: ViewIndex, @propagate_inbounds, tail

immutable Value{T}
val::T
tol::T
end
Value(v,t) = Value(promote(v,t)...)
atvalue(x; rtol=Base.rtoldefault(typeof(x)), atol=zero(x)) = Value(x, atol+rtol*abs(x))

# Defer IndexStyle to the wrapped array
@compat Base.IndexStyle{T,N,D,Ax}(::Type{AxisArray{T,N,D,Ax}}) = IndexStyle(D)

Expand Down Expand Up @@ -170,14 +177,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)
return 0
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I can defeat this with the following:

julia> using AxisArrays, OffsetArrays

julia> A = AxisArray(OffsetArray([1 2; 3 4], 0:1, 1:2), Axis{:x}([1.0,4.0]), Axis{:y}([2.0,6.1]))
2-dimensional AxisArray{Int64,2,...} with axes:
    :x, [1.0, 4.0]
    :y, [2.0, 6.1]
And data, a OffsetArrays.OffsetArray{Int64,2,Array{Int64,2}} with indices 0:1×1:2:
 1  2
 3  4

julia> A[atvalue(5.0)]
1-dimensional AxisArray{Int64,1,...} with axes:
    :y, [2.0, 6.1]
And data, a OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}} with indices 1:2:
 1
 2

Perhaps you can kill two birds with one stone (see your comment below about throwing a more informative BoundsError) and throw the error from here?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I fixed the error you identified and added the ability to throw a BoundsError with a value index. The issue isn't entirely resolved though, as indicated by a @test_broken in the tests.

I've found a few parts of AxisArrays that I haven't touched that don't work as expected when using OffsetArrays, so I'm not sure how much I should attempt to resolve in this PR. Consider the following:

julia> using AxisArrays,OffsetArrays

julia> B = AxisArray(OffsetArray([1,2,3], 2:4), Axis{:junk}([:a, :b, :c]))
1-dimensional AxisArray{Int64,1,...} with axes:
    :junk, Symbol[:a, :b, :c]
And data, a OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}} with indices 2:4:
 1
 2
 3

julia> B[:a]
ERROR: BoundsError: attempt to access OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}} with indices 2:4 at index [1]
Stacktrace:
 [1] throw_boundserror(::OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}}, ::Tuple{Int64}) at ./abstractarray.jl:433
 [2] checkbounds at ./abstractarray.jl:362 [inlined]
 [3] getindex at /Users/ajkeller/.julia/v0.6/OffsetArrays/src/OffsetArrays.jl:94 [inlined]
 [4] getindex at /Users/ajkeller/.julia/v0.6/AxisArrays/src/indexing.jl:26 [inlined]
 [5] getindex(::AxisArrays.AxisArray{Int64,1,OffsetArrays.OffsetArray{Int64,1,Array{Int64,1}},Tuple{AxisArrays.Axis{:junk,Array{Symbol,1}}}}, ::Symbol) at /Users/ajkeller/.julia/v0.6/AxisArrays/src/indexing.jl:111

Perhaps there should be a separate PR that tries to fix these issues, perhaps via translating from standard indices to offset indices using Base.indices? I'm not terribly familiar with the non-standard indexing schemes in Julia so perhaps there's a better way to go.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, don't worry too much about that case yet, I suspect #90 is relevant.

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)
Expand Down
22 changes: 22 additions & 0 deletions test/indexing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -165,3 +165,25 @@ 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)] # more suitable error?
@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 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]
@test_throws ArgumentError A[:,6.1]