BCnearest bug

[ghost]

Here's an minimal example. Sorry for the format:

1d linear interpolation + BCnearest extrapolation. Values higher than the grid look fine but values below the bottom value of the grid are all wrong.

Any suggestions?

Edit: Version 0.3.2 (2014-10-21 22:05 UTC)

[timholy]

Sorry for the slow response, and of course for the bug. I'm traveling now and for the next week, and am unlikely to have time to debug this myself anytime soon. It presumably comes down to a couple of functions in src/interp.jl. Do you know about @which? That should make it more straightforward to find the call chain.

[ghost]

Here's a self-contained example and some @which output. I hope this might help. I have a hard time reading/understanding Julia polymorphism so I'm not sure which of the methods gets called next. Sorry!

julia> using Grid

julia> x = -1.0:0.2:1.0
-1.0:0.2:1.0

julia> y = vec([ 1.0 2.0 3.0 3.0 3.0 3.0 5.0 7.0 5.0 5.0 6.0 ])
11-element Array{Float64,1}:
 1.0
 2.0
 3.0
 3.0
 3.0
 3.0
 5.0
 7.0
 5.0
 5.0
 6.0

julia> ff = CoordInterpGrid(x, y, BCnearest, InterpLinear)
11-element CoordInterpGrid{Float64,1,BCnearest,InterpLinear,(FloatRange{Float64},)}:
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef
 #undef

julia> @which CoordInterpGrid(x, y, BCnearest, InterpLinear)
CoordInterpGrid{R<:Range{T},T<:FloatingPoint (coord::R<:Range{T},A::Array{T<:FloatingPoint,1},args...) at /home/gmihalac/.julia/v0.3/Grid/src/coord.jl:24

julia> ff[-2.3]
0.5

julia> @which ff[-2.3]
getindex{T}(C::CoordInterpGrid{T,1,BC<:BoundaryCondition,IT<:InterpType,R},x::Real) at /home/gmihalac/.julia/v0.3/Grid/src/coord.jl:36

julia> ff[ -3.0:0.1:-2.0  ]
11-element Array{Float64,1}:
 1.0
 0.5
 1.0
 0.5
 1.0
 0.5
 1.0
 0.5
 1.0
 0.5
 1.0

julia> @which ff[ -3.0:0.1:-2.0  ]
getindex{T,R<:Number}(G::AbstractInterpGrid{T,1,BC<:BoundaryCondition,IT<:InterpType},x::AbstractArray{R<:Number,1}) at /home/gmihalac/.julia/v0.3/Grid/src/interp.jl:274

[ghost]

Whatever the problem is, it seems to happen exclusively when extrapolating to the left of the grid. I also tried a case where both bounds of the grid were positive and the problem persists.

[timholy]

What I do is create variables with the same names as the arguments of the functions holding the inputs I want to test, then start copy/pasting the code into the REPL. That way you're effectively executing the function line-by-line.

Of course, when Keno's debugger gets merged, this will get easier.

[tomasaschan]

We've seen this before (see plots in #31) but I don't think we really investigated it though (at least the discussion in that issue doesn't indicate we did). It's likely that there's some subtle bug in the implementation of the wrap function (or one of the ones used with it), but I think a main reason we didn't invest any time in this was that we discovered it right around the time we started focusing on Interpolations.jl instead.

Interpolations.jl is, unfortunately, 0.4-only, and it also doesn't have anything corresponding to CoordInterpGrid yet (although that is trivial to implement, especially if you don't need derivatives), so it's understandable if you don't want to use that just yet. However, there are solutions to these problems in the pipeline.

[ghost]

@tlycken Thank you for the reply! I don't know whether the cluster people here would be open to installing a 0.4 snapshop, otherwise I would gladly try it. At this point I'm contemplating a C or Fortran rewrite (horrible, I know) since I suspect this extrapolation behavior is messing with convergence.

[tomasaschan]

Instead of a complete rewrite in a different language, you could quite easily wrap InterpGrid in a type that works for your use case (constant extrapolation):

type InterpGridBugfix{T<:FloatingPoint, N, BC<:BoundaryCondition, IT<:InterpType} <: AbstractInterpGrid{T, N}
    g::InterpGrid{T, N, BC, IT}
end

import Base.getindex
function getindex(itp::InterpGridBugfix, xs...)
    # these lines can probably be done more efficiently, but this will
    # hopefully clarify the idea enough for you to optimize if you need to
    clamped = similar(xs)
    for i in 1:length(xs)
        clamped[i] = clamp(xs[i], 1, size(itp.g, i))
    end
    return itp.g[clamped...]
end

The idea is simply that you wrap the indexing operation with a method that makes sure you never actually extrapolate, which lets you avoid the bug in Grid.jl. If you, for example, only need 1D interpolation objects, you can make it even simpler:

Base.getindex(itp::InterpGridBugfix, x) = itp[clamp(x, 1, size(itp.g,1))

Constructing one of these objects is really simple - you just pass it the InterpGrid instance. However, you might need to use an explicit constructor of CoordInterpGrid to be able to work it that way.

[agbristol]

I've also noticed this bug. Although it seems from the discussion that effort is being focused on Interpolations.jl and 0.4, I thought I'd record things here.

I'm using BCfill, and I see the same problem. An additional detail I've noticed is that if the negative(/left-of-the-grid) values are expressible as integers, the bug disappears and the call to getindex() behaves as expected. Perhaps that might shed some light?

[ghost]

@agbristol Thank you! I'm still holding out hope for some fix here, since 0.4 and its libraries are far from "production ready" and it's doubtful I can get the cluster people here to install it for me.

[tomasaschan]

@agbristol For BCfill, it's possible to construct a similar workaround as what I proposed for BCnearest, so you could implement that yourself if you're in a hurry. The basic principle is that you wrap the InterpGrid in your own type, define a new indexing method that handles the extrapolation for you, and only intex into the InterpGrid with indices that are inside the domain.

However, the info you provided about indexing working for integer indices is valuable. I'll see if I can find the root cause and fix it, rather than force verbose workarounds on you :)

[mbauman]

I was looking into this a little bit over the weekend. I was getting close to a fix, but I didn't have much time. To save you some digging, there's one issue here at interp.jl#611 (needs to also check ix==1 && ifx < 1 - I think iswrap would be more accurately named isoob (out-of-bounds)). But that doesn't entirely fix it - there is still a bug in how this information is used in set_position (interp.jl#408). I may have a chance to finish this up tonight if you don't beat me to it.

[tomasaschan]

@mbauman Thanks a whole bunch! It's dinner time in Sweden, so you have an hour or two on me. What do you say, shall we race to a PR? ;)

[ghost]

@mbauman @tlycken any hope for this? Thank you for all your time and effort!

[mbauman]

I don't really have time to work on this, but I've pushed my incomplete WIP to a branch on my fork. It at least points a direction on where some of the problems are. Anyone is welcome to move the branch to this repo, add commits, etc.

[tomasaschan]

See #57 for at least a partial fix.

[ghost]

Will Pkg.update() update with the fix? I can't seem to get the package to update. Or do we need timholy to do something on the server-side? Thank you very much for fixing this!

[tomasaschan]

@gaabriel #57 hasn't been merged, so the fix isn't actually pushed or released yet. I don't think there are any big hurdles to clear to get it done, though, so I hope to have time for pushing these today.

[tomasaschan]

@gaabriel I just merged the other PRs and tagged a new version, so as of now, Pkg.update() should give you these fixes.

[ghost]

Great! Thank you very much for all the work, all! Now I just hope this helps with my convergence issues.

[timholy]

@tlycken, thanks a bunch. Can you also do git push --tags?

[tomasaschan]

@timholy Ah, sorry - Pkg.publish() failed, so I just pushed the changes to METADATA manually, but forgot to do that. Should be fixed now.

[timholy]

Thanks!