# The `BSplineBasis`

type and related functions

To work with B-spline bases, this package defines the parametric `BSplineBasis{T}`

type which represents B-spline bases with breakpoint vector of type `T`

.

A B-spline basis is completely characterized by its order $k$ and knot vector $t$. In the case of the `BSplineBasis`

type, **the knot vector of a basis is generated from its breakpoint vector by repeating the first and last breakpoints so that they appear $k$ times.**

In this package, $k$ always refers to the *order* of a B-spline. Some B-spline libraries, including DIERCKX and scipy.interpolate, use the symbol $k$ to refer to the *degree* of a B-spline. The relationship between the two is

\[\mathrm{order} = \mathrm{degree} + 1.\]

For example, using $k=4$ in this package is equivalent to $k=3$ (cubic splines) in scipy.interpolate.

Knot sequences where the first and last knot do not appear $k$ times are not supported by the `BSplineBasis`

type. The reason for this is that it simplifies the evaluation of B-splines, since it means that exactly $k$ B-splines are non-zero on each interval. If the first or last knot would appear less than $k$ times, this would not be the case.

## Properties of B-spline bases – `order`

, `breakpoints`

, and `knots`

A `BSplineBasis`

is constructed from its order and breakpoint vector. The `order`

and `breakpoints`

functions return these properties. The `length`

function returns the number of B-splines in the basis:

```
julia> basis = BSplineBasis(4, 0:5)
8-element BSplineBasis{UnitRange{Int64}}:
order: 4
breakpoints: 0:5
julia> order(basis)
4
julia> breakpoints(basis)
0:5
julia> length(basis)
8
```

To obtain a valid B-spline basis, the breakpoint vector must be sorted in ascending order. This is not checked by the `BSplineBasis`

constructor.

The `knots`

function returns the knot vector that is generated from the breakpoints. In order to not allocate memory for a new array, a wrapper type around the original breakpoint vector is used:

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

The interval over which a `BSplineBasis`

is defined (i.e., the interval between the first and the last breakpoint) can be obtained (as a `Tuple`

) with the `support`

function:

```
julia> support(basis)
(0, 5)
```

## Evaluating B-Splines and their derivatives

The `bsplines`

and `bsplines!`

functions can be used to obtain the values (or derivatives) of all B-splines that are non-zero at a given point.

### Evaluating B-splines

If `x`

is within the support of the B-spline basis, `bsplines(basis, x)`

returns an `OffsetArray`

which contains the value of the `i`

-th B-spline at the index `i`

. The array always has the length `order(basis)`

:

```
julia> basis = BSplineBasis(4, 0:5);
julia> bsplines(basis, 3.2)
4-element OffsetArray(::Vector{Float64}, 4:7) with eltype Float64 with indices 4:7:
0.08533333333333327
0.6306666666666667
0.28200000000000014
0.0020000000000000052
julia> bsplines(basis, 7//3)
4-element OffsetArray(::Vector{Rational{Int64}}, 3:6) with eltype Rational{Int64} with indices 3:6:
4//81
31//54
10//27
1//162
```

The type of the elements of the array depends on the type of the breakpoints and the type of `x`

. If the point `x`

is outside of the support of `basis`

, `nothing`

is returned instead:

`julia> bsplines(basis, 6) # returns `nothing`, which is not printed in the REPL`

### Evaluating derivatives

The `bsplines`

function can also calculate derivatives of the B-splines. The `N`

-th derivative is specified via an optional third argument of the singleton type `Derivative{N}`

. Instead of `Derivative{N}()`

, the shorter constructor `Derivative(N)`

can be used:

```
julia> basis = BSplineBasis(4, 0:5);
julia> bsplines(basis, 7//3, Derivative(1)) # calculate first derivatives of non-zero B-splines
4-element OffsetArray(::Vector{Rational{Int64}}, 3:6) with eltype Rational{Int64} with indices 3:6:
-2//9
-1//2
2//3
1//18
```

To calculate the value of the non-zero B-splines as well as all of their derivatives up to the `N-1`

-th, the `AllDerivatives{N}`

singleton type can be used instead. `AllDerivatives(N)`

is equivalent to `AllDerivatives{N}()`

. If `x`

is within the support of the B-spline basis, `bsplines(basis, x, AllDerivatives(N))`

returns an `OffsetArray`

which contains the `j`

-th derivative of the `i`

-th B-spline at the index `i, j`

(the zeroth derivative is the value itself):

```
julia> bsplines(basis, 4, AllDerivatives(3)) # calculate zeroth, first and second derivatives
4×3 OffsetArray(::Matrix{Float64}, 5:8, 0:2) with eltype Float64 with indices 5:8×0:2:
0.166667 -0.5 1.0
0.583333 -0.25 -2.5
0.25 0.75 1.5
0.0 0.0 0.0
```

### Pre-allocating output arrays

The `bsplines`

function allocates a new array to return the values (except when it returns `nothing`

). In order to use a pre-allocated array instead, the `bsplines!`

function can be used: `bsplines!(dest, args...)`

behaves like `bsplines(args...)`

, but the calculations are done in `dest`

and the returned `OffsetArray`

is a wrapper around `dest`

.

```
julia> basis = BSplineBasis(4, 0:5);
julia> destvec = zeros(4);
julia> bsplines!(destvec, basis, 3.2)
4-element OffsetArray(::Vector{Float64}, 4:7) with eltype Float64 with indices 4:7:
0.08533333333333327
0.6306666666666667
0.28200000000000014
0.0020000000000000052
julia> parent(ans) === destvec
true
julia> bsplines!(destvec, basis, 7//3, Derivative(2))
4-element OffsetArray(::Vector{Float64}, 3:6) with eltype Float64 with indices 3:6:
0.6666666666666665
-0.9999999999999996
-4.440892098500626e-16
0.3333333333333335
julia> parent(ans) === destvec
true
julia> destmat = zeros(Rational{Int}, 4, 2);
julia> bsplines!(destmat, basis, 4, AllDerivatives(2))
4×2 OffsetArray(::Matrix{Rational{Int64}}, 5:8, 0:1) with eltype Rational{Int64} with indices 5:8×0:1:
1//6 -1//2
7//12 -1//4
1//4 3//4
0//1 0//1
julia> parent(ans) === destmat
true
```

When calculating values of B-splines or their derivatives via the `Derivative{N}`

argument, `dest`

must be a vector of length `order(basis)`

. In the case of the `AllDerivatives{N}`

argument, it must be a matrix of size `(order(basis), N)`

. In any case, `dest`

must not have offset axes itself.

### Improving performance

The `bsplines`

and `bsplines!`

functions accept two optional keyword arguments (one of them only when evaluating derivatives) can be used to speed up the evaluation of splines:

`derivspace`

: When evaluating derivatives, coefficients are stored in a matrix of size`(order(basis), order(basis))`

. In order to avoid allocating a new matrix every time, a pre-allocated matrix can be supplied with the`derivspace`

argument. It can only be used when calculating derivatives, i.e., with`Derivative(N)`

where`N ≥ 1`

or`AllDerivatives(N)`

where`N ≥ 2`

.`leftknot`

: Evaluating B-splines at a point`x`

requires finding the largest index`i`

such that`t[i] ≤ x`

and`t[i] < t[end]`

where`t`

is the knot vector. The`bsplines`

and`bsplines!`

functions use the`intervalindex`

function to find this index. If the index is already known, it can be specified with the`leftknot`

keyword argument.

```
julia> using BenchmarkTools
julia> basis = BSplineBasis(4, 0:5);
julia> dest = zeros(order(basis));
julia> space = zeros(order(basis), order(basis));
julia> left = intervalindex(basis, 2.5);
julia> @btime bsplines!($dest, $basis, 2.5, Derivative(1));
541.320 ns (2 allocations: 240 bytes)
julia> @btime bsplines!($dest, $basis, 2.5, Derivative(1), derivspace=$space);
554.131 ns (1 allocation: 32 bytes)
julia> @btime bsplines!($dest, $basis, 2.5, Derivative(1), derivspace=$space, leftknot=$left);
547.200 ns (0 allocations: 0 bytes)
```

## Constructing knot vectors

As stated above, the `knots`

of a `BSplineBasis`

are generated from its `breakpoints`

by repeating the first and last breakpoints so that they appear `order(basis)`

times. However, in many situations it may be desirable to

- repeat other knots than the first and last or
- have the first and/or last knot repeated less than $k$ times (e.g. to implement certain boundary conditions).

In order to have repeated interior knots, it is sufficient to include it in the breakpoint vector multiple times when constructing the basis:

```
julia> basis = BSplineBasis(3, [1, 2, 3, 4, 5, 5, 6, 7, 8]) # 5 appears twice
10-element BSplineBasis{Vector{Int64}}:
order: 3
breakpoints: [1, 2, 3, 4, 5, 5, 6, 7, 8]
julia> knots(basis)
13-element BSplines.KnotVector{Int64, Vector{Int64}}:
1
1
1
2
3
4
5
5
6
7
8
8
8
```

There is currently no way to create a `BSplineBasis`

where the first and last knots appear less than $k$ times. However, some functions provided by this package (see Higher-level functions) provide a `indices`

keyword argument to select only a certain range of B-splines from a basis, thus achieving the same result as if the knot vector had been shortened.

Because of repeated knots, not every pair of knots `(t[i], t[i+1])`

describes an actual interval. The `intervalindices`

function helps with finding all indices which describe relevant intervals. See its docstrings for more information.