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

Fix interval indexing with offset axes #90

Closed
wants to merge 7 commits into from
Closed

Fix interval indexing with offset axes #90

wants to merge 7 commits into from

Conversation

mbauman
Copy link
Member

@mbauman mbauman commented Jun 8, 2017

This one is terrifying. We were only testing against axes of the form step:step:last. Requires PainterQubits/Unitful.jl#90 to pass tests.

This one is terrifying. We were only testing against axes of the form `step:step:last`. Requires PainterQubits/Unitful.jl#90.
@timholy
Copy link
Member

timholy commented Jun 9, 2017

Whew, glad you caught this.

Maybe add some tests for potential roundoff error? 0.0:0.1:0.3 is often a good test case since

julia> 3*0.1
0.30000000000000004

Or we can just encourage users to use a bit of padding in the intervals they supply.

@mbauman
Copy link
Member Author

mbauman commented Jun 9, 2017

Yeah, I'm ashamed I didn't notice this any sooner.

That's a very good point about floating point errors. I don't want to automatically add slop in the library, but we should be doing this calculation in double-precision when possible:

julia> A = AxisArray(-10:10, -1.0:0.1:1.0);

julia> A[-.3 .. .3]
1-dimensional AxisArray{Int64,1,...} with axes:
    :row, -0.3:0.1:0.3
And data, a 7-element UnitRange{Int64}:
 -3
 -2
 -1
  0
  1
  2
  3

julia> A[(-.3 .. .3) + [0]]
2-dimensional AxisArray{Int64,2,...} with axes:
    :row_sub, -0.2:0.1:0.2
    :row_rep, [0.0]
And data, a 5×1 Array{Int64,2}:
 -2
 -1
  0
  1
  2

In the first iteration of this fix, I had kept the use of unsafe_searchsorted and constructed a temporary range with step(r):step(r):step(r):

julia> AxisArrays.unsafe_searchsorted(.1:.1:.1, -.3 .. .3)
-3:3

That works, but it's not dimensionally correct. I suppose I could add zero(eltype(r)) to the endpoints, but at that point it seemed simpler to just do the division and multiplication directly. I'm not sure what the best approach here would be.

@mbauman
Copy link
Member Author

mbauman commented Jun 10, 2017

Actually, on second thought -- I think a proxy like step(r):step(r):step(r) is exactly the right thing. It may solve dates, too. The only tricky thing is if it will preserve the identity A[(x..y)] == vec(A[(x..y)+[0]]) for floating point ranges.

@timholy
Copy link
Member

timholy commented Jun 10, 2017

There's an advantage in leveraging range methods, since they exploit the TwicePrecision type to handle roundoff. How about basing it on something like this?

nsteps(x, step) = floor(Int, abs(x / step))
function nsteps{T}(x, step::TwicePrecision{T})
    # this is basically a hack because Base hasn't defined x/step at TwicePrecision resolution
    nf = abs(x / convert(T, step))
    nc = ceil(Int, nf)
    return abs(convert(T, nc*step)) <= abs(x) ? nc : floor(Int, nf)
end

julia> r = -1.0:0.1:1.0
-1.0:0.1:1.0

julia> step(r)
0.1

julia> nsteps(0.3, step(r))
2

julia> nsteps(0.3, r.step)
3

Why this works:

julia> r.step
Base.TwicePrecision{Float64}(0.09999999999999987, 1.3322676295501878e-16)

julia> 3*r.step
Base.TwicePrecision{Float64}(0.3, 1.1102230246251526e-17)

julia> convert(Float64, 3*r.step)
0.3

By calling nsteps for both ends of the interval, you can construct the UnitRange{Int} that offsets the indices.

@mbauman
Copy link
Member Author

mbauman commented Jun 11, 2017

Yes, that was the first thing I thought of myself, but I didn't particularly care to think through the required math to define that division. It's a great workaround — thanks! Unfortunately, it doesn't solve the problem for 0.5. I tried my "phony" range idea, but that too runs into issues with 0.5.

@timholy
Copy link
Member

timholy commented Jun 11, 2017

So, make 0.6 required? Alternatively I can try to backport JuliaLang/julia#18777 via compat, but that's a fair amount of work and might introduce problems.

@timholy
Copy link
Member

timholy commented Jun 23, 2017

I really like the commit to just support more TwicePrecision operations (d080469). OK if I steal it and contribute it to Base (with proper attribution, of course)?

I'd love to see this cross the finish line. If it's necessary, though I'd rather not I'm hereby offering to look into backporting JuliaLang/julia#18777 (or at least the pieces needed to get this to pass).

@timholy
Copy link
Member

timholy commented Aug 6, 2017

Rebased and merged in #102

@timholy timholy closed this Aug 6, 2017
@mbauman mbauman deleted the mb/offset branch August 3, 2018 19:33
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

Successfully merging this pull request may close these issues.

2 participants