Source code for GPy.examples.dimensionality_reduction

# Copyright (c) 2012-2014, GPy authors (see AUTHORS.txt).
# Licensed under the BSD 3-clause license (see LICENSE.txt)
import numpy as _np

default_seed = 123344

# default_seed = _np.random.seed(123344)


[docs]def bgplvm_test_model(optimize=False, verbose=1, plot=False, output_dim=200, nan=False): """ model for testing purposes. Samples from a GP with rbf kernel and learns the samples with a new kernel. Normally not for optimization, just model cheking """ import GPy num_inputs = 13 num_inducing = 5 if plot: output_dim = 1 input_dim = 3 else: input_dim = 2 output_dim = output_dim # generate GPLVM-like data X = _np.random.rand(num_inputs, input_dim) lengthscales = _np.random.rand(input_dim) k = GPy.kern.RBF(input_dim, 0.5, lengthscales, ARD=True) K = k.K(X) Y = _np.random.multivariate_normal(_np.zeros(num_inputs), K, (output_dim,)).T # k = GPy.kern.RBF_inv(input_dim, .5, _np.ones(input_dim) * 2., ARD=True) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim) # k = GPy.kern.linear(input_dim)# + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001) # k = GPy.kern.RBF(input_dim, ARD = False) + GPy.kern.white(input_dim, 0.00001) # k = GPy.kern.RBF(input_dim, .5, _np.ones(input_dim) * 2., ARD=True) + GPy.kern.RBF(input_dim, .3, _np.ones(input_dim) * .2, ARD=True) # k = GPy.kern.RBF(input_dim, .5, 2., ARD=0) + GPy.kern.RBF(input_dim, .3, .2, ARD=0) # k = GPy.kern.RBF(input_dim, .5, _np.ones(input_dim) * 2., ARD=True) + GPy.kern.linear(input_dim, _np.ones(input_dim) * .2, ARD=True) p = 0.3 m = GPy.models.BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing) if nan: m.inference_method = ( GPy.inference.latent_function_inference.var_dtc.VarDTCMissingData() ) m.Y[_np.random.binomial(1, p, size=(Y.shape)).astype(bool)] = _np.nan m.parameters_changed() # =========================================================================== # randomly obstruct data with percentage p # =========================================================================== # m2 = GPy.models.BayesianGPLVMWithMissingData(Y_obstruct, input_dim, kernel=k, num_inducing=num_inducing) # m.lengthscales = lengthscales if plot: import matplotlib.pyplot as pb m.plot() pb.title("PCA initialisation") # m2.plot() # pb.title('PCA initialisation') if optimize: m.optimize("scg", messages=verbose) # m2.optimize('scg', messages=verbose) if plot: m.plot() pb.title("After optimisation") # m2.plot() # pb.title('After optimisation') return m
[docs]def gplvm_oil_100(optimize=True, verbose=1, plot=True): import GPy import pods data = pods.datasets.oil_100() Y = data["X"] # create simple GP model kernel = GPy.kern.RBF(6, ARD=True) + GPy.kern.Bias(6) m = GPy.models.GPLVM(Y, 6, kernel=kernel) m.data_labels = data["Y"].argmax(axis=1) if optimize: m.optimize("scg", messages=verbose) if plot: m.plot_latent(labels=m.data_labels) return m
[docs]def sparse_gplvm_oil( optimize=True, verbose=0, plot=True, N=100, Q=6, num_inducing=15, max_iters=50 ): import GPy import pods _np.random.seed(0) data = pods.datasets.oil() Y = data["X"][:N] Y = Y - Y.mean(0) Y /= Y.std(0) # Create the model kernel = GPy.kern.RBF(Q, ARD=True) + GPy.kern.Bias(Q) m = GPy.models.SparseGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing) m.data_labels = data["Y"][:N].argmax(axis=1) if optimize: m.optimize("scg", messages=verbose, max_iters=max_iters) if plot: m.plot_latent(labels=m.data_labels) m.kern.plot_ARD() return m
[docs]def swiss_roll( optimize=True, verbose=1, plot=True, N=1000, num_inducing=25, Q=4, sigma=0.2 ): import GPy from pods.datasets import swiss_roll_generated from GPy.models import BayesianGPLVM data = swiss_roll_generated(num_samples=N, sigma=sigma) Y = data["Y"] Y -= Y.mean() Y /= Y.std() t = data["t"] c = data["colors"] try: from sklearn.manifold.isomap import Isomap iso = Isomap().fit(Y) X = iso.embedding_ if Q > 2: X = _np.hstack((X, _np.random.randn(N, Q - 2))) except ImportError: X = _np.random.randn(N, Q) if plot: import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # @UnusedImport fig = plt.figure("Swiss Roll Data") ax = fig.add_subplot(121, projection="3d") ax.scatter(*Y.T, c=c) ax.set_title("Swiss Roll") ax = fig.add_subplot(122) ax.scatter(*X.T[:2], c=c) ax.set_title("BGPLVM init") var = 0.5 S = ( var * _np.ones_like(X) + _np.clip(_np.random.randn(N, Q) * var ** 2, -(1 - var), (1 - var)) ) + 0.001 Z = _np.random.permutation(X)[:num_inducing] kernel = ( GPy.kern.RBF(Q, ARD=True) + GPy.kern.Bias(Q, _np.exp(-2)) + GPy.kern.White(Q, _np.exp(-2)) ) m = BayesianGPLVM( Y, Q, X=X, X_variance=S, num_inducing=num_inducing, Z=Z, kernel=kernel ) m.data_colors = c m.data_t = t if optimize: m.optimize("bfgs", messages=verbose, max_iters=2e3) if plot: fig = plt.figure("fitted") ax = fig.add_subplot(111) s = m.input_sensitivity().argsort()[::-1][:2] ax.scatter(*m.X.mean.T[s], c=c) return m
[docs]def bgplvm_oil( optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40, max_iters=1000, **k ): import GPy from matplotlib import pyplot as plt import numpy as np _np.random.seed(0) try: import pods data = pods.datasets.oil() except ImportError: data = GPy.util.datasets.oil() kernel = GPy.kern.RBF( Q, 1.0, 1.0 / _np.random.uniform(0, 1, (Q,)), ARD=True ) # + GPy.kern.Bias(Q, _np.exp(-2)) Y = data["X"][:N] m = GPy.models.BayesianGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing, **k) m.data_labels = data["Y"][:N].argmax(axis=1) if optimize: m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: fig, (latent_axes, sense_axes) = plt.subplots(1, 2) m.plot_latent(ax=latent_axes, labels=m.data_labels) data_show = GPy.plotting.matplot_dep.visualize.vector_show((m.Y[0, :])) lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect( m.X.mean.values[0:1, :], # @UnusedVariable m, data_show, latent_axes=latent_axes, sense_axes=sense_axes, labels=m.data_labels, ) input("Press enter to finish") plt.close(fig) return m
[docs]def ssgplvm_oil( optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40, max_iters=1000, **k ): import GPy from matplotlib import pyplot as plt import pods _np.random.seed(0) data = pods.datasets.oil() kernel = GPy.kern.RBF( Q, 1.0, 1.0 / _np.random.uniform(0, 1, (Q,)), ARD=True ) # + GPy.kern.Bias(Q, _np.exp(-2)) Y = data["X"][:N] m = GPy.models.SSGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing, **k) m.data_labels = data["Y"][:N].argmax(axis=1) if optimize: m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: fig, (latent_axes, sense_axes) = plt.subplots(1, 2) m.plot_latent(ax=latent_axes, labels=m.data_labels) data_show = GPy.plotting.matplot_dep.visualize.vector_show((m.Y[0, :])) lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect( m.X.mean.values[0:1, :], # @UnusedVariable m, data_show, latent_axes=latent_axes, sense_axes=sense_axes, labels=m.data_labels, ) input("Press enter to finish") plt.close(fig) return m
def _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim=False): """Simulate some data drawn from a matern covariance and a periodic exponential for use in MRD demos.""" Q_signal = 4 import GPy import numpy as np np.random.seed(3000) k = GPy.kern.Matern32( Q_signal, 1.0, lengthscale=(np.random.uniform(1, 6, Q_signal)), ARD=1 ) for i in range(Q_signal): k += GPy.kern.PeriodicExponential( 1, variance=1.0, active_dims=[i], period=3.0, lower=-2, upper=6 ) t = np.c_[[np.linspace(-1, 5, N) for _ in range(Q_signal)]].T K = k.K(t) s2, s1, s3, sS = np.random.multivariate_normal(np.zeros(K.shape[0]), K, size=(4))[ :, :, None ] Y1, Y2, Y3, S1, S2, S3 = _generate_high_dimensional_output( D1, D2, D3, s1, s2, s3, sS ) slist = [sS, s1, s2, s3] slist_names = ["sS", "s1", "s2", "s3"] Ylist = [Y1, Y2, Y3] if plot_sim: from matplotlib import pyplot as plt import matplotlib.cm as cm import itertools fig = plt.figure("MRD Simulation Data", figsize=(8, 6)) fig.clf() ax = fig.add_subplot(2, 1, 1) labls = slist_names for S, lab in zip(slist, labls): ax.plot(S, label=lab) ax.legend() for i, Y in enumerate(Ylist): ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i) ax.imshow(Y, aspect="auto", cmap=cm.gray) # @UndefinedVariable ax.set_title("Y{}".format(i + 1)) plt.draw() plt.tight_layout() return slist, [S1, S2, S3], Ylist def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False): """Simulate some data drawn from sine and cosine for use in demos of MRD""" _np.random.seed(1234) x = _np.linspace(0, 4 * _np.pi, N)[:, None] s1 = _np.vectorize(lambda x: _np.sin(x)) s2 = _np.vectorize(lambda x: _np.cos(x)) s3 = _np.vectorize(lambda x: -_np.exp(-_np.cos(2 * x))) sS = _np.vectorize(lambda x: _np.cos(x)) s1 = s1(x) s2 = s2(x) s3 = s3(x) sS = sS(x) s1 -= s1.mean() s1 /= s1.std(0) s2 -= s2.mean() s2 /= s2.std(0) s3 -= s3.mean() s3 /= s3.std(0) sS -= sS.mean() sS /= sS.std(0) Y1, Y2, Y3, S1, S2, S3 = _generate_high_dimensional_output( D1, D2, D3, s1, s2, s3, sS ) slist = [sS, s1, s2, s3] slist_names = ["sS", "s1", "s2", "s3"] Ylist = [Y1, Y2, Y3] if plot_sim: from matplotlib import pyplot as plt import matplotlib.cm as cm import itertools fig = plt.figure("MRD Simulation Data", figsize=(8, 6)) fig.clf() ax = fig.add_subplot(2, 1, 1) labls = slist_names for S, lab in zip(slist, labls): ax.plot(S, label=lab) ax.legend() for i, Y in enumerate(Ylist): ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i) ax.imshow(Y, aspect="auto", cmap=cm.gray) # @UndefinedVariable ax.set_title("Y{}".format(i + 1)) plt.draw() plt.tight_layout() return slist, [S1, S2, S3], Ylist def _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS): S1 = _np.hstack([s1, sS]) S2 = _np.hstack([sS]) S3 = _np.hstack([s1, s3, sS]) Y1 = S1.dot(_np.random.randn(S1.shape[1], D1)) Y2 = S2.dot(_np.random.randn(S2.shape[1], D2)) Y3 = S3.dot(_np.random.randn(S3.shape[1], D3)) Y1 += 0.3 * _np.random.randn(*Y1.shape) Y2 += 0.2 * _np.random.randn(*Y2.shape) Y3 += 0.25 * _np.random.randn(*Y3.shape) Y1 -= Y1.mean(0) Y2 -= Y2.mean(0) Y3 -= Y3.mean(0) Y1 /= Y1.std(0) Y2 /= Y2.std(0) Y3 /= Y3.std(0) return Y1, Y2, Y3, S1, S2, S3
[docs]def bgplvm_simulation( optimize=True, verbose=1, plot=True, plot_sim=False, max_iters=2e4, ): from GPy import kern from GPy.models import BayesianGPLVM D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9 _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim) Y = Ylist[0] k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) # k = kern.RBF(Q, ARD=True, lengthscale=10.) m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k) m.X.variance[:] = _np.random.uniform(0, 0.01, m.X.shape) m.likelihood.variance = 0.1 if optimize: print("Optimizing model:") m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: m.X.plot("BGPLVM Latent Space 1D") m.kern.plot_ARD() return m
[docs]def gplvm_simulation( optimize=True, verbose=1, plot=True, plot_sim=False, max_iters=2e4, ): from GPy import kern from GPy.models import GPLVM D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9 _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim) Y = Ylist[0] k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) # k = kern.RBF(Q, ARD=True, lengthscale=10.) m = GPLVM(Y, Q, init="PCA", kernel=k) m.likelihood.variance = 0.1 if optimize: print("Optimizing model:") m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: m.X.plot("BGPLVM Latent Space 1D") m.kern.plot_ARD() return m
[docs]def ssgplvm_simulation( optimize=True, verbose=1, plot=True, plot_sim=False, max_iters=2e4, useGPU=False ): from GPy import kern from GPy.models import SSGPLVM D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9 _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim) Y = Ylist[0] k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) # k = kern.RBF(Q, ARD=True, lengthscale=10.) m = SSGPLVM( Y, Q, init="rand", num_inducing=num_inducing, kernel=k, group_spike=True ) m.X.variance[:] = _np.random.uniform(0, 0.01, m.X.shape) m.likelihood.variance = 0.01 if optimize: print("Optimizing model:") m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: m.X.plot("SSGPLVM Latent Space 1D") m.kern.plot_ARD() return m
[docs]def bgplvm_simulation_missing_data( optimize=True, verbose=1, plot=True, plot_sim=False, max_iters=2e4, percent_missing=0.1, d=13, ): from GPy import kern from GPy.models.bayesian_gplvm_minibatch import BayesianGPLVMMiniBatch D1, D2, D3, N, num_inducing, Q = d, 5, 8, 400, 3, 4 _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim) Y = Ylist[0] k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) inan = _np.random.binomial(1, percent_missing, size=Y.shape).astype( bool ) # 80% missing data Ymissing = Y.copy() Ymissing[inan] = _np.nan m = BayesianGPLVMMiniBatch( Ymissing, Q, init="random", num_inducing=num_inducing, kernel=k, missing_data=True, ) m.Yreal = Y if optimize: print("Optimizing model:") m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: m.X.plot("BGPLVM Latent Space 1D") m.kern.plot_ARD() return m
[docs]def bgplvm_simulation_missing_data_stochastics( optimize=True, verbose=1, plot=True, plot_sim=False, max_iters=2e4, percent_missing=0.1, d=13, batchsize=2, ): from GPy import kern from GPy.models.bayesian_gplvm_minibatch import BayesianGPLVMMiniBatch D1, D2, D3, N, num_inducing, Q = d, 5, 8, 400, 3, 4 _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim) Y = Ylist[0] k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) inan = _np.random.binomial(1, percent_missing, size=Y.shape).astype( bool ) # 80% missing data Ymissing = Y.copy() Ymissing[inan] = _np.nan m = BayesianGPLVMMiniBatch( Ymissing, Q, init="random", num_inducing=num_inducing, kernel=k, missing_data=True, stochastic=True, batchsize=batchsize, ) m.Yreal = Y if optimize: print("Optimizing model:") m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05) if plot: m.X.plot("BGPLVM Latent Space 1D") m.kern.plot_ARD() return m
[docs]def mrd_simulation(optimize=True, verbose=True, plot=True, plot_sim=True, **kw): from GPy import kern from GPy.models import MRD D1, D2, D3, N, num_inducing, Q = 60, 20, 36, 60, 6, 5 _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim) k = kern.Linear(Q, ARD=True) + kern.White(Q, variance=1e-4) m = MRD( Ylist, input_dim=Q, num_inducing=num_inducing, kernel=k, initx="PCA_concat", initz="permute", **kw ) m[".*noise"] = [Y.var() / 40.0 for Y in Ylist] if optimize: print("Optimizing Model:") m.optimize(messages=verbose, max_iters=8e3) if plot: m.X.plot("MRD Latent Space 1D") m.plot_scales() return m
[docs]def mrd_simulation_missing_data( optimize=True, verbose=True, plot=True, plot_sim=True, **kw ): from GPy import kern from GPy.models import MRD D1, D2, D3, N, num_inducing, Q = 60, 20, 36, 60, 6, 5 _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim) k = kern.Linear(Q, ARD=True) + kern.White(Q, variance=1e-4) inanlist = [] for Y in Ylist: inan = _np.random.binomial(1, 0.6, size=Y.shape).astype(bool) inanlist.append(inan) Y[inan] = _np.nan m = MRD( Ylist, input_dim=Q, num_inducing=num_inducing, kernel=k, inference_method=None, initx="random", initz="permute", **kw ) if optimize: print("Optimizing Model:") m.optimize("bfgs", messages=verbose, max_iters=8e3, gtol=0.1) if plot: m.X.plot("MRD Latent Space 1D") m.plot_scales() return m
[docs]def brendan_faces(optimize=True, verbose=True, plot=True): import GPy import pods data = pods.datasets.brendan_faces() Q = 2 Y = data["Y"] Yn = Y - Y.mean() Yn /= Yn.std() m = GPy.models.BayesianGPLVM(Yn, Q, num_inducing=20) # optimize if optimize: m.optimize("bfgs", messages=verbose, max_iters=1000) if plot: ax = m.plot_latent(which_indices=(0, 1)) y = m.Y[0, :] data_show = GPy.plotting.matplot_dep.visualize.image_show( y[None, :], dimensions=(20, 28), transpose=True, order="F", invert=False, scale=False, ) lvm = GPy.plotting.matplot_dep.visualize.lvm( m.X.mean[0, :].copy(), m, data_show, ax ) input("Press enter to finish") return m
[docs]def olivetti_faces(optimize=True, verbose=True, plot=True): import GPy import pods data = pods.datasets.olivetti_faces() Q = 2 Y = data["Y"] Yn = Y - Y.mean() Yn /= Yn.std() m = GPy.models.BayesianGPLVM(Yn, Q, num_inducing=20) if optimize: m.optimize("bfgs", messages=verbose, max_iters=1000) if plot: ax = m.plot_latent(which_indices=(0, 1)) y = m.Y[0, :] data_show = GPy.plotting.matplot_dep.visualize.image_show( y[None, :], dimensions=(112, 92), transpose=False, invert=False, scale=False ) lvm = GPy.plotting.matplot_dep.visualize.lvm( m.X.mean[0, :].copy(), m, data_show, ax ) input("Press enter to finish") return m
[docs]def stick_play(range=None, frame_rate=15, optimize=False, verbose=True, plot=True): import GPy import pods data = pods.datasets.osu_run1() # optimize if range == None: Y = data["Y"].copy() else: Y = data["Y"][range[0] : range[1], :].copy() if plot: y = Y[0, :] data_show = GPy.plotting.matplot_dep.visualize.stick_show( y[None, :], connect=data["connect"] ) GPy.plotting.matplot_dep.visualize.data_play(Y, data_show, frame_rate) return Y
[docs]def stick(kernel=None, optimize=True, verbose=True, plot=True): from matplotlib import pyplot as plt import GPy import pods data = pods.datasets.osu_run1() # optimize m = GPy.models.GPLVM(data["Y"], 2, kernel=kernel) if optimize: m.optimize("bfgs", messages=verbose, max_f_eval=10000) if plot: plt.clf ax = m.plot_latent() y = m.Y[0, :] data_show = GPy.plotting.matplot_dep.visualize.stick_show( y[None, :], connect=data["connect"] ) lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm( m.X[:1, :].copy(), m, data_show, latent_axes=ax ) input("Press enter to finish") lvm_visualizer.close() data_show.close() return m
[docs]def bcgplvm_linear_stick(kernel=None, optimize=True, verbose=True, plot=True): from matplotlib import pyplot as plt import GPy import pods data = pods.datasets.osu_run1() # optimize mapping = GPy.mappings.Linear(data["Y"].shape[1], 2) m = GPy.models.BCGPLVM(data["Y"], 2, kernel=kernel, mapping=mapping) if optimize: m.optimize(messages=verbose, max_f_eval=10000) if plot and GPy.plotting.matplot_dep.visualize.visual_available: plt.clf ax = m.plot_latent() y = m.likelihood.Y[0, :] data_show = GPy.plotting.matplot_dep.visualize.stick_show( y[None, :], connect=data["connect"] ) GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax) input("Press enter to finish") return m
[docs]def bcgplvm_stick(kernel=None, optimize=True, verbose=True, plot=True): from matplotlib import pyplot as plt import GPy import pods data = pods.datasets.osu_run1() # optimize back_kernel = GPy.kern.RBF(data["Y"].shape[1], lengthscale=5.0) mapping = GPy.mappings.Kernel(X=data["Y"], output_dim=2, kernel=back_kernel) m = GPy.models.BCGPLVM(data["Y"], 2, kernel=kernel, mapping=mapping) if optimize: m.optimize(messages=verbose, max_f_eval=10000) if plot and GPy.plotting.matplot_dep.visualize.visual_available: plt.clf ax = m.plot_latent() y = m.likelihood.Y[0, :] data_show = GPy.plotting.matplot_dep.visualize.stick_show( y[None, :], connect=data["connect"] ) GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax) # input('Press enter to finish') return m
[docs]def robot_wireless(optimize=True, verbose=True, plot=True): from matplotlib import pyplot as plt import GPy import pods data = pods.datasets.robot_wireless() # optimize m = GPy.models.BayesianGPLVM(data["Y"], 4, num_inducing=25) if optimize: m.optimize(messages=verbose, max_f_eval=10000) if plot: m.plot_latent() return m
[docs]def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True): """Interactive visualisation of the Stick Man data from Ohio State University with the Bayesian GPLVM.""" from GPy.models import BayesianGPLVM from matplotlib import pyplot as plt import numpy as np import GPy import pods data = pods.datasets.osu_run1() Q = 6 kernel = GPy.kern.RBF(Q, lengthscale=np.repeat(0.5, Q), ARD=True) m = BayesianGPLVM(data["Y"], Q, init="PCA", num_inducing=20, kernel=kernel) m.data = data m.likelihood.variance = 0.001 # optimize try: if optimize: m.optimize("bfgs", messages=verbose, max_iters=5e3, bfgs_factor=10) except KeyboardInterrupt: print("Keyboard interrupt, continuing to plot and return") if plot: fig, (latent_axes, sense_axes) = plt.subplots(1, 2) plt.sca(latent_axes) m.plot_latent(ax=latent_axes) y = m.Y[:1, :].copy() data_show = GPy.plotting.matplot_dep.visualize.stick_show( y, connect=data["connect"] ) dim_select = GPy.plotting.matplot_dep.visualize.lvm_dimselect( m.X.mean[:1, :].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes, ) fig.canvas.draw() # Canvas.show doesn't work on OSX. # fig.canvas.show() input("Press enter to finish") return m
[docs]def cmu_mocap( subject="35", motion=["01"], in_place=True, optimize=True, verbose=True, plot=True ): import matplotlib.pyplot as plt import GPy import pods data = pods.datasets.cmu_mocap(subject, motion) if in_place: # Make figure move in place. data["Y"][:, 0:3] = 0.0 Y = data["Y"] m = GPy.models.GPLVM(Y, 2, normalizer=True) if optimize: m.optimize(messages=verbose, max_f_eval=10000) if plot: fig, _ = plt.subplots(figsize=(8, 5)) latent_axes = fig.add_subplot(131) sense_axes = fig.add_subplot(132) viz_axes = fig.add_subplot(133, projection="3d") m.plot_latent(ax=latent_axes) latent_axes.set_aspect("equal") y = m.Y[0, :] data_show = GPy.plotting.matplot_dep.visualize.skeleton_show( y[None, :], data["skel"], viz_axes ) lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm( m.X[0].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes ) input("Press enter to finish") lvm_visualizer.close() data_show.close() return m
[docs]def ssgplvm_simulation_linear(): import numpy as np import GPy N, D, Q = 1000, 20, 5 pi = 0.2 def sample_X(Q, pi): x = np.empty(Q) dies = np.random.rand(Q) for q in range(Q): if dies[q] < pi: x[q] = np.random.randn() else: x[q] = 0.0 return x Y = np.empty((N, D)) X = np.empty((N, Q)) # Generate data from random sampled weight matrices for n in range(N): X[n] = sample_X(Q, pi) w = np.random.randn(D, Q) Y[n] = np.dot(w, X[n])