Non-uniform interpolation in 1D

[Uwe Fechner]

Hello,

I am trying to port the following function from python to julia:

# -*- coding: utf-8 -*-
from scipy.interpolate import InterpolatedUnivariateSpline
import numpy as np
from pylab import plot

P_NOM = [1.5, 2.2, 3.7, 5.6, 7.5, 11.2, 14.9]
ETA   = [93., 94., 94., 95., 95., 95.5, 95.5]

calc_eta = InterpolatedUnivariateSpline(P_NOM, ETA, k=1)

# plotting code, only for testing
if __name__ == "__main__":
    X = np.linspace(1.5, 14.9, 1024, endpoint=True)
    ETA = []
    for alpha in X:
        eta = calc_eta(alpha)
        ETA.append(eta)
    plot(X, ETA)

The resulting plot is shown at the end of this posting.

How can I port this to Julia?

I am trying to use the package "Interpolations.jl", but I do not see any
example, that shows the interpolation on a non-uniform grid.

For now I need only linear interpolation, but I want to use B-Splines
later.

Any hint appreciated!

Uwe Fechner

[Cedric St-Jean]

Hi Uwe,

Have you tried Grid.jl? I haven't tried it, but this example looks like it might work with a non-uniform grid.

# Let's define a quadratic function in one dimension, and evaluate it on an evenly-spaced grid of 5 points:
c = 2.3  # center
a = 8.1  # quadratic coefficient
o = 1.6  # vertical offset
qfunc = x -> a*(x-c).^2 + o
xg = Float64[1:5]
y = qfunc(xg)
yi = InterpGrid(y, BCnil, InterpQuadratic)

[Cedric St-Jean]

Actually, nevermind, I misread the documentation, Grid is for a regularly-spaced grid of points.

[Uwe Fechner]

Hello,

I don't think, that this works on a non-uniform grid. The array xg is evenly spaced, and it
is NOT passed to the function InterpGrid.

Uwe

[Yichao Yu]

On Sat, Feb 27, 2016 at 9:40 AM, Uwe Fechner <uwe.fec...@gmail.com> wrote:

Hello,

I don't think, that this works on a non-uniform grid. The array xg is evenly spaced, and it
is NOT passed to the function InterpGrid.

I've recently tried Dierckx which support non-uniform interpolation. I only need very basic functions so I don't know if it has all the flexibility you need but it's probably worth a look if you haven't.

[Uwe Fechner]

Thanks. The following code works:

using Dierckx

P_NOM = [1.5, 2.2, 3.7, 5.6, 7.5, 11.2, 14.9]
ETA   = [93., 94., 94., 95., 95., 95.5, 95.5]

calc_eta = Spline1D(P_NOM, ETA, k=1)

println(calc_eta(3.5))

Nevertheless I would be interested how to do that with Interpolations.jl. Sometimes you don't have Fortran available.

Best regards:

Uwe

[Matt Bauman]

Interpolations is very similar, but it currently only supports linear and nearest-neighbor schemes for gridded interpolations:

using Interpolations

itp = interpolate((P_NOM,), ETA, Gridded(Linear())) # You pass the x-values as a tuple, since this generalizes to multi-dimensional coordinates
println(itp[3.5])

x = linspace(1.5, 14.9, 1024)
y = itp[x]

plot(x,y)

[Tim Holy]

FYI this is covered in the README:

https://github.com/tlycken/Interpolations.jl#gridded-interpolation

--Tim

> >>>>> <https://lh3.googleusercontent.com/-8OofwCQWohg/VtGwKR-1BOI/AAAAAAAAAQ
> >>>>> I/UTLksCCMIPo/s1600/LinearInterpolation.png>

[Eric Forgy]

Quick note...

Grid.jl has some support for irregular grids, but is not documented (If I recall). I'm using it with linear interpolation on irregular grids, but don't remember if it can do irregular splines.

[Tomas Lycken]

Higher than linear order interpolation on irregular grids is not supported in either Grid.jl or Interpolations.jl (and since the latter is a replacement for the former, it probably never will be supported in Grid). If you have a good scheme that you'd like to see implemented, do file an issue with some pointers to relevant resources that can be used as a baseline for an implementation.

Interpolations.jl now has a superset of the functionality in Grid.jl, but is (a lot!) faster, so if you're on Julia 0.4 (Interpolations.jl is not supported on 0.3) there's no reason to use Grid.jl.

// T

[Uwe Fechner]

Hello,

thanks for your explanations.

Nevertheless I like the syntax of the package Dierckx much more.
The expression Gridded(Linear()) is very confusing for me. What is gridded,
in this example?

Furthermore the Readme file of Interpolations is missing the information, how
to use it for the given problem.

Best regards:

Uwe

[Tomas Lycken]

Gridded here is the interpolation scheme, as opposed to (implicitly uniform) BSpline interpolation; it’s simply the way Interpolations.jl lets you specify which algorithm you want to use. Compare the following incantations:

itp = interpolate(A, BSpline(Linear()), OnGrid()) # linear b-spline interpolation; x is implicitly 1:length(A)
itp = interpolate((x,), A, Gridded(Linear())) # linear, "gridded" interpolation; x can be irregular
itp = interpolate(A, BSpline(Cubic(Flat())), OnCell()) # cubic b-spline with free boundary condition

That is; you give the data to be interpolated (A and, where applicable, x) as well as one or more arguments specifying the algortihm you want to use (for details on OnGrid/OnCell, see the readme). Gridded is what just we’ve called the family of algorithms that support irregular grids.

This is all documented in the readme, but documentation is not just about putting the information in writing - it’s also about putting it in the correct place, where it seems obvious to look for it. If you have suggestions on how this information can be made easier to find and digest, please file an issue or PR. All feedback is most welcome! :)

// T

[jmarcell...@ufpi.edu.br]

I did not understand how to use the "A". "A" is a vector tuple ...?

[Tim Holy]

A is the array of values you want to interpolate. x are the positions.
Effectively, `A = [f(xx) for xx in x]` assuming you want to interpolate the
function `f`.

--Tim

On Wednesday, March 09, 2016 02:25:24 PM jmarcell...@ufpi.edu.br wrote:
> I did not understand how to use the "A". "A" is a vector tuple ...?
>
> Em domingo, 28 de fevereiro de 2016 18:32:23 UTC-3, Tomas Lycken escreveu:

> > Gridded here is the *interpolation scheme*, as opposed to (implicitly

> >>>>>>>> <https://lh3.googleusercontent.com/-8OofwCQWohg/VtGwKR-1BOI/AAAAAAA
> >>>>>>>> AAQI/UTLksCCMIPo/s1600/LinearInterpolation.png>>>>>>
> >>>>> ​

[jmarcell...@ufpi.edu.br]

Hello Tim, good morning
In this case, where we have discrete points (P_NOM = [1.5, 2.2, 3.7, 5.6, 7.5, 11.2, 14.9] and ETA = [93, 94, 94, 95, 95, 95.5, 95.5 ]), as would be "A" for P_NOM(2.8) value ? Since function f (x) is more intuitive.

[Tim Holy]

I don't know what P_NOM and ETA "mean," so I can't answer your question.
Bottom line: if you have a function f, and x is a vector of locations, then

f_itp = interpolate((x,), [f(xx) for xx in x], Gridded(Linear()))

will construct an object so that f_itp[xx] interpolates f at xx.

--Tim