GPR: Unit tests covers GPR kernels.

docs/source/surrogates/gpr.rst

@@ -144,7 +144,7 @@ The Matern kernel takes the following form:
 
 .. math:: K(x, s, \theta) = \sigma^2 \frac{2^{1-\nu}}{\Gamma(\nu)} \bigg( \sqrt{2 \nu} \frac{d}{l} \bigg)^{\nu} K_v \bigg(\sqrt{2 \nu} \frac{d}{l} \bigg)
 
-where :math:`d = ||x-s||_2^{1/2}` is the euclidean distance and :math:`\theta` is consist of lengthscale (:math:`l`), process variance (:math:`\sigma^2`) and smoothing parameter (:math:`\nu`). Also, :math:`\Gamma` is tha gamma function and :math:`K_v` is the modified Bessel function. This kernel concides with exponential and RBF kernel for :math:`\nu=0.5` and :math:`\infty`, respectively.
+where :math:`d = ||x-s||_2^{1/2}` is the euclidean distance and :math:`\theta` is consist of lengthscale (:math:`l`), process variance (:math:`\sigma^2`) and smoothing parameter (:math:`\nu`). Also, :math:`\Gamma` is tha gamma function and :math:`K_v` is the modified Bessel function. This kernel concides with absolute exponential and RBF kernel for :math:`\nu=0.5` and :math:`\infty`, respectively.
 
 User-Defined Kernel
 """"""""""""""""""""""""""""

src/UQpy/surrogates/gaussian_process/GaussianProcessRegression.py

@@ -91,17 +91,20 @@ def __init__(
         self.mu = 0
 
         if bounds is None:
-            if self.optimizer._bounds is None:
-                bounds_ = [[10**-3, 10**3]] * self.hyperparameters.shape[0]
-                if self.noise:
-                    bounds_[-1] = [10**-10, 10**-1]
-                self.bounds = bounds_
+            if isinstance(self.optimizer, type(None)) or isinstance(self.optimizer._bounds, type(None)):
+                self._define_bounds()
             else:
                 self.bounds = self.optimizer._bounds
 
         self.jac = False
         self.random_state = process_random_state(random_state)
 
+    def _define_bounds(self):
+        bounds_ = [[10 ** -3, 10 ** 3]] * self.hyperparameters.shape[0]
+        if self.noise:
+            bounds_[-1] = [10 ** -10, 10 ** -1]
+        self.bounds = bounds_
+
     def fit(
             self,
             samples,

src/UQpy/surrogates/gaussian_process/kernels/Matern.py

@@ -13,20 +13,20 @@ def __init__(self, nu=1.5):
         """
         self.nu = nu
 
-    def c(self, x, s, params, dt=False, dx=False):
+    def c(self, x, s, params):
         l, sigma = params[:-1], params[-1]
         stack = cdist(x/l, s/l, metric='euclidean')
         if self.nu == 0.5:
-            return np.exp(-stack)
+            return sigma**2 * np.exp(-np.abs(stack))
         elif self.nu == 1.5:
-            return (1+np.sqrt(3)*stack)*np.exp(-np.sqrt(3)*stack)
+            return sigma**2 * (1+np.sqrt(3)*stack)*np.exp(-np.sqrt(3)*stack)
         elif self.nu == 2.5:
-            return (1+np.sqrt(5)*stack+5*(stack**2)/3)*np.exp(-np.sqrt(5)*stack)
+            return sigma**2 * (1+np.sqrt(5)*stack+5*(stack**2)/3)*np.exp(-np.sqrt(5)*stack)
         elif self.nu == np.inf:
-            return np.exp(-(stack**2)/2)
+            return sigma**2 * np.exp(-(stack**2)/2)
         else:
             stack *= np.sqrt(2*self.nu)
             tmp = 1/(gamma(self.nu)*(2**(self.nu-1)))
             tmp1 = stack ** self.nu
-            tmp2 = kv(tmp1)
+            tmp2 = kv(self.nu, stack)
             return sigma**2 * tmp * tmp1 * tmp2

src/UQpy/surrogates/gaussian_process/kernels/baseclass/Kernel.py

@@ -24,24 +24,24 @@ def check_samples_and_return_stack(x, s):
         ) - np.tile(s_, (np.size(x_, 0), 1, 1))
         return stack
 
-    @staticmethod
-    def derivatives(x_, s_, params):
-        stack = Kernel.check_samples_and_return_stack(x_, s_)
-        # Taking stack and creating array of all thetaj*dij
-        after_parameters = params * abs(stack)
-        # Create matrix of all ones to compare
-        comp_ones = np.ones((np.size(x_, 0), np.size(s_, 0), np.size(s_, 1)))
-        # zeta_matrix has all values min{1,theta*dij}
-        zeta_matrix_ = np.minimum(after_parameters, comp_ones)
-        # Copy zeta_matrix to another matrix that will used to find where derivative should be zero
-        indices = zeta_matrix_.copy()
-        # If value of min{1,theta*dij} is 1, the derivative should be 0.
-        # So, replace all values of 1 with 0, then perform the .astype(bool).astype(int)
-        # operation like in the linear example, so you end up with an array of 1's where
-        # the derivative should be caluclated and 0 where it should be zero
-        indices[indices == 1] = 0
-        # Create matrix of all |dij| (where non zero) to be used in calculation of dR/dtheta
-        dtheta_derivs_ = indices.astype(bool).astype(int) * abs(stack)
-        # Same as above, but for matrix of all thetaj where non-zero
-        dx_derivs_ = indices.astype(bool).astype(int) * params * np.sign(stack)
-        return zeta_matrix_, dtheta_derivs_, dx_derivs_
+    # @staticmethod
+    # def derivatives(x_, s_, params):
+    #     stack = Kernel.check_samples_and_return_stack(x_, s_)
+    #     # Taking stack and creating array of all thetaj*dij
+    #     after_parameters = params * abs(stack)
+    #     # Create matrix of all ones to compare
+    #     comp_ones = np.ones((np.size(x_, 0), np.size(s_, 0), np.size(s_, 1)))
+    #     # zeta_matrix has all values min{1,theta*dij}
+    #     zeta_matrix_ = np.minimum(after_parameters, comp_ones)
+    #     # Copy zeta_matrix to another matrix that will used to find where derivative should be zero
+    #     indices = zeta_matrix_.copy()
+    #     # If value of min{1,theta*dij} is 1, the derivative should be 0.
+    #     # So, replace all values of 1 with 0, then perform the .astype(bool).astype(int)
+    #     # operation like in the linear example, so you end up with an array of 1's where
+    #     # the derivative should be caluclated and 0 where it should be zero
+    #     indices[indices == 1] = 0
+    #     # Create matrix of all |dij| (where non zero) to be used in calculation of dR/dtheta
+    #     dtheta_derivs_ = indices.astype(bool).astype(int) * abs(stack)
+    #     # Same as above, but for matrix of all thetaj where non-zero
+    #     dx_derivs_ = indices.astype(bool).astype(int) * params * np.sign(stack)
+    #     return zeta_matrix_, dtheta_derivs_, dx_derivs_

tests/unit_tests/surrogates/test_gpr.py

@@ -15,48 +15,139 @@
 #
 samples = np.linspace(0, 5, 20).reshape(-1, 1)
 values = np.cos(samples)
-optimizer = MinimizeOptimizer(method="L-BFGS-B")
-gpr = GaussianProcessRegression(kernel=RBF(), optimizer=optimizer,
-                                hyperparameters=[0.14], optimize=False, random_state=1)
-gpr.fit(samples=samples, values=values, hyperparameters=[0.3])
+optimizer = MinimizeOptimizer(method="L-BFGS-B", bounds=[[0.1, 5], [0.1, 3]])
+gpr = GaussianProcessRegression(kernel=RBF(), hyperparameters=[2.852, 2.959], random_state=1)
+gpr.fit(samples=samples, values=values)
 
-optimizer = MinimizeOptimizer(method="L-BFGS-B")
-gpr2 = GaussianProcessRegression(kernel=Matern(nu=0.5), optimizer=optimizer,
-                                 hyperparameters=[0.3], bounds=[[0.01, 5]], optimize=False,
-                                 normalize=False, random_state=2)
+optimizer1 = MinimizeOptimizer(method="L-BFGS-B")
+gpr2 = GaussianProcessRegression(kernel=Matern(nu=0.5), hyperparameters=[10, 4], optimizer=optimizer1,
+                                 bounds=[[0.01, 100], [0.1, 10]], optimizations_number=10, random_state=2)
 gpr2.fit(samples=samples, values=values)
 
+points11 = np.array([[0], [1], [2]])
+points12 = np.array([[-1], [0], [-1]])
+hyperparameters1 = np.array([1.5, 10])
+
+points21 = np.array([[0, 0], [1, 1], [2, 2]])
+points22 = np.array([[-1, 0], [1, 2], [0, -2]])
+hyperparameters2 = np.array([3, 1, 2])
+rbf_kernel = RBF()
+matern_kernel_1_2 = Matern(nu=0.5)
+matern_kernel_3_2 = Matern(nu=1.5)
+matern_kernel_5_2 = Matern(nu=2.5)
+matern_kernel_inf = Matern(nu=np.inf)
+matern_kernel_2_1 = Matern(nu=2)
+
+
+def test_predict():
+    prediction = np.round(gpr.predict([[1], [np.pi/2], [np.pi]], True), 3)
+    expected_prediction = np.array([[0.54, 0., -1.], [0.,  0.,  0.]])
+    assert (expected_prediction == prediction).all()
+
+
+def test_predict1():
+    prediction = np.round(gpr2.predict([[1], [2*np.pi], [np.pi]], True), 3)
+    expected_prediction = np.array([[0.537,  0.238, -0.998], [0.077,  0.39,  0.046]])
+    assert (expected_prediction == prediction).all()
+
+
+def test_rbf():
+    """
+    Test RBF kernel for 1-d input model
+    """
+    cov = rbf_kernel.c(points11, points12, hyperparameters1)
+    covariance_matrix = np.round(cov, 3)
+    expected_covariance_matrix = np.array([[80.074, 100.,  80.074], [41.111,  80.074,  41.111],
+                                           [13.534,  41.111,  13.534]])
+    assert (expected_covariance_matrix == covariance_matrix).all()
+
+
+def test_rbf2():
+    """
+    Test RBF kernel for 2-d input model
+    """
+    cov = rbf_kernel.c(points21, points22, hyperparameters2)
+    covariance_matrix2 = np.round(cov, 3)
+    expected_covariance_matrix2 = np.array([[3.784e+00, 5.120e-01, 5.410e-01], [1.943e+00, 2.426e+00, 4.200e-02],
+                                            [3.280e-01, 3.784e+00, 1.000e-03]])
+    assert (expected_covariance_matrix2 == covariance_matrix2).all()
+
+
+def test_matern_inf():
+    """
+    Test Matern kernel (nu=inf) for 1-d input model, should concide with RBF kernel
+    """
+    cov = matern_kernel_inf.c(points11, points12, hyperparameters1)
+    covariance_matrix = np.round(cov, 3)
+    expected_covariance_matrix = np.array([[80.074, 100.,  80.074], [41.111,  80.074,  41.111],
+                                           [13.534,  41.111,  13.534]])
+    assert (expected_covariance_matrix == covariance_matrix).all()
+
+
+def tets_matern_inf2():
+    """
+    Test Matern kernel (nu=inf) for 2-d input model
+    """
+    cov = matern_kernel_inf.c(points21, points22, hyperparameters2)
+    covariance_matrix2 = np.round(cov, 3)
+    expected_covariance_matrix2 = np.array([[3.784e+00, 5.120e-01, 5.410e-01], [1.943e+00, 2.426e+00, 4.200e-02],
+                                            [3.280e-01, 3.784e+00, 1.000e-03]])
+    assert (expected_covariance_matrix2 == covariance_matrix2).all()
+
+
+def test_matern_1_2():
+    """
+    Test Matern kernel (nu=0.5) for 1-d input model, should concide with absolute exponential kernel
+    """
+    cov = matern_kernel_1_2.c(points11, points12, hyperparameters1)
+    covariance_matrix = np.round(cov, 3)
+    expected_covariance_matrix = np.array([[51.342, 100.,  51.342], [26.36,  51.342,  26.36],
+                                           [13.534,  26.36,  13.534]])
+    assert (expected_covariance_matrix == covariance_matrix).all()
+
+
+def test_matern_1_2_2():
+    """
+    Test Matern kernel (nu=0.5) for 2-d input model
+    """
+    cov = matern_kernel_1_2.c(points21, points22, hyperparameters2)
+    covariance_matrix2 = np.round(cov, 3)
+    expected_covariance_matrix2 = np.array([[2.866, 0.527, 0.541], [1.203, 1.472, 0.196], [0.428, 2.866, 0.069]])
+    assert (expected_covariance_matrix2 == covariance_matrix2).all()
+
+
+def test_matern_3_2():
+    """
+    Test Matern kernel (nu=1.5) for 1-d input model
+    """
+    cov = matern_kernel_3_2.c(points11, points12, hyperparameters1)
+    covariance_matrix = np.round(cov, 3)
+    expected_covariance_matrix = np.array([[67.906, 100.,  67.906], [32.869,  67.906,  32.869],
+                                           [13.973,  32.869,  13.973]])
+    assert (expected_covariance_matrix == covariance_matrix).all()
+
+
+def test_matern_5_2():
+    """
+    Test Matern kernel (nu=2.5) for 1-d input model
+    """
+    cov = matern_kernel_5_2.c(points11, points12, hyperparameters1)
+    covariance_matrix = np.round(cov, 3)
+    expected_covariance_matrix = np.array([[72.776, 100.,  72.776], [35.222,  72.776,  35.222],
+                                           [13.866,  35.222,  13.866]])
+    assert (expected_covariance_matrix == covariance_matrix).all()
+
+
+def test_matern_2_1():
+    """
+    Test Matern kernel (nu=2) for 1-d input model
+    """
+    cov = matern_kernel_2_1.c(points11, points12, hyperparameters1)
+    covariance_matrix = np.round(cov, 3)
+    expected_covariance_matrix = np.array([[70.894, np.nan,  70.894], [34.25, 70.894, 34.25],
+                                           [13.921, 34.25, 13.921]])
+    assert np.array_equal(expected_covariance_matrix, covariance_matrix, equal_nan=True)
 
-# # Using the in-built linear regression model as a function
-# # linear_regression_model = Kriging(regression_model=Linear(), correlation_model=Gaussian(), optimizer=optimizer,
-# #                                   correlation_model_parameters=[1]).regression_model
-# # optimizer = MinimizeOptimizer(method="L-BFGS-B")
-# rbf_kernel = GaussianProcessRegression(correlation_model=RBF(), optimizer=optimizer,
-#                                       correlation_model_parameters=[1]).kernel
-#
-# optimizer = MinimizeOptimizer(method="L-BFGS-B")
-# gpr3 = GaussianProcessRegression(correlation_model=rbf_kernel, optimizer=optimizer,
-#                                 correlation_model_parameters=[1], optimize=False, normalize=False, random_state=0)
-# gpr3.fit(samples=samples, values=values)
-#
-#
-# def test_predict():
-#     prediction = np.round(gpr.predict([[1], [np.pi/2], [np.pi]], True), 3)
-#     expected_prediction = np.array([[0.54,  0., -1.], [0.,  0.,  0.]])
-#     assert (expected_prediction == prediction).all()
-#
-#
-# def test_predict1():
-#     prediction = np.round(gpr2.predict([[1], [2*np.pi], [np.pi]], True), 3)
-#     expected_prediction = np.array([[0.54,  1.009, -1.], [0.,  0.031,  0.]])
-#     assert (expected_prediction == prediction).all()
-#
-#
-# def test_predict2():
-#     prediction = np.round(gpr3.predict([[1], [np.pi/2], [np.pi]]), 3)
-#     expected_prediction = np.array([[0.54, -0., -1.]])
-#     assert (expected_prediction == prediction).all()
-#
 #
 # def test_jacobian():
 #     jacobian = np.round(gpr.jacobian([[np.pi], [np.pi/2]]), 3)