25 days ago
Julia library for 1-d and 2-d splines
This is a Julia wrapper for the dierckx Fortran library, the same library underlying the spline classes in scipy.interpolate. Some of the functionality here overlaps with Interpolations.jl, a pure-Julia interpolation package. This package provides is intended to complement Interpolations.jl and to serve as a benchmark in cases of overlapping functionality.
The Fortran library source code is distributed with the package, so
you need a Fortran compiler on OSX or Linux. On Ubuntu,
sudo apt-get install gfortran will do it.
On Windows, a compiled dll will be downloaded.
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
Equivalent to the above:
evaluate(spl, [1.5, 2.5]) evaluate(spl, 1.5)
Evaluate derivative, integral, or roots:
derivative(spl, 1.5) # derivate at x=1.5; result is 5.75 integrate(spl, 0., 4.) # integrate from x=0 to x=4; result is 60.0 roots(spl) # result is [1.0]
roots() only works for cubic splines (k=3).
Fit a 2-d spline to data on a (possibly irregular) grid:
x = [0.5, 2., 3., 4., 5.5, 8.] y = [0.5, 2., 3., 4.] z = [1. 2. 1. 2.; # size is (length(x), length(y)) 1. 2. 1. 2.; 1. 2. 3. 2.; 1. 2. 2. 2.; 1. 2. 1. 2.; 1. 2. 3. 1.] spline = Spline2D(x, y, z)
Note that if you consider
z as a matrix,
x refers to row
y refers to column coordinates.
Evaluate at element-wise points:
xi = [1., 1.5, 2.3, 4.5, 3.3, 3.2, 3.] yi = [1., 2.3, 5.3, 0.5, 3.3, 1.2, 3.] zi = spline(xi, yi) # 1-d array of length 7 zi = evaluate(spline, xi, yi) # equivalent to previous line
Evaluate at grid spanned by input arrays:
xi = [1., 1.5, 2.3, 4.5] yi = [1., 2.3, 5.3] zi = evalgrid(spline, xi, yi) # 2-d array of size (4, 3)
Spline1D(x, y; w=ones(length(x)), k=3, bc="nearest", s=0.0) Spline1D(x, y, xknots; w=ones(length(x)), k=3, bc="nearest")
k(1 = linear, 2 = quadratic, 3 = cubic, up to 5) from 1-d arrays
y. The parameter
bcspecifies the behavior when evaluating the spline outside the support domain, which is
(minimum(x), maximum(x)). The allowed values are
In the first form, the number and positions of knots are chosen
automatically. The smoothness of the spline is then achieved by
minimalizing the discontinuity jumps of the
kth derivative of the
spline at the knots. The amount of smoothness is determined by the
sum((w[i]*(y[i]-spline(x[i])))**2) <= s, with
given non-negative constant, called the smoothing factor. The number
of knots is increased until the condition is satisfied. By means of
this parameter, the user can control the tradeoff between closeness
of fit and smoothness of fit of the approximation. if
s is too
large, the spline will be too smooth and signal will be lost ; if
s is too small the spline will pick up too much noise. in the
extreme cases the program will return an interpolating spline if
s=0.0 and the weighted least-squares polynomial of degree
s is very large.
In the second form, the knots are supplied by the user. There is no
smoothing parameter in this form. The program simply minimizes the
discontinuity jumps of the
kth derivative of the spline at the
Evalute the 1-d spline
spl at points given in
x, which can be a
1-d array or scalar. If a 1-d array, the values must be monotonically
derivative(spl, x; nu=1)
nu-th derivative of the spline at points in
integrate(spl, a, b)
Definite integral of spline between
For cubic splines (
k=3) only, find roots. Only up to
are returned. A warning is issued if the spline has more roots than
the number returned.
These are the B-spline representation of a curve through N-dimensional space.
ParametricSpline(X; s=0.0, ...) ParametricSpline(u, X; s=0.0, ...) ParametricSpline(X, knots, ...) ParametricSpline(u, X, knots, ...)
X is a 2-d array with size
N is the number of dimensions
of the space (must be between 1 and 10) and
m is the number of points.
X[:, i] gives the coordinates of the
u is a 1-d array giving parameter values at each of the
m points. If not
given, values are calculated automatically.
knots is a 1-d array giving user-specified knots, if desired.
Keyword arguemnts common to all constructor methods:
w: weight applied to each point (length
k: Spline order (between 1 and 5; default 3).
bc: Boundary condition (default
periodic: Treat curve as periodic? (Default is
Spline2D(x, y, z; w=ones(length(x)), kx=3, ky=3, s=0.0) Spline2D(x, y, z; kx=3, ky=3, s=0.0)
ymust be 1-d arrays.
z is also a 1-d array, the inputs are assumed to represent
unstructured data, with
z[i] being the function value at point
(x[i], y[i]). In this case, the lengths of all inputs must match.
z is a 2-d array, the data are assumed to be gridded:
is the function value at
(x[i], y[j]). In this case, it is
size(z) == (length(x), length(y)). (Note that when
z as a matrix,
x gives the row coordinates and
gives the column coordinates.)
evaluate(spl, x, y)
Evalute the 2-d spline
spl at points
(x[i], y[i]). Inputs can be
Vectors or scalars. Points outside the domain of the spline are set to
the values at the boundary.
evalgrid(spl, x, y)
Evaluate the 2-d spline
spl at the grid points spanned by the
y. The input arrays must be
monotonically increasing. The output is a 2-d array with size
output[i, j] is the function value at
(x[i], y[j]). In other words, when interpreting the result as a
x gives the row coordinates and
y gives the column
integral of a 2-d spline over the domain
[x0, x1]*[y0, y1]
integrate(spl, x0, x1, y0, y1)
Spline classes in scipy.interpolate are also thin wrappers
for the Dierckx Fortran library. The performance of Dierckx.jl should
be similar or better than the scipy.interpolate classes. (Better for
small arrays where Python overhead is more significant.) The
equivalent of a specific classes in scipy.interpolate:
|scipy.interpolate class||Dierckx.jl constructor method|
|scipy.interpolate function||Dierckx.jl constructor method|
Dierckx.jl is distributed under a 3-clause BSD license. See LICENSE.md for details. The real*8 version of the Dierckx Fortran library as well as some test cases and error messages are copied from the scipy package, which is distributed under this license.