README.md

@@ -5,7 +5,7 @@ Modification of the classical vkoga code to solve PDEs in a greedy kernel way.
 
 ## Installation
 
-pip install git+https://gitlab.mathematik.uni-stuttgart.de/pub/ians-anm/pde-vkoga@v0.1.0
+pip install git+https://gitlab.mathematik.uni-stuttgart.de/pub/ians-anm/pde-vkoga@v0.1.1
 
 
 ## Usage

examples/example.py

@@ -6,7 +6,7 @@
 import numpy as np
 from matplotlib import pyplot as plt
 
-from vkoga_pde.kernels_PDE import Gaussian_laplace
+from vkoga_pde.kernels_PDE import Gaussian_laplace, wendland32_laplace, cubicMatern_laplace
 from vkoga_pde.vkoga_PDE import VKOGA_PDE
 
 np.random.seed(1)
@@ -28,17 +28,19 @@ f = lambda x: -6 * x[:, [0]] * x[:, [1]]
 
 # Define right hand side values
 y1 = f(X1)
+y1_sol = u(X1)
 y2 = u(X2)
 
 # Initialize and run models
-kernel = Gaussian_laplace(dim=2)
+# kernel = Gaussian_laplace(dim=2)
+kernel = cubicMatern_laplace(dim=2)
 maxIter = 200
 
 model_f = VKOGA_PDE(kernel=kernel, greedy_type='f_greedy')
 model_p = VKOGA_PDE(kernel=kernel, greedy_type='p_greedy')
 
-_ = model_f.fit(X1, y1, X2, y2, maxIter=maxIter)
-_ = model_p.fit(X1, y1, X2, y2, maxIter=maxIter)
+_ = model_f.fit(X1, y1, X2, y2, maxIter=maxIter, y1_sol=y1_sol)
+_ = model_p.fit(X1, y1, X2, y2, maxIter=maxIter, y1_sol=y1_sol)
 
 # Plot some results
 plt.figure(1)
@@ -48,24 +50,34 @@ plt.plot(np.array(model_p.train_hist['f'])**.5, 'b')
 plt.yscale('log')
 plt.xscale('log')
 plt.legend(['f-greedy', 'p-greedy'])
-plt.title('residuals')
+plt.title('L(u-s_n) residuals')
 plt.draw()
 
 plt.figure(2)
 plt.clf()
+plt.plot(np.array(model_f.train_hist['f sol']), 'r')
+plt.plot(np.array(model_p.train_hist['f sol']), 'b')
+plt.yscale('log')
+plt.xscale('log')
+plt.legend(['f-greedy', 'p-greedy'])
+plt.title('(u-s_n) residuals')
+plt.draw()
+
+plt.figure(3)
+plt.clf()
 plt.plot(model_f.ctrs_[model_f.ind_ctrs1, 0], model_f.ctrs_[model_f.ind_ctrs1, 1], 'rx')
 plt.plot(model_f.ctrs_[model_f.ind_ctrs2, 0], model_f.ctrs_[model_f.ind_ctrs2, 1], 'ro')
 plt.title('Selected centers of f-greedy')
 plt.draw()
 
-plt.figure(3)
+plt.figure(4)
 plt.clf()
 plt.plot(model_p.ctrs_[model_p.ind_ctrs1, 0], model_p.ctrs_[model_p.ind_ctrs1, 1], 'bx')
 plt.plot(model_p.ctrs_[model_p.ind_ctrs2, 0], model_p.ctrs_[model_p.ind_ctrs2, 1], 'bo')
 plt.title('Selected centers of P-greedy')
 plt.draw()
 
-plt.figure(4)
+plt.figure(5)
 plt.clf()
 plt.plot(X1[:, 0], X1[:, 1], 'm.')
 plt.plot(X2[:, 0], X2[:, 1], 'r.')

setup.py

@@ -2,7 +2,7 @@ from setuptools import setup
 
 setup(
     name='PDE-VKOGA',
-    version='0.1.0',    
+    version='0.1.1',    
     description='Python implementation of PDE-VKOGA to approximate the solution of linear PDEs.',
     url='https://gitlab.mathematik.uni-stuttgart.de/pub/ians-anm/pde-vkoga',
     author='Tizian Wenzel',

vkoga_pde/kernels_PDE.py

@@ -89,34 +89,130 @@ class L_RBF(Kernel):
         # Since the LaPlace is radial, the formula is the same as d2_eval
         return self.d2_eval(x, y)
 
-    def mixed_eval_L(self, x, ctrs, list_ind_interior):
-        # Kernel evaluation L^x k(x, y), whereby y1 := y[ind_interior, :] are centers in the interior,
-        # y2 := y[~ind_interior, :] are centers on the boundary.
-        # Therefore ddeval needs to be used with respect to y1, but deval with respect to y2.
+    def mixed_L_eval(self, x, ctrs, n_ctrs, list_ind_L, list_ind_k, list_ind_n):
+        # L evaluations in x of Riesz representers which are located in ctrs.
+        # list_ind_L, list_ind_k and list_ind_n indicate to which set of functionals the Riesz
+        # representer belong to:
+        # y1 := y[list_ind_L, :] corresponds to the interior, y2 := y[list_ind_k, :] to Dirichlet
+        # boundary centers and y3 := y[list_ind_n, :] to Neumann boundary centers.
+        # Different eval-methods need to be used with respect to y1, y2, y3.
 
-        list_not_ind_interior = [idx for idx in range(x.shape[0]) if idx not in list_ind_interior]
+        x = np.atleast_2d(x)
+        ctrs = np.atleast_2d(ctrs)
 
         array_result = np.zeros((x.shape[0], ctrs.shape[0]))
 
-        array_result[list_ind_interior, :] = self.dd_eval(x[list_ind_interior, :], ctrs)
-        array_result[list_not_ind_interior, :] = self.d1_eval(x[list_not_ind_interior, :], ctrs)
+        array_result[:, list_ind_L] = self.dd_eval(x, ctrs[list_ind_L])
+        array_result[:, list_ind_k] = self.d1_eval(x, ctrs[list_ind_k])
+        array_result[:, list_ind_n] = self.mixed_L_n(x, ctrs[list_ind_n], n_ctrs[list_ind_n])
 
         return array_result
 
-    def mixed_eval(self, x, ctrs, list_ind_interior):
-        # Kernel evaluation k(x, y), whereby y1 := y[ind_interior, :] are centers in the interior,
-        # y2 := y[~ind_interior, :] are centers on the boundary.
-        # Therefore deval needs to be used with respect to y1, but eval with respect to y2.
+    def mixed_k_eval(self, x, ctrs, n_ctrs, list_ind_L, list_ind_k, list_ind_n):
+        # Point evaluations in x of Riesz representers which are located in ctrs.
+        # list_ind_L, list_ind_k and list_ind_n indicate to which set of functionals the Riesz
+        # representer belong to:
+        # y1 := y[list_ind_L, :] corresponds to the interior, y2 := y[list_ind_k, :] to Dirichlet
+        # boundary centers and y3 := y[list_ind_n, :] to Neumann boundary centers.
+        # Different eval-methods need to be used with respect to y1, y2, y3.
 
-        list_not_ind_interior = [idx for idx in range(x.shape[0]) if idx not in list_ind_interior]
+        x = np.atleast_2d(x)
+        ctrs = np.atleast_2d(ctrs)
 
         array_result = np.zeros((x.shape[0], ctrs.shape[0]))
 
-        array_result[list_ind_interior, :] = self.d2_eval(x[list_ind_interior, :], ctrs)
-        array_result[list_not_ind_interior, :] = self.eval(x[list_not_ind_interior, :], ctrs)
+        array_result[:, list_ind_L] = self.d2_eval(x, ctrs[list_ind_L])
+        array_result[:, list_ind_k] = self.eval(x, ctrs[list_ind_k])
+        array_result[:, list_ind_n] = self.mixed_k_n(x, ctrs[list_ind_n], n_ctrs[list_ind_n])
 
         return array_result
 
+    def mixed_n_eval(self, x, n_x, ctrs, n_ctrs, list_ind_L, list_ind_k, list_ind_n):
+        # Neumann derivative evaluations in x with normal vector n_x of Riesz
+        # representers which are located in ctrs. n_ctrs contains normal vectors
+        # to those centers listed in list_ind_n - if this list is empty, then the
+        # array n_ctrs should be a zero array.
+        # list_ind_L, list_ind_k and list_ind_n indicate to which set of functionals the Riesz
+        # representer belong to:
+        # y1 := y[list_ind_L, :] corresponds to the interior, y2 := y[list_ind_k, :] to Dirichlet
+        # boundary centers and y3 := y[list_ind_n, :] to Neumann boundary centers.
+        # Different eval-methods need to be used with respect to y1, y2, y3.
+
+        x = np.atleast_2d(x)
+        ctrs = np.atleast_2d(ctrs)
+
+        array_result = np.zeros((x.shape[0], ctrs.shape[0]))
+
+        array_result[:, list_ind_L] = self.mixed_L_n(ctrs[list_ind_L], x, n_x).T
+        array_result[:, list_ind_k] = self.mixed_k_n(ctrs[list_ind_k], x, n_x).T
+        array_result[:, list_ind_n] = self.mixed_n_n(
+            x, ctrs[list_ind_n], n_x, n_ctrs[list_ind_n])
+
+        return array_result
+
+    # Now implement some functions related to Neumann boundary conditions
+    def nn_diagonal(self, X):
+        # Evaluation of < d/dn k(*, x), d/dn k(*, y) > for x = y,
+        # i.e. dot product between two Neumann boundary functionals.
+        # Value is independent of direction on n, therefore it is not provided/required.
+
+        X = np.atleast_2d(X)
+
+        return - np.ones(X.shape[0]) * self.rbf1_divided_by_r(self.ep, 0)
+
+    def mixed_L_n(self, x, y, n_y):
+        # Evaluation of < L^{2}k(*, x), d/dn k(*, y) >,  i.e. dot product of L evaluation and Neumann boundary
+        # See 14.04.22\C, especially \C4.
+        # n_y is the normal vectors for y.
+
+        x = np.atleast_2d(x)
+        y = np.atleast_2d(y)
+        n_y = np.atleast_2d(n_y)
+
+        array_dists = distance_matrix(x, y)
+
+        array_result1 = x @ n_y.T - np.sum(y * n_y, axis=1, keepdims=True).T
+        array_result2 = (self.dim - 1) * self.rbf_util_neuman(self.ep, array_dists) \
+                        + self.rbf3_divided_by_r(self.ep, array_dists)
+
+        return array_result1 * array_result2
+
+    def mixed_k_n(self, x, y, n_y):
+        # Evaluation of < k(*, x), d/dn k(*, y) >,  i.e. dot product of k evaluation and Neumann boundary
+        # See 14.04.22\C, especially \C2.
+        # n_y is the normal vectors for y.
+
+        x = np.atleast_2d(x)
+        y = np.atleast_2d(y)
+        n_y = np.atleast_2d(n_y)
+
+        array_dists = distance_matrix(x, y)
+
+        array_result1 = np.sum(y * n_y, axis=1, keepdims=True).T - x @ n_y.T
+        array_result2 = self.rbf1_divided_by_r(self.ep, array_dists)
+
+        return array_result1 * array_result2
+
+    def mixed_n_n(self, x, y, n_x, n_y):
+        # Evaluation of < d/dn k(*, x), d/dn k(*, y) >,  i.e. dot product between two Neumann boundary functionals
+        # See 14.04.22\C, especially \C2.
+        # n_y is the normal vectors for y, same for n_x.
+
+        x = np.atleast_2d(x)
+        y = np.atleast_2d(y)
+        n_x = np.atleast_2d(n_x)
+        n_y = np.atleast_2d(n_y)
+
+        array_dists = distance_matrix(x, y)
+
+        array_1 = - self.rbf1_divided_by_r(self.ep, array_dists) * (n_x @ n_y.T)
+
+        array_2 = - self.rbf_util_neuman(self.ep, array_dists) * \
+                  (x @ n_y.T - np.sum(y * n_y, axis=1, keepdims=True).T) \
+                  * (np.sum(x * n_x, axis=1, keepdims=True) - n_x @ y.T)
+
+        return array_1 + array_2
+
 
 # Implementation of concrete RBFs
 class Gaussian_laplace(L_RBF):
@@ -139,10 +235,15 @@ class Gaussian_laplace(L_RBF):
 
         self.rbf2 = lambda ep, r: (4 * r ** 2 - 2) * np.exp(-r**2)
 
+        self.rbf3_divided_by_r = lambda ep, r: -4 * (2 * r**2 - 3) * np.exp(-r**2)
+
         # Utility RBF for dd_eval: Double LaPlacian
         self.rbf_utility = lambda ep, r: \
                 np.exp(-r ** 2) * (16 * r ** 4 - (16 * self.dim + 32) * r ** 2 + (4 * self.dim ** 2 + 8 * self.dim))
 
+        # Utility RBF for neumann: See 14.04.22\C4, this is Phi'' / r^2 - Phi' / r^3
+        self.rbf_util_neuman = lambda ep, r: 4 * np.exp(-r ** 2)
+
 
 class wendland32_laplace(L_RBF):
     # Wendland kernel 3, 2: phi(r) = (1-r)_+^6 * (35r**2 + 18r + 3)
@@ -174,6 +275,10 @@ class wendland32_laplace(L_RBF):
                    + r**2*(10080*self.dim**2 + 60480*self.dim + 80640)
                    + r*(-6720*self.dim**2 - 26880*self.dim - 20160)) * (r < 1)
 
+        # Utility RBF for neumann: See 14.04.22\C4, this is Phi'' / r^2 - Phi' / r^3
+        self.rbf3_divided_by_r = lambda ep, r: np.zeros_like(r)     # ToDo: Not (yet) implemented
+        self.rbf_util_neuman = lambda ep, r: np.zeros_like(r)       # ToDo: Not (yet) implemented
+
 
 class wendland33_laplace(L_RBF):
     # Wendland kernel 3, 3: phi(r) = (1-r)_+^8 * (32r**3+ 25r^2 + 8r + 1)
@@ -197,6 +302,8 @@ class wendland33_laplace(L_RBF):
 
         self.rbf2 = lambda ep, r: np.maximum(1 - r, 0) ** 6 * 22 * (160 * r ** 3 + 15 * r ** 2 - 6 * r - 1)
 
+        self.rbf3_divided_by_r = lambda ep, r: -np.maximum(1-r, 0) ** 5 * 1584 * (20 * r ** 2 - 5 * r - 1)
+
         # Utility RBF for dd_eval: Double LaPlacian
         self.rbf_utility = lambda ep, r: \
             (528 * self.dim ** 2 + 1056 * self.dim + r ** 7 * (3168 * self.dim ** 2 + 50688 * self.dim + 199584)
@@ -206,6 +313,8 @@ class wendland33_laplace(L_RBF):
              + r ** 3 * (36960 * self.dim ** 2 + 295680 * self.dim + 554400)
              + r ** 2 * (-11088 * self.dim ** 2 - 66528 * self.dim - 88704)) * (r < 1)
 
+        self.rbf_util_neuman = lambda ep, r: 528 * np.maximum(1 - r, 0) ** 6 * (6*r + 1)
+
 
 
 
@@ -232,6 +341,221 @@ class quadrMatern_laplace(L_RBF):
         self.rbf_utility = lambda ep, r: \
             np.exp(-r) * (r**2 - (2 * self.dim + 3) * r + self.dim ** 2 + 2 * self.dim)
 
+        # Utility RBF for neumann: See 14.04.22\C4, this is Phi'' / r^2 - Phi' / r^3
+        self.rbf3_divided_by_r = lambda ep, r: np.zeros_like(r)  # ToDo: Not (yet) implemented
+        self.rbf_util_neuman = lambda ep, r: np.zeros_like(r)       # ToDo: Not (yet) implemented
+
+
+class cubicMatern_laplace(L_RBF):
+    # Quadratic matern kernel phi(r) = (3 + 3r + r^2)*exp(-r),
+    # used in conjunction with L = -LaPlace
+
+    def __init__(self, dim, shift=0, ep=1):
+        super().__init__()
+        self.dim = dim
+        self.shift = shift
+        self.ep = 1  # ep ToDo: Implement general case!
+        self.name = 'cubic Matern'
+
+        # Implement not only the rbf, but also its derivatives
+        self.rbf = lambda ep, r: np.exp(-r) * (15 + 15 * r + 6 * r ** 2 + r ** 3)
+
+        self.rbf1 = lambda ep, r: np.exp(-r) * (-r) * (r ** 2 + 3 * r + 3)
+        self.rbf1_divided_by_r = lambda ep, r: np.exp(-r) * (-1) * (r ** 2 + 3 * r + 3)
+
+        self.rbf2 = lambda ep, r: np.exp(-r) * (r ** 3 - 3 * r - 3)
+
+        # Utility RBF for dd_eval: Double LaPlacian
+        self.rbf_utility = lambda ep, r: \
+            np.exp(-r) * ((self.dim ** 2 + 2 * self.dim)
+                          + (self.dim ** 2 + 2 * self.dim) * r
+                          - (4 + 2 * self.dim) * r ** 2 + r ** 3)
+
+
+        # Utility RBF for neumann: See 14.04.22\C4, this is Phi'' / r^2 - Phi' / r^3
+        self.rbf3_divided_by_r = lambda ep, r: np.zeros_like(r)  # ToDo: Not (yet) implemented
+        self.rbf_util_neuman = lambda ep, r: np.zeros_like(r)       # ToDo: Not (yet) implemented
+
+
+
+
+class L_RBF_derivative(Kernel):
+    # L = d/dx, i.e. only first derivative
+
+    @abstractmethod
+    def __init__(self):
+        super(L_RBF_derivative, self).__init__()
+
+    def eval(self, x, y):
+        return self.rbf(self.ep, distance_matrix(np.atleast_2d(x), np.atleast_2d(y)))
+
+    def diagonal(self, X):
+        return np.ones(X.shape[0]) * self.rbf(self.ep, 0)
+
+    def __str__(self):
+        return self.name + ' [gamma = %2.2e]' % self.ep
+
+    def set_params(self, par):
+        self.ep = par
+
+    def dd_diagonal(self, X):
+        # Evaluation of (L^x - shift) (L^y - shift) k(x,y) for x=y: Likely some constant
+
+        return np.ones(X.shape[0]) * self.dd_eval(np.zeros((1, 1)), np.zeros((1, 1))).item()
+
+    def dd_eval(self, x, y):
+        # Evaluation of (L^x - shift) (L^y - shift) k(x,y)
+
+        array_dists = distance_matrix(np.atleast_2d(x), np.atleast_2d(y))
+
+        return - self.rbf2(self.ep, array_dists)
+
+    def d2_eval(self, x, y):
+        # Evaluation of (L^z - shift * Id) k(*,z)|_{x,y}
+        # Required exactly like this for "compute the nth basis", see also 14.12.21\A1
+
+        array_dists = distance_matrix(np.atleast_2d(x), np.atleast_2d(y))
+        array_x_minus_y = np.atleast_2d(x) - np.atleast_2d(y).transpose()
+
+        # Introduce minus since derivative wrt y instead of x
+        array_result = - np.sign(array_x_minus_y) * self.rbf1(self.ep, array_dists)
+
+        return array_result
+
+    def d1_eval(self, x, y):
+        # Evaluation of (L - shift * Id) k(*,y)|_{x}
+        # Required exactly like this for "compute the nth basis", see also 14.12.21\A1
+
+        # Need a minus here since first derivative depends (in sign) on whether x or y
+        return - self.d2_eval(x, y)
+
+    def mixed_L_eval(self, x, ctrs, n_ctrs, list_ind_L, list_ind_k, list_ind_n):
+        # L evaluations in x of Riesz representers which are located in ctrs.
+        # list_ind_L, list_ind_k and list_ind_n indicate to which set of functionals the Riesz
+        # representer belong to:
+        # y1 := y[list_ind_L, :] corresponds to the interior, y2 := y[list_ind_k, :] to Dirichlet
+        # boundary centers and y3 := y[list_ind_n, :] to Neumann boundary centers.
+        # Different eval-methods need to be used with respect to y1, y2, y3.
+
+        x = np.atleast_2d(x)
+        ctrs = np.atleast_2d(ctrs)
+
+        array_result = np.zeros((x.shape[0], ctrs.shape[0]))
+
+        array_result[:, list_ind_L] = self.dd_eval(x, ctrs[list_ind_L])
+        array_result[:, list_ind_k] = self.d1_eval(x, ctrs[list_ind_k])
+        array_result[:, list_ind_n] = self.mixed_L_n(x, ctrs[list_ind_n], n_ctrs[list_ind_n])
+
+        return array_result
+
+    def mixed_k_eval(self, x, ctrs, n_ctrs, list_ind_L, list_ind_k, list_ind_n):
+        # Point evaluations in x of Riesz representers which are located in ctrs.
+        # list_ind_L, list_ind_k and list_ind_n indicate to which set of functionals the Riesz
+        # representer belong to:
+        # y1 := y[list_ind_L, :] corresponds to the interior, y2 := y[list_ind_k, :] to Dirichlet
+        # boundary centers and y3 := y[list_ind_n, :] to Neumann boundary centers.
+        # Different eval-methods need to be used with respect to y1, y2, y3.
+
+        x = np.atleast_2d(x)
+        ctrs = np.atleast_2d(ctrs)
+
+        array_result = np.zeros((x.shape[0], ctrs.shape[0]))
+
+        array_result[:, list_ind_L] = self.d2_eval(x, ctrs[list_ind_L])
+        array_result[:, list_ind_k] = self.eval(x, ctrs[list_ind_k])
+        array_result[:, list_ind_n] = self.mixed_k_n(x, ctrs[list_ind_n], n_ctrs[list_ind_n])
+
+        return array_result
+
+    def mixed_n_eval(self, x, n_x, ctrs, n_ctrs, list_ind_L, list_ind_k, list_ind_n):
+        # not required therefore just return zeros
+
+        return np.zeros((x.shape[0], ctrs.shape[0]))
+
+    def nn_diagonal(self, X):
+        # not required therefore just return zeros
+
+        return np.zeros(X.shape[0])
+
+    def mixed_L_n(self, x, y, n_y):
+        # not required therefore just return zeros
+
+        return np.zeros((x.shape[0], y.shape[0]))
+
+    def mixed_k_n(self, x, y, n_y):
+        # not required, therefore just return zeros
+
+        return np.zeros((x.shape[0], y.shape[0]))
+
+    def mixed_n_n(self, x, y, n_x, n_y):
+        # not required, therefore just return zeros
+
+        return np.zeros((x.shape[0], y.shape[0]))
+
+
+class Gaussian_derivative(L_RBF_derivative):
+    # # Gaussian kernel phi(r) = exp(-r^2),
+    # used in conjunction with L = d/dx in 1D
+
+    def __init__(self, dim, shift=0, ep=1):
+        super().__init__()
+        self.dim = dim
+        self.shift = shift
+        self.ep = 1     # ToDo: Generalize this one!!!
+        self.name = 'gauss'
+
+        # Implement not only the rbf, but also its derivatives
+        self.rbf = lambda ep, r: np.exp(-r**2)
+        self.rbf1 = lambda ep, r: -2 * r * np.exp(-r**2)
+        self.rbf2 = lambda ep, r: (4 * r ** 2 - 2) * np.exp(-r**2)
+
+
+class LinearMatern_derivative(L_RBF_derivative):
+    # used in conjunction with L = d/dx in 1D
+
+    def __init__(self, dim, shift=0, ep=1):
+        super().__init__()
+        self.dim = dim
+        self.shift = shift
+        self.ep = 1     # ToDo: Generalize this one!!!
+        self.name = 'quadratic Matern'
+
+        # Implement not only the rbf, but also its derivatives
+        self.rbf = lambda ep, r:  np.exp(-r) * (1 + r)
+        self.rbf1 = lambda ep, r: - np.exp(-r) * r
+        self.rbf2 = lambda ep, r: np.exp(-r) * (r - 1)
+
+class QuadraticMatern_derivative(L_RBF_derivative):
+    # used in conjunction with L = d/dx in 1D
+
+    def __init__(self, dim, shift=0, ep=1):
+        super().__init__()
+        self.dim = dim
+        self.shift = shift
+        self.ep = 1     # ToDo: Generalize this one!!!
+        self.name = 'quadratic Matern'
+
+        # Implement not only the rbf, but also its derivatives
+        self.rbf = lambda ep, r: np.exp(-r) * (3 + 3 * r + r ** 2)
+        self.rbf1 = lambda ep, r: np.exp(-r) * (-r) * (r+1)
+        self.rbf2 = lambda ep, r: np.exp(-r) * (r ** 2 - r - 1)
+
+
+
+class CubicMatern_derivative(L_RBF_derivative):
+    # used in conjunction with L = d/dx in 1D
+
+    def __init__(self, dim, shift=0, ep=1):
+        super().__init__()
+        self.dim = dim
+        self.shift = shift
+        self.ep = 1     # ToDo: Generalize this one!!!
+        self.name = 'cubic Matern'
+
+        # Implement not only the rbf, but also its derivatives
+        self.rbf = lambda ep, r: np.exp(-r) * (15 + 15 * r + 6 * r ** 2 + r ** 3)
+        self.rbf1 = lambda ep, r: np.exp(-r) * (-r) * (r ** 2 + 3 * r + 3)
+        self.rbf2 = lambda ep, r: np.exp(-r) * (r ** 3 - 3 * r - 3)