diff --git a/docs/code/dimension_reduction/grassmann/plot_grassmann_kernel.py b/docs/code/dimension_reduction/grassmann/plot_grassmann_kernel.py index baea30bac..fa546b567 100644 --- a/docs/code/dimension_reduction/grassmann/plot_grassmann_kernel.py +++ b/docs/code/dimension_reduction/grassmann/plot_grassmann_kernel.py @@ -20,7 +20,7 @@ from UQpy.dimension_reduction.grassmann_manifold.projections.SVDProjection import SVDProjection from UQpy.dimension_reduction import GrassmannOperations from UQpy.utilities import GrassmannPoint -from UQpy.utilities.kernels import Kernel, ProjectionKernel +from UQpy.utilities.kernels import ProjectionKernel import sys # %% md @@ -73,10 +73,10 @@ # %% projection_kernel = ProjectionKernel() -projection_kernel.calculate_kernel_matrix(points=manifold_projection.u) +projection_kernel.calculate_kernel_matrix(x=manifold_projection.u, s=manifold_projection.u) kernel_psi = projection_kernel.kernel_matrix -projection_kernel.calculate_kernel_matrix(points=manifold_projection.v) +projection_kernel.calculate_kernel_matrix(x=manifold_projection.v, s=manifold_projection.v) kernel_phi = projection_kernel.kernel_matrix fig, (ax1, ax2) = plt.subplots(1, 2) @@ -92,9 +92,9 @@ # %% -projection_kernel.calculate_kernel_matrix(points=[manifold_projection.u[0], - manifold_projection.u[1], - manifold_projection.u[2]]) +projection_kernel.\ + calculate_kernel_matrix(x=[manifold_projection.u[0], manifold_projection.u[1], manifold_projection.u[2]], + s=[manifold_projection.u[0], manifold_projection.u[1], manifold_projection.u[2]]) kernel_01 = projection_kernel.kernel_matrix fig = plt.figure() @@ -111,17 +111,18 @@ class UserKernel(GrassmannianKernel): - def kernel_entry(self, xi: GrassmannPoint, xj: GrassmannPoint): - r = np.dot(xi.data.T, xj.data) + def element_wise_operation(self, xi_j) -> float: + xi, xj = xi_j + r = np.dot(xi.T, xj) det = np.linalg.det(r) return det * det user_kernel = UserKernel() -user_kernel.calculate_kernel_matrix(points=manifold_projection.u) +user_kernel.calculate_kernel_matrix(x=manifold_projection.u, s=manifold_projection.u) kernel_user_psi = user_kernel.kernel_matrix -user_kernel.calculate_kernel_matrix(points=manifold_projection.v) +user_kernel.calculate_kernel_matrix(x=manifold_projection.v, s=manifold_projection.v) kernel_user_phi = user_kernel.kernel_matrix fig, (ax1, ax2) = plt.subplots(1, 2) diff --git a/docs/code/reliability/form/plot_FORM_linear function_2d.py b/docs/code/reliability/form/FORM_linear function_2d.py similarity index 100% rename from docs/code/reliability/form/plot_FORM_linear function_2d.py rename to docs/code/reliability/form/FORM_linear function_2d.py diff --git a/docs/code/surrogates/gpr/plot_gpr_custom2D.py b/docs/code/surrogates/gpr/plot_gpr_custom2D.py index 764130894..b78166e35 100644 --- a/docs/code/surrogates/gpr/plot_gpr_custom2D.py +++ b/docs/code/surrogates/gpr/plot_gpr_custom2D.py @@ -32,13 +32,14 @@ import matplotlib.pyplot as plt from matplotlib import cm from matplotlib.ticker import LinearLocator, FormatStrFormatter -from UQpy.surrogates import GaussianProcessRegression, Matern +from UQpy.surrogates import GaussianProcessRegression # %% md # # Create a distribution object. # %% +from UQpy.utilities import Matern marginals = [Uniform(loc=0., scale=1.), Uniform(loc=0., scale=1.)] diff --git a/docs/code/surrogates/gpr/plot_gpr_no_noise.py b/docs/code/surrogates/gpr/plot_gpr_no_noise.py index 13c1240bc..01ccb3333 100644 --- a/docs/code/surrogates/gpr/plot_gpr_no_noise.py +++ b/docs/code/surrogates/gpr/plot_gpr_no_noise.py @@ -19,16 +19,17 @@ # %% -import numpy as np -import matplotlib.pyplot as plt import warnings +import matplotlib.pyplot as plt +import numpy as np + from UQpy.surrogates.gaussian_process.regression_models.LinearRegression import LinearRegression +from UQpy.utilities import RBF warnings.filterwarnings('ignore') from UQpy.utilities.MinimizeOptimizer import MinimizeOptimizer -from UQpy.utilities.FminCobyla import FminCobyla -from UQpy.surrogates import GaussianProcessRegression, NonNegative, RBF +from UQpy.surrogates import GaussianProcessRegression # %% md diff --git a/docs/code/surrogates/gpr/plot_gpr_noisy.py b/docs/code/surrogates/gpr/plot_gpr_noisy.py index 46d758274..a809da8ba 100644 --- a/docs/code/surrogates/gpr/plot_gpr_noisy.py +++ b/docs/code/surrogates/gpr/plot_gpr_noisy.py @@ -22,10 +22,13 @@ import numpy as np import matplotlib.pyplot as plt import warnings + +from UQpy.utilities import RBF + warnings.filterwarnings('ignore') from UQpy.utilities.MinimizeOptimizer import MinimizeOptimizer from UQpy.surrogates.gaussian_process.regression_models.LinearRegression import LinearRegression -from UQpy.surrogates import GaussianProcessRegression, RBF +from UQpy.surrogates import GaussianProcessRegression # %% md diff --git a/docs/code/surrogates/gpr/plot_gpr_sine.py b/docs/code/surrogates/gpr/plot_gpr_sine.py index 9056a7eaa..13c17c280 100644 --- a/docs/code/surrogates/gpr/plot_gpr_sine.py +++ b/docs/code/surrogates/gpr/plot_gpr_sine.py @@ -29,13 +29,14 @@ from UQpy.distributions import Gamma import numpy as np import matplotlib.pyplot as plt -from UQpy.surrogates import GaussianProcessRegression, RBF +from UQpy.surrogates import GaussianProcessRegression # %% md # # Create a distribution object. # %% +from UQpy.utilities import RBF marginals = [Gamma(a=2., loc=1., scale=3.)] diff --git a/docs/code/surrogates/pce/plot_pce_camel.py b/docs/code/surrogates/pce/plot_pce_camel.py index 7e0a0cafb..d15336958 100644 --- a/docs/code/surrogates/pce/plot_pce_camel.py +++ b/docs/code/surrogates/pce/plot_pce_camel.py @@ -9,14 +9,15 @@ .. math:: f(x) = \Big(4-2.1x_1^2 + \frac{x_1^4}{3} \Big)x_1^2 + x_1x_2 + (-4 + 4x_2^2)x_2^2 + **Description:** Dimensions: 2 **Input Domain:** This function is evaluated on the hypercube :math:`x_1 \in [-3, 3], x_2 \in [-2, 2]`. -**Global minimum:** :math:`f(x^*)=-1.0316,` at :math:`x^* = (0.0898, -0.7126)` and -:math:`(-0.0898, 0.7126)`. +**Global minimum:** :math:`f(x^*)=-1.0316,` at :math:`x^* = (0.0898, -0.7126)` and :math:`(-0.0898, 0.7126)`. **Reference:** Molga, M., & Smutnicki, C. Test functions for optimization needs (2005). Retrieved June 2013, from http://www.zsd.ict.pwr.wroc.pl/files/docs/functions.pdf. + """ # %% md diff --git a/docs/code/transformations/nataf/plot_nataf.py b/docs/code/transformations/nataf/nataf.py similarity index 100% rename from docs/code/transformations/nataf/plot_nataf.py rename to docs/code/transformations/nataf/nataf.py diff --git a/docs/requirements.txt b/docs/requirements.txt index fda2af357..0b93f7a8d 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -2,4 +2,5 @@ sphinx_autodoc_typehints sphinx_rtd_theme sphinx_gallery - sphinxcontrib_bibtex \ No newline at end of file + sphinxcontrib_bibtex + Sphinx==6.1.3 \ No newline at end of file diff --git a/docs/source/architecture.rst b/docs/source/architecture.rst index 3e4a926ed..03cb0d23b 100644 --- a/docs/source/architecture.rst +++ b/docs/source/architecture.rst @@ -129,7 +129,7 @@ Surrogates Another module that has extensively restructured in v4 is the surrogates. Apart from the :class:`.SROM` method which was retained as an independent algorithm, the previous Kriging functionality was removed. It is now replaced with :class:`.GaussianProcessRegression`. The functionality of the Gaussian is constructed using object composition, -and the specific implementation of :class:`.Regression` and :class:`.Kernel` abstract base classes. An additional +and the specific implementation of :class:`.Regression` and :class:`Kernel` abstract base classes. An additional functionality of constrained surrogates is added by implementing the :class:`.ConstraintsGPR` abstract class. The functionality of :class:`.PolynomialChaosExpansion` was rewritten from scratch to address some performance issues of v3. The Strategy Design pattern was used here as well, with three abstract base classes :class:`.Polynomials`, diff --git a/docs/source/bibliography.bib b/docs/source/bibliography.bib index c3f85d8b7..39d445a81 100644 --- a/docs/source/bibliography.bib +++ b/docs/source/bibliography.bib @@ -829,12 +829,19 @@ @article{saltelli_2002 keywords = {Sensitivity analysis, Sensitivity measures, Sensitivity indices, Importance measures}, } -@article{PTMCMC1, -title = {Parallel Tempering: Theory, Applications, and New Perspectives}, -author = {David J. Earl, Michael W. Deem}, -year = {2005}, -doi = {https://doi.org/10.48550/arXiv.physics/0508111} -} +@Article{PTMCMC1, +author ="Earl, David J. and Deem, Michael W.", +title ="Parallel tempering: Theory{,} applications{,} and new perspectives", +journal ="Phys. Chem. Chem. Phys.", +year ="2005", +volume ="7", +issue ="23", +pages ="3910-3916", +publisher ="The Royal Society of Chemistry", +doi ="10.1039/B509983H", +url ="http://dx.doi.org/10.1039/B509983H", +abstract ="We review the history of the parallel tempering simulation method. From its origins in data analysis{,} the parallel tempering method has become a standard workhorse of physicochemical simulations. We discuss the theory behind the method and its various generalizations. We mention a selected set of the many applications that have become possible with the introduction of parallel tempering{,} and we suggest several promising avenues for future research."} + @inproceedings{PTMCMC2, title = {Using Thermodynamic Integration to Calculate the Posterior Probability in Bayesian Model Selection Problems}, @@ -844,13 +851,22 @@ @inproceedings{PTMCMC2 author = {Paul M. Goggans and Ying Chi} } + @article{STMCMC_ChingChen, +author = {Jianye Ching and Yi-Chu Chen }, title = {Transitional Markov Chain Monte Carlo Method for Bayesian Model Updating, Model Class Selection, and Model Averaging}, -author = {Jianye Ching, Yi-Chu Chen}, +journal = {Journal of Engineering Mechanics}, +volume = {133}, +number = {7}, +pages = {816-832}, year = {2007}, -doi = {https://doi.org/10.1061/(ASCE)0733-9399(2007)133:7(816)} +doi = {10.1061/(ASCE)0733-9399(2007)133:7(816)}, +URL = {https://ascelibrary.org/doi/abs/10.1061/%28ASCE%290733-9399%282007%29133%3A7%28816%29}, +eprint = {https://ascelibrary.org/doi/pdf/10.1061/%28ASCE%290733-9399%282007%29133%3A7%28816%29} } + + @article{Kle2D, title = {Simulation of multi-dimensional random fields by Karhunen–Loève expansion}, journal = {Computer Methods in Applied Mechanics and Engineering}, diff --git a/docs/source/conf.py b/docs/source/conf.py index 4d44183fb..b4891bf87 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -25,7 +25,7 @@ ) # The full version, including alpha/beta/rc tags -release = "v4.0.0" +release = "v4.1.0" # -- General configuration --------------------------------------------------- @@ -173,7 +173,7 @@ # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. # -language = None +language = "en" pygments_style = None html_theme = "sphinx_rtd_theme" @@ -181,7 +181,7 @@ html_theme_options = { "logo_only": True, "style_nav_header_background": "#F0F0F0", - "vcs_pageview_mode": "view", + # "vcs_pageview_mode": "view", } github_url = "https://github.com/SURGroup/UQpy" diff --git a/docs/source/sensitivity/sobol.rst b/docs/source/sensitivity/sobol.rst index f2ecb537f..81b97e77b 100644 --- a/docs/source/sensitivity/sobol.rst +++ b/docs/source/sensitivity/sobol.rst @@ -81,8 +81,8 @@ Attributes .. autoattribute:: UQpy.sensitivity.SobolSensitivity.second_order_confidence_interval .. autoattribute:: UQpy.sensitivity.SobolSensitivity.total_order_confidence_interval .. autoattribute:: UQpy.sensitivity.SobolSensitivity.n_samples -.. autoattribute:: UQpy.sensitivity.Sobol.n_variables -.. autoattribute:: UQpy.sensitivity.Sobol.is_multi_output +.. autoattribute:: UQpy.sensitivity.SobolSensitivity.n_variables +.. autoattribute:: UQpy.sensitivity.SobolSensitivity.is_multi_output Examples diff --git a/docs/source/stochastic_process/karhunen_loeve_1d.rst b/docs/source/stochastic_process/karhunen_loeve_1d.rst index a952d7c84..11422269d 100644 --- a/docs/source/stochastic_process/karhunen_loeve_1d.rst +++ b/docs/source/stochastic_process/karhunen_loeve_1d.rst @@ -29,4 +29,4 @@ Examples .. toctree:: - Karhunen Loeve Examples <../auto_examples/stochastic_processes/karhunen_loeve/index> \ No newline at end of file + Karhunen Loeve Examples <../auto_examples/stochastic_processes/karhunen_loeve_1d/index> \ No newline at end of file diff --git a/docs/source/stochastic_process/karhunen_loeve_2d.rst b/docs/source/stochastic_process/karhunen_loeve_2d.rst index 5e23d00a6..3bb981a61 100644 --- a/docs/source/stochastic_process/karhunen_loeve_2d.rst +++ b/docs/source/stochastic_process/karhunen_loeve_2d.rst @@ -1,5 +1,5 @@ Karhunen Loève Expansion for Multi-Dimensional Fields ----------------------------- +----------------------------------------------------- The Karhunen Loève Expansion expands the stochastic field as follows: @@ -9,7 +9,7 @@ where :math:`\eta_{nk}(\theta)` are uncorrelated standardized normal random vari at :cite:`Kle2D` KarhunenLoeve2D Class -^^^^^^^^^^^^^^^^^^^^ +^^^^^^^^^^^^^^^^^^^^^ The :class:`.KarhunenLoeve2D` class is imported using the following command: @@ -23,7 +23,6 @@ Methods Attributes """""""""" .. autoattribute:: UQpy.stochastic_process.KarhunenLoeveExpansion2D.samples -.. autoattribute:: UQpy.stochastic_process.KarhunenLoeveExpansion2D.xi Examples """""""""" diff --git a/docs/source/surrogates/gpr.rst b/docs/source/surrogates/gpr.rst index 98af231bc..ba4267b47 100644 --- a/docs/source/surrogates/gpr.rst +++ b/docs/source/surrogates/gpr.rst @@ -150,16 +150,13 @@ User-Defined Kernel """""""""""""""""""""""""""" Adding a new kernel to the :class:`.GaussianProcessRegression` class is straightforward. This is done by creating a new class -that extends the :class:`UQpy.surrogates.gaussian_process.kernels.baseclass.Kernel` abstract base class. +that extends the :class:`.Kernel` abstract base class. This new class must have a method ``c(self, x, s, params)`` that takes as input the new points, training points and hyperparameters. Notice that the input ``params`` include lengthscales and process standard deviation, not noise standard deviation (even for noisy data). -The :class:`UQpy.surrogates.gaussian_process.kernels.baseclass.Kernel` class is imported using the following command: +The :class:`.Kernel` class is imported using the following command: ->>> from UQpy.surrogates.gaussian_process.kernels.baseclass.Kernel import Kernel - -.. autoclass:: UQpy.surrogates.gaussian_process.Kernel - :members: +>>> from UQpy.utilities.kernels.baseclass.Kernel import Kernel The method should return covariance matrix, i.e. a 2-D array with first dimension being the number of new points and second dimension being the number of training points. diff --git a/docs/source/utilities/kernels/euclidean_kernels.rst b/docs/source/utilities/kernels/euclidean_kernels.rst index 82d6737f7..26407eb02 100644 --- a/docs/source/utilities/kernels/euclidean_kernels.rst +++ b/docs/source/utilities/kernels/euclidean_kernels.rst @@ -29,7 +29,7 @@ Methods ~~~~~~~~~ .. autoclass:: UQpy.utilities.kernels.GaussianKernel - :members: kernel_entry, optimize_parameters + :members: optimize_parameters Attributes ~~~~~~~~~~ diff --git a/docs/source/utilities/kernels/index.rst b/docs/source/utilities/kernels/index.rst index 5cbc1fae6..5972171b6 100644 --- a/docs/source/utilities/kernels/index.rst +++ b/docs/source/utilities/kernels/index.rst @@ -20,7 +20,7 @@ The :class:`UQpy.utilities.kernels.baseclass.Kernel` class is imported using the >>> from UQpy.utilities.kernels.baseclass.Kernel import Kernel .. autoclass:: UQpy.utilities.kernels.baseclass.Kernel - :members: kernel_entry, optimize_parameters, calculate_kernel_matrix + :members: calculate_kernel_matrix Types of Kernels ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/src/UQpy/dimension_reduction/diffusion_maps/DiffusionMaps.py b/src/UQpy/dimension_reduction/diffusion_maps/DiffusionMaps.py index 5d47354e5..38762ed31 100644 --- a/src/UQpy/dimension_reduction/diffusion_maps/DiffusionMaps.py +++ b/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") diff --git a/src/UQpy/dimension_reduction/grassmann_manifold/projections/baseclass/GrassmannProjection.py b/src/UQpy/dimension_reduction/grassmann_manifold/projections/baseclass/GrassmannProjection.py index b32af3a63..acc3ac6cd 100644 --- a/src/UQpy/dimension_reduction/grassmann_manifold/projections/baseclass/GrassmannProjection.py +++ b/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): """ diff --git a/src/UQpy/sampling/mcmc/tempering_mcmc/baseclass/TemperingMCMC.py b/src/UQpy/sampling/mcmc/tempering_mcmc/baseclass/TemperingMCMC.py index 00d0921ba..6dac7172c 100644 --- a/src/UQpy/sampling/mcmc/tempering_mcmc/baseclass/TemperingMCMC.py +++ b/src/UQpy/sampling/mcmc/tempering_mcmc/baseclass/TemperingMCMC.py @@ -13,7 +13,7 @@ def __init__(self, pdf_intermediate=None, log_pdf_intermediate=None, args_pdf_in least two inputs :code:`x` (ndarray, point(s) at which to evaluate the function), and :code:`temper_param` (float, tempering parameter). Eit her `pdf_intermediate` or `log_pdf_intermediate` must be provided (`log_pdf_intermediate` is preferred). Within the code, the `log_pdf_intermediate` is evaluated as: - :code:`log_pdf_intermediate(x, temper_param, *args_pdf_intermediate)` + :code:`log_pdf_intermediate(x, temper_param, *args_pdf_intermediate)` where `args_pdf_intermediate` are additional positional arguments that are provided to the class via its `args_pdf_intermediate` input :param log_pdf_intermediate: see `pdf_intermediate` diff --git a/src/UQpy/surrogates/gaussian_process/GaussianProcessRegression.py b/src/UQpy/surrogates/gaussian_process/GaussianProcessRegression.py index 1d4741e58..41c114ae3 100755 --- a/src/UQpy/surrogates/gaussian_process/GaussianProcessRegression.py +++ b/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,7 +294,10 @@ 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 @@ -298,7 +305,9 @@ def predict(self, points, return_std: bool = False, hyperparameters: list = None 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] diff --git a/src/UQpy/surrogates/gaussian_process/__init__.py b/src/UQpy/surrogates/gaussian_process/__init__.py index 83b0f56c8..d9439caf9 100644 --- a/src/UQpy/surrogates/gaussian_process/__init__.py +++ b/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 * diff --git a/src/UQpy/surrogates/gaussian_process/kernels/Matern.py b/src/UQpy/surrogates/gaussian_process/kernels/Matern.py deleted file mode 100644 index 78b2edb97..000000000 --- a/src/UQpy/surrogates/gaussian_process/kernels/Matern.py +++ /dev/null @@ -1,32 +0,0 @@ -from UQpy.surrogates.gaussian_process.kernels.baseclass.Kernel import * -from scipy.spatial.distance import pdist, cdist, squareform -from scipy.special import gamma, kv - - -class Matern(Kernel): - def __init__(self, nu=1.5): - """ - Matern Kernel is a generalization of Radial Basis Function kernel. - - :params nu: Shape parameter. For nu=0.5, 1.5, 2.5 and infinity, matern coincides with the exponential, - matern-3/2, matern-5/2 and RBF covariance function, respectively. - """ - self.nu = nu - - 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 sigma**2 * np.exp(-np.abs(stack)) - elif self.nu == 1.5: - return sigma**2 * (1+np.sqrt(3)*stack)*np.exp(-np.sqrt(3)*stack) - elif self.nu == 2.5: - return sigma**2 * (1+np.sqrt(5)*stack+5*(stack**2)/3)*np.exp(-np.sqrt(5)*stack) - elif self.nu == np.inf: - 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(self.nu, stack) - return sigma**2 * tmp * tmp1 * tmp2 diff --git a/src/UQpy/surrogates/gaussian_process/kernels/RBF.py b/src/UQpy/surrogates/gaussian_process/kernels/RBF.py deleted file mode 100644 index bd56c2b91..000000000 --- a/src/UQpy/surrogates/gaussian_process/kernels/RBF.py +++ /dev/null @@ -1,15 +0,0 @@ -from UQpy.surrogates.gaussian_process.kernels.baseclass.Kernel import * - - -class RBF(Kernel): - def c(self, x, s, params): - """ - This method compute the RBF kernel on sample points 'x' and 's'. - - :params x: An array containing input samples. - :params s: An array containing input samples. - :params params: A list/array of hyperparameters containing length scales and the process variance. - """ - stack = Kernel.check_samples_and_return_stack(x/params[:-1], s/params[:-1]) - k = params[-1] ** 2 * np.exp(np.sum(-0.5 * (stack ** 2), axis=2)) - return k diff --git a/src/UQpy/surrogates/gaussian_process/kernels/__init__.py b/src/UQpy/surrogates/gaussian_process/kernels/__init__.py deleted file mode 100644 index c9c6bd495..000000000 --- a/src/UQpy/surrogates/gaussian_process/kernels/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -from UQpy.surrogates.gaussian_process.kernels.baseclass import * -from UQpy.surrogates.gaussian_process.kernels.RBF import RBF -from UQpy.surrogates.gaussian_process.kernels.Matern import Matern diff --git a/src/UQpy/surrogates/gaussian_process/kernels/baseclass/Kernel.py b/src/UQpy/surrogates/gaussian_process/kernels/baseclass/Kernel.py deleted file mode 100644 index 2336930bc..000000000 --- a/src/UQpy/surrogates/gaussian_process/kernels/baseclass/Kernel.py +++ /dev/null @@ -1,47 +0,0 @@ -from abc import ABC, abstractmethod -import numpy as np - - -class Kernel(ABC): - """ - Abstract base class of all Kernels. Serves as a template for creating new Gaussian Process covariance - functions. - """ - - @abstractmethod - def c(self, x, s, params): - """ - Abstract method that needs to be implemented by the user when creating a new Covariance function. - """ - pass - - @staticmethod - def check_samples_and_return_stack(x, s): - x_, s_ = np.atleast_2d(x), np.atleast_2d(s) - # Create stack matrix, where each block is x_i with all s - stack = np.tile( - np.swapaxes(np.atleast_3d(x_), 1, 2), (1, np.size(s_, 0), 1) - ) - 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_ diff --git a/src/UQpy/surrogates/gaussian_process/kernels/baseclass/__init__.py b/src/UQpy/surrogates/gaussian_process/kernels/baseclass/__init__.py deleted file mode 100644 index 305c426fb..000000000 --- a/src/UQpy/surrogates/gaussian_process/kernels/baseclass/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from UQpy.surrogates.gaussian_process.kernels.baseclass.Kernel import Kernel diff --git a/src/UQpy/utilities/kernels/BinetCauchyKernel.py b/src/UQpy/utilities/kernels/BinetCauchyKernel.py deleted file mode 100644 index 0d10fcb96..000000000 --- a/src/UQpy/utilities/kernels/BinetCauchyKernel.py +++ /dev/null @@ -1,25 +0,0 @@ -import numpy as np - -from UQpy.utilities.GrassmannPoint import GrassmannPoint -from UQpy.utilities.kernels import GrassmannianKernel - - -class BinetCauchyKernel(GrassmannianKernel): - """ - A class to calculate the Binet-Cauchy kernel. - - """ - def apply_method(self, points): - points.evaluate_matrix(self, self.calculate_kernel_matrix) - - def kernel_entry(self, xi: GrassmannPoint, xj: GrassmannPoint) -> float: - """ - Compute the Binet-Cauchy kernel entry for two points on the Grassmann manifold. - - :param xi: Orthonormal matrix representing the first point. - :param xj: Orthonormal matrix representing the second point. - - """ - r = np.dot(xi.data.T, xj.data) - det = np.linalg.det(r) - return det * det diff --git a/src/UQpy/utilities/kernels/GaussianKernel.py b/src/UQpy/utilities/kernels/GaussianKernel.py index 3eb0c03be..46326be26 100644 --- a/src/UQpy/utilities/kernels/GaussianKernel.py +++ b/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, diff --git a/src/UQpy/utilities/kernels/ProjectionKernel.py b/src/UQpy/utilities/kernels/ProjectionKernel.py deleted file mode 100644 index c23353baf..000000000 --- a/src/UQpy/utilities/kernels/ProjectionKernel.py +++ /dev/null @@ -1,25 +0,0 @@ -import numpy as np - -from UQpy.utilities.GrassmannPoint import GrassmannPoint -from UQpy.utilities.kernels import GrassmannianKernel - - -class ProjectionKernel(GrassmannianKernel): - """ - A class to calculate the Projection kernel - - """ - def apply_method(self, points): - points.evaluate_matrix(self, self.calculate_kernel_matrix) - - def kernel_entry(self, xi: GrassmannPoint, xj: GrassmannPoint) -> float: - """ - Compute the Projection kernel entry for two points on the Grassmann manifold. - - :param xi: Orthonormal matrix representing the first point. - :param xj: Orthonormal matrix representing the second point. - """ - r = np.dot(xi.data.T, xj.data) - n = np.linalg.norm(r, "fro") - kij = n * n - return kij diff --git a/src/UQpy/utilities/kernels/__init__.py b/src/UQpy/utilities/kernels/__init__.py index 360dfe64f..6e3877dab 100644 --- a/src/UQpy/utilities/kernels/__init__.py +++ b/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 diff --git a/src/UQpy/utilities/kernels/baseclass/EuclideanKernel.py b/src/UQpy/utilities/kernels/baseclass/EuclideanKernel.py index 1fb0e62a5..a71f7997a 100644 --- a/src/UQpy/utilities/kernels/baseclass/EuclideanKernel.py +++ b/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 diff --git a/src/UQpy/utilities/kernels/baseclass/GrassmannianKernel.py b/src/UQpy/utilities/kernels/baseclass/GrassmannianKernel.py index ec4249fe9..d2d219196 100644 --- a/src/UQpy/utilities/kernels/baseclass/GrassmannianKernel.py +++ b/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 diff --git a/src/UQpy/utilities/kernels/baseclass/Kernel.py b/src/UQpy/utilities/kernels/baseclass/Kernel.py index f24fd78ce..4b390b494 100644 --- a/src/UQpy/utilities/kernels/baseclass/Kernel.py +++ b/src/UQpy/utilities/kernels/baseclass/Kernel.py @@ -1,52 +1,40 @@ -import itertools from abc import ABC, abstractmethod from typing import Union import numpy as np -from UQpy.utilities import GrassmannPoint -from UQpy.utilities.ValidationTypes import NumpyFloatArray, Numpy2DFloatArray - class Kernel(ABC): """ - This is the baseclass for all kernels in :py:mod:`UQpy`. - - This serves as a blueprint to show the methods for kernels implemented in the :py:mod:`.kernels` module . + Abstract base class of all Kernels. Serves as a template for creating new Gaussian Process covariance + functions. """ - def __init__(self): - self.kernel_matrix: np.ndarray = None - """Kernel matrix defining the similarity between the points""" - def calculate_kernel_matrix(self, points: Union[Numpy2DFloatArray, list[GrassmannPoint]]): - """ - Using the kernel-specific :py:meth:`.kernel_entry` method, this function assembles the kernel matrix. + def __init__(self, kernel_parameter: Union[int, float]): + self.__kernel_parameter = kernel_parameter + self.kernel_matrix=None - :param points: Set of data points. Depending on the type of kernel these should be either :class:`numpy.ndarray` - or of type :class:`.GrassmannPoint`. - """ - pass + @property + def kernel_parameter(self): + return self.__kernel_parameter - def optimize_parameters(self, data: Union[Numpy2DFloatArray, list[GrassmannPoint]], **kwargs_optimization): - """ - This serves as a blueprint function in case a kernel provides the ability to optimize its parameters. In that - case, the kernel will override this method, and store the optimized parameters in the kernel's attributes. + @kernel_parameter.setter + def kernel_parameter(self, value): + self.__kernel_parameter = value - :param data: Set of data points. - :param kwargs_optimization: Keyword arguments containing any extra parameters needed to perform the - optimization. These will be unique to the specific kernel. - """ - pass @abstractmethod - def kernel_entry(self, xi: Union[Numpy2DFloatArray, GrassmannPoint], - xj: Union[Numpy2DFloatArray, GrassmannPoint]): + def calculate_kernel_matrix(self, x, s): """ - Given two points, this method calculates the respective kernel entry. Each concrete kernel implementation must - override this method and provide its own implementation. - - :param xi: First point. - :param xj: Second point. - :return: Float representing the kernel entry. + Abstract method that needs to be implemented by the user when creating a new Covariance function. """ pass + + @staticmethod + def check_samples_and_return_stack(x, s): + x_, s_ = np.atleast_2d(x), np.atleast_2d(s) + # Create stack matrix, where each block is x_i with all s + stack = np.tile( + np.swapaxes(np.atleast_3d(x_), 1, 2), (1, np.size(s_, 0), 1) + ) - np.tile(s_, (np.size(x_, 0), 1, 1)) + return stack diff --git a/src/UQpy/utilities/kernels/baseclass/__init__.py b/src/UQpy/utilities/kernels/baseclass/__init__.py index b069d28fe..20f0ff402 100644 --- a/src/UQpy/utilities/kernels/baseclass/__init__.py +++ b/src/UQpy/utilities/kernels/baseclass/__init__.py @@ -1,3 +1,3 @@ -from UQpy.utilities.kernels.baseclass.Kernel import Kernel from UQpy.utilities.kernels.baseclass.EuclideanKernel import EuclideanKernel from UQpy.utilities.kernels.baseclass.GrassmannianKernel import GrassmannianKernel +from UQpy.utilities.kernels.baseclass.Kernel import Kernel diff --git a/src/UQpy/utilities/kernels/euclidean_kernels/Matern.py b/src/UQpy/utilities/kernels/euclidean_kernels/Matern.py new file mode 100644 index 000000000..e21a72d1a --- /dev/null +++ b/src/UQpy/utilities/kernels/euclidean_kernels/Matern.py @@ -0,0 +1,38 @@ +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(EuclideanKernel): + def __init__(self, kernel_parameter: Union[int, float] = 1, nu=1.5): + """ + Matern Kernel is a generalization of Radial Basis Function kernel. + + :params nu: Shape parameter. For nu=0.5, 1.5, 2.5 and infinity, matern coincides with the exponential, + matern-3/2, matern-5/2 and RBF covariance function, respectively. + """ + super().__init__(kernel_parameter) + self.nu = nu + + def calculate_kernel_matrix(self, x, s): + l = self.kernel_parameter + stack = cdist(x / l, s / l, metric='euclidean') + if self.nu == 0.5: + self.kernel_matrix = np.exp(-np.abs(stack)) + elif self.nu == 1.5: + self.kernel_matrix = (1 + np.sqrt(3) * stack) * np.exp(-np.sqrt(3) * stack) + elif self.nu == 2.5: + self.kernel_matrix = (1 + np.sqrt(5) * stack + 5 * (stack ** 2) / 3) * np.exp(-np.sqrt(5) * stack) + elif self.nu == np.inf: + 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) + self.kernel_matrix = tmp * tmp1 * tmp2 + return self.kernel_matrix diff --git a/src/UQpy/utilities/kernels/euclidean_kernels/RBF.py b/src/UQpy/utilities/kernels/euclidean_kernels/RBF.py new file mode 100644 index 000000000..768e4372b --- /dev/null +++ b/src/UQpy/utilities/kernels/euclidean_kernels/RBF.py @@ -0,0 +1,22 @@ +from typing import Union + +import numpy as np + +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) + + def calculate_kernel_matrix(self, x, s): + """ + This method compute the RBF kernel on sample points 'x' and 's'. + + :params x: An array containing training points. + :params s: An array containing input points. + """ + stack = Kernel.check_samples_and_return_stack(x / self.kernel_parameter, s / self.kernel_parameter) + self.kernel_matrix = np.exp(np.sum(-0.5 * (stack ** 2), axis=2)) + return self.kernel_matrix diff --git a/src/UQpy/utilities/kernels/euclidean_kernels/__init__.py b/src/UQpy/utilities/kernels/euclidean_kernels/__init__.py new file mode 100644 index 000000000..e53d45504 --- /dev/null +++ b/src/UQpy/utilities/kernels/euclidean_kernels/__init__.py @@ -0,0 +1,2 @@ +from UQpy.utilities.kernels.euclidean_kernels.Matern import Matern +from UQpy.utilities.kernels.euclidean_kernels.RBF import RBF \ No newline at end of file diff --git a/src/UQpy/utilities/kernels/grassmannian_kernels/BinetCauchyKernel.py b/src/UQpy/utilities/kernels/grassmannian_kernels/BinetCauchyKernel.py new file mode 100644 index 000000000..91b9eeec7 --- /dev/null +++ b/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 diff --git a/src/UQpy/utilities/kernels/grassmannian_kernels/ProjectionKernel.py b/src/UQpy/utilities/kernels/grassmannian_kernels/ProjectionKernel.py new file mode 100644 index 000000000..a98fa46b4 --- /dev/null +++ b/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 diff --git a/src/UQpy/utilities/kernels/grassmannian_kernels/__init__.py b/src/UQpy/utilities/kernels/grassmannian_kernels/__init__.py new file mode 100644 index 000000000..43b4f509f --- /dev/null +++ b/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 diff --git a/tests/unit_tests/dimension_reduction/test_dmaps.py b/tests/unit_tests/dimension_reduction/test_dmaps.py index edc2e6bb9..a677d4b40 100644 --- a/tests/unit_tests/dimension_reduction/test_dmaps.py +++ b/tests/unit_tests/dimension_reduction/test_dmaps.py @@ -21,7 +21,7 @@ def test_dmaps_swiss_roll(): Y0 = 1. / 6 * (phi + sigma * xi) * np.cos(phi) swiss_roll = np.array([X0, Y0, Z0]).transpose() - dmaps = DiffusionMaps(data=swiss_roll, kernel=GaussianKernel(epsilon=0.5), + dmaps = DiffusionMaps(data=swiss_roll, kernel=GaussianKernel(kernel_parameter=0.5), alpha=0.5, n_eigenvectors=3, is_sparse=True, n_neighbors=100) evals = dmaps.eigenvalues @@ -54,7 +54,7 @@ def test_dmaps_circular(): X = np.array([x, y, z]).transpose() dmaps = DiffusionMaps(data=X, alpha=1, n_eigenvectors=3, - kernel=GaussianKernel(epsilon=0.3)) + kernel=GaussianKernel(kernel_parameter=0.3)) evals = dmaps.eigenvalues assert np.round(evals[0], 5) == 1.0 diff --git a/tests/unit_tests/dimension_reduction/test_kernel.py b/tests/unit_tests/dimension_reduction/test_kernel.py index c48f012a2..53434368a 100644 --- a/tests/unit_tests/dimension_reduction/test_kernel.py +++ b/tests/unit_tests/dimension_reduction/test_kernel.py @@ -1,9 +1,7 @@ -import sys - +from UQpy.utilities import ProjectionKernel from UQpy.utilities.GrassmannPoint import GrassmannPoint from UQpy.dimension_reduction.grassmann_manifold.projections.SVDProjection import SVDProjection -from UQpy.utilities.kernels.ProjectionKernel import ProjectionKernel -from UQpy.utilities.kernels.BinetCauchyKernel import BinetCauchyKernel +from UQpy.utilities.kernels.grassmannian_kernels.BinetCauchyKernel import BinetCauchyKernel from UQpy.utilities.kernels.GaussianKernel import GaussianKernel import numpy as np @@ -16,12 +14,11 @@ def test_kernel_projection(): [-0.63305978, 0.51850616]])) points = [xi, xj, xk] k = ProjectionKernel() - k.calculate_kernel_matrix(points) + k.calculate_kernel_matrix(points, points) kernel = np.matrix.round(k.kernel_matrix, 4) assert np.allclose(kernel, np.array([[2, 1.0063, 1.2345], [1.0063, 2, 1.0101], [1.2345, 1.0101, 2]])) - def test_kernel_binet_cauchy(): xi = GrassmannPoint(np.array([[-np.sqrt(2) / 2, -np.sqrt(2) / 4], [np.sqrt(2) / 2, -np.sqrt(2) / 4], [0, -np.sqrt(3) / 2]])) @@ -30,7 +27,7 @@ def test_kernel_binet_cauchy(): points = [xi, xj, xk] kernel = BinetCauchyKernel() - kernel.calculate_kernel_matrix(points) + kernel.calculate_kernel_matrix(points, points) kernel = np.matrix.round(kernel.kernel_matrix, 4) assert np.allclose(kernel, np.array([[1, 0.0063, 0.2345], [0.0063, 1, 0.0101], [0.2345, 0.0101, 1]])) @@ -41,13 +38,13 @@ def test_kernel_gaussian_1d(): xj = np.array([0.2, -1, 3, 5]) xk = np.array([1, 2, 3, 4]) points = [xi, xj, xk] - gaussian = GaussianKernel(epsilon=2.0) - gaussian.calculate_kernel_matrix(points) + gaussian = GaussianKernel(kernel_parameter=2.0) + gaussian.calculate_kernel_matrix(points, points) assert np.allclose(np.matrix.round(gaussian.kernel_matrix, 4), np.array([[1., 0.26447726, 1.], [0.26447726, 1., 0.26447726], [1, 0.26447726, 1]]), atol=1e-04) - assert np.round(gaussian.epsilon, 4) == 2 + assert np.round(gaussian.kernel_parameter, 4) == 2 def test_kernel_gaussian_2d(): @@ -56,12 +53,12 @@ def test_kernel_gaussian_2d(): xk = np.array([[-0.69535592, -0.0546034], [-0.34016974, -0.85332868], [-0.63305978, 0.51850616]]) points = [xi, xj, xk] gaussian = GaussianKernel() - gaussian.calculate_kernel_matrix(points) + gaussian.calculate_kernel_matrix(points, points) assert np.allclose(np.matrix.round(gaussian.kernel_matrix, 4), np.array([[1., 0.39434829, 0.15306655], [0.39434829, 1., 0.06422136], [0.15306655, 0.06422136, 1.]]), atol=1e-4) - assert np.round(gaussian.epsilon, 4) == 1.0 + assert np.round(gaussian.kernel_parameter, 4) == 1.0 sol0 = np.array([[0.61415, 1.03029, 1.02001, 0.57327, 0.79874, 0.73274], @@ -104,6 +101,6 @@ def test_kernel(): manifold_projection = SVDProjection(Solutions, p="max") kernel = ProjectionKernel() - kernel.calculate_kernel_matrix(manifold_projection.u) + kernel.calculate_kernel_matrix(manifold_projection.u, manifold_projection.u) assert np.round(kernel.kernel_matrix[0, 1], 8) == 6.0 diff --git a/tests/unit_tests/sampling/test_adaptive_kriging.py b/tests/unit_tests/sampling/test_adaptive_kriging.py index 4fcf3b1e6..e191feb32 100644 --- a/tests/unit_tests/sampling/test_adaptive_kriging.py +++ b/tests/unit_tests/sampling/test_adaptive_kriging.py @@ -1,6 +1,7 @@ import pytest -from UQpy import GaussianProcessRegression, RBF, LinearRegression +from UQpy import GaussianProcessRegression, LinearRegression +from UQpy.utilities.kernels.euclidean_kernels.RBF import RBF from UQpy.run_model.model_execution.PythonModel import PythonModel from UQpy.utilities.MinimizeOptimizer import MinimizeOptimizer diff --git a/tests/unit_tests/sampling/test_refined_stratified.py b/tests/unit_tests/sampling/test_refined_stratified.py index 13afea15f..03e7921d4 100644 --- a/tests/unit_tests/sampling/test_refined_stratified.py +++ b/tests/unit_tests/sampling/test_refined_stratified.py @@ -1,7 +1,8 @@ import pytest from beartype.roar import BeartypeCallHintPepParamException -from UQpy import GaussianProcessRegression, RBF, LinearRegression +from UQpy import GaussianProcessRegression, LinearRegression +from UQpy.utilities.kernels.euclidean_kernels.RBF import RBF from UQpy.run_model.model_execution.PythonModel import PythonModel from UQpy.utilities.MinimizeOptimizer import MinimizeOptimizer from UQpy.sampling.stratified_sampling.refinement.GradientEnhancedRefinement import GradientEnhancedRefinement diff --git a/tests/unit_tests/surrogates/test_gpr.py b/tests/unit_tests/surrogates/test_gpr.py index 504369b97..9a8a1be8f 100644 --- a/tests/unit_tests/surrogates/test_gpr.py +++ b/tests/unit_tests/surrogates/test_gpr.py @@ -1,18 +1,14 @@ import pytest from beartype.roar import BeartypeCallHintPepParamException +from UQpy.utilities.kernels.euclidean_kernels import RBF, Matern from UQpy.utilities.MinimizeOptimizer import MinimizeOptimizer from UQpy.utilities.FminCobyla import FminCobyla from UQpy.surrogates.gaussian_process.GaussianProcessRegression import GaussianProcessRegression -# # from UQpy.utilities.strata.Rectangular import Rectangular -# # from UQpy.sampling.StratifiedSampling import StratifiedSampling -# # from UQpy.RunModel import RunModel -# # from UQpy.distributions.collection.Uniform import Uniform import numpy as np import shutil from UQpy.surrogates.gaussian_process.constraints import NonNegative from UQpy.surrogates.gaussian_process.regression_models import LinearRegression, ConstantRegression, QuadraticRegression -from UQpy.surrogates.gaussian_process.kernels import RBF, Matern samples = np.linspace(0, 5, 20).reshape(-1, 1) @@ -49,35 +45,34 @@ gpr3 = GaussianProcessRegression(regression_model=constant_reg, kernel=RBF(), hyperparameters=[2.852, 2.959], random_state=1, optimizer=optimizer3, optimize_constraints=non_negative, bounds=[[0.1, 5], [0.1, 3]]) -gpr3.fit(samples=samples, values=values+1.05) +gpr3.fit(samples=samples, values=values + 1.05) gpr4 = GaussianProcessRegression(kernel=RBF(), hyperparameters=[2.852, 2.959, 0.001], random_state=1, optimizer=optimizer3, bounds=[[0.1, 5], [0.1, 3], [1e-6, 1e-1]], noise=True) gpr4.fit(samples=samples, values=values) - def test_predict1(): - prediction = np.round(gpr.predict([[1], [np.pi/2], [np.pi]], True), 3) - expected_prediction = np.array([[0.54, 0., -1.], [0., 0., 0.]]) + 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_predict2(): - 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]]) + 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_predict3(): - prediction = np.round(gpr3.predict([[1], [2*np.pi], [np.pi]], True), 3) + prediction = np.round(gpr3.predict([[1], [2 * np.pi], [np.pi]], True), 3) expected_prediction = np.array([[1.59, 2.047, 0.05], [0., 0.006, 0.]]) assert (expected_prediction == prediction).all() def test_predict4(): prediction = np.round(gpr4.predict([[1], [2 * np.pi], [np.pi]], True), 3) - expected_prediction = np.array([[0.54, 0.983, -1.], [0., 0.04, 0.]]) + expected_prediction = np.array([[0.54, 0.983, -1.], [0., 0.04, 0.]]) assert np.isclose(prediction, expected_prediction, atol=0.05).all() @@ -85,10 +80,12 @@ def test_rbf(): """ Test RBF kernel for 1-d input model """ - cov = rbf_kernel.c(points11, points12, hyperparameters1) + rbf_kernel.kernel_parameter = hyperparameters1[:-1] + sigma = hyperparameters1[-1] + cov = sigma ** 2 * rbf_kernel.calculate_kernel_matrix(points11, points12) 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]]) + 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() @@ -96,7 +93,9 @@ def test_rbf2(): """ Test RBF kernel for 2-d input model """ - cov = rbf_kernel.c(points21, points22, hyperparameters2) + rbf_kernel.kernel_parameter = hyperparameters2[:-1] + sigma = hyperparameters2[-1] + cov = sigma ** 2 * rbf_kernel.calculate_kernel_matrix(points21, points22) 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]]) @@ -107,18 +106,22 @@ 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) + matern_kernel_inf.kernel_parameter = hyperparameters1[:-1] + sigma = hyperparameters1[-1] + cov = sigma ** 2 * matern_kernel_inf.calculate_kernel_matrix(points11, points12) 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]]) + 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(): +def test_matern_inf2(): """ Test Matern kernel (nu=inf) for 2-d input model """ - cov = matern_kernel_inf.c(points21, points22, hyperparameters2) + matern_kernel_inf.kernel_parameter = hyperparameters2[:-1] + sigma = hyperparameters2[-1] + cov = sigma ** 2 * matern_kernel_inf.calculate_kernel_matrix(points21, points22) 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]]) @@ -129,10 +132,12 @@ 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) + matern_kernel_1_2.kernel_parameter = hyperparameters1[:-1] + sigma = hyperparameters1[-1] + cov = sigma ** 2 * matern_kernel_1_2.calculate_kernel_matrix(points11, points12) 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]]) + 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() @@ -140,7 +145,9 @@ 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) + matern_kernel_1_2.kernel_parameter = hyperparameters2[:-1] + sigma = hyperparameters2[-1] + cov = sigma ** 2 * matern_kernel_1_2.calculate_kernel_matrix(points21, points22) 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() @@ -150,10 +157,12 @@ 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) + matern_kernel_3_2.kernel_parameter = hyperparameters1[:-1] + sigma = hyperparameters1[-1] + cov = sigma ** 2 * matern_kernel_3_2.calculate_kernel_matrix(points11, points12) 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]]) + 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() @@ -161,10 +170,12 @@ 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) + matern_kernel_5_2.kernel_parameter = hyperparameters1[:-1] + sigma = hyperparameters1[-1] + cov = sigma ** 2 * matern_kernel_5_2.calculate_kernel_matrix(points11, points12) 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]]) + 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() @@ -172,9 +183,11 @@ def test_matern_2_1(): """ Test Matern kernel (nu=2) for 1-d input model """ - cov = matern_kernel_2_1.c(points11, points12, hyperparameters1) + matern_kernel_2_1.kernel_parameter = hyperparameters1[:-1] + sigma = hyperparameters1[-1] + cov = sigma ** 2 * matern_kernel_2_1.calculate_kernel_matrix(points11, points12) covariance_matrix = np.round(cov, 3) - expected_covariance_matrix = np.array([[70.894, np.nan, 70.894], [34.25, 70.894, 34.25], + 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) @@ -426,4 +439,4 @@ def test_hyperparameters_length(): # # Put a legend to the right of the current axis # ax.legend(loc='center left', bbox_to_anchor=(1, 0.5)) # plt.grid() -# plt.show() \ No newline at end of file +# plt.show()