%pip install gpy

    \begin{center}
      \begin{tikzpicture}
        % Define nodes
        \draw node[latent] (XS1) {$\mathbf{Z}^{(1)}$};
        \draw node[obs, below right=of XS1] (Y1) {$\dataMatrix^{(1)}$};
        \draw node[latent, above right=of Y1] (X) {$\latentMatrix$};
        \draw node[obs, below right=of X] (Y2) {$\dataMatrix^{(2)}$};
        \draw node[latent, above right=of Y2] (XS2) {$\mathbf{Z}^{(2)}$};
        
        % Connect the nodes
        \draw [->] (XS1) to (Y1);%
        \draw [->] (XS2) to (Y2);%
        \draw [->] (X) to (Y1);%
        \draw [->] (X) to (Y2);%
        
      \end{tikzpicture}
    \end{center}

from IPython.lib.display import YouTubeVideo
YouTubeVideo('UvLI8o8z4IU')

import numpy as np
import GPy
import string

from matplotlib import pyplot as plt
import mlai.plot as plot
import mlai

colors = ["#3FCC94", "#DD4F23", "#C6D63B", "#D44271", 
          "#E4A42C", "#4F9139", "#6DDA4C", "#85831F", 
          "#B36A29", "#CF4E4A"]

def plot_model(X, which_dims, labels):
    fig, ax = plt.subplots(figsize=plot.big_figsize)
    X = X[:,which_dims]
    ulabs = []
    for lab in labels:
        if not lab in ulabs:
            ulabs.append(lab)
            pass
        pass
    for i, lab in enumerate(ulabs):
        ax.scatter(*X[labels==lab].T,marker='o',color=colors[i],label=lab)
        pass
    pass

which = [0,1,2,6,7,9] # which digits to work on
data = pods.datasets.decampos_digits(which_digits=which)
Y = data['Y']
labels = data['str_lbls']


# Model optimization
input_dim = 5 # How many latent dimensions to use
kern = GPy.kern.RBF(input_dim,ARD=True) # ARD kernel
m = GPy.models.BayesianGPLVM(Yn, input_dim=input_dim, kernel=kern, num_inducing=25)

# initialize noise as 1% of variance in data
#m.likelihood.variance = m.Y.var()/100.
m.optimize(messages=1)

# Plotting the model
plot_model(m.X.mean, m.rbf.lengthscale.argsort()[:2], labels.flatten())
pb.legend()
m.kern.plot_ARD()
# Saving the model:
m.pickle('digit_bgplvm_rbf.pickle')

m = GPy.examples.dimensionality_reduction.mrd_simulation(optimize = False, plot=False)
m.optimize(messages = True, max_iters=3e3, optimizer = 'bfgs')

_ = m.X.plot()
m.plot_scales()

import urllib2, os, sys

model_path =  'digit_bgplvm_demo.pickle' # local name for model file
status = ""

re = 0
if len(sys.argv) == 2:
    re = 1

if re or not os.path.exists(model_path): # only download the model new, if it was not already
    url = 'http://staffwww.dcs.sheffield.ac.uk/people/M.Zwiessele/gpss/lab3/digit_bgplvm_demo.pickle'
    with open(model_path, 'wb') as f:
        u = urllib2.urlopen(url)
        meta = u.info()
        file_size = int(meta.getheaders("Content-Length")[0])
        print "Downloading: %s" % (model_path)

        file_size_dl = 0
        block_sz = 8192
        while True:
            buff = u.read(block_sz)
            if not buff:
                break
            file_size_dl += len(buff)
            f.write(buff)
            sys.stdout.write(" "*(len(status)) + "\r")
            status = r"{:7.3f}/{:.3f}MB [{: >7.2%}]".format(file_size_dl/(1.*1e6), file_size/(1.*1e6), float(file_size_dl)/file_size)
            sys.stdout.write(status)
            sys.stdout.flush()
        sys.stdout.write(" "*(len(status)) + "\r")
        print status
else:
    print("Already cached, to reload run with 'reload' as the only argument")

import cPickle as pickle
with open('./digit_bgplvm_demo.pickle', 'rb') as f:
    m = pickle.load(f)

fig = pb.figure('Latent Space & Scales', figsize=(16,6))
ax_latent = fig.add_subplot(121)
ax_scales = fig.add_subplot(122)

fig_out = pb.figure('Output', figsize=(1,1))
ax_image  = fig_out.add_subplot(111)
fig_out.tight_layout(pad=0)

data_show = GPy.plotting.matplot_dep.visualize.image_show(m.Y[0:1, :], dimensions=(16, 16), transpose=0, invert=0, scale=False, axes=ax_image)
lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X.mean.copy(), m, data_show, ax_latent, ax_scales, labels=labels.flatten())

GPy.examples.dimensionality_reduction.bgplvm_simulation()

import pods

data = pods.datasets.singlecell()
Y = data['Y']
Y.describe

import numpy as np

# obtain a centred version of data.
centredY = Y - Y.mean()
# compute inner product matrix
C = np.dot(centredY,centredY.T)
# perform eigendecomposition
V, U = np.linalg.eig(C)
# sort eigenvalues and vectors according to size
ind = V.argsort()
ev = V[ind[::-1]]
U = U[:, ind[::-1]]

import GPy
import matplotlib.pyplot as plt

def plot_labels(ax, x, y, labels, symbols):
    """A small helper function for plotting with labels"""
    # make sure labels are in order of input:
    ulabels = []
    for lab in labels:
        if not lab in ulabels:
            ulabels.append(lab)
    for i, label in enumerate(ulabels):
        symbol = symbols[i % len(symbols)]
        ind = labels == label
        ax.plot(x[ind], y[ind], symbol)
    ax.legend(ulabels)

import matplotlib.pyplot as plt
import mlai.plot as plot
import mlai

fig, ax = plt.subplots(figsize=plot.big_figsize)
plot_labels(ax, U[:, 0], U[:, 1], data['labels'], '<>^vsd')

mlai.write_figure('singlecell-data-pca.svg', directory='./datasets')

import GPy

kernel=GPy.kern.RBF(5,ARD=1)+GPy.kern.Bias(5)
model = GPy.models.BayesianGPLVM(Y.values, 5, num_inducing=15, kernel=kernel)
model.optimize(messages=True)

import matplotlib.pyplot as plt
import mlai.plot as plot
import mlai

fig, ax = plt.subplots(figsize=plot.big_figsize)
model.plot_latent(ax=ax, labels=data['labels'], marker='<>^vsd')

mlai.write_figure('singlecell-bayes-gplvm.svg', directory='./gplvm')

fig, ax = plt.subplots(figsize=plot.big_wide_figsize)
model.kern.plot_ARD(ax=ax)

mlai.write_figure('singlecell-bayes-gplvm-ard.svg', directory='./gplvm')