1# Copyright (c) 2012-2014, GPy authors (see AUTHORS.txt). 2# Licensed under the BSD 3-clause license (see LICENSE.txt) 3import numpy as _np 4default_seed = 123344 5 6# default_seed = _np.random.seed(123344) 7 8def bgplvm_test_model(optimize=False, verbose=1, plot=False, output_dim=200, nan=False): 9 """ 10 model for testing purposes. Samples from a GP with rbf kernel and learns 11 the samples with a new kernel. Normally not for optimization, just model cheking 12 """ 13 import GPy 14 15 num_inputs = 13 16 num_inducing = 5 17 if plot: 18 output_dim = 1 19 input_dim = 3 20 else: 21 input_dim = 2 22 output_dim = output_dim 23 24 # generate GPLVM-like data 25 X = _np.random.rand(num_inputs, input_dim) 26 lengthscales = _np.random.rand(input_dim) 27 k = GPy.kern.RBF(input_dim, .5, lengthscales, ARD=True) 28 K = k.K(X) 29 Y = _np.random.multivariate_normal(_np.zeros(num_inputs), K, (output_dim,)).T 30 31 # 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) 32 # k = GPy.kern.linear(input_dim)# + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001) 33 # k = GPy.kern.RBF(input_dim, ARD = False) + GPy.kern.white(input_dim, 0.00001) 34 # 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) 35 # k = GPy.kern.RBF(input_dim, .5, 2., ARD=0) + GPy.kern.RBF(input_dim, .3, .2, ARD=0) 36 # 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) 37 38 p = .3 39 40 m = GPy.models.BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing) 41 42 if nan: 43 m.inference_method = GPy.inference.latent_function_inference.var_dtc.VarDTCMissingData() 44 m.Y[_np.random.binomial(1, p, size=(Y.shape)).astype(bool)] = _np.nan 45 m.parameters_changed() 46 47 #=========================================================================== 48 # randomly obstruct data with percentage p 49 #=========================================================================== 50 # m2 = GPy.models.BayesianGPLVMWithMissingData(Y_obstruct, input_dim, kernel=k, num_inducing=num_inducing) 51 # m.lengthscales = lengthscales 52 53 if plot: 54 import matplotlib.pyplot as pb 55 m.plot() 56 pb.title('PCA initialisation') 57 # m2.plot() 58 # pb.title('PCA initialisation') 59 60 if optimize: 61 m.optimize('scg', messages=verbose) 62 # m2.optimize('scg', messages=verbose) 63 if plot: 64 m.plot() 65 pb.title('After optimisation') 66 # m2.plot() 67 # pb.title('After optimisation') 68 69 return m 70 71def gplvm_oil_100(optimize=True, verbose=1, plot=True): 72 import GPy 73 import pods 74 data = pods.datasets.oil_100() 75 Y = data['X'] 76 # create simple GP model 77 kernel = GPy.kern.RBF(6, ARD=True) + GPy.kern.Bias(6) 78 m = GPy.models.GPLVM(Y, 6, kernel=kernel) 79 m.data_labels = data['Y'].argmax(axis=1) 80 if optimize: m.optimize('scg', messages=verbose) 81 if plot: m.plot_latent(labels=m.data_labels) 82 return m 83 84def sparse_gplvm_oil(optimize=True, verbose=0, plot=True, N=100, Q=6, num_inducing=15, max_iters=50): 85 import GPy 86 import pods 87 88 _np.random.seed(0) 89 data = pods.datasets.oil() 90 Y = data['X'][:N] 91 Y = Y - Y.mean(0) 92 Y /= Y.std(0) 93 # Create the model 94 kernel = GPy.kern.RBF(Q, ARD=True) + GPy.kern.Bias(Q) 95 m = GPy.models.SparseGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing) 96 m.data_labels = data['Y'][:N].argmax(axis=1) 97 98 if optimize: m.optimize('scg', messages=verbose, max_iters=max_iters) 99 if plot: 100 m.plot_latent(labels=m.data_labels) 101 m.kern.plot_ARD() 102 return m 103 104def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=25, Q=4, sigma=.2): 105 import GPy 106 from pods.datasets import swiss_roll_generated 107 from GPy.models import BayesianGPLVM 108 109 data = swiss_roll_generated(num_samples=N, sigma=sigma) 110 Y = data['Y'] 111 Y -= Y.mean() 112 Y /= Y.std() 113 114 t = data['t'] 115 c = data['colors'] 116 117 try: 118 from sklearn.manifold.isomap import Isomap 119 iso = Isomap().fit(Y) 120 X = iso.embedding_ 121 if Q > 2: 122 X = _np.hstack((X, _np.random.randn(N, Q - 2))) 123 except ImportError: 124 X = _np.random.randn(N, Q) 125 126 if plot: 127 import matplotlib.pyplot as plt 128 from mpl_toolkits.mplot3d import Axes3D # @UnusedImport 129 fig = plt.figure("Swiss Roll Data") 130 ax = fig.add_subplot(121, projection='3d') 131 ax.scatter(*Y.T, c=c) 132 ax.set_title("Swiss Roll") 133 134 ax = fig.add_subplot(122) 135 ax.scatter(*X.T[:2], c=c) 136 ax.set_title("BGPLVM init") 137 138 var = .5 139 S = (var * _np.ones_like(X) + _np.clip(_np.random.randn(N, Q) * var ** 2, 140 - (1 - var), 141 (1 - var))) + .001 142 Z = _np.random.permutation(X)[:num_inducing] 143 144 kernel = GPy.kern.RBF(Q, ARD=True) + GPy.kern.Bias(Q, _np.exp(-2)) + GPy.kern.White(Q, _np.exp(-2)) 145 146 m = BayesianGPLVM(Y, Q, X=X, X_variance=S, num_inducing=num_inducing, Z=Z, kernel=kernel) 147 m.data_colors = c 148 m.data_t = t 149 150 if optimize: 151 m.optimize('bfgs', messages=verbose, max_iters=2e3) 152 153 if plot: 154 fig = plt.figure('fitted') 155 ax = fig.add_subplot(111) 156 s = m.input_sensitivity().argsort()[::-1][:2] 157 ax.scatter(*m.X.mean.T[s], c=c) 158 159 return m 160 161def bgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40, max_iters=1000, **k): 162 import GPy 163 from matplotlib import pyplot as plt 164 import numpy as np 165 _np.random.seed(0) 166 try: 167 import pods 168 data = pods.datasets.oil() 169 except ImportError: 170 data = GPy.util.datasets.oil() 171 172 173 kernel = GPy.kern.RBF(Q, 1., 1. / _np.random.uniform(0, 1, (Q,)), ARD=True) # + GPy.kern.Bias(Q, _np.exp(-2)) 174 Y = data['X'][:N] 175 m = GPy.models.BayesianGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing, **k) 176 m.data_labels = data['Y'][:N].argmax(axis=1) 177 178 if optimize: 179 m.optimize('bfgs', messages=verbose, max_iters=max_iters, gtol=.05) 180 181 if plot: 182 fig, (latent_axes, sense_axes) = plt.subplots(1, 2) 183 m.plot_latent(ax=latent_axes, labels=m.data_labels) 184 data_show = GPy.plotting.matplot_dep.visualize.vector_show((m.Y[0, :])) 185 lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X.mean.values[0:1, :], # @UnusedVariable 186 m, data_show, latent_axes=latent_axes, sense_axes=sense_axes, labels=m.data_labels) 187 input('Press enter to finish') 188 plt.close(fig) 189 return m 190 191def ssgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40, max_iters=1000, **k): 192 import GPy 193 from matplotlib import pyplot as plt 194 import pods 195 196 _np.random.seed(0) 197 data = pods.datasets.oil() 198 199 kernel = GPy.kern.RBF(Q, 1., 1. / _np.random.uniform(0, 1, (Q,)), ARD=True) # + GPy.kern.Bias(Q, _np.exp(-2)) 200 Y = data['X'][:N] 201 m = GPy.models.SSGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing, **k) 202 m.data_labels = data['Y'][:N].argmax(axis=1) 203 204 if optimize: 205 m.optimize('bfgs', messages=verbose, max_iters=max_iters, gtol=.05) 206 207 if plot: 208 fig, (latent_axes, sense_axes) = plt.subplots(1, 2) 209 m.plot_latent(ax=latent_axes, labels=m.data_labels) 210 data_show = GPy.plotting.matplot_dep.visualize.vector_show((m.Y[0, :])) 211 lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X.mean.values[0:1, :], # @UnusedVariable 212 m, data_show, latent_axes=latent_axes, sense_axes=sense_axes, labels=m.data_labels) 213 input('Press enter to finish') 214 plt.close(fig) 215 return m 216 217def _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim=False): 218 """Simulate some data drawn from a matern covariance and a periodic exponential for use in MRD demos.""" 219 Q_signal = 4 220 import GPy 221 import numpy as np 222 np.random.seed(3000) 223 224 k = GPy.kern.Matern32(Q_signal, 1., lengthscale=(np.random.uniform(1, 6, Q_signal)), ARD=1) 225 for i in range(Q_signal): 226 k += GPy.kern.PeriodicExponential(1, variance=1., active_dims=[i], period=3., lower=-2, upper=6) 227 t = np.c_[[np.linspace(-1, 5, N) for _ in range(Q_signal)]].T 228 K = k.K(t) 229 s2, s1, s3, sS = np.random.multivariate_normal(np.zeros(K.shape[0]), K, size=(4))[:, :, None] 230 231 Y1, Y2, Y3, S1, S2, S3 = _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS) 232 233 slist = [sS, s1, s2, s3] 234 slist_names = ["sS", "s1", "s2", "s3"] 235 Ylist = [Y1, Y2, Y3] 236 237 if plot_sim: 238 from matplotlib import pyplot as plt 239 import matplotlib.cm as cm 240 import itertools 241 fig = plt.figure("MRD Simulation Data", figsize=(8, 6)) 242 fig.clf() 243 ax = fig.add_subplot(2, 1, 1) 244 labls = slist_names 245 for S, lab in zip(slist, labls): 246 ax.plot(S, label=lab) 247 ax.legend() 248 for i, Y in enumerate(Ylist): 249 ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i) 250 ax.imshow(Y, aspect='auto', cmap=cm.gray) # @UndefinedVariable 251 ax.set_title("Y{}".format(i + 1)) 252 plt.draw() 253 plt.tight_layout() 254 255 return slist, [S1, S2, S3], Ylist 256 257def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False): 258 """Simulate some data drawn from sine and cosine for use in demos of MRD""" 259 _np.random.seed(1234) 260 261 x = _np.linspace(0, 4 * _np.pi, N)[:, None] 262 s1 = _np.vectorize(lambda x: _np.sin(x)) 263 s2 = _np.vectorize(lambda x: _np.cos(x)) 264 s3 = _np.vectorize(lambda x:-_np.exp(-_np.cos(2 * x))) 265 sS = _np.vectorize(lambda x: _np.cos(x)) 266 267 s1 = s1(x) 268 s2 = s2(x) 269 s3 = s3(x) 270 sS = sS(x) 271 272 s1 -= s1.mean(); s1 /= s1.std(0) 273 s2 -= s2.mean(); s2 /= s2.std(0) 274 s3 -= s3.mean(); s3 /= s3.std(0) 275 sS -= sS.mean(); sS /= sS.std(0) 276 277 Y1, Y2, Y3, S1, S2, S3 = _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS) 278 279 slist = [sS, s1, s2, s3] 280 slist_names = ["sS", "s1", "s2", "s3"] 281 Ylist = [Y1, Y2, Y3] 282 283 if plot_sim: 284 from matplotlib import pyplot as plt 285 import matplotlib.cm as cm 286 import itertools 287 fig = plt.figure("MRD Simulation Data", figsize=(8, 6)) 288 fig.clf() 289 ax = fig.add_subplot(2, 1, 1) 290 labls = slist_names 291 for S, lab in zip(slist, labls): 292 ax.plot(S, label=lab) 293 ax.legend() 294 for i, Y in enumerate(Ylist): 295 ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i) 296 ax.imshow(Y, aspect='auto', cmap=cm.gray) # @UndefinedVariable 297 ax.set_title("Y{}".format(i + 1)) 298 plt.draw() 299 plt.tight_layout() 300 301 return slist, [S1, S2, S3], Ylist 302 303def _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS): 304 S1 = _np.hstack([s1, sS]) 305 S2 = _np.hstack([sS]) 306 S3 = _np.hstack([s1, s3, sS]) 307 Y1 = S1.dot(_np.random.randn(S1.shape[1], D1)) 308 Y2 = S2.dot(_np.random.randn(S2.shape[1], D2)) 309 Y3 = S3.dot(_np.random.randn(S3.shape[1], D3)) 310 Y1 += .3 * _np.random.randn(*Y1.shape) 311 Y2 += .2 * _np.random.randn(*Y2.shape) 312 Y3 += .25 * _np.random.randn(*Y3.shape) 313 Y1 -= Y1.mean(0) 314 Y2 -= Y2.mean(0) 315 Y3 -= Y3.mean(0) 316 Y1 /= Y1.std(0) 317 Y2 /= Y2.std(0) 318 Y3 /= Y3.std(0) 319 return Y1, Y2, Y3, S1, S2, S3 320 321def bgplvm_simulation(optimize=True, verbose=1, 322 plot=True, plot_sim=False, 323 max_iters=2e4, 324 ): 325 from GPy import kern 326 from GPy.models import BayesianGPLVM 327 328 D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9 329 _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim) 330 Y = Ylist[0] 331 k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) 332 # k = kern.RBF(Q, ARD=True, lengthscale=10.) 333 m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k) 334 m.X.variance[:] = _np.random.uniform(0, .01, m.X.shape) 335 m.likelihood.variance = .1 336 337 if optimize: 338 print("Optimizing model:") 339 m.optimize('bfgs', messages=verbose, max_iters=max_iters, 340 gtol=.05) 341 if plot: 342 m.X.plot("BGPLVM Latent Space 1D") 343 m.kern.plot_ARD() 344 return m 345 346def gplvm_simulation(optimize=True, verbose=1, 347 plot=True, plot_sim=False, 348 max_iters=2e4, 349 ): 350 from GPy import kern 351 from GPy.models import GPLVM 352 353 D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9 354 _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim) 355 Y = Ylist[0] 356 k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) 357 # k = kern.RBF(Q, ARD=True, lengthscale=10.) 358 m = GPLVM(Y, Q, init="PCA", kernel=k) 359 m.likelihood.variance = .1 360 361 if optimize: 362 print("Optimizing model:") 363 m.optimize('bfgs', messages=verbose, max_iters=max_iters, 364 gtol=.05) 365 if plot: 366 m.X.plot("BGPLVM Latent Space 1D") 367 m.kern.plot_ARD() 368 return m 369def ssgplvm_simulation(optimize=True, verbose=1, 370 plot=True, plot_sim=False, 371 max_iters=2e4, useGPU=False 372 ): 373 from GPy import kern 374 from GPy.models import SSGPLVM 375 376 D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9 377 _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim) 378 Y = Ylist[0] 379 k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) 380 # k = kern.RBF(Q, ARD=True, lengthscale=10.) 381 m = SSGPLVM(Y, Q, init="rand", num_inducing=num_inducing, kernel=k, group_spike=True) 382 m.X.variance[:] = _np.random.uniform(0, .01, m.X.shape) 383 m.likelihood.variance = .01 384 385 if optimize: 386 print("Optimizing model:") 387 m.optimize('bfgs', messages=verbose, max_iters=max_iters, 388 gtol=.05) 389 if plot: 390 m.X.plot("SSGPLVM Latent Space 1D") 391 m.kern.plot_ARD() 392 return m 393 394def bgplvm_simulation_missing_data(optimize=True, verbose=1, 395 plot=True, plot_sim=False, 396 max_iters=2e4, percent_missing=.1, d=13, 397 ): 398 from GPy import kern 399 from GPy.models.bayesian_gplvm_minibatch import BayesianGPLVMMiniBatch 400 401 D1, D2, D3, N, num_inducing, Q = d, 5, 8, 400, 3, 4 402 _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim) 403 Y = Ylist[0] 404 k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) 405 406 inan = _np.random.binomial(1, percent_missing, size=Y.shape).astype(bool) # 80% missing data 407 Ymissing = Y.copy() 408 Ymissing[inan] = _np.nan 409 410 m = BayesianGPLVMMiniBatch(Ymissing, Q, init="random", num_inducing=num_inducing, 411 kernel=k, missing_data=True) 412 413 m.Yreal = Y 414 415 if optimize: 416 print("Optimizing model:") 417 m.optimize('bfgs', messages=verbose, max_iters=max_iters, 418 gtol=.05) 419 if plot: 420 m.X.plot("BGPLVM Latent Space 1D") 421 m.kern.plot_ARD() 422 return m 423 424def bgplvm_simulation_missing_data_stochastics(optimize=True, verbose=1, 425 plot=True, plot_sim=False, 426 max_iters=2e4, percent_missing=.1, d=13, batchsize=2, 427 ): 428 from GPy import kern 429 from GPy.models.bayesian_gplvm_minibatch import BayesianGPLVMMiniBatch 430 431 D1, D2, D3, N, num_inducing, Q = d, 5, 8, 400, 3, 4 432 _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim) 433 Y = Ylist[0] 434 k = kern.Linear(Q, ARD=True) # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) 435 436 inan = _np.random.binomial(1, percent_missing, size=Y.shape).astype(bool) # 80% missing data 437 Ymissing = Y.copy() 438 Ymissing[inan] = _np.nan 439 440 m = BayesianGPLVMMiniBatch(Ymissing, Q, init="random", num_inducing=num_inducing, 441 kernel=k, missing_data=True, stochastic=True, batchsize=batchsize) 442 443 m.Yreal = Y 444 445 if optimize: 446 print("Optimizing model:") 447 m.optimize('bfgs', messages=verbose, max_iters=max_iters, 448 gtol=.05) 449 if plot: 450 m.X.plot("BGPLVM Latent Space 1D") 451 m.kern.plot_ARD() 452 return m 453 454 455def mrd_simulation(optimize=True, verbose=True, plot=True, plot_sim=True, **kw): 456 from GPy import kern 457 from GPy.models import MRD 458 459 D1, D2, D3, N, num_inducing, Q = 60, 20, 36, 60, 6, 5 460 _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim) 461 462 k = kern.Linear(Q, ARD=True) + kern.White(Q, variance=1e-4) 463 m = MRD(Ylist, input_dim=Q, num_inducing=num_inducing, kernel=k, initx="PCA_concat", initz='permute', **kw) 464 465 m['.*noise'] = [Y.var() / 40. for Y in Ylist] 466 467 if optimize: 468 print("Optimizing Model:") 469 m.optimize(messages=verbose, max_iters=8e3) 470 if plot: 471 m.X.plot("MRD Latent Space 1D") 472 m.plot_scales() 473 return m 474 475def mrd_simulation_missing_data(optimize=True, verbose=True, plot=True, plot_sim=True, **kw): 476 from GPy import kern 477 from GPy.models import MRD 478 479 D1, D2, D3, N, num_inducing, Q = 60, 20, 36, 60, 6, 5 480 _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim) 481 482 k = kern.Linear(Q, ARD=True) + kern.White(Q, variance=1e-4) 483 inanlist = [] 484 485 for Y in Ylist: 486 inan = _np.random.binomial(1, .6, size=Y.shape).astype(bool) 487 inanlist.append(inan) 488 Y[inan] = _np.nan 489 490 m = MRD(Ylist, input_dim=Q, num_inducing=num_inducing, 491 kernel=k, inference_method=None, 492 initx="random", initz='permute', **kw) 493 494 if optimize: 495 print("Optimizing Model:") 496 m.optimize('bfgs', messages=verbose, max_iters=8e3, gtol=.1) 497 if plot: 498 m.X.plot("MRD Latent Space 1D") 499 m.plot_scales() 500 return m 501 502def brendan_faces(optimize=True, verbose=True, plot=True): 503 import GPy 504 import pods 505 506 data = pods.datasets.brendan_faces() 507 Q = 2 508 Y = data['Y'] 509 Yn = Y - Y.mean() 510 Yn /= Yn.std() 511 512 m = GPy.models.BayesianGPLVM(Yn, Q, num_inducing=20) 513 514 # optimize 515 516 if optimize: m.optimize('bfgs', messages=verbose, max_iters=1000) 517 518 if plot: 519 ax = m.plot_latent(which_indices=(0, 1)) 520 y = m.Y[0, :] 521 data_show = GPy.plotting.matplot_dep.visualize.image_show(y[None, :], dimensions=(20, 28), transpose=True, order='F', invert=False, scale=False) 522 lvm = GPy.plotting.matplot_dep.visualize.lvm(m.X.mean[0, :].copy(), m, data_show, ax) 523 input('Press enter to finish') 524 525 return m 526 527def olivetti_faces(optimize=True, verbose=True, plot=True): 528 import GPy 529 import pods 530 531 data = pods.datasets.olivetti_faces() 532 Q = 2 533 Y = data['Y'] 534 Yn = Y - Y.mean() 535 Yn /= Yn.std() 536 537 m = GPy.models.BayesianGPLVM(Yn, Q, num_inducing=20) 538 539 if optimize: m.optimize('bfgs', messages=verbose, max_iters=1000) 540 if plot: 541 ax = m.plot_latent(which_indices=(0, 1)) 542 y = m.Y[0, :] 543 data_show = GPy.plotting.matplot_dep.visualize.image_show(y[None, :], dimensions=(112, 92), transpose=False, invert=False, scale=False) 544 lvm = GPy.plotting.matplot_dep.visualize.lvm(m.X.mean[0, :].copy(), m, data_show, ax) 545 input('Press enter to finish') 546 547 return m 548 549def stick_play(range=None, frame_rate=15, optimize=False, verbose=True, plot=True): 550 import GPy 551 import pods 552 553 data = pods.datasets.osu_run1() 554 # optimize 555 if range == None: 556 Y = data['Y'].copy() 557 else: 558 Y = data['Y'][range[0]:range[1], :].copy() 559 if plot: 560 y = Y[0, :] 561 data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect']) 562 GPy.plotting.matplot_dep.visualize.data_play(Y, data_show, frame_rate) 563 return Y 564 565def stick(kernel=None, optimize=True, verbose=True, plot=True): 566 from matplotlib import pyplot as plt 567 import GPy 568 import pods 569 570 data = pods.datasets.osu_run1() 571 # optimize 572 m = GPy.models.GPLVM(data['Y'], 2, kernel=kernel) 573 if optimize: m.optimize('bfgs', messages=verbose, max_f_eval=10000) 574 if plot: 575 plt.clf 576 ax = m.plot_latent() 577 y = m.Y[0, :] 578 data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect']) 579 lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm(m.X[:1, :].copy(), m, data_show, latent_axes=ax) 580 input('Press enter to finish') 581 lvm_visualizer.close() 582 data_show.close() 583 return m 584 585def bcgplvm_linear_stick(kernel=None, optimize=True, verbose=True, plot=True): 586 from matplotlib import pyplot as plt 587 import GPy 588 import pods 589 590 data = pods.datasets.osu_run1() 591 # optimize 592 mapping = GPy.mappings.Linear(data['Y'].shape[1], 2) 593 m = GPy.models.BCGPLVM(data['Y'], 2, kernel=kernel, mapping=mapping) 594 if optimize: m.optimize(messages=verbose, max_f_eval=10000) 595 if plot and GPy.plotting.matplot_dep.visualize.visual_available: 596 plt.clf 597 ax = m.plot_latent() 598 y = m.likelihood.Y[0, :] 599 data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect']) 600 GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax) 601 input('Press enter to finish') 602 603 return m 604 605def bcgplvm_stick(kernel=None, optimize=True, verbose=True, plot=True): 606 from matplotlib import pyplot as plt 607 import GPy 608 import pods 609 610 data = pods.datasets.osu_run1() 611 # optimize 612 back_kernel = GPy.kern.RBF(data['Y'].shape[1], lengthscale=5.) 613 mapping = GPy.mappings.Kernel(X=data['Y'], output_dim=2, kernel=back_kernel) 614 m = GPy.models.BCGPLVM(data['Y'], 2, kernel=kernel, mapping=mapping) 615 if optimize: m.optimize(messages=verbose, max_f_eval=10000) 616 if plot and GPy.plotting.matplot_dep.visualize.visual_available: 617 plt.clf 618 ax = m.plot_latent() 619 y = m.likelihood.Y[0, :] 620 data_show = GPy.plotting.matplot_dep.visualize.stick_show(y[None, :], connect=data['connect']) 621 GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax) 622 # input('Press enter to finish') 623 624 return m 625 626def robot_wireless(optimize=True, verbose=True, plot=True): 627 from matplotlib import pyplot as plt 628 import GPy 629 import pods 630 631 data = pods.datasets.robot_wireless() 632 # optimize 633 m = GPy.models.BayesianGPLVM(data['Y'], 4, num_inducing=25) 634 if optimize: m.optimize(messages=verbose, max_f_eval=10000) 635 if plot: 636 m.plot_latent() 637 638 return m 639 640def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True): 641 """Interactive visualisation of the Stick Man data from Ohio State University with the Bayesian GPLVM.""" 642 from GPy.models import BayesianGPLVM 643 from matplotlib import pyplot as plt 644 import numpy as np 645 import GPy 646 import pods 647 648 data = pods.datasets.osu_run1() 649 Q = 6 650 kernel = GPy.kern.RBF(Q, lengthscale=np.repeat(.5, Q), ARD=True) 651 m = BayesianGPLVM(data['Y'], Q, init="PCA", num_inducing=20, kernel=kernel) 652 653 m.data = data 654 m.likelihood.variance = 0.001 655 656 # optimize 657 try: 658 if optimize: m.optimize('bfgs', messages=verbose, max_iters=5e3, bfgs_factor=10) 659 except KeyboardInterrupt: 660 print("Keyboard interrupt, continuing to plot and return") 661 662 if plot: 663 fig, (latent_axes, sense_axes) = plt.subplots(1, 2) 664 plt.sca(latent_axes) 665 m.plot_latent(ax=latent_axes) 666 y = m.Y[:1, :].copy() 667 data_show = GPy.plotting.matplot_dep.visualize.stick_show(y, connect=data['connect']) 668 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) 669 fig.canvas.draw() 670 # Canvas.show doesn't work on OSX. 671 #fig.canvas.show() 672 input('Press enter to finish') 673 674 return m 675 676 677def cmu_mocap(subject='35', motion=['01'], in_place=True, optimize=True, verbose=True, plot=True): 678 import GPy 679 import pods 680 681 data = pods.datasets.cmu_mocap(subject, motion) 682 if in_place: 683 # Make figure move in place. 684 data['Y'][:, 0:3] = 0.0 685 Y = data['Y'] 686 Y_mean = Y.mean(0) 687 Y_std = Y.std(0) 688 m = GPy.models.GPLVM((Y - Y_mean) / Y_std, 2) 689 690 if optimize: m.optimize(messages=verbose, max_f_eval=10000) 691 if plot: 692 ax = m.plot_latent() 693 y = m.Y[0, :] 694 data_show = GPy.plotting.matplot_dep.visualize.skeleton_show(y[None, :], data['skel']) 695 lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm(m.X[0].copy(), m, data_show, latent_axes=ax) 696 input('Press enter to finish') 697 lvm_visualizer.close() 698 data_show.close() 699 700 return m 701 702def ssgplvm_simulation_linear(): 703 import numpy as np 704 import GPy 705 N, D, Q = 1000, 20, 5 706 pi = 0.2 707 708 def sample_X(Q, pi): 709 x = np.empty(Q) 710 dies = np.random.rand(Q) 711 for q in range(Q): 712 if dies[q] < pi: 713 x[q] = np.random.randn() 714 else: 715 x[q] = 0. 716 return x 717 718 Y = np.empty((N, D)) 719 X = np.empty((N, Q)) 720 # Generate data from random sampled weight matrices 721 for n in range(N): 722 X[n] = sample_X(Q, pi) 723 w = np.random.randn(D, Q) 724 Y[n] = np.dot(w, X[n]) 725 726