Output format about contour.jl

[tomasaschan]

Actually, it might be better to introduce two types:

type ContourVertex{T<:Number}
    x::T
    y::T
end

type ContourLine{T<:Number}
    level::Number
    vertices::Vector{ContourVertex{T}}
end

We might even want to parametrize the type of level as well, to make sure we store an unboxed value (I need to read up some more on Julia's type system to be sure I know what I'm talking about here...).

It might also make sense to define some helper methods for these types, e.g.

import Base.length
length{T<:Number}(c::ContourLine{T}) = length(c.vertices)

[tomasaschan]

We could also define, then,

xs{T}(c::ContourLine{T}) = [v.x for v in c.vertices]
ys{T}(c::ContourLine{T}) = [v.y for y in c.vertices]

to enable e.g. Gadfly to plot them using

plot(x=xs(c), y=ys(c), Geom.line)

[tomasaschan]

I asked a question about this kind of stuff on the dev list and it seems there are already a few versions of this kind of types, and there also seems to be consensus that it's probably a good idea to develop a standalone package GeometryBase or something that has types for things like points, lines and line segments. If/when that is done, I think we should use that package, to make sure it's as easy as possible to re-use the output from Contour.jl in other packages.

[tomasaschan]

...in the meantime, maybe ImmutableArrays.jl is a good choice of point representation?

[darwindarak]

Do you mean something like this?

typealias Vertex Vector2{Float64}

type ContourLine{T<:Number}
    level::Number
    vertices::Vector{Vertex}
end

[tomasaschan]

Yeah, that sounds good. I can't really make up my mind, however, on exactly how a set of contours over a field is best represented. I can think of at least two ways to do it - either by a vector of the ContourLines using the type you just defined, or by something like

typealias Vertex Vector2{Float64}

type ContourLevel
    level::Float64
    lines::Vector{Curve2}
end

type Curve2{T}
    vertices::Vector{Vector2{T}}
end

contours would in this case return a Vector{ContourLevel}. The reason to do this distinction is that Curve2{T} is probably generic enough to be included in a GeometryBase package, so then we could easily make the switch in the future. What do you think?

[darwindarak]

I think the ContourLevel approach makes more sense. It gives us the option of changing how we want to store the vertices without significantly changing the output from contours. I'll make the change tonight.

[darwindarak]

I'm running into a lot or problems trying to get the parametric types to work. On Julia 0.3, the line contours = ContourLevel(h, []) seems to work just fine. But when testing on Julia 0.2, the message
ERROR: no method ContourLevel(Float64,Array{None,1})
pops up unless I do something like contours = ContourLevel(h, Array(Curve2,0))

I am not familiar enough with Julia's parametric types to make a clean fix. Do you mind taking a look at it? In particular, lines L127,129, and 177 in 72733ac.

[darwindarak]

I was able to resolve the type errors by forcing z::Matrix{T<:FloatingPoint} in the arguments list of trace_contour. That might just be taking the easy way out, but I don't think it is too restrictive. What do you think?

[tomasaschan]

I think I was able to fix this by adding a couple of constructors to the Curve2 and ContourLine types, as well as assuming that regardless of eltype(z), the contour curves will always be Curve{Float64}. This is somewhat restrictive - for example, it doesn't allow for arbitrary precision floats to be used - but I think it is quite easy to extend the API to fix this if we need to. Anyway, I think it's more important to allow e.g. eltype(z) <: Integer, which we do now.

[tomasaschan]

It seems we have a decision for now, so I'm closing this for housekeeping.

[tomasaschan]

Actually, I just found something that might not have been obvious, but that needs fixing: the relative scaling of the vertices' positions are correct, but the actual values are not. I didn't realize this when I looked at the tests, but the problem is really obvious when you think about it - the current output format requires the user to manually keep track of interpolating the vertices back to the original xy-grid (just as the test script does).

It's also obvious that this will happen if you look at the entry method contour:

function contour(x, y, z, level::Number)
    trace_contour(z,level,get_level_cells(z,level)) # x and y are never used!
end

I tried to fix this by keeping track of the actual x and y grid along the way, but I must have gotten lost in the algorithm somewhere because my output was garbage. But I can see two ways forward here, and I don't know which one of them is better:

1. Keep track of the user-supplied x- and y-grids and make sure to add vertices to the correct places on these grids (i.e. do what I tried but failed to do), or
2. Make Grid.jl a dependency not only for testing, and use something similar as in the test script to interpolate the vertices back to the original xy-grid before returning. In this case, we could maybe leverage Grid.jl also for the interpolation we're currently doing.

@darwindarak, what do you think?

[darwindarak]

Oh wow! Completely forgot about that!

I think rather than interpolating twice (once to find the position relative to the indices, then using Grid.jl to map the indices to the axes), it might be better to just do the entire interpolation in interpolate. I'll add a PR to implement the first approach, and we can hash out the details there? Of course, it would also be great if we can leverage Grid.jl as well, I just haven't quite figured out a clean way to do so.

On a not so related note, I've been attempting to add a contour statistic to Gadfly.jl and our output format makes it slightly cumbersome. Ultimately, I think the output of Stat.contour would include the vectors, x, y, and color to be used with Geom.line. So far, I've been doing something like

for curve in contourlevel.lines
    append!(x,  [vertex.x for vertex in curve.vertices])
    append!(y,  [vertex.y for vertex in curve.vertices])
end

That seems to undo our efforts to pack the coordinates into vertices. What do you think?

[tomasaschan]

I've thought about that too - I'm in a similar pickle, but because I'm using the contours as vertices in an algorithm depending on polygons.

I tried (once, but never again), to inverse-zip the vertex list using xs,ys = zip(vertices...) but it's excruciatingly slow... But it's basically an inverse of Base.zip that we want here.

Eventually, we're hoping that the geometry aspects of this will be abstracted away in a separate library, meaning the problem "here's a curve in the plane - give me the coordinates of the vertices in two lists" will be solved independently of Contour.jl. In the meantime, we could provide methods like

xlist(c::Curve2) = [x for (x,y) in c.vertices]
ylist(c::Curve2) = [y for (x,y) in c.vertices]

in the module, but I don't really like the solution - for one thing, it requires iterating over the vertices twice. Maybe it would be worth defining a single method that does it in one pass?

function coordinates(c::Curve2)
    N = length(c.vertices)
    xlist = Array(Float64,N)
    ylist = Array(Float64,N)

    for (i,v) in enumerate(c.vertices)
        xlist[i] = v[1]
        ylist[i] = v[2]
    end
    xlist, ylist
end

to be used as

xs, ys = coordinates(curve)
append!(x,xs)
append!(y,ys)

in Stat.contour.

[darwindarak]

Let's take the coordinates route. I think we'll get a better idea on how to modify this when/if more people use this package. I can add in the changes later today.

[tomasaschan]

Now that #7 is merged, should we close this? It seems the decisions on output format is more or less done; adding functions to make accessing stuff on these types can be seen as new features.

[darwindarak]

Sounds good!

[tomasaschan]

Nice! I'm going to take a stab at making JuliaLang/julia#7393 merge-ready - hopefully that lends itself to having a clean solution to #6 (via JuliaCI/Coverage.jl#9). I think that's the only thing left for #3 now.

(Moar crosslinkz!)