JuliaMath / Interpolations.jl

Fast, continuous interpolation of discrete datasets in Julia
http://juliamath.github.io/Interpolations.jl/
Other
534 stars 110 forks source link

Simple, concrete example in docs #403

Open OliverEvans96 opened 3 years ago

OliverEvans96 commented 3 years ago

Hello,

There's a whole page of benchmarks front-and-center in the README, but I've been searching for 10 minutes now to figure out how to do a 1D linear interpolation of a set of points similar to scipy's interp1d.

This page discusses (but doesn't give an example of) the first argument to interpolate but not the required options.

This page discusses the options but not the data (until the section on Parametric Splines, which is clear but not what I'm trying to do).

Also, The tests are abstracted such that they do not contain a simple, concrete example either.

The equivalent scipy documentation is abundantly clear and easy to find.

from scipy.interpolate import interp1d

x = np.linspace(0, 10, num=11, endpoint=True)
y = np.cos(-x**2/9.0)
f = interp1d(x, y)
f2 = interp1d(x, y, kind='cubic')
xnew = np.linspace(0, 10, num=41, endpoint=True)

import matplotlib.pyplot as plt
plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')
plt.legend(['data', 'linear', 'cubic'], loc='best')
plt.show()

plot

Sorry to complain, but hopefully it's helpful complaining :) Oliver

OliverEvans96 commented 3 years ago

The first Dierckx README is also great (but I'd rather avoid a fortran wrapper).

Example Usage

using Dierckx

Fit a 1-d spline to some input data (points can be unevenly spaced):

x = [0., 1., 2., 3., 4.]
y = [-1., 0., 7., 26., 63.]  # x.^3 - 1.
spl = Spline1D(x, y)

Evaluate the spline at some new points:

spl([1.5, 2.5])  # result = [2.375, 14.625]
spl(1.5)  # result = 2.375
OliverEvans96 commented 3 years ago

This post provides a good but outdated example:


using Interpolations

function interp1(xpt, ypt, x; method="linear", extrapvalue=nothing)

    if extrapvalue == nothing
        y = zeros(x)
        idx = trues(x)
    else
        y = extrapvalue*ones(x)
        idx = (x .>= xpt[1]) .& (x .<= xpt[end])
    end

    if method == "linear"
        intf = interpolate((xpt,), ypt, Gridded(Linear()))
        y[idx] = intf[x[idx]]

    elseif method == "cubic"
        itp = interpolate(ypt, BSpline(Cubic(Natural())), OnGrid())
        intf = scale(itp, xpt)
        y[idx] = [intf[xi] for xi in x[idx]]
    end

    return y
end

x = 0:pi/4:2*pi
v = sin.(x)
xq = 0:pi/16:1.5*2*pi
vq = interp1(x, v, xq, method="cubic")
mkitti commented 3 years ago

Something like this maybe? http://juliamath.github.io/Interpolations.jl/latest/convenience-construction/#Convenience-notation-1

f(x) = log(x)
xs = 1:0.2:5
A = [f(x) for x in xs]

# linear interpolation
interp_linear = LinearInterpolation(xs, A)
interp_linear(3) # exactly log(3)
interp_linear(3.1) # approximately log(3.1)

# cubic spline interpolation
interp_cubic = CubicSplineInterpolation(xs, A)
interp_cubic(3) # exactly log(3)
interp_cubic(3.1) # approximately log(3.1)
OliverEvans96 commented 3 years ago

Yes, exactly like that :facepalm: Thank you very much. Apparently "Convenience Constructors" didn't pop out to me, although I now see that it's also linked at the bottom of the General Usage page.

vini-fda commented 3 years ago

I'd also love if there was a plotting example. It took hours for me to discover that if I create an interpolation object itp, I can't just treat it as a function input to plots, e.g. plot(itp, x) like in the basic examples of Plots.jl. So, I created an example(directly mirroring the scipy docs mentioned by @OliverEvans96), which I think helps with visual learners like me:

using Interpolations, Plots

# Lower and higher bound of interval
a = 1.0
b = 10.0
# Interval definition
x = a:1.0:b
# Function application by broadcasting
y = @. cos(x^2 / 9.0)
# Interpolations
itp_linear = LinearInterpolation(x, y)
itp_cubic = CubicSplineInterpolation(x, y)
# Interpolation functions
f_linear(x) = itp_linear(x)
f_cubic(x) = itp_cubic(x)
# Plots
width, height = 1500, 800 # not strictly necessary
x_new = a:0.1:b # smoother interval, necessary for cubic spline

scatter(x, y, markersize=10,label="Data points")
plot!(f_linear, x_new, w=3,label="Linear interpolation")
plot!(f_cubic, x_new, linestyle=:dash, w=3, label="Cubic Spline interpolation")
plot!(size = (width, height))
plot!(legend = :bottomleft)

And the generated plot is: image

mkitti commented 3 years ago

Could you submit it as a pull request?

timholy commented 3 years ago

Arguably Plots should support plot(itp, x). What's broken about it? Maybe they just put ::Function in places where it might be better not to constrain the type.

vini-fda commented 3 years ago

Done! @mkitti this is my first PR here, hope I haven't done anything wonky, :sweat_smile:. Feel free to modify as you please!