{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Gaussian Process Latent Variable Model (Bayesian GPLVM)\n",
    "This notebook shows how to use the Bayesian GPLVM model. This is an unsupervised learning method usually used for dimensionality reduction. For an in-depth overview of GPLVMs,see *[1, 2]*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T09:41:37.159594Z",
     "start_time": "2018-06-19T09:41:36.178948Z"
    }
   },
   "outputs": [],
   "source": [
    "import gpflow\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data\n",
    "We are using the \"three phase oil flow\" dataset used initially for demonstrating the Generative Topographic mapping from *[3]*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T09:41:37.164944Z",
     "start_time": "2018-06-19T09:41:37.160737Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of points: 100 and Number of dimensions: 12\n"
     ]
    }
   ],
   "source": [
    "data = np.load('./data/three_phase_oil_flow.npz')\n",
    "Y = data['Y']  # following the GPflow notation we assume this dataset has size [num_data, output_dim]\n",
    "labels = data['labels']  # integer in [0, 2] indicating to which class the datapoint belongs [num_data,]. Not used for model fitting, only for plotting afterwards.\n",
    "\n",
    "print('Number of points: {} and Number of dimensions: {}'.format(Y.shape[0], Y.shape[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model construction\n",
    "\n",
    "We start by initialising the required variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T09:42:12.018601Z",
     "start_time": "2018-06-19T09:41:37.165915Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "latent_dim = 2  # number of latent dimensions\n",
    "num_inducing = 20  # number of inducing pts\n",
    "num_data = Y.shape[0]  # number of data points\n",
    "X_mean_init = gpflow.models.PCA_reduce(Y, latent_dim)  # Initialise via PCA\n",
    "X_var_init = np.ones((num_data, latent_dim))\n",
    "Z_inducing_inputs_init = np.random.permutation(X_mean_init.copy())[:num_inducing]  # Pick inducing inputs randomly from dataset initialisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We construct a Squared Exponential (SE) kernel operating on the two-dimensional latent space. \n",
    "The `ARD` parameter stands for Automatic Relevance Determination, which in practice means that\n",
    "we learn a different lengthscale for each of the input dimensions. See [Manipulating kernels](../advanced/kernels.ipynb) for more information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel = gpflow.kernels.RBF(latent_dim, ARD=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have all the necessary ingredients to construct the model. GPflow contains an implementation of the Bayesian GPLVM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = gpflow.models.BayesianGPLVM(\n",
    "    X_mean=X_mean_init,\n",
    "    X_var=X_var_init,\n",
    "    Y=Y,\n",
    "    kern=kernel,\n",
    "    M=num_inducing,\n",
    "    Z=Z_inducing_inputs_init\n",
    ")\n",
    "\n",
    "# we change the default likelihood variance, which is 1, to 0.01.\n",
    "m.likelihood.variance = 0.01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we optimise the created model. Given that this model has a deterministic evidence lower bound (ELBO), we can use SciPy's L-BFGS-B optimiser."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T09:42:12.018601Z",
     "start_time": "2018-06-19T09:41:37.165915Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO:tensorflow:Optimization terminated with:\n",
      "  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n",
      "  Objective function value: -173.650389\n",
      "  Number of iterations: 1309\n",
      "  Number of functions evaluations: 1364\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:tensorflow:Optimization terminated with:\n",
      "  Message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n",
      "  Objective function value: -173.650389\n",
      "  Number of iterations: 1309\n",
      "  Number of functions evaluations: 1364\n"
     ]
    }
   ],
   "source": [
    "opt = gpflow.train.ScipyOptimizer()\n",
    "opt.minimize(m, maxiter=gpflow.test_util.notebook_niter(10000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model analysis\n",
    "GPflow allows you to inspect the learned model hyperparameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T09:42:12.050803Z",
     "start_time": "2018-06-19T09:42:12.027308Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>class</th>\n",
       "      <th>prior</th>\n",
       "      <th>transform</th>\n",
       "      <th>trainable</th>\n",
       "      <th>shape</th>\n",
       "      <th>fixed_shape</th>\n",
       "      <th>value</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>BayesianGPLVM/X_mean</th>\n",
       "      <td>Parameter</td>\n",
       "      <td>None</td>\n",
       "      <td>(none)</td>\n",
       "      <td>True</td>\n",
       "      <td>(100, 2)</td>\n",
       "      <td>True</td>\n",
       "      <td>[[0.5335858726494092, -3.585462463168241], [-0...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BayesianGPLVM/X_var</th>\n",
       "      <td>Parameter</td>\n",
       "      <td>None</td>\n",
       "      <td>+ve</td>\n",
       "      <td>True</td>\n",
       "      <td>(100, 2)</td>\n",
       "      <td>True</td>\n",
       "      <td>[[0.00025938337775691946, 0.002987260607075345...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BayesianGPLVM/feature/Z</th>\n",
       "      <td>Parameter</td>\n",
       "      <td>None</td>\n",
       "      <td>(none)</td>\n",
       "      <td>True</td>\n",
       "      <td>(20, 2)</td>\n",
       "      <td>True</td>\n",
       "      <td>[[0.723448722354364, 0.0216113992401794], [-1....</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BayesianGPLVM/kern/lengthscales</th>\n",
       "      <td>Parameter</td>\n",
       "      <td>None</td>\n",
       "      <td>+ve</td>\n",
       "      <td>True</td>\n",
       "      <td>(2,)</td>\n",
       "      <td>True</td>\n",
       "      <td>[0.5724092125532642, 2.7044551950083004]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BayesianGPLVM/kern/variance</th>\n",
       "      <td>Parameter</td>\n",
       "      <td>None</td>\n",
       "      <td>+ve</td>\n",
       "      <td>True</td>\n",
       "      <td>()</td>\n",
       "      <td>True</td>\n",
       "      <td>0.7480714466964397</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BayesianGPLVM/likelihood/variance</th>\n",
       "      <td>Parameter</td>\n",
       "      <td>None</td>\n",
       "      <td>+ve</td>\n",
       "      <td>True</td>\n",
       "      <td>()</td>\n",
       "      <td>True</td>\n",
       "      <td>0.005255421275529302</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                       class prior transform  trainable  \\\n",
       "BayesianGPLVM/X_mean               Parameter  None    (none)       True   \n",
       "BayesianGPLVM/X_var                Parameter  None       +ve       True   \n",
       "BayesianGPLVM/feature/Z            Parameter  None    (none)       True   \n",
       "BayesianGPLVM/kern/lengthscales    Parameter  None       +ve       True   \n",
       "BayesianGPLVM/kern/variance        Parameter  None       +ve       True   \n",
       "BayesianGPLVM/likelihood/variance  Parameter  None       +ve       True   \n",
       "\n",
       "                                      shape  fixed_shape  \\\n",
       "BayesianGPLVM/X_mean               (100, 2)         True   \n",
       "BayesianGPLVM/X_var                (100, 2)         True   \n",
       "BayesianGPLVM/feature/Z             (20, 2)         True   \n",
       "BayesianGPLVM/kern/lengthscales        (2,)         True   \n",
       "BayesianGPLVM/kern/variance              ()         True   \n",
       "BayesianGPLVM/likelihood/variance        ()         True   \n",
       "\n",
       "                                                                               value  \n",
       "BayesianGPLVM/X_mean               [[0.5335858726494092, -3.585462463168241], [-0...  \n",
       "BayesianGPLVM/X_var                [[0.00025938337775691946, 0.002987260607075345...  \n",
       "BayesianGPLVM/feature/Z            [[0.723448722354364, 0.0216113992401794], [-1....  \n",
       "BayesianGPLVM/kern/lengthscales             [0.5724092125532642, 2.7044551950083004]  \n",
       "BayesianGPLVM/kern/variance                                       0.7480714466964397  \n",
       "BayesianGPLVM/likelihood/variance                               0.005255421275529302  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.as_pandas_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting vs. Principle Component Analysis (PCA)\n",
    "The reduction of the dimensionality of the dataset to two dimensions allows us to visualise the learned manifold.\n",
    "We compare the Bayesian GPLVM's latent space to the deterministic PCA's one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T09:42:12.370753Z",
     "start_time": "2018-06-19T09:42:12.184055Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlUAAAF1CAYAAADBQh8ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X94XVd95/vP17IcKYKRgKRIlm2CHxi3IdbEvgJSJ6VtTGVAJDEBHFJ+hNJ73Q70IhjG1IbUPWSYicEzZJSh81AXOoRLSuMa4yQ94XFCwq8kTYqxg5wQDEEFYkUGQ5BCZDmW5XX/2OfIR9LeR+fHPmf/OO/X8+iRtM7W1tKRztZ3r/Vd32XOOQEAAKA6i6LuAAAAQBoQVAEAAISAoAoAACAEBFUAAAAhIKgCAAAIAUEVAABACAiqAAANx8x+z8yORN0PpAtBFapmZj8xs0kze9bMfm5mnzez5+Ue22Bm3zKz35jZcTP7ppldOefr/8DMnJn9ZTQ/AYBam3Od+LWZZc1seVT9cc592zm3qlbnN7O3mdnDZjZhZr/IffxeM7Pc4583s1O55+NpM7vHzH4791jGzL7oc86vmdl2n/Y3m9mImTWZ2Rdz19P+Ocf8r1z7O2r1M4OgCuG5wjn3PElrJfVKut7M3iLpnyR9QdIySS+WtF3SFXO+9jpJT0t6V/26CyAC+etEl6SfS/pfEfenJszsQ5IGJe2U1Cnv2vfnki6VtKTg0E/mno9lkn4h6fMLnPoWSe/0aX+npC8656Zzn/9QBddTM2uW9BZJw+X+LCgPQRVC5ZwbkfRVSaslfUrSf3HOfdY5N+6cO+Oc+6Zz7v/JH29mbfJe7O+T9HIz642k4wDqxjl3UtIeSRfm28ys38wOmdkzZvakmWUKHsua2f9beA4zGzKzN+U+/u3cSM/TZnbEzDYVHPcGM/t+brR8xMz+c679D8zsaMFxW83sx7njvp8/d+6xd5vZ/Wb233OjbP9mZq/3+9nMrF3SDZLe65zb45z7jfMccs693Tn3nM/zcULSP0i6aIGnbq+kTjNbV/D9XiTpDfJuXvP2SfqDXF8kqV/SAUnHFzg/qkRQhVDlhvPfIOmEpOXyLpzFXC3pWXkjWvvljVoBSDEzO1fSNZIeKmiekDe60iEvCPiPZrYx99gtkt5R8PX/QVK3pGzuxuweeUHJb0l6m6T/bWb5gO1zkv7MOfd8eUHLfQHd+rGk35PULuljkr5oZl0Fj79a0hFJ50n6pKTP5afy5vhdSedIun2Bp2FGLl3i7ZIOFTvOOTch75paOKr/NklDzrnHCtomJWUl5YPLd2l20IUaIahCWPaZ2Zik+yV9U9L/zLWPLvB110m6LTds/Q+S3pYbqgaQPvnrxLikP5I3PSZJcs59wzl3ODeiPSTpS5J+P/fwHZL+vZm9PPf5O+VdN05JeqOknzjn/o9z7rRz7pCkL0t6a+7YKUkXmtm/c8792jl30K9jzrl/cs49lfv+t0n6kaRXFRzyU+fc3+WuVbfIm8J8sc+pzpP0S+fc6XyDmT1oZmO5nLLXFBz7n3PPxxOSnifp3cWevJxbJG0ys3Nyn78r1zbXFyS9y8xeKGmdvOcQNUZQhbBsdM51OOde4px7r6Rf5dq7gr4gN6r1h5JuzTXdLqlF3l0qgPTZ6JzrkPc6/wtJ3zSzTkkys1eb2ddzC1rG5eUgnSfNTBfeJukdZrZI0rWS/r/cOV8i6dW5oGUsF6S8XV4ukyS9Wd7o+U9zC2V+169jZvYuM3uk4BwX5b9/zrH8B7npOskLhOb6laTzzGxxwfHrcj/3rzT7/+5/z103O51zVzrnflzkucv7pqRnJF1hZv9e0hp5AajfccskbZN0u9+0I8JHUIVaOSLpSXkXtCDvlPc3eKeZHZOXRNkipgCBVHPOTTvn9kqalnRZrvkf5I2mLHfOtUv6jKTC6bVb5AVL6yWdcM79S679SUnfzAUn+bfnOef+Y+57fcc5d5W8qcF9knbP7Y+ZvUTS38kL9F6UC4AenfP9S/Uvkp6TdFUFX7sg55xTbhRK3jX0LufcLwOOu1XSh8TUX90QVKEmci/o/yTpr8zsT8zs35nZIjO7zMx25Q67Tl7uwsUFb2+W9IZc8iWAFDLPVZJeIOnxXPPzJT3tnDtpZq+S9MeFX5MLos5I+h86O0olSf8sb2rwnWbWnHt7pZn9jpktMbO3m1m7c25K3gjPGZ8utUlyyiVym9mfaOGkcV/OuTF517X/bWZvMbPn5659F+e+T6kWmVlLwds5BY99QdLrJL1H/lN/eTdJ+iPn3APl/hyoDEEVasY5t0deMup7JD0lbwn1xyXdbmaXyBu2/xvn3LGCtzvk5RdcG1W/AdTMnWb2rLzg5r9Kuq4gwfq9km4ws9/IK70yb0RJXjCxWtJMDSfn3G8k9clL2H5K3jTdJ+Qli0veaM5PzOwZeVOKb597Uufc9+UFa/8i7zq1WlLFgYhz7pPybio/nDvfzyX9raS/lPRgiae5Vl7Cef5tZmrQOfeEpH+V9zNmi/TjV865eyv4EVAh8wYUAACINzN7l6TNzrnLFjwYiAAjVQCA2MuVYXivpF0LHQtEhaAKABBrZrZBXr7Tz+UltAOxxPQfAABACBipAgAACAFBFQAAQAgWL3xI+M477zx3wQUXRPGtAUTku9/97i+dc+dH3Y9qcf0CGk+p169IgqoLLrhABw4ciOJbA4iImf006j6EgesX0HhKvX4x/QcAABACgioAAIAQEFQBAACEgKAKAAAgBARVAAAAISCoAgAACAFBFQAAQAgIqgAAAEJAUAUAABACgioAAIAQEFQBAACEgKAKAFAz2eGs+vb0qeeWHvXt6VN2OBt1l4CaqTqoMrPlZvZ1M/u+mT1mZgNhdAyYi4szkCzZ4awyD2Y0OjEqJ6fRiVFlHszw2kVqhTFSdVrSh5xzF0q6RNL7zOzCEM4LzODiDCTP4MFBnZw+Oavt5PRJDR4cjKhHQG1VHVQ550adcwdzH/9G0uOSuqs9L1CIizOQPMcmjpXVDiRdqDlVZnaBpDWSHg7zvAAXZyB5Ots6y2oHki60oMrMnifpy5I+4Jx7xufxzWZ2wMwOHD9+PKxviwbBxRlInoG1A2ppapnV1tLUooG1pN4inUIJqsysWV5Adatzbq/fMc65Xc65Xudc7/nnnx/Gt0UD4eJcQ0O7pZsukjId3vuh3VH3CCnRv7JfmXUZdbV1yWTqautSZl1G/Sv7o+4aUBOLqz2BmZmkz0l63Dn3qeq7BMyXvwgPHhzUsYlj6mzr1MDaAS7O1RraLd35fmlq0vt8/Envc0nq2RRdv5Aa/Sv7eZ2iYVQdVEm6VNI7JR02s0dybR9xzt0VwrmBGVyca+DeG84GVHlTk147QRUAlKXqoMo5d78kC6EvAOpt/Gh57UAJssNZRpXRkKioDjSy9mXltQMLoKYcGhlBFRCBfYdGdOmO+/TSrVlduuM+7Ts0Ek1H1m+XmltntzW3eu1ABagph0YWRk4VgDLsOzSibXsPa3JqWpI0MjapbXsPS5I2rqlz3dx83tS9N3hTfu3LvICKfCpUiJpyaGQEVUCd7dx/ZCagypucmtbO/UfqH1RJXgBFEIWQdLZ1anRi1LcdSDum/4A6e2pssqx2IEmoKYdGxkgVYqURVg0t7WjViE8AtbSj1edoIFmoKYdGRlCF2MivGsonueZXDUlK1QV5y4ZVs3KqJKm1uUlbNqyqzTcc2k3OFOoq0TXleL2gCkz/ITYaZdXQxjXduvHq1eruaJVJ6u5o1Y1Xr65NPlW+Yvr4k5Lc2YrpbEUDzMfrBVVipAqx0Uirhjau6a5PUjoV04HS8XpBlRipQmwErQ5i1VAVqJgOlI7XC6pEUIXYSMuqodgU9pSomA6Ug9cLqkRQhdjoX9mvzLqMutq6ZDJ1tXUpsy6TqITXfGHPkbFJOZ0t7EnFdCABeL2gSuRUIVYSvWpIMS3sKbGaCXWR+JIovF5QJYIqRCLxF98AsSzsScV01EFqSqLwekEVmP5D3aV5F/ugAp4U9kTaNUpJFKAYgirUXZovvls2rFJrc9OstpoW9gRiopFKogBBCKpQd2m++Na1sCcQI5REAcipQgTSvot93Qp7IhRmtlzSFyS9WJKTtMs5l/xh0zobWDswK6dKSmZJFKAajFSh7tJSjwqpcVrSh5xzF0q6RNL7zOzCiPuUOGkoiQJUi5Eq1B272CNOnHOjkkZzH//GzB6X1C3p+5F2LIGSXhIFqBZBFSLBxRdxZGYXSFoj6eE57ZslbZakFStW1L1fAJKB6T8AkGRmz5P0ZUkfcM49U/iYc26Xc67XOdd7/vnnR9NBALFHUAWg4ZlZs7yA6lbn3N6o+wMgmQiqADQ0MzNJn5P0uHPuU1H3B0ByEVQBaHSXSnqnpMvN7JHc2xui7lSSZIez6tvTp55betS3py+ZuyMM7ZZuukjKdHjvh3ZH3SMkEInqABqac+5+SRZ1P5IqFXv+De2W7ny/NJXbo3P8Se9ziX0AURZGqgAAFUvFtlP33nA2oMqbmvTagTIQVCF1UjEVASREKradGj9aXjsQgKAKqZKfihidGJWTm5mKILACaiMVe/61LyuvHQhAUIVUScVUBJAgqdh2av12qbl1dltzq9cOlIFEdaRKKqYi4mJot5dTMn7Uu2Nfv52kXcyTim2n8n/X/L2jSgRVSJXOtk6NToz6tqMMrIZCGVKx7VTPJv62UTWm/5AqqZiKkKKvmcNqKAAoGyNVSJVUTEXEYZSI1VAAUDaCKqRO4qciio0S1Suoal/mBXN+7UDakD+IkDD9B8RNHEaJWA2FRpEfGR5/UpI7OzLMNjWoAEEVEDdxqJnTs0m64mapfbkk895fcTN370gf8gcRIqb/gLhZv312TpUUzSgRq6HQCOIwMozUYKQKiBtGiYD6icPIMFKDkSogjhglAuojLiPDSAVGqgAAjYuRYYSIkSoAQGNjZBghYaQKAAAgBARVAAAAISCoAgAACAE5VUAd7Ts0op37j+ipsUkt7WjVlg2rtHFNd9TdAgCEgJGqiGSHs+rb06eeW3rUt6dP2eFs1F1Cje07NKJtew9rZGxSTtLI2KS27T2sfYdGou4a6ozXP5BOBFURyA5nlXkwo9GJUTk5jU6MKvNghgtryu3cf0STU9Oz2ianprVz/5GIeoQo8PoH0ougKgKDBwd1cvrkrLaT0yc1eHAwoh6hHp4amyyrHenE6x9IL3KqInBs4lhZ7YincvOjlna0asQngFra0VrLbiJmeP1HYGi3t0Hy+FFv+5n126lLhZpgpCoCnW2dZbUjfirJj9qyYZVam5tmtbU2N2nLhlU17i3ihNd/nQ3t9rahGX9SkvPe3/l+rx0IGUFVBAbWDqilqWVWW0tTiwbWDkTUI5SrkvyojWu6dePVq9Xd0SqT1N3RqhuvXs3qvwbD67/O7r1h9r5+kvf5vTdE0x+kGtN/Eehf2S/Jy604NnFMnW2dGlg7MNOO+Ks0P2rjmm6CqAbH67/Oxo+W1w5UgaAqIv0r+7mIJhj5UagGr/86al+Wm/rzaQdCxvQfUAHyo4CEWL9dap5zs9Pc6rX7Gdot3XSRlOnw3pN7hTIwUgVUID+FR3V0IObyq/xKWf2XT2rP52Dlk9oLzwMUQVAFVCjR+VEsMa+77HBWO/51h8aeG5MktS9p17ZXb2MasB56NpX2910sqZ3XB0pAUAU0Gu7G6y47nNVfPfBXmjozNdM2fmpc199/vSQRWMUFSe2oEjlVQKP56l+yxLyOssNZfeT+j8wKqPJOu9NUUo+ToOR1ktpRIoIqoJEM7ZYmn/Z/jLvx0OX3+TvjzgQeQyX1GCk3qR2Yg6AKaCTFRqO4Gw+d3z5/c1FJPUZ6NklX3Cy1L5dk3vsrbmZaHCUjpwpoJMVGo7gbD91Co1CLbTGV1OOm1KR2wEcoI1Vm9vdm9gszezSM8wGokaDRqNYX8o+kBoqNQrUvadfHL/s4SepAioQ1/fd5Sa8L6VwAaiUoZ+T1n4imPykXtM/fjt/bofuvvZ+ACkiZUIIq59y3JAVkvwKIDXJG5qnlSHv/yn5l1mXU1dYlk6mrrUuZdRmCKSCl6pZTZWabJW2WpBUrVtTr2yZGdjjLBqs4q5bFOckZmevzkj4t6Qu1ODn7/AGNo26r/5xzu5xzvc653vPPP79e3zYR8suuRydG5eQ0OjGqzIMZZYezUXcNUcgX5xx/UpI7W5yTPchqgpF2AGGhpEIM+C27Pjl9kqKAjarYVhkAgNgiqIqBoGXXFAVsUGyVETtmttnMDpjZgePHj0fdHQAxFVZJhS9J+hdJq8zsqJn9aRjnTbPscFZ9e/rUc0uPzMz3GIoCNii2yogd0hcAlCKURHXn3LVhnKdR5HOo8lN+zrl5x7Q0tVAUsFGt3z57w2OJrTIAIAGY/otA0NYVi2wRy65B2YM6Y6QdQFjYpiYCQblSzjkNXTdU594glih7UDeMtAMICyNVEQjKleps65yVa9W3p4+yCgAAJARBVQSCtq54zbLXUK8KAICEIqiKQNDWFd86+i3qVQEAkFDkVEXEb+uKbd/e5nss9aoAAIg/RqpipFiuFQAAiDeCqhgJyrWiXhUAAPHH9F+M5KcDBw8O6tjEMXW2dWpg7QD1qgAASACCqpjxy7VCOuw7NKKd+4/oqbFJLe1o1ZYNq7RxTXfU3QIAhISgCqiD79zxt3rldz+pb+uXemrJefrkM5u0be8pSSKwAoCUIKgCam1oty46+FdqteckScvsl9rR/FlpStq5fwlBFRCmod3SvTdI40e9TcjXb2+I3QkYCY8Hgiqg1u69Qa16blbTuXZKH168W783dllEnQJSaGj37M3Ix5/0PpdiGViFFQjtOzSibXsPa3JqWpI0MjapbXsPS2IkvN5Y/QfU2vhR3+al9ist7Wgt7RxDu6WbLpIyHd77od0hdhBIiXtvOBtQ5U1Neu0xkw+ERsYm5XQ2ENp3aKTsc+3cf2QmoMqbnJrWzv1HQuotSsVIFVBr7cu8O+Y5RvUibdmwauGvT9jdNxCZgBuYwPYILRQIlTOC9dTYZFntoWvQKVc/jFQBtbZ+u9Q8e0RqUufoqf/rw6UNzVdy983IFhpR+7Ly2iMUFPDkR6zKGcEKGvEueSS8GvmbvvEnJbmzN30Nes0hqAJCtu/QiC7dcZ9eujWrS3fcp33Tl0pX3Cy1L5dkUvtytV79ab3yyj8r7YTl3n1zkUOj8rmBUXOr1x4zQQFPk5nvCNYHbnvEu574BFdbNqxSa3PTrLbW5qbSRsKrVcpNn99NXkpv/Jj+A0IUmDB69aXa+MFHKztpwPRh4N13sYtcgw7Jo0Hk/74TMBX1h799vr740M/mtU87F/g1QQno+Y8jWf230E2fX/rCvvdKZtL0qbNtKUlpIKgCQlQsT6LiC9z67bMvSlLxu+8E5ZUAoevZlIh/zF//wXHf9iazooFV0PVk45ruaFb6LXTT53eTd2Zq/vEpufFj+g8IUU0SRns2zZs+1BU3B198EpRXAoQqQVNKQdeEaefmTeWV+rWRWGjKtZybuRTc+BFUASGqWcJozybpg49KmTHvfbG7uQTllQChSVguYdA1obujVTdevVrdRa4ZdUlAL9VCN33l3Mzlj01QcDwXQRUQokgTRvPKHdkC0iBBNaqk4teKjWu69cDWy/U/r7l43jEmLx8rVord9Pnd5C1qlpqWzG7L3/j5Bcd7N0v//J9q/mOEgZwqIETlJIzWdFuJhOSVAKFJWC5hKdeKjWu6deCnT+vWh36mfJaVk3TrQz/TFx/6mbqTsB1N0OIBv7aeTd7I1NzgWE468PfSiktif10zVyQhrlZ6e3vdgQMH6v59gbiYu0pQ8u5Sb7x6dbwvkFUws+8653qj7ke1uH7F1E0XBSRML/dGT+ZKSMHKS3fcp5EiOVQmL9BKRIBVikyHpIC4JOh3WQelXr+Y/gMiwLYSQMjKySWMWf7VvNp2BbWoFkpKz4cfI2OT2rLne7r4Y3f7nicxiuVgxXTUsRBBFRCByLeVANKmnFzCGOVfLbQHYDlJ6VPTTmOTU1XvJRip9dvljb/5SMAKZoIqoALF7ixLEem2EkBalbpKNkb5VwuNWvsltJcqktHvalfu9WySet+jeYFVQlYwE1QBZQpjd/lYrBIEGlWMarktNGq9cU33giUWKjl/TYQ1rfrGT0lX70rkCmZW/wFlCqNqeqTbSgCNrtxdCmpoaUerbyJ64ah1YbX0iz92t8YmfSqSFzl/kNBXIIe5RVZCVzATVAFlCisfquRtJRKySglIjBjtEbhlwyrflcBBo9bjZQRUxc4TuE+pVHlgFaNp1agQVAFlKuXOMjR+m5GmZONRIFIxGQlZaNR67mhSe2tzySNVxUq0BI24f2j392b1qyzlbv6eQgRVQJnKvbOsaog9zOF0ALEUNGrtN5rU3BSwMm6O7o7WoteZYnsPVjxiFaNp1aiQqA6UqTBx1HR2r65iF8WKk9oZTgcalt9o0tS0U9uSpqCiA5JKW/RSbGS94lWDbJHFSBVQiVLzoapOamc4HWhYQZXUT5ya1k3XXDwzAt7e2iwzaezEVMmj4X4j7oUqXjUY9rRqwnJKUztSlR3Oqm9Pn3pu6VHfnj5lh7NRdwkNqOqk9nKqRANIjX2HRgJHo5bmpvYe2Hq5brrmYrWds7isgEo6O+LeZP7fJRY182JW+b4UqRypyg5nlXkwo5PTJyVJoxOjyjyYkST1r+yPsGdoNFUntcdolRKA2rh+32F96eEnNe2cmsx07auX6+s/OB60A54mnjutl27NquPcZj178rSmznhHlruCL39MOTmidZXAnNJUBlWDBwdnAqq8k9MnNXhwkKAKtVcwXH1Pa6e2L3mz9pxaN/Nw2ResmKxSmpGw4Xggzq7fd1hffOhnM59POzfrcz/51X+/PjF/FWCqauYlMKc0FUFVdjirwYODOjZxTJ1tnRqdGPU97tjEsTr3DA1nTgmEcydHtaP5s3reksW65dlXlXfBimPwQokHIFRfetgnZ7JKNauZV28JzClNfE5VfqpvdGJUTi4woJKkzrbOma8h3wo14TNcvXj6pDJtX9a/7ejXA1svLz2gimMuQYw2ogXSYNoFTfKp4j3/YpEPFYaX95XXHgOJD6r8pvr8tDS1aGDtgG8QlnkwQ2CFcIQ1XB3X4CWBw/Fpx01isgUlijeZzSvd8oJzmxc8X3OTxSMfKgw/uru89hhIfFBVbEpvkXk/XldblzLrMupf2V803wqoWlgbtcY1eInRRrTwH6nnJjFZrn31ct/2JYu9YOuBrZfPjHL/9RWvWHj0KnjgK3nieh0sIvFBVX5Kz88Zd2ZmhCqfoB4UhJFvhVCEVQIhquBlaLd000VSpsN7P3e6MaUlHszsdWZ2xMyeMLOtUfenVNwkJt/HN67WOy5ZoUVzBqwmp87MKxQ8t/Cw3yjX1BlXWeHOOAq63tmi6FMhAiQ+qBpYO6CWppbAx+deYIKCsGLBGVCysCoKRxG8lJLHlcKKyWbWJOlvJL1e0oWSrjWzC6PtVWm4SUyHj29cra72+XlQhZXN9x0a0aU77tMHb3tEknTTNRfrTEA+VsWFO+PG7zooSW46HjmmPhIfVPWv7FdmXUZdbV2BxxReYPyCsPxoFhCKnk3SBx/Vvqse06XP3ayX/kObLt1xX+lb0+TPUe/gpdQ8rtzPp8yY9z7BAVXOqyQ94Zwbds6dkvSPkq6KuE9F5fOoXMBcDzeJyVOsUHDQdlcdATlWqUlUz18HzWfKMw45pj5SUVKhf2W/+lf2q29Pn+/qv8ILTH4asLAEQ+H0IBAGv41Qy96ktN71qRKYvxCSbkmF67aPSnp1RH1Z0NzixnM1L2rmJjGBihUKDtru6pzFi9Ta3BTPwp1h6dkk7d3s/1gMr02JH6kqVOooVP/Kft39lrs1dN2Q7n7L3QRUCF2xPf9iiyT0QGa22cwOmNmB48ePR9qXhVY8uyJL9BFfWzasmpeEng+QgkaxxienSt7cPdESdG1KxUhVXjmjUHMLhjJahTBVvedfFNZvn13YU0pFEnoJRiQVLsFalmub4ZzbJWmXJPX29kYatRSrxSdJp93pxt09Io4Fc0tUrLL5zv1HAkexSi3cue/QSLyqppfzu0rQtSlVQZV0diqwGPYGRK0FDeW3ty5cZyYyjbvP4HckvdzMXiovmHqbpD+OtkvBFtkinXFnih7TkInqKaj2HxQgbdmwqqr9+fYdGtGWPd/T1PTZPQK37PnezPesu3J/Vwm6NqVq+q9ULENGrW3ZsErNc9dIS5o4dbq8hPV6S18S+oKcc6cl/YWk/ZIel7TbOfdYtL0KtlBAJTVoonpcC+aGYG4phXKn+T5252MzAVXe1LTTx+6M6M+8kt9VQq5NqRupKgXLkFFrG9d062N3PjZvw9OpaVfWZqeBEjzNEUfOubsk3RV1P0rR1da14BRgQyaqp3yhRTX78/ltvFysveZS/LtKdVCVz5sanRidGTLvautS+zntGntubN7xDXl3h1D45SuMBVywqs6rSsE0Byr3mmWv0W1Hbgt8/JpV1zRmGkMCN99tWCn+XaV2+q9w+wbp7JD56MSonj31rJoXzc5toVYVKlX3GjIpnuZAcdnhrG5/4vaix1x/yfV16k3MpLTafxg6AnI5g9prLsW/q9QGVcWWHZ92p3Xu4nPV1dYlk83aGxAoV1D5hOempgOXSJetcPsYvzs8KRVD5yhuoXIKxYogp4rfdkoprPYflsyVr5iX49m8yJS58hXRdCjFv6vUTv8tlB/1zKlndP+199epN0izoOm8E1Nn9I5LVujrPzhe3TLmudN9QVIwdI7iiuVSNcxo+0LT3yn4x1yNYqUTYlVSIaW/q9QGVZ1tnUUvQH75U9SuQiWCyidI0td/cFwPbL28um/gN90316LmVAydw1/+2lTMVS+7qjGuV8Wmv1P4T7ocC+3kkLqioDGU2um/Yhst+93RFeZgObmZ2lXZ4Ww9uosEKzadF0qxz1Km9Xx2q0c6zM0PDfKto9+qU48iluKVY9VK5E4OKZPaoGp50cPJAAAeaElEQVTuRsuLbNHM+3xNqsKAidpVqNTGNd2BCZ+LzPTSrdnyN1QuVMq03vSp4onqfjkoSISF8qjyGqYkTKlblsTsb37foRFduuO+6q8HRSRyJ4eUSW1QJZ3d4+/wdYf13y77b2ppapm1CrBwJIraVahG5spXzEtKl6Rp52atCKzoQuq3UsZP0J16Pgdl/ElJ7mwOCoFVrGWHs4GbxPtpmJIwpawcG9ot3f6+2X/zt78vsr/5oBXCYQdWQSuLq15xjJKlOqgqtNBIVNAFqWEuVKjK3IrHTT7TcRUPw89dKWPzgzdJwXfwlGBInFKn/PIaJkldKm3l2Ff/0hu9LTR9ymv3U+NRrXpNyxXblBn1kdpE9bkWGokaWDswaz9AqcEuVKhaYSLoS7f65+JVPAxfuFLGbzVgsRov5KAkTqlTfpKX0tBwJWEWWjk2+XTp7XUopluvablYrvKrRIJ3jGiYoCpoNWD7Oe2Szm6kzOo/hCFoRWAow/Dlbi4aVL249QXV9wU1UWraQUtTS+MFVGGrw2rCcq4HxUoi5B8bGZtUk5mmnVP3nGMSv8ov4TtGhDL9Z2avM7MjZvaEmW0N45xhG1g7MK+KuiQ9e+rZmbyqfA7W0HVDuvstd3OhQsVqPgxfzuai67dLTUvmtz/3G/KqYqqUtAOTJSqgqkei9ozWF5beXoeR3FKvB8Vyrwofk7x8Tal2+VmRSXi6QtVBlZk1SfobSa+XdKGka83swmrPm5dP1uy5pUd9e/oqLnHQv7Jf5y4+d177aXeaFX4IXbW7yoeqZ5O05Hnz289MJeZC1WhKTTtIUkBVj0TtGa//hFe7rdCiZq99rlJXE1ah1OtBsdwrv8fmHpMKCU9XCGP671WSnnDODUuSmf2jpKskfb/aE+eTNfO5BfkVe1JlF5NnTj3j284KP9SC3zB8saH9mpr8tX97Qi5UjaZ/Zb+2frv4oH+SFtEUCxZq8vdfzhT5+u3l5ShWqJRpuWpyr1JTNiHhmy2HMf3XLanwGTiaa6vaQiv2yh3FYoUfolT3u/VCdbgbR7jal7QHPpa0RTSR1E8KmiKfu9JPis0+dMVKIiyUj5masgkJ32y5biUVzGyzmR0wswPHjx8v6WuKrdirpAK6X5X1pF2ckFyRVjtO+IWq0WSHszpx+oTvYx3ndCQql0qqU/2kUsoiBNVsk0rPUayhYrlXfo/NPSYVEr7ZchjTfyOSlhd8vizXNotzbpekXZLU29vrSjlx0Iq9zrbOoqNYQRcbVvih3gqn+4L+6OsybF/uikFEavDgoKbOTM1r7zinQ99+27cj6FF1tmxYNWtPOinkQKDUFWMx3zewlJIIC63+S4UEb7ZszpUU3wSfwGyxpB9KWi8vmPqOpD92zj0W9DW9vb3uwIEDC557bk6VdHYJ8bZvb5ML+DfV1dZF0ITIzd3cNEh3R2v1my4ngJl91znXG3U/qlXq9asaq29ZHfjY4esO1/R710r+BqP3mXu0bck/6cX6pSys4P6miwLycJZ7I095mQ7J9/+GeaNUiF5Ma1SVev2qevrPOXda0l9I2i/pcUm7iwVU5Sjcv89k6mrrmhn2LpYHxabIiINiq3Xyit6tx2zvMtRPfq/SUtuTYOOabj3whl9qsO3/qFPHZWFumVTqijFyC+MtBVtqhVL80zl3l6S7wjjXXP0r+31HmvwqoPs5OX1SNz58I1N+qLti03omFV/9l/ACeKhOfo/SUtsTo1bTb8VWjBWOfLS+wCutUDi1Sm5hfMR8erYUib3t8RvFCjJ+apzRK9RdUBJud0er/m1Hvx7YenlwHkTCC+ChOkHXs2LXuUSoVQ2ioIUYL++bPfIx+bRklisCmrwk6NRLeI0qKcFBlTS/AnqpF5zCsgxArVRVVb3SiwtThqmQ2pXKtZp+C1ox9qO759+cTJ+SlrRFvtIPPlIwPZvooGouvwtREAp+otaqqqpeycUlBfkI8BTLJ020Wpb28KtLFaORj7pu05NUKSj9kqoNlf1KJkyentTYc/NXdVDwE/VQ8eamlVR5TkE+As4KyidNtHqX9ohJde65K4HzhX8lpasUQrVSUPolMUFVdjhbUrL53AtRUFmGxA+jI90qubjE6K4cCFTPGkR12oJmIXXfpidpYlpGoRKJCKqq2QOQgp+IStX7/JX7zycmd+VAbMRk5COSbXqSImUrnRMRVFVSPb1QKofREWuRDPfH5K4ciJUYVOde2tGqEZ8AKjX79VUjZWkLiUhUL7YHIBAn+WTUD9z2SP33+Uv4nllAWlW1EjjtUpa2kIiRqmJ7AAJxse/QiLbs+Z6mpoO3fqr5cH8M7soBzFbKnn4NK2VpC4kIqvyqp5Nsjrj52J2PFQ2oJIb7gUZV8UrgtEtZ2kIigiqSzRF3+w6N6Ncnpooew3A/AMwRk8UEYUlEUCWRbI74yielF9PNcD8A+EtR2kJigiogrvxq0BTqaG3WA1svr2OPACRF1aVXULo61MNKRFBVauFPIArFks+bF5kyV76ijr0BkBRUWq+jOtXDin1JhXzhz9GJUTm5mcKf2eFs1F0DJAUnnzeZaedb/wMXRwC+ilVaR8iK1cMKUeyDqqDCnx+5/yMEVoiFoBo0/2NTCAHV0G7ppoukTIf3ns2RgdSg0nod1akeVuyDqqACn2fcGUasEAsb13TrxqtXq7ujVSYvKf3Gq1eHE1Dd+f5cDRd3driawApIhaBRbkqv1EBQ3auQ62HFPqgqVuAzv1UNELWNa7r1wNbL9W87+vXA1svDmfKr03A1gGhQab2O1m/36l8VqkE9rNgHVQNrB9TS1BL4OFvVILVStn0DgNlqNsqN+eq0jVfsV//lV/l95P6P6Iw7M+9xtqpBaqVs+4Y4MrO3SspI+h1Jr3LOHYi2R2g0VFqvozrUw4p9UCWdDazYqgZJUU7tmcBjU7Z9Q0w9KulqSX8bdUcAJF8igiqJrWqQHOXUnil+bLq2b4gj59zjkmRmUXcFQAokJqiS2KoGyVCs9szcoGrBY1O0fUOSmdlmSZslacWKFRH3BkBcJSqoApJgJKDGjF87dWpqz8y+Jskv+fKjzrnbSzmHc26XpF2S1Nvb60LsHoAUIagCQtZkpmk3//9uk88U09KOVt9gizo14XHOvTbqPgBoDLEvqQAkjV9AFdROnRoASA+CKiBk3QGjTDPtBVvPbPzGBn3hlT+lTk1EzOxNZnZU0u9KyprZ/qj7BCC5mP4DQrZlw6pZK/qkgtEnn53SX3n4r/VADYrQYWHOua9I+krU/QCQDgRVQMjyo0y+tadu8t965sRXt+uP7jqvpLpWAIB4IqgCaiCwSnLAFjMtJ45p5Dkv2CpW1wpIq3IK5gJxRU4VUE8BW8w85V406/N8rSqgEeSL4I6MTcrp7I3FvkMjUXcNKAtBFVBPPjuln3BL9MnT8/OpqFWFRlGsCC6QJEz/AfXUM3/rmU9OvFl3PPeqeYdSqwqNgiK4SAuCKqBKZeeCzNl65uJDI2oNWi0INACK4CItCKqAKgRtiHzgp0/r6z84XlKgVXS1INAAipYhARKEoAqoQlAuyK0P/Uz5+umlrOYLXC0INABuLJAWJKoDVQjK+Zi7IU1VSbcFFdh100Xe50DKbFzTrS0bVmlpR6ueGpvUzv1HWP2HxGGkCqhCUC6In4qSbn0qsOvO93sfU4EdKRI0lS5Rrw3JwUgVUAW/DZEt4NiKkm7v9a/ArntvKP9cQIxRVgFpQFAFVGHjmm7dePXqWRsiv/2SFfMCrYqTbgMqsAe2AwlFWQWkAdN/QJX8ksx7X/LCypNuh3afrWNliyQ3Pf+YgMrsQFJRVgFpQFAF1EDFq/nm5lD5BVTNrV5ldiBFKKuANCCoAuLEL4eqUPtyL6AiSR0pQ1kFpAFBFRAnxXKl8iNUBFRIKeq1IelIVAfipFiuFKv+ACDWCKqAOFm/3RuRCsKqPwCILab/gDjJT+195c9Z9QcACcNIFRA3PZukN31m/ogVq/4AINYIqoA46tkkXXGzt9pP5r2/4maS1AEgxpj+A+KqZxNBFAAkCCNVAAAAISCoAgAoO5xV354+9dzSo749fcoOZ6PuEpA4TP8BQIPLDmeVeTCjk9MnJUmjE6PKPJiRJPWv7I+wZ0CyMFIFAA1u8ODgTECVd3L6pAYPDkbUIyCZGKkCKrTv0Aj7lCHRssNZ3fjwjRo/Ne77+LGJY3XuEZBsjFQBFdh3aETb9h7WyNiknKSRsUl94LZHtOaGu7Xv0EjU3QMWlB3O6vr7rw8MqCSps62zjj0Cko+gCqjAzv1HNDk1v+L5r09MadvewwRWiL3Bg4M67U4HPt7S1KKBtQN17BGQfARVQAWeGpsMfGxyalo79x+pY2+A8i00tZdZlyFJHSgTQRXKwrJrz9KOIpseq3jQBcTBQlN7BFRA+QiqULL8suvRiVE5uZll140YWG3ZsEqtzU2Bjy8UdAFRKza1176kvY49AdKDoAolY9n1WRvXdOvGq1ero7V53mOtzU3asmFVBL0CSte/sl/XrLpmXvtiW6xtr94WQY+A5KOkAkoWlIPRCMuug8onbFzT7f9Y0wPSTTdI40el9mXS+u3s44fYuf6S67Xmt9Zo8OCgjk0cU2dbpwbWDjD1B1SoqqDKzN4qKSPpdyS9yjl3IIxOIZ462zo1OjHq255m+fIJ+dV+I2OT2rb3sCTNBFaz6lMN7ZbufL80lcurGn/S+1wisELs9K/sJ4gCQlLt9N+jkq6W9K0Q+oKYe82y18xra4Rl137lE4qu8Lv3hrMBVd7UpNcOxASLToDwVTVS5Zx7XJLMLJzeILayw1nd/sTt89qvetlVqb/LDVrJF7jCb/xoee1AnbHXH1AbdUtUN7PNZnbAzA4cP368Xt8WIfFLUpekbx09O0iZ1jvfoJV8gSv82peV1w7UGYtOgNpYMKgys6+Z2aM+b1eV842cc7ucc73Oud7zzz+/8h4jEgslqae53IJf+YSiK/zWb5ea5wRcza1eO2LFzHaa2Q/MbMjMvmJmHVH3qR4aedEJUEsLBlXOudc65y7yeZs/F4TUCkpGz7en+c43Xz6hu6NVJqm7o1U3Xr3ad/PkfYdGdOld52lg4k90TOfLyaT25dIVN5OkHk/3SLrIOdcj6YeSGqKWwEKvZwCVoaQCSjKwdmBWDkZePnk97Xe+81b4+ShcJTiiy3T7ycvU2tykG/9gtTb2FP9aRMM5d3fBpw9JektUfaknv9dzIyw6AWqtqpwqM3uTmR2V9LuSsma2P5xuIW76V/brqpfNn/G9/YnblR3OcuerClYJIm7eI+mrUXeiHvpX9iuzLqOuti6ZTF1tXez1B4Sg2tV/X5H0lZD6gpgrTErPy0/xcedbwSpB1IWZfU2SX3T/0Xwag5l9VNJpSbcGnGOzpM2StGLFihr1tL6oTwWEj+k/lKzYFF/+4tzIlZmXdrRqxCeAYh/AaDnnXlvscTN7t6Q3SlrvnHMB59glaZck9fb2+h4DAARVKNlCFdUb/c53y4ZVsyqvS+wDGHdm9jpJH5b0+865E1H3B0CysaEySjawdkAtTS2z2hptii/Q0G5t/MYGfb/pbXqoZUBXLrpfTWYzOVX7Do1E3UP4+7Sk50u6x8weMbPPRN0hAMnFSBVKxhRfgIK9/kxSp45rR/NnpSnpDneZRsYmtWXP9yRpwRWEqC/n3Mui7gOA9CCoQlkafYrPl89ef+faKX148W7dceoySdLUtNPH7nyMoAoAUozpP6BaAXv6LbVfzfr81yem6tEbAEBECKqAagXs6feUe1GdOwIAiBJBFVAtn73+Trgl+uTp2dvSdLQ217NXAIA6I6cKqFZ+T797b5DGj+pEa6c++purdceZS2cOaV5kylz5iog6CACoB4IqIAw9m2aCq3Ml/f6hEX3jzsdm8qjazuGlBgBpx/QfUCMnp87MfDw2OaVtew9TrwoAUoygCqgBNlcGgMZDUAXUQDWbK2eHs+rb06eeW3rUt6dP2eFs2N0DANQAQRVQA0GbKC+0uXJ2OKvMgxmNTozKyWl0YlSZBzMEVgCQAARVKBkjKKXbsmGVWpubZrWVsrny4MFBnZw+Oavt5PRJDR4cDL2PAIBwsSQJJcmPoOT/4edHUCSxbY2P/HY0O/cf0VNjk1ra0aotG1YtuE3NsYljZbUDAOKDoAolKTaCQlDlb+Oa7rL3+uts69ToxKhvOwAg3pj+Q0kYQamPgbUDamlqmdXW0tSigbUDEfUIAFAqRqpQEkZQ6iM/6jd4cFDHJo6ps61TA2sHGA0EgAQgqEJJBtYOzMqpkhhBqZX+lf0EUQCQQARVKAkjKAAAFEdQhZIxggIAQDAS1QEAAEJAUAUAABACgioAAIAQEFQBAACEgKAKAAAgBARVQC0N7ZZuukjKdHjvh3ZH3SMAQI1QUgGolaHd0p3vl6Ymvc/Hn/Q+l6SeTdH1CwBQE4xUAbVy7w1nA6q8qUmvHQCQOgRVQK2MHy2vHQCQaARVQK20LyuvHUDsZYez6tvTp55betS3p0/Z4WzUXUKMEFQBtbJ+u9TcOrutudVrB5A42eGsMg9mNDoxKien0YlRZR7MEFhhBkEVUCs9m6Qrbpbal0sy7/0VN5OkDiTU4MFBnZw+Oavt5PRJDR4cjKhHiBtW/wG11LOJIApIiWMTx8pqR+NhpAoAgBJ0tnWW1Y7GQ1AFAEAJBtYOqKWpZVZb86JmnZg6QeI6JDH9BwBASfpX9kvycquOTRxT+zntevbUsxo/NS5JM4nrhceisTBSBQBAifpX9uvut9ytoeuG1Lq4Vafd6VmPk7je2AiqAACoAInrmIugCgCACpC4jrkIqgAAqIBf4npLU4sG1g5E1CNEjUR1AAAqMDdxvbOtUwNrB0hSb2AEVQAalpn9F0lXSToj6ReS3u2ceyraXiFJ+lf2E0RhBtN/ABrZTudcj3PuYkn/LImNGQFUjKAKiRK4Q/zQbummi6RMh/d+aHe0HUUiOOeeKfi0TZKLqi8Ako/pPyRGfof4/IamM4X2fvaQ+h/4O2lq0jtw/Enpzvd7H7PvHhZgZv9V0rskjUv6w4BjNkvaLEkrVqyoX+cAJAojVUiMwB3ih79yNqDKm5qU7r2hjr1DXJnZ18zsUZ+3qyTJOfdR59xySbdK+gu/czjndjnnep1zveeff349uw8gQRipQmIEFtoLujUYP1q7ziAxnHOvLfHQWyXdJemva9gdACnGSBUSI7DQ3pmAL2hfVrvOIBXM7OUFn14l6QdR9QVA8hFUITECC+2tfJPU3Dr74OZWaT0LubCgHbmpwCFJfZKo2gigYkz/ITGKFtp74Wovh2r8qDdCtX47SepYkHPuzVH3AUB6EFQhUQIL7fVsilUQte/QiHbuP6Knxia1tKNVWzas0sY13VF3CwBQQwRVQMj2HRrRtr2HNTk1LUkaGZvUtr2HJYnACgBSjJwqIEBgodEF7Nx/ZCagypucmtbO/Udq0U0AQEwwUgX4CCw0Ki24z9dTY5NltQMA0oGRKsBHYKHRg4MLfu3Sjtay2gEA6UBQBfgILDQa0F5oy4ZVam1umtXW2tykLRtWhdI3AEA8EVQBPgILjQa0F9q4pls3Xr1a3R2tMkkdrc1qaV6kD972iC7dcZ/2HRoJubcAgDggqAJ8BBYaXVtabciNa7r1wNbLddM1F+u502f06xNTcjq7EpDACgDSh6AK8NG/sl+ZdRl1tXXJZOpq61JmXWbBJPW5WAkIAI2jqtV/ZrZT0hWSTkn6saQ/cc6NhdExIGqBhUbLwEpAAGgc1Y5U3SPpIudcj6QfStpWfZeA9GAlIAA0jqqCKufc3c6507lPH5K0rPouAenBSkAAaBxhFv98j6TbQjwfkHj5bWnYBxAA0m/BoMrMvibJbx35R51zt+eO+aik05JuLXKezZI2S9KKFSsq6iyQRBvXdBNEAUADWDCocs69ttjjZvZuSW+UtN4554qcZ5ekXZLU29sbeBwAAEASVbv673WSPizp951zJ8LpEgAAQPJUu/rv05KeL+keM3vEzD4TQp8AAAASp6qRKufcy8LqCAAAQJJRUR0AACAEBFUAAMRIdjirvj196rmlR317+pQdzkbdJZQozDpVAEqw79AIdasA+MoOZ5V5MKOT0yclSaMTo8o8mJGkqrfNQu0xUgXU0b5DI9q297BGxiblJI2MTWrb3sPad2gk6q4BiIHBg4MzAVXeyemTGjw4GFGPUA6CKqCOdu4/osmp6Vltk1PT2rn/SEQ9AhAnxyaOldWOeCGoAuroqbHJstoBNJbONr8NTILbES8EVUAdLe1oLasdQGMZWDuglqaWWW0tTS0aWDsQUY9QDoIqoI62bFil1uamWW2tzU3asmFVRD0CECf9K/uVWZdRV1uXTKauti5l1mVIUk8IVv8BdZRf5cfqPwBB+lf2E0QlFEEVUGcb13QTRAFACjH9BwAAEAKCKgAAgBAQVAEAAISAoAoAACAEBFUAAAAhIKgCAAAIAUEVAABACAiqAAAAQkBQBQAAEAKCKgAAgBAQVAEAAITAnHP1/6ZmxyX9VNJ5kn5Z9w6UJq59o1/li2vf4tovqTZ9e4lz7vyQz1l3BdevNIvz32Y9NPLPz8/ur6TrVyRB1cw3NzvgnOuNrANFxLVv9Kt8ce1bXPslxbtvqL1G//038s/Pz17dz870HwAAQAgIqgAAAEIQdVC1K+LvX0xc+0a/yhfXvsW1X1K8+4baa/TffyP//PzsVYg0pwoAACAtoh6pAgAASIW6BlVmttPMfmBmQ2b2FTPrCDjudWZ2xMyeMLOtderbW83sMTM7Y2aB2f9m9hMzO2xmj5jZgRj1q67PmZm90MzuMbMf5d6/IOC46dxz9YiZ3VHjPhV9DszsHDO7Lff4w2Z2QS37U0a/3m1mxwuep/+7Tv36ezP7hZk9GvC4mdnNuX4PmdnaevQL8VDqtSdNovjfExcLXQ/SzMyWm9nXzez7ub/5gYpP5pyr25ukPkmLcx9/QtInfI5pkvRjSSslLZH0PUkX1qFvvyNplaRvSOotctxPJJ1Xx+dswX5F8ZxJ+qSkrbmPt/r9LnOPPVun52nB50DSeyV9Jvfx2yTdFpN+vVvSp+v1N1XwfV8jaa2kRwMef4Okr0oySZdIerjefeQturdSr4lpeYvqf09c3ha6HqT5TVKXpLW5j58v6YeV/u7rOlLlnLvbOXc69+lDkpb5HPYqSU8454adc6ck/aOkq+rQt8edc0dq/X3KVWK/onjOrpJ0S+7jWyRtrPH3W0gpz0Fhn/dIWm9mFoN+RcI59y1JTxc55CpJX3CehyR1mFlXfXqHqMX1mlhDsX2t1kMJ14PUcs6NOucO5j7+jaTHJXVXcq4oc6reI+8ueK5uSU8WfH5UFf5wNeIk3W1m3zWzzVF3JieK5+zFzrnR3MfHJL044LgWMztgZg+ZWS0Dr1Keg5ljcsH9uKQX1bBPpfZLkt6cm2LbY2bLa9ynUsX9tQiEib93KJcWskbSw5V8/eIwOyNJZvY1SZ0+D33UOXd77piPSjot6dawv3+1fSvBZc65ETP7LUn3mNkPchF+1P0KXbF+FX7inHNmFrSM9CW552ulpPvM7LBz7sdh9zXh7pT0Jefcc2b2Z/JG0y6PuE9oAHG99gBRMLPnSfqypA84556p5ByhB1XOudcWe9zM3i3pjZLWu9wE5hwjkgrv1Jfl2mretxLPMZJ7/wsz+4q8IeOqgqoQ+lWT56xYv8zs52bW5ZwbzU0J/SLgHPnna9jMviHvDqAWQVUpz0H+mKNmtlhSu6Rf1aAvZfXLOVfYh8/Ky1eLg5q9FhEPYVwTU4S/9wZmZs3yAqpbnXN7Kz1PvVf/vU7ShyVd6Zw7EXDYdyS93MxeamZL5CUU13TVWKnMrM3Mnp//WF7ifRxWSkTxnN0h6brcx9dJmndXa2YvMLNzch+fJ+lSSd+vUX9KeQ4K+/wWSfcFBPZ17decPKUr5c3nx8Edkt6VWwV4iaTxgilfIG1i+78HtZXLrf2cpMedc5+q6mR1zrB/Qt6c9SO5t/xKrKWS7io47g3ysu9/LG8Yuh59e5O8OfTnJP1c0v65fZO3KuR7ubfH6tG3UvoVxXMmLxfpXkk/kvQ1SS/MtfdK+mzu43WSDueer8OS/rTGfZr3HEi6QV4QL0ktkv4p93f4r5JW1ulva6F+3Zj7e/qepK9L+u069etLkkYlTeX+xv5U0p9L+vPc4ybpb3L9PqwGWAHG26y/D99rT5rfovjfE5c3v+tB1H2q489+mbx86aGC+OQNlZyLiuoAAAAhoKI6AABACAiqAAAAQkBQBQAAEAKCKgAAgBAQVAEAAISAoAoAACAEBFUAAAAhIKgCAAAIwf8PPJnXHUkvLk8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "XPCAplot = gpflow.models.PCA_reduce(Y, 2)\n",
    "# when a model is trained we can access the values from the tensor as a `np.ndarray` using `read_value()`.\n",
    "GPLVM_X_mean = m.X_mean.read_value()\n",
    "\n",
    "f, ax = plt.subplots(1, 2, figsize=(10,6))\n",
    "\n",
    "for i in np.unique(labels):\n",
    "    ax[0].scatter(XPCAplot[labels==i, 0], XPCAplot[labels==i, 1], label=i)\n",
    "    ax[1].scatter(GPLVM_X_mean[labels==i, 0], GPLVM_X_mean[labels==i, 1], label=i)\n",
    "    ax[0].set_title('PCA')\n",
    "    ax[1].set_title('Bayesian GPLVM')\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\\[1\\] Lawrence, Neil D. 'Gaussian process latent variable models for visualisation of high dimensional data'. *Advances in Neural Information Processing Systems*. 2004.\n",
    "\n",
    "\\[2\\] Titsias, Michalis, and Neil D. Lawrence. 'Bayesian Gaussian process latent variable model'. *Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics*. 2010.\n",
    "\n",
    "\\[3\\] Bishop, Christopher M., and Gwilym D. James. 'Analysis of multiphase flows using dual-energy gamma densitometry and neural networks'. *Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment* 327.2-3 (1993): 580-593."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}