Refactors Dimension Reduction kernels to match GP kernels signatures

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:

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/__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,7 +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

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

@@ -1,5 +1,6 @@
 import itertools
-from abc import ABC
+from abc import ABC, abstractmethod
+from typing import Union, Tuple
 
 import numpy as np
 
@@ -10,36 +11,21 @@
 class GrassmannianKernel(Kernel, ABC):
     """The parent class for Grassmannian kernels implemented in the :py:mod:`kernels` module ."""
 
-    def __init__(self):
-        super().__init__()
-
-    def calculate_kernel_matrix(self, points: list[GrassmannPoint], p: int = None):
+    def __init__(self, kernel_parameter: Union[int, float] = None):
         """
-        Compute the kernel matrix given a list of points on the Grassmann manifold.
-
-        :param points: Points on the Grassmann manifold
-        :param p: Number of independent p-planes of each Grassmann point.
-        :return: :class:`ndarray` The kernel matrix.
+        :param kernel_parameter: Number of independent p-planes of each Grassmann point.
         """
-        nargs = len(points)
-        # Define the pairs of points to compute the entries of the kernel matrix.
-        indices = range(nargs)
-        pairs = list(itertools.combinations_with_replacement(indices, 2))
-
-        # Estimate entries of the kernel matrix.
-        kernel = np.zeros((nargs, nargs))
-        for id_pair in range(np.shape(pairs)[0]):
-            i = pairs[id_pair][0]  # Point i
-            j = pairs[id_pair][1]  # Point j
-            if not p:
-                xi = points[i]
-                xj = points[j]
-            else:
-                xi = GrassmannPoint(points[i].data[:, :p])
-                xj = GrassmannPoint(points[j].data[:, :p])
-
-            # RiemannianDistance.check_rows(xi, xj)
-            kernel[i, j] = self.kernel_entry(xi, xj)
-            kernel[j, i] = kernel[i, j]
-
-        self.kernel_matrix = kernel
+        super().__init__(kernel_parameter)
+
+    def calculate_kernel_matrix(self, x: list[GrassmannPoint], s: list[GrassmannPoint]):
+        p = self.kernel_parameter
+        list1 = [point.data if not p else point.data[:, :p] for point in x]
+        list2 = [point.data if not p else point.data[:, :p] for point in s]
+        product = [self.element_wise_operation(point_pair)
+                   for point_pair in list(itertools.product(list1, list2))]
+        self.kernel_matrix = np.array(product).reshape(len(list1), len(list2))
+        return self.kernel_matrix
+
+    @abstractmethod
+    def element_wise_operation(self, xi_j: Tuple) -> float:
+        pass

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

@@ -12,6 +12,7 @@ class Kernel(ABC):
 
     def __init__(self, kernel_parameter: Union[int, float]):
         self.__kernel_parameter = kernel_parameter
+        self.kernel_matrix=None
 
     @property
     def kernel_parameter(self):

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

@@ -1,3 +1,2 @@
-from UQpy.utilities.kernels.baseclass.DRKernel import Kernel
 from UQpy.utilities.kernels.baseclass.EuclideanKernel import EuclideanKernel
 from UQpy.utilities.kernels.baseclass.GrassmannianKernel import GrassmannianKernel

src/UQpy/utilities/kernels/euclidean_kernels/Matern.py

@@ -1,9 +1,13 @@
-from UQpy.utilities.kernels.baseclass.Kernel import *
+from typing import Union
+
+import numpy as np
+
+from UQpy.utilities.kernels.baseclass.EuclideanKernel import EuclideanKernel
 from scipy.spatial.distance import cdist
 from scipy.special import gamma, kv
 
 
-class Matern(Kernel):
+class Matern(EuclideanKernel):
     def __init__(self, kernel_parameter: Union[int, float] = 1, nu=1.5):
         """
         Matern Kernel is a generalization of Radial Basis Function kernel.
@@ -18,16 +22,17 @@ def calculate_kernel_matrix(self, x, s):
         l = self.kernel_parameter
         stack = cdist(x / l, s / l, metric='euclidean')
         if self.nu == 0.5:
-            return np.exp(-np.abs(stack))
+            self.kernel_matrix = np.exp(-np.abs(stack))
         elif self.nu == 1.5:
-            return (1 + np.sqrt(3) * stack) * np.exp(-np.sqrt(3) * stack)
+            self.kernel_matrix = (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)
+            self.kernel_matrix = (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)
+            self.kernel_matrix = 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(self.nu, stack)
-            return tmp * tmp1 * tmp2
+            self.kernel_matrix = tmp * tmp1 * tmp2
+        return self.kernel_matrix

src/UQpy/utilities/kernels/euclidean_kernels/RBF.py

@@ -1,7 +1,12 @@
-from UQpy.utilities.kernels.baseclass.Kernel import *
+from typing import Union
 
+import numpy as np
 
-class RBF(Kernel):
+from UQpy.utilities.kernels.baseclass.EuclideanKernel import EuclideanKernel
+from UQpy.utilities.kernels.baseclass.Kernel import Kernel
+
+
+class RBF(EuclideanKernel):
     def __init__(self, kernel_parameter: Union[int, float] = 1.0):
         super().__init__(kernel_parameter)
 
@@ -13,4 +18,5 @@ def calculate_kernel_matrix(self, x, s):
         :params s: An array containing input points.
         """
         stack = Kernel.check_samples_and_return_stack(x / self.kernel_parameter, s / self.kernel_parameter)
-        return np.exp(np.sum(-0.5 * (stack ** 2), axis=2))
+        self.kernel_matrix = np.exp(np.sum(-0.5 * (stack ** 2), axis=2))
+        return self.kernel_matrix

src/UQpy/utilities/kernels/grassmannian_kernels/BinetCauchyKernel.py

@@ -0,0 +1,22 @@
+from typing import Tuple
+
+import numpy as np
+
+from UQpy.utilities.kernels import GrassmannianKernel
+
+
+class BinetCauchyKernel(GrassmannianKernel):
+    """
+    A class to calculate the Binet-Cauchy kernel.
+
+    """
+    def element_wise_operation(self, xi_j: Tuple) -> float:
+        """
+        Compute the Projection kernel entry for a tuple of points on the Grassmann manifold.
+
+        :param xi_j: Tuple of orthonormal matrices representing the grassmann points.
+        """
+        xi, xj = xi_j
+        r = np.dot(xi.T, xj)
+        det = np.linalg.det(r)
+        return det * det

src/UQpy/utilities/kernels/grassmannian_kernels/ProjectionKernel.py

@@ -0,0 +1,25 @@
+from typing import Union, Tuple
+
+import numpy as np
+
+from UQpy.utilities.kernels.baseclass.GrassmannianKernel import GrassmannianKernel
+
+
+class ProjectionKernel(GrassmannianKernel):
+
+    def __init__(self, kernel_parameter: Union[int, float] = None):
+        """
+        :param kernel_parameter: Number of independent p-planes of each Grassmann point.
+        """
+        super().__init__(kernel_parameter)
+
+    def element_wise_operation(self, xi_j: Tuple) -> float:
+        """
+        Compute the Projection kernel entry for a tuple of points on the Grassmann manifold.
+
+        :param xi_j: Tuple of orthonormal matrices representing the grassmann points.
+        """
+        xi, xj = xi_j
+        r = np.dot(xi.T, xj)
+        n = np.linalg.norm(r, "fro")
+        return n * n

src/UQpy/utilities/kernels/grassmannian_kernels/__init__.py

@@ -0,0 +1,2 @@
+from UQpy.utilities.kernels.grassmannian_kernels.ProjectionKernel import ProjectionKernel
+from UQpy.utilities.kernels.grassmannian_kernels.BinetCauchyKernel import BinetCauchyKernel