Mismatched coordinate scales

[brendanjohnharris]

Nice package, and easy to use!

I come across some weird behavior, though, when the x and y coordinates of the grid points have different scales.
I wonder if this could be easily fixed by having a normalization step before constructing the interpolant, which would scale the given points to an e.g. unit box (then, when the interpolant is called, the coordinates used for the new interpolated data can be transformed with the same normalization)?

An example of the current problem below, and let me know if I've missed something obvious in the docs that solves this 😅

---

When the y indices have a scale that is 50-times smaller than the x indices:

   using NaturalNeighbours
    using CairoMakie
    N = 50
    f = (xy) -> sinc(sqrt((xy[1])^2 + (xy[2] * N)^2) .+ 0.1 * randn())
    x = Float64.(-6:0.1:6)
    y = Float64.((-6:0.1:6) ./ N)
    xy = Iterators.product(x, y)
    z = f.(xy)

    itp = interpolate(first.(xy)[:], last.(xy)[:], z[:])
    _x = range(-6, 6, length=100)
    _y = range(-6 / N, 6 / N, length=100)
    _xy = Iterators.product(_x, _y)
    _z = itp(first.(_xy)[:], last.(_xy)[:])
    Z = reshape(_z, length.([_x, _y])...)

    f = Figure(size=(600, 300))
    ax = Axis(f[1, 1]; title="Input")
    p = heatmap!(ax, x, y, z)
    ax1 = Axis(f[1, 2]; title="Output")
    heatmap!(ax1, _x, _y, Z)
    display(f)

It works fine for N = 1 though (matched scales):

[DanielVandH]

Interesting. I tried this for some other functions and, it of course is fine for e.g. f(x, y) = x^2 + y^2, but I can find it for others like e.g.

Nx = 1
Ny = 50
f = (x, y) -> tanh(x * Nx + y * Ny - (y * Ny)^2)
xg = collect(-6:0.1:6) ./ Nx
yg = collect(-6:0.1:6) ./ Ny
x = [x for x in xg, y in yg] |> vec
y = [y for x in xg, y in yg] |> vec
z = f.(x, y)
method = Hiyoshi(2)
itp = interpolate(x, y, z, method=method)
xvg = collect(-6:0.04:6.0) ./ Nx
yvg = collect(-6:0.04:6.0) ./ Ny
xv = [x for x in xvg, y in yvg] |> vec
yv = [y for x in xvg, y in yvg] |> vec
zv = itp(xv, yv)
f = Figure(size=(600, 300))
ax = Axis(f[1, 1]; title="Input")
p = heatmap!(ax, x, y, z)
ax1 = Axis(f[1, 2]; title="Output")
heatmap!(ax1, xv, yv, zv)
display(f)

I tried looking online for other such issues but can't really find any discussion of it. I suppose your idea would probably work, Here's an idea:

struct ScaledInterpolant{I}
    itp::I
    xmin::Float64
    xmax::Float64
    ymin::Float64
    ymax::Float64
end
function scaled_interpolant(x, y, z; kwargs...)
    xmin, xmax = extrema(x)
    ymin, ymax = extrema(y)
    x̂ = (x .- xmin) ./ (xmax - xmin)
    ŷ = (y .- ymin) ./ (ymax - ymin)
    itp = interpolate(x̂, ŷ, z; kwargs...)
    return ScaledInterpolant(itp, xmin, xmax, ymin, ymax)
end
function (itp::ScaledInterpolant)(x, y; kwargs...)
    x̂ = (x .- itp.xmin) ./ (itp.xmax - itp.xmin)
    ŷ = (y .- itp.ymin) ./ (itp.ymax - itp.ymin)
    return itp.itp(x̂, ŷ; kwargs...)
end

Nx = 1
Ny = 50
f = (x, y) -> sinc(sqrt((x * Nx)^2 + (y * Ny)^2) .+ 0.1 * randn())
xg = collect(-6:0.1:6) ./ Nx
yg = collect(-6:0.1:6) ./ Ny
x = [x for x in xg, y in yg] |> vec
y = [y for x in xg, y in yg] |> vec
z = f.(x, y)
method = Hiyoshi(2)
itp = interpolate(x, y, z, method=method)
xvg = collect(-6:0.04:6.0) ./ Nx
yvg = collect(-6:0.04:6.0) ./ Ny
xv = [x for x in xvg, y in yvg] |> vec
yv = [y for x in xvg, y in yvg] |> vec
zv = itp(xv, yv)
f = Figure(size=(600, 300))
ax = Axis(f[1, 1]; title="Input")
p = heatmap!(ax, x, y, z)
ax1 = Axis(f[1, 2]; title="Output w/o Scaling")
heatmap!(ax1, xv, yv, zv)
ax2 = Axis(f[1, 3]; title="Output w/ Scaling")
itp = scaled_interpolant(x, y, z, method=method)
zv = itp(xv, yv)
heatmap!(ax2, xv, yv, zv)
display(f)

What do you think?

I won't be able to dedicate the time to work on this (in fact, I'll be putting a notice on this package soon that I unfortunately won't have the time to push new content into this package), but if you are interested I would be happy to assist with any PRs / review. Care would have to be taken for how we differentiate this interpolant, of course.

[DanielVandH]

I also note that, I don't want this to be a step that's done automatically within the interpolator since it could cause more problems / won't be obvious to the user. I reckon either having it as a keyword or as its own struct (the former idea, having it as a keyword, would just redirect into a struct like the one I defined above) is more appropriate.

[DanielVandH]

I tried looking online for other such issues but can't really find any discussion of it

My first guess was that this was a triangulation issue, and actually it is (if you don't know, the interpolant is pretty reliant on the quality of the Delaunay triangulation of the points). I found some other discussion about it: https://blogs.mathworks.com/loren/2019/02/04/data-scaling-for-scattered-interpolation/. The normalisation in this blog post might be more appropriate.

[brendanjohnharris]

Ah, it makes sense that the triangulation isn't normalized, and in general it shouldn't be; it could be desirable for the quality of the meshing could depend on the absolute scale of the in each dimension, or, as in my case, it could not. I guess it's not too much work to normalize the data before interpolating, but that means keeping track of the parameters externally which could be inconvenient when working with the interpolant on the fly.

I'll see if I can put together a small PR from your example above over the weekend. Thanks for your help!

[DanielVandH]

Did you ever have a shot at trying any of this @brendanjohnharris? No problem if not, just curious if the normalisation step gave the results you expected. Could be a nice extra tutorial in the documentation if not a full feature (which would be obviously more time-consuming to implement).

[brendanjohnharris]

In the end I went with pre-normalizing the data, as below:

using NaturalNeighbours, Normalization, TimeseriesTools
function interpolate(X::TimeseriesTools.AbstractTimeSeries; kwargs...)
    xy = [Float64.(x) for x in lookup(X)]
    N = fit.([MinMax], xy) # Could just be (x-minimum(x))/(maximum(x) - minimum(x))
    points = Iterators.product((xy .|> N)...) |> collect |> vec
    itp = interpolate(Float64.(first.(points)), Float64.(last.(points)), X.data[:];
                      kwargs...)
    return itp, N
end

MinMax normalization worked well, but you could also use any of the others (z-score, sigmoid, etc.) from Normalization.jl
Yeah I think having something in the docs to flag when and how you might need to pre-normalize sounds like a good idea

[DanielVandH]

Oh nice, glad that works. I will add that to the documentation soon, thanks! I think some care might be needed for differentiation too which I'll look into.