mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-06-11 15:15:15 +02:00
Merge remote-tracking branch 'origin/devel' into devel
This commit is contained in:
Ricardo
commit
1d4c64e17f
8 changed files with 81 additions and 31 deletions
|
|
@ -237,7 +237,7 @@ class Model(Parameterised):
|
|||
try:
|
||||
self._set_params_transformed(x)
|
||||
self._fail_count = 0
|
||||
except (LinAlgError, ZeroDivisionError) as e:
|
||||
except (LinAlgError, ZeroDivisionError, ValueError) as e:
|
||||
if self._fail_count >= self._allowed_failures:
|
||||
raise e
|
||||
self._fail_count += 1
|
||||
|
|
@ -255,7 +255,7 @@ class Model(Parameterised):
|
|||
try:
|
||||
self._set_params_transformed(x)
|
||||
self._fail_count = 0
|
||||
except (LinAlgError, ZeroDivisionError) as e:
|
||||
except (LinAlgError, ZeroDivisionError, ValueError) as e:
|
||||
if self._fail_count >= self._allowed_failures:
|
||||
raise e
|
||||
self._fail_count += 1
|
||||
|
|
@ -267,7 +267,7 @@ class Model(Parameterised):
|
|||
self._set_params_transformed(x)
|
||||
obj_f = -self.log_likelihood() - self.log_prior()
|
||||
self._fail_count = 0
|
||||
except (LinAlgError, ZeroDivisionError) as e:
|
||||
except (LinAlgError, ZeroDivisionError, ValueError) as e:
|
||||
if self._fail_count >= self._allowed_failures:
|
||||
raise e
|
||||
self._fail_count += 1
|
||||
|
|
|
|||
|
|
@ -63,6 +63,7 @@ class SparseGP(GPBase):
|
|||
|
||||
def _computations(self):
|
||||
|
||||
|
||||
# factor Kmm
|
||||
self.Lm = jitchol(self.Kmm)
|
||||
|
||||
|
|
@ -89,17 +90,18 @@ class SparseGP(GPBase):
|
|||
self.B = np.eye(self.num_inducing) + self.A
|
||||
self.LB = jitchol(self.B)
|
||||
|
||||
# TODO: make a switch for either first compute psi1V, or VV.T
|
||||
self.psi1V = np.dot(self.psi1.T, self.likelihood.V)
|
||||
#VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency!
|
||||
self.psi1Vf = np.dot(self.psi1.T, self.likelihood.VVT_factor)
|
||||
|
||||
# back substutue C into psi1V
|
||||
tmp, info1 = dtrtrs(self.Lm, np.asfortranarray(self.psi1V), lower=1, trans=0)
|
||||
self._LBi_Lmi_psi1V, _ = dtrtrs(self.LB, np.asfortranarray(tmp), lower=1, trans=0)
|
||||
# back substutue C into psi1Vf
|
||||
tmp, info1 = dtrtrs(self.Lm, np.asfortranarray(self.psi1Vf), lower=1, trans=0)
|
||||
self._LBi_Lmi_psi1Vf, _ = dtrtrs(self.LB, np.asfortranarray(tmp), lower=1, trans=0)
|
||||
tmp, info2 = dpotrs(self.LB, tmp, lower=1)
|
||||
self.Cpsi1V, info3 = dtrtrs(self.Lm, tmp, lower=1, trans=1)
|
||||
self.Cpsi1Vf, info3 = dtrtrs(self.Lm, tmp, lower=1, trans=1)
|
||||
|
||||
# Compute dL_dKmm
|
||||
tmp = tdot(self._LBi_Lmi_psi1V)
|
||||
tmp = tdot(self._LBi_Lmi_psi1Vf)
|
||||
self.data_fit = np.trace(tmp)
|
||||
self.DBi_plus_BiPBi = backsub_both_sides(self.LB, self.output_dim * np.eye(self.num_inducing) + tmp)
|
||||
tmp = -0.5 * self.DBi_plus_BiPBi
|
||||
tmp += -0.5 * self.B * self.output_dim
|
||||
|
|
@ -108,7 +110,7 @@ class SparseGP(GPBase):
|
|||
|
||||
# Compute dL_dpsi # FIXME: this is untested for the heterscedastic + uncertain inputs case
|
||||
self.dL_dpsi0 = -0.5 * self.output_dim * (self.likelihood.precision * np.ones([self.num_data, 1])).flatten()
|
||||
self.dL_dpsi1 = np.dot(self.Cpsi1V, self.likelihood.V.T).T
|
||||
self.dL_dpsi1 = np.dot(self.likelihood.VVT_factor, self.Cpsi1Vf.T)
|
||||
dL_dpsi2_beta = 0.5 * backsub_both_sides(self.Lm, self.output_dim * np.eye(self.num_inducing) - self.DBi_plus_BiPBi)
|
||||
|
||||
if self.likelihood.is_heteroscedastic:
|
||||
|
|
@ -138,18 +140,18 @@ class SparseGP(GPBase):
|
|||
# likelihood is not heterscedatic
|
||||
self.partial_for_likelihood = -0.5 * self.num_data * self.output_dim * self.likelihood.precision + 0.5 * self.likelihood.trYYT * self.likelihood.precision ** 2
|
||||
self.partial_for_likelihood += 0.5 * self.output_dim * (self.psi0.sum() * self.likelihood.precision ** 2 - np.trace(self.A) * self.likelihood.precision)
|
||||
self.partial_for_likelihood += self.likelihood.precision * (0.5 * np.sum(self.A * self.DBi_plus_BiPBi) - np.sum(np.square(self._LBi_Lmi_psi1V)))
|
||||
self.partial_for_likelihood += self.likelihood.precision * (0.5 * np.sum(self.A * self.DBi_plus_BiPBi) - self.data_fit)
|
||||
|
||||
def log_likelihood(self):
|
||||
""" Compute the (lower bound on the) log marginal likelihood """
|
||||
if self.likelihood.is_heteroscedastic:
|
||||
A = -0.5 * self.num_data * self.output_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.likelihood.precision)) - 0.5 * np.sum(self.likelihood.V * self.likelihood.Y)
|
||||
A = -0.5 * self.num_data * self.output_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.likelihood.precision)) - 0.5 * np.sum(self.likelihood.V*self.likelihood.Y)
|
||||
B = -0.5 * self.output_dim * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A))
|
||||
else:
|
||||
A = -0.5 * self.num_data * self.output_dim * (np.log(2.*np.pi) - np.log(self.likelihood.precision)) - 0.5 * self.likelihood.precision * self.likelihood.trYYT
|
||||
B = -0.5 * self.output_dim * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A))
|
||||
C = -self.output_dim * (np.sum(np.log(np.diag(self.LB)))) # + 0.5 * self.num_inducing * np.log(sf2))
|
||||
D = 0.5 * np.sum(np.square(self._LBi_Lmi_psi1V))
|
||||
D = 0.5 * self.data_fit
|
||||
return A + B + C + D + self.likelihood.Z
|
||||
|
||||
def _set_params(self, p):
|
||||
|
|
@ -158,6 +160,7 @@ class SparseGP(GPBase):
|
|||
self.likelihood._set_params(p[self.Z.size + self.kern.num_params:])
|
||||
self._compute_kernel_matrices()
|
||||
self._computations()
|
||||
self.Cpsi1V = None
|
||||
|
||||
def _get_params(self):
|
||||
return np.hstack([self.Z.flatten(), self.kern._get_params_transformed(), self.likelihood._get_params()])
|
||||
|
|
@ -224,6 +227,14 @@ class SparseGP(GPBase):
|
|||
symmetrify(Bi)
|
||||
Kmmi_LmiBLmi = backsub_both_sides(self.Lm, np.eye(self.num_inducing) - Bi)
|
||||
|
||||
if self.Cpsi1V is None:
|
||||
psi1V = np.dot(self.psi1.T,self.likelihood.V)
|
||||
tmp, _ = dtrtrs(self.Lm, np.asfortranarray(psi1V), lower=1, trans=0)
|
||||
tmp, _ = dpotrs(self.LB, tmp, lower=1)
|
||||
self.Cpsi1V, _ = dtrtrs(self.Lm, tmp, lower=1, trans=1)
|
||||
|
||||
|
||||
|
||||
if X_variance_new is None:
|
||||
Kx = self.kern.K(self.Z, Xnew, which_parts=which_parts)
|
||||
mu = np.dot(Kx.T, self.Cpsi1V)
|
||||
|
|
|
|||
|
|
@ -2,7 +2,7 @@
|
|||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
import numpy as np
|
||||
from matplotlib import pyplot as plt
|
||||
from matplotlib import pyplot as plt, cm
|
||||
|
||||
import GPy
|
||||
from GPy.core.transformations import logexp
|
||||
|
|
@ -182,7 +182,7 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False):
|
|||
sS = sS(x)
|
||||
|
||||
S1 = np.hstack([s1, sS])
|
||||
S2 = np.hstack([s2, sS])
|
||||
S2 = np.hstack([s2, s3, sS])
|
||||
S3 = np.hstack([s3, sS])
|
||||
|
||||
Y1 = S1.dot(np.random.randn(S1.shape[1], D1))
|
||||
|
|
@ -216,7 +216,7 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False):
|
|||
ax.legend()
|
||||
for i, Y in enumerate(Ylist):
|
||||
ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i)
|
||||
ax.imshow(Y)
|
||||
ax.imshow(Y, aspect='auto', cmap=cm.gray) # @UndefinedVariable
|
||||
ax.set_title("Y{}".format(i + 1))
|
||||
pylab.draw()
|
||||
pylab.tight_layout()
|
||||
|
|
@ -261,6 +261,7 @@ def bgplvm_simulation(optimize='scg',
|
|||
|
||||
k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q)
|
||||
m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k, _debug=True)
|
||||
|
||||
# m.constrain('variance|noise', logexp_clipped())
|
||||
m.ensure_default_constraints()
|
||||
m['noise'] = Y.var() / 100.
|
||||
|
|
@ -298,7 +299,7 @@ def mrd_simulation(optimize=True, plot=True, plot_sim=True, **kw):
|
|||
|
||||
if optimize:
|
||||
print "Optimizing Model:"
|
||||
m.optimize('scg', messages=1, max_iters=8e3, max_f_eval=8e3, gtol=.1)
|
||||
m.optimize(messages=1, max_iters=8e3, max_f_eval=8e3, gtol=.1)
|
||||
if plot:
|
||||
m.plot_X_1d("MRD Latent Space 1D")
|
||||
m.plot_scales("MRD Scales")
|
||||
|
|
@ -327,28 +328,56 @@ def brendan_faces():
|
|||
data_show = GPy.util.visualize.image_show(y[None, :], dimensions=(20, 28), transpose=True, invert=False, scale=False)
|
||||
lvm_visualizer = GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
|
||||
raw_input('Press enter to finish')
|
||||
lvm_visualizer.close()
|
||||
|
||||
return m
|
||||
def stick_play(range=None, frame_rate=15):
|
||||
data = GPy.util.datasets.stick()
|
||||
# optimize
|
||||
if range==None:
|
||||
Y = data['Y'].copy()
|
||||
else:
|
||||
Y = data['Y'][range[0]:range[1], :].copy()
|
||||
y = Y[0, :]
|
||||
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
|
||||
GPy.util.visualize.data_play(Y, data_show, frame_rate)
|
||||
return Y
|
||||
|
||||
def stick():
|
||||
data = GPy.util.datasets.stick()
|
||||
m = GPy.models.GPLVM(data['Y'], 2)
|
||||
|
||||
# optimize
|
||||
m = GPy.models.GPLVM(data['Y'], 2)
|
||||
m.ensure_default_constraints()
|
||||
m.optimize(messages=1, max_f_eval=10000)
|
||||
m._set_params(m._get_params())
|
||||
|
||||
plt.clf
|
||||
ax = m.plot_latent()
|
||||
y = m.likelihood.Y[0, :]
|
||||
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
|
||||
lvm_visualizer = GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
|
||||
raw_input('Press enter to finish')
|
||||
lvm_visualizer.close()
|
||||
|
||||
return m
|
||||
|
||||
def stick_bgplvm(model=None):
|
||||
data = GPy.util.datasets.stick()
|
||||
Q = 6
|
||||
kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2))
|
||||
m = BayesianGPLVM(data['Y'], Q, init="PCA", num_inducing=20,kernel=kernel)
|
||||
# optimize
|
||||
m.ensure_default_constraints()
|
||||
m.optimize(messages=1, max_f_eval=3000,xtol=1e-300,ftol=1e-300)
|
||||
m._set_params(m._get_params())
|
||||
plt.clf, (latent_axes, sense_axes) = plt.subplots(1, 2)
|
||||
plt.sca(latent_axes)
|
||||
m.plot_latent()
|
||||
y = m.likelihood.Y[0, :].copy()
|
||||
data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect'])
|
||||
lvm_visualizer = GPy.util.visualize.lvm_dimselect(m.X[0, :].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
|
||||
raw_input('Press enter to finish')
|
||||
|
||||
return m
|
||||
|
||||
|
||||
def cmu_mocap(subject='35', motion=['01'], in_place=True):
|
||||
|
||||
data = GPy.util.datasets.cmu_mocap(subject, motion)
|
||||
|
|
|
|||
|
|
@ -21,7 +21,7 @@ class Brownian(Kernpart):
|
|||
def __init__(self,input_dim,variance=1.):
|
||||
self.input_dim = input_dim
|
||||
assert self.input_dim==1, "Brownian motion in 1D only"
|
||||
self.num_params = 1.
|
||||
self.num_params = 1
|
||||
self.name = 'Brownian'
|
||||
self._set_params(np.array([variance]).flatten())
|
||||
|
||||
|
|
|
|||
|
|
@ -34,6 +34,8 @@ class EP(likelihood):
|
|||
self.Z = 0
|
||||
self.YYT = None
|
||||
self.V = self.precision * self.Y
|
||||
self.VVT_factor = self.V
|
||||
self.trYYT = 0.
|
||||
|
||||
def restart(self):
|
||||
self.tau_tilde = np.zeros(self.N)
|
||||
|
|
@ -44,6 +46,8 @@ class EP(likelihood):
|
|||
self.Z = 0
|
||||
self.YYT = None
|
||||
self.V = self.precision * self.Y
|
||||
self.VVT_factor = self.V
|
||||
self.trYYT = 0.
|
||||
|
||||
def predictive_values(self,mu,var,full_cov):
|
||||
if full_cov:
|
||||
|
|
@ -71,6 +75,8 @@ class EP(likelihood):
|
|||
self.covariance_matrix = np.diag(1./self.tau_tilde)
|
||||
self.precision = self.tau_tilde[:,None]
|
||||
self.V = self.precision * self.Y
|
||||
self.VVT_factor = self.V
|
||||
self.trYYT = np.trace(self.YYT)
|
||||
|
||||
def fit_full(self,K):
|
||||
"""
|
||||
|
|
|
|||
|
|
@ -40,9 +40,11 @@ class Gaussian(likelihood):
|
|||
if D > self.N:
|
||||
self.YYT = np.dot(self.Y, self.Y.T)
|
||||
self.trYYT = np.trace(self.YYT)
|
||||
self.YYT_factor = jitchol(self.YYT)
|
||||
else:
|
||||
self.YYT = None
|
||||
self.trYYT = np.sum(np.square(self.Y))
|
||||
self.YYT_factor = self.Y
|
||||
|
||||
def _get_params(self):
|
||||
return np.asarray(self._variance)
|
||||
|
|
@ -53,12 +55,13 @@ class Gaussian(likelihood):
|
|||
def _set_params(self, x):
|
||||
x = np.float64(x)
|
||||
if np.all(self._variance != x):
|
||||
if x == 0.:
|
||||
if x == 0.:#special case of zero noise
|
||||
self.precision = np.inf
|
||||
self.V = None
|
||||
else:
|
||||
self.precision = 1. / x
|
||||
self.V = (self.precision) * self.Y
|
||||
self.VVT_factor = self.precision * self.YYT_factor
|
||||
self.covariance_matrix = np.eye(self.N) * x
|
||||
self._variance = x
|
||||
|
||||
|
|
|
|||
|
|
@ -193,7 +193,7 @@ class lvm_dimselect(lvm):
|
|||
GPy.examples.dimensionality_reduction.BGPVLM_oil()
|
||||
|
||||
"""
|
||||
def __init__(self, vals, Model, data_visualize, latent_axes=None, sense_axes=None, latent_index=[0, 1]):
|
||||
def __init__(self, vals, Model, data_visualize, latent_axes=None, sense_axes=None, latent_index=[0, 1], labels=None):
|
||||
if latent_axes==None and sense_axes==None:
|
||||
self.fig,(latent_axes,self.sense_axes) = plt.subplots(1,2)
|
||||
elif sense_axes==None:
|
||||
|
|
@ -201,8 +201,9 @@ class lvm_dimselect(lvm):
|
|||
self.sense_axes = fig.add_subplot(111)
|
||||
else:
|
||||
self.sense_axes = sense_axes
|
||||
|
||||
self.labels = labels
|
||||
lvm.__init__(self,vals,Model,data_visualize,latent_axes,sense_axes,latent_index)
|
||||
self.show_sensitivities()
|
||||
print "use left and right mouse butons to select dimensions"
|
||||
|
||||
|
||||
|
|
@ -221,7 +222,7 @@ class lvm_dimselect(lvm):
|
|||
|
||||
self.latent_axes.cla()
|
||||
self.Model.plot_latent(which_indices=self.latent_index,
|
||||
ax=self.latent_axes)
|
||||
ax=self.latent_axes, labels=self.labels)
|
||||
self.latent_handle = self.latent_axes.plot([0],[0],'rx',mew=2)[0]
|
||||
self.modify(self.latent_values)
|
||||
|
||||
|
|
@ -506,5 +507,5 @@ def data_play(Y, visualizer, frame_rate=30):
|
|||
|
||||
|
||||
for y in Y:
|
||||
visualizer.modify(y)
|
||||
visualizer.modify(y[None, :])
|
||||
time.sleep(1./float(frame_rate))
|
||||
|
|
|
|||
4
setup.py
4
setup.py
|
|
@ -5,7 +5,7 @@ import os
|
|||
from setuptools import setup
|
||||
|
||||
# Version number
|
||||
version = '0.3.2'
|
||||
version = '0.4.5'
|
||||
|
||||
def read(fname):
|
||||
return open(os.path.join(os.path.dirname(__file__), fname)).read()
|
||||
|
|
@ -23,7 +23,7 @@ setup(name = 'GPy',
|
|||
package_data = {'GPy': ['GPy/examples']},
|
||||
py_modules = ['GPy.__init__'],
|
||||
long_description=read('README.md'),
|
||||
install_requires=['sympy', 'numpy>=1.6', 'scipy>=0.9','matplotlib>=1.1', 'nose'],
|
||||
install_requires=['numpy>=1.6', 'scipy>=0.9','matplotlib>=1.1', 'nose'],
|
||||
extras_require = {
|
||||
'docs':['Sphinx', 'ipython'],
|
||||
},
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue