API documentation

AllDerivatives{N}

A singleton type that represents all m-th derivatives where 0 ≤ m < N.

AllDerivatives(N)

A shortcut for AllDerivatives{N}(), representing all m-th derivatives where 0 ≤ m < N.

BSpline{B} <: Spline{B}

Type for a B-spline from a B-spline basis of type B. BSplines can be obtained by indexing into a B-spline basis.

BSpline(basis::BSplineBasis, index)

Return the index-th B-spline of basis.

Example

julia> basis = BSplineBasis(5, 1:10);

julia> BSpline(basis, 3)
BSpline{BSplineBasis{UnitRange{Int64}}}:
 basis: 13-element BSplineBasis{UnitRange{Int64}}:
  order: 5
  breakpoints: 1:10
 index: 3 (knots: [1, 1, 1, 2, 3, 4])

BSplineBasis{T<:AbstractVector{<:Real}}

Type for a B-spline basis with breakpoint vector of type T.

Here, a B-spline basis is completely specified by its order $k$ and breakpoint sequence. The knot sequence is derived from the breakpoint sequence by duplicating the first and last breakpoints so they each appear $k$ times. Knot sequences where the first and last breakpoints do not appear $k$ times are not supported by this data type.

BSplineBasis(order, breakpoints)

Create a B-spline basis with order order and breakpoint vector breakpoints. The breakpoint vector is assumed to be sorted.

Derivative{N}

A singleton type that represents the N-th derivative.

Derivative(N)

A shortcut for Derivative{N}(), representing the N-th derivative.

Spline{B<:BSplineBasis, C<:AbstractVector{<:Real}}

Type for a spline based on a B-spline basis of type B and coefficient vector of type C.

Spline(basis, coeffs)

Create a spline from a B-spline basis and a vector of coefficients.

Function(spline::Spline, [deriv::Derivative], [onlysupport=true])

Wrap spline in an object that is a subtype of Function. Calling the returned function with a single argument x will evaluate the spline (or one of its derivatives as specified by deriv) at x.

If the optional argument onlysupport is set to true (the default), the returned function will return NaN if evaluated outside of the support of spline. If onlysupport is set to false, it will return zero there (as does calling the spline directly).

Note that a Spline can be called without wrapping them as described here, although they are not a subtype of Function. Wrapping a spline in a Function object is mainly intended to aid in plotting, which is the rationale behind the onlysupport=true default: when using Plots.jl, this will cause the spline to not be drawn outside of its support.

Examples

julia> spline = approximate(sin, BSplineBasis(5, 0:5)); # create a Spline

julia> f = Function(spline, false); # f is zero outside of the interval [0,5]

julia> g = Function(spline, Derivative(1)); # g is NaN outside of the interval [0,5]

julia> f(1.5) == spline(1.5)
true

julia> g(1.5) == spline(1.5, Derivative(1))
true

julia> f(-1), g(-1)
(0.0, NaN)

approximate(f, basis::BSplineBasis; indices=eachindex(basis)) -> Spline

Approximate the function f in the B-spline basis basis. If indices is supplied, only the basis functions at the given indices of basis are used.

The approximation is calculated by interpolating samples of f at the knotaverages of the basis.

See also: interpolate

Examples

julia> basis = BSplineBasis(6, 0:0.25:1);

julia> spl = approximate(sin, basis, indices=2:length(basis))
Spline{BSplineBasis{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{Float64}}:
 basis: 9-element BSplineBasis{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}:
  order: 6
  breakpoints: 0.0:0.25:1.0
 coeffs: [0.0, 0.05, 0.15, 0.298438, 0.486979, 0.651299, 0.755166, 0.814456, 0.841471]

julia> spl(π/4) ≈ sqrt(1/2)
true

averagebasis(order, datapoints) -> BSplineBasis

Returns a B-spline basis with the specified order that is well-suited for interpolation on the given datapoints. The datapoints vector is assumed to be sorted.

The calculated breakpoints are described in ^[deBoor1978], p. 219, as a “reasonable alternative” to the optimal breakpoint sequence since they are “often very close to the optimum” and are computationally inexpensive.

Examples

julia> averagebasis(5, 0:10)
11-element BSplineBasis{Vector{Float64}}:
 order: 5
 breakpoints: [0.0, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 10.0]

basis(spline::Spline)

Return the B-spline basis on which the spline is defined.

basismatrix(basis::BSplineBasis, xvalues; indices=eachindex(basis))

Calculate the matrix [basis[i](x) for x=xvalues, i=indices].

Examples

julia> basis = BSplineBasis(3, 0:5);

julia> x = [0.3, 1.5, 3.2, 4.5];

julia> basismatrix(basis, x)
4×7 Matrix{Float64}:
 0.49  0.465  0.045  0.0    0.0    0.0    0.0
 0.0   0.125  0.75   0.125  0.0    0.0    0.0
 0.0   0.0    0.0    0.32   0.66   0.02   0.0
 0.0   0.0    0.0    0.0    0.125  0.625  0.25

basismatrix!(dest, basis::BSplineBasis, xvalues; indices=eachindex(basis))

Calculate the matrix [basis[i](x) for x=xvalues, i=indices] and store it in dest.

Examples

julia> basis = BSplineBasis(3, 0:5);

julia> x = [0.3, 1.5, 3.2, 4.5];

julia> dest = Array{Float64}(undef, length(x), length(basis));

julia> basismatrix!(dest, basis, x)
4×7 Matrix{Float64}:
 0.49  0.465  0.045  0.0    0.0    0.0    0.0
 0.0   0.125  0.75   0.125  0.0    0.0    0.0
 0.0   0.0    0.0    0.32   0.66   0.02   0.0
 0.0   0.0    0.0    0.0    0.125  0.625  0.25

breakpoints(basis::BSplineBasis)

Return the breakpoint sequence of the B-spline basis.

Examples

julia> breakpoints(BSplineBasis(3, 0:5))
0:5

julia> breakpoints(BSplineBasis(4, [1.0, 1.5, 2.5, 4.0]))
4-element Vector{Float64}:
 1.0
 1.5
 2.5
 4.0

bsplines!(dest, basis, x[, derivative]; kwargs...)

Calculate the values and/or derivatives of all non-zero B-splines of basis at x in-place (i.e., in dest).

If all B-splines of basis are zero at x, nothing is returned. If any B-splines are non-zero at x, an OffsetArray that wraps dest is returned. Its size and contents depend on the optional derivative argument:

• If derivative is not given, the returned OffsetArray contains the value of the i-th B-spline at the index i. In this case, dest must be a vector of length order(basis).
• If derivative == Derivative(N), the returned OffsetArray contains the N-th derivative of the i-th B-spline at the index i. Again, dest must be a vector of length order(basis).
• If derivative == AllDerivatives(N), the returned OffsetArray contains the m-th derivative (0 ≤ m < N) of the i-th B-spline at the index i, m. In this case, dest must be a matrix of size (order(basis), N).

Two optional keyword arguments can be used to increase performance:

• derivspace: When calculating derivatives, some coefficients are stored in a matrix of size (order(basis), order(basis)). By default, the function allocates a new matrix. To avoid this, a pre-allocated matrix can be supplied with the derivspace keyword. It can only be used when calculating derivatives, i.e., with Derivative(N) where N ≥ 1 or AllDerivatives(N) where N ≥ 2.
• leftknot: If the index of the relevant interval (i.e., intervalindex(basis(spline), x)) is already known, it can be supplied with the leftknot keyword.

Examples

julia> basis = BSplineBasis(4, 0:5);

julia> dest = zeros(4);

julia> bsplines!(dest, basis, 2.4)
4-element OffsetArray(::Vector{Float64}, 3:6) with eltype Float64 with indices 3:6:
 0.03600000000000002
 0.5386666666666667
 0.41466666666666663
 0.01066666666666666

julia> parent(ans) === dest
true

julia> bsplines!(dest, basis, -1.0) # returns nothing

julia> dest = zeros(4, 3);

julia> bsplines!(dest, basis, 3.75, AllDerivatives(3))
4×3 OffsetArray(::Matrix{Float64}, 4:7, 0:2) with eltype Float64 with indices 4:7×0:2:
 0.00260417  -0.03125    0.25
 0.315104    -0.65625    0.25
 0.576823     0.265625  -1.625
 0.105469     0.421875   1.125

julia> parent(ans) === dest
true

bsplines(basis, x[, derivative]; kwargs...)

Calculate the values and/or derivatives of all non-zero B-splines of basis at x.

If all B-splines of basis are zero at x, nothing is returned. If any B-splines are non-zero at x, an OffsetArray is returned. Its size and contents depend on the optional derivative argument:

• If derivative is not given, the returned OffsetArray contains the value of the i-th B-spline at the index i.
• If derivative == Derivative(N), the returned OffsetArray contains the N-th derivative of the i-th B-spline at the index i.
• If derivative == AllDerivatives(N), the returned OffsetArray contains the m-th derivative (0 ≤ m < N) of the i-th B-spline at the index i, m.

Two optional keyword arguments can be used to increase performance:

• derivspace: When calculating derivatives, some coefficients are stored in a matrix of size (order(basis), order(basis)). By default, the function allocates a new matrix. To avoid this, a pre-allocated matrix can be supplied with the derivspace keyword. It can only be used when calculating derivatives, i.e., with Derivative(N) where N ≥ 1 or AllDerivatives(N) where N ≥ 2. To also pre-allocate the array that contains the result, use the bsplines! function instead.
• leftknot: If the index of the relevant interval (i.e., intervalindex(basis(spline), x)) is already known, it can be supplied with the leftknot keyword.

Examples

julia> basis = BSplineBasis(4, 0:5);

julia> bsplines(basis, 2.4)
4-element OffsetArray(::Vector{Float64}, 3:6) with eltype Float64 with indices 3:6:
 0.03600000000000002
 0.5386666666666667
 0.41466666666666663
 0.01066666666666666

julia> bsplines(basis, 2.4, Derivative(1), derivspace=zeros(4,4))
4-element OffsetArray(::Vector{Float64}, 3:6) with eltype Float64 with indices 3:6:
 -0.18000000000000005
 -0.5599999999999999
  0.66
  0.07999999999999996

julia> bsplines(basis, 6) # returns nothing

julia> bsplines(basis, 17//5, leftknot=7)
4-element OffsetArray(::Vector{Rational{Int64}}, 4:7) with eltype Rational{Int64} with indices 4:7:
   9//250
 202//375
 307//750
   2//125

julia> bsplines(basis, 2.4, AllDerivatives(3))
4×3 OffsetArray(::Matrix{Float64}, 3:6, 0:2) with eltype Float64 with indices 3:6×0:2:
 0.036      -0.18   0.6
 0.538667   -0.56  -0.8
 0.414667    0.66  -0.2
 0.0106667   0.08   0.4

coeffs(spline::Spline)

Return the coefficient vector of spline.

interpolate(basis::BSplineBasis, xvalues, yvalues; indices=eachindex(basis)) -> Spline

Interpolate the data given by xvalues and yvalues in the B-spline basis basis. If indices is supplied, only the basis functions at the given indices of basis are used.

The spline interpolation is calculated by creating the matrix B = [basis[i](x) for x=xvalues, i=indices] and then calculating B\yvalues.

See also: approximate

Examples

julia> basis = BSplineBasis(5, 1:10);

julia> xs = range(1, stop=10, length=length(basis)); ys = log.(xs);

julia> spl = interpolate(basis, xs, ys)
Spline{BSplineBasis{UnitRange{Int64}}, Vector{Float64}}:
 basis: 13-element BSplineBasis{UnitRange{Int64}}:
  order: 5
  breakpoints: 1:10
 coeffs: [0.0, 0.248019, 0.596872, 0.946671, 1.26894, 1.51405, 1.71149, 1.87666, 2.01856, 2.14292, 2.22592, 2.27758, 2.30259]

julia> spl(float(ℯ))
0.9999766059171411

intervalindex(vec, x[, start])

If v is an AbstractVector, return the largest index i so that vec[i] ≤ x ≤ vec[end] and vec[i] < vec[end]. Return nothing if no such index exists. The vector vec is assumed to be sorted in ascending order.

If vec is a BSplineBasis, return intervalindex(knots(vec), x[, start]).

If start is given, a linear search is performed, starting from the index start going forward or backward. If start is not given, a binary search is performed.

Examples

julia> intervalindex([1,1,2,3,4,4,4,5,6,6], 2.5)
3

julia> intervalindex([1,1,2,3,4,4,4,5,6,6], 7) # returns nothing

julia> intervalindex([1,1,2,3,4,4,4,5,6,6], 1)
2

julia> intervalindex([1,1,2,3,4,4,4,5,6,6], 4)
7

julia> intervalindex([1,1,2,3,4,4,4,5,6,6], 6.0)
8

intervalindices(basis::BSplineBasis, indices=eachindex(basis))

Return an iterator that yields the indices of all intervals on which basis is defined, i.e., it produces all indices ind (in ascending order) for which (knots(basis)[ind], knots(basis)[ind+1]) is such an interval.

If a range of indices is supplied, the iterator yields only those intervals on which at least one of the B-splines basis[j] for j=indices is non-zero.

Examples

julia> intervalindices(BSplineBasis(3, 0:5))
3:7

julia> intervalindices(BSplineBasis(3, 0:5), 1:4)
3:6

julia> intervalindices(BSplineBasis(4, [1,2,3,4,4,4,5,6]))
BSplines.IntervalIndices{Vector{Int64}}([1, 2, 3, 4, 4, 4, 5, 6], 1:8, 3)

julia> collect(ans)
5-element Vector{Int64}:
  4
  5
  6
  9
 10

intervalindices(basis::BSplineBasis, i, j, ...)

For integers i, j, …, return an iterator that yields the indices of all intervals on which all of the B-splines basis[i], basis[j], … are non-zero, i.e., it produces all indices ind (in ascending order) for which (knots(basis)[ind], knots(basis)[ind+1]) is such an interval.

Examples

julia> intervalindices(BSplineBasis(3, 0:5), 3)
3:5

julia> intervalindices(BSplineBasis(3, 0:5), 4, 5)
5:6

julia> intervalindices(BSplineBasis(3, 0:5), 2, 6) # B-splines do not overlap
6:5

julia> intervalindices(BSplineBasis(3, 0:5), 3, 5, 4)
5:5

julia> intervalindices(BSplineBasis(4, [1,2,3,4,4,4,5,6]), 3, 5)
BSplines.IntervalIndices{Vector{Int64}}([1, 2, 3, 4, 4, 4, 5, 6], 2:4, 3)

julia> collect(ans)
2-element Vector{Int64}:
 5
 6

knotaverages!(dest, basis::BSplineBasis; indices=eachindex(basis))

Calculate the knot averages τ[i] = mean(knots[i+1:i+order-1]) for i ∈ indices and the knots and order of basis and store the result in dest. The knot averages are recommended in ^[deBoor1978] (p. 214) as data points for interpolation.

See also: knotaverages

Examples

julia> dest = Vector{Float64}(undef, 7);

julia> knotaverages!(dest, BSplineBasis(3, 0:5))
7-element Vector{Float64}:
 0.0
 0.5
 1.5
 2.5
 3.5
 4.5
 5.0

julia> dest = Vector{Rational{Int}}(undef, 5);

julia> knotaverages!(dest, BSplineBasis(3, 0:5), indices=2:6)
5-element Vector{Rational{Int64}}:
 1//2
 3//2
 5//2
 7//2
 9//2

knotaverages(basis::BSplineBasis; indices=eachindex(basis))

Return the knot averages τ[i] = mean(knots[i+1:i+order-1]) for i ∈ indices and the knots and order of basis. The knot averages are recommended in ^[deBoor1978] (p. 214) as data points for interpolation.

See also: knotaverages!

Examples

julia> knotaverages(BSplineBasis(3, 0:5))
7-element Vector{Float64}:
 0.0
 0.5
 1.5
 2.5
 3.5
 4.5
 5.0

julia> knotaverages(BSplineBasis(4, [1, 3//2, 5//2, 4]), indices=2:6)
5-element Vector{Rational{Int64}}:
 7//6
 5//3
 8//3
 7//2
 4//1

knots(basis::BSplineBasis)

Return the knot sequence of the B-spline basis.

The knot sequence is the breakpoint sequence except that the first and last values are duplicated so they appear order(basis) times.

Examples

julia> knots(BSplineBasis(3, 0:5))
10-element BSplines.KnotVector{Int64, UnitRange{Int64}}:
 0
 0
 0
 1
 2
 3
 4
 5
 5
 5

order(spline::Spline)
order(basis::BSplineBasis)

Return the order of a spline or a B-spline basis.

Examples

julia> order(BSplineBasis(3, 0:5))
3

julia> order(BSplineBasis(4, [1.0, 1.5, 2.5, 4.0]))
4

splinevalue(spline::Spline, x[, ::Derivative{N}]; kwargs...)

Calculate the value of spline, or its N-th derivative, at x.

Two optional keyword arguments can be used to increase performance:

• workspace: By default, the function allocates a vector of length order(spline) in which the calculation is performed. To avoid this, a pre-allocated vector can be supplied with the workspace keyword. In this case, the returned value is always of type eltype(workspace).
• leftknot: If the index of the relevant interval (i.e., intervalindex(basis(spline), x)) is already known, it can be supplied with the leftknot keyword.

Instead of calling splinevalue, a spline object can be called directly: spline(x[, Derivative(N)]; kwargs...) is equivalent to splinevalue(spline, x[, Derivative(N)]; kwargs...).

Examples

julia> spl = Spline(BSplineBasis(4, 0:5), 1:8);

julia> splinevalue(spl, 1.7)
3.69775

julia> splinevalue(spl, 1.7, Derivative(1))
1.0225

julia> splinevalue(spl, 3.6, leftknot=7)
5.618

julia> spl(18//5)
2809//500

julia> spl(3.6, Derivative(3), leftknot=7)
0.5

support(basis::BSplineBasis) -> a, b

Return the interval $[a,b]$ on which the B-spline basis is defined, i.e., a is the first and b the last breakpoint of the basis.

Examples

julia> support(BSplineBasis(3, 0:5))
(0, 5)

julia> support(BSplineBasis(4, [1.0, 1.5, 2.5, 4.0]))
(1.0, 4.0)

support(spline::Spline) -> a, b

If spline is a BSpline, return the interval $[a,b]$ on which the B-spline is non-zero. Otherwise, return support(basis(spline)).

Examples

julia> basis = BSplineBasis(3, 0:5);

julia> spline = Spline(basis, ones(7));

julia> zerospline = Spline(basis, zeros(7));

julia> support(spline)
(0, 5)

julia> support(zerospline) # even though the spline is zero everywhere
(0, 5)

julia> support(basis[4]) # for BSplines, return their actual support
(1, 4)