mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-24 14:15:14 +02:00
Merge branch 'mrd' into devel
This commit is contained in:
Max Zwiessele
commit
3bb42481be
3 changed files with 34 additions and 168 deletions
|
|
@ -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,94 +258,22 @@ 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
|
||||
slist, Slist, Ylist = _simulate_sincos(D1, D2, D3, N, M, Q, plot_sim)
|
||||
|
||||
|
|
@ -386,50 +283,16 @@ 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)
|
||||
|
||||
print "initializing beta"
|
||||
cstr = "noise"
|
||||
m.unconstrain(cstr); m.constrain_fixed(cstr)
|
||||
m.optimize('scg', messages=1, max_f_eval=2e3, gtol=100)
|
||||
|
||||
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
|
||||
|
|
|
|||
|
|
@ -154,5 +154,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
|
||||
|
|
|
|||
|
|
@ -208,23 +208,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 +265,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 +413,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 +485,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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue