contours gives “key not found” when n is increased about contour.jl

[darwindarak]

Both Matplotlib and Matlab were able to create contours for that dataset, so it's definitely a bug in Contour.jl. I probably won't have time to dig deeper into this until I finish my prelim exams though (~ a month).

[tomasaschan]

Thanks for the bug report! I started looking into this, and it seems there's some odd edge case stuff going on for these parameters - a larger number of levels sometimes works for me (8, 9, 11, 12, 13 and 16 all work, but 10, 14 and 15 don't. I haven't tested higher).

I'm not sure I'm familiar enough with the algorithm to be able to pin this down, but I'll take a stab at it - and if I can't fix it, I'll leave my debugging findings here for @darwindarak when you have more time on your hands (good luck with the exams!).

[tomasaschan]

The problem is our handling of saddle points. I haven't figured out exactly what's wrong yet, but I have been able to produce a simpler example that reproduces the error:

julia> x = float([1:3]); y = copy(x); z = eye(3,3);
julia> using Contour
julia> contours(x,y,z)
ERROR: key not found: (1,2)
 in chase at /home/tlycken/.julia/v0.4/Contour/src/Contour.jl:148
 in trace_contour at /home/tlycken/.julia/v0.4/Contour/src/Contour.jl:224
 in contours at /home/tlycken/.julia/v0.4/Contour/src/Contour.jl:23
 in contours at /home/tlycken/.julia/v0.4/Contour/src/Contour.jl:28
 in contours at /home/tlycken/.julia/v0.4/Contour/src/Contour.jl:30

I'll see if I can investigate this further.

[tomasaschan]

OK, the problem seems to be these lines, which handle the ambiguous cases. There's an implicit assumption in this code that the two contour lines that cross ambiguous cells are parsed in a specific order (for example, case 16 has two contour lines running from top/left to right/bottom, and assumes that the one from left to bottom is processed first). When this is not the case, the lookup tables here will yield incorrect exit faces and tracing directions, which in turn results in the next call to chase here (or here) updating xi and yi to the wrong cells. In cases where these cells don't have any crossing contour lines, a KeyNotFoundError is thrown the next time cells[xi,yi] is called.

It seems to me that we need to re-think the implementation of these ambiguous cases - fun! ;)

[darwindarak]

That's great that you were able to track down the problem!

What do you think about using something like an adaptive grid method? I haven't had much time to think about it in detail, but what if we use a bilinear interpolation to find the contour crossing whenever we identify an ambiguous cell? We can then split the cell around the crossing point. We then modify the chase function to work recursively through different layers of cell sizes.

[tomasaschan]

An adaptive grid method is probably a nice idea, but I'm not sure it solves the current problem (I haven't given it a lot of thought either...).

A quick fix to solve this would be to keep track of from which edge chase enters a cell, and simplify the correct way. Something like:

simplify(case,inedge)
    # cases that are already simple are returned unchanged
    case <= 15 && return case

    case == 16 &&
        if inedge == up || inedge == left
             return 4
        else
             return 7
        end
   case == 17 && 
        # etc
   # no warranty that the above is actually correct for case 16...
   # if we get to the end, we did something wrong:
   error("cannot simplify case $case")
end

That could probably be made more efficient by making use of a lookup table, but you get the general idea. Now, instead of doing this, we can do something like

simplecase = simplify(case, inedge)
if case != simplecase
    case = simplecase
else
   delete!(cells, (xi,yi))
end

That will ensure that cells always has a correct mapping from cell indices to contour edges we still haven't traced.

[tomasaschan]

OK, I may have cheered too early on this. After more or less re-implementing the entire algorithm (mostly to increase my own understanding) and verifying that the saddle points are handled correctly, I still get a similar KeyNotFoundError for the sample data.

Surprisingly, it only occurs for the unevenly spaced grid; this works fine on my machine:

contours(float([1:length(x)]), float([1:length(y)]), z, 7) # or any other n I've tested

I've tried to use Debug.jl to track down exactly what's going on, but I haven't been able to put my finger on it yet. (I figured I might as well share my findings so far with you anyway...)

@darwindarak Is there something in the algorithm that assumes evenly spaced grids? If so, how should we work around it?

@dcjones I've overwritten my Julia 0.3 instance with 0.4, under which Gadfly doesn't work. Could you check the result of plotting versus index of x and y rather than of the actual values, and get back with the results?

I know Darwin still has much to do elsewhere, so I'll keep digging here, but if any of you has any input it is most welcome.

[darwindarak]

Thanks for looking into this!

If I understand it correctly, I think it only requires Cartesian grids, not necessarily evenly spaced.

[darwindarak]

I think I might have a fix for this. I'll aim for a PR by the end of this weekend.

[darwindarak]

Not quite fixed yet, the problem still shows up using the dataset provided here

https://gist.github.com/dlfivefifty/17f72c2e306f4a78094b