From 73ef24a9222390925cfc2e64743c9b4a555f5825 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
Date: Wed, 30 Oct 2024 14:09:01 +0100
Subject: [PATCH] Add interface for general bases (#100)

* add initial interface for general bases

* Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* fix least squares in tutorial

* fix SpatialDiscretization

* format

* fix docstring rendering

* fix TemporalDiscretization

* fix unit tests

* add unit tests

* Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>

* fix docstring rendering

* use StandardBasis in some examples

---------

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
---
 docs/src/ref.md                               |   7 ++
 docs/src/tutorial_noisy_data.md               |   2 +-
 .../interpolation/interpolation_2d_sphere.jl  |   3 +-
 examples/interpolation/least_squares_2d.jl    |   2 +-
 src/KernelInterpolation.jl                    |   4 +-
 src/basis.jl                                  |  73 ++++++++++++
 src/discretization.jl                         |  49 +++++---
 src/interpolation.jl                          | 105 +++++++++++-------
 src/kernel_matrices.jl                        |  91 +++++++++------
 src/visualization.jl                          |  19 +++-
 test/test_unit.jl                             |  28 ++++-
 test/test_visualization.jl                    |   4 +
 12 files changed, 287 insertions(+), 100 deletions(-)
 create mode 100644 src/basis.jl

diff --git a/docs/src/ref.md b/docs/src/ref.md
index 883add03..60701d89 100644
--- a/docs/src/ref.md
+++ b/docs/src/ref.md
@@ -23,6 +23,13 @@ Modules = [KernelInterpolation]
 Pages = ["nodes.jl"]
 ```
 
+## Bases
+
+```@autodocs
+Modules = [KernelInterpolation]
+Pages = ["basis.jl"]
+```
+
 ## Interpolation
 
 ```@autodocs
diff --git a/docs/src/tutorial_noisy_data.md b/docs/src/tutorial_noisy_data.md
index 966d1fc7..f6b15a64 100644
--- a/docs/src/tutorial_noisy_data.md
+++ b/docs/src/tutorial_noisy_data.md
@@ -128,7 +128,7 @@ approximation and the polynomial augmentation is not changed. In KernelInterpola
 ```@example noisy-itp
 M = 81
 centers = random_hypercube(M; dim = 2)
-ls = interpolate(nodeset, centers, values_noisy, kernel)
+ls = interpolate(centers, nodeset, values_noisy, kernel)
 ```
 
 We plot the least-squares approximation and, again, see a better fit to the underlying target function.
diff --git a/examples/interpolation/interpolation_2d_sphere.jl b/examples/interpolation/interpolation_2d_sphere.jl
index 4a341bf3..88f3257c 100644
--- a/examples/interpolation/interpolation_2d_sphere.jl
+++ b/examples/interpolation/interpolation_2d_sphere.jl
@@ -13,7 +13,8 @@ merge!(nodeset, nodeset_boundary)
 values = f.(nodeset)
 
 kernel = InverseMultiquadricKernel{dim(nodeset)}()
-itp = interpolate(nodeset, values, kernel)
+basis = StandardBasis(nodeset, kernel)
+itp = interpolate(basis, values)
 
 N = 500
 many_nodes = random_hypersphere(N, r, center)
diff --git a/examples/interpolation/least_squares_2d.jl b/examples/interpolation/least_squares_2d.jl
index cee18fcf..31ab0746 100644
--- a/examples/interpolation/least_squares_2d.jl
+++ b/examples/interpolation/least_squares_2d.jl
@@ -16,7 +16,7 @@ M = 81
 centers = random_hypercube(M; dim = 2)
 
 kernel = ThinPlateSplineKernel{dim(nodeset)}()
-ls = interpolate(nodeset, centers, values, kernel)
+ls = interpolate(StandardBasis(centers, kernel), values, nodeset)
 itp = interpolate(nodeset, values, kernel)
 
 N = 40
diff --git a/src/KernelInterpolation.jl b/src/KernelInterpolation.jl
index 06940920..3d10499d 100644
--- a/src/KernelInterpolation.jl
+++ b/src/KernelInterpolation.jl
@@ -31,10 +31,11 @@ using WriteVTK: WriteVTK, vtk_grid, paraview_collection, MeshCell, VTKCellTypes,
 
 include("kernels/kernels.jl")
 include("nodes.jl")
+include("basis.jl")
+include("regularization.jl")
 include("differential_operators.jl")
 include("equations.jl")
 include("kernel_matrices.jl")
-include("regularization.jl")
 include("interpolation.jl")
 include("discretization.jl")
 include("callbacks_step/callbacks_step.jl")
@@ -48,6 +49,7 @@ export GaussKernel, MultiquadricKernel, InverseMultiquadricKernel,
        RadialCharacteristicKernel, MaternKernel, Matern12Kernel, Matern32Kernel,
        Matern52Kernel, Matern72Kernel, RieszKernel,
        TransformationKernel, ProductKernel, SumKernel
+export StandardBasis
 export phi, Phi, order
 export PartialDerivative, Gradient, Laplacian, EllipticOperator
 export PoissonEquation, EllipticEquation, AdvectionEquation, HeatEquation,
diff --git a/src/basis.jl b/src/basis.jl
new file mode 100644
index 00000000..693b0618
--- /dev/null
+++ b/src/basis.jl
@@ -0,0 +1,73 @@
+"""
+    AbstractBasis
+
+Abstract type for a basis of a kernel function space. Every basis represents a
+set of functions, which can be obtained by indexing the basis object. Every basis
+object holds a kernel function and a [`NodeSet`](@ref) of centers and potentially
+more fields depending on the concrete basis type.
+"""
+abstract type AbstractBasis end
+
+function (basis::AbstractBasis)(x)
+    return [basis[i](x) for i in eachindex(basis)]
+end
+
+"""
+    interpolation_kernel(basis)
+
+Return the kernel from a basis.
+"""
+interpolation_kernel(basis::AbstractBasis) = basis.kernel
+
+"""
+	centers(basis)
+
+Return the centers from a basis object.
+"""
+centers(basis::AbstractBasis) = basis.centers
+
+"""
+    order(basis)
+
+Return the order ``m`` of the polynomial, which is needed by this `basis` for
+the interpolation, i.e., the polynomial degree plus 1. If ``m = 0``,
+no polynomial is added.
+"""
+order(basis::AbstractBasis) = order(interpolation_kernel(basis))
+dim(basis::AbstractBasis) = dim(basis.centers)
+Base.length(basis::AbstractBasis) = length(centers(basis))
+Base.eachindex(basis::AbstractBasis) = Base.OneTo(length(basis))
+function Base.iterate(basis::AbstractBasis, state = 1)
+    state > length(basis) ? nothing : (basis[state], state + 1)
+end
+Base.collect(basis::AbstractBasis) = Function[basis[i] for i in 1:length(basis)]
+
+function Base.show(io::IO, basis::AbstractBasis)
+    return print(io,
+                 "$(nameof(typeof(basis))) with $(length(centers(basis))) centers and kernel $(interpolation_kernel(basis)).")
+end
+
+@doc raw"""
+    StandardBasis(centers, kernel)
+
+The standard basis for a function space defined by a kernel and a [`NodeSet`](@ref) of `centers`.
+The basis functions are given by
+
+```math
+    b_j(x) = K(x, x_j)
+```
+
+where `K` is the kernel and `x_j` are the nodes in `centers`.
+"""
+struct StandardBasis{Kernel} <: AbstractBasis
+    centers::NodeSet
+    kernel::Kernel
+    function StandardBasis(centers::NodeSet, kernel::Kernel) where {Kernel}
+        if dim(kernel) != dim(centers)
+            throw(DimensionMismatch("The dimension of the kernel and the centers must be the same"))
+        end
+        new{typeof(kernel)}(centers, kernel)
+    end
+end
+
+Base.getindex(basis::StandardBasis, i) = x -> basis.kernel(x, centers(basis)[i])
diff --git a/src/discretization.jl b/src/discretization.jl
index ca102747..17dfef1f 100644
--- a/src/discretization.jl
+++ b/src/discretization.jl
@@ -1,4 +1,5 @@
 """
+    SpatialDiscretization(equations, nodeset_inner, boundary_condition, nodeset_boundary, basis)
     SpatialDiscretization(equations, nodeset_inner, boundary_condition, nodeset_boundary,
                           [centers,] kernel = GaussKernel{dim(nodeset_inner)}())
 
@@ -10,27 +11,35 @@ is a function describing the Dirichlet boundary conditions. The `centers` are th
 See also [`Semidiscretization`](@ref), [`solve_stationary`](@ref).
 """
 struct SpatialDiscretization{Dim, RealT, Equations, BoundaryCondition,
-                             Kernel <: AbstractKernel{Dim}}
+                             Basis <: AbstractBasis}
     equations::Equations
     nodeset_inner::NodeSet{Dim, RealT}
     boundary_condition::BoundaryCondition
     nodeset_boundary::NodeSet{Dim, RealT}
-    centers::NodeSet{Dim, RealT}
-    kernel::Kernel
+    basis::Basis
 
     function SpatialDiscretization(equations, nodeset_inner::NodeSet{Dim, RealT},
                                    boundary_condition,
                                    nodeset_boundary::NodeSet{Dim, RealT},
-                                   centers::NodeSet{Dim, RealT},
-                                   kernel = GaussKernel{Dim}()) where {Dim,
-                                                                       RealT}
+                                   basis::AbstractBasis) where {Dim,
+                                                                RealT}
         new{Dim, RealT, typeof(equations), typeof(boundary_condition),
-            typeof(kernel)}(equations, nodeset_inner,
-                            boundary_condition, nodeset_boundary,
-                            centers, kernel)
+            typeof(basis)}(equations, nodeset_inner,
+                           boundary_condition, nodeset_boundary,
+                           basis)
     end
 end
 
+function SpatialDiscretization(equations, nodeset_inner::NodeSet{Dim, RealT},
+                               boundary_condition,
+                               nodeset_boundary::NodeSet{Dim, RealT},
+                               centers::NodeSet{Dim, RealT},
+                               kernel = GaussKernel{Dim}()) where {Dim,
+                                                                   RealT}
+    SpatialDiscretization(equations, nodeset_inner, boundary_condition,
+                          nodeset_boundary, StandardBasis(centers, kernel))
+end
+
 function SpatialDiscretization(equations, nodeset_inner::NodeSet{Dim, RealT},
                                boundary_condition,
                                nodeset_boundary::NodeSet{Dim, RealT},
@@ -42,8 +51,11 @@ function SpatialDiscretization(equations, nodeset_inner::NodeSet{Dim, RealT},
 end
 
 function Base.show(io::IO, sd::SpatialDiscretization)
+    N_i = length(sd.nodeset_inner)
+    N_b = length(sd.nodeset_boundary)
+    k = interpolation_kernel(sd.basis)
     print(io,
-          "SpatialDiscretization with $(dim(sd)) dimensions, $(length(sd.nodeset_inner)) inner nodes, $(length(sd.nodeset_boundary)) boundary nodes, and kernel $(sd.kernel)")
+          "SpatialDiscretization with $(dim(sd)) dimensions, $N_i inner nodes, $N_b boundary nodes, and kernel $k")
 end
 
 dim(::SpatialDiscretization{Dim}) where {Dim} = Dim
@@ -60,7 +72,8 @@ function solve_stationary(spatial_discretization::SpatialDiscretization{Dim, Rea
                                                                                             Dim,
                                                                                             RealT
                                                                                             }
-    @unpack equations, nodeset_inner, boundary_condition, nodeset_boundary, centers, kernel = spatial_discretization
+    @unpack equations, nodeset_inner, boundary_condition, nodeset_boundary, basis = spatial_discretization
+    @unpack centers, kernel = basis
 
     system_matrix = pde_boundary_matrix(equations, nodeset_inner, nodeset_boundary, centers,
                                         kernel)
@@ -71,7 +84,7 @@ function solve_stationary(spatial_discretization::SpatialDiscretization{Dim, Rea
     xx = polyvars(Dim)
     ps = monomials(xx, 0:-1)
     nodeset = merge(nodeset_inner, nodeset_boundary)
-    return Interpolation(kernel, nodeset, centers, c, system_matrix,
+    return Interpolation(basis, nodeset, c, system_matrix,
                          ps, xx)
 end
 
@@ -96,10 +109,11 @@ end
 
 function Semidiscretization(spatial_discretization::SpatialDiscretization,
                             initial_condition)
-    @unpack equations, nodeset_inner, boundary_condition, nodeset_boundary, centers, kernel = spatial_discretization
+    @unpack equations, nodeset_inner, boundary_condition, nodeset_boundary, basis = spatial_discretization
+    @unpack centers, kernel = basis
     @assert length(centers)==length(nodeset_inner) + length(nodeset_boundary) "The number of centers must be equal to the number of inner and boundary nodes."
-    k_matrix_inner = kernel_matrix(nodeset_inner, centers, kernel)
-    k_matrix_boundary = kernel_matrix(nodeset_boundary, centers, kernel)
+    k_matrix_inner = kernel_matrix(centers, nodeset_inner, kernel)
+    k_matrix_boundary = kernel_matrix(centers, nodeset_boundary, kernel)
     # whole kernel matrix is not needed for rhs, but for initial condition
     k_matrix = [k_matrix_inner
                 k_matrix_boundary]
@@ -133,8 +147,11 @@ function Semidiscretization(equations, nodeset_inner::NodeSet{Dim, RealT},
 end
 
 function Base.show(io::IO, semi::Semidiscretization)
+    N_i = length(semi.spatial_discretization.nodeset_inner)
+    N_b = length(semi.spatial_discretization.nodeset_boundary)
+    k = interpolation_kernel(semi.spatial_discretization.basis)
     print(io,
-          "Semidiscretization with $(dim(semi)) dimensions, $(length(semi.spatial_discretization.nodeset_inner)) inner nodes, $(length(semi.spatial_discretization.nodeset_boundary)) boundary nodes, and kernel $(semi.spatial_discretization.kernel)")
+          "Semidiscretization with $(dim(semi)) dimensions, $N_i inner nodes, $N_b boundary nodes, and kernel $k")
 end
 
 dim(semi::Semidiscretization) = dim(semi.spatial_discretization)
diff --git a/src/interpolation.jl b/src/interpolation.jl
index 8079f10a..48a04ca4 100644
--- a/src/interpolation.jl
+++ b/src/interpolation.jl
@@ -1,25 +1,25 @@
-abstract type AbstractInterpolation{Kernel, Dim, RealT} end
+abstract type AbstractInterpolation{Basis, Dim, RealT} end
 
 @doc raw"""
     Interpolation
 
 Interpolation object that can be evaluated at a node and represents a kernel interpolation of the form
 ```math
-    s(x) = \sum_{j = 1}^N c_jK(x, x_j) + \sum_{k = 1}^Q d_kp_k(x),
+    s(x) = \sum_{j = 1}^N c_jb_j(x) + \sum_{k = 1}^Q d_kp_k(x),
 ```
-where ``x_j`` are the nodes in the nodeset and ``s(x)`` the interpolant satisfying ``s(x_j) = f(x_j)``, where
-``f(x_j)`` are given by `values` in [`interpolate`](@ref) and ``p_k`` is a basis of the `Q`-dimensional space
-of multivariate polynomials of order [`order`](@ref). The additional conditions
+where ``b_j`` are the basis functions and ``p_k`` is a basis of the `Q`-dimensional space of multivariate
+polynomials of order [`order`](@ref). The additional conditions
 ```math
     \sum_{j = 1}^N c_jp_k(x_j) = 0, \quad k = 1,\ldots, Q
 ```
 are enforced.
+
+See also [`interpolate`](@ref).
 """
-struct Interpolation{Kernel, Dim, RealT, A, Monomials, PolyVars} <:
-       AbstractInterpolation{Kernel, Dim, RealT}
-    kernel::Kernel
+struct Interpolation{Basis, Dim, RealT, A, Monomials, PolyVars} <:
+       AbstractInterpolation{Basis, Dim, RealT}
+    basis::Basis
     nodeset::NodeSet{Dim, RealT}
-    centers::NodeSet{Dim, RealT}
     c::Vector{RealT}
     system_matrix::A
     ps::Monomials
@@ -36,14 +36,21 @@ end
 
 Return the dimension of the input variables of the interpolation.
 """
-dim(::Interpolation{Kernel, Dim, RealT, A}) where {Kernel, Dim, RealT, A} = Dim
+dim(::Interpolation{Basis, Dim}) where {Basis, Dim} = Dim
+
+"""
+    basis(itp)
+
+Return the basis from an interpolation object.
+"""
+basis(itp::Interpolation) = itp.basis
 
 """
     interpolation_kernel(itp)
 
 Return the kernel from an interpolation object.
 """
-interpolation_kernel(itp::AbstractInterpolation) = itp.kernel
+interpolation_kernel(itp::AbstractInterpolation) = interpolation_kernel(basis(itp))
 
 """
 	nodeset(itp)
@@ -52,6 +59,13 @@ Return the node set from an interpolation object.
 """
 nodeset(itp::AbstractInterpolation) = itp.nodeset
 
+"""
+    centers(itp::Interpolation)
+
+Return the centers from the basis of an interpolation object.
+"""
+centers(itp::Interpolation) = centers(basis(itp))
+
 """
     coefficients(itp::Interpolation)
 
@@ -70,7 +84,7 @@ interpolant.
 
 See also [`coefficients`](@ref) and [`polynomial_coefficients`](@ref).
 """
-kernel_coefficients(itp::Interpolation) = itp.c[eachindex(itp.centers)]
+kernel_coefficients(itp::Interpolation) = itp.c[eachindex(centers(itp))]
 
 """
     polynomial_coefficients(itp::Interpolation)
@@ -80,7 +94,7 @@ interpolant.
 
 See also [`coefficients`](@ref) and [`kernel_coefficients`](@ref).
 """
-polynomial_coefficients(itp::Interpolation) = itp.c[(length(itp.centers) + 1):end]
+polynomial_coefficients(itp::Interpolation) = itp.c[(length(centers(itp)) + 1):end]
 
 """
     polynomial_basis(itp::Interpolation)
@@ -121,59 +135,68 @@ of known values. The exact form of ``A`` differs depending on which method is us
 system_matrix(itp::Interpolation) = itp.system_matrix
 
 @doc raw"""
-    interpolate(nodeset, [centers,] values, kernel = GaussKernel{dim(nodeset)}();
+    interpolate(basis, values, nodeset = centers(basis); m = order(basis),
+                regularization = NoRegularization())
+    interpolate(centers, [nodeset,] values, kernel = GaussKernel{dim(nodeset)}();
                 m = order(kernel), regularization = NoRegularization())
 
 Interpolate the `values` evaluated at the nodes in the `nodeset` to a function using the kernel `kernel`
 and polynomials up to a order `m` (i.e. degree - 1), i.e., determine the coefficients ``c_j`` and ``d_k`` in the expansion
 ```math
-    s(x) = \sum_{j = 1}^N c_jK(x, x_j) + \sum_{k = 1}^Q d_kp_k(x),
+    s(x) = \sum_{j = 1}^N c_jb_j(x) + \sum_{k = 1}^Q d_kp_k(x),
 ```
-where ``x_j`` are the nodes in the nodeset and ``s(x)`` the interpolant ``s(x_j) = f(x_j)``, where ``f(x_j)``
-are given by `values` and ``p_k`` is a basis of the ``Q``-dimensional space of multivariate polynomials with
-maximum degree of `m - 1`. If `m = 0`, no polynomial is added. The additional conditions
+where ``b_j`` are the basis functions in the `basis` and ``s(x)`` the interpolant ``s(x_j) = f(x_j)``, where ``f(x_j)``
+are given by `values`, ``x_j`` are the nodes in the `nodeset`, and ``p_k`` is a basis of the ``Q``-dimensional
+space of multivariate polynomials with maximum degree of `m - 1`. If `m = 0`, no polynomial is added.
+The additional conditions
 ```math
     \sum_{j = 1}^N c_jp_k(x_j) = 0, \quad k = 1,\ldots, Q = \begin{pmatrix}m - 1 + d\\d\end{pmatrix}
 ```
 are enforced. Returns an [`Interpolation`](@ref) object.
 
-If `centers` is provided, the interpolant is a least squares approximation with the centers used for the basis.
-Otherwise, `centers` is set to `nodeset`.
+If `nodeset` is provided, the interpolant is a least squares approximation with a different set of nodes as the centers
+used for the basis.
+Otherwise, `nodeset` is set to `centers(basis)` or `centers`.
 
 A regularization can be applied to the kernel matrix using the `regularization` argument, cf. [`regularize!`](@ref).
 """
-function interpolate(nodeset::NodeSet{Dim, RealT}, centers::NodeSet{Dim, RealT},
-                     values::Vector{RealT}, kernel = GaussKernel{Dim}();
-                     m = order(kernel),
+function interpolate(basis::AbstractBasis, values::Vector{RealT},
+                     nodeset::NodeSet{Dim} = centers(basis);
+                     m = order(basis),
                      regularization = NoRegularization()) where {Dim, RealT}
-    @assert dim(kernel) == Dim
+    @assert dim(basis) == Dim
     n = length(nodeset)
     @assert length(values) == n
     xx = polyvars(Dim)
     ps = monomials(xx, 0:(m - 1))
     q = length(ps)
 
-    if nodeset == centers
-        system_matrix = interpolation_matrix(nodeset, kernel, ps, regularization)
+    if nodeset == centers(basis)
+        system_matrix = interpolation_matrix(basis, ps, regularization)
     else
-        system_matrix = least_squares_matrix(nodeset, centers, kernel, ps, regularization)
+        system_matrix = least_squares_matrix(basis, nodeset, ps, regularization)
     end
     b = [values; zeros(q)]
     c = system_matrix \ b
-    return Interpolation(kernel, nodeset, centers, c, system_matrix, ps, xx)
+    return Interpolation(basis, nodeset, c, system_matrix, ps, xx)
+end
+function interpolate(centers::NodeSet{Dim, RealT}, nodeset::NodeSet{Dim, RealT},
+                     values::AbstractVector{RealT}, kernel = GaussKernel{Dim}();
+                     kwargs...) where {Dim, RealT}
+    interpolate(StandardBasis(centers, kernel), values, nodeset; kwargs...)
 end
 
-function interpolate(nodeset::NodeSet{Dim, RealT},
-                     values::Vector{RealT}, kernel = GaussKernel{Dim}();
+function interpolate(centers::NodeSet{Dim, RealT},
+                     values::AbstractVector{RealT}, kernel = GaussKernel{Dim}();
                      kwargs...) where {Dim, RealT}
-    interpolate(nodeset, nodeset, values, kernel; kwargs...)
+    interpolate(StandardBasis(centers, kernel), values; kwargs...)
 end
 
 # Evaluate interpolant
 function (itp::Interpolation)(x)
     s = 0
     kernel = interpolation_kernel(itp)
-    xis = itp.centers
+    xis = centers(itp)
     c = kernel_coefficients(itp)
     d = polynomial_coefficients(itp)
     ps = polynomial_basis(itp)
@@ -197,7 +220,7 @@ end
 function (diff_op_or_pde::Union{AbstractDifferentialOperator, AbstractStationaryEquation})(itp::Interpolation,
                                                                                            x)
     kernel = interpolation_kernel(itp)
-    xis = itp.centers
+    xis = centers(itp)
     c = kernel_coefficients(itp)
     s = zero(eltype(x))
     for j in eachindex(c)
@@ -208,7 +231,7 @@ end
 
 function (g::Gradient)(itp::Interpolation, x)
     kernel = interpolation_kernel(itp)
-    xis = itp.centers
+    xis = centers(itp)
     c = kernel_coefficients(itp)
     s = zero(x)
     for j in eachindex(c)
@@ -236,8 +259,8 @@ function kernel_inner_product(itp1, itp2)
     @assert kernel == interpolation_kernel(itp2)
     c_f = kernel_coefficients(itp1)
     c_g = kernel_coefficients(itp2)
-    xs = itp1.centers
-    xis = itp2.centers
+    xs = centers(itp1)
+    xis = centers(itp2)
     s = 0
     for i in eachindex(c_f)
         for j in eachindex(c_g)
@@ -275,20 +298,24 @@ end
 function Base.show(io::IO, titp::TemporalInterpolation)
     sd = titp.ode_sol.prob.p.spatial_discretization
     tspan = titp.ode_sol.prob.tspan
+    N_i = length(sd.nodeset_inner)
+    N_b = length(sd.nodeset_boundary)
+    k = interpolation_kernel(sd.basis)
     return print(io,
-                 "Temporal interpolation with $(length(sd.nodeset_inner)) inner nodes, $(length(sd.nodeset_boundary)) boundary nodes, kernel $(sd.kernel), and time span $tspan")
+                 "Temporal interpolation with $N_i inner nodes, $N_b boundary nodes, kernel $k, and time span $tspan")
 end
 
 function (titp::TemporalInterpolation)(t)
     ode_sol = titp.ode_sol
     semi = ode_sol.prob.p
-    @unpack kernel, nodeset_inner, boundary_condition, nodeset_boundary, centers = semi.spatial_discretization
+    @unpack nodeset_inner, boundary_condition, nodeset_boundary, basis = semi.spatial_discretization
+    @unpack centers, kernel = basis
     c = ode_sol(t)
     # Do not support additional polynomial basis for now
     xx = polyvars(dim(semi))
     ps = monomials(xx, 0:-1)
     nodeset = merge(nodeset_inner, nodeset_boundary)
-    itp = Interpolation(kernel, nodeset, centers, c,
+    itp = Interpolation(StandardBasis(centers, kernel), nodeset, c,
                         semi.cache.mass_matrix, ps, xx)
     return itp
 end
diff --git a/src/kernel_matrices.jl b/src/kernel_matrices.jl
index df4cdacd..4f2a8cd6 100644
--- a/src/kernel_matrices.jl
+++ b/src/kernel_matrices.jl
@@ -1,35 +1,37 @@
 @doc raw"""
-    kernel_matrix(nodeset1, nodeset2, kernel)
+    kernel_matrix(basis, nodeset = centers(basis))
+    kernel_matrix(nodeset1[, nodeset2], kernel)
 
-Return the kernel matrix for the nodeset and kernel. The kernel matrix is defined as
+Return the kernel matrix for the `nodes` and `kernel`. The kernel matrix is defined as
 ```math
-    A_{ij} = K(x_i, \xi_j),
+    A_{ij} = b_j(x_i),
 ```
-where ``x_i`` are the nodes in the `nodeset1`, ``\xi_j`` are the nodes in the `nodeset2`,
-and ``K`` the `kernel`.
+where ``b_i`` are the basis function in the `basis` and `x_i` are the nodes in the `nodeset`.
+If two nodesets and a `kernel` are given, the kernel matrix is computed for the [`StandardBasis`](@ref) meaning
+```math
+    A_{ij} = K(\xi_j, x_i),
+```
+where ``\xi_j`` are the nodes/centers in `nodeset1`, ``x_i`` are the nodes in `nodeset2`, and `K` is the `kernel`.
+If `nodeset2` is not given, it defaults to `nodeset1`.
 """
-function kernel_matrix(nodeset1, nodeset2, kernel)
-    n = length(nodeset1)
-    m = length(nodeset2)
-    A = Matrix{eltype(nodeset1)}(undef, n, m)
+function kernel_matrix(basis::AbstractBasis, nodeset::NodeSet = centers(basis))
+    n = length(nodeset)
+    m = length(basis)
+    A = Matrix{eltype(nodeset)}(undef, n, m)
     for i in 1:n
         for j in 1:m
-            A[i, j] = kernel(nodeset1[i], nodeset2[j])
+            A[i, j] = basis[j](nodeset[i])
         end
     end
     return A
 end
 
-"""
-    kernel_matrix(nodeset, kernel)
+function kernel_matrix(nodeset1::NodeSet{Dim}, nodeset2::NodeSet{Dim},
+                       kernel::AbstractKernel{Dim}) where {Dim}
+    kernel_matrix(StandardBasis(nodeset1, kernel), nodeset2)
+end
 
-Return the kernel matrix for the nodeset and kernel. The kernel matrix is defined as
-```math
-    A_{ij} = K(x_i, x_j),
-```
-where ``x_i`` are the nodes in the `nodeset` and ``K`` the `kernel`.
-"""
-function kernel_matrix(nodeset, kernel)
+function kernel_matrix(nodeset::NodeSet, kernel::AbstractKernel)
     kernel_matrix(nodeset, nodeset, kernel)
 end
 
@@ -42,7 +44,7 @@ Return the polynomial matrix for the nodeset and polynomials. The polynomial mat
 ```
 where ``x_i`` are the nodes in the `nodeset` and ``p_j`` the polynomials.
 """
-function polynomial_matrix(nodeset, ps)
+function polynomial_matrix(nodeset::NodeSet, ps)
     n = length(nodeset)
     q = length(ps)
     A = Matrix{eltype(nodeset)}(undef, n, q)
@@ -55,48 +57,65 @@ function polynomial_matrix(nodeset, ps)
     return A
 end
 
-"""
-    interpolation_matrix(nodeset, kernel, ps, regularization)
+@doc raw"""
+    interpolation_matrix(centers, kernel, ps, regularization = NoRegularization())
+    interpolation_matrix(basis, ps, regularization)
 
-Return the interpolation matrix for the `nodeset`, `kernel`, polynomials `ps`, and `regularization`.
+Return the interpolation matrix for the `basis`, polynomials `ps`, and `regularization`.
 The interpolation matrix is defined as
 ```math
     A = \begin{pmatrix}K & P\\P^T & 0\end{pmatrix},
 ```
-where ``K`` is the [`regularize!`](@ref)d [`kernel_matrix`](@ref) and ``P`` the [`polynomial_matrix`](@ref)`.
+where ``K`` is the [`regularize!`](@ref)d [`kernel_matrix`](@ref) and ``P`` the [`polynomial_matrix`](@ref).
+If a node set of `centers` and a `kernel` are given, the interpolation matrix is computed for the [`StandardBasis`](@ref).
 """
-function interpolation_matrix(nodeset, kernel, ps, regularization)
+function interpolation_matrix(basis::AbstractBasis, ps,
+                              regularization::AbstractRegularization = NoRegularization())
     q = length(ps)
-    k_matrix = kernel_matrix(nodeset, kernel)
+    k_matrix = kernel_matrix(basis)
     regularize!(k_matrix, regularization)
-    p_matrix = polynomial_matrix(nodeset, ps)
+    p_matrix = polynomial_matrix(centers(basis), ps)
     system_matrix = [k_matrix p_matrix
                      p_matrix' zeros(q, q)]
     return Symmetric(system_matrix)
 end
 
-"""
-    least_squares_matrix(nodeset, centers, kernel, ps, regularization)
+function interpolation_matrix(centers::NodeSet, kernel::AbstractKernel, ps,
+                              regularization::AbstractRegularization = NoRegularization())
+    interpolation_matrix(StandardBasis(centers, kernel), ps, regularization)
+end
 
-Return the least squares matrix for the `nodeset`, `centers`, `kernel`, polynomials `ps`, and `regularization`.
+@doc raw"""
+    least_squares_matrix(basis, nodeset, ps, regularization = NoRegularization())
+    least_squares_matrix(centers, nodeset, kernel, ps, regularization = NoRegularization())
+
+Return the least squares matrix for the `basis`, `nodeset`, polynomials `ps`, and `regularization`.
 The least squares matrix is defined as
 ```math
-    A = \begin{pmatrix}K & P_1\\P_2' & 0\end{pmatrix},
+    A = \begin{pmatrix}K & P_1\\P_2^T & 0\end{pmatrix},
 ```
-where ``K`` is the [`regularize!`](@ref)d [`kernel_matrix`](@ref), ``P_1`` the [`polynomial_matrix`](@ref)`
+where ``K`` is the [`regularize!`](@ref)d [`kernel_matrix`](@ref), ``P_1`` the [`polynomial_matrix`](@ref)
 for the `nodeset` and ``P_2`` the [`polynomial_matrix`](@ref)` for the `centers`.
+If a `nodeset` and `kernel` are given, the least squares matrix is computed for the [`StandardBasis`](@ref).
 """
-function least_squares_matrix(nodeset, centers, kernel, ps, regularization)
+function least_squares_matrix(basis::AbstractBasis, nodeset::NodeSet, ps,
+                              regularization::AbstractRegularization = NoRegularization())
     q = length(ps)
-    k_matrix = kernel_matrix(nodeset, centers, kernel)
+    k_matrix = kernel_matrix(basis, nodeset)
     regularize!(k_matrix, regularization)
     p_matrix1 = polynomial_matrix(nodeset, ps)
-    p_matrix2 = polynomial_matrix(centers, ps)
+    p_matrix2 = polynomial_matrix(centers(basis), ps)
     system_matrix = [k_matrix p_matrix1
                      p_matrix2' zeros(q, q)]
     return system_matrix
 end
 
+function least_squares_matrix(centers::NodeSet, nodeset::NodeSet, kernel::AbstractKernel,
+                              ps,
+                              regularization::AbstractRegularization = NoRegularization())
+    least_squares_matrix(StandardBasis(centers, kernel), nodeset, ps, regularization)
+end
+
 @doc raw"""
     pde_matrix(diff_op_or_pde, nodeset1, nodeset2, kernel)
 
@@ -142,7 +161,7 @@ See also [`pde_matrix`](@ref) and [`kernel_matrix`](@ref).
 function pde_boundary_matrix(diff_op_or_pde, nodeset_inner, nodeset_boundary, centers,
                              kernel)
     pd_matrix = pde_matrix(diff_op_or_pde, nodeset_inner, centers, kernel)
-    b_matrix = kernel_matrix(nodeset_boundary, centers, kernel)
+    b_matrix = kernel_matrix(centers, nodeset_boundary, kernel)
     return [pd_matrix
             b_matrix]
 end
diff --git a/src/visualization.jl b/src/visualization.jl
index 412622f1..3608fafc 100644
--- a/src/visualization.jl
+++ b/src/visualization.jl
@@ -106,13 +106,13 @@ end
     elseif dim(nodeset) == 2
         if training_nodes
             @series begin
-                x = values_along_dim(itp.centers, 1)
-                y = values_along_dim(itp.centers, 2)
+                x = values_along_dim(centers(itp), 1)
+                y = values_along_dim(centers(itp), 2)
                 seriestype := :scatter
                 markershape --> :star
                 markersize --> 10
                 label --> "training nodes"
-                x, y, itp.(itp.centers)
+                x, y, itp.(centers(itp))
             end
         end
         @series begin
@@ -177,3 +177,16 @@ end
 @recipe function f(nodeset::NodeSet, f::Function)
     nodeset, f.(nodeset)
 end
+
+@recipe function f(basis::AbstractBasis, nodeset::NodeSet = centers(basis))
+    if dim(nodeset) != dim(basis)
+        throw(DimensionMismatch("The dimension of the node set and the basis must be the same"))
+    end
+    for i in 1:length(basis)
+        @series begin
+            label --> nothing
+            seriestype --> :surface
+            nodeset, basis[i]
+        end
+    end
+end
diff --git a/test/test_unit.jl b/test/test_unit.jl
index b9c3d9b9..4da1fcb6 100644
--- a/test/test_unit.jl
+++ b/test/test_unit.jl
@@ -609,6 +609,26 @@ end
     @test PointSet(NodeSet(ps)) == ps
 end
 
+@testitem "Basis" setup=[Setup, AdditionalImports] begin
+    nodeset = NodeSet([0.0 0.0
+                       1.0 0.0
+                       0.0 1.0
+                       1.0 1.0])
+    kernel = GaussKernel{dim(nodeset)}(shape_parameter = 0.5)
+    basis = @test_nowarn StandardBasis(nodeset, kernel)
+    @test_throws DimensionMismatch StandardBasis(nodeset,
+                                                 GaussKernel{1}(shape_parameter = 0.5))
+    @test_nowarn println(basis)
+    @test_nowarn display(basis)
+    A = kernel_matrix(basis)
+    @test isapprox(stack(basis.(nodeset)), A)
+    @test isapprox(stack(basis.(nodeset)), kernel.(distance_matrix(nodeset, nodeset)))
+    basis_functions = collect(basis)
+    for (i, b) in enumerate(basis)
+        @test b.(nodeset) == basis_functions[i].(nodeset)
+    end
+end
+
 @testitem "Interpolation" setup=[Setup, AdditionalImports] begin
     nodes = NodeSet([0.0 0.0
                      1.0 0.0
@@ -624,6 +644,8 @@ end
     @test nodeset(itp) == nodes
     @test dim(itp) == dim(kernel)
     @test dim(itp) == dim(nodes)
+    @test system_matrix(itp) ==
+          KernelInterpolation.interpolation_matrix(nodes, kernel, itp.ps)
     # Saving the interpolation and the function to a VTK file
     @test_nowarn vtk_save("itp", nodes, f, itp, ff; keys = ["f", "itp", "f2"])
     nodes2, point_data = @test_nowarn vtk_read("itp.vtu")
@@ -703,7 +725,7 @@ end
     centers = NodeSet([0.0 0.0
                        1.0 0.0
                        0.0 1.0])
-    itp = @test_nowarn interpolate(nodes, centers, ff, kernel)
+    itp = @test_nowarn interpolate(centers, nodes, ff, kernel)
     expected_coefficients = [
         0.0,
         0.0,
@@ -712,12 +734,14 @@ end
         1.0,
         1.0]
     coeffs = coefficients(itp)
+    @test system_matrix(itp) ==
+          KernelInterpolation.least_squares_matrix(centers, nodes, kernel, itp.ps)
     @test length(coeffs) == length(expected_coefficients)
     for i in eachindex(coeffs)
         @test isapprox(coeffs[i], expected_coefficients[i], atol = 1e-15)
     end
     @test order(itp) == order(kernel)
-    @test length(kernel_coefficients(itp)) == length(itp.centers)
+    @test length(kernel_coefficients(itp)) == length(KernelInterpolation.centers(itp))
     @test length(polynomial_coefficients(itp)) == order(itp) + 1
     @test length(polynomial_basis(itp)) ==
           binomial(order(itp) - 1 + dim(nodes), dim(nodes))
diff --git a/test/test_visualization.jl b/test/test_visualization.jl
index 66889e97..31728cef 100644
--- a/test/test_visualization.jl
+++ b/test/test_visualization.jl
@@ -26,6 +26,10 @@
                 nodes2d = homogeneous_hypercube(5; dim = 2)
                 @test_nowarn plot!(nodes2d)
             end
+            kernel1 = GaussKernel{dim}(shape_parameter = 0.5)
+            basis = StandardBasis(nodes, kernel1)
+            @test_nowarn plot(basis, nodes_fine)
+            @test_throws DimensionMismatch plot(basis, random_hypercube(5; dim = dim + 1))
         end
     end
 end