Merge pull request #215 from SURGroup/feature/unify_kernels

unify kernels

src/UQpy/dimension_reduction/diffusion_maps/DiffusionMaps.py

@@ -16,7 +16,7 @@
 from typing import Annotated, Union
 from beartype.vale import Is
 from UQpy.utilities.ValidationTypes import Numpy2DFloatArray, NumpyFloatArray
-from UQpy.utilities.kernels.baseclass import Kernel
+from UQpy.utilities.kernels.baseclass.Kernel import Kernel
 
 
 class DiffusionMaps:
@@ -68,7 +68,7 @@ def __init__(
         if kernel_matrix is not None:
             self.kernel_matrix = kernel_matrix
         elif data is not None and kernel is not None:
-            kernel.calculate_kernel_matrix(points=data)
+            kernel.calculate_kernel_matrix(x=data, s=data)
             self.kernel_matrix = kernel.kernel_matrix
         else:
             raise ValueError("Either `kernel_matrix` or both `data` and `kernel` must be provided")

src/UQpy/dimension_reduction/grassmann_manifold/projections/baseclass/GrassmannProjection.py

@@ -1,7 +1,5 @@
 from abc import ABC, abstractmethod
 
-from UQpy.utilities.kernels.baseclass.Kernel import Kernel
-
 
 class GrassmannProjection(ABC):
     """

src/UQpy/surrogates/gaussian_process/GaussianProcessRegression.py

@@ -4,13 +4,11 @@
 
 from beartype import beartype
 
-
 from UQpy.utilities.Utilities import process_random_state
 from UQpy.surrogates.baseclass.Surrogate import Surrogate
 from UQpy.utilities.ValidationTypes import RandomStateType
-from UQpy.surrogates.gaussian_process.kernels.baseclass.Kernel import Kernel
+from UQpy.utilities.kernels.baseclass.Kernel import Kernel
 from UQpy.surrogates.gaussian_process.constraints.baseclass.Constraints import ConstraintsGPR
-from UQpy.surrogates.gaussian_process.regression_models.baseclass.Regression import Regression
 
 
 class GaussianProcessRegression(Surrogate):
@@ -204,14 +202,18 @@ def fit(
                 raise NotImplementedError("Maximum likelihood estimator failed: Choose different starting point or "
                                           "increase nopt")
             t = np.argmin(fun_value)
-            self.hyperparameters = 10**minimizer[t, :]
+            self.hyperparameters = 10 ** minimizer[t, :]
 
         # Updated Correlation matrix corresponding to MLE estimates of hyperparameters
         if self.noise:
-            self.K = self.kernel.c(x=s_, s=s_, params=self.hyperparameters[:-1]) + \
-                     np.eye(nsamples)*(self.hyperparameters[-1])**2
+            self.kernel.kernel_parameter = self.hyperparameters[:-2]
+            sigma = self.hyperparameters[-2]
+            self.K = sigma ** 2 * self.kernel.calculate_kernel_matrix(x=s_, s=s_) + \
+                     np.eye(nsamples) * (self.hyperparameters[-1]) ** 2
         else:
-            self.K = self.kernel.c(x=s_, s=s_, params=self.hyperparameters)
+            self.kernel.kernel_parameter = self.hyperparameters[:-1]
+            sigma = self.hyperparameters[-1]
+            self.K = sigma ** 2 * self.kernel.calculate_kernel_matrix(x=s_, s=s_)
 
         self.cc = cholesky(self.K + 1e-10 * np.eye(nsamples), lower=True)
         self.alpha_ = cho_solve((self.cc, True), y_)
@@ -264,7 +266,9 @@ def predict(self, points, return_std: bool = False, hyperparameters: list = None
 
         if hyperparameters is not None:
             # This is used for MLE constraints, if constraints call 'predict' method.
-            K = self.kernel.c(x=s_, s=s_, params=kernelparameters) + \
+            self.kernel.kernel_parameter = kernelparameters[:-1]
+            sigma = kernelparameters[-1]
+            K = sigma ** 2 * self.kernel.calculate_kernel_matrix(x=s_, s=s_) + \
                 np.eye(self.samples.shape[0]) * (noise_std ** 2)
             cc = np.linalg.cholesky(K + 1e-10 * np.eye(self.samples.shape[0]))
             mu = 0
@@ -278,7 +282,7 @@ def predict(self, points, return_std: bool = False, hyperparameters: list = None
                 # Design parameters (beta: regression coefficient)
                 beta = np.linalg.solve(g_, np.matmul(np.transpose(q_), y_dash))
                 mu = np.einsum("ij,jk->ik", self.F, beta)
-            alpha_ = cho_solve((cc, True), y_-mu)
+            alpha_ = cho_solve((cc, True), y_ - mu)
         else:
             cc, alpha_ = self.cc, self.alpha_
 
@@ -290,15 +294,20 @@ def predict(self, points, return_std: bool = False, hyperparameters: list = None
             else:
                 mu1 = np.einsum("ij,jk->ik", fx, self.beta)
 
-        k = self.kernel.c(x=x_, s=s_, params=kernelparameters)
+        self.kernel.kernel_parameter = kernelparameters[:-1]
+        sigma = kernelparameters[-1]
+
+        k = sigma**2*self.kernel.calculate_kernel_matrix(x=x_, s=s_)
         y = mu1 + k @ alpha_
         if self.normalize:
             y = self.value_mean + y * self.value_std
         if x_.shape[1] == 1:
             y = y.flatten()
 
         if return_std:
-            k1 = self.kernel.c(x=x_, s=x_, params=kernelparameters)
+            self.kernel.kernel_parameter = kernelparameters[:-1]
+            sigma = kernelparameters[-1]
+            k1 = sigma**2*self.kernel.calculate_kernel_matrix(x=x_, s=x_)
             var = (k1 - k @ cho_solve((cc, True), k.T)).diagonal()
             mse = np.sqrt(var)
             if self.normalize:
@@ -326,9 +335,13 @@ def log_likelihood(p0, k_, s, y, ind_noise, fx_):
         m = s.shape[0]
 
         if ind_noise:
-            k__ = k_.c(x=s, s=s, params=10 ** p0[:-1]) + np.eye(m) * (10 ** p0[-1]) ** 2
+            k_.kernel_parameter = 10 ** p0[:-2]
+            sigma = 10 ** p0[-2]
+            k__ = sigma ** 2 * k_.calculate_kernel_matrix(x=s, s=s) + np.eye(m) * (10 ** p0[-1]) ** 2
         else:
-            k__ = k_.c(x=s, s=s, params=10 ** p0)
+            k_.kernel_parameter = 10 ** p0[:-1]
+            sigma = 10 ** p0[-1]
+            k__ = sigma ** 2 * k_.calculate_kernel_matrix(x=s, s=s)
         cc = cholesky(k__ + 1e-10 * np.eye(m), lower=True)
 
         mu = 0
@@ -343,7 +356,7 @@ def log_likelihood(p0, k_, s, y, ind_noise, fx_):
             beta = np.linalg.solve(g_, np.matmul(np.transpose(q_), y_dash))
             mu = np.einsum("ij,jk->ik", fx_, beta)
 
-        term1 = (y-mu).T @ (cho_solve((cc, True), y-mu))
+        term1 = (y - mu).T @ (cho_solve((cc, True), y - mu))
         term2 = 2 * np.sum(np.log(np.abs(np.diag(cc))))
 
         return 0.5 * (term1 + term2 + m * np.log(2 * np.pi))[0, 0]

src/UQpy/surrogates/gaussian_process/__init__.py

@@ -1,5 +1,4 @@
 from UQpy.surrogates.gaussian_process.GaussianProcessRegression import GaussianProcessRegression
 
-from UQpy.surrogates.gaussian_process.kernels import *
 from UQpy.surrogates.gaussian_process.regression_models import *
 from UQpy.surrogates.gaussian_process.constraints import *

src/UQpy/utilities/kernels/GaussianKernel.py

@@ -1,10 +1,14 @@
+import itertools
+from typing import Tuple
+
 import numpy as np
 import scipy
-import scipy.spatial.distance as sd
+from scipy.spatial.distance import cdist
 
 from UQpy.utilities.ValidationTypes import RandomStateType, Numpy2DFloatArray
 from UQpy.utilities.kernels import EuclideanKernel
 from scipy.spatial.distance import pdist
+import scipy.spatial.distance as sd
 
 
 class GaussianKernel(EuclideanKernel):
@@ -15,26 +19,26 @@ class GaussianKernel(EuclideanKernel):
         k(x_j, x_i) = \exp[-(x_j - xj)^2/4\epsilon]
 
     """
-    def __init__(self, epsilon: float = 1.0):
+    def __init__(self, kernel_parameter: float = 1.0):
         """
         :param epsilon: Scale parameter of the Gaussian kernel
         """
-        super().__init__()
-        self.epsilon = epsilon
+        super().__init__(kernel_parameter=kernel_parameter)
 
-    def kernel_entry(self, xi: Numpy2DFloatArray, xj: Numpy2DFloatArray):
-        """
-        Given two points, this method computes the Gaussian kernel value between those two points
+    def calculate_kernel_matrix(self, x, s):
+        product = [self.element_wise_operation(point_pair)
+                   for point_pair in list(itertools.product(x, s))]
+        self.kernel_matrix = np.array(product).reshape(len(x), len(s))
+        return self.kernel_matrix
+
+    def element_wise_operation(self, xi_j: Tuple) -> float:
+        xi, xj = xi_j
 
-        :param xi: First point.
-        :param xj: Second point.
-        :return: Float representing the kernel entry.
-        """
         if len(xi.shape) == 1:
             d = pdist(np.array([xi, xj]), "sqeuclidean")
         else:
-            d = np.linalg.norm(xi-xj, 'fro') ** 2
-        return np.exp(-d / (2*self.epsilon**2))
+            d = np.linalg.norm(xi - xj, 'fro') ** 2
+        return np.exp(-d / (2 * self.kernel_parameter ** 2))
 
     def optimize_parameters(self, data: np.ndarray, tolerance: float,
                             n_nearest_neighbors: int,

src/UQpy/utilities/kernels/__init__.py

@@ -1,6 +1,7 @@
 from UQpy.utilities.kernels.baseclass import *
+from UQpy.utilities.kernels.euclidean_kernels import *
+from UQpy.utilities.kernels.grassmannian_kernels import *
 
-from UQpy.utilities.kernels.BinetCauchyKernel import BinetCauchyKernel
+from UQpy.utilities.kernels.grassmannian_kernels.BinetCauchyKernel import BinetCauchyKernel
 from UQpy.utilities.kernels.GaussianKernel import GaussianKernel
-from UQpy.utilities.kernels.ProjectionKernel import ProjectionKernel
 

src/UQpy/utilities/kernels/baseclass/EuclideanKernel.py

@@ -1,39 +1,7 @@
-import itertools
-from abc import ABC, abstractmethod
-from typing import Union
-import scipy.spatial.distance as sd
-import numpy as np
+from abc import ABC
 
-from UQpy.utilities import GrassmannPoint
-from UQpy.utilities.ValidationTypes import NumpyFloatArray, Numpy2DFloatArray
 from UQpy.utilities.kernels.baseclass.Kernel import Kernel
 
 
 class EuclideanKernel(Kernel, ABC):
     """This is a blueprint for Euclidean kernels implemented in the :py:mod:`kernels` module ."""
-
-    def __init__(self):
-        super().__init__()
-
-    def calculate_kernel_matrix(self, points: Numpy2DFloatArray):
-        """
-        Using the kernel-specific :py:meth:`.Kernel.kernel_entry` method, this function assembles the kernel matrix.
-
-        :param points: Set of data points in the Euclidean space
-
-        """
-        nargs = len(points)
-        indices = range(nargs)
-        pairs = list(itertools.combinations_with_replacement(indices, 2))
-        kernel = np.zeros((nargs, nargs))
-        for id_pair in range(np.shape(pairs)[0]):
-            i = pairs[id_pair][0]
-            j = pairs[id_pair][1]
-
-            xi = points[i]
-            xj = points[j]
-
-            kernel[i, j] = self.kernel_entry(xi, xj)
-            kernel[j, i] = kernel[i, j]
-
-        self.kernel_matrix = kernel