### Initialize Submatrix Slice

Source: https://jipolanco.github.io/BSplineKit.jl/dev/tensors

Example of initializing a specific submatrix slice.

```julia
A[:, :, k] = rand(2b + 1, 2b + 1)
```

--------------------------------

### Initialize and Set BandedTensor3D Elements

Source: https://jipolanco.github.io/BSplineKit.jl/dev/tensors

Example of initializing a tensor and setting its slices.

```julia
A = BandedTensor3D(undef, 20, Val(2))  # tensor of size 20×20×20 and band width b = 2
for k in axes(A, 3)
    A[:, :, k] = rand(5, 5)
end
```

--------------------------------

### Construct RecombinedBSplineBasis with multiple operators

Source: https://jipolanco.github.io/BSplineKit.jl/dev/recombination

Examples of applying multiple boundary conditions and natural conditions.

```julia
julia> ops = (Derivative(0), Derivative(1));


julia> R1 = RecombinedBSplineBasis(B, ops)
9-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{0}, D{1})
 BCs right: (D{0}, D{1})

julia> ops = (Derivative(0), Derivative(1) - 4Derivative(2));


julia> R2 = RecombinedBSplineBasis(B, ops)
9-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{0}, D{1} + -4 * D{2})
 BCs right: (D{0}, D{1} + -4 * D{2})

julia> R3 = RecombinedBSplineBasis(B, Natural())
11-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{2},)
 BCs right: (D{2},)
```

--------------------------------

### Perform spline interpolation with examples

Source: https://jipolanco.github.io/BSplineKit.jl/dev/interpolation

Demonstrates creating standard and natural spline interpolations, and evaluating derivatives.

```julia
julia> xs = -1:0.1:1;

julia> ys = cospi.(xs);

julia> S = interpolate(xs, ys, BSplineOrder(4))
SplineInterpolation containing the 21-element Spline{Float64}:
 basis: 21-element BSplineBasis of order 4, domain [-1.0, 1.0]
 order: 4
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3  …  0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0, 1.0, 1.0, 1.0]
 coefficients: [-1.0, -1.00111, -0.8975, -0.597515, -0.314147, 1.3265e-6, 0.314142, 0.597534, 0.822435, 0.96683  …  0.96683, 0.822435, 0.597534, 0.314142, 1.3265e-6, -0.314147, -0.597515, -0.8975, -1.00111, -1.0]
 interpolation points: -1.0:0.1:1.0

julia> S(-1)
-1.0

julia> (Derivative(1) * S)(-1)
-0.01663433622896893

julia> (Derivative(2) * S)(-1)
10.527273287554928

julia> Snat = interpolate(xs, ys, BSplineOrder(4), Natural())
SplineInterpolation containing the 21-element Spline{Float64}:
 basis: 21-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0], BCs {left => (D{2},), right => (D{2},)}
 order: 4
 knots: [-1.0, -1.0, -1.0, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4  …  0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.0, 1.0, 1.0]
 coefficients: [-0.833333, -0.647516, -0.821244, -0.597853, -0.314057, -2.29076e-5, 0.314148, 0.597532, 0.822435, 0.96683  …  0.96683, 0.822435, 0.597532, 0.314148, -2.29076e-5, -0.314057, -0.597853, -0.821244, -0.647516, -0.833333]
 interpolation points: -1.0:0.1:1.0

julia> Snat(-1)
-1.0

julia> (Derivative(1) * Snat)(-1)
0.2872618670889516

julia> (Derivative(2) * Snat)(-1)
-3.33066907387547e-14
```

--------------------------------

### SplineApproximation output structure

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/approximation

Example output showing the structure of a SplineApproximation object resulting from the L2 minimization.

```text
SplineApproximation containing the 17-element Spline{Float64}:
 basis: 17-element BSplineBasis of order 4, domain [0.0, 1.0]
 order: 4
 knots: [0.0, 0.0, 0.0, 0.0, 0.0714286, 0.142857, 0.214286, 0.285714, 0.357143, 0.428571  …  0.571429, 0.642857, 0.714286, 0.785714, 0.857143, 0.928571, 1.0, 1.0, 1.0, 1.0]
 coefficients: [0.97272, 1.0617, -0.457591, -1.31138, 0.91464, 0.769569, -1.11789, -0.20053, 1.05012, -0.259743, -0.804373, 0.560952, 0.460536, -0.66901, -0.138778, 0.409463, 0.358481]
 approximation method: MinimiseL2Error()
```

--------------------------------

### Example usage of PeriodicBSplineBasis

Source: https://jipolanco.github.io/BSplineKit.jl/dev/bsplines

Demonstrates creating a periodic basis and querying its properties, including period, length, boundaries, and knot values.

```julia
julia> ts = range(-1, 1; length = 11)
-1.0:0.2:1.0

julia> B = PeriodicBSplineBasis(BSplineOrder(4), ts)
10-element PeriodicBSplineBasis of order 4, domain [-1.0, 1.0), period 2.0
 knots: [..., -1.4, -1.2, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, ...]

julia> period(B)
2.0

julia> length(B)
10

julia> boundaries(B)
(-1.0, 1.0)

julia> B(-0.42)
(5, (0.12150000000000002, 0.6571666666666667, 0.22116666666666668, 0.00016666666666666563))

julia> B(-0.42 + 2)
(15, (0.12150000000000015, 0.6571666666666667, 0.22116666666666657, 0.00016666666666666674))

julia> knots(B)
14-element PeriodicKnots{Float64, 4, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}:
 -1.4
 -1.2
 -1.0
 -0.8
 -0.6
 -0.4
 -0.2
  0.0
  0.2
  0.4
  0.6
  0.8
  1.0
  1.2
```

--------------------------------

### Periodic boundary condition interpolation

Source: https://jipolanco.github.io/BSplineKit.jl/dev/interpolation

Example of using Periodic boundary conditions for periodic signals.

```julia
julia> xp = -1:0.1:0.9;

julia> yp = cospi.(xp);

julia> Sper = interpolate(xp, yp, BSplineOrder(4), Periodic(2))
SplineInterpolation containing the 20-element Spline{Float64}:
 basis: 20-element PeriodicBSplineBasis of order 4, domain [-1.0, 1.0), period 2.0
 order: 4
 knots: [..., -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3  …  0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, ...]
 coefficients: [..., -1.01659, -0.96683, -0.822435, -0.597534, -0.314142, 1.10589e-17, 0.314142, 0.597534, 0.822435, 0.96683, 1.01659, 0.96683, 0.822435, 0.597534, 0.314142, 1.51788e-17, -0.314142, -0.597534, -0.822435, -0.96683, ...]
 interpolation points: -1.0:0.1:0.9
```

--------------------------------

### Construct RecombinedBSplineBasis with differential operators

Source: https://jipolanco.github.io/BSplineKit.jl/dev/recombination

Examples of creating recombined bases using Neumann and Robin boundary conditions.

```julia
julia> B = BSplineBasis(BSplineOrder(4), -1:0.2:1)
13-element BSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]

julia> R_neumann = RecombinedBSplineBasis(B, Derivative(1))
11-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{1},)
 BCs right: (D{1},)

julia> R_robin = RecombinedBSplineBasis(B, Derivative(0) + 3Derivative(1))
11-element RecombinedBSplineBasis of order 4, domain [-1.0, 1.0]
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 BCs left:  (D{0} + 3 * D{1},)
 BCs right: (D{0} + 3 * D{1},)
```

--------------------------------

### Get Band Shifts

Source: https://jipolanco.github.io/BSplineKit.jl/dev/tensors

Retrieve the band shifts along each dimension.

```julia
bandshift(A::BandedTensor3D) -> (δi, δj, δk)
```

--------------------------------

### Spline Interpolation Object Structure

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/interpolation

Example output showing the internal structure of a SplineInterpolation object with natural boundary conditions applied.

```text
SplineInterpolation containing the 20-element Spline{Float64}:
 basis: 20-element RecombinedBSplineBasis of order 8, domain [0.0, 1.0], BCs {left => (D{2}, D{3}, D{4}), right => (D{2}, D{3}, D{4})}
 order: 8
 knots: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0518416, 0.109298  …  0.900587, 0.933845, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
 coefficients: [0.0137583, 0.142842, 0.381383, 0.482671, 0.653978, 0.576594, 1.10858, 0.718314, 1.09114, 0.975321, 0.940988, 1.12046, 0.863254, 0.931811, 0.613254, 0.83853, 0.206394, 0.679701, 0.046989, -0.0095833]
 interpolation points: [0.0, 0.0518416, 0.109298, 0.160794, 0.209856, 0.269298, 0.325959, 0.367882, 0.415349, 0.465107, 0.524201, 0.570411, 0.635975, 0.671061, 0.738611, 0.786926, 0.824332, 0.900587, 0.933845, 1.0]
```

--------------------------------

### Get number of constraints for a basis

Source: https://jipolanco.github.io/BSplineKit.jl/dev/recombination

Returns the number of boundary conditions satisfied on each boundary for a given basis or matrix.

```julia
num_constraints(R::AbstractBSplineBasis) -> (Int, Int)
num_constraints(A::RecombineMatrix) -> (Int, Int)
```

--------------------------------

### Get Basis Function Support

Source: https://jipolanco.github.io/BSplineKit.jl/dev/bsplines

Retrieves the knot range supported by a basis function or a specific index in a basis.

```julia
support(b::BasisFunction) -> UnitRange{Int}
```

```julia
support(B::AbstractBSplineBasis, i::Integer) -> UnitRange{Int}
```

--------------------------------

### Get Band Width

Source: https://jipolanco.github.io/BSplineKit.jl/dev/tensors

Retrieve the band width b of the tensor.

```julia
bandwidth(A::BandedTensor3D)
```

--------------------------------

### max_order

Source: https://jipolanco.github.io/BSplineKit.jl/dev/diffops

Get the maximum derivative order of one or more differential operators.

```APIDOC
## max_order(op::AbstractDifferentialOp)
## max_order(ops...)

### Description
Returns the maximum derivative order present in the provided differential operator or collection of operators.

### Signature
`max_order(op::AbstractDifferentialOp) -> Int`
`max_order(ops...) -> Int`
```

--------------------------------

### Get length of B-spline basis

Source: https://jipolanco.github.io/BSplineKit.jl/dev/bsplines

Returns the number of B-splines composing a spline.

```julia
length(g::BSplineBasis) -> Int
```

--------------------------------

### Approximate functions with splines in Julia

Source: https://jipolanco.github.io/BSplineKit.jl/dev/approximation

Demonstrates creating a spline basis and approximating functions using different methods, including in-place updates for efficiency.

```julia
julia> B = BSplineBasis(BSplineOrder(3), -1:0.4:1);

julia> S_interp = approximate(sin, B)
SplineApproximation containing the 7-element Spline{Float64}:
 basis: 7-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.6, -0.2, 0.2, 0.6, 1.0, 1.0, 1.0]
 coefficients: [-0.841471, -0.731727, -0.39727, 2.85767e-17, 0.39727, 0.731727, 0.841471]
 approximation method: interpolation at [-1.0, -0.8, -0.4, 0.0, 0.4, 0.8, 1.0]

julia> sin(0.3), S_interp(0.3)
(0.29552020666133955, 0.2959895327282942)

julia> approximate!(exp, S_interp)
SplineApproximation containing the 7-element Spline{Float64}:
 basis: 7-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.6, -0.2, 0.2, 0.6, 1.0, 1.0, 1.0]
 coefficients: [0.367879, 0.440373, 0.65701, 0.980127, 1.46223, 2.18111, 2.71828]
 approximation method: interpolation at [-1.0, -0.8, -0.4, 0.0, 0.4, 0.8, 1.0]

julia> exp(0.3), S_interp(0.3)
(1.3498588075760032, 1.3491015490105396)

julia> S_fast = approximate(exp, B, VariationDiminishing())
SplineApproximation containing the 7-element Spline{Float64}:
 basis: 7-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.6, -0.2, 0.2, 0.6, 1.0, 1.0, 1.0]
 coefficients: [0.367879, 0.449329, 0.67032, 1.0, 1.49182, 2.22554, 2.71828]
 approximation method: VariationDiminishing()

julia> S_opt = approximate(exp, B, MinimiseL2Error())
SplineApproximation containing the 7-element Spline{Float64}:
 basis: 7-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.6, -0.2, 0.2, 0.6, 1.0, 1.0, 1.0]
 coefficients: [0.368074, 0.440342, 0.657077, 0.980279, 1.46216, 2.18201, 2.71669]
 approximation method: MinimiseL2Error()

julia> x = 0.34; exp(x), S_opt(x), S_interp(x), S_fast(x)
(1.4049475905635938, 1.4044530324752076, 1.4044149581073813, 1.4328668494041878)
```

--------------------------------

### Constructing Derivative Ranges in Julia

Source: https://jipolanco.github.io/BSplineKit.jl/dev/diffops

Demonstrates different syntax options for creating a range of derivatives and converting them to a tuple.

```julia
julia> Derivative(2):Derivative(4)
Derivative(2:4)

julia> Derivative(2:4)
Derivative(2:4)

julia> Tuple(Derivative(2:4))
(D{2}, D{3}, D{4})
```

--------------------------------

### Accessing the BoundaryConditions Module

Source: https://jipolanco.github.io/BSplineKit.jl/dev/boundary_conditions

Reference the module containing boundary condition definitions.

```julia
BSplineKit.BoundaryConditions
```

--------------------------------

### Get number of recombined functions

Source: https://jipolanco.github.io/BSplineKit.jl/dev/recombination

Returns the count of basis functions modified by recombination on each boundary.

```julia
num_recombined(R::AbstractBSplineBasis) -> (Int, Int)
num_recombined(A::RecombineMatrix) -> (Int, Int)
```

--------------------------------

### Supported operator configurations

Source: https://jipolanco.github.io/BSplineKit.jl/dev/recombination

Syntax for defining various boundary condition operator sets.

```julia
ops = (Derivative(0), ..., Derivative(m))
```

```julia
ops = (Derivative(0), ..., Derivative(m), D(n))
```

```julia
ops = Natural()
```

```julia
ops = ()
ops = nothing
```

--------------------------------

### Get Band Indices

Source: https://jipolanco.github.io/BSplineKit.jl/dev/tensors

Return the range of indices where values may be non-zero for a given slice.

```julia
band_indices(A::BandedTensor3D, k)
```

--------------------------------

### Visualize Approximation

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/heat

Plots the exact initial condition against its spline approximation and displays the difference between them.

```julia
fig = Figure(size = (800, 400))
let ax = Axis(fig[1, 1]; xlabel = rich("x"; font = :italic), ylabel = rich("θ"; font = :italic))
    lines!(ax, -1..1, θ₀; label = rich("θ", subscript("0"), rich("(x)"; offset = (0.1, 0.0)); font = :italic), color = :blue)
    lines!(ax, -1..1, θ₀_spline; label = "Approximation", color = :orange, linestyle = :dash)
    axislegend(ax; position = :cb)
end
let ax = Axis(fig[1, 2]; xlabel = rich("x"; font = :italic), ylabel = "Difference")
    lines!(ax, -1..1, x -> θ₀(x) - θ₀_spline(x))
    plot_knots!(ax, knots(R); knot_offset = 0)
end
fig
```

--------------------------------

### Get recombination matrix

Source: https://jipolanco.github.io/BSplineKit.jl/dev/recombination

Retrieves the RecombineMatrix associated with a recombined basis, or the identity matrix for non-recombined bases.

```julia
recombination_matrix(R::AbstractBSplineBasis)
```

--------------------------------

### Construct Galerkin matrices

Source: https://jipolanco.github.io/BSplineKit.jl/dev

Generates mass and stiffness matrices for the Galerkin method using a recombined basis.

```julia
# By default, M and L are Hermitian banded matrices
M = galerkin_matrix(R)
L = galerkin_matrix(R, (Derivative(1), Derivative(1)))
```

--------------------------------

### Initialize ApproxByInterpolation with basis

Source: https://jipolanco.github.io/BSplineKit.jl/dev/approximation

Automatically determines interpolation points based on the provided B-spline basis.

```julia
ApproxByInterpolation(B::AbstractBSplineBasis)
```

--------------------------------

### Time integration and visualization

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/heat

Define the right-hand side function for the ODE solver, set up the problem, and visualize the solution.

```julia
using LinearAlgebra
using OrdinaryDiffEqTsit5

function heat_rhs!(du, u, params, t)
    mul!(du, params.L, u)    # du = ν * L * u
    ldiv!(du, params.A, du)  # du = A \ (ν * L * u)
    du
end

# Solver parameters
params_col = (
    A = lu(Acol),  # we pass the factorised matrix A for performance
    L = Lcol,
)

tspan = (0.0, 10.0)
prob = ODEProblem(heat_rhs!, u_init, tspan, params_col)
prob = ODEProblem{true}(heat_rhs!, u_init, tspan, params_col)
sol_collocation = solve(prob, Tsit5(); saveat = 0.5)

function plot_heat_solution(sol, R)
    fig = Figure()
    ax = Axis(fig[1, 1]; xlabel = rich("x"; font = :italic), ylabel = rich("θ(x,t)"; font = :italic))
    colormap = cgrad(:viridis)
    tspan = sol.prob.tspan
    Δt = tspan[2] - tspan[1]
    for (u, t) in tuples(sol)
        S = Spline(R, u)
        color = colormap[(t - tspan[1]) / Δt]
        lines!(ax, -1..1, S; label = string(t), color, linewidth = 2)
    end
    Colorbar(fig[1, 2]; colormap, limits = tspan, label = "Time")
    fig
end

# NOTE: there's an issue in CairoMakie 0.11.10 when saving SVGs with colourbars, so we fall
# back to PNG output.
# See https://github.com/MakieOrg/Makie.jl/issues/3016
CairoMakie.activate!(type = "png", px_per_unit = 2.0)
plot_heat_solution(sol_collocation, R)
```

--------------------------------

### Create a recombined basis with boundary conditions

Source: https://jipolanco.github.io/BSplineKit.jl/dev

Generates a basis that satisfies homogeneous Robin boundary conditions.

```julia
bc = Derivative(0) + 3Derivative(1)
R = RecombinedBSplineBasis(B, bc)  # satisfies u ∓ 3u' = 0 on the left/right boundary
```

--------------------------------

### Construct the stiffness matrix R

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/heat

Construct the matrix R representing the inner product of first derivatives.

```julia
Rgal = galerkin_matrix(R, (Derivative(1), Derivative(1)))
```

--------------------------------

### Get Non-zero Basis Functions in Segment

Source: https://jipolanco.github.io/BSplineKit.jl/dev/bsplines

Returns the range of indices for basis functions that are non-zero within a specific knot segment.

```julia
nonzero_in_segment(B::AbstractBSplineBasis, n::Int) -> UnitRange{Int}
```

--------------------------------

### Visualize extrapolation methods

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/interpolation

Plotting the different extrapolation strategies using Makie.

```julia
fig = Figure(size = (600, 400))
ax = Axis(fig[1, 1])
scatter!(ax, xs, ys; label = "Data", color = :black)
lines!(ax, -0.5..1.5, S; label = "No extrapolation", linewidth = 2)
lines!(ax, -0.5..1.5, E_smooth; label = "Smooth", linestyle = :dash, linewidth = 2)
lines!(ax, -0.5..1.5, E_linear; label = "Linear", linestyle = :dashdot, linewidth = 2)
lines!(ax, -0.5..1.5, E_flat; label = "Flat", linestyle = :dot, linewidth = 2)
axislegend(ax)
fig
```

--------------------------------

### Initialize ApproxByInterpolation with points

Source: https://jipolanco.github.io/BSplineKit.jl/dev/approximation

Specifies interpolation points for the approximation.

```julia
ApproxByInterpolation(xs::AbstractVector)
```

--------------------------------

### RecombinedBSplineBasis constructor signatures

Source: https://jipolanco.github.io/BSplineKit.jl/dev/recombination

Standard constructor signatures for creating a recombined basis.

```julia
RecombinedBSplineBasis(B::BSplineBasis, ops)
RecombinedBSplineBasis(B::BSplineBasis, ops_left, ops_right)
```

--------------------------------

### Define and Plot the Target Function

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/approximation

Sets up the CairoMakie environment and defines the target function f(x) = e^-x * cos(8πx) for visualization.

```julia
using CairoMakie
CairoMakie.activate!(type = "svg", pt_per_unit = 2.0)

x_interval = 0..1
f(x) = exp(-x) * cospi(8x)

fig = Figure()
ax = Axis(fig[1, 1]; xlabel = rich("x"; font = :italic))
lines!(ax, x_interval, f)
fig
```

--------------------------------

### constraints(R::AbstractBSplineBasis)

Source: https://jipolanco.github.io/BSplineKit.jl/dev/recombination

Retrieves the boundary condition constraints for a basis.

```APIDOC
## constraints(R::AbstractBSplineBasis)

### Description
Returns a tuple containing the left and right boundary constraints applied to the basis.

### Returns
- **(left, right)** (Tuple) - The constraints applied to the left and right boundaries respectively.
```

--------------------------------

### Visualize spline approximation methods

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/approximation

Uses Makie.jl to plot the original function against three different spline approximation methods and their respective errors.

```julia
fig = Figure(size = (950, 750))
colours = theme(fig.scene).palette.color[]
style_vd = (color = colours[3], label = "Variation diminishing")
style_interp = (color = colours[2], label = "Interpolation")
style_minL2 = (color = colours[1], label = "L² minimisation")
let ax = Axis(fig[1:2, 1]; xlabel = rich("x"; font = :italic), ylabel = "Approximation")
    plot_knots!(ax, knots(B); knot_offset = nothing)
    lines!(ax, x_interval, f; color = :black, linewidth = 2, label = "Original")
    lines!(ax, x_interval, S_vd; style_vd...)
    lines!(ax, x_interval, S_interp; style_interp...)
    lines!(ax, x_interval, S_minL2; style_minL2...)
    axislegend(ax)
end
let ax = Axis(fig[1, 2]; ylabel = "Difference with original")
    plot_knots!(ax, knots(B); knot_offset = nothing)
    lines!(ax, x_interval, x -> S_interp(x) - f(x); style_interp...)
    lines!(ax, x_interval, x -> S_minL2(x) - f(x); style_minL2...)
    hidexdecorations!(ax; grid = false)
    axislegend(ax; position = :rt, orientation = :horizontal)
end
let ax = Axis(fig[2, 2]; xlabel = rich("x"; font = :italic), ylabel = "Squared difference", yscale = log10)
    ylims!(1e-8, 1e-2)
    plot_knots!(ax, knots(B); knot_offset = nothing, ybase = 1e-6)
    lines!(ax, x_interval, x -> abs2(S_interp(x) - f(x)); style_interp...)
    lines!(ax, x_interval, x -> abs2(S_minL2(x) - f(x)); style_minL2...)
end
fig
```

--------------------------------

### galerkin_matrix!

Source: https://jipolanco.github.io/BSplineKit.jl/dev/galerkin

Fills a preallocated matrix with Galerkin mass or stiffness values.

```APIDOC
## galerkin_matrix!(f, A, B, deriv; quadrature)

### Description
Fills an existing preallocated matrix A with the computed Galerkin matrix values. This is useful for performance when the matrix structure is known.

### Parameters
- **f** (Function) - Optional - A function to include in the inner product integral.
- **A** (AbstractMatrix) - Required - The preallocated matrix to fill.
- **B** (AbstractBSplineBasis or Tuple) - Required - The B-spline basis or a tuple of two bases.
- **deriv** (Tuple) - Optional - A tuple of Derivative objects specifying the derivative order.
- **quadrature** (Any) - Optional - The quadrature rule to use for integration.
```

--------------------------------

### Construct the mass matrix

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/heat

Construct the mass matrix A using the galerkin_matrix function.

```julia
Agal = galerkin_matrix(R)
```

--------------------------------

### Fill Preallocated Galerkin Matrix

Source: https://jipolanco.github.io/BSplineKit.jl/dev/galerkin

Fills an existing matrix with Galerkin integration results, useful for performance-critical applications where memory allocation should be minimized.

```julia
galerkin_matrix!(
    [f::Function], A::AbstractMatrix, B::AbstractBSplineBasis, [deriv = (Derivative(0), Derivative(0))];
    [quadrature],
)
```

--------------------------------

### Visualize solution comparison with Makie

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/heat

Uses Makie to plot the final state of the heat equation solutions and their difference for both collocation and Galerkin methods.

```julia
fig = Figure(size = (800, 400))
let ax = Axis(fig[1, 1]; xlabel = rich("x"; font = :italic), ylabel = rich("θ(x, t = $(tspan[end]))"; font = :italic))
    for pair in (
            "Collocation" => sol_collocation,
            "Galerkin" => sol_galerkin,
        )
        label, sol = pair
        u = last(sol.u)
        S = Spline(R, u)
        lines!(ax, -1..1, S; label, linewidth = 2)
    end
    axislegend(ax; position = :cb)
end
let ax = Axis(fig[1, 2]; xlabel = rich("x"; font = :italic), ylabel = "Difference")
    Sc = Spline(R, last(sol_collocation.u))
    Sg = Spline(R, last(sol_galerkin.u))
    lines!(ax, -1..1, x -> Sc(x) - Sg(x); linewidth = 2)
end
fig
```

--------------------------------

### Plot Basis Function Derivatives

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/heat

Visualizes the derivatives of the recombined basis functions to verify boundary conditions.

```julia
fig = Figure()
ax = Axis(
    fig[1, 1];
    xlabel = rich("x"; font = :italic),
    ylabel = rich("ϕ′", subscript("i"; offset = (-0.3, 0.0)), rich("(x)"; offset = (0.1, 0.0)); font = :italic),
)
plot_basis!(ax, R; eval_args = (Derivative(1), ), knot_offset = 0.4)
fig
```

--------------------------------

### Visualize Basis Functions and Knots

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/approximation

Defines helper functions to plot knots and basis functions, then renders the spline space visualization.

```julia
function plot_knots!(ax, ts; ybase = 0, knot_offset = 0.03, kws...)
    ys = zero(ts) .+ ybase
    # Add offset to distinguish knots with multiplicity > 1
    if knot_offset !== nothing
        for i in eachindex(ts)[(begin + 1):end]
            if ts[i] == ts[i - 1]
                ys[i] = ys[i - 1] + knot_offset
            end
        end
    end
    scatter!(ax, ts, ys; marker = '×', color = :gray, markersize = 24, kws...)
    ax
end

function plot_basis!(ax, B; eval_args = (), kws...)
    cmap = cgrad(:tab20)
    N = length(B)
    ts = knots(B)
    hlines!(ax, 0; color = :gray)
    for (n, bi) in enumerate(B)
        color = cmap[(n - 1) / (N - 1)]
        i, j = extrema(support(bi))
        lines!(ax, ts[i]..ts[j], x -> bi(x, eval_args...); color, linewidth = 2.5)
    end
    plot_knots!(ax, ts; kws...)
    ax
end

fig = Figure()
ax = Axis(
    fig[1, 1];
    xlabel = rich("x"; font = :italic),
    ylabel = rich("b", subscript("i"), rich("(x)"; offset = (0.1, 0.0)); font = :italic),
)
plot_basis!(ax, B; knot_offset = 0.05)
fig
```

--------------------------------

### Utility methods for RecombinedBSplineBasis

Source: https://jipolanco.github.io/BSplineKit.jl/dev/recombination

Methods to retrieve the original basis, length, or constraints.

```julia
parent(R::RecombinedBSplineBasis)
```

```julia
length(R::RecombinedBSplineBasis)
```

```julia
constraints(R::AbstractBSplineBasis) -> (left, right)
constraints(A::RecombineMatrix) -> (left, right)
```

--------------------------------

### RecombinedBSplineBasis Constructors

Source: https://jipolanco.github.io/BSplineKit.jl/dev/recombination

Constructors for creating a recombined basis, either with a single differential operator applied to both boundaries or with distinct operators for the left and right boundaries.

```julia
RecombinedBSplineBasis(B::BSplineBasis, op::AbstractDifferentialOp)
RecombinedBSplineBasis(B::BSplineBasis, op_left, op_right)
```

--------------------------------

### Defining Periodic Boundary Conditions

Source: https://jipolanco.github.io/BSplineKit.jl/dev/boundary_conditions

Represents periodic boundary conditions with a specified period.

```julia
Periodic <: BoundaryCondition
```

```julia
Periodic(L::Real)
```

--------------------------------

### Plot interpolation and fitting results

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/interpolation

Visualize the data, periodic interpolation, and smoothing splines using Makie.

```julia
fig = Figure()
ax = Axis(fig[1, 1])
scatter!(ax, xs, ys; label = "Data")
lines!(ax, -1..1, S_interp; label = "Interpolation", color = (:grey, 0.5))
lines!(ax, -1..1, S_fit_periodic; linewidth = 3, label = "Smoothing (periodic)")
lines!(ax, -1..1, S_fit_natural; linewidth = 3, linestyle = :dash, label = "Smoothing (natural)")
axislegend(ax)
fig
```

--------------------------------

### Visualize B-spline basis

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/heat

Plot the knot locations and the generated B-spline basis functions using CairoMakie.

```julia
using CairoMakie
using LaTeXStrings
CairoMakie.activate!(type = "svg", pt_per_unit = 2.0)

function plot_knots!(ax, ts; knot_offset = 0.05, kws...)
    ys = zero(ts)
    # Add offset to distinguish knots with multiplicity > 1
    for i in eachindex(ts)[(begin + 1):end]
        if ts[i] == ts[i - 1]
            ys[i] = ys[i - 1] + knot_offset
        end
    end
    scatter!(ax, ts, ys; marker = '×', markersize = 24, color = :gray, kws...)
    ax
end

function plot_basis!(ax, B; eval_args = (), kws...)
    cmap = cgrad(:tab20)
    N = length(B)
    ts = knots(B)
    hlines!(ax, 0; color = :gray)
    for (n, bi) in enumerate(B)
        color = cmap[(n - 1) / (N - 1)]
        i, j = extrema(support(bi))
        lines!(ax, ts[i]..ts[j], x -> bi(x, eval_args...); color, linewidth = 2.5)
    end
    plot_knots!(ax, ts; kws...)
    ax
end

fig = Figure()
ax = Axis(
    fig[1, 1];
    xlabel = rich("x"; font = :italic),
    ylabel = rich("b", subscript("i"), rich("(x)"; offset = (0.1, 0.0)); font = :italic),
)
plot_basis!(ax, B)
fig
```

--------------------------------

### Initialize SplineInterpolation

Source: https://jipolanco.github.io/BSplineKit.jl/dev/interpolation

Initializes the interpolation structure using a B-spline basis and interpolation points. The number of points must match the number of B-splines.

```julia
SplineInterpolation(undef, B::AbstractBSplineBasis, x::AbstractVector, [T = eltype(x)])
```

--------------------------------

### Construct a Spline

Source: https://jipolanco.github.io/BSplineKit.jl/dev/splines

Create a spline instance using a B-spline basis and a vector of coefficients.

```julia
julia> B = BSplineBasis(BSplineOrder(4), -1:0.2:1);

julia> coefs = rand(length(B));

julia> S = Spline(B, coefs)
13-element Spline{Float64}:
 basis: 13-element BSplineBasis of order 4, domain [-1.0, 1.0]
 order: 4
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 coefficients: [0.173575, 0.321662, 0.258585, 0.166439, 0.527015, 0.483022, 0.390663, 0.802763, 0.721983, 0.372347, 0.0301856, 0.0793339, 0.663758]
```

--------------------------------

### Apply extrapolation strategies

Source: https://jipolanco.github.io/BSplineKit.jl/dev/generated/interpolation

Enable extrapolation by wrapping the spline with a specific strategy like Flat, Linear, or Smooth.

```julia
E_flat   = extrapolate(S, Flat())
E_linear = extrapolate(S, Linear())
E_smooth = extrapolate(S, Smooth())
```

--------------------------------

### Fit a cubic smoothing spline in Julia

Source: https://jipolanco.github.io/BSplineKit.jl/dev/interpolation

Demonstrates fitting a cubic spline to noisy data using a specified smoothing parameter.

```julia
julia> xs = (0:0.01:1).^2;

julia> ys = @. cospi(2 * xs) + 0.1 * sinpi(200 * xs);  # smooth + highly fluctuating components

julia> λ = 0.001;  # smoothing parameter

julia> S = fit(BSplineOrder(4), xs, ys, λ)
101-element Spline{Float64}:
 basis: 101-element RecombinedBSplineBasis of order 4, domain [0.0, 1.0], BCs {left => (D{2},), right => (D{2},)}
 order: 4
 knots: [0.0, 0.0, 0.0, 0.0, 0.0001, 0.0004, 0.0009, 0.0016, 0.0025, 0.0036  …  0.8836, 0.9025, 0.9216, 0.9409, 0.9604, 0.9801, 1.0, 1.0, 1.0, 1.0]
 coefficients: [0.946872, 0.631018, 1.05101, 1.04986, 1.04825, 1.04618, 1.04366, 1.04067, 1.03722, 1.03331  …  0.437844, 0.534546, 0.627651, 0.716043, 0.798813, 0.875733, 0.947428, 1.01524, 0.721199, 0.954231]
```

--------------------------------

### Construct BasisFunction

Source: https://jipolanco.github.io/BSplineKit.jl/dev/bsplines

Creates the i-th basis function for a given basis.

```julia
BasisFunction(basis::AbstractBSplineBasis, i::Int, [T = Float64])
```

--------------------------------

### AvgKnots selection method

Source: https://jipolanco.github.io/BSplineKit.jl/dev/collocation

Implements Greville sites by calculating collocation points as sliding averages over k-1 knots.

```julia
AvgKnots <: SelectionMethod
```

--------------------------------

### Defining Natural Boundary Conditions

Source: https://jipolanco.github.io/BSplineKit.jl/dev/boundary_conditions

Represents generalized natural boundary conditions for even-order splines.

```julia
Natural <: BoundaryCondition
```

--------------------------------

### SameAsKnots selection method

Source: https://jipolanco.github.io/BSplineKit.jl/dev/collocation

Selects collocation points to match B-spline knots, applicable when the basis length equals the number of knots.

```julia
SameAsKnots <: SelectionMethod
```