From 0d6bc35661d1c370fbe57a5d869ee15a07a4de5b Mon Sep 17 00:00:00 2001 From: Brian Delhaisse Date: Wed, 13 Nov 2019 20:11:30 +0100 Subject: [PATCH] correct GP/GPR to work with Python3 and add example --- examples/models/gpr.py | 98 +++++++++++++++++++ pyrobolearn/models/gp/gp.py | 189 ++++++++++++++++++++++++------------ 2 files changed, 227 insertions(+), 60 deletions(-) create mode 100644 examples/models/gpr.py diff --git a/examples/models/gpr.py b/examples/models/gpr.py new file mode 100644 index 0000000..98139c9 --- /dev/null +++ b/examples/models/gpr.py @@ -0,0 +1,98 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +"""Provide some examples using GPR. +""" + +import sys # to check if Python 3 or 2 +import numpy as np +import matplotlib.pyplot as plt + +from pyrobolearn.models.gp import GPR +from pyrobolearn.utils.converter import torch_to_numpy + + +# create data: Generate random sample following a sine curve +# Ref: https://scikit-learn.org/stable/auto_examples/mixture/plot_gmm_sin.html#sphx-glr-auto-examples-mixture-\ +# plot-gmm-sin-py +n_samples = 100 +np.random.seed(0) +X = np.zeros((n_samples, 2)) +step = 4. * np.pi / n_samples + +for i in range(X.shape[0]): + x = i * step - 6. + X[i, 0] = x + np.random.normal(0, 0.1) + X[i, 1] = 3. * (np.sin(x) + np.random.normal(0, .2)) + +x, y = X[:, 0], X[:, 1] +xlim, ylim = [-8, 8], [-8, 8] + +# plot data +plt.title('Training data') +plt.scatter(x, y) +plt.xlim(xlim) +plt.ylim(ylim) +plt.show() + +# create GPR +model = GPR() + +# plot prior possible functions +f = model.sample(x, num_samples=10, to_numpy=True) +plt.title('Sampled functions from prior distribution') +plt.scatter(x, y, alpha=0.3) +plt.plot(x, f.T) +plt.xlim(xlim) +plt.ylim(ylim) +plt.show() + +# compute log likelihoods +print("\nBefore training:") +print("Log likelihood: {}".format(model.log_likelihood(x, y, to_numpy=True))) +print("Log marginal likelihood: {}".format(model.log_marginal_likelihood(x, y, to_numpy=True))) + +# fit the data +optimizer = 'adam' # 'lbfgs' +model.fit(x, y, num_iters=100, optimizer=optimizer, verbose=True) + +# compute log likelihoods +print("\nAfter training:") +print("Log likelihood: {}".format(model.log_likelihood(x, y, to_numpy=True))) +print("Log marginal likelihood: {}".format(model.log_marginal_likelihood(x, y, to_numpy=True))) + +# sample function and plot it +f = model.sample(x, num_samples=1, to_numpy=True) +plt.title('one sampled function after training') +plt.scatter(x, y) +plt.plot(x, f.T, 'k', linewidth=2.) +plt.xlim(xlim) +plt.ylim(ylim) +plt.show() + +# predict prob +x_test = np.linspace(-10, 10, 101) +xlim, ylim = [-10, 10], [-10, 10] +mean_y, var_y = model.predict_prob(x_test, to_numpy=True) +std_y = np.sqrt(var_y) +plt.title('Prediction') +plt.scatter(x, y) +plt.plot(x_test, mean_y, 'b') +plt.fill_between(x_test, mean_y-2*std_y, mean_y+2*std_y, facecolor='green', alpha=0.3) +plt.xlim(xlim) +plt.ylim(ylim) +plt.show() + +# Another way to predict +pred = model.forward(x_test) +lower, upper = pred.confidence_region() + +plt.title("Another way to predict (see code)") +plt.plot(x, y, 'k*') +if sys.version_info[0] < 3: # Python 2 + plt.plot(x_test, torch_to_numpy(pred.mean()), 'b') +else: # Python 3 + plt.plot(x_test, torch_to_numpy(pred.mean), 'b') +plt.fill_between(x_test, torch_to_numpy(lower), torch_to_numpy(upper), alpha=0.5) +plt.xlim(xlim) +plt.ylim(ylim) +plt.show() diff --git a/pyrobolearn/models/gp/gp.py b/pyrobolearn/models/gp/gp.py index 0d25d01..c4b2cbd 100644 --- a/pyrobolearn/models/gp/gp.py +++ b/pyrobolearn/models/gp/gp.py @@ -21,6 +21,9 @@ import torch import gpytorch # import GPy +# to check Python version (if sys.version_info[0] < 3, then python 2) +import sys + # from pyrobolearn.models.model import Model __author__ = "Brian Delhaisse" @@ -32,6 +35,7 @@ __maintainer__ = "Brian Delhaisse" __email__ = "briandelhaisse@gmail.com" __status__ = "Development" +echo '# -*- coding: utf-8 -*-\n#!/usr/bin/env python' class GP(object): r"""Gaussian Process model @@ -51,11 +55,11 @@ class GP(object): - `GPFlow` (which uses TensorFlow) [5] References: - [1] "Gaussian Processes for Machine Learning", Rasmussen and Williams, 2006 - [2] GPyTorch: https://github.com/cornellius-gp/gpytorch - [3] GPyTorch examples: https://github.com/cornellius-gp/gpytorch/tree/master/examples - [4] GPy: https://gpy.readthedocs.io/en/deploy/ - [5] GPFlow: http://gpflow.readthedocs.io/en/latest/intro.html + - [1] "Gaussian Processes for Machine Learning", Rasmussen and Williams, 2006 + - [2] GPyTorch: https://github.com/cornellius-gp/gpytorch + - [3] GPyTorch examples: https://github.com/cornellius-gp/gpytorch/tree/master/examples + - [4] GPy: https://gpy.readthedocs.io/en/deploy/ + - [5] GPFlow: http://gpflow.readthedocs.io/en/latest/intro.html """ def fit(self, *args, **kwargs): @@ -66,10 +70,11 @@ class GPC(GP): r"""Gaussian Process Classification References: - [1] "Gaussian Processes for Machine Learning", Rasmussen and Williams, 2006 - [2] GPyTorch: https://github.com/cornellius-gp/gpytorch - [3] GPyTorch examples: https://github.com/cornellius-gp/gpytorch/tree/master/examples + - [1] "Gaussian Processes for Machine Learning", Rasmussen and Williams, 2006 + - [2] GPyTorch: https://github.com/cornellius-gp/gpytorch + - [3] GPyTorch examples: https://github.com/cornellius-gp/gpytorch/tree/master/examples """ + # TODO pass @@ -140,12 +145,20 @@ class ExactGPModel(gpytorch.models.ExactGP): # Methods # ########### - def forward(self, x): - r"""Return the prior probability density function :math:`p(f|x) = \mathcal{N}(\cdot | \mu(x), K(x,x))`.""" - mean_x = self.mean(x) - covar_x = self.kernel(x) - # return gpytorch.distributions.MultivariateNormal(mean_x, covar_x) - return gpytorch.random_variables.GaussianRandomVariable(mean_x, covar_x) + if sys.version_info[0] < 3: # Python 2.7 + def forward(self, x): + r"""Return the prior probability density function :math:`p(f|x) = \mathcal{N}(\cdot | \mu(x), K(x,x))`.""" + mean_x = self.mean(x) + covar_x = self.kernel(x) + # return gpytorch.distributions.MultivariateNormal(mean_x, covar_x) + return gpytorch.random_variables.GaussianRandomVariable(mean_x, covar_x) + else: # Python 3.* + def forward(self, x): + r"""Return the prior probability density function :math:`p(f|x) = \mathcal{N}(\cdot | \mu(x), K(x,x))`.""" + mean_x = self.mean(x) + covar_x = self.kernel(x) + # return gpytorch.random_variables.GaussianRandomVariable(mean_x, covar_x) + return gpytorch.distributions.MultivariateNormal(mean_x, covar_x) class GPR(GP): @@ -214,10 +227,10 @@ class GPR(GP): References: - [1] "Gaussian Processes for Machine Learning", Rasmussen and Williams, 2006 - [2] GPy: https://gpy.readthedocs.io/en/deploy/ - [3] GPyTorch: https://github.com/cornellius-gp/gpytorch - [4] GPFlow: http://gpflow.readthedocs.io/en/latest/intro.html + - [1] "Gaussian Processes for Machine Learning", Rasmussen and Williams, 2006 + - [2] GPy: https://gpy.readthedocs.io/en/deploy/ + - [3] GPyTorch: https://github.com/cornellius-gp/gpytorch + - [4] GPFlow: http://gpflow.readthedocs.io/en/latest/intro.html """ def __init__(self, mean=None, kernel=None, model=None, likelihood=None): @@ -231,7 +244,7 @@ class GPR(GP): model (None, gpytorch.module.Module): the prior GP model. If None, it will create `ExactGPModel()`, a GP model using the provided mean, kernel, and likelihood. likelihood (None, gpytorch.likelihoods.Likelihood): the likelihood pdf. If None, it will use the - `gpytorch.likelihoods.GaussianLikelihood()` + `gpytorch.likelihoods.GaussianLikelihood()`. """ # check model if model is None: @@ -428,26 +441,44 @@ class GPR(GP): likelihood = torch.exp(self.log_likelihood(x, y, to_numpy=False)) return self._convert(likelihood, to_numpy=to_numpy) - def log_likelihood(self, x, y, to_numpy=False): - r"""Evaluate the log likelihood: log p(y|f,x).""" - x = self._convert_to_torch(x) - y = self._convert_to_torch(y) - f = self.model(x) - log_likelihood = self.likelihood_prob.log_probability(f, y) - return self._convert(log_likelihood, to_numpy=to_numpy) + if sys.version_info[0] < 3: # Python 2.7 + def log_likelihood(self, x, y, to_numpy=False): + r"""Evaluate the log likelihood: log p(y|f,x).""" + x = self._convert_to_torch(x) + y = self._convert_to_torch(y) + f = self.model(x) + log_likelihood = self.likelihood_prob.log_probability(f, y) + return self._convert(log_likelihood, to_numpy=to_numpy) + else: + def log_likelihood(self, x, y, to_numpy=False): + r"""Evaluate the log likelihood: log p(y|f,x).""" + x = self._convert_to_torch(x) + y = self._convert_to_torch(y) + f = self.model(x) + log_likelihood = self.likelihood_prob.variational_log_probability(f, y) + return self._convert(log_likelihood, to_numpy=to_numpy) def marginal_likelihood(self, x, y, to_numpy=False): r"""Evaluate the marginal likelihood: p(y|x).""" ml = torch.exp(self.log_marginal_likelihood(x, y, to_numpy=False)) return self._convert(ml, to_numpy=to_numpy) - def log_marginal_likelihood(self, x, y, to_numpy=False): - r"""Evaluate the log marginal likelihood: log p(y|x).""" - x = self._convert_to_torch(x) - y = self._convert_to_torch(y) - f = self.model(x) - mll = self.mll(f, y) - return self._convert(mll[0], to_numpy=to_numpy) + if sys.version_info[0] < 3: # Python 2.7 + def log_marginal_likelihood(self, x, y, to_numpy=False): + r"""Evaluate the log marginal likelihood: log p(y|x).""" + x = self._convert_to_torch(x) + y = self._convert_to_torch(y) + f = self.model(x) + mll = self.mll(f, y) + return self._convert(mll[0], to_numpy=to_numpy) + else: + def log_marginal_likelihood(self, x, y, to_numpy=False): + r"""Evaluate the log marginal likelihood: log p(y|x).""" + x = self._convert_to_torch(x) + y = self._convert_to_torch(y) + f = self.model(x) + mll = self.mll(f, y) + return self._convert(mll, to_numpy=to_numpy) def fit(self, x, y, num_iters=100, tolerance=1e-5, optimizer=None, verbose=False): r"""Fit the input and output data; find optimal model hyperparameters. @@ -546,30 +577,56 @@ class GPR(GP): return self._convert_to_numpy(y.mean()) return y.mean() - def predict_prob(self, x, to_numpy=True): - r""" - Predict p(y|x) by returning the mean and the covariance arrays. + if sys.version_info[0] < 3: # Python 2.7 + def predict_prob(self, x, to_numpy=True): + r""" + Predict p(y|x) by returning the mean and the covariance arrays. - Args: - x (np.ndarray, torch.Tensor): input array - to_numpy (bool): if True, return a np.array + Args: + x (np.ndarray, torch.Tensor): input array + to_numpy (bool): if True, return a np.array - Returns: - np.ndarray, torch.Tensor: output mean array - np.ndarray, torch.Tensor: output covariance array - """ - x = self._convert_to_torch(x) + Returns: + np.ndarray, torch.Tensor: output mean array + np.ndarray, torch.Tensor: output covariance array + """ + x = self._convert_to_torch(x) - # compute p(f|x) - f = self.model(x) + # compute p(f|x) + f = self.model(x) - # compute p(y|f,x) - y = self.likelihood_prob(f) + # compute p(y|f,x) + y = self.likelihood_prob(f) - # return mean and covariance - if to_numpy: - return self._convert_to_numpy(y.mean()), self._convert_to_numpy(y.var()) # y.covar()) - return y.mean(), y.var() # y.covar() + # return mean and covariance + if to_numpy: + return self._convert_to_numpy(y.mean()), self._convert_to_numpy(y.var()) # y.covar()) + return y.mean(), y.var() # y.covar() + else: # Python 3.5 + def predict_prob(self, x, to_numpy=True): + r""" + Predict p(y|x) by returning the mean and the covariance arrays. + + Args: + x (np.ndarray, torch.Tensor): input array + to_numpy (bool): if True, return a np.array + + Returns: + np.ndarray, torch.Tensor: output mean array + np.ndarray, torch.Tensor: output covariance array + """ + x = self._convert_to_torch(x) + + # compute p(f|x) + f = self.model(x) + + # compute p(y|f,x) + y = self.likelihood_prob(f) + + # return mean and covariance + if to_numpy: + return self._convert_to_numpy(y.mean), self._convert_to_numpy(y.variance) # y.covariance) + return y.mean, y.variance # y.covariance def forward(self, x): r""" @@ -589,13 +646,25 @@ class GPR(GP): # return p(y|f,x) return self.likelihood_prob(f) - def sample(self, x, num_samples=1, to_numpy=True): - """Sample the function vector from the GP; i.e. f ~ p(f|x).""" - x = self._convert_to_torch(x) - f = self.model(x) - if to_numpy: - return self._convert_to_numpy(f.sample(num_samples)) - return f.sample(num_samples) + if sys.version_info[0] < 3: # Python 2.7 + def sample(self, x, num_samples=1, to_numpy=True): + """Sample the function vector from the GP; i.e. f ~ p(f|x).""" + x = self._convert_to_torch(x) + f = self.model(x) + if to_numpy: + return self._convert_to_numpy(f.sample(num_samples)) + return f.sample(num_samples) + + else: # Python 3.* + def sample(self, x, num_samples=1, to_numpy=True): + """Sample the function vector from the GP; i.e. f ~ p(f|x).""" + x = self._convert_to_torch(x) + f = self.model(x) + if isinstance(num_samples, int): + num_samples = (num_samples,) + if to_numpy: + return self._convert_to_numpy(f.sample(num_samples)) + return f.sample(num_samples) ############# # Operators # @@ -679,7 +748,7 @@ if __name__ == '__main__': plt.plot(x, f.T, 'k', linewidth=2.) # predict prob - x_test = torch.linspace(0, 1, 51).numpy() + x_test = torch.linspace(-0.5, 1.5, 51).numpy() mean_y, var_y = model.predict_prob(x_test, to_numpy=True) std_y = np.sqrt(var_y) plt.plot(x_test, mean_y, 'b')