Run black on examples.

GPy/examples/classification.py

@@ -7,39 +7,47 @@ import GPy
 
 default_seed = 10000
 
+
 def oil(num_inducing=50, max_iters=100, kernel=None, optimize=True, plot=True):
     """
     Run a Gaussian process classification on the three phase oil data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood.
 
     """
-    try:import pods
+    try:
-    except ImportError:raise ImportWarning('Need pods for example datasets. See https://github.com/sods/ods, or pip install pods.')
+        import pods
+    except ImportError:
+        raise ImportWarning(
+            "Need pods for example datasets. See https://github.com/sods/ods, or pip install pods."
+        )
     data = pods.datasets.oil()
-    X = data['X']
+    X = data["X"]
-    Xtest = data['Xtest']
+    Xtest = data["Xtest"]
-    Y = data['Y'][:, 0:1]
+    Y = data["Y"][:, 0:1]
-    Ytest = data['Ytest'][:, 0:1]
+    Ytest = data["Ytest"][:, 0:1]
-    Y[Y.flatten()==-1] = 0
+    Y[Y.flatten() == -1] = 0
-    Ytest[Ytest.flatten()==-1] = 0
+    Ytest[Ytest.flatten() == -1] = 0
 
     # Create GP model
-    m = GPy.models.SparseGPClassification(X, Y, kernel=kernel, num_inducing=num_inducing)
+    m = GPy.models.SparseGPClassification(
+        X, Y, kernel=kernel, num_inducing=num_inducing
+    )
     m.Ytest = Ytest
 
     # Contrain all parameters to be positive
-    #m.tie_params('.*len')
+    # m.tie_params('.*len')
-    m['.*len'] = 10.
+    m[".*len"] = 10.0
 
     # Optimize
     if optimize:
         m.optimize(messages=1)
     print(m)
 
-    #Test
+    # Test
     probs = m.predict(Xtest)[0]
     GPy.util.classification.conf_matrix(probs, Ytest)
     return m
 
+
 def toy_linear_1d_classification(seed=default_seed, optimize=True, plot=True):
     """
     Simple 1D classification example using EP approximation
@@ -48,26 +56,31 @@ def toy_linear_1d_classification(seed=default_seed, optimize=True, plot=True):
     :type seed: int
 
     """
-    try:import pods
+    try:
-    except ImportError:raise ImportWarning('Need pods for example datasets. See https://github.com/sods/ods, or pip install pods.')
+        import pods
+    except ImportError:
+        raise ImportWarning(
+            "Need pods for example datasets. See https://github.com/sods/ods, or pip install pods."
+        )
     data = pods.datasets.toy_linear_1d_classification(seed=seed)
-    Y = data['Y'][:, 0:1]
+    Y = data["Y"][:, 0:1]
     Y[Y.flatten() == -1] = 0
 
     # Model definition
-    m = GPy.models.GPClassification(data['X'], Y)
+    m = GPy.models.GPClassification(data["X"], Y)
 
     # Optimize
     if optimize:
-        #m.update_likelihood_approximation()
+        # m.update_likelihood_approximation()
         # Parameters optimization:
         m.optimize()
-        #m.update_likelihood_approximation()
+        # m.update_likelihood_approximation()
-        #m.pseudo_EM()
+        # m.pseudo_EM()
 
     # Plot
     if plot:
         from matplotlib import pyplot as plt
+
         fig, axes = plt.subplots(2, 1)
         m.plot_f(ax=axes[0])
         m.plot(ax=axes[1])
@@ -75,6 +88,7 @@ def toy_linear_1d_classification(seed=default_seed, optimize=True, plot=True):
     print(m)
     return m
 
+
 def toy_linear_1d_classification_laplace(seed=default_seed, optimize=True, plot=True):
     """
     Simple 1D classification example using Laplace approximation
@@ -84,10 +98,12 @@ def toy_linear_1d_classification_laplace(seed=default_seed, optimize=True, plot=
 
     """
 
-    try:import pods
+    try:
-    except ImportError:print('pods unavailable, see https://github.com/sods/ods for example datasets')
+        import pods
+    except ImportError:
+        print("pods unavailable, see https://github.com/sods/ods for example datasets")
     data = pods.datasets.toy_linear_1d_classification(seed=seed)
-    Y = data['Y'][:, 0:1]
+    Y = data["Y"][:, 0:1]
     Y[Y.flatten() == -1] = 0
 
     likelihood = GPy.likelihoods.Bernoulli()
@@ -95,18 +111,21 @@ def toy_linear_1d_classification_laplace(seed=default_seed, optimize=True, plot=
     kernel = GPy.kern.RBF(1)
 
     # Model definition
-    m = GPy.core.GP(data['X'], Y, kernel=kernel, likelihood=likelihood, inference_method=laplace_inf)
+    m = GPy.core.GP(
+        data["X"], Y, kernel=kernel, likelihood=likelihood, inference_method=laplace_inf
+    )
 
     # Optimize
     if optimize:
         try:
-            m.optimize('scg', messages=1)
+            m.optimize("scg", messages=1)
         except Exception as e:
             return m
 
     # Plot
     if plot:
         from matplotlib import pyplot as plt
+
         fig, axes = plt.subplots(2, 1)
         m.plot_f(ax=axes[0])
         m.plot(ax=axes[1])
@@ -114,7 +133,10 @@ def toy_linear_1d_classification_laplace(seed=default_seed, optimize=True, plot=
     print(m)
     return m
 
-def sparse_toy_linear_1d_classification(num_inducing=10, seed=default_seed, optimize=True, plot=True):
+
+def sparse_toy_linear_1d_classification(
+    num_inducing=10, seed=default_seed, optimize=True, plot=True
+):
     """
     Sparse 1D classification example
 
@@ -123,15 +145,17 @@ def sparse_toy_linear_1d_classification(num_inducing=10, seed=default_seed, opti
 
     """
 
-    try:import pods
+    try:
-    except ImportError:print('pods unavailable, see https://github.com/sods/ods for example datasets')
+        import pods
+    except ImportError:
+        print("pods unavailable, see https://github.com/sods/ods for example datasets")
     data = pods.datasets.toy_linear_1d_classification(seed=seed)
-    Y = data['Y'][:, 0:1]
+    Y = data["Y"][:, 0:1]
     Y[Y.flatten() == -1] = 0
 
     # Model definition
-    m = GPy.models.SparseGPClassification(data['X'], Y, num_inducing=num_inducing)
+    m = GPy.models.SparseGPClassification(data["X"], Y, num_inducing=num_inducing)
-    m['.*len'] = 4.
+    m[".*len"] = 4.0
 
     # Optimize
     if optimize:
@@ -140,6 +164,7 @@ def sparse_toy_linear_1d_classification(num_inducing=10, seed=default_seed, opti
     # Plot
     if plot:
         from matplotlib import pyplot as plt
+
         fig, axes = plt.subplots(2, 1)
         m.plot_f(ax=axes[0])
         m.plot(ax=axes[1])
@@ -147,7 +172,10 @@ def sparse_toy_linear_1d_classification(num_inducing=10, seed=default_seed, opti
     print(m)
     return m
 
-def sparse_toy_linear_1d_classification_uncertain_input(num_inducing=10, seed=default_seed, optimize=True, plot=True):
+
+def sparse_toy_linear_1d_classification_uncertain_input(
+    num_inducing=10, seed=default_seed, optimize=True, plot=True
+):
     """
     Sparse 1D classification example
 
@@ -156,18 +184,23 @@ def sparse_toy_linear_1d_classification_uncertain_input(num_inducing=10, seed=de
 
     """
 
-    try:import pods
+    try:
-    except ImportError:print('pods unavailable, see https://github.com/sods/ods for example datasets')
+        import pods
+    except ImportError:
+        print("pods unavailable, see https://github.com/sods/ods for example datasets")
     import numpy as np
+
     data = pods.datasets.toy_linear_1d_classification(seed=seed)
-    Y = data['Y'][:, 0:1]
+    Y = data["Y"][:, 0:1]
     Y[Y.flatten() == -1] = 0
-    X = data['X']
+    X = data["X"]
-    X_var = np.random.uniform(0.3,0.5,X.shape)
+    X_var = np.random.uniform(0.3, 0.5, X.shape)
 
     # Model definition
-    m = GPy.models.SparseGPClassificationUncertainInput(X, X_var, Y, num_inducing=num_inducing)
+    m = GPy.models.SparseGPClassificationUncertainInput(
-    m['.*len'] = 4.
+        X, X_var, Y, num_inducing=num_inducing
+    )
+    m[".*len"] = 4.0
 
     # Optimize
     if optimize:
@@ -176,6 +209,7 @@ def sparse_toy_linear_1d_classification_uncertain_input(num_inducing=10, seed=de
     # Plot
     if plot:
         from matplotlib import pyplot as plt
+
         fig, axes = plt.subplots(2, 1)
         m.plot_f(ax=axes[0])
         m.plot(ax=axes[1])
@@ -183,6 +217,7 @@ def sparse_toy_linear_1d_classification_uncertain_input(num_inducing=10, seed=de
     print(m)
     return m
 
+
 def toy_heaviside(seed=default_seed, max_iters=100, optimize=True, plot=True):
     """
     Simple 1D classification example using a heavy side gp transformation
@@ -192,29 +227,41 @@ def toy_heaviside(seed=default_seed, max_iters=100, optimize=True, plot=True):
 
     """
 
-    try:import pods
+    try:
-    except ImportError:print('pods unavailable, see https://github.com/sods/ods for example datasets')
+        import pods
+    except ImportError:
+        print("pods unavailable, see https://github.com/sods/ods for example datasets")
     data = pods.datasets.toy_linear_1d_classification(seed=seed)
-    Y = data['Y'][:, 0:1]
+    Y = data["Y"][:, 0:1]
     Y[Y.flatten() == -1] = 0
 
     # Model definition
     kernel = GPy.kern.RBF(1)
-    likelihood = GPy.likelihoods.Bernoulli(gp_link=GPy.likelihoods.link_functions.Heaviside())
+    likelihood = GPy.likelihoods.Bernoulli(
+        gp_link=GPy.likelihoods.link_functions.Heaviside()
+    )
     ep = GPy.inference.latent_function_inference.expectation_propagation.EP()
-    m = GPy.core.GP(X=data['X'], Y=Y, kernel=kernel, likelihood=likelihood, inference_method=ep, name='gp_classification_heaviside')
+    m = GPy.core.GP(
-    #m = GPy.models.GPClassification(data['X'], likelihood=likelihood)
+        X=data["X"],
+        Y=Y,
+        kernel=kernel,
+        likelihood=likelihood,
+        inference_method=ep,
+        name="gp_classification_heaviside",
+    )
+    # m = GPy.models.GPClassification(data['X'], likelihood=likelihood)
 
     # Optimize
     if optimize:
         # Parameters optimization:
         for _ in range(5):
-            m.optimize(max_iters=int(max_iters/5))
+            m.optimize(max_iters=int(max_iters / 5))
         print(m)
 
     # Plot
     if plot:
         from matplotlib import pyplot as plt
+
         fig, axes = plt.subplots(2, 1)
         m.plot_f(ax=axes[0])
         m.plot(ax=axes[1])
@@ -222,7 +269,15 @@ def toy_heaviside(seed=default_seed, max_iters=100, optimize=True, plot=True):
     print(m)
     return m
 
-def crescent_data(model_type='Full', num_inducing=10, seed=default_seed, kernel=None, optimize=True, plot=True):
+
+def crescent_data(
+    model_type="Full",
+    num_inducing=10,
+    seed=default_seed,
+    kernel=None,
+    optimize=True,
+    plot=True,
+):
     """
     Run a Gaussian process classification on the crescent data. The demonstration calls the basic GP classification model and uses EP to approximate the likelihood.
 
@@ -234,22 +289,28 @@ def crescent_data(model_type='Full', num_inducing=10, seed=default_seed, kernel=
     :param kernel: kernel to use in the model
     :type kernel: a GPy kernel
     """
-    try:import pods
+    try:
-    except ImportError:print('pods unavailable, see https://github.com/sods/ods for example datasets')
+        import pods
+    except ImportError:
+        print("pods unavailable, see https://github.com/sods/ods for example datasets")
     data = pods.datasets.crescent_data(seed=seed)
-    Y = data['Y']
+    Y = data["Y"]
-    Y[Y.flatten()==-1] = 0
+    Y[Y.flatten() == -1] = 0
 
-    if model_type == 'Full':
+    if model_type == "Full":
-        m = GPy.models.GPClassification(data['X'], Y, kernel=kernel)
+        m = GPy.models.GPClassification(data["X"], Y, kernel=kernel)
 
-    elif model_type == 'DTC':
+    elif model_type == "DTC":
-        m = GPy.models.SparseGPClassification(data['X'], Y, kernel=kernel, num_inducing=num_inducing)
+        m = GPy.models.SparseGPClassification(
-        m['.*len'] = 10.
+            data["X"], Y, kernel=kernel, num_inducing=num_inducing
+        )
+        m[".*len"] = 10.0
 
-    elif model_type == 'FITC':
+    elif model_type == "FITC":
-        m = GPy.models.FITCClassification(data['X'], Y, kernel=kernel, num_inducing=num_inducing)
+        m = GPy.models.FITCClassification(
-        m['.*len'] = 3.
+            data["X"], Y, kernel=kernel, num_inducing=num_inducing
+        )
+        m[".*len"] = 3.0
     if optimize:
         m.optimize(messages=1)
 

GPy/examples/dimensionality_reduction.py

@@ -1,10 +1,12 @@
 # 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)
 
+
 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
@@ -24,7 +26,7 @@ def bgplvm_test_model(optimize=False, verbose=1, plot=False, output_dim=200, nan
     # generate GPLVM-like data
     X = _np.random.rand(num_inputs, input_dim)
     lengthscales = _np.random.rand(input_dim)
-    k = GPy.kern.RBF(input_dim, .5, lengthscales, ARD=True)
+    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
 
@@ -35,88 +37,102 @@ def bgplvm_test_model(optimize=False, verbose=1, plot=False, output_dim=200, nan
     # 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 = .3
+    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.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')
+        pb.title("PCA initialisation")
         # m2.plot()
         # pb.title('PCA initialisation')
 
     if optimize:
-        m.optimize('scg', messages=verbose)
+        m.optimize("scg", messages=verbose)
         # m2.optimize('scg', messages=verbose)
         if plot:
             m.plot()
-            pb.title('After optimisation')
+            pb.title("After optimisation")
             # m2.plot()
             # pb.title('After optimisation')
 
     return m
 
+
 def gplvm_oil_100(optimize=True, verbose=1, plot=True):
     import GPy
     import pods
+
     data = pods.datasets.oil_100()
-    Y = data['X']
+    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)
+    m.data_labels = data["Y"].argmax(axis=1)
-    if optimize: m.optimize('scg', messages=verbose)
+    if optimize:
+        m.optimize("scg", messages=verbose)
     if plot:
         m.plot_latent(labels=m.data_labels)
     return m
 
-def sparse_gplvm_oil(optimize=True, verbose=0, plot=True, N=100, Q=6, num_inducing=15, max_iters=50):
+
+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 = 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)
+    m.data_labels = data["Y"][:N].argmax(axis=1)
 
-    if optimize: m.optimize('scg', messages=verbose, max_iters=max_iters)
+    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
 
-def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=25, Q=4, sigma=.2):
+
+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 = data["Y"]
     Y -= Y.mean()
     Y /= Y.std()
 
-    t = data['t']
+    t = data["t"]
-    c = data['colors']
+    c = data["colors"]
 
     try:
         from sklearn.manifold.isomap import Isomap
+
         iso = Isomap().fit(Y)
         X = iso.embedding_
         if Q > 2:
@@ -127,8 +143,9 @@ def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=25, Q=4
     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 = fig.add_subplot(121, projection="3d")
         ax.scatter(*Y.T, c=c)
         ax.set_title("Swiss Roll")
 
@@ -136,60 +153,96 @@ def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=25, Q=4
         ax.scatter(*X.T[:2], c=c)
         ax.set_title("BGPLVM init")
 
-    var = .5
+    var = 0.5
-    S = (var * _np.ones_like(X) + _np.clip(_np.random.randn(N, Q) * var ** 2,
+    S = (
-                                         - (1 - var),
+        var * _np.ones_like(X)
-                                         (1 - var))) + .001
+        + _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))
+    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 = 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)
+        m.optimize("bfgs", messages=verbose, max_iters=2e3)
 
     if plot:
-        fig = plt.figure('fitted')
+        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
 
-def bgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40, max_iters=1000, **k):
+
+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(
-    kernel = GPy.kern.RBF(Q, 1., 1. / _np.random.uniform(0, 1, (Q,)), ARD=True)  # + GPy.kern.Bias(Q, _np.exp(-2))
+        Q, 1.0, 1.0 / _np.random.uniform(0, 1, (Q,)), ARD=True
-    Y = data['X'][:N]
+    )  # + 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)
+    m.data_labels = data["Y"][:N].argmax(axis=1)
 
     if optimize:
-        m.optimize('bfgs', messages=verbose, max_iters=max_iters, gtol=.05)
+        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
+        lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(
-            m, data_show, latent_axes=latent_axes, sense_axes=sense_axes, labels=m.data_labels)
+            m.X.mean.values[0:1, :],  # @UnusedVariable
-        input('Press enter to finish')
+            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 ssgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40, max_iters=1000, **k):
+
+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
@@ -197,39 +250,57 @@ def ssgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40
     _np.random.seed(0)
     data = pods.datasets.oil()
 
-    kernel = GPy.kern.RBF(Q, 1., 1. / _np.random.uniform(0, 1, (Q,)), ARD=True)  # + GPy.kern.Bias(Q, _np.exp(-2))
+    kernel = GPy.kern.RBF(
-    Y = data['X'][:N]
+        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)
+    m.data_labels = data["Y"][:N].argmax(axis=1)
 
     if optimize:
-        m.optimize('bfgs', messages=verbose, max_iters=max_iters, gtol=.05)
+        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
+        lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(
-            m, data_show, latent_axes=latent_axes, sense_axes=sense_axes, labels=m.data_labels)
+            m.X.mean.values[0:1, :],  # @UnusedVariable
-        input('Press enter to finish')
+            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., lengthscale=(np.random.uniform(1, 6, Q_signal)), ARD=1)
+    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., active_dims=[i], period=3., lower=-2, upper=6)
+        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]
+    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)
+    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"]
@@ -239,6 +310,7 @@ def _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim=False):
         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)
@@ -248,13 +320,14 @@ def _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim=False):
         ax.legend()
         for i, Y in enumerate(Ylist):
             ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i)
-            ax.imshow(Y, aspect='auto', cmap=cm.gray)  # @UndefinedVariable
+            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)
@@ -262,7 +335,7 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False):
     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)))
+    s3 = _np.vectorize(lambda x: -_np.exp(-_np.cos(2 * x)))
     sS = _np.vectorize(lambda x: _np.cos(x))
 
     s1 = s1(x)
@@ -270,12 +343,18 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False):
     s3 = s3(x)
     sS = sS(x)
 
-    s1 -= s1.mean(); s1 /= s1.std(0)
+    s1 -= s1.mean()
-    s2 -= s2.mean(); s2 /= s2.std(0)
+    s1 /= s1.std(0)
-    s3 -= s3.mean(); s3 /= s3.std(0)
+    s2 -= s2.mean()
-    sS -= sS.mean(); sS /= sS.std(0)
+    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)
+    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"]
@@ -285,6 +364,7 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False):
         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)
@@ -294,13 +374,14 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, plot_sim=False):
         ax.legend()
         for i, Y in enumerate(Ylist):
             ax = fig.add_subplot(2, len(Ylist), len(Ylist) + 1 + i)
-            ax.imshow(Y, aspect='auto', cmap=cm.gray)  # @UndefinedVariable
+            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])
@@ -308,9 +389,9 @@ def _generate_high_dimensional_output(D1, D2, D3, s1, s2, 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 += .3 * _np.random.randn(*Y1.shape)
+    Y1 += 0.3 * _np.random.randn(*Y1.shape)
-    Y2 += .2 * _np.random.randn(*Y2.shape)
+    Y2 += 0.2 * _np.random.randn(*Y2.shape)
-    Y3 += .25 * _np.random.randn(*Y3.shape)
+    Y3 += 0.25 * _np.random.randn(*Y3.shape)
     Y1 -= Y1.mean(0)
     Y2 -= Y2.mean(0)
     Y3 -= Y3.mean(0)
@@ -319,10 +400,10 @@ def _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS):
     Y3 /= Y3.std(0)
     return Y1, Y2, Y3, S1, S2, S3
 
-def bgplvm_simulation(optimize=True, verbose=1,
+
-                      plot=True, plot_sim=False,
+def bgplvm_simulation(
-                      max_iters=2e4,
+    optimize=True, verbose=1, plot=True, plot_sim=False, max_iters=2e4,
-                      ):
+):
     from GPy import kern
     from GPy.models import BayesianGPLVM
 
@@ -332,22 +413,21 @@ def bgplvm_simulation(optimize=True, verbose=1,
     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, .01, m.X.shape)
+    m.X.variance[:] = _np.random.uniform(0, 0.01, m.X.shape)
-    m.likelihood.variance = .1
+    m.likelihood.variance = 0.1
 
     if optimize:
         print("Optimizing model:")
-        m.optimize('bfgs', messages=verbose, max_iters=max_iters,
+        m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05)
-                   gtol=.05)
     if plot:
         m.X.plot("BGPLVM Latent Space 1D")
         m.kern.plot_ARD()
     return m
 
-def gplvm_simulation(optimize=True, verbose=1,
+
-                      plot=True, plot_sim=False,
+def gplvm_simulation(
-                      max_iters=2e4,
+    optimize=True, verbose=1, plot=True, plot_sim=False, max_iters=2e4,
-                      ):
+):
     from GPy import kern
     from GPy.models import GPLVM
 
@@ -357,20 +437,20 @@ def gplvm_simulation(optimize=True, verbose=1,
     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 = .1
+    m.likelihood.variance = 0.1
 
     if optimize:
         print("Optimizing model:")
-        m.optimize('bfgs', messages=verbose, max_iters=max_iters,
+        m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05)
-                   gtol=.05)
     if plot:
         m.X.plot("BGPLVM Latent Space 1D")
         m.kern.plot_ARD()
     return m
-def ssgplvm_simulation(optimize=True, verbose=1,
+
-                      plot=True, plot_sim=False,
+
-                      max_iters=2e4, useGPU=False
+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
 
@@ -379,23 +459,30 @@ def ssgplvm_simulation(optimize=True, verbose=1,
     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 = SSGPLVM(
-    m.X.variance[:] = _np.random.uniform(0, .01, m.X.shape)
+        Y, Q, init="rand", num_inducing=num_inducing, kernel=k, group_spike=True
-    m.likelihood.variance = .01
+    )
+    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,
+        m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05)
-                   gtol=.05)
     if plot:
         m.X.plot("SSGPLVM Latent Space 1D")
         m.kern.plot_ARD()
     return m
 
-def bgplvm_simulation_missing_data(optimize=True, verbose=1,
+
-                      plot=True, plot_sim=False,
+def bgplvm_simulation_missing_data(
-                      max_iters=2e4, percent_missing=.1, d=13,
+    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
 
@@ -404,28 +491,42 @@ def bgplvm_simulation_missing_data(optimize=True, verbose=1,
     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
+    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,
+    m = BayesianGPLVMMiniBatch(
-                      kernel=k, missing_data=True)
+        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,
+        m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05)
-                   gtol=.05)
     if plot:
         m.X.plot("BGPLVM Latent Space 1D")
         m.kern.plot_ARD()
     return m
 
-def bgplvm_simulation_missing_data_stochastics(optimize=True, verbose=1,
+
-                      plot=True, plot_sim=False,
+def bgplvm_simulation_missing_data_stochastics(
-                      max_iters=2e4, percent_missing=.1, d=13, batchsize=2,
+    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
 
@@ -434,19 +535,28 @@ def bgplvm_simulation_missing_data_stochastics(optimize=True, verbose=1,
     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
+    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,
+    m = BayesianGPLVMMiniBatch(
-                      kernel=k, missing_data=True, stochastic=True, batchsize=batchsize)
+        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,
+        m.optimize("bfgs", messages=verbose, max_iters=max_iters, gtol=0.05)
-                   gtol=.05)
     if plot:
         m.X.plot("BGPLVM Latent Space 1D")
         m.kern.plot_ARD()
@@ -461,9 +571,17 @@ def mrd_simulation(optimize=True, verbose=True, plot=True, plot_sim=True, **kw):
     _, _, 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 = MRD(
+        Ylist,
+        input_dim=Q,
+        num_inducing=num_inducing,
+        kernel=k,
+        initx="PCA_concat",
+        initz="permute",
+        **kw
+    )
 
-    m['.*noise'] = [Y.var() / 40. for Y in Ylist]
+    m[".*noise"] = [Y.var() / 40.0 for Y in Ylist]
 
     if optimize:
         print("Optimizing Model:")
@@ -473,7 +591,10 @@ def mrd_simulation(optimize=True, verbose=True, plot=True, plot_sim=True, **kw):
         m.plot_scales()
     return m
 
-def mrd_simulation_missing_data(optimize=True, verbose=True, plot=True, plot_sim=True, **kw):
+
+def mrd_simulation_missing_data(
+    optimize=True, verbose=True, plot=True, plot_sim=True, **kw
+):
     from GPy import kern
     from GPy.models import MRD
 
@@ -484,29 +605,37 @@ def mrd_simulation_missing_data(optimize=True, verbose=True, plot=True, plot_sim
     inanlist = []
 
     for Y in Ylist:
-        inan = _np.random.binomial(1, .6, size=Y.shape).astype(bool)
+        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,
+    m = MRD(
-            kernel=k, inference_method=None,
+        Ylist,
-            initx="random", initz='permute', **kw)
+        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=.1)
+        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
 
+
 def brendan_faces(optimize=True, verbose=True, plot=True):
     import GPy
     import pods
 
     data = pods.datasets.brendan_faces()
     Q = 2
-    Y = data['Y']
+    Y = data["Y"]
     Yn = Y - Y.mean()
     Yn /= Yn.std()
 
@@ -514,39 +643,56 @@ def brendan_faces(optimize=True, verbose=True, plot=True):
 
     # optimize
 
-    if optimize: m.optimize('bfgs', messages=verbose, max_iters=1000)
+    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)
+        data_show = GPy.plotting.matplot_dep.visualize.image_show(
-        lvm = GPy.plotting.matplot_dep.visualize.lvm(m.X.mean[0, :].copy(), m, data_show, ax)
+            y[None, :],
-        input('Press enter to finish')
+            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
 
+
 def olivetti_faces(optimize=True, verbose=True, plot=True):
     import GPy
     import pods
 
     data = pods.datasets.olivetti_faces()
     Q = 2
-    Y = data['Y']
+    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 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)
+        data_show = GPy.plotting.matplot_dep.visualize.image_show(
-        lvm = GPy.plotting.matplot_dep.visualize.lvm(m.X.mean[0, :].copy(), m, data_show, ax)
+            y[None, :], dimensions=(112, 92), transpose=False, invert=False, scale=False
-        input('Press enter to finish')
+        )
+        lvm = GPy.plotting.matplot_dep.visualize.lvm(
+            m.X.mean[0, :].copy(), m, data_show, ax
+        )
+        input("Press enter to finish")
 
     return m
 
+
 def stick_play(range=None, frame_rate=15, optimize=False, verbose=True, plot=True):
     import GPy
     import pods
@@ -554,15 +700,18 @@ def stick_play(range=None, frame_rate=15, optimize=False, verbose=True, plot=Tru
     data = pods.datasets.osu_run1()
     # optimize
     if range == None:
-        Y = data['Y'].copy()
+        Y = data["Y"].copy()
     else:
-        Y = data['Y'][range[0]:range[1], :].copy()
+        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'])
+        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
 
+
 def stick(kernel=None, optimize=True, verbose=True, plot=True):
     from matplotlib import pyplot as plt
     import GPy
@@ -570,19 +719,25 @@ def stick(kernel=None, optimize=True, verbose=True, plot=True):
 
     data = pods.datasets.osu_run1()
     # optimize
-    m = GPy.models.GPLVM(data['Y'], 2, kernel=kernel)
+    m = GPy.models.GPLVM(data["Y"], 2, kernel=kernel)
-    if optimize: m.optimize('bfgs', messages=verbose, max_f_eval=10000)
+    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'])
+        data_show = GPy.plotting.matplot_dep.visualize.stick_show(
-        lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm(m.X[:1, :].copy(), m, data_show, latent_axes=ax)
+            y[None, :], connect=data["connect"]
-        input('Press enter to finish')
+        )
+        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
 
+
 def bcgplvm_linear_stick(kernel=None, optimize=True, verbose=True, plot=True):
     from matplotlib import pyplot as plt
     import GPy
@@ -590,19 +745,23 @@ def bcgplvm_linear_stick(kernel=None, optimize=True, verbose=True, plot=True):
 
     data = pods.datasets.osu_run1()
     # optimize
-    mapping = GPy.mappings.Linear(data['Y'].shape[1], 2)
+    mapping = GPy.mappings.Linear(data["Y"].shape[1], 2)
-    m = GPy.models.BCGPLVM(data['Y'], 2, kernel=kernel, mapping=mapping)
+    m = GPy.models.BCGPLVM(data["Y"], 2, kernel=kernel, mapping=mapping)
-    if optimize: m.optimize(messages=verbose, max_f_eval=10000)
+    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'])
+        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')
+        input("Press enter to finish")
 
     return m
 
+
 def bcgplvm_stick(kernel=None, optimize=True, verbose=True, plot=True):
     from matplotlib import pyplot as plt
     import GPy
@@ -610,20 +769,24 @@ def bcgplvm_stick(kernel=None, optimize=True, verbose=True, plot=True):
 
     data = pods.datasets.osu_run1()
     # optimize
-    back_kernel = GPy.kern.RBF(data['Y'].shape[1], lengthscale=5.)
+    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)
+    mapping = GPy.mappings.Kernel(X=data["Y"], output_dim=2, kernel=back_kernel)
-    m = GPy.models.BCGPLVM(data['Y'], 2, kernel=kernel, mapping=mapping)
+    m = GPy.models.BCGPLVM(data["Y"], 2, kernel=kernel, mapping=mapping)
-    if optimize: m.optimize(messages=verbose, max_f_eval=10000)
+    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'])
+        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
 
+
 def robot_wireless(optimize=True, verbose=True, plot=True):
     from matplotlib import pyplot as plt
     import GPy
@@ -631,13 +794,15 @@ def robot_wireless(optimize=True, verbose=True, plot=True):
 
     data = pods.datasets.robot_wireless()
     # optimize
-    m = GPy.models.BayesianGPLVM(data['Y'], 4, num_inducing=25)
+    m = GPy.models.BayesianGPLVM(data["Y"], 4, num_inducing=25)
-    if optimize: m.optimize(messages=verbose, max_f_eval=10000)
+    if optimize:
+        m.optimize(messages=verbose, max_f_eval=10000)
     if plot:
         m.plot_latent()
 
     return m
 
+
 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
@@ -648,15 +813,16 @@ def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True):
 
     data = pods.datasets.osu_run1()
     Q = 6
-    kernel = GPy.kern.RBF(Q, lengthscale=np.repeat(.5, Q), ARD=True)
+    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 = 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)
+        if optimize:
+            m.optimize("bfgs", messages=verbose, max_iters=5e3, bfgs_factor=10)
     except KeyboardInterrupt:
         print("Keyboard interrupt, continuing to plot and return")
 
@@ -665,17 +831,27 @@ def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True):
         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'])
+        data_show = GPy.plotting.matplot_dep.visualize.stick_show(
-        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)
+            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()
+        # fig.canvas.show()
-        input('Press enter to finish')
+        input("Press enter to finish")
 
     return m
 
 
-def cmu_mocap(subject='35', motion=['01'], in_place=True, optimize=True, verbose=True, plot=True):
+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
@@ -683,34 +859,40 @@ def cmu_mocap(subject='35', motion=['01'], in_place=True, optimize=True, verbose
     data = pods.datasets.cmu_mocap(subject, motion)
     if in_place:
         # Make figure move in place.
-        data['Y'][:, 0:3] = 0.0
+        data["Y"][:, 0:3] = 0.0
-    Y = data['Y']
+    Y = data["Y"]
     m = GPy.models.GPLVM(Y, 2, normalizer=True)
 
-    if optimize: m.optimize(messages=verbose, max_f_eval=10000)
+    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')
+        viz_axes = fig.add_subplot(133, projection="3d")
-
 
         m.plot_latent(ax=latent_axes)
-        latent_axes.set_aspect('equal')
+        latent_axes.set_aspect("equal")
 
         y = m.Y[0, :]
-        data_show = GPy.plotting.matplot_dep.visualize.skeleton_show(y[None, :], data['skel'], viz_axes)
+        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)        
+        lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm(
-        input('Press enter to finish')
+            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
 
+
 def ssgplvm_simulation_linear():
     import numpy as np
     import GPy
+
     N, D, Q = 1000, 20, 5
     pi = 0.2
 
@@ -721,7 +903,7 @@ def ssgplvm_simulation_linear():
             if dies[q] < pi:
                 x[q] = np.random.randn()
             else:
-                x[q] = 0.
+                x[q] = 0.0
         return x
 
     Y = np.empty((N, D))
@@ -731,4 +913,3 @@ def ssgplvm_simulation_linear():
         X[n] = sample_X(Q, pi)
         w = np.random.randn(D, Q)
         Y[n] = np.dot(w, X[n])
-

GPy/examples/non_gaussian.py

@@ -4,38 +4,40 @@
 import GPy
 import numpy as np
 from GPy.util import datasets
+
 try:
     import matplotlib.pyplot as plt
 except:
     pass
 
+
 def student_t_approx(optimize=True, plot=True):
     """
     Example of regressing with a student t likelihood using Laplace
     """
     real_std = 0.1
-    #Start a function, any function
+    # Start a function, any function
-    X = np.linspace(0.0, np.pi*2, 100)[:, None]
+    X = np.linspace(0.0, np.pi * 2, 100)[:, None]
-    Y = np.sin(X) + np.random.randn(*X.shape)*real_std
+    Y = np.sin(X) + np.random.randn(*X.shape) * real_std
-    Y = Y/Y.max()
+    Y = Y / Y.max()
     Yc = Y.copy()
 
-    X_full = np.linspace(0.0, np.pi*2, 500)[:, None]
+    X_full = np.linspace(0.0, np.pi * 2, 500)[:, None]
     Y_full = np.sin(X_full)
-    Y_full = Y_full/Y_full.max()
+    Y_full = Y_full / Y_full.max()
 
-    #Slightly noisy data
+    # Slightly noisy data
     Yc[75:80] += 1
 
-    #Very noisy data
+    # Very noisy data
-    #Yc[10] += 100
+    # Yc[10] += 100
-    #Yc[25] += 10
+    # Yc[25] += 10
-    #Yc[23] += 10
+    # Yc[23] += 10
-    #Yc[26] += 1000
+    # Yc[26] += 1000
-    #Yc[24] += 10
+    # Yc[24] += 10
-    #Yc = Yc/Yc.max()
+    # Yc = Yc/Yc.max()
 
-    #Add student t random noise to datapoints
+    # Add student t random noise to datapoints
     deg_free = 1
     print("Real noise: ", real_std)
     initial_var_guess = 0.5
@@ -47,45 +49,50 @@ def student_t_approx(optimize=True, plot=True):
     kernel3 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
     kernel4 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
 
-    #Gaussian GP model on clean data
+    # Gaussian GP model on clean data
     m1 = GPy.models.GPRegression(X, Y.copy(), kernel=kernel1)
     # optimize
-    m1['.*white'].constrain_fixed(1e-5)
+    m1[".*white"].constrain_fixed(1e-5)
     m1.randomize()
 
-    #Gaussian GP model on corrupt data
+    # Gaussian GP model on corrupt data
     m2 = GPy.models.GPRegression(X, Yc.copy(), kernel=kernel2)
-    m2['.*white'].constrain_fixed(1e-5)
+    m2[".*white"].constrain_fixed(1e-5)
     m2.randomize()
 
-    #Student t GP model on clean data
+    # Student t GP model on clean data
     t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd)
     laplace_inf = GPy.inference.latent_function_inference.Laplace()
-    m3 = GPy.core.GP(X, Y.copy(), kernel3, likelihood=t_distribution, inference_method=laplace_inf)
+    m3 = GPy.core.GP(
-    m3['.*t_scale2'].constrain_bounded(1e-6, 10.)
+        X, Y.copy(), kernel3, likelihood=t_distribution, inference_method=laplace_inf
-    m3['.*white'].constrain_fixed(1e-5)
+    )
+    m3[".*t_scale2"].constrain_bounded(1e-6, 10.0)
+    m3[".*white"].constrain_fixed(1e-5)
     m3.randomize()
 
-    #Student t GP model on corrupt data
+    # Student t GP model on corrupt data
     t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd)
     laplace_inf = GPy.inference.latent_function_inference.Laplace()
-    m4 = GPy.core.GP(X, Yc.copy(), kernel4, likelihood=t_distribution, inference_method=laplace_inf)
+    m4 = GPy.core.GP(
-    m4['.*t_scale2'].constrain_bounded(1e-6, 10.)
+        X, Yc.copy(), kernel4, likelihood=t_distribution, inference_method=laplace_inf
-    m4['.*white'].constrain_fixed(1e-5)
+    )
+    m4[".*t_scale2"].constrain_bounded(1e-6, 10.0)
+    m4[".*white"].constrain_fixed(1e-5)
     m4.randomize()
     print(m4)
-    debug=True
+    debug = True
     if debug:
         m4.optimize(messages=1)
         from matplotlib import pyplot as pb
+
         pb.plot(m4.X, m4.inference_method.f_hat)
-        pb.plot(m4.X, m4.Y, 'rx')
+        pb.plot(m4.X, m4.Y, "rx")
         m4.plot()
         print(m4)
         return m4
 
     if optimize:
-        optimizer='scg'
+        optimizer = "scg"
         print("Clean Gaussian")
         m1.optimize(optimizer, messages=1)
         print("Corrupt Gaussian")
@@ -97,77 +104,91 @@ def student_t_approx(optimize=True, plot=True):
 
     if plot:
         plt.figure(1)
-        plt.suptitle('Gaussian likelihood')
+        plt.suptitle("Gaussian likelihood")
         ax = plt.subplot(211)
         m1.plot(ax=ax)
         plt.plot(X_full, Y_full)
         plt.ylim(-1.5, 1.5)
-        plt.title('Gaussian clean')
+        plt.title("Gaussian clean")
 
         ax = plt.subplot(212)
         m2.plot(ax=ax)
         plt.plot(X_full, Y_full)
         plt.ylim(-1.5, 1.5)
-        plt.title('Gaussian corrupt')
+        plt.title("Gaussian corrupt")
 
         plt.figure(2)
-        plt.suptitle('Student-t likelihood')
+        plt.suptitle("Student-t likelihood")
         ax = plt.subplot(211)
         m3.plot(ax=ax)
         plt.plot(X_full, Y_full)
         plt.ylim(-1.5, 1.5)
-        plt.title('Student-t rasm clean')
+        plt.title("Student-t rasm clean")
 
         ax = plt.subplot(212)
         m4.plot(ax=ax)
         plt.plot(X_full, Y_full)
         plt.ylim(-1.5, 1.5)
-        plt.title('Student-t rasm corrupt')
+        plt.title("Student-t rasm corrupt")
 
     return m1, m2, m3, m4
 
+
 def boston_example(optimize=True, plot=True):
     raise NotImplementedError("Needs updating")
     import sklearn
     from sklearn.cross_validation import KFold
-    optimizer='bfgs'
+
-    messages=0
+    optimizer = "bfgs"
+    messages = 0
     data = datasets.boston_housing()
     degrees_freedoms = [3, 5, 8, 10]
-    X = data['X'].copy()
+    X = data["X"].copy()
-    Y = data['Y'].copy()
+    Y = data["Y"].copy()
-    X = X-X.mean(axis=0)
+    X = X - X.mean(axis=0)
-    X = X/X.std(axis=0)
+    X = X / X.std(axis=0)
-    Y = Y-Y.mean()
+    Y = Y - Y.mean()
-    Y = Y/Y.std()
+    Y = Y / Y.std()
     num_folds = 10
     kf = KFold(len(Y), n_folds=num_folds, indices=True)
-    num_models = len(degrees_freedoms) + 3 #3 for baseline, gaussian, gaussian laplace approx
+    num_models = (
+        len(degrees_freedoms) + 3
+    )  # 3 for baseline, gaussian, gaussian laplace approx
     score_folds = np.zeros((num_models, num_folds))
     pred_density = score_folds.copy()
 
     def rmse(Y, Ystar):
-        return np.sqrt(np.mean((Y-Ystar)**2))
+        return np.sqrt(np.mean((Y - Ystar) ** 2))
 
     for n, (train, test) in enumerate(kf):
         X_train, X_test, Y_train, Y_test = X[train], X[test], Y[train], Y[test]
         print("Fold {}".format(n))
 
-        noise = 1e-1 #np.exp(-2)
+        noise = 1e-1  # np.exp(-2)
         rbf_len = 0.5
         data_axis_plot = 4
-        kernelstu = GPy.kern.RBF(X.shape[1]) + GPy.kern.white(X.shape[1]) + GPy.kern.bias(X.shape[1])
+        kernelstu = (
-        kernelgp = GPy.kern.RBF(X.shape[1]) + GPy.kern.white(X.shape[1]) + GPy.kern.bias(X.shape[1])
+            GPy.kern.RBF(X.shape[1])
+            + GPy.kern.white(X.shape[1])
+            + GPy.kern.bias(X.shape[1])
+        )
+        kernelgp = (
+            GPy.kern.RBF(X.shape[1])
+            + GPy.kern.white(X.shape[1])
+            + GPy.kern.bias(X.shape[1])
+        )
 
-        #Baseline
+        # Baseline
         score_folds[0, n] = rmse(Y_test, np.mean(Y_train))
 
-        #Gaussian GP
+        # Gaussian GP
         print("Gauss GP")
-        mgp = GPy.models.GPRegression(X_train.copy(), Y_train.copy(), kernel=kernelgp.copy())
+        mgp = GPy.models.GPRegression(
-        mgp.constrain_fixed('.*white', 1e-5)
+            X_train.copy(), Y_train.copy(), kernel=kernelgp.copy()
-        mgp['.*len'] = rbf_len
+        )
-        mgp['.*noise'] = noise
+        mgp.constrain_fixed(".*white", 1e-5)
+        mgp[".*len"] = rbf_len
+        mgp[".*noise"] = noise
         print(mgp)
         if optimize:
             mgp.optimize(optimizer=optimizer, messages=messages)
@@ -179,13 +200,20 @@ def boston_example(optimize=True, plot=True):
 
         print("Gaussian Laplace GP")
         N, D = Y_train.shape
-        g_distribution = GPy.likelihoods.noise_model_constructors.gaussian(variance=noise, N=N, D=D)
+        g_distribution = GPy.likelihoods.noise_model_constructors.gaussian(
+            variance=noise, N=N, D=D
+        )
         g_likelihood = GPy.likelihoods.Laplace(Y_train.copy(), g_distribution)
-        mg = GPy.models.GPRegression(X_train.copy(), Y_train.copy(), kernel=kernelstu.copy(), likelihood=g_likelihood)
+        mg = GPy.models.GPRegression(
-        mg.constrain_positive('noise_variance')
+            X_train.copy(),
-        mg.constrain_fixed('.*white', 1e-5)
+            Y_train.copy(),
-        mg['rbf_len'] = rbf_len
+            kernel=kernelstu.copy(),
-        mg['noise'] = noise
+            likelihood=g_likelihood,
+        )
+        mg.constrain_positive("noise_variance")
+        mg.constrain_fixed(".*white", 1e-5)
+        mg["rbf_len"] = rbf_len
+        mg["noise"] = noise
         print(mg)
         if optimize:
             mg.optimize(optimizer=optimizer, messages=messages)
@@ -196,101 +224,108 @@ def boston_example(optimize=True, plot=True):
         print(mg)
 
         for stu_num, df in enumerate(degrees_freedoms):
-            #Student T
+            # Student T
             print("Student-T GP {}df".format(df))
-            t_distribution = GPy.likelihoods.noise_model_constructors.student_t(deg_free=df, sigma2=noise)
+            t_distribution = GPy.likelihoods.noise_model_constructors.student_t(
+                deg_free=df, sigma2=noise
+            )
             stu_t_likelihood = GPy.likelihoods.Laplace(Y_train.copy(), t_distribution)
-            mstu_t = GPy.models.GPRegression(X_train.copy(), Y_train.copy(), kernel=kernelstu.copy(), likelihood=stu_t_likelihood)
+            mstu_t = GPy.models.GPRegression(
-            mstu_t.constrain_fixed('.*white', 1e-5)
+                X_train.copy(),
-            mstu_t.constrain_bounded('.*t_scale2', 0.0001, 1000)
+                Y_train.copy(),
-            mstu_t['rbf_len'] = rbf_len
+                kernel=kernelstu.copy(),
-            mstu_t['.*t_scale2'] = noise
+                likelihood=stu_t_likelihood,
+            )
+            mstu_t.constrain_fixed(".*white", 1e-5)
+            mstu_t.constrain_bounded(".*t_scale2", 0.0001, 1000)
+            mstu_t["rbf_len"] = rbf_len
+            mstu_t[".*t_scale2"] = noise
             print(mstu_t)
             if optimize:
                 mstu_t.optimize(optimizer=optimizer, messages=messages)
             Y_test_pred = mstu_t.predict(X_test)
-            score_folds[3+stu_num, n] = rmse(Y_test, Y_test_pred[0])
+            score_folds[3 + stu_num, n] = rmse(Y_test, Y_test_pred[0])
-            pred_density[3+stu_num, n] = np.mean(mstu_t.log_predictive_density(X_test, Y_test))
+            pred_density[3 + stu_num, n] = np.mean(
+                mstu_t.log_predictive_density(X_test, Y_test)
+            )
             print(pred_density)
             print(mstu_t)
 
     if plot:
         plt.figure()
         plt.scatter(X_test[:, data_axis_plot], Y_test_pred[0])
-        plt.scatter(X_test[:, data_axis_plot], Y_test, c='r', marker='x')
+        plt.scatter(X_test[:, data_axis_plot], Y_test, c="r", marker="x")
-        plt.title('GP gauss')
+        plt.title("GP gauss")
 
         plt.figure()
         plt.scatter(X_test[:, data_axis_plot], Y_test_pred[0])
-        plt.scatter(X_test[:, data_axis_plot], Y_test, c='r', marker='x')
+        plt.scatter(X_test[:, data_axis_plot], Y_test, c="r", marker="x")
-        plt.title('Lap gauss')
+        plt.title("Lap gauss")
 
         plt.figure()
         plt.scatter(X_test[:, data_axis_plot], Y_test_pred[0])
-        plt.scatter(X_test[:, data_axis_plot], Y_test, c='r', marker='x')
+        plt.scatter(X_test[:, data_axis_plot], Y_test, c="r", marker="x")
-        plt.title('Stu t {}df'.format(df))
+        plt.title("Stu t {}df".format(df))
 
     print("Average scores: {}".format(np.mean(score_folds, 1)))
     print("Average pred density: {}".format(np.mean(pred_density, 1)))
 
     if plot:
-        #Plotting
+        # Plotting
-        stu_t_legends = ['Student T, df={}'.format(df) for df in degrees_freedoms]
+        stu_t_legends = ["Student T, df={}".format(df) for df in degrees_freedoms]
-        legends = ['Baseline', 'Gaussian', 'Laplace Approx Gaussian'] + stu_t_legends
+        legends = ["Baseline", "Gaussian", "Laplace Approx Gaussian"] + stu_t_legends
 
-        #Plot boxplots for RMSE density
+        # Plot boxplots for RMSE density
         fig = plt.figure()
-        ax=fig.add_subplot(111)
+        ax = fig.add_subplot(111)
-        plt.title('RMSE')
+        plt.title("RMSE")
-        bp = ax.boxplot(score_folds.T, notch=0, sym='+', vert=1, whis=1.5)
+        bp = ax.boxplot(score_folds.T, notch=0, sym="+", vert=1, whis=1.5)
-        plt.setp(bp['boxes'], color='black')
+        plt.setp(bp["boxes"], color="black")
-        plt.setp(bp['whiskers'], color='black')
+        plt.setp(bp["whiskers"], color="black")
-        plt.setp(bp['fliers'], color='red', marker='+')
+        plt.setp(bp["fliers"], color="red", marker="+")
         xtickNames = plt.setp(ax, xticklabels=legends)
         plt.setp(xtickNames, rotation=45, fontsize=8)
-        ax.set_ylabel('RMSE')
+        ax.set_ylabel("RMSE")
-        ax.set_xlabel('Distribution')
+        ax.set_xlabel("Distribution")
-        #Make grid and put it below boxes
+        # Make grid and put it below boxes
-        ax.yaxis.grid(True, linestyle='-', which='major', color='lightgrey',
+        ax.yaxis.grid(True, linestyle="-", which="major", color="lightgrey", alpha=0.5)
-                alpha=0.5)
         ax.set_axisbelow(True)
 
-        #Plot boxplots for predictive density
+        # Plot boxplots for predictive density
         fig = plt.figure()
-        ax=fig.add_subplot(111)
+        ax = fig.add_subplot(111)
-        plt.title('Predictive density')
+        plt.title("Predictive density")
-        bp = ax.boxplot(pred_density[1:,:].T, notch=0, sym='+', vert=1, whis=1.5)
+        bp = ax.boxplot(pred_density[1:, :].T, notch=0, sym="+", vert=1, whis=1.5)
-        plt.setp(bp['boxes'], color='black')
+        plt.setp(bp["boxes"], color="black")
-        plt.setp(bp['whiskers'], color='black')
+        plt.setp(bp["whiskers"], color="black")
-        plt.setp(bp['fliers'], color='red', marker='+')
+        plt.setp(bp["fliers"], color="red", marker="+")
         xtickNames = plt.setp(ax, xticklabels=legends[1:])
         plt.setp(xtickNames, rotation=45, fontsize=8)
-        ax.set_ylabel('Mean Log probability P(Y*|Y)')
+        ax.set_ylabel("Mean Log probability P(Y*|Y)")
-        ax.set_xlabel('Distribution')
+        ax.set_xlabel("Distribution")
-        #Make grid and put it below boxes
+        # Make grid and put it below boxes
-        ax.yaxis.grid(True, linestyle='-', which='major', color='lightgrey',
+        ax.yaxis.grid(True, linestyle="-", which="major", color="lightgrey", alpha=0.5)
-                alpha=0.5)
         ax.set_axisbelow(True)
     return mstu_t
 
-#def precipitation_example():
-    #import sklearn
-    #from sklearn.cross_validation import KFold
-    #data = datasets.boston_housing()
-    #X = data['X'].copy()
-    #Y = data['Y'].copy()
-    #X = X-X.mean(axis=0)
-    #X = X/X.std(axis=0)
-    #Y = Y-Y.mean()
-    #Y = Y/Y.std()
-    #import ipdb; ipdb.set_trace()  # XXX BREAKPOINT
-    #num_folds = 10
-    #kf = KFold(len(Y), n_folds=num_folds, indices=True)
-    #score_folds = np.zeros((4, num_folds))
-    #def rmse(Y, Ystar):
-        #return np.sqrt(np.mean((Y-Ystar)**2))
-    ##for train, test in kf:
-    #for n, (train, test) in enumerate(kf):
-        #X_train, X_test, Y_train, Y_test = X[train], X[test], Y[train], Y[test]
-        #print "Fold {}".format(n)
 
+# def precipitation_example():
+# import sklearn
+# from sklearn.cross_validation import KFold
+# data = datasets.boston_housing()
+# X = data['X'].copy()
+# Y = data['Y'].copy()
+# X = X-X.mean(axis=0)
+# X = X/X.std(axis=0)
+# Y = Y-Y.mean()
+# Y = Y/Y.std()
+# import ipdb; ipdb.set_trace()  # XXX BREAKPOINT
+# num_folds = 10
+# kf = KFold(len(Y), n_folds=num_folds, indices=True)
+# score_folds = np.zeros((4, num_folds))
+# def rmse(Y, Ystar):
+# return np.sqrt(np.mean((Y-Ystar)**2))
+##for train, test in kf:
+# for n, (train, test) in enumerate(kf):
+# X_train, X_test, Y_train, Y_test = X[train], X[test], Y[train], Y[test]
+# print "Fold {}".format(n)

GPy/examples/regression.py

@@ -11,88 +11,112 @@ except:
 import numpy as np
 import GPy
 
+
 def olympic_marathon_men(optimize=True, plot=True):
     """Run a standard Gaussian process regression on the Olympic marathon data."""
-    try:import pods
+    try:
+        import pods
     except ImportError:
-        print('pods unavailable, see https://github.com/sods/ods for example datasets')
+        print("pods unavailable, see https://github.com/sods/ods for example datasets")
         return
     data = pods.datasets.olympic_marathon_men()
 
     # create simple GP Model
-    m = GPy.models.GPRegression(data['X'], data['Y'])
+    m = GPy.models.GPRegression(data["X"], data["Y"])
 
     # set the lengthscale to be something sensible (defaults to 1)
-    m.kern.lengthscale = 10.
+    m.kern.lengthscale = 10.0
 
     if optimize:
-        m.optimize('bfgs', max_iters=200)
+        m.optimize("bfgs", max_iters=200)
     if plot:
         m.plot(plot_limits=(1850, 2050))
 
     return m
 
+
 def coregionalization_toy(optimize=True, plot=True):
     """
     A simple demonstration of coregionalization on two sinusoidal functions.
     """
-    #build a design matrix with a column of integers indicating the output
+    # build a design matrix with a column of integers indicating the output
     X1 = np.random.rand(50, 1) * 8
     X2 = np.random.rand(30, 1) * 5
 
-    #build a suitable set of observed variables
+    # build a suitable set of observed variables
     Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
-    Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2.
+    Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2.0
 
-    m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2])
+    m = GPy.models.GPCoregionalizedRegression(X_list=[X1, X2], Y_list=[Y1, Y2])
 
     if optimize:
-        m.optimize('bfgs', max_iters=100)
+        m.optimize("bfgs", max_iters=100)
 
     if plot:
-        slices = GPy.util.multioutput.get_slices([X1,X2])
+        slices = GPy.util.multioutput.get_slices([X1, X2])
-        m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0})
+        m.plot(
-        m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca())
+            fixed_inputs=[(1, 0)],
+            which_data_rows=slices[0],
+            Y_metadata={"output_index": 0},
+        )
+        m.plot(
+            fixed_inputs=[(1, 1)],
+            which_data_rows=slices[1],
+            Y_metadata={"output_index": 1},
+            ax=pb.gca(),
+        )
     return m
 
+
 def coregionalization_sparse(optimize=True, plot=True):
     """
     A simple demonstration of coregionalization on two sinusoidal functions using sparse approximations.
     """
-    #build a design matrix with a column of integers indicating the output
+    # build a design matrix with a column of integers indicating the output
     X1 = np.random.rand(50, 1) * 8
     X2 = np.random.rand(30, 1) * 5
 
-    #build a suitable set of observed variables
+    # build a suitable set of observed variables
     Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
-    Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2.
+    Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2.0
 
-    m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2])
+    m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1, X2], Y_list=[Y1, Y2])
 
     if optimize:
-        m.optimize('bfgs', max_iters=100)
+        m.optimize("bfgs", max_iters=100)
 
     if plot:
-        slices = GPy.util.multioutput.get_slices([X1,X2])
+        slices = GPy.util.multioutput.get_slices([X1, X2])
-        m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0})
+        m.plot(
-        m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca())
+            fixed_inputs=[(1, 0)],
+            which_data_rows=slices[0],
+            Y_metadata={"output_index": 0},
+        )
+        m.plot(
+            fixed_inputs=[(1, 1)],
+            which_data_rows=slices[1],
+            Y_metadata={"output_index": 1},
+            ax=pb.gca(),
+        )
         pb.ylim(-3,)
 
     return m
 
+
 def epomeo_gpx(max_iters=200, optimize=True, plot=True):
     """
     Perform Gaussian process regression on the latitude and longitude data
     from the Mount Epomeo runs. Requires gpxpy to be installed on your system
     to load in the data.
     """
-    try:import pods
+    try:
+        import pods
     except ImportError:
-        print('pods unavailable, see https://github.com/sods/ods for example datasets')
+        print("pods unavailable, see https://github.com/sods/ods for example datasets")
         return
     data = pods.datasets.epomeo_gpx()
     num_data_list = []
-    for Xpart in data['X']:
+    for Xpart in data["X"]:
         num_data_list.append(Xpart.shape[0])
 
     num_data_array = np.array(num_data_list)
@@ -100,29 +124,43 @@ def epomeo_gpx(max_iters=200, optimize=True, plot=True):
     Y = np.zeros((num_data, 2))
     t = np.zeros((num_data, 2))
     start = 0
-    for Xpart, index in zip(data['X'], range(len(data['X']))):
+    for Xpart, index in zip(data["X"], range(len(data["X"]))):
-        end = start+Xpart.shape[0]
+        end = start + Xpart.shape[0]
-        t[start:end, :] = np.hstack((Xpart[:, 0:1],
+        t[start:end, :] = np.hstack(
-                                    index*np.ones((Xpart.shape[0], 1))))
+            (Xpart[:, 0:1], index * np.ones((Xpart.shape[0], 1)))
+        )
         Y[start:end, :] = Xpart[:, 1:3]
 
     num_inducing = 200
-    Z = np.hstack((np.linspace(t[:,0].min(), t[:, 0].max(), num_inducing)[:, None],
+    Z = np.hstack(
-                   np.random.randint(0, 4, num_inducing)[:, None]))
+        (
+            np.linspace(t[:, 0].min(), t[:, 0].max(), num_inducing)[:, None],
+            np.random.randint(0, 4, num_inducing)[:, None],
+        )
+    )
 
     k1 = GPy.kern.RBF(1)
     k2 = GPy.kern.Coregionalize(output_dim=5, rank=5)
-    k = k1**k2
+    k = k1 ** k2
 
     m = GPy.models.SparseGPRegression(t, Y, kernel=k, Z=Z, normalize_Y=True)
-    m.constrain_fixed('.*variance', 1.)
+    m.constrain_fixed(".*variance", 1.0)
     m.inducing_inputs.constrain_fixed()
     m.Gaussian_noise.variance.constrain_bounded(1e-3, 1e-1)
-    m.optimize(max_iters=max_iters,messages=True)
+    m.optimize(max_iters=max_iters, messages=True)
 
     return m
 
-def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=10000, max_iters=300, optimize=True, plot=True):
+
+def multiple_optima(
+    gene_number=937,
+    resolution=80,
+    model_restarts=10,
+    seed=10000,
+    max_iters=300,
+    optimize=True,
+    plot=True,
+):
     """
     Show an example of a multimodal error surface for Gaussian process
     regression. Gene 939 has bimodal behaviour where the noisy mode is
@@ -130,25 +168,30 @@ def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=1000
     """
 
     # Contour over a range of length scales and signal/noise ratios.
-    length_scales = np.linspace(0.1, 60., resolution)
+    length_scales = np.linspace(0.1, 60.0, resolution)
-    log_SNRs = np.linspace(-3., 4., resolution)
+    log_SNRs = np.linspace(-3.0, 4.0, resolution)
 
-    try:import pods
+    try:
+        import pods
     except ImportError:
-        print('pods unavailable, see https://github.com/sods/ods for example datasets')
+        print("pods unavailable, see https://github.com/sods/ods for example datasets")
         return
-    data = pods.datasets.della_gatta_TRP63_gene_expression(data_set='della_gatta',gene_number=gene_number)
+    data = pods.datasets.della_gatta_TRP63_gene_expression(
+        data_set="della_gatta", gene_number=gene_number
+    )
     # data['Y'] = data['Y'][0::2, :]
     # data['X'] = data['X'][0::2, :]
 
-    data['Y'] = data['Y'] - np.mean(data['Y'])
+    data["Y"] = data["Y"] - np.mean(data["Y"])
 
-    lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.RBF)
+    lls = GPy.examples.regression._contour_data(
+        data, length_scales, log_SNRs, GPy.kern.RBF
+    )
     if plot:
         pb.contour(length_scales, log_SNRs, np.exp(lls), 20, cmap=pb.cm.jet)
         ax = pb.gca()
-        pb.xlabel('length scale')
+        pb.xlabel("length scale")
-        pb.ylabel('log_10 SNR')
+        pb.ylabel("log_10 SNR")
 
         xlim = ax.get_xlim()
         ylim = ax.get_ylim()
@@ -160,28 +203,41 @@ def multiple_optima(gene_number=937, resolution=80, model_restarts=10, seed=1000
     np.random.seed(seed=seed)
     for i in range(0, model_restarts):
         # kern = GPy.kern.RBF(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.))
-        kern = GPy.kern.RBF(1, variance=np.random.uniform(1e-3, 1), lengthscale=np.random.uniform(5, 50))
+        kern = GPy.kern.RBF(
+            1, variance=np.random.uniform(1e-3, 1), lengthscale=np.random.uniform(5, 50)
+        )
 
-        m = GPy.models.GPRegression(data['X'], data['Y'], kernel=kern)
+        m = GPy.models.GPRegression(data["X"], data["Y"], kernel=kern)
         m.likelihood.variance = np.random.uniform(1e-3, 1)
         optim_point_x[0] = m.rbf.lengthscale
-        optim_point_y[0] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);
+        optim_point_y[0] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance)
 
         # optimize
         if optimize:
-            m.optimize('scg', xtol=1e-6, ftol=1e-6, max_iters=max_iters)
+            m.optimize("scg", xtol=1e-6, ftol=1e-6, max_iters=max_iters)
 
         optim_point_x[1] = m.rbf.lengthscale
-        optim_point_y[1] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance);
+        optim_point_y[1] = np.log10(m.rbf.variance) - np.log10(m.likelihood.variance)
 
         if plot:
-            pb.arrow(optim_point_x[0], optim_point_y[0], optim_point_x[1] - optim_point_x[0], optim_point_y[1] - optim_point_y[0], label=str(i), head_length=1, head_width=0.5, fc='k', ec='k')
+            pb.arrow(
+                optim_point_x[0],
+                optim_point_y[0],
+                optim_point_x[1] - optim_point_x[0],
+                optim_point_y[1] - optim_point_y[0],
+                label=str(i),
+                head_length=1,
+                head_width=0.5,
+                fc="k",
+                ec="k",
+            )
         models.append(m)
 
     if plot:
         ax.set_xlim(xlim)
         ax.set_ylim(ylim)
-    return m # (models, lls)
+    return m  # (models, lls)
+
 
 def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.RBF):
     """
@@ -195,19 +251,19 @@ def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.RBF):
     """
 
     lls = []
-    total_var = np.var(data['Y'])
+    total_var = np.var(data["Y"])
-    kernel = kernel_call(1, variance=1., lengthscale=1.)
+    kernel = kernel_call(1, variance=1.0, lengthscale=1.0)
-    model = GPy.models.GPRegression(data['X'], data['Y'], kernel=kernel)
+    model = GPy.models.GPRegression(data["X"], data["Y"], kernel=kernel)
     for log_SNR in log_SNRs:
-        SNR = 10.**log_SNR
+        SNR = 10.0 ** log_SNR
-        noise_var = total_var / (1. + SNR)
+        noise_var = total_var / (1.0 + SNR)
         signal_var = total_var - noise_var
-        model.kern['.*variance'] = signal_var
+        model.kern[".*variance"] = signal_var
         model.likelihood.variance = noise_var
         length_scale_lls = []
 
         for length_scale in length_scales:
-            model['.*lengthscale'] = length_scale
+            model[".*lengthscale"] = length_scale
             length_scale_lls.append(model.log_likelihood())
 
         lls.append(length_scale_lls)
@@ -217,86 +273,97 @@ def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.RBF):
 
 def olympic_100m_men(optimize=True, plot=True):
     """Run a standard Gaussian process regression on the Rogers and Girolami olympics data."""
-    try:import pods
+    try:
+        import pods
     except ImportError:
-        print('pods unavailable, see https://github.com/sods/ods for example datasets')
+        print("pods unavailable, see https://github.com/sods/ods for example datasets")
         return
     data = pods.datasets.olympic_100m_men()
 
     # create simple GP Model
-    m = GPy.models.GPRegression(data['X'], data['Y'])
+    m = GPy.models.GPRegression(data["X"], data["Y"])
 
     # set the lengthscale to be something sensible (defaults to 1)
     m.rbf.lengthscale = 10
 
     if optimize:
-        m.optimize('bfgs', max_iters=200)
+        m.optimize("bfgs", max_iters=200)
 
     if plot:
         m.plot(plot_limits=(1850, 2050))
     return m
 
+
 def toy_rbf_1d(optimize=True, plot=True):
     """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
-    try:import pods
+    try:
+        import pods
     except ImportError:
-        print('pods unavailable, see https://github.com/sods/ods for example datasets')
+        print("pods unavailable, see https://github.com/sods/ods for example datasets")
         return
     data = pods.datasets.toy_rbf_1d()
 
     # create simple GP Model
-    m = GPy.models.GPRegression(data['X'], data['Y'])
+    m = GPy.models.GPRegression(data["X"], data["Y"])
 
     if optimize:
-        m.optimize('bfgs')
+        m.optimize("bfgs")
     if plot:
         m.plot()
 
     return m
 
+
 def toy_rbf_1d_50(optimize=True, plot=True):
     """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
-    try:import pods
+    try:
+        import pods
     except ImportError:
-        print('pods unavailable, see https://github.com/sods/ods for example datasets')
+        print("pods unavailable, see https://github.com/sods/ods for example datasets")
         return
     data = pods.datasets.toy_rbf_1d_50()
 
     # create simple GP Model
-    m = GPy.models.GPRegression(data['X'], data['Y'])
+    m = GPy.models.GPRegression(data["X"], data["Y"])
 
     if optimize:
-        m.optimize('bfgs')
+        m.optimize("bfgs")
     if plot:
         m.plot()
 
     return m
 
+
 def toy_poisson_rbf_1d_laplace(optimize=True, plot=True):
     """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
-    optimizer='scg'
+    optimizer = "scg"
     x_len = 100
     X = np.linspace(0, 10, x_len)[:, None]
     f_true = np.random.multivariate_normal(np.zeros(x_len), GPy.kern.RBF(1).K(X))
-    Y = np.array([np.random.poisson(np.exp(f)) for f in f_true])[:,None]
+    Y = np.array([np.random.poisson(np.exp(f)) for f in f_true])[:, None]
 
     kern = GPy.kern.RBF(1)
     poisson_lik = GPy.likelihoods.Poisson()
     laplace_inf = GPy.inference.latent_function_inference.Laplace()
 
     # create simple GP Model
-    m = GPy.core.GP(X, Y, kernel=kern, likelihood=poisson_lik, inference_method=laplace_inf)
+    m = GPy.core.GP(
+        X, Y, kernel=kern, likelihood=poisson_lik, inference_method=laplace_inf
+    )
 
     if optimize:
         m.optimize(optimizer)
     if plot:
         m.plot()
         # plot the real underlying rate function
-        pb.plot(X, np.exp(f_true), '--k', linewidth=2)
+        pb.plot(X, np.exp(f_true), "--k", linewidth=2)
 
     return m
 
-def toy_ARD(max_iters=1000, kernel_type='linear', num_samples=300, D=4, optimize=True, plot=True):
+
+def toy_ARD(
+    max_iters=1000, kernel_type="linear", num_samples=300, D=4, optimize=True, plot=True
+):
     # Create an artificial dataset where the values in the targets (Y)
     # only depend in dimensions 1 and 3 of the inputs (X). Run ARD to
     # see if this dependency can be recovered
@@ -310,14 +377,14 @@ def toy_ARD(max_iters=1000, kernel_type='linear', num_samples=300, D=4, optimize
     Y2 = np.asarray(4 * (X[:, 2] - 1.5 * X[:, 0])).reshape(-1, 1)
     Y = np.hstack((Y1, Y2))
 
-    Y = np.dot(Y, np.random.rand(2, D));
+    Y = np.dot(Y, np.random.rand(2, D))
     Y = Y + 0.2 * np.random.randn(Y.shape[0], Y.shape[1])
     Y -= Y.mean()
     Y /= Y.std()
 
-    if kernel_type == 'linear':
+    if kernel_type == "linear":
         kernel = GPy.kern.Linear(X.shape[1], ARD=1)
-    elif kernel_type == 'rbf_inv':
+    elif kernel_type == "rbf_inv":
         kernel = GPy.kern.RBF_inv(X.shape[1], ARD=1)
     else:
         kernel = GPy.kern.RBF(X.shape[1], ARD=1)
@@ -327,14 +394,17 @@ def toy_ARD(max_iters=1000, kernel_type='linear', num_samples=300, D=4, optimize
     # m.set_prior('.*lengthscale',len_prior)
 
     if optimize:
-        m.optimize(optimizer='scg', max_iters=max_iters)
+        m.optimize(optimizer="scg", max_iters=max_iters)
 
     if plot:
         m.kern.plot_ARD()
 
     return m
 
-def toy_ARD_sparse(max_iters=1000, kernel_type='linear', num_samples=300, D=4, optimize=True, plot=True):
+
+def toy_ARD_sparse(
+    max_iters=1000, kernel_type="linear", num_samples=300, D=4, optimize=True, plot=True
+):
     # Create an artificial dataset where the values in the targets (Y)
     # only depend in dimensions 1 and 3 of the inputs (X). Run ARD to
     # see if this dependency can be recovered
@@ -348,69 +418,73 @@ def toy_ARD_sparse(max_iters=1000, kernel_type='linear', num_samples=300, D=4, o
     Y2 = np.asarray(4 * (X[:, 2] - 1.5 * X[:, 0]))[:, None]
     Y = np.hstack((Y1, Y2))
 
-    Y = np.dot(Y, np.random.rand(2, D));
+    Y = np.dot(Y, np.random.rand(2, D))
     Y = Y + 0.2 * np.random.randn(Y.shape[0], Y.shape[1])
     Y -= Y.mean()
     Y /= Y.std()
 
-    if kernel_type == 'linear':
+    if kernel_type == "linear":
         kernel = GPy.kern.Linear(X.shape[1], ARD=1)
-    elif kernel_type == 'rbf_inv':
+    elif kernel_type == "rbf_inv":
         kernel = GPy.kern.RBF_inv(X.shape[1], ARD=1)
     else:
         kernel = GPy.kern.RBF(X.shape[1], ARD=1)
-    #kernel += GPy.kern.Bias(X.shape[1])
+    # kernel += GPy.kern.Bias(X.shape[1])
     X_variance = np.ones(X.shape) * 0.5
     m = GPy.models.SparseGPRegression(X, Y, kernel, X_variance=X_variance)
     # len_prior = GPy.priors.inverse_gamma(1,18) # 1, 25
     # m.set_prior('.*lengthscale',len_prior)
 
     if optimize:
-        m.optimize(optimizer='scg', max_iters=max_iters)
+        m.optimize(optimizer="scg", max_iters=max_iters)
 
     if plot:
         m.kern.plot_ARD()
 
     return m
 
+
 def robot_wireless(max_iters=100, kernel=None, optimize=True, plot=True):
     """Predict the location of a robot given wirelss signal strength readings."""
-    try:import pods
+    try:
+        import pods
     except ImportError:
-        print('pods unavailable, see https://github.com/sods/ods for example datasets')
+        print("pods unavailable, see https://github.com/sods/ods for example datasets")
         return
     data = pods.datasets.robot_wireless()
 
     # create simple GP Model
-    m = GPy.models.GPRegression(data['Y'], data['X'], kernel=kernel)
+    m = GPy.models.GPRegression(data["Y"], data["X"], kernel=kernel)
 
     # optimize
     if optimize:
         m.optimize(max_iters=max_iters)
 
-    Xpredict = m.predict(data['Ytest'])[0]
+    Xpredict = m.predict(data["Ytest"])[0]
     if plot:
-        pb.plot(data['Xtest'][:, 0], data['Xtest'][:, 1], 'r-')
+        pb.plot(data["Xtest"][:, 0], data["Xtest"][:, 1], "r-")
-        pb.plot(Xpredict[:, 0], Xpredict[:, 1], 'b-')
+        pb.plot(Xpredict[:, 0], Xpredict[:, 1], "b-")
-        pb.axis('equal')
+        pb.axis("equal")
-        pb.title('WiFi Localization with Gaussian Processes')
+        pb.title("WiFi Localization with Gaussian Processes")
-        pb.legend(('True Location', 'Predicted Location'))
+        pb.legend(("True Location", "Predicted Location"))
 
-    sse = ((data['Xtest'] - Xpredict)**2).sum()
+    sse = ((data["Xtest"] - Xpredict) ** 2).sum()
 
-    print(('Sum of squares error on test data: ' + str(sse)))
+    print(("Sum of squares error on test data: " + str(sse)))
     return m
 
+
 def silhouette(max_iters=100, optimize=True, plot=True):
     """Predict the pose of a figure given a silhouette. This is a task from Agarwal and Triggs 2004 ICML paper."""
-    try:import pods
+    try:
+        import pods
     except ImportError:
-        print('pods unavailable, see https://github.com/sods/ods for example datasets')
+        print("pods unavailable, see https://github.com/sods/ods for example datasets")
         return
     data = pods.datasets.silhouette()
 
     # create simple GP Model
-    m = GPy.models.GPRegression(data['X'], data['Y'])
+    m = GPy.models.GPRegression(data["X"], data["Y"])
 
     # optimize
     if optimize:
@@ -419,10 +493,18 @@ def silhouette(max_iters=100, optimize=True, plot=True):
     print(m)
     return m
 
-def sparse_GP_regression_1D(num_samples=400, num_inducing=5, max_iters=100, optimize=True, plot=True, checkgrad=False):
+
+def sparse_GP_regression_1D(
+    num_samples=400,
+    num_inducing=5,
+    max_iters=100,
+    optimize=True,
+    plot=True,
+    checkgrad=False,
+):
     """Run a 1D example of a sparse GP regression."""
     # sample inputs and outputs
-    X = np.random.uniform(-3., 3., (num_samples, 1))
+    X = np.random.uniform(-3.0, 3.0, (num_samples, 1))
     Y = np.sin(X) + np.random.randn(num_samples, 1) * 0.05
     # construct kernel
     rbf = GPy.kern.RBF(1)
@@ -433,20 +515,23 @@ def sparse_GP_regression_1D(num_samples=400, num_inducing=5, max_iters=100, opti
         m.checkgrad()
 
     if optimize:
-        m.optimize('tnc', max_iters=max_iters)
+        m.optimize("tnc", max_iters=max_iters)
 
     if plot:
         m.plot()
 
     return m
 
-def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100, optimize=True, plot=True, nan=False):
+
+def sparse_GP_regression_2D(
+    num_samples=400, num_inducing=50, max_iters=100, optimize=True, plot=True, nan=False
+):
     """Run a 2D example of a sparse GP regression."""
     np.random.seed(1234)
-    X = np.random.uniform(-3., 3., (num_samples, 2))
+    X = np.random.uniform(-3.0, 3.0, (num_samples, 2))
     Y = np.sin(X[:, 0:1]) * np.sin(X[:, 1:2]) + np.random.randn(num_samples, 1) * 0.05
     if nan:
-        inan = np.random.binomial(1,.2,size=Y.shape)
+        inan = np.random.binomial(1, 0.2, size=Y.shape)
         Y[inan] = np.nan
 
     # construct kernel
@@ -456,13 +541,13 @@ def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100, opt
     m = GPy.models.SparseGPRegression(X, Y, kernel=rbf, num_inducing=num_inducing)
 
     # contrain all parameters to be positive (but not inducing inputs)
-    m['.*len'] = 2.
+    m[".*len"] = 2.0
 
     m.checkgrad()
 
     # optimize
     if optimize:
-        m.optimize('tnc', messages=1, max_iters=max_iters)
+        m.optimize("tnc", messages=1, max_iters=max_iters)
 
     # plot
     if plot:
@@ -471,76 +556,79 @@ def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100, opt
     print(m)
     return m
 
+
 def uncertain_inputs_sparse_regression(max_iters=200, optimize=True, plot=True):
     """Run a 1D example of a sparse GP regression with uncertain inputs."""
     fig, axes = pb.subplots(1, 2, figsize=(12, 5), sharex=True, sharey=True)
 
     # sample inputs and outputs
     S = np.ones((20, 1))
-    X = np.random.uniform(-3., 3., (20, 1))
+    X = np.random.uniform(-3.0, 3.0, (20, 1))
     Y = np.sin(X) + np.random.randn(20, 1) * 0.05
     # likelihood = GPy.likelihoods.Gaussian(Y)
-    Z = np.random.uniform(-3., 3., (7, 1))
+    Z = np.random.uniform(-3.0, 3.0, (7, 1))
 
     k = GPy.kern.RBF(1)
     # create simple GP Model - no input uncertainty on this one
     m = GPy.models.SparseGPRegression(X, Y, kernel=k, Z=Z)
 
     if optimize:
-        m.optimize('scg', messages=1, max_iters=max_iters)
+        m.optimize("scg", messages=1, max_iters=max_iters)
 
     if plot:
         m.plot(ax=axes[0])
-        axes[0].set_title('no input uncertainty')
+        axes[0].set_title("no input uncertainty")
     print(m)
 
     # the same Model with uncertainty
     m = GPy.models.SparseGPRegression(X, Y, kernel=GPy.kern.RBF(1), Z=Z, X_variance=S)
     if optimize:
-        m.optimize('scg', messages=1, max_iters=max_iters)
+        m.optimize("scg", messages=1, max_iters=max_iters)
     if plot:
         m.plot(ax=axes[1])
-        axes[1].set_title('with input uncertainty')
+        axes[1].set_title("with input uncertainty")
         fig.canvas.draw()
 
     print(m)
     return m
 
+
 def simple_mean_function(max_iters=100, optimize=True, plot=True):
     """
     The simplest possible mean function. No parameters, just a simple Sinusoid.
     """
-    #create  simple mean function
+    # create  simple mean function
-    mf = GPy.core.Mapping(1,1)
+    mf = GPy.core.Mapping(1, 1)
     mf.f = np.sin
-    mf.update_gradients = lambda a,b: None
+    mf.update_gradients = lambda a, b: None
 
-    X = np.linspace(0,10,50).reshape(-1,1)
+    X = np.linspace(0, 10, 50).reshape(-1, 1)
-    Y = np.sin(X) + 0.5*np.cos(3*X) + 0.1*np.random.randn(*X.shape)
+    Y = np.sin(X) + 0.5 * np.cos(3 * X) + 0.1 * np.random.randn(*X.shape)
 
-    k =GPy.kern.RBF(1)
+    k = GPy.kern.RBF(1)
     lik = GPy.likelihoods.Gaussian()
     m = GPy.core.GP(X, Y, kernel=k, likelihood=lik, mean_function=mf)
     if optimize:
         m.optimize(max_iters=max_iters)
     if plot:
-        m.plot(plot_limits=(-10,15))
+        m.plot(plot_limits=(-10, 15))
     return m
 
+
 def parametric_mean_function(max_iters=100, optimize=True, plot=True):
     """
     A linear mean function with parameters that we'll learn alongside the kernel
     """
-    #create  simple mean function
+    # create  simple mean function
-    mf = GPy.core.Mapping(1,1)
+    mf = GPy.core.Mapping(1, 1)
     mf.f = np.sin
 
-    X = np.linspace(0,10,50).reshape(-1,1)
+    X = np.linspace(0, 10, 50).reshape(-1, 1)
-    Y = np.sin(X) + 0.5*np.cos(3*X) + 0.1*np.random.randn(*X.shape) + 3*X
+    Y = np.sin(X) + 0.5 * np.cos(3 * X) + 0.1 * np.random.randn(*X.shape) + 3 * X
 
-    mf = GPy.mappings.Linear(1,1)
+    mf = GPy.mappings.Linear(1, 1)
 
-    k =GPy.kern.RBF(1)
+    k = GPy.kern.RBF(1)
     lik = GPy.likelihoods.Gaussian()
     m = GPy.core.GP(X, Y, kernel=k, likelihood=lik, mean_function=mf)
     if optimize:
@@ -556,23 +644,27 @@ def warped_gp_cubic_sine(max_iters=100):
     Snelson's paper.
     """
     X = (2 * np.pi) * np.random.random(151) - np.pi
-    Y = np.sin(X) + np.random.normal(0,0.2,151)
+    Y = np.sin(X) + np.random.normal(0, 0.2, 151)
-    Y = np.array([np.power(abs(y),float(1)/3) * (1,-1)[y<0] for y in Y])
+    Y = np.array([np.power(abs(y), float(1) / 3) * (1, -1)[y < 0] for y in Y])
     X = X[:, None]
     Y = Y[:, None]
 
     warp_k = GPy.kern.RBF(1)
     warp_f = GPy.util.warping_functions.TanhFunction(n_terms=2)
     warp_m = GPy.models.WarpedGP(X, Y, kernel=warp_k, warping_function=warp_f)
-    warp_m['.*\.d'].constrain_fixed(1.0)
+    warp_m[".*\.d"].constrain_fixed(1.0)
     m = GPy.models.GPRegression(X, Y)
-    m.optimize_restarts(parallel=False, robust=True, num_restarts=5, max_iters=max_iters)
+    m.optimize_restarts(
-    warp_m.optimize_restarts(parallel=False, robust=True, num_restarts=5, max_iters=max_iters)
+        parallel=False, robust=True, num_restarts=5, max_iters=max_iters
-    #m.optimize(max_iters=max_iters)
+    )
-    #warp_m.optimize(max_iters=max_iters)
+    warp_m.optimize_restarts(
+        parallel=False, robust=True, num_restarts=5, max_iters=max_iters
+    )
+    # m.optimize(max_iters=max_iters)
+    # warp_m.optimize(max_iters=max_iters)
 
     print(warp_m)
-    print(warp_m['.*warp.*'])
+    print(warp_m[".*warp.*"])
 
     warp_m.predict_in_warped_space = False
     warp_m.plot(title="Warped GP - Latent space")
@@ -584,55 +676,88 @@ def warped_gp_cubic_sine(max_iters=100):
     return warp_m
 
 
-
 def multioutput_gp_with_derivative_observations():
-    def plot_gp_vs_real(m, x, yreal, size_inputs, title, fixed_input=1, xlim=[0,11], ylim=[-1.5,3]):
+    def plot_gp_vs_real(
+        m, x, yreal, size_inputs, title, fixed_input=1, xlim=[0, 11], ylim=[-1.5, 3]
+    ):
         fig, ax = pb.subplots()
         ax.set_title(title)
-        pb.plot(x, yreal, "r", label='Real function')
+        pb.plot(x, yreal, "r", label="Real function")
-        rows = slice(0, size_inputs[0]) if fixed_input == 0 else slice(size_inputs[0], size_inputs[0]+size_inputs[1])
+        rows = (
-        m.plot(fixed_inputs=[(1, fixed_input)], which_data_rows=rows, xlim=xlim, ylim=ylim, ax=ax)
+            slice(0, size_inputs[0])
-    f = lambda x: np.sin(x)+0.1*(x-2.)**2-0.005*x**3
+            if fixed_input == 0
-    fd = lambda x: np.cos(x)+0.2*(x-2.)-0.015*x**2
+            else slice(size_inputs[0], size_inputs[0] + size_inputs[1])
-    N=10 # Number of observations
+        )
-    M=10 # Number of derivative observations
+        m.plot(
-    Npred=100 # Number of prediction points
+            fixed_inputs=[(1, fixed_input)],
-    sigma = 0.05 # Noise of observations
+            which_data_rows=rows,
-    sigma_der = 0.05 # Noise of derivative observations
+            xlim=xlim,
-    x = np.array([np.linspace(1,10,N)]).T
+            ylim=ylim,
-    y = f(x) + np.array(sigma*np.random.normal(0,1,(N,1)))
+            ax=ax,
+        )
 
-    xd = np.array([np.linspace(2,8,M)]).T
+    f = lambda x: np.sin(x) + 0.1 * (x - 2.0) ** 2 - 0.005 * x ** 3
-    yd = fd(xd) + np.array(sigma_der*np.random.normal(0,1,(M,1)))
+    fd = lambda x: np.cos(x) + 0.2 * (x - 2.0) - 0.015 * x ** 2
+    N = 10  # Number of observations
+    M = 10  # Number of derivative observations
+    Npred = 100  # Number of prediction points
+    sigma = 0.05  # Noise of observations
+    sigma_der = 0.05  # Noise of derivative observations
+    x = np.array([np.linspace(1, 10, N)]).T
+    y = f(x) + np.array(sigma * np.random.normal(0, 1, (N, 1)))
 
-    xpred = np.array([np.linspace(0,11,Npred)]).T
+    xd = np.array([np.linspace(2, 8, M)]).T
+    yd = fd(xd) + np.array(sigma_der * np.random.normal(0, 1, (M, 1)))
+
+    xpred = np.array([np.linspace(0, 11, Npred)]).T
     ypred_true = f(xpred)
     ydpred_true = fd(xpred)
 
     # squared exponential kernel:
-    se = GPy.kern.RBF(input_dim = 1, lengthscale=1.5, variance=0.2)
+    se = GPy.kern.RBF(input_dim=1, lengthscale=1.5, variance=0.2)
     # We need to generate separate kernel for the derivative observations and give the created kernel as an input:
     se_der = GPy.kern.DiffKern(se, 0)
 
-    #Then
+    # Then
-    gauss = GPy.likelihoods.Gaussian(variance=sigma**2)
+    gauss = GPy.likelihoods.Gaussian(variance=sigma ** 2)
-    gauss_der = GPy.likelihoods.Gaussian(variance=sigma_der**2)
+    gauss_der = GPy.likelihoods.Gaussian(variance=sigma_der ** 2)
 
     # Then create the model, we give everything in lists, the order of the inputs indicates the order of the outputs
     # Now we have the regular observations first and derivative observations second, meaning that the kernels and
     # the likelihoods must follow the same order. Crosscovariances are automatically taken care of
-    m = GPy.models.MultioutputGP(X_list=[x, xd], Y_list=[y, yd],
+    m = GPy.models.MultioutputGP(
-                                 kernel_list=[se, se_der],
+        X_list=[x, xd],
-                                 likelihood_list=[gauss, gauss_der])
+        Y_list=[y, yd],
+        kernel_list=[se, se_der],
+        likelihood_list=[gauss, gauss_der],
+    )
 
     # Optimize the model
     m.optimize(messages=0, ipython_notebook=False)
 
-    #Plot the model, the syntax is same as for multioutput models:
+    # Plot the model, the syntax is same as for multioutput models:
-    plot_gp_vs_real(m, xpred, ydpred_true, [x.shape[0], xd.shape[0]], title='Latent function derivatives', fixed_input=1, xlim=[0,11], ylim=[-1.5,3])
+    plot_gp_vs_real(
-    plot_gp_vs_real(m, xpred, ypred_true, [x.shape[0], xd.shape[0]], title='Latent function', fixed_input=0, xlim=[0,11], ylim=[-1.5,3])
+        m,
+        xpred,
+        ydpred_true,
+        [x.shape[0], xd.shape[0]],
+        title="Latent function derivatives",
+        fixed_input=1,
+        xlim=[0, 11],
+        ylim=[-1.5, 3],
+    )
+    plot_gp_vs_real(
+        m,
+        xpred,
+        ypred_true,
+        [x.shape[0], xd.shape[0]],
+        title="Latent function",
+        fixed_input=0,
+        xlim=[0, 11],
+        ylim=[-1.5, 3],
+    )
 
-    #making predictions for the values:
+    # making predictions for the values:
-    mu, var = m.predict_noiseless(Xnew=[xpred, np.empty((0,1))])
+    mu, var = m.predict_noiseless(Xnew=[xpred, np.empty((0, 1))])
 
     return m

GPy/examples/state_space.py

@@ -4,26 +4,26 @@ import matplotlib.pyplot as plt
 
 import GPy.models.state_space_model as SS_model
 
+
 def state_space_example():
     X = np.linspace(0, 10, 2000)[:, None]
-    Y = np.sin(X) + np.random.randn(*X.shape)*0.1
+    Y = np.sin(X) + np.random.randn(*X.shape) * 0.1
 
     kernel1 = GPy.kern.Matern32(X.shape[1])
-    m1  = GPy.models.GPRegression(X,Y, kernel1)
+    m1 = GPy.models.GPRegression(X, Y, kernel1)
 
     print(m1)
-    m1.optimize(optimizer='bfgs',messages=True)
+    m1.optimize(optimizer="bfgs", messages=True)
 
     print(m1)
 
     kernel2 = GPy.kern.sde_Matern32(X.shape[1])
-    #m2  = SS_model.StateSpace(X,Y, kernel2)
+    # m2  = SS_model.StateSpace(X,Y, kernel2)
-    m2 = GPy.models.StateSpace(X,Y, kernel2)
+    m2 = GPy.models.StateSpace(X, Y, kernel2)
     print(m2)
 
-    m2.optimize(optimizer='bfgs',messages=True)
+    m2.optimize(optimizer="bfgs", messages=True)
 
     print(m2)
 
     return m1, m2
-