mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-27 14:25:16 +02:00
Merge branch 'devel' of github.com:SheffieldML/GPy into devel
This commit is contained in:
Teo de Campos
commit
056d68251c
19 changed files with 346 additions and 409 deletions
|
|
@ -66,7 +66,7 @@ class model(parameterised):
|
|||
|
||||
|
||||
# check constraints are okay
|
||||
if isinstance(what, (priors.gamma, priors.log_Gaussian)):
|
||||
if isinstance(what, (priors.gamma, priors.inverse_gamma, priors.log_Gaussian)):
|
||||
constrained_positive_indices = [i for i, t in zip(self.constrained_indices, self.constraints) if t.domain == 'positive']
|
||||
if len(constrained_positive_indices):
|
||||
constrained_positive_indices = np.hstack(constrained_positive_indices)
|
||||
|
|
|
|||
|
|
@ -251,7 +251,18 @@ class parameterised(object):
|
|||
|
||||
def _set_params_transformed(self, x):
|
||||
""" takes the vector x, which is then modified (by untying, reparameterising or inserting fixed values), and then call self._set_params"""
|
||||
self._set_params(self._untransform_params(x))
|
||||
|
||||
def _untransform_params(self,x):
|
||||
"""
|
||||
The transformation required for _set_params_transformed.
|
||||
|
||||
This moves the vector x seen by the optimiser (unconstrained) to the
|
||||
valid parameter vector seen by the model
|
||||
|
||||
Note:
|
||||
- This function is separate from _set_params_transformed for downstream flexibility
|
||||
"""
|
||||
# work out how many places are fixed, and where they are. tricky logic!
|
||||
fix_places = self.fixed_indices + [t[1:] for t in self.tied_indices]
|
||||
if len(fix_places):
|
||||
|
|
@ -272,7 +283,8 @@ class parameterised(object):
|
|||
[np.put(xx,i,t.f(xx[i])) for i,t in zip(self.constrained_indices, self.constraints)]
|
||||
if hasattr(self,'debug'):
|
||||
stop
|
||||
self._set_params(xx)
|
||||
|
||||
return xx
|
||||
|
||||
def _get_param_names_transformed(self):
|
||||
"""
|
||||
|
|
|
|||
|
|
@ -26,7 +26,6 @@ class Gaussian(prior):
|
|||
:param mu: mean
|
||||
:param sigma: standard deviation
|
||||
|
||||
|
||||
.. Note:: Bishop 2006 notation is used throughout the code
|
||||
|
||||
"""
|
||||
|
|
@ -144,7 +143,6 @@ def gamma_from_EV(E,V):
|
|||
b = E/V
|
||||
return gamma(a,b)
|
||||
|
||||
|
||||
class gamma(prior):
|
||||
"""
|
||||
Implementation of the Gamma probability function, coupled with random variables.
|
||||
|
|
@ -155,7 +153,6 @@ class gamma(prior):
|
|||
.. Note:: Bishop 2006 notation is used throughout the code
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self,a,b):
|
||||
self.a = float(a)
|
||||
self.b = float(b)
|
||||
|
|
@ -183,3 +180,30 @@ class gamma(prior):
|
|||
|
||||
def rvs(self,n):
|
||||
return np.random.gamma(scale=1./self.b,shape=self.a,size=n)
|
||||
|
||||
class inverse_gamma(prior):
|
||||
"""
|
||||
Implementation of the inverse-Gamma probability function, coupled with random variables.
|
||||
|
||||
:param a: shape parameter
|
||||
:param b: rate parameter (warning: it's the *inverse* of the scale)
|
||||
|
||||
.. Note:: Bishop 2006 notation is used throughout the code
|
||||
|
||||
"""
|
||||
def __init__(self,a,b):
|
||||
self.a = float(a)
|
||||
self.b = float(b)
|
||||
self.constant = -gammaln(self.a) + a*np.log(b)
|
||||
|
||||
def __str__(self):
|
||||
return "iGa("+str(np.round(self.a))+', '+str(np.round(self.b))+')'
|
||||
|
||||
def lnpdf(self,x):
|
||||
return self.constant - (self.a+1)*np.log(x) - self.b/x
|
||||
|
||||
def lnpdf_grad(self,x):
|
||||
return -(self.a+1.)/x + self.b/x**2
|
||||
|
||||
def rvs(self,n):
|
||||
return 1./np.random.gamma(scale=1./self.b,shape=self.a,size=n)
|
||||
|
|
|
|||
|
|
@ -39,8 +39,8 @@ class logexp(transformation):
|
|||
return '(+ve)'
|
||||
|
||||
class logexp_clipped(transformation):
|
||||
max_bound = 1e300
|
||||
min_bound = 1e-10
|
||||
max_bound = 1e250
|
||||
min_bound = 1e-9
|
||||
log_max_bound = np.log(max_bound)
|
||||
log_min_bound = np.log(min_bound)
|
||||
def __init__(self, lower=1e-6):
|
||||
|
|
@ -49,11 +49,13 @@ class logexp_clipped(transformation):
|
|||
def f(self, x):
|
||||
exp = np.exp(np.clip(x, self.log_min_bound, self.log_max_bound))
|
||||
f = np.log(1. + exp)
|
||||
if np.isnan(f).any():
|
||||
import ipdb;ipdb.set_trace()
|
||||
return f
|
||||
def finv(self, f):
|
||||
return np.log(np.exp(np.clip(f, self.min_bound, self.max_bound)) - 1.)
|
||||
def gradfactor(self, f):
|
||||
ef = np.exp(f)
|
||||
ef = np.exp(f) # np.clip(f, self.min_bound, self.max_bound))
|
||||
gf = (ef - 1.) / ef
|
||||
return np.where(f < self.lower, 0, gf)
|
||||
def initialize(self, f):
|
||||
|
|
|
|||
|
|
@ -9,8 +9,8 @@ import pylab as pb
|
|||
import numpy as np
|
||||
import GPy
|
||||
|
||||
default_seed=10000
|
||||
def crescent_data(seed=default_seed): #FIXME
|
||||
default_seed = 10000
|
||||
def crescent_data(seed=default_seed): # FIXME
|
||||
"""Run a Gaussian process classification on the crescent data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood.
|
||||
|
||||
:param model_type: type of model to fit ['Full', 'FITC', 'DTC'].
|
||||
|
|
@ -27,10 +27,10 @@ def crescent_data(seed=default_seed): #FIXME
|
|||
|
||||
# Likelihood object
|
||||
distribution = GPy.likelihoods.likelihood_functions.probit()
|
||||
likelihood = GPy.likelihoods.EP(data['Y'],distribution)
|
||||
likelihood = GPy.likelihoods.EP(data['Y'], distribution)
|
||||
|
||||
|
||||
m = GPy.models.GP(data['X'],likelihood,kernel)
|
||||
m = GPy.models.GP(data['X'], likelihood, kernel)
|
||||
m.ensure_default_constraints()
|
||||
|
||||
m.update_likelihood_approximation()
|
||||
|
|
@ -54,10 +54,10 @@ def oil():
|
|||
|
||||
# Likelihood object
|
||||
distribution = GPy.likelihoods.likelihood_functions.probit()
|
||||
likelihood = GPy.likelihoods.EP(data['Y'][:, 0:1],distribution)
|
||||
likelihood = GPy.likelihoods.EP(data['Y'][:, 0:1], distribution)
|
||||
|
||||
# Create GP model
|
||||
m = GPy.models.GP(data['X'],likelihood=likelihood,kernel=kernel)
|
||||
m = GPy.models.GP(data['X'], likelihood=likelihood, kernel=kernel)
|
||||
|
||||
# Contrain all parameters to be positive
|
||||
m.constrain_positive('')
|
||||
|
|
@ -85,17 +85,17 @@ def toy_linear_1d_classification(seed=default_seed):
|
|||
|
||||
# Likelihood object
|
||||
distribution = GPy.likelihoods.likelihood_functions.probit()
|
||||
likelihood = GPy.likelihoods.EP(Y,distribution)
|
||||
likelihood = GPy.likelihoods.EP(Y, distribution)
|
||||
|
||||
# Model definition
|
||||
m = GPy.models.GP(data['X'],likelihood=likelihood,kernel=kernel)
|
||||
m = GPy.models.GP(data['X'], likelihood=likelihood, kernel=kernel)
|
||||
m.ensure_default_constraints()
|
||||
|
||||
# Optimize
|
||||
m.update_likelihood_approximation()
|
||||
# Parameters optimization:
|
||||
m.optimize()
|
||||
#m.pseudo_EM() #FIXME
|
||||
# m.pseudo_EM() #FIXME
|
||||
|
||||
# Plot
|
||||
pb.subplot(211)
|
||||
|
|
@ -121,20 +121,20 @@ def sparse_toy_linear_1d_classification(seed=default_seed):
|
|||
|
||||
# Likelihood object
|
||||
distribution = GPy.likelihoods.likelihood_functions.probit()
|
||||
likelihood = GPy.likelihoods.EP(Y,distribution)
|
||||
likelihood = GPy.likelihoods.EP(Y, distribution)
|
||||
|
||||
Z = np.random.uniform(data['X'].min(),data['X'].max(),(10,1))
|
||||
Z = np.random.uniform(data['X'].min(), data['X'].max(), (10, 1))
|
||||
|
||||
# Model definition
|
||||
m = GPy.models.sparse_GP(data['X'],likelihood=likelihood,kernel=kernel,Z=Z,normalize_X=False)
|
||||
m.set('len',2.)
|
||||
m = GPy.models.sparse_GP(data['X'], likelihood=likelihood, kernel=kernel, Z=Z, normalize_X=False)
|
||||
m.set('len', 2.)
|
||||
|
||||
m.ensure_default_constraints()
|
||||
# Optimize
|
||||
m.update_likelihood_approximation()
|
||||
# Parameters optimization:
|
||||
m.optimize()
|
||||
#m.EPEM() #FIXME
|
||||
# m.EPEM() #FIXME
|
||||
|
||||
# Plot
|
||||
pb.subplot(211)
|
||||
|
|
@ -162,15 +162,15 @@ def sparse_crescent_data(inducing=10, seed=default_seed):
|
|||
|
||||
# Likelihood object
|
||||
distribution = GPy.likelihoods.likelihood_functions.probit()
|
||||
likelihood = GPy.likelihoods.EP(data['Y'],distribution)
|
||||
likelihood = GPy.likelihoods.EP(data['Y'], distribution)
|
||||
|
||||
sample = np.random.randint(0,data['X'].shape[0],inducing)
|
||||
Z = data['X'][sample,:]
|
||||
sample = np.random.randint(0, data['X'].shape[0], inducing)
|
||||
Z = data['X'][sample, :]
|
||||
|
||||
# create sparse GP EP model
|
||||
m = GPy.models.sparse_GP(data['X'],likelihood=likelihood,kernel=kernel,Z=Z)
|
||||
m = GPy.models.sparse_GP(data['X'], likelihood=likelihood, kernel=kernel, Z=Z)
|
||||
m.ensure_default_constraints()
|
||||
m.set('len',10.)
|
||||
m.set('len', 10.)
|
||||
|
||||
m.update_likelihood_approximation()
|
||||
|
||||
|
|
|
|||
|
|
@ -17,11 +17,11 @@ def BGPLVM(seed=default_seed):
|
|||
D = 4
|
||||
# generate GPLVM-like data
|
||||
X = np.random.rand(N, Q)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
k = GPy.kern.rbf(Q) + GPy.kern.white(Q, 0.00001)
|
||||
K = k.K(X)
|
||||
Y = np.random.multivariate_normal(np.zeros(N), K, D).T
|
||||
|
||||
k = GPy.kern.linear(Q, ARD=True) + GPy.kern.white(Q)
|
||||
k = GPy.kern.rbf(Q, ARD=True) + GPy.kern.linear(Q, ARD=True) + GPy.kern.rbf(Q, ARD=True) + GPy.kern.white(Q)
|
||||
# k = GPy.kern.rbf(Q) + GPy.kern.rbf(Q) + GPy.kern.white(Q)
|
||||
# k = GPy.kern.rbf(Q) + GPy.kern.bias(Q) + GPy.kern.white(Q, 0.00001)
|
||||
# k = GPy.kern.rbf(Q, ARD = False) + GPy.kern.white(Q, 0.00001)
|
||||
|
|
@ -118,9 +118,9 @@ def swiss_roll(optimize=True, N=1000, M=15, Q=4, sigma=.2, plot=False):
|
|||
return m
|
||||
|
||||
def BGPLVM_oil(optimize=True, N=100, Q=5, M=25, max_f_eval=4e3, plot=False, **k):
|
||||
np.random.seed(0)
|
||||
data = GPy.util.datasets.oil()
|
||||
from GPy.core.transformations import logexp_clipped
|
||||
np.random.seed(0)
|
||||
|
||||
# create simple GP model
|
||||
kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2))
|
||||
|
|
@ -131,8 +131,7 @@ def BGPLVM_oil(optimize=True, N=100, Q=5, M=25, max_f_eval=4e3, plot=False, **k)
|
|||
m = GPy.models.Bayesian_GPLVM(Yn, Q, kernel=kernel, M=M, **k)
|
||||
m.data_labels = data['Y'][:N].argmax(axis=1)
|
||||
|
||||
# m.constrain('variance', logexp_clipped())
|
||||
# m.constrain('length', logexp_clipped())
|
||||
m.constrain('variance|leng', logexp_clipped())
|
||||
m['lengt'] = m.X.var(0).max() / m.X.var(0)
|
||||
m['noise'] = Yn.var() / 100.
|
||||
|
||||
|
|
@ -140,10 +139,6 @@ def BGPLVM_oil(optimize=True, N=100, Q=5, M=25, max_f_eval=4e3, plot=False, **k)
|
|||
|
||||
# optimize
|
||||
if optimize:
|
||||
# m.unconstrain('noise'); m.constrain_fixed('noise')
|
||||
# m.optimize('scg', messages=1, max_f_eval=200)
|
||||
# m.unconstrain('noise')
|
||||
# m.constrain('noise', logexp_clipped())
|
||||
m.optimize('scg', messages=1, max_f_eval=max_f_eval)
|
||||
|
||||
if plot:
|
||||
|
|
@ -155,11 +150,6 @@ def BGPLVM_oil(optimize=True, N=100, Q=5, M=25, max_f_eval=4e3, plot=False, **k)
|
|||
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')
|
||||
plt.close('all')
|
||||
# # plot
|
||||
# print(m)
|
||||
# m.plot_latent(labels=m.data_labels)
|
||||
# pb.figure()
|
||||
# pb.bar(np.arange(m.kern.D), 1. / m.input_sensitivity())
|
||||
return m
|
||||
|
||||
def oil_100():
|
||||
|
|
@ -189,15 +179,6 @@ def _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim=False):
|
|||
s3 = s3(x)
|
||||
sS = sS(x)
|
||||
|
||||
# s1 -= s1.mean()
|
||||
# s2 -= s2.mean()
|
||||
# s3 -= s3.mean()
|
||||
# sS -= sS.mean()
|
||||
# s1 /= .5 * (np.abs(s1).max() - np.abs(s1).min())
|
||||
# s2 /= .5 * (np.abs(s2).max() - np.abs(s2).min())
|
||||
# s3 /= .5 * (np.abs(s3).max() - np.abs(s3).min())
|
||||
# sS /= .5 * (np.abs(sS).max() - np.abs(sS).min())
|
||||
|
||||
S1 = np.hstack([s1, sS])
|
||||
S2 = np.hstack([s2, sS])
|
||||
S3 = np.hstack([s3, sS])
|
||||
|
|
@ -217,16 +198,17 @@ def _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim=False):
|
|||
Y2 /= Y2.std(0)
|
||||
Y3 /= Y3.std(0)
|
||||
|
||||
slist = [s1, s2, s3, sS]
|
||||
slist = [sS, s1, s2, s3]
|
||||
slist_names = ["sS", "s1", "s2", "s3"]
|
||||
Ylist = [Y1, Y2, Y3]
|
||||
|
||||
if plot_sim:
|
||||
import pylab
|
||||
import itertools
|
||||
fig = pylab.figure("MRD Simulation", figsize=(8, 6))
|
||||
fig = pylab.figure("MRD Simulation Data", figsize=(8, 6))
|
||||
fig.clf()
|
||||
ax = fig.add_subplot(2, 1, 1)
|
||||
labls = sorted(filter(lambda x: x.startswith("s"), locals()))
|
||||
labls = slist_names
|
||||
for S, lab in itertools.izip(slist, labls):
|
||||
ax.plot(S, label=lab)
|
||||
ax.legend()
|
||||
|
|
@ -250,7 +232,6 @@ def bgplvm_simulation_matlab_compare():
|
|||
from GPy.models import mrd
|
||||
from GPy import kern
|
||||
reload(mrd); reload(kern)
|
||||
# k = kern.rbf(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
|
||||
k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
|
||||
m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k,
|
||||
# X=mu,
|
||||
|
|
@ -260,26 +241,14 @@ def bgplvm_simulation_matlab_compare():
|
|||
m.auto_scale_factor = True
|
||||
m['noise'] = Y.var() / 100.
|
||||
m['linear_variance'] = .01
|
||||
|
||||
# lscstr = 'X_variance'
|
||||
# m[lscstr] = .01
|
||||
# m.unconstrain(lscstr); m.constrain_fixed(lscstr, .1)
|
||||
|
||||
# cstr = 'white'
|
||||
# m.unconstrain(cstr); m.constrain_bounded(cstr, .01, 1.)
|
||||
|
||||
# cstr = 'noise'
|
||||
# m.unconstrain(cstr); m.constrain_bounded(cstr, .01, 1.)
|
||||
return m
|
||||
|
||||
def bgplvm_simulation(burnin='scg', plot_sim=False,
|
||||
max_burnin=100, true_X=False,
|
||||
do_opt=True,
|
||||
max_f_eval=1000):
|
||||
def bgplvm_simulation(optimize='scg',
|
||||
plot=True,
|
||||
max_f_eval=2e4):
|
||||
from GPy.core.transformations import logexp_clipped
|
||||
|
||||
D1, D2, D3, N, M, Q = 15, 8, 8, 350, 3, 6
|
||||
slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim)
|
||||
D1, D2, D3, N, M, Q = 15, 8, 8, 100, 3, 5
|
||||
slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot)
|
||||
|
||||
from GPy.models import mrd
|
||||
from GPy import kern
|
||||
|
|
@ -289,95 +258,23 @@ def bgplvm_simulation(burnin='scg', plot_sim=False,
|
|||
Y = Ylist[0]
|
||||
|
||||
k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q)
|
||||
# k = kern.white(Q, .00001) + kern.bias(Q)
|
||||
m = Bayesian_GPLVM(Y, Q, init="PCA", M=M, kernel=k, _debug=True)
|
||||
# m.set('noise',)
|
||||
m.constrain('variance', logexp_clipped())
|
||||
|
||||
m.ensure_default_constraints()
|
||||
m.constrain('variance|noise', logexp_clipped())
|
||||
# m.ensure_default_constraints()
|
||||
m['noise'] = Y.var() / 100.
|
||||
m['linear_variance'] = .001
|
||||
# m.auto_scale_factor = True
|
||||
# m.scale_factor = 1.
|
||||
m['linear_variance'] = .01
|
||||
|
||||
|
||||
if burnin:
|
||||
print "initializing beta"
|
||||
cstr = "noise"
|
||||
m.unconstrain(cstr); m.constrain_fixed(cstr, Y.var() / 70.)
|
||||
m.optimize(burnin, messages=1, max_f_eval=max_burnin)
|
||||
|
||||
print "releasing beta"
|
||||
cstr = "noise"
|
||||
m.unconstrain(cstr); m.constrain_positive(cstr)
|
||||
|
||||
if true_X:
|
||||
true_X = np.hstack((slist[0], slist[3], 0. * np.ones((N, Q - 2))))
|
||||
m.set('X_\d', true_X)
|
||||
m.constrain_fixed("X_\d")
|
||||
|
||||
cstr = 'X_variance'
|
||||
# m.unconstrain(cstr), m.constrain_fixed(cstr, .0001)
|
||||
m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-7, .1)
|
||||
|
||||
# cstr = 'X_variance'
|
||||
# m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-3, 1.)
|
||||
|
||||
# m['X_var'] = np.ones(N * Q) * .5 + np.random.randn(N * Q) * .01
|
||||
|
||||
# cstr = "iip"
|
||||
# m.unconstrain(cstr); m.constrain_fixed(cstr)
|
||||
|
||||
# cstr = 'variance'
|
||||
# m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-10, 1.)
|
||||
# cstr = 'X_\d'
|
||||
# m.unconstrain(cstr), m.constrain_bounded(cstr, -10., 10.)
|
||||
#
|
||||
# cstr = 'noise'
|
||||
# m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-5, 1.)
|
||||
#
|
||||
# cstr = 'white'
|
||||
# m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-6, 1.)
|
||||
#
|
||||
# cstr = 'linear_variance'
|
||||
# m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-10, 10.)
|
||||
|
||||
# cstr = 'variance'
|
||||
# m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-10, 10.)
|
||||
|
||||
# np.seterr(all='call')
|
||||
# def ipdbonerr(errtype, flags):
|
||||
# import ipdb; ipdb.set_trace()
|
||||
# np.seterrcall(ipdbonerr)
|
||||
|
||||
if do_opt and burnin:
|
||||
try:
|
||||
m.optimize(burnin, messages=1, max_f_eval=max_f_eval)
|
||||
except:
|
||||
pass
|
||||
finally:
|
||||
return m
|
||||
if optimize:
|
||||
print "Optimizing model:"
|
||||
m.optimize('scg', max_iters=max_f_eval, max_f_eval=max_f_eval, messages=True)
|
||||
if plot:
|
||||
import pylab
|
||||
m.plot_X_1d()
|
||||
pylab.figure(); pylab.axis(); m.kern.plot_ARD()
|
||||
return m
|
||||
|
||||
def mrd_simulation(plot_sim=False):
|
||||
# num = 2
|
||||
# ard1 = np.array([1., 1, 0, 0], dtype=float)
|
||||
# ard2 = np.array([0., 1, 1, 0], dtype=float)
|
||||
# ard1[ard1 == 0] = 1E-10
|
||||
# ard2[ard2 == 0] = 1E-10
|
||||
|
||||
# ard1i = 1. / ard1
|
||||
# ard2i = 1. / ard2
|
||||
|
||||
# k = GPy.kern.rbf(Q, ARD=True, lengthscale=ard1i) + GPy.kern.bias(Q, 0) + GPy.kern.white(Q, 0.0001)
|
||||
# Y1 = np.random.multivariate_normal(np.zeros(N), k.K(X), D1).T
|
||||
# Y1 -= Y1.mean(0)
|
||||
#
|
||||
# k = GPy.kern.rbf(Q, ARD=True, lengthscale=ard2i) + GPy.kern.bias(Q, 0) + GPy.kern.white(Q, 0.0001)
|
||||
# Y2 = np.random.multivariate_normal(np.zeros(N), k.K(X), D2).T
|
||||
# Y2 -= Y2.mean(0)
|
||||
# make_params = lambda ard: np.hstack([[1], ard, [1, .3]])
|
||||
D1, D2, D3, N, M, Q = 150, 250, 300, 700, 3, 7
|
||||
def mrd_simulation(optimize=True, plot_sim=False):
|
||||
D1, D2, D3, N, M, Q = 150, 250, 30, 300, 3, 7
|
||||
slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim)
|
||||
|
||||
from GPy.models import mrd
|
||||
|
|
@ -386,50 +283,23 @@ def mrd_simulation(plot_sim=False):
|
|||
|
||||
reload(mrd); reload(kern)
|
||||
|
||||
# k = kern.rbf(2, ARD=True) + kern.bias(2) + kern.white(2)
|
||||
# Y1 = np.random.multivariate_normal(np.zeros(N), k.K(S1), D1).T
|
||||
# Y2 = np.random.multivariate_normal(np.zeros(N), k.K(S2), D2).T
|
||||
# Y3 = np.random.multivariate_normal(np.zeros(N), k.K(S3), D3).T
|
||||
|
||||
# Ylist = Ylist[0:2]
|
||||
|
||||
# k = kern.rbf(Q, ARD=True) + kern.bias(Q) + kern.white(Q)
|
||||
|
||||
k = kern.linear(Q, [0.01] * Q, True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2))
|
||||
m = mrd.MRD(*Ylist, Q=Q, M=M, kernel=k, initx="concat", initz='permute', _debug=False)
|
||||
m = mrd.MRD(*Ylist, Q=Q, M=M, kernel=k, initx="concat", initz='permute')
|
||||
|
||||
for i, Y in enumerate(Ylist):
|
||||
m['{}_noise'.format(i + 1)] = Y.var() / 100.
|
||||
|
||||
m.constrain('variance', logexp_clipped())
|
||||
m.constrain('variance|noise', logexp_clipped())
|
||||
m.ensure_default_constraints()
|
||||
# m.auto_scale_factor = True
|
||||
|
||||
# cstr = 'variance'
|
||||
# m.unconstrain(cstr), m.constrain_bounded(cstr, 1e-12, 1.)
|
||||
#
|
||||
# cstr = 'linear_variance'
|
||||
# m.unconstrain(cstr), m.constrain_positive(cstr)
|
||||
# DEBUG
|
||||
np.seterr("raise")
|
||||
|
||||
print "initializing beta"
|
||||
cstr = "noise"
|
||||
m.unconstrain(cstr); m.constrain_fixed(cstr)
|
||||
m.optimize('scg', messages=1, max_f_eval=2e3, gtol=100)
|
||||
if optimize:
|
||||
print "Optimizing Model:"
|
||||
m.optimize('scg', messages=1, max_iters=3e3)
|
||||
|
||||
print "releasing beta"
|
||||
cstr = "noise"
|
||||
m.unconstrain(cstr); m.constrain(cstr, logexp_clipped())
|
||||
|
||||
# np.seterr(all='call')
|
||||
# def ipdbonerr(errtype, flags):
|
||||
# import ipdb; ipdb.set_trace()
|
||||
# np.seterrcall(ipdbonerr)
|
||||
|
||||
return m # , mtest
|
||||
|
||||
def mrd_silhouette():
|
||||
|
||||
pass
|
||||
return m
|
||||
|
||||
def brendan_faces():
|
||||
from GPy import kern
|
||||
|
|
|
|||
|
|
@ -36,8 +36,14 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
|
|||
Returns
|
||||
x the optimal value for x
|
||||
flog : a list of all the objective values
|
||||
|
||||
function_eval number of fn evaluations
|
||||
status: string describing convergence status
|
||||
"""
|
||||
|
||||
if display:
|
||||
print " SCG"
|
||||
print ' {0:{mi}s} {1:11s} {2:11s} {3:11s}'.format("I", "F", "Scale", "|g|", mi=len(str(maxiters)))
|
||||
|
||||
if xtol is None:
|
||||
xtol = 1e-6
|
||||
if ftol is None:
|
||||
|
|
@ -45,18 +51,18 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
|
|||
if gtol is None:
|
||||
gtol = 1e-5
|
||||
sigma0 = 1.0e-4
|
||||
fold = f(x, *optargs) # Initial function value.
|
||||
fold = f(x, *optargs) # Initial function value.
|
||||
function_eval = 1
|
||||
fnow = fold
|
||||
gradnew = gradf(x, *optargs) # Initial gradient.
|
||||
gradnew = gradf(x, *optargs) # Initial gradient.
|
||||
current_grad = np.dot(gradnew, gradnew)
|
||||
gradold = gradnew.copy()
|
||||
d = -gradnew # Initial search direction.
|
||||
success = True # Force calculation of directional derivs.
|
||||
nsuccess = 0 # nsuccess counts number of successes.
|
||||
beta = 1.0 # Initial scale parameter.
|
||||
betamin = 1.0e-15 # Lower bound on scale.
|
||||
betamax = 1.0e100 # Upper bound on scale.
|
||||
d = -gradnew # Initial search direction.
|
||||
success = True # Force calculation of directional derivs.
|
||||
nsuccess = 0 # nsuccess counts number of successes.
|
||||
beta = 1.0 # Initial scale parameter.
|
||||
betamin = 1.0e-15 # Lower bound on scale.
|
||||
betamax = 1.0e100 # Upper bound on scale.
|
||||
status = "Not converged"
|
||||
|
||||
flog = [fold]
|
||||
|
|
@ -106,12 +112,12 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
|
|||
fnow = fold
|
||||
|
||||
# Store relevant variables
|
||||
flog.append(fnow) # Current function value
|
||||
flog.append(fnow) # Current function value
|
||||
|
||||
iteration += 1
|
||||
if display:
|
||||
print '\r',
|
||||
print 'Iter: {0:>0{mi}g} Obj:{1:> 12e} Scale:{2:> 12e} |g|:{3:> 12e}'.format(iteration, float(fnow), float(beta), float(current_grad), mi=len(str(maxiters))),
|
||||
print '{0:>0{mi}g} {1:> 12e} {2:> 12e} {3:> 12e}'.format(iteration, float(fnow), float(beta), float(current_grad), mi=len(str(maxiters))),
|
||||
# print 'Iteration:', iteration, ' Objective:', fnow, ' Scale:', beta, '\r',
|
||||
sys.stdout.flush()
|
||||
|
||||
|
|
@ -153,5 +159,6 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto
|
|||
# iterations.
|
||||
status = "maxiter exceeded"
|
||||
|
||||
print ""
|
||||
if display:
|
||||
print ""
|
||||
return x, flog, function_eval, status
|
||||
|
|
|
|||
|
|
@ -18,7 +18,7 @@ class opt_SGD(Optimizer):
|
|||
|
||||
"""
|
||||
|
||||
def __init__(self, start, iterations = 10, learning_rate = 1e-4, momentum = 0.9, model = None, messages = False, batch_size = 1, self_paced = False, center = True, iteration_file = None, **kwargs):
|
||||
def __init__(self, start, iterations = 10, learning_rate = 1e-4, momentum = 0.9, model = None, messages = False, batch_size = 1, self_paced = False, center = True, iteration_file = None, learning_rate_adaptation=None, **kwargs):
|
||||
self.opt_name = "Stochastic Gradient Descent"
|
||||
|
||||
self.model = model
|
||||
|
|
@ -33,6 +33,13 @@ class opt_SGD(Optimizer):
|
|||
self.center = center
|
||||
self.param_traces = [('noise',[])]
|
||||
self.iteration_file = iteration_file
|
||||
self.learning_rate_adaptation = learning_rate_adaptation
|
||||
if self.learning_rate_adaptation != None:
|
||||
if self.learning_rate_adaptation == 'annealing':
|
||||
self.learning_rate_0 = self.learning_rate
|
||||
else:
|
||||
self.learning_rate_0 = self.learning_rate.mean()
|
||||
|
||||
# if len([p for p in self.model.kern.parts if p.name == 'bias']) == 1:
|
||||
# self.param_traces.append(('bias',[]))
|
||||
# if len([p for p in self.model.kern.parts if p.name == 'linear']) == 1:
|
||||
|
|
@ -204,6 +211,7 @@ class opt_SGD(Optimizer):
|
|||
|
||||
ci = self.shift_constraints(j)
|
||||
f, fp = f_fp(self.x_opt[j])
|
||||
|
||||
step[j] = self.momentum * step[j] + self.learning_rate[j] * fp
|
||||
self.x_opt[j] -= step[j]
|
||||
self.restore_constraints(ci)
|
||||
|
|
@ -216,9 +224,53 @@ class opt_SGD(Optimizer):
|
|||
|
||||
return f, step, self.model.N
|
||||
|
||||
def adapt_learning_rate(self, t):
|
||||
if self.learning_rate_adaptation == 'adagrad':
|
||||
if t > 5:
|
||||
g = np.array(self.grads)
|
||||
l2_g = np.sqrt(np.square(g).sum(0))
|
||||
self.learning_rate = 0.001/l2_g
|
||||
else:
|
||||
self.learning_rate = np.zeros_like(self.learning_rate)
|
||||
elif self.learning_rate_adaptation == 'annealing':
|
||||
self.learning_rate = self.learning_rate_0/(1+float(t+1)/10)
|
||||
elif self.learning_rate_adaptation == 'semi_pesky':
|
||||
if self.model.__class__.__name__ == 'Bayesian_GPLVM':
|
||||
if t == 0:
|
||||
N = self.model.N
|
||||
Q = self.model.Q
|
||||
M = self.model.M
|
||||
|
||||
iip_pos = np.arange(2*N*Q,2*N*Q+M*Q)
|
||||
mu_pos = np.arange(0,N*Q)
|
||||
S_pos = np.arange(N*Q,2*N*Q)
|
||||
self.vbparam_dict = {'iip': [iip_pos],
|
||||
'mu': [mu_pos],
|
||||
'S': [S_pos]}
|
||||
|
||||
for k in self.vbparam_dict.keys():
|
||||
hbar_t = 0.0
|
||||
tau_t = 1000.0
|
||||
gbar_t = 0.0
|
||||
self.vbparam_dict[k].append(hbar_t)
|
||||
self.vbparam_dict[k].append(tau_t)
|
||||
self.vbparam_dict[k].append(gbar_t)
|
||||
|
||||
g_t = self.model.grads
|
||||
|
||||
for k in self.vbparam_dict.keys():
|
||||
pos, hbar_t, tau_t, gbar_t = self.vbparam_dict[k]
|
||||
|
||||
gbar_t = (1-1/tau_t)*gbar_t + 1/tau_t * g_t[pos]
|
||||
hbar_t = (1-1/tau_t)*hbar_t + 1/tau_t * np.dot(g_t[pos].T, g_t[pos])
|
||||
self.learning_rate[pos] = np.dot(gbar_t.T, gbar_t) / hbar_t
|
||||
tau_t = tau_t*(1-self.learning_rate[pos]) + 1
|
||||
self.vbparam_dict[k] = [pos, hbar_t, tau_t, gbar_t]
|
||||
|
||||
|
||||
def opt(self, f_fp=None, f=None, fp=None):
|
||||
self.x_opt = self.model._get_params_transformed()
|
||||
self.model.grads = np.zeros_like(self.x_opt)
|
||||
self.grads = []
|
||||
|
||||
X, Y = self.model.X.copy(), self.model.likelihood.Y.copy()
|
||||
|
||||
|
|
@ -235,6 +287,7 @@ class opt_SGD(Optimizer):
|
|||
|
||||
step = np.zeros_like(num_params)
|
||||
for it in range(self.iterations):
|
||||
self.model.grads = np.zeros_like(self.x_opt) # TODO this is ugly
|
||||
|
||||
if it == 0 or self.self_paced is False:
|
||||
features = np.random.permutation(Y.shape[1])
|
||||
|
|
@ -272,16 +325,17 @@ class opt_SGD(Optimizer):
|
|||
sys.stdout.write(status)
|
||||
sys.stdout.flush()
|
||||
self.param_traces['noise'].append(noise)
|
||||
NLL.append(f)
|
||||
|
||||
self.fopt_trace.append(f)
|
||||
NLL.append(f)
|
||||
self.fopt_trace.append(NLL[-1])
|
||||
# fig = plt.figure('traces')
|
||||
# plt.clf()
|
||||
# plt.plot(self.param_traces['noise'])
|
||||
|
||||
# for k in self.param_traces.keys():
|
||||
# self.param_traces[k].append(self.model.get(k)[0])
|
||||
|
||||
self.grads.append(self.model.grads.tolist())
|
||||
self.adapt_learning_rate(it)
|
||||
# should really be a sum(), but earlier samples in the iteration will have a very crappy ll
|
||||
self.f_opt = np.mean(NLL)
|
||||
self.model.N = N
|
||||
|
|
@ -303,6 +357,6 @@ class opt_SGD(Optimizer):
|
|||
|
||||
if self.messages != 0:
|
||||
sys.stdout.write('\r' + ' '*len(status)*2 + ' \r')
|
||||
status = "SGD Iteration: {0: 3d}/{1: 3d} f: {2: 2.3f}\n".format(it+1, self.iterations, self.f_opt)
|
||||
status = "SGD Iteration: {0: 3d}/{1: 3d} f: {2: 2.3f} max eta: {3: 1.5f}\n".format(it+1, self.iterations, self.f_opt, self.learning_rate.max())
|
||||
sys.stdout.write(status)
|
||||
sys.stdout.flush()
|
||||
|
|
|
|||
|
|
@ -55,7 +55,8 @@ class bias(kernpart):
|
|||
target += self.variance
|
||||
|
||||
def psi1(self, Z, mu, S, target):
|
||||
target += self.variance
|
||||
self._psi1 = self.variance
|
||||
target += self._psi1
|
||||
|
||||
def psi2(self, Z, mu, S, target):
|
||||
target += self.variance**2
|
||||
|
|
|
|||
130
GPy/kern/kern.py
130
GPy/kern/kern.py
|
|
@ -315,31 +315,20 @@ class kern(parameterised):
|
|||
|
||||
# compute the "cross" terms
|
||||
# TODO: input_slices needed
|
||||
crossterms = 0
|
||||
|
||||
for p1, p2 in itertools.combinations(self.parts, 2):
|
||||
# white doesn;t combine with anything
|
||||
if p1.name == 'white' or p2.name == 'white':
|
||||
pass
|
||||
# rbf X bias
|
||||
elif p1.name == 'bias' and p2.name == 'rbf':
|
||||
target += p1.variance * (p2._psi1[:, :, None] + p2._psi1[:, None, :])
|
||||
elif p2.name == 'bias' and p1.name == 'rbf':
|
||||
target += p2.variance * (p1._psi1[:, :, None] + p1._psi1[:, None, :])
|
||||
# linear X bias
|
||||
elif p1.name == 'bias' and p2.name == 'linear':
|
||||
tmp = np.zeros((mu.shape[0], Z.shape[0]))
|
||||
p2.psi1(Z, mu, S, tmp)
|
||||
target += p1.variance * (tmp[:, :, None] + tmp[:, None, :])
|
||||
elif p2.name == 'bias' and p1.name == 'linear':
|
||||
tmp = np.zeros((mu.shape[0], Z.shape[0]))
|
||||
p1.psi1(Z, mu, S, tmp)
|
||||
target += p2.variance * (tmp[:, :, None] + tmp[:, None, :])
|
||||
# rbf X linear
|
||||
elif p1.name == 'linear' and p2.name == 'rbf':
|
||||
raise NotImplementedError # TODO
|
||||
elif p2.name == 'linear' and p1.name == 'rbf':
|
||||
raise NotImplementedError # TODO
|
||||
else:
|
||||
raise NotImplementedError, "psi2 cannot be computed for this kernel"
|
||||
|
||||
# TODO psi1 this must be faster/better/precached/more nice
|
||||
tmp1 = np.zeros((mu.shape[0], Z.shape[0]))
|
||||
p1.psi1(Z, mu, S, tmp1)
|
||||
tmp2 = np.zeros((mu.shape[0], Z.shape[0]))
|
||||
p2.psi1(Z, mu, S, tmp2)
|
||||
|
||||
prod = np.multiply(tmp1, tmp2)
|
||||
crossterms += prod[:,:,None] + prod[:, None, :]
|
||||
|
||||
target += crossterms
|
||||
return target
|
||||
|
||||
def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S):
|
||||
|
|
@ -348,71 +337,34 @@ class kern(parameterised):
|
|||
|
||||
# compute the "cross" terms
|
||||
# TODO: better looping, input_slices
|
||||
for i1, i2 in itertools.combinations(range(len(self.parts)), 2):
|
||||
for i1, i2 in itertools.permutations(range(len(self.parts)), 2):
|
||||
p1, p2 = self.parts[i1], self.parts[i2]
|
||||
# ipsl1, ipsl2 = self.input_slices[i1], self.input_slices[i2]
|
||||
ps1, ps2 = self.param_slices[i1], self.param_slices[i2]
|
||||
|
||||
# white doesn;t combine with anything
|
||||
if p1.name == 'white' or p2.name == 'white':
|
||||
pass
|
||||
# rbf X bias
|
||||
elif p1.name == 'bias' and p2.name == 'rbf':
|
||||
p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1.variance * 2., Z, mu, S, target[ps2])
|
||||
p1.dpsi1_dtheta(dL_dpsi2.sum(1) * p2._psi1 * 2., Z, mu, S, target[ps1])
|
||||
elif p2.name == 'bias' and p1.name == 'rbf':
|
||||
p1.dpsi1_dtheta(dL_dpsi2.sum(1) * p2.variance * 2., Z, mu, S, target[ps1])
|
||||
p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1._psi1 * 2., Z, mu, S, target[ps2])
|
||||
# linear X bias
|
||||
elif p1.name == 'bias' and p2.name == 'linear':
|
||||
p2.dpsi1_dtheta(dL_dpsi2.sum(1) * p1.variance * 2., Z, mu, S, target[ps2]) # [ps1])
|
||||
psi1 = np.zeros((mu.shape[0], Z.shape[0]))
|
||||
p2.psi1(Z, mu, S, psi1)
|
||||
p1.dpsi1_dtheta(dL_dpsi2.sum(1) * psi1 * 2., Z, mu, S, target[ps1])
|
||||
elif p2.name == 'bias' and p1.name == 'linear':
|
||||
p1.dpsi1_dtheta(dL_dpsi2.sum(1) * p2.variance * 2., Z, mu, S, target[ps1])
|
||||
psi1 = np.zeros((mu.shape[0], Z.shape[0]))
|
||||
p1.psi1(Z, mu, S, psi1)
|
||||
p2.dpsi1_dtheta(dL_dpsi2.sum(1) * psi1 * 2., Z, mu, S, target[ps2])
|
||||
# rbf X linear
|
||||
elif p1.name == 'linear' and p2.name == 'rbf':
|
||||
raise NotImplementedError # TODO
|
||||
elif p2.name == 'linear' and p1.name == 'rbf':
|
||||
raise NotImplementedError # TODO
|
||||
else:
|
||||
raise NotImplementedError, "psi2 cannot be computed for this kernel"
|
||||
tmp = np.zeros((mu.shape[0], Z.shape[0]))
|
||||
p1.psi1(Z, mu, S, tmp)
|
||||
p2.dpsi1_dtheta((tmp[:,None,:]*dL_dpsi2).sum(1)*2., Z, mu, S, target[ps2])
|
||||
|
||||
return self._transform_gradients(target)
|
||||
|
||||
def dpsi2_dZ(self, dL_dpsi2, Z, mu, S):
|
||||
target = np.zeros_like(Z)
|
||||
[p.dpsi2_dZ(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
|
||||
#target *= 2
|
||||
|
||||
# compute the "cross" terms
|
||||
# TODO: we need input_slices here.
|
||||
for p1, p2 in itertools.combinations(self.parts, 2):
|
||||
# white doesn;t combine with anything
|
||||
if p1.name == 'white' or p2.name == 'white':
|
||||
pass
|
||||
# rbf X bias
|
||||
elif p1.name == 'bias' and p2.name == 'rbf':
|
||||
p2.dpsi1_dX(dL_dpsi2.sum(1).T * p1.variance, Z, mu, S, target)
|
||||
elif p2.name == 'bias' and p1.name == 'rbf':
|
||||
p1.dpsi1_dZ(dL_dpsi2.sum(1).T * p2.variance, Z, mu, S, target)
|
||||
# linear X bias
|
||||
elif p1.name == 'bias' and p2.name == 'linear':
|
||||
p2.dpsi1_dZ(dL_dpsi2.sum(1).T * p1.variance, Z, mu, S, target)
|
||||
elif p2.name == 'bias' and p1.name == 'linear':
|
||||
p1.dpsi1_dZ(dL_dpsi2.sum(1).T * p2.variance, Z, mu, S, target)
|
||||
# rbf X linear
|
||||
elif p1.name == 'linear' and p2.name == 'rbf':
|
||||
raise NotImplementedError # TODO
|
||||
elif p2.name == 'linear' and p1.name == 'rbf':
|
||||
raise NotImplementedError # TODO
|
||||
else:
|
||||
raise NotImplementedError, "psi2 cannot be computed for this kernel"
|
||||
for p1, p2 in itertools.permutations(self.parts, 2):
|
||||
if p1.name == 'linear' and p2.name == 'linear':
|
||||
raise NotImplementedError("We don't handle linear/linear cross-terms")
|
||||
tmp = np.zeros((mu.shape[0], Z.shape[0]))
|
||||
p1.psi1(Z, mu, S, tmp)
|
||||
tmp2 = np.zeros_like(target)
|
||||
p2.dpsi1_dZ((tmp[:,None,:]*dL_dpsi2).sum(1).T, Z, mu, S, tmp2)
|
||||
target += tmp2
|
||||
|
||||
return target * 2.
|
||||
return target * 2
|
||||
|
||||
def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S):
|
||||
target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1]))
|
||||
|
|
@ -420,27 +372,13 @@ class kern(parameterised):
|
|||
|
||||
# compute the "cross" terms
|
||||
# TODO: we need input_slices here.
|
||||
for p1, p2 in itertools.combinations(self.parts, 2):
|
||||
# white doesn;t combine with anything
|
||||
if p1.name == 'white' or p2.name == 'white':
|
||||
pass
|
||||
# rbf X bias
|
||||
elif p1.name == 'bias' and p2.name == 'rbf':
|
||||
p2.dpsi1_dmuS(dL_dpsi2.sum(1).T * p1.variance * 2., Z, mu, S, target_mu, target_S)
|
||||
elif p2.name == 'bias' and p1.name == 'rbf':
|
||||
p1.dpsi1_dmuS(dL_dpsi2.sum(1).T * p2.variance * 2., Z, mu, S, target_mu, target_S)
|
||||
# linear X bias
|
||||
elif p1.name == 'bias' and p2.name == 'linear':
|
||||
p2.dpsi1_dmuS(dL_dpsi2.sum(1).T * p1.variance * 2., Z, mu, S, target_mu, target_S)
|
||||
elif p2.name == 'bias' and p1.name == 'linear':
|
||||
p1.dpsi1_dmuS(dL_dpsi2.sum(1).T * p2.variance * 2., Z, mu, S, target_mu, target_S)
|
||||
# rbf X linear
|
||||
elif p1.name == 'linear' and p2.name == 'rbf':
|
||||
raise NotImplementedError # TODO
|
||||
elif p2.name == 'linear' and p1.name == 'rbf':
|
||||
raise NotImplementedError # TODO
|
||||
else:
|
||||
raise NotImplementedError, "psi2 cannot be computed for this kernel"
|
||||
for p1, p2 in itertools.permutations(self.parts, 2):
|
||||
if p1.name == 'linear' and p2.name == 'linear':
|
||||
raise NotImplementedError("We don't handle linear/linear cross-terms")
|
||||
|
||||
tmp = np.zeros((mu.shape[0], Z.shape[0]))
|
||||
p1.psi1(Z, mu, S, tmp)
|
||||
p2.dpsi1_dmuS((tmp[:,None,:]*dL_dpsi2).sum(1).T*2., Z, mu, S, target_mu, target_S)
|
||||
|
||||
return target_mu, target_S
|
||||
|
||||
|
|
|
|||
|
|
@ -54,5 +54,3 @@ class kernpart(object):
|
|||
raise NotImplementedError
|
||||
def dK_dX(self,X,X2,target):
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -18,6 +18,7 @@ class white(kernpart):
|
|||
self.Nparam = 1
|
||||
self.name = 'white'
|
||||
self._set_params(np.array([variance]).flatten())
|
||||
self._psi1 = 0 # TODO: more elegance here
|
||||
|
||||
def _get_params(self):
|
||||
return self.variance
|
||||
|
|
@ -81,4 +82,3 @@ class white(kernpart):
|
|||
|
||||
def dpsi2_dmuS(self,dL_dpsi2,Z,mu,S,target_mu,target_S):
|
||||
pass
|
||||
|
||||
|
|
|
|||
|
|
@ -69,6 +69,7 @@ class Gaussian(likelihood):
|
|||
# Note. for D>1, we need to re-normalise all the outputs independently.
|
||||
# This will mess up computations of diag(true_var), below.
|
||||
# note that the upper, lower quantiles should be the same shape as mean
|
||||
# Augment the output variance with the likelihood variance and rescale.
|
||||
true_var = (var + np.eye(var.shape[0]) * self._variance) * self._scale ** 2
|
||||
_5pc = mean - 2.*np.sqrt(np.diag(true_var))
|
||||
_95pc = mean + 2.*np.sqrt(np.diag(true_var))
|
||||
|
|
|
|||
|
|
@ -171,9 +171,6 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
self.dbound_dZtheta = sparse_GP._log_likelihood_gradients(self)
|
||||
return np.hstack((self.dbound_dmuS.flatten(), self.dbound_dZtheta))
|
||||
|
||||
def _log_likelihood_normal_gradients(self):
|
||||
Si, _, _, _ = pdinv(self.X_variance)
|
||||
|
||||
def plot_latent(self, which_indices=None, *args, **kwargs):
|
||||
|
||||
if which_indices is None:
|
||||
|
|
@ -208,23 +205,25 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
else:
|
||||
colors = iter(colors)
|
||||
plots = []
|
||||
x = np.arange(self.X.shape[0])
|
||||
for i in range(self.X.shape[1]):
|
||||
if axes is None:
|
||||
ax = fig.add_subplot(self.X.shape[1], 1, i + 1)
|
||||
else:
|
||||
ax = axes[i]
|
||||
ax.plot(self.X, c='k', alpha=.3)
|
||||
plots.extend(ax.plot(self.X.T[i], c=colors.next(), label=r"$\mathbf{{X_{{{}}}}}$".format(i)))
|
||||
ax.fill_between(np.arange(self.X.shape[0]),
|
||||
plots.extend(ax.plot(x, self.X.T[i], c=colors.next(), label=r"$\mathbf{{X_{{{}}}}}$".format(i)))
|
||||
ax.fill_between(x,
|
||||
self.X.T[i] - 2 * np.sqrt(self.X_variance.T[i]),
|
||||
self.X.T[i] + 2 * np.sqrt(self.X_variance.T[i]),
|
||||
facecolor=plots[-1].get_color(),
|
||||
alpha=.3)
|
||||
ax.legend(borderaxespad=0.)
|
||||
ax.set_xlim(x.min(), x.max())
|
||||
if i < self.X.shape[1] - 1:
|
||||
ax.set_xticklabels('')
|
||||
pylab.draw()
|
||||
fig.tight_layout(h_pad=.01) # , rect=(0, 0, 1, .95))
|
||||
fig.tight_layout(h_pad=.01) # , rect=(0, 0, 1, .95))
|
||||
return fig
|
||||
|
||||
def _debug_filter_params(self, x):
|
||||
|
|
@ -263,7 +262,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
kllls = np.array(self._savedklll)
|
||||
LL, = ax1.plot(kllls[:, 0], kllls[:, 1] - kllls[:, 2], '-', label=r'$\log p(\mathbf{Y})$', mew=1.5)
|
||||
KL, = ax1.plot(kllls[:, 0], kllls[:, 2], '-', label=r'$\mathcal{KL}(p||q)$', mew=1.5)
|
||||
L, = ax1.plot(kllls[:, 0], kllls[:, 1], '-', label=r'$L$', mew=1.5) # \mathds{E}_{q(\mathbf{X})}[p(\mathbf{Y|X})\frac{p(\mathbf{X})}{q(\mathbf{X})}]
|
||||
L, = ax1.plot(kllls[:, 0], kllls[:, 1], '-', label=r'$L$', mew=1.5) # \mathds{E}_{q(\mathbf{X})}[p(\mathbf{Y|X})\frac{p(\mathbf{X})}{q(\mathbf{X})}]
|
||||
|
||||
param_dict = dict(self._savedparams)
|
||||
gradient_dict = dict(self._savedgradients)
|
||||
|
|
@ -411,7 +410,7 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
|
||||
# parameter changes
|
||||
# ax2 = pylab.subplot2grid((4, 1), (1, 0), 3, 1, projection='3d')
|
||||
button_options = [0, 0] # [0]: clicked -- [1]: dragged
|
||||
button_options = [0, 0] # [0]: clicked -- [1]: dragged
|
||||
|
||||
def update_plots(event):
|
||||
if button_options[0] and not button_options[1]:
|
||||
|
|
@ -483,4 +482,4 @@ class Bayesian_GPLVM(sparse_GP, GPLVM):
|
|||
cidp = figs[0].canvas.mpl_connect('button_press_event', onclick)
|
||||
cidd = figs[0].canvas.mpl_connect('motion_notify_event', motion)
|
||||
|
||||
return ax1, ax2, ax3, ax4, ax5 # , ax6, ax7
|
||||
return ax1, ax2, ax3, ax4, ax5 # , ax6, ax7
|
||||
|
|
|
|||
|
|
@ -3,7 +3,7 @@
|
|||
|
||||
import numpy as np
|
||||
import pylab as pb
|
||||
from ..util.linalg import mdot, jitchol, tdot, symmetrify, backsub_both_sides
|
||||
from ..util.linalg import mdot, jitchol, tdot, symmetrify, backsub_both_sides,chol_inv
|
||||
from ..util.plot import gpplot
|
||||
from .. import kern
|
||||
from GP import GP
|
||||
|
|
@ -16,9 +16,9 @@ class sparse_GP(GP):
|
|||
:param X: inputs
|
||||
:type X: np.ndarray (N x Q)
|
||||
:param likelihood: a likelihood instance, containing the observed data
|
||||
:type likelihood: GPy.likelihood.(Gaussian | EP)
|
||||
:param kernel : the kernel/covariance function. See link kernels
|
||||
:type kernel: a GPy kernel
|
||||
:type likelihood: GPy.likelihood.(Gaussian | EP | Laplace)
|
||||
:param kernel : the kernel (covariance function). See link kernels
|
||||
:type kernel: a GPy.kern.kern instance
|
||||
:param X_variance: The uncertainty in the measurements of X (Gaussian variance)
|
||||
:type X_variance: np.ndarray (N x Q) | None
|
||||
:param Z: inducing inputs (optional, see note)
|
||||
|
|
@ -30,8 +30,6 @@ class sparse_GP(GP):
|
|||
"""
|
||||
|
||||
def __init__(self, X, likelihood, kernel, Z, X_variance=None, normalize_X=False):
|
||||
# self.scale_factor = 100.0 # a scaling factor to help keep the algorithm stable
|
||||
# self.auto_scale_factor = False
|
||||
self.Z = Z
|
||||
self.M = Z.shape[0]
|
||||
self.likelihood = likelihood
|
||||
|
|
@ -63,49 +61,29 @@ class sparse_GP(GP):
|
|||
self.psi2 = None
|
||||
|
||||
def _computations(self):
|
||||
# sf = self.scale_factor
|
||||
# sf2 = sf ** 2
|
||||
|
||||
# factor Kmm
|
||||
self.Lm = jitchol(self.Kmm)
|
||||
|
||||
# The rather complex computations of self.A
|
||||
if self.likelihood.is_heteroscedastic:
|
||||
assert self.likelihood.D == 1 # TODO: what if the likelihood is heterscedatic and there are multiple independent outputs?
|
||||
if self.has_uncertain_inputs:
|
||||
# psi2_beta_scaled = (self.psi2 * (self.likelihood.precision.flatten().reshape(self.N, 1, 1) / sf2)).sum(0)
|
||||
psi2_beta_scaled = (self.psi2 * (self.likelihood.precision.flatten().reshape(self.N, 1, 1))).sum(0)
|
||||
evals, evecs = linalg.eigh(psi2_beta_scaled)
|
||||
clipped_evals = np.clip(evals, 0., 1e6) # TODO: make clipping configurable
|
||||
if not np.allclose(evals, clipped_evals):
|
||||
print "Warning: clipping posterior eigenvalues"
|
||||
tmp = evecs * np.sqrt(clipped_evals)
|
||||
tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
|
||||
self.A = tdot(tmp)
|
||||
if self.has_uncertain_inputs:
|
||||
if self.likelihood.is_heteroscedastic:
|
||||
psi2_beta = (self.psi2 * (self.likelihood.precision.flatten().reshape(self.N, 1, 1))).sum(0)
|
||||
else:
|
||||
# tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.N)) / sf)
|
||||
tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.N)))
|
||||
tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
|
||||
self.A = tdot(tmp)
|
||||
psi2_beta = self.psi2.sum(0) * self.likelihood.precision
|
||||
evals, evecs = linalg.eigh(psi2_beta)
|
||||
clipped_evals = np.clip(evals, 0., 1e6) # TODO: make clipping configurable
|
||||
tmp = evecs * np.sqrt(clipped_evals)
|
||||
else:
|
||||
if self.has_uncertain_inputs:
|
||||
# psi2_beta_scaled = (self.psi2 * (self.likelihood.precision / sf2)).sum(0)
|
||||
psi2_beta_scaled = (self.psi2 * (self.likelihood.precision)).sum(0)
|
||||
evals, evecs = linalg.eigh(psi2_beta_scaled)
|
||||
clipped_evals = np.clip(evals, 0., 1e15) # TODO: make clipping configurable
|
||||
if not np.allclose(evals, clipped_evals):
|
||||
print "Warning: clipping posterior eigenvalues"
|
||||
tmp = evecs * np.sqrt(clipped_evals)
|
||||
tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
|
||||
self.A = tdot(tmp)
|
||||
if self.likelihood.is_heteroscedastic:
|
||||
tmp = self.psi1 * (np.sqrt(self.likelihood.precision.flatten().reshape(1, self.N)))
|
||||
else:
|
||||
# tmp = self.psi1 * (np.sqrt(self.likelihood.precision) / sf)
|
||||
tmp = self.psi1 * (np.sqrt(self.likelihood.precision))
|
||||
tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
|
||||
self.A = tdot(tmp)
|
||||
tmp, _ = linalg.lapack.flapack.dtrtrs(self.Lm, np.asfortranarray(tmp), lower=1)
|
||||
self.A = tdot(tmp)
|
||||
|
||||
|
||||
# factor B
|
||||
# self.B = np.eye(self.M) / sf2 + self.A
|
||||
self.B = np.eye(self.M) + self.A
|
||||
self.LB = jitchol(self.B)
|
||||
|
||||
|
|
@ -121,8 +99,6 @@ class sparse_GP(GP):
|
|||
# Compute dL_dKmm
|
||||
tmp = tdot(self._LBi_Lmi_psi1V)
|
||||
self.DBi_plus_BiPBi = backsub_both_sides(self.LB, self.D * np.eye(self.M) + tmp)
|
||||
# tmp = -0.5 * self.DBi_plus_BiPBi / sf2
|
||||
# tmp += -0.5 * self.B * sf2 * self.D
|
||||
tmp = -0.5 * self.DBi_plus_BiPBi
|
||||
tmp += -0.5 * self.B * self.D
|
||||
tmp += self.D * np.eye(self.M)
|
||||
|
|
@ -132,9 +108,10 @@ class sparse_GP(GP):
|
|||
self.dL_dpsi0 = -0.5 * self.D * (self.likelihood.precision * np.ones([self.N, 1])).flatten()
|
||||
self.dL_dpsi1 = np.dot(self.Cpsi1V, self.likelihood.V.T)
|
||||
dL_dpsi2_beta = 0.5 * backsub_both_sides(self.Lm, self.D * np.eye(self.M) - self.DBi_plus_BiPBi)
|
||||
|
||||
if self.likelihood.is_heteroscedastic:
|
||||
if self.has_uncertain_inputs:
|
||||
self.dL_dpsi2 = self.likelihood.precision[:, None, None] * dL_dpsi2_beta[None, :, :]
|
||||
self.dL_dpsi2 = self.likelihood.precision.flatten()[:, None, None] * dL_dpsi2_beta[None, :, :]
|
||||
else:
|
||||
self.dL_dpsi1 += 2.*np.dot(dL_dpsi2_beta, self.psi1 * self.likelihood.precision.reshape(1, self.N))
|
||||
self.dL_dpsi2 = None
|
||||
|
|
@ -158,7 +135,6 @@ class sparse_GP(GP):
|
|||
else:
|
||||
# likelihood is not heterscedatic
|
||||
self.partial_for_likelihood = -0.5 * self.N * self.D * self.likelihood.precision + 0.5 * self.likelihood.trYYT * self.likelihood.precision ** 2
|
||||
# self.partial_for_likelihood += 0.5 * self.D * (self.psi0.sum() * self.likelihood.precision ** 2 - np.trace(self.A) * self.likelihood.precision * sf2)
|
||||
self.partial_for_likelihood += 0.5 * self.D * (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)))
|
||||
|
||||
|
|
@ -166,16 +142,12 @@ class sparse_GP(GP):
|
|||
|
||||
def log_likelihood(self):
|
||||
""" Compute the (lower bound on the) log marginal likelihood """
|
||||
# sf2 = self.scale_factor ** 2
|
||||
if self.likelihood.is_heteroscedastic:
|
||||
A = -0.5 * self.N * self.D * 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.D * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A) * sf2)
|
||||
B = -0.5 * self.D * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A))
|
||||
else:
|
||||
A = -0.5 * self.N * self.D * (np.log(2.*np.pi) - np.log(self.likelihood.precision)) - 0.5 * self.likelihood.precision * self.likelihood.trYYT
|
||||
# B = -0.5 * self.D * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A) * sf2)
|
||||
B = -0.5 * self.D * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A))
|
||||
# C = -self.D * (np.sum(np.log(np.diag(self.LB))) + 0.5 * self.M * np.log(sf2))
|
||||
C = -self.D * (np.sum(np.log(np.diag(self.LB)))) # + 0.5 * self.M * np.log(sf2))
|
||||
D = 0.5 * np.sum(np.square(self._LBi_Lmi_psi1V))
|
||||
return A + B + C + D
|
||||
|
|
@ -185,14 +157,6 @@ class sparse_GP(GP):
|
|||
self.kern._set_params(p[self.Z.size:self.Z.size + self.kern.Nparam])
|
||||
self.likelihood._set_params(p[self.Z.size + self.kern.Nparam:])
|
||||
self._compute_kernel_matrices()
|
||||
# if self.auto_scale_factor:
|
||||
# self.scale_factor = np.sqrt(self.psi2.sum(0).mean()*self.likelihood.precision)
|
||||
# if self.auto_scale_factor:
|
||||
# if self.likelihood.is_heteroscedastic:
|
||||
# self.scale_factor = max(100,np.sqrt(self.psi2_beta_scaled.sum(0).mean()))
|
||||
# else:
|
||||
# self.scale_factor = np.sqrt(self.psi2.sum(0).mean()*self.likelihood.precision)
|
||||
# self.scale_factor = 100.
|
||||
self._computations()
|
||||
|
||||
def _get_params(self):
|
||||
|
|
@ -205,11 +169,17 @@ class sparse_GP(GP):
|
|||
"""
|
||||
Approximates a non-gaussian likelihood using Expectation Propagation
|
||||
|
||||
For a Gaussian (or direct: TODO) likelihood, no iteration is required:
|
||||
For a Gaussian likelihood, no iteration is required:
|
||||
this function does nothing
|
||||
"""
|
||||
if self.has_uncertain_inputs:
|
||||
raise NotImplementedError, "EP approximation not implemented for uncertain inputs"
|
||||
|
||||
Lmi = chol_inv(self.Lm)
|
||||
Kmmi = tdot(Lmi.T)
|
||||
diag_tr_psi2Kmmi = np.array([np.trace(psi2_Kmmi) for psi2_Kmmi in np.dot(self.psi2,Kmmi)])
|
||||
|
||||
self.likelihood.fit_FITC(self.Kmm,self.psi1,diag_tr_psi2Kmmi) #This uses the fit_FITC code, but does not perfomr a FITC-EP.#TODO solve potential confusion
|
||||
#raise NotImplementedError, "EP approximation not implemented for uncertain inputs"
|
||||
else:
|
||||
self.likelihood.fit_DTC(self.Kmm, self.psi1)
|
||||
# self.likelihood.fit_FITC(self.Kmm,self.psi1,self.psi0)
|
||||
|
|
|
|||
|
|
@ -5,6 +5,8 @@ import GPy
|
|||
import scipy.sparse
|
||||
import scipy.io
|
||||
import cPickle as pickle
|
||||
import urllib2 as url
|
||||
|
||||
data_path = os.path.join(os.path.dirname(__file__), 'datasets')
|
||||
default_seed = 10000
|
||||
|
||||
|
|
@ -15,6 +17,18 @@ def sample_class(f):
|
|||
c = np.where(c, 1, -1)
|
||||
return c
|
||||
|
||||
def fetch_dataset(resource, file_name, messages = True):
|
||||
if messages:
|
||||
print "Downloading resource: " , resource, " ... "
|
||||
response = url.urlopen(resource)
|
||||
# TODO: Some error checking...
|
||||
html = response.read()
|
||||
response.close()
|
||||
with open(file_name, "w") as text_file:
|
||||
text_file.write("%s"%html)
|
||||
if messages:
|
||||
print "Done!"
|
||||
|
||||
def della_gatta_TRP63_gene_expression(gene_number=None):
|
||||
mat_data = scipy.io.loadmat(os.path.join(data_path, 'DellaGattadata.mat'))
|
||||
X = np.double(mat_data['timepoints'])
|
||||
|
|
|
|||
|
|
@ -236,7 +236,7 @@ def tdot(*args, **kwargs):
|
|||
else:
|
||||
return tdot_numpy(*args,**kwargs)
|
||||
|
||||
def DSYR(A,x,alpha=1.):
|
||||
def DSYR_blas(A,x,alpha=1.):
|
||||
"""
|
||||
Performs a symmetric rank-1 update operation:
|
||||
A <- A + alpha * np.dot(x,x.T)
|
||||
|
|
@ -258,6 +258,26 @@ def DSYR(A,x,alpha=1.):
|
|||
x_, byref(INCX), A_, byref(LDA))
|
||||
symmetrify(A,upper=True)
|
||||
|
||||
def DSYR_numpy(A,x,alpha=1.):
|
||||
"""
|
||||
Performs a symmetric rank-1 update operation:
|
||||
A <- A + alpha * np.dot(x,x.T)
|
||||
|
||||
Arguments
|
||||
---------
|
||||
:param A: Symmetric NxN np.array
|
||||
:param x: Nx1 np.array
|
||||
:param alpha: scalar
|
||||
"""
|
||||
A += alpha*np.dot(x[:,None],x[None,:])
|
||||
|
||||
|
||||
def DSYR(*args, **kwargs):
|
||||
if _blas_available:
|
||||
return DSYR_blas(*args,**kwargs)
|
||||
else:
|
||||
return DSYR_numpy(*args,**kwargs)
|
||||
|
||||
def symmetrify(A,upper=False):
|
||||
"""
|
||||
Take the square matrix A and make it symmetrical by copting elements from the lower half to the upper
|
||||
|
|
|
|||
13
GPy/util/mocap_fetch.py
Normal file
13
GPy/util/mocap_fetch.py
Normal file
|
|
@ -0,0 +1,13 @@
|
|||
import GPy
|
||||
import urllib2
|
||||
|
||||
# TODO...
|
||||
class mocap_fetch(base_url = 'http://mocap.cs.cmu.edu:8080/subjects/', skel_store_dir = './', motion_store_dir = './'):
|
||||
def __init__(self):
|
||||
self.base_url = base_url
|
||||
self.store_dir = store_dir
|
||||
self.motion_dict = []
|
||||
|
||||
def fetch_motions(self, motion_dict = None):
|
||||
response = urllib2.urlopen(...)
|
||||
html = response.read()
|
||||
|
|
@ -44,7 +44,7 @@ class vector_show(data_show):
|
|||
|
||||
|
||||
class lvm(data_show):
|
||||
def __init__(self, vals, model, data_visualize, latent_axes=None, latent_index=[0,1]):
|
||||
def __init__(self, vals, model, data_visualize, latent_axes=None, sense_axes=None, latent_index=[0,1]):
|
||||
"""Visualize a latent variable model
|
||||
|
||||
:param model: the latent variable model to visualize.
|
||||
|
|
@ -71,7 +71,7 @@ class lvm(data_show):
|
|||
self.data_visualize = data_visualize
|
||||
self.model = model
|
||||
self.latent_axes = latent_axes
|
||||
|
||||
self.sense_axes = sense_axes
|
||||
self.called = False
|
||||
self.move_on = False
|
||||
self.latent_index = latent_index
|
||||
|
|
@ -81,10 +81,12 @@ class lvm(data_show):
|
|||
self.latent_values = vals
|
||||
self.latent_handle = self.latent_axes.plot([0],[0],'rx',mew=2)[0]
|
||||
self.modify(vals)
|
||||
self.show_sensitivities()
|
||||
|
||||
def modify(self, vals):
|
||||
"""When latent values are modified update the latent representation and ulso update the output visualization."""
|
||||
y = self.model.predict(vals)[0]
|
||||
print y
|
||||
self.data_visualize.modify(y)
|
||||
self.latent_handle.set_data(vals[self.latent_index[0]], vals[self.latent_index[1]])
|
||||
self.axes.figure.canvas.draw()
|
||||
|
|
@ -99,6 +101,7 @@ class lvm(data_show):
|
|||
if event.inaxes!=self.latent_axes: return
|
||||
self.move_on = not self.move_on
|
||||
self.called = True
|
||||
|
||||
def on_move(self, event):
|
||||
if event.inaxes!=self.latent_axes: return
|
||||
if self.called and self.move_on:
|
||||
|
|
@ -107,22 +110,54 @@ class lvm(data_show):
|
|||
self.latent_values[self.latent_index[1]]=event.ydata
|
||||
self.modify(self.latent_values)
|
||||
|
||||
def show_sensitivities(self):
|
||||
# A click in the bar chart axis for selection a dimension.
|
||||
if self.sense_axes != None:
|
||||
self.sense_axes.cla()
|
||||
self.sense_axes.bar(np.arange(self.model.Q),1./self.model.input_sensitivity(),color='b')
|
||||
|
||||
if self.latent_index[1] == self.latent_index[0]:
|
||||
self.sense_axes.bar(np.array(self.latent_index[0]),1./self.model.input_sensitivity()[self.latent_index[0]],color='y')
|
||||
self.sense_axes.bar(np.array(self.latent_index[1]),1./self.model.input_sensitivity()[self.latent_index[1]],color='y')
|
||||
|
||||
else:
|
||||
self.sense_axes.bar(np.array(self.latent_index[0]),1./self.model.input_sensitivity()[self.latent_index[0]],color='g')
|
||||
self.sense_axes.bar(np.array(self.latent_index[1]),1./self.model.input_sensitivity()[self.latent_index[1]],color='r')
|
||||
|
||||
self.sense_axes.figure.canvas.draw()
|
||||
|
||||
|
||||
class lvm_subplots(lvm):
|
||||
"""
|
||||
latent_axes is a np array of dimension np.ceil(Q/2) + 1,
|
||||
one for each pair of the axes, and the last one for the sensitiity bar chart
|
||||
latent_axes is a np array of dimension np.ceil(Q/2),
|
||||
one for each pair of the latent dimensions.
|
||||
"""
|
||||
def __init__(self, vals, model, data_visualize, latent_axes=None, latent_index=[0,1]):
|
||||
lvm.__init__(self, vals, model,data_visualize,latent_axes,[0,1])
|
||||
def __init__(self, vals, model, data_visualize, latent_axes=None, sense_axes=None):
|
||||
self.nplots = int(np.ceil(model.Q/2.))+1
|
||||
lvm.__init__(self,model,data_visualize,latent_axes,latent_index)
|
||||
self.latent_values = np.zeros(2*np.ceil(self.model.Q/2.)) # possibly an extra dimension on this
|
||||
assert latent_axes.size == self.nplots
|
||||
assert len(latent_axes)==self.nplots
|
||||
if vals==None:
|
||||
vals = model.X[0, :]
|
||||
self.latent_values = vals
|
||||
|
||||
for i, axis in enumerate(latent_axes):
|
||||
if i == self.nplots-1:
|
||||
if self.nplots*2!=model.Q:
|
||||
latent_index = [i*2, i*2]
|
||||
lvm.__init__(self, self.latent_vals, model, data_visualize, axis, sense_axes, latent_index=latent_index)
|
||||
else:
|
||||
latent_index = [i*2, i*2+1]
|
||||
lvm.__init__(self, self.latent_vals, model, data_visualize, axis, latent_index=latent_index)
|
||||
|
||||
|
||||
|
||||
class lvm_dimselect(lvm):
|
||||
"""
|
||||
A visualizer for latent variable models which allows selection of the latent dimensions to use by clicking on a bar chart of their length scales.
|
||||
|
||||
For an example of the visualizer's use try:
|
||||
|
||||
GPy.examples.dimensionality_reduction.BGPVLM_oil()
|
||||
|
||||
"""
|
||||
def __init__(self, vals, model, data_visualize, latent_axes=None, sense_axes=None, latent_index=[0, 1]):
|
||||
if latent_axes==None and sense_axes==None:
|
||||
|
|
@ -133,24 +168,9 @@ class lvm_dimselect(lvm):
|
|||
else:
|
||||
self.sense_axes = sense_axes
|
||||
|
||||
lvm.__init__(self,vals,model,data_visualize,latent_axes,latent_index)
|
||||
self.show_sensitivities()
|
||||
lvm.__init__(self,vals,model,data_visualize,latent_axes,sense_axes,latent_index)
|
||||
print "use left and right mouse butons to select dimensions"
|
||||
|
||||
def show_sensitivities(self):
|
||||
# A click in the bar chart axis for selection a dimension.
|
||||
self.sense_axes.cla()
|
||||
self.sense_axes.bar(np.arange(self.model.Q),1./self.model.input_sensitivity(),color='b')
|
||||
|
||||
if self.latent_index[1] == self.latent_index[0]:
|
||||
self.sense_axes.bar(np.array(self.latent_index[0]),1./self.model.input_sensitivity()[self.latent_index[0]],color='y')
|
||||
self.sense_axes.bar(np.array(self.latent_index[1]),1./self.model.input_sensitivity()[self.latent_index[1]],color='y')
|
||||
|
||||
else:
|
||||
self.sense_axes.bar(np.array(self.latent_index[0]),1./self.model.input_sensitivity()[self.latent_index[0]],color='g')
|
||||
self.sense_axes.bar(np.array(self.latent_index[1]),1./self.model.input_sensitivity()[self.latent_index[1]],color='r')
|
||||
|
||||
self.sense_axes.figure.canvas.draw()
|
||||
|
||||
def on_click(self, event):
|
||||
|
||||
|
|
@ -177,12 +197,6 @@ class lvm_dimselect(lvm):
|
|||
self.called = True
|
||||
|
||||
|
||||
def on_move(self, event):
|
||||
if event.inaxes!=self.latent_axes: return
|
||||
if self.called and self.move_on:
|
||||
self.latent_values[self.latent_index[0]]=event.xdata
|
||||
self.latent_values[self.latent_index[1]]=event.ydata
|
||||
self.modify(self.latent_values)
|
||||
|
||||
def on_leave(self,event):
|
||||
latent_values = self.latent_values.copy()
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue