apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-06-11 15:15:15 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

more massive and destructive changes everywhere

Browse source
This commit is contained in:
James Hensman 2013-12-05 15:09:31 -05:00
parent 881800126f
commit 5ec64d2279
18 changed files with 202 additions and 166 deletions
Download patch file Download diff file Expand all files Collapse all files

28
GPy/core/gp.py
View file

@ -3,6 +3,9 @@
import numpy as np
from gp_base import GPBase
from ..util.linalg import dtrtrs
from ..inference.latent_function_inference import exact_gaussian_inference, expectation_propagation
from .. import likelihoods
class GP(GPBase):
"""
@ -18,14 +21,25 @@ class GP(GPBase):
.. Note:: Multiple independent outputs are allowed using columns of Y
"""
def __init__(self, X, likelihood, kernel, normalize_X=False):
super(GP, self).__init__(X, likelihood, kernel, normalize_X=normalize_X)
self.Posterior = self.inference_method.inference(K, likelihood, self.Y)
def __init__(self, X, Y, kernel, likelihood, inference_method=None, name='gp'):
if inference_method is None:
if isinstance(likelihood, likelihoods.Gaussian):
inference_method = exact_gaussian_inference.ExactGaussianInference()
else:
inference_method = expectation_propagation
print "defaulting to ", inference_method, "for latent function inference"
super(GP, self).__init__(X, Y, kernel, likelihood, inference_method, name)
self.parameters_changed()
def parameters_changed(self):
super(GP, self).parameters_changed()
self.K = self.kern.K(self.X)
self.Posterior = self.inference_method.inference(K, likelihood, self.Y)
self.posterior = self.inference_method.inference(self.K, self.likelihood, self.Y)
def dL_dtheta_K(self):
return self.kern.dK_dtheta(self.posterior.dL_dK, self.X)
def log_likelihood(self):
return self.posterior.log_marginal
@ -41,15 +55,15 @@ class GP(GPBase):
"""
Kx = self.kern.K(_Xnew, self.X, which_parts=which_parts).T
LiKx, _ = dptrrs(self.posterior._woodbury_chol, np.asfortranarray(Kx), lower=1)
LiKx, _ = dtrtrs(self.posterior._woodbury_chol, np.asfortranarray(Kx), lower=1)
mu = np.dot(Kx.T, self.posterior._woodbury_vector)
if full_cov:
Kxx = self.kern.K(_Xnew, which_parts=which_parts)
var = Kxx - tdot(LiKx.T)
else:
Kxx = self.kern.Kdiag(_Xnew, which_parts=which_parts)
var = Kxx - np.sum(LiKx.T*LiKx, 0)
var = var.reshape(self.num_data, 1)
var = Kxx - np.sum(LiKx*LiKx, 0)
var = var.reshape(-1, 1)
return mu, var
def predict(self, Xnew, which_parts='all', full_cov=False, **likelihood_args):

24
GPy/core/gp_base.py
View file

@ -4,34 +4,36 @@ import warnings
from .. import kern
from ..util.plot import gpplot, Tango, x_frame1D, x_frame2D
from model import Model
from ..core.parameter import ObservableArray
from parameter import ObservableArray
from .. import likelihoods
class GPBase(Model):
"""
Gaussian process base model for holding shared behaviour between
sparse_GP and GP models.
"""
def __init__(self, X, Y, kernel, normalize_X=False, inference_method=None, name=''):
def __init__(self, X, Y, kernel, likelihood, inference_method, name=''):
super(GPBase, self).__init__(name)
assert X.ndim == 2
self.X = ObservableArray(X)
self.num_data, self.input_dim = self.X.shape
assert Y.ndim == 2
self.Y = ObservableArray(Y)
assert Y.shape[0] == self.num_data
_, self.output_dim = self.Y.shape
assert isinstance(kernel, kern.kern)
self.kern = kernel
assert isinstance(likelihood, likelihoods.Likelihood)
self.likelihood = likelihood
if inference_method is None:
self.inference_method = self.likelihood.preferred_inference_method
print "defaulting to ", inference_method, "for latent function inference"
else:
self.inference_method = inference_method
self.inference_method = inference_method
assert self.X.shape[0] == Y.shape[0]
self.num_data, self.output_dim = self.Y.shape
# reinstate this later TODO
normalize_X = False
if normalize_X:
self._Xoffset = X.mean(0)[None, :]
@ -41,8 +43,8 @@ class GPBase(Model):
self._Xoffset = np.zeros((1, self.input_dim))
self._Xscale = np.ones((1, self.input_dim))
self.add_parameter(self.kern, gradient=lambda:self.kern.dK_dtheta(self.dL_dK, self.X))
self.add_parameter(self.likelihood, gradient=lambda:self.likelihood._gradients(partial=np.diag(self.dL_dK)))
self.add_parameter(self.kern, gradient=self.dL_dtheta_K)
self.add_parameter(self.likelihood, gradient=lambda:self.posterior.dL_dtheta_lik)
def posterior_samples_f(self,X,size=10,which_parts='all',full_cov=True):
"""

1
GPy/core/sparse_gp.py
View file

@ -5,7 +5,6 @@ import numpy as np
import pylab as pb
from ..util.linalg import mdot, jitchol, tdot, symmetrify, backsub_both_sides, chol_inv, dtrtrs, dpotrs, dpotri
from scipy import linalg
from ..likelihoods import Gaussian, EP,EP_Mixed_Noise
from gp_base import GPBase
from GPy.core.parameter import Param

1
GPy/core/svigp.py
View file

@ -5,7 +5,6 @@ import numpy as np
import pylab as pb
from .. import kern
from ..util.linalg import pdinv, mdot, tdot, dpotrs, dtrtrs, jitchol, backsub_both_sides
from ..likelihoods import EP
from gp_base import GPBase
from model import Model
import time

3
GPy/inference/latent_function_inference/__init__.py
View file

@ -22,3 +22,6 @@ use this posterior object for making predictions, optimizing hyper-parameters,
etc.
"""
from exact_gaussian_inference import ExactGaussianInference
expectation_propagation = 'foo' # TODO

17
GPy/inference/latent_function_inference/exact_gaussian_inference.py
View file

@ -2,11 +2,12 @@
# Licensed under the BSD 3-clause license (see LICENSE.txt)
from posterior import Posterior
from .../util.linalg import pdinv, dpotrs, tdot
from ...util.linalg import pdinv, dpotrs, tdot
import numpy as np
log_2_pi = np.log(2*np.pi)
class exact_gaussian_inference(object):
class ExactGaussianInference(object):
"""
An object for inference when the likelihood is Gaussian.
@ -17,9 +18,9 @@ class exact_gaussian_inference(object):
"""
def __init__(self):
self._YYTfactor_cache = caching.cache()
pass#self._YYTfactor_cache = caching.cache()
def get_YYTfactor(self, Y):
def get_YYTfactor(self, Y):
"""
find a matrix L which satisfies LLT = YYT.
@ -38,16 +39,16 @@ class exact_gaussian_inference(object):
"""
YYT_factor = self.get_YYTfactor(Y)
Wi, LW, LWi, W_logdet = pdinv(K + likelhood.covariance(Y, Y_metadata))
Wi, LW, LWi, W_logdet = pdinv(K + likelihood.covariance_matrix(Y, Y_metadata))
alpha, _ = dpotrs(LW, YYT_factor, lower=1)
dL_dK = 0.5 * (tdot(alpha) - Y.shape[1] * Wi)
log_marginal = 0.5*(-Y.size * log_2_pi - Y.shape[1] * W_logdet - np.sum(alpha * YYT_factor.T)
log_marginal = 0.5*(-Y.size * log_2_pi - Y.shape[1] * W_logdet - np.sum(alpha * YYT_factor))
dL_dtheta_lik = likelihood.dL_dtheta(np.diag(dL_dK))
dL_dtheta_lik = likelihood._gradients(np.diag(dL_dK))
return Posterior(log_marginal, dL_DK, dL_dtheta_lik, LW, alpha, K):
return Posterior(log_marginal, dL_dK, dL_dtheta_lik, LW, alpha, K)

22
GPy/inference/latent_function_inference/posterior.py
View file

@ -6,27 +6,27 @@ import numpy as np
class Posterior(object):
"""
An object to represent a Gaussian posterior over latent function values.
this may be computed exactly for Gaussian likelihoods, or approximated for
This may be computed exactly for Gaussian likelihoods, or approximated for
non-Gaussian likelihoods.
The purpose of this clas is to serve as an interface between the inference
The purpose of this class is to serve as an interface between the inference
schemes and the model classes.
"""
def __init__(self, log_marginal, dLM_DK, dLM_dtheta_lik, woodbury_chol, woodbury_mean, K):
def __init__(self, log_marginal, dL_dK, dL_dtheta_lik, woodbury_chol, woodbury_vector, K):
"""
log_marginal: log p(Y|X)
DLM_dK: d/dK log p(Y|X)
dLM_dtheta_lik : d/dtheta log p(Y|X) (where theta are the parameters of the likelihood)
DL_dK: d/dK log p(Y|X)
dL_dtheta_lik : d/dtheta log p(Y|X) (where theta are the parameters of the likelihood)
woodbury_chol : a lower triangular matrix L that satisfies posterior_covariance = K - K L^{-T} L^{-1} K
woodbury_mean : a matrix (or vector, as Nx1 matrix) M which satisfies posterior_mean = K M
woodbury_vector : a matrix (or vector, as Nx1 matrix) M which satisfies posterior_mean = K M
K : the proir covariance (required for lazy computation of various quantities)
"""
self.log_marginal = log_marginal
self.dLM_DK = dLM_DK
self.dLM_dtheta_lik = _dLM_dtheta_lik
self.dL_dK = dL_dK
self.dL_dtheta_lik = dL_dtheta_lik
self._woodbury_chol = woodbury_chol
self._woodbury_mean = woodbury_mean
self._woodbury_vector = woodbury_vector
self._K = K
#these are computed lazily below
@ -37,13 +37,13 @@ class Posterior(object):
@property
def mean(self):
if self._mean is None:
self._mean = np.dot(self._K, self._woodbury_mean)
self._mean = np.dot(self._K, self._woodbury_vector)
return self._mean
@property
def covariance(self):
if self._covariance is None:
LiK, _ = dpotrs
LiK, _ = dpotrs(self._woodbury_chol, self._K)
self._covariance = self._K - tdot(LiK.T)
return self._covariance

2
GPy/inference/optimization/__init__.py
View file

@ -0,0 +1,2 @@
from scg import SCG
from optimization import *

3
GPy/likelihoods/__init__.py
View file

@ -3,4 +3,5 @@ from exponential import Exponential
from gaussian import Gaussian
from gamma import Gamma
from poisson import Poisson
from student_t import Student_t
from student_t import StudentT
from likelihood import Likelihood

46
GPy/likelihoods/bernoulli.py
View file

@ -2,10 +2,10 @@
# Licensed under the BSD 3-clause license (see LICENSE.txt)
import numpy as np
from scipy import stats,special
from scipy import stats, special
import scipy as sp
from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
import gp_transformations
from GPy.util.univariate_Gaussian import std_norm_pdf, std_norm_cdf
import link_functions
from likelihood import Likelihood
class Bernoulli(Likelihood):
@ -16,29 +16,29 @@ class Bernoulli(Likelihood):
p(y_{i}|\\lambda(f_{i})) = \\lambda(f_{i})^{y_{i}}(1-f_{i})^{1-y_{i}}
.. Note::
Y is expected to take values in {-1,1}
Y is expected to take values in {-1, 1}
Probit likelihood usually used
"""
def __init__(self,gp_link=None,analytical_mean=False,analytical_variance=False):
super(Bernoulli, self).__init__(gp_link,analytical_mean,analytical_variance)
if isinstance(gp_link , (gp_transformations.Heaviside, gp_transformations.Probit)):
def __init__(self, gp_link=None, analytical_mean=False, analytical_variance=False):
super(Bernoulli, self).__init__(gp_link, analytical_mean, analytical_variance)
if isinstance(gp_link , (link_functions.Heaviside, link_functions.Probit)):
self.log_concave = True
def _preprocess_values(self,Y):
def _preprocess_values(self, Y):
"""
Check if the values of the observations correspond to the values
assumed by the likelihood function.
..Note:: Binary classification algorithm works better with classes {-1,1}
..Note:: Binary classification algorithm works better with classes {-1, 1}
"""
Y_prep = Y.copy()
Y1 = Y[Y.flatten()==1].size
Y2 = Y[Y.flatten()==0].size
assert Y1 + Y2 == Y.size, 'Bernoulli likelihood is meant to be used only with outputs in {0,1}.'
assert Y1 + Y2 == Y.size, 'Bernoulli likelihood is meant to be used only with outputs in {0, 1}.'
Y_prep[Y.flatten() == 0] = -1
return Y_prep
def _moments_match_analytical(self,data_i,tau_i,v_i):
def _moments_match_analytical(self, data_i, tau_i, v_i):
"""
Moments match of the marginal approximation in EP algorithm
@ -51,15 +51,15 @@ class Bernoulli(Likelihood):
elif data_i == 0:
sign = -1
else:
raise ValueError("bad value for Bernouilli observation (0,1)")
if isinstance(self.gp_link,gp_transformations.Probit):
raise ValueError("bad value for Bernouilli observation (0, 1)")
if isinstance(self.gp_link, link_functions.Probit):
z = sign*v_i/np.sqrt(tau_i**2 + tau_i)
Z_hat = std_norm_cdf(z)
phi = std_norm_pdf(z)
mu_hat = v_i/tau_i + sign*phi/(Z_hat*np.sqrt(tau_i**2 + tau_i))
sigma2_hat = 1./tau_i - (phi/((tau_i**2+tau_i)*Z_hat))*(z+phi/Z_hat)
elif isinstance(self.gp_link,gp_transformations.Heaviside):
elif isinstance(self.gp_link, link_functions.Heaviside):
a = sign*v_i/np.sqrt(tau_i)
Z_hat = std_norm_cdf(a)
N = std_norm_pdf(a)
@ -68,24 +68,24 @@ class Bernoulli(Likelihood):
if np.any(np.isnan([Z_hat, mu_hat, sigma2_hat])):
stop
else:
raise ValueError("Exact moment matching not available for link {}".format(self.gp_link.gp_transformations.__name__))
raise ValueError("Exact moment matching not available for link {}".format(self.gp_link.__name__))
return Z_hat, mu_hat, sigma2_hat
def _predictive_mean_analytical(self,mu,variance):
def _predictive_mean_analytical(self, mu, variance):
if isinstance(self.gp_link,gp_transformations.Probit):
if isinstance(self.gp_link, link_functions.Probit):
return stats.norm.cdf(mu/np.sqrt(1+variance))
elif isinstance(self.gp_link,gp_transformations.Heaviside):
elif isinstance(self.gp_link, link_functions.Heaviside):
return stats.norm.cdf(mu/np.sqrt(variance))
else:
raise NotImplementedError
def _predictive_variance_analytical(self,mu,variance, pred_mean):
def _predictive_variance_analytical(self, mu, variance, pred_mean):
if isinstance(self.gp_link,gp_transformations.Heaviside):
if isinstance(self.gp_link, link_functions.Heaviside):
return 0.
else:
raise NotImplementedError
@ -106,7 +106,7 @@ class Bernoulli(Likelihood):
:rtype: float
.. Note:
Each y_i must be in {0,1}
Each y_i must be in {0, 1}
"""
assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape
objective = (link_f**y) * ((1.-link_f)**(1.-y))
@ -195,13 +195,13 @@ class Bernoulli(Likelihood):
d3logpdf_dlink3 = 2*(y/(link_f**3) - (1-y)/((1-link_f)**3))
return d3logpdf_dlink3
def _mean(self,gp):
def _mean(self, gp):
"""
Mass (or density) function
"""
return self.gp_link.transf(gp)
def _variance(self,gp):
def _variance(self, gp):
"""
Mass (or density) function
"""

4
GPy/likelihoods/exponential.py
View file

@ -6,10 +6,10 @@ import numpy as np
from scipy import stats,special
import scipy as sp
from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
import gp_transformations
import link_functions
from likelihood import Likelihood
class Exponential(NoiseDistribution):
class Exponential(Likelihood):
"""
Expoential likelihood
Y is expected to take values in {0,1,2,...}

4
GPy/likelihoods/gamma.py
View file

@ -6,10 +6,10 @@ import numpy as np
from scipy import stats,special
import scipy as sp
from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
import gp_transformations
import link_functions
from likelihood import Likelihood
class Gamma(NoiseDistribution):
class Gamma(Likelihood):
"""
Gamma likelihood

90
GPy/likelihoods/gaussian.py
View file

@ -2,13 +2,13 @@
# Licensed under the BSD 3-clause license (see LICENSE.txt)
import numpy as np
from scipy import stats,special
import scipy as sp
from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
import gp_transformations
from scipy import stats, special
from GPy.util.univariate_Gaussian import std_norm_pdf, std_norm_cdf
import link_functions
from likelihood import Likelihood
from ..core.parameter import Param
class Gaussian(NoiseDistribution):
class Gaussian(Likelihood):
"""
Gaussian likelihood
@ -19,40 +19,39 @@ class Gaussian(NoiseDistribution):
:param N: Number of data points
:type N: int
"""
def __init__(self,gp_link=None,analytical_mean=False,analytical_variance=False,variance=1., D=None, N=None):
self.variance = variance
self.N = N
self._set_params(np.asarray(variance))
super(Gaussian, self).__init__(gp_link,analytical_mean,analytical_variance)
if isinstance(gp_link , gp_transformations.Identity):
def __init__(self, gp_link=None, variance=1., name='Gaussian_noise'):
if gp_link is None:
gp_link = link_functions.Identity()
if isinstance(gp_link, link_functions.Identity):
analytical_variance = True
analytical_mean = True
else:
analytical_variance = False
analytical_mean = False
super(Gaussian, self).__init__(gp_link, analytical_mean, analytical_variance, name=name)
self.variance = Param('variance', variance)
self.add_parameter(self.variance)
if isinstance(gp_link , link_functions.Identity):
self.log_concave = True
def _get_params(self):
return np.array([self.variance])
def covariance_matrix(self, Y, Y_metadata=None):
return np.eye(Y.shape[0]) * self.variance
def _get_param_names(self):
return ['noise_model_variance']
def _gradients(self, partial):
return np.sum(partial)
def _set_params(self, p):
self.variance = float(p)
self.I = np.eye(self.N)
self.covariance_matrix = self.I * self.variance
self.Ki = self.I*(1.0 / self.variance)
#self.ln_det_K = np.sum(np.log(np.diag(self.covariance_matrix)))
self.ln_det_K = self.N*np.log(self.variance)
def _gradients(self,partial):
return np.zeros(1)
#return np.sum(partial)
def _preprocess_values(self,Y):
def _preprocess_values(self, Y):
"""
Check if the values of the observations correspond to the values
assumed by the likelihood function.
"""
return Y
def _moments_match_analytical(self,data_i,tau_i,v_i):
def _moments_match_analytical(self, data_i, tau_i, v_i):
"""
Moments match of the marginal approximation in EP algorithm
@ -66,36 +65,13 @@ class Gaussian(NoiseDistribution):
Z_hat = 1./np.sqrt(2.*np.pi*sum_var)*np.exp(-.5*(data_i - v_i/tau_i)**2./sum_var)
return Z_hat, mu_hat, sigma2_hat
def _predictive_mean_analytical(self,mu,sigma):
new_sigma2 = self.predictive_variance(mu,sigma)
def _predictive_mean_analytical(self, mu, sigma):
new_sigma2 = self.predictive_variance(mu, sigma)
return new_sigma2*(mu/sigma**2 + self.gp_link.transf(mu)/self.variance)
def _predictive_variance_analytical(self,mu,sigma,predictive_mean=None):
def _predictive_variance_analytical(self, mu, sigma, predictive_mean=None):
return 1./(1./self.variance + 1./sigma**2)
def _mass(self, link_f, y, extra_data=None):
NotImplementedError("Deprecated, now doing chain in noise_model.py for link function evaluation\
Please negate your function and use pdf in noise_model.py, if implementing a likelihood\
rederivate the derivative without doing the chain and put in logpdf, dlogpdf_dlink or\
its derivatives")
def _nlog_mass(self, link_f, y, extra_data=None):
NotImplementedError("Deprecated, now doing chain in noise_model.py for link function evaluation\
Please negate your function and use logpdf in noise_model.py, if implementing a likelihood\
rederivate the derivative without doing the chain and put in logpdf, dlogpdf_dlink or\
its derivatives")
def _dnlog_mass_dgp(self, link_f, y, extra_data=None):
NotImplementedError("Deprecated, now doing chain in noise_model.py for link function evaluation\
Please negate your function and use dlogpdf_df in noise_model.py, if implementing a likelihood\
rederivate the derivative without doing the chain and put in logpdf, dlogpdf_dlink or\
its derivatives")
def _d2nlog_mass_dgp2(self, link_f, y, extra_data=None):
NotImplementedError("Deprecated, now doing chain in noise_model.py for link function evaluation\
Please negate your function and use d2logpdf_df2 in noise_model.py, if implementing a likelihood\
rederivate the derivative without doing the chain and put in logpdf, dlogpdf_dlink or\
its derivatives")
def pdf_link(self, link_f, y, extra_data=None):
"""
Likelihood function given link(f)
@ -270,7 +246,7 @@ class Gaussian(NoiseDistribution):
d2logpdf_dlink2_dvar = self.d2logpdf_dlink2_dvar(f, y, extra_data=extra_data)
return d2logpdf_dlink2_dvar
def _mean(self,gp):
def _mean(self, gp):
"""
Expected value of y under the Mass (or density) function p(y|f)
@ -279,7 +255,7 @@ class Gaussian(NoiseDistribution):
"""
return self.gp_link.transf(gp)
def _variance(self,gp):
def _variance(self, gp):
"""
Variance of y under the Mass (or density) function p(y|f)

29
GPy/likelihoods/likelihood.py
View file

@ -1,19 +1,19 @@
# Copyright (c) 2012, 2013 Ricardo Andrade
# Licensed under the BSD 3-clause license (see LICENSE.txt)
import numpy as np
from scipy import stats,special
import scipy as sp
import pylab as pb
from GPy.util.plot import gpplot
from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
import gp_transformations
from GPy.util.misc import chain_1, chain_2, chain_3
from ..util.plot import gpplot
from ..util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
import link_functions
from ..util.misc import chain_1, chain_2, chain_3
from scipy.integrate import quad
import warnings
from ..core.parameterized import Parameterized
class Likelihood(object):
class Likelihood(Parameterized):
"""
Likelihood base class
@ -21,8 +21,12 @@ class Likelihood(object):
The minimum required funciotnality is... TODO
"""
def __init__(self,gp_link,analytical_mean=False,analytical_variance=False):
assert isinstance(gp_link,gp_transformations.GPTransformation), "gp_link is not a valid GPTransformation."
def __init__(self,gp_link,analytical_mean=False,analytical_variance=False, name='likelihood_base'):
"""
What are analytical_mean, analyitical_variance? TODO
"""
super(Likelihood, self).__init__(name)
assert isinstance(gp_link,link_functions.GPTransformation), "gp_link is not a valid GPTransformation."
self.gp_link = gp_link
self.analytical_mean = analytical_mean
self.analytical_variance = analytical_variance
@ -39,15 +43,6 @@ class Likelihood(object):
self.log_concave = False
def _get_params(self):
return np.zeros(0)
def _get_param_names(self):
return []
def _set_params(self,p):
pass
def _gradients(self,partial):
return np.zeros(0)

2
GPy/likelihoods/poisson.py
View file

@ -6,7 +6,7 @@ import numpy as np
from scipy import stats,special
import scipy as sp
from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
import gp_transformations
import link_functions
from likelihood import Likelihood
class Poisson(Likelihood):

2
GPy/likelihoods/student_t.py
View file

@ -4,7 +4,7 @@
import numpy as np
from scipy import stats, special
import scipy as sp
import gp_transformations
import link_functions
from scipy import stats, integrate
from scipy.special import gammaln, gamma
from likelihood import Likelihood

13
GPy/models/gp_regression.py
View file

@ -16,28 +16,21 @@ class GPRegression(GP):
:param X: input observations
:param Y: observed values
:param kernel: a GPy kernel, defaults to rbf
:param normalize_X: whether to normalize the input data before computing (predictions will be in original scales)
:type normalize_X: False|True
:param normalize_Y: whether to normalize the input data before computing (predictions will be in original scales)
:type normalize_Y: False|True
.. Note:: Multiple independent outputs are allowed using columns of Y
"""
def __init__(self, X, Y, kernel=None, normalize_X=False, normalize_Y=False):
def __init__(self, X, Y, kernel=None):
if kernel is None:
kernel = kern.rbf(X.shape[1])
likelihood = likelihoods.Gaussian(Y, normalize=normalize_Y)
likelihood = likelihoods.Gaussian()
super(GPRegression, self).__init__(X, likelihood, kernel, normalize_X=normalize_X)
self.ensure_default_constraints()
super(GPRegression, self).__init__(X, Y, kernel, likelihood, name='gp regression')
def getstate(self):
return GP.getstate(self)
def setstate(self, state):
return GP.setstate(self, state)
pass

77
GPy/util/misc.py
View file

@ -3,6 +3,34 @@
import numpy as np
from scipy import weave
from config import *
def chain_1(df_dg, dg_dx):
"""
Generic chaining function for first derivative
.. math::
\\frac{d(f . g)}{dx} = \\frac{df}{dg} \\frac{dg}{dx}
"""
return df_dg * dg_dx
def chain_2(d2f_dg2, dg_dx, df_dg, d2g_dx2):
"""
Generic chaining function for second derivative
.. math::
\\frac{d^{2}(f . g)}{dx^{2}} = \\frac{d^{2}f}{dg^{2}}(\\frac{dg}{dx})^{2} + \\frac{df}{dg}\\frac{d^{2}g}{dx^{2}}
"""
return d2f_dg2*(dg_dx**2) + df_dg*d2g_dx2
def chain_3(d3f_dg3, dg_dx, d2f_dg2, d2g_dx2, df_dg, d3g_dx3):
"""
Generic chaining function for third derivative
.. math::
\\frac{d^{3}(f . g)}{dx^{3}} = \\frac{d^{3}f}{dg^{3}}(\\frac{dg}{dx})^{3} + 3\\frac{d^{2}f}{dg^{2}}\\frac{dg}{dx}\\frac{d^{2}g}{dx^{2}} + \\frac{df}{dg}\\frac{d^{3}g}{dx^{3}}
"""
return d3f_dg3*(dg_dx**3) + 3*d2f_dg2*dg_dx*d2g_dx2 + df_dg*d3g_dx3
def opt_wrapper(m, **kwargs):
"""
@ -57,11 +85,18 @@ def kmm_init(X, m = 10):
return X[inducing]
def fast_array_equal(A, B):
if config.getboolean('parallel', 'openmp'):
pragma_string = '#pragma omp parallel for private(i, j)'
else:
pragma_string = ''
code2="""
int i, j;
return_val = 1;
// #pragma omp parallel for private(i, j)
%s
for(i=0;i<N;i++){
for(j=0;j<D;j++){
if(A(i, j) != B(i, j)){
@ -70,13 +105,18 @@ def fast_array_equal(A, B):
}
}
}
"""
""" % pragma_string
if config.getboolean('parallel', 'openmp'):
pragma_string = '#pragma omp parallel for private(i, j, z)'
else:
pragma_string = ''
code3="""
int i, j, z;
return_val = 1;
// #pragma omp parallel for private(i, j, z)
%s
for(i=0;i<N;i++){
for(j=0;j<D;j++){
for(z=0;z<Q;z++){
@ -87,20 +127,33 @@ def fast_array_equal(A, B):
}
}
}
"""
""" % pragma_string
if config.getboolean('parallel', 'openmp'):
pragma_string = '#include <omp.h>'
else:
pragma_string = ''
support_code = """
// #include <omp.h>
%s
#include <math.h>
"""
""" % pragma_string
weave_options = {'headers' : ['<omp.h>'],
'extra_compile_args': ['-fopenmp -O3'],
'extra_link_args' : ['-lgomp']}
weave_options_openmp = {'headers' : ['<omp.h>'],
'extra_compile_args': ['-fopenmp -O3'],
'extra_link_args' : ['-lgomp'],
'libraries': ['gomp']}
weave_options_noopenmp = {'extra_compile_args': ['-O3']}
if config.getboolean('parallel', 'openmp'):
weave_options = weave_options_openmp
else:
weave_options = weave_options_noopenmp
value = False
if (A == None) and (B == None):
return True
elif ((A == None) and (B != None)) or ((A != None) and (B == None)):
@ -110,14 +163,12 @@ def fast_array_equal(A, B):
N, D = [int(i) for i in A.shape]
value = weave.inline(code2, support_code=support_code,
arg_names=['A', 'B', 'N', 'D'],
type_converters=weave.converters.blitz)
# libraries=['gomp'], **weave_options)
type_converters=weave.converters.blitz, **weave_options)
elif A.ndim == 3:
N, D, Q = [int(i) for i in A.shape]
value = weave.inline(code3, support_code=support_code,
arg_names=['A', 'B', 'N', 'D', 'Q'],
type_converters=weave.converters.blitz)
#libraries=['gomp'], **weave_options)
type_converters=weave.converters.blitz, **weave_options)
else:
value = np.array_equal(A,B)