{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bayesian GPLVM\n",
    "--\n",
    "This notebook shows how to use the Bayesian GPLVM model."
   ]
  },
  {
   "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",
    "from gpflow import kernels\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "%matplotlib inline\n",
    "np.random.seed(42)\n",
    "gpflow.settings.numerics.quadrature = 'error'  # throw error if quadrature is used for kernel expectations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data\n",
    "\n",
    "We are using the \"three phase oil flow\" data set used initially for demonstrating the Generative Topographic mapping from *Bishop, C. M. and James, G. D. (1993): Analysis of multiphase flows using dual-energy gamma densitometry and neural networks. Nuclear Instruments and Methods in Physics Research A327, 580-593*."
   ]
  },
  {
   "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 x Number of dimensions (100, 12)\n"
     ]
    }
   ],
   "source": [
    "data = np.load('data/three_phase_oil_flow.npz')\n",
    "Y = data['Y']\n",
    "labels = data['labels']\n",
    "\n",
    "print('Number of points x Number of dimensions', Y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model construction\n",
    "Create Bayesian GPLVM model using additive kernel."
   ]
  },
  {
   "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": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "lib/python3.5/site-packages/tensorflow/python/ops/gradients_impl.py:100: UserWarning: Converting sparse IndexedSlices to a dense Tensor of unknown shape. This may consume a large amount of memory.\n",
      "  \"Converting sparse IndexedSlices to a dense Tensor of unknown shape. \"\n"
     ]
    },
    {
     "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: -148.020057\n",
      "  Number of iterations: 5402\n",
      "  Number of functions evaluations: 5551\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: -148.020057\n",
      "  Number of iterations: 5402\n",
      "  Number of functions evaluations: 5551\n"
     ]
    }
   ],
   "source": [
    "Q = 5\n",
    "M = 20  # number of inducing pts\n",
    "N = Y.shape[0]\n",
    "X_mean = gpflow.models.PCA_reduce(Y, Q) # Initialise via PCA\n",
    "Z = np.random.permutation(X_mean.copy())[:M]\n",
    "\n",
    "fHmmm = False\n",
    "if(fHmmm):\n",
    "    k = (kernels.RBF(3, ARD=True, active_dims=slice(0,3)) + \n",
    "         kernels.Linear(2, ARD=False, active_dims=slice(3,5)))\n",
    "else:\n",
    "    k = (kernels.RBF(3, ARD=True, active_dims=[0,1,2]) + \n",
    "         kernels.Linear(2, ARD=False, active_dims=[3, 4]))\n",
    "    \n",
    "m = gpflow.models.BayesianGPLVM(X_mean=X_mean, X_var=0.1*np.ones((N, Q)), Y=Y,\n",
    "                                kern=k, M=M, Z=Z)\n",
    "m.likelihood.variance = 0.01\n",
    "\n",
    "opt = gpflow.train.ScipyOptimizer()\n",
    "m.compile()\n",
    "opt.minimize(m, maxiter=gpflow.test_util.notebook_niter(10000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute and sensitivity to input\n",
    "Sensitivity is a measure of the importance of each latent dimension. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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/kern/kernels/0/lengthscales</th>\n",
       "      <td>Parameter</td>\n",
       "      <td>None</td>\n",
       "      <td>+ve</td>\n",
       "      <td>True</td>\n",
       "      <td>(3,)</td>\n",
       "      <td>True</td>\n",
       "      <td>[2.92768860478, 7.13288180325, 1.45483099534]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BayesianGPLVM/kern/kernels/0/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.6787500924312497</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BayesianGPLVM/kern/kernels/1/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.027179996202905617</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                               class prior transform  \\\n",
       "BayesianGPLVM/kern/kernels/0/lengthscales  Parameter  None       +ve   \n",
       "BayesianGPLVM/kern/kernels/0/variance      Parameter  None       +ve   \n",
       "BayesianGPLVM/kern/kernels/1/variance      Parameter  None       +ve   \n",
       "\n",
       "                                           trainable shape  fixed_shape  \\\n",
       "BayesianGPLVM/kern/kernels/0/lengthscales       True  (3,)         True   \n",
       "BayesianGPLVM/kern/kernels/0/variance           True    ()         True   \n",
       "BayesianGPLVM/kern/kernels/1/variance           True    ()         True   \n",
       "\n",
       "                                                                                   value  \n",
       "BayesianGPLVM/kern/kernels/0/lengthscales  [2.92768860478, 7.13288180325, 1.45483099534]  \n",
       "BayesianGPLVM/kern/kernels/0/variance                                 0.6787500924312497  \n",
       "BayesianGPLVM/kern/kernels/1/variance                               0.027179996202905617  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m.kern.as_pandas_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T09:42:12.057723Z",
     "start_time": "2018-06-19T09:42:12.052705Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.28140387  0.11550211  0.56629458]\n"
     ]
    }
   ],
   "source": [
    "kern = m.kern.kernels[0]\n",
    "sens = np.sqrt(kern.variance.read_value())/kern.lengthscales.read_value()\n",
    "print(sens)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T09:42:12.178982Z",
     "start_time": "2018-06-19T09:42:12.059247Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEWCAYAAACdaNcBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAFEdJREFUeJzt3X+0XWV95/H3hwSsbZRWk1okgTAlYxt/gQa0a5SmLXWAVnAN2ELFitqmzppMcVVbcXBYgnas1Fq1zVojDgrlhwFp60qnaUNrTfEnTUBEQ8RmEEgAS0QQQQQC3/lj74tn7tybe+7Nudzk8f1a667s/exnP/t7zj3rc/Z99jk7qSokSW3Zb64LkCSNnuEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw10jleSBJP9uN9u3JFk5xDivSXL1SIuboSRnJPnsXNcxkSSH9M/5vLmuRXsXw71hSV6W5PNJvpPk20k+l+So2TxmVS2oqlv641+U5N3jtj+3qjYOMc5lVfWKsfUkleTwmdSUZGWSHTPZdwbHemeSS0c43sYkvzXZ9qq6vX/OHxvVMXdTy61Jjp3t42g05s91AZodSZ4O/G/gPwNXAgcALwcensu6JD1JqsqfBn+AFcB9U/R5A7AVuBfYABw6sK2ANwH/CtwHrAHSbzsc+GfgO8C3gCvG7Xc4sAp4FHgEeAD4m377rcCxwLOBh4BnDOx7ZD/e/sAZwGf79mv6cR/sx/p14KvAKwf23b/f98hxj/HH+uM83u/7QH/spwAfAO7sfz4APGWS5+mJWvr1DwLbgfuB64CX9+3H9Y/30f44X+7bDwQuBO4C7gDeDcwbHBt4X/97+AZwfL/tD4HHgO/34/35BLUt7Z+b+f36RuBdwOeA7wJXAwvH9V3VP+a7gLcOjHUR8O6B9ZXAjn75kv45fKiv5Q+AHwEuBe6he41sAp411699f7ofp2Xa9XXgsSQXJzk+yU8MbkxyEvDfgP8ELAI+A3x83Bi/ChwFvAD4NeA/9u3voguNnwAWA382/uBVdQFwGXB+ddMGrxy3/U7gC8DJA82/AVxVVY+O63tMv/jCfqwrgL8ATh/odgJwV1V9ady+DwLHA3f2+y7oj3028FLgCOCFwNHAO8Y/jkls6vd7BnA58IkkP1JVfw/8D7o3uwVV9cK+/0XALro3vSOBVwCDUy0vAW4GFgLnAxcmSVWdTfd7Wd2Pt3rI+n4DeD3wk3R/sb113PZfAJb1dbxtmKmWqnotcDvdG+qCqjofeB3dG9cS4Jl0JwMPDVmjZpnh3qiquh94Gd2Z2keAnUnWJXlW3+VNwHuqamtV7aILpSOSHDowzB9V1X1VdTvwabpAg+7M9FDg2VX1/aqa6cXGy4HTAJIEOLVvG8alwAn99BPAa+nOLof1GuC8qrq7qnYC5/ZjTKmqLq2qe6pqV1X9Cd1fAc+ZqG//fJ8AvLmqHqyqu4E/pXusY26rqo9UN29+MXAQ8Kz/f7Shfayqvl5VD9FNyR0xbvu5fS1fAT5G/zuYgUfpQv3wqnqsqq7rX3faCxjuDeuD+4yqWgw8j2464gP95kOBDya5L8l9wLeBAAcPDPHNgeXvAQv65T/o+/5L/+mXN8ywxL8Efi7JQcAxdH/2f2aYHfuz788BJyf5cbqz88umcexnA7cNrN/Wt00pyVuTbO0vVN9Hd/a6cJLuh9JNGd018Fx/mO6seswTz3NVfa9fXMDMTfZ7G7N9YHnoxz2BS+im89YmuTPJ+Un2n+FYGjEvqP6QqKqvJbkI+J2+aTvwh1U1nUAcG+ubwG9D94kc4B+TXFNV28Z3nWKce/uPO/468LPA2qqazm1KL6ab3pgPfKGq7pjsUBO03UkXvFv69UP6tt1K8nK6N7dfArZU1eNJ7qV7s5voWNvpLmIv7P9Cmq7ZuG3rEuBr/fLg434Q+NGBfj+1u1r66bNzgXOTLAXW000vXTjacjUTnrk3KsnPJHlLksX9+hK6P7+/2Hf5n8Dbkzy3335gklcPOfarx8aluwhYdGfd4/0bMOln3nuXA78JnMLup2QmGuuTwIuAM+nm4He37zOTHDjQ9nHgHUkWJVkInEM31TOVp9HNn+8E5ic5B3j6wPZ/A5Ym2Q+gqu6iuz7xJ0menmS/JD+d5OeHONbYeFM9h9P135P8aP+7fz1wRd9+A91U1zOS/BTw5t3VkuQXkjy//4z9/XTTNBO9DjQHDPd2fZfuQt21SR6kC/WvAm8BqKq/Bt5L9yf1/f2244cc+6h+3AeAdcCZ1X+2fZwLgeX9dMQnJxlrHd3FvW9W1Zd3c8x3Ahf3Y/1a/xgeopvaOQz4q8l2rKqv0YX5Lf3+z6b7xMpm4EbgK8D1fdtUNgB/T3fB+ja6T7IMTnN8ov/3niTX98u/SXdh8ya6N8Or6ObVh/FB4JQk9yb50JD7TOWfgW3Ap4D3VdXYl8UuAb5M94mmq/lB6I95D90b4n1J3kp3Zn8VXbBv7cedznUPzaJM769gae/Snzn/+6o6fcrOP+T6qZNvAPvPcIpI+xDn3LXPSvIM4I0M+SkX6YeJ0zLaJyX5bbrpkL+rqmvmuh5pb+O0jCQ1yDN3SWrQnM25L1y4sJYuXTpXh5ekfdJ11133rapaNFW/OQv3pUuXsnnz5rk6vCTtk5LcNnUvp2UkqUmGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalB3vJXkmZg48ZM3WkSK1fO/g0bPXOXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDhgr3JMcluTnJtiRnTbD9jCQ7k9zQ//zW6EuVJA1r/lQdkswD1gC/DOwANiVZV1U3jet6RVWtnoUaJUnTNMyZ+9HAtqq6paoeAdYCJ81uWZKkPTFMuB8MbB9Y39G3jXdykhuTXJVkyUQDJVmVZHOSzTt37pxBuZKkYYzqgurfAEur6gXAPwAXT9Spqi6oqhVVtWLRokUjOrQkabxhwv0OYPBMfHHf9oSquqeqHu5X/xfw4tGUJ0maiWHCfROwLMlhSQ4ATgXWDXZIctDA6onA1tGVKEmarik/LVNVu5KsBjYA84CPVtWWJOcBm6tqHfC7SU4EdgHfBs6YxZolSVOYMtwBqmo9sH5c2zkDy28H3j7a0iRJM+U3VCWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBg0V7kmOS3Jzkm1JztpNv5OTVJIVoytRkjRdU4Z7knnAGuB4YDlwWpLlE/R7GnAmcO2oi5QkTc8wZ+5HA9uq6paqegRYC5w0Qb93Ae8Fvj/C+iRJMzBMuB8MbB9Y39G3PSHJi4AlVfW3I6xNkjRDe3xBNcl+wPuBtwzRd1WSzUk279y5c08PLUmaxDDhfgewZGB9cd825mnA84CNSW4FXgqsm+iialVdUFUrqmrFokWLZl61JGm3hgn3TcCyJIclOQA4FVg3trGqvlNVC6tqaVUtBb4InFhVm2elYknSlKYM96raBawGNgBbgSurakuS85KcONsFSpKmb/4wnapqPbB+XNs5k/RduedlSZL2hN9QlaQGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoOGCvckxyW5Ocm2JGdNsP1NSb6S5IYkn02yfPSlSpKGNWW4J5kHrAGOB5YDp00Q3pdX1fOr6gjgfOD9I69UkjS0Yc7cjwa2VdUtVfUIsBY4abBDVd0/sPpjQI2uREnSdM0fos/BwPaB9R3AS8Z3SvJfgN8DDgB+caKBkqwCVgEccsgh061VkjSkkV1Qrao1VfXTwNuAd0zS54KqWlFVKxYtWjSqQ0uSxhkm3O8AlgysL+7bJrMWeNWeFCVJ2jPDhPsmYFmSw5IcAJwKrBvskGTZwOqvAP86uhIlSdM15Zx7Ve1KshrYAMwDPlpVW5KcB2yuqnXA6iTHAo8C9wKvm82iJUm7N8wFVapqPbB+XNs5A8tnjrguSdIe8BuqktSgoc7c9zYbN2bG+65c6UfwJbXPM3dJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNGirckxyX5OYk25KcNcH230tyU5Ibk3wqyaGjL1WSNKwpwz3JPGANcDywHDgtyfJx3b4ErKiqFwBXAeePulBJ0vCGOXM/GthWVbdU1SPAWuCkwQ5V9emq+l6/+kVg8WjLlCRNxzDhfjCwfWB9R982mTcCfzfRhiSrkmxOsnnnzp3DVylJmpaRXlBNcjqwAvjjibZX1QVVtaKqVixatGiUh5YkDZg/RJ87gCUD64v7tv9HkmOBs4Gfr6qHR1OeJGkmhjlz3wQsS3JYkgOAU4F1gx2SHAl8GDixqu4efZmSpOmYMtyrahewGtgAbAWurKotSc5LcmLf7Y+BBcAnktyQZN0kw0mSngTDTMtQVeuB9ePazhlYPnbEdUmS9oDfUJWkBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBs2f6wKkubBxY2a878qVNcJKpNnhmbskNWiocE9yXJKbk2xLctYE249Jcn2SXUlOGX2ZkqTpmDLck8wD1gDHA8uB05IsH9ftduAM4PJRFyhJmr5h5tyPBrZV1S0ASdYCJwE3jXWoqlv7bY/PQo2SpGkaZlrmYGD7wPqOvm3akqxKsjnJ5p07d85kCEnSEJ7UC6pVdUFVraiqFYsWLXoyDy1JP1SGCfc7gCUD64v7NknSXmqYcN8ELEtyWJIDgFOBdbNbliRpT0wZ7lW1C1gNbAC2AldW1ZYk5yU5ESDJUUl2AK8GPpxky2wWLUnavaG+oVpV64H149rOGVjeRDddI0naC/gNVUlqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDVoqHBPclySm5NsS3LWBNufkuSKfvu1SZaOulBJ0vCmDPck84A1wPHAcuC0JMvHdXsjcG9VHQ78KfDeURcqSRreMGfuRwPbquqWqnoEWAucNK7PScDF/fJVwC8lyejKlCRNx/wh+hwMbB9Y3wG8ZLI+VbUryXeAZwLfGuyUZBWwql99IMnNMyl6CAvHH3ugilk6pBqym9cP+BrSEGbzNXToMJ2GCfeRqaoLgAtm+zhJNlfVitk+jtrk60d7am94DQ0zLXMHsGRgfXHfNmGfJPOBA4F7RlGgJGn6hgn3TcCyJIclOQA4FVg3rs864HX98inAP1VVja5MSdJ0TDkt08+hrwY2APOAj1bVliTnAZurah1wIXBJkm3At+neAObSrE/9qGm+frSn5vw1FE+wJak9fkNVkhpkuEtSg5oK96lukyDtTpKPJrk7yVfnuhbtm5IsSfLpJDcl2ZLkzDmrpZU59/42CV8Hfpnui1abgNOq6qY5LUz7jCTHAA8Af1FVz5vrerTvSXIQcFBVXZ/kacB1wKvmIodaOnMf5jYJ0qSq6hq6T3tJM1JVd1XV9f3yd4GtdN/gf9K1FO4T3SZhTp5UServjnskcO1cHL+lcJekvUKSBcBfAm+uqvvnooaWwn2Y2yRI0qxKsj9dsF9WVX81V3W0FO7D3CZBkmZNf6vzC4GtVfX+uaylmXCvql3A2G0StgJXVtWWua1K+5IkHwe+ADwnyY4kb5zrmrTP+Q/Aa4FfTHJD/3PCXBTSzEchJUk/0MyZuyTpBwx3SWqQ4S5JDTLcJalBhrskNehJ/Q+ypVFL8k66m309Hbimqv7xST7+icDyqvqjJ/O40lT8KKT2aWPhXlXvm+tapL2J0zLa5yQ5O8nXk3wWeE7fdlGSU/rlW5O8p/8CyeYkL0qyIcn/SfKmgXF+P8mmJDcmObdvW5pka5KP9PfjvjrJU/ttv9vfp/vGJGv7tjOS/PnAvv/Ub/9UkkMGavtQks8nuWWsTmk2Ge7apyR5Md2tJY4ATgCOmqTr7VV1BPAZ4CLgFOClwFiIvwJYRner6COAF/f3c6dvX1NVzwXuA07u288CjqyqFwBPvEkM+DPg4n77ZcCHBrYdBLwM+FXAKRzNOsNd+5qXA39dVd/r77Y32f2Dxtq/AlxbVd+tqp3Aw0l+HHhF//Ml4HrgZ+hCHeAbVXVDv3wdsLRfvhG4LMnpwK4JjvlzwOX98iV0YT7mk1X1eP+fNjxr6EcrzZAXVNWqh/t/Hx9YHlufDwR4T1V9eHCn/h7cg/0fA57aL/8KcAzwSuDsJM+fQT30x5ZmlWfu2tdcA7wqyVP7/8bslTMcZwPwhv6+2yQ5OMlPTtY5yX7Akqr6NPA24EBgwbhun6ebMgJ4Dd2UkDQnPHPXPqX/vymvAL4M3E13q+eZjHN1kp8FvtDdpZUHgNPpztQnMg+4NMmBdGfeH6qq+/p9x/xX4GNJfh/YCbx+JrVJo+BHISWpQU7LSFKDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoP8Lxxd2p5BjWzEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7cde4a7208>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "dims = np.arange(len(sens))\n",
    "ax.bar(dims, sens, 0.1, color='y')\n",
    "ax.set_xticks(dims)\n",
    "ax.set_xlabel('dimension')\n",
    "ax.set_title('Sensitivity to latent inputs');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting vs PCA\n",
    "We see that using the 2 more relevant dimensions, the Bayesian GPLVM is able to seperate the\n",
    "three classes while PCA cannot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-19T09:42:12.182381Z",
     "start_time": "2018-06-19T09:42:12.180135Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dim1, dim2 = sens.argsort()[::-1][:2]  # the two dimensions with highest sensitivity"
   ]
  },
  {
   "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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3X98ZHd93/v3R7JstNiIYG+Bu0ZS0jhKCGoSri4JmLR0lzZmWWMgJIU7GANpVfBNH+Zx0wcl0S0NzlXvI6RNvU1CXRF+LDAhNGBjvF4nMVYClxRo1sQwxs5iw10tVggs6yBsVsWy9nP/OHNWI+mc+Xl+z+v5eMxjNGeOZr46mjnzme/38/18zd0FAACAwYzk3QAAAIAqIKgCAABIAEEVAABAAgiqAAAAEkBQBQAAkACCKgAAgAQQVAEAho6Z/ayZnci7HagWgioMzMxOmtm6mT1mZt80s/eb2cXN+37OzD5tZo+a2Wkz+5SZvWzH77/IzNzM/k0+fwGAtO04T/ydmd1hZs/Kqz3u/v+6+0xaj29mrzazz5vZ98zsW82frzcza97/fjN7vHk8HjGzu8zsR5v3/bqZfSjiMT9pZm+P2P7zZrZqZqNm9qHm+fSlO/b5neb216b1N4OgCsm52t0vlvRcSXOS/i8ze5WkP5L0AUmXS3q6pLdLunrH714n6RFJr8uuuQByEJ4nninpm5J+J+f2pMLMfkXSYUm/JekZCs59b5J0paQLW3Z9Z/N4XC7pW5Le3+Ghj0i6NmL7tZI+5O6bzdtfUcv51MzGJL1K0td6/VvQG4IqJMrdVyXdKWlW0m9L+g13/313X3P3c+7+KXf/F+H+ZvZkBW/2/0PSFWY2l0vDAWTG3f+npI9Kena4zcxeamZ/ZWbfNbOvm9mvt9x3h5n9q9bHMLMvmdkrmj//aLOn5xEzO2Fmv9iy30Ezu7/ZW75qZv+6uf1FZvZwy35vM7OvNve7P3zs5n2vN7PPmNl/aPay/X9m9pKov83MJiTdKOl6d/+ouz/qgb9y95q7fz/ieJyV9AeSntPh0N0i6Rlm9oKW57tU0kEFX15DH5f0omZbJOmlko5LOt3h8TEggiokqtmdf1DSWUnPUnDibOeVkh5T0KP1Jwp6rQBUmJntkfTPJH2uZfP3FPSuPFVBEPBmM3t5874jkl7b8vs/IWmfpDuaX8zuUhCU/D1Jr5b0LjMLA7b3SPqX7n6JgqBlOaZZX5X0s5ImJL1D0ofM7Jkt9/+0pBOSLpP0TknvCYfydni+pIsk3dbhMJzXTJeoSfqrdvu5+/cUnFNbe/VfLelL7v7llm3rku6QFAaXr9P2oAspIahCUj5uZt+R9BlJn5J0U3P7Nzr83nWSPtLstv4DSa9udlUDqJ7wPLEm6Z8oGB6TJLn7n7t7o9mj/SVJH5b0j5p3f0LSj5jZFc3b1yo4bzwu6ZCkk+7+Pnd/wt3/StLHJP1Cc98NSc82s6e4+9+5+xeiGubuf+Tuf9N8/o9IelDS81p2WXH3dzfPVUcUDGE+PeKhLpP0bXd/ItxgZv/dzL7TzCn7hy37/uvm8XhI0sWSXt/u4DUdkfSLZnZR8/brmtt2+oCk15nZ0yS9QMExRMoIqpCUl7v7U919yt2vl3Smuf2Zcb/Q7NX6x5LqzU23SXqSgm+pAKrn5e7+VAXv81+W9Ckze4YkmdlPm9mfNSe0rCnIQbpMOj9c+BFJrzWzEUmvkfTB5mNOSfrpZtDynWaQUlOQyyRJP6+g93ylOVHm+VENM7PXmdm9LY/xnPD5m/42/KE5XCcFgdBOZyRdZmYXtOz/gubffUbbP3f/Q/O8+Qx3f5m7f7XNsQt9StJ3JV1tZj8i6acUBKBR+10u6Vcl3RY17IjkEVQhLSckfV3BCS3OtQpeg7eb2d8qSKJ8khgCBCrN3Tfd/RZJm5Je2Nz8Bwp6U57l7hOSbpbUOrx2REGwdEDSWXf/bHP71yV9qhmchJeL3f3Nzef6S3e/RsHQ4Mcl/bed7TGzKUnvVhDoXdoMgO7b8fzd+qyk70u6po/f7cjdXc1eKAXn0GPu/u2Y/eqSfkUM/WWGoAqpaL6h/09J/9bM3mBmTzGzETN7oZktNXe7TkHuwk+2XH5e0sFm8iWACrLANZJ+QNIDzc2XSHrE3f+nmT1P0v/e+jvNIOqcpP+orV4qSTqqYGjwWjMba17+NzP7MTO70MxqZjbh7hsKenjORTTpyZJczURuM3uDOieNR3L37yg4r73LzF5lZpc0z30/2Xyebo2Y2ZNaLhe13PcBSVdJeqOih/5C/0nSP3H3v+j170B/CKqQGnf/qIJk1DdK+hsFU6j/b0m3mdnPKOi2/z13/9uWyycU5Be8Jq92A0jN7Wb2mILgZlHSdS0J1tdLutHMHlVQemVXj5KCYGJW0vkaTu7+qKR/qiBh+28UDNP9poJkcSnozTlpZt9VMKRY2/mg7n6/gmDtswrOU7OS+g5E3P2dCr5UvrX5eN+U9F8l/RtJ/73Lh3mNgoTz8HJ+aNDdH5L0PxT8jXe0accZd7+7jz8BfbKgQwEAgGIzs9dJmnf3F3bcGcgBPVUAgMJrlmG4XtJSp32BvBBUAQAKzcx+TkG+0zcVJLQDhcTwH4Ch1izt8QEFNYdc0pK7H863VQDKiKAKwFBrVs1+prt/wcwukXSPgnpK9+fcNAAlw/AfgKHm7t8Iq2w3Z5I9oGAJFADoyQWdd0neZZdd5tPT03k8NYCc3HPPPd929715t6MdM5tWUKH683H7cP4Chk+3569cgqrp6WkdP348j6cGkBMzW8m7De00F7X9mKS3uPt3d9w3L2lekiYnJzl/AUOm2/MXw38Ahl5zEe+PSao3l0/Zxt2X3H3O3ef27i10ZxuAHBFUARhqZmaS3iPpAXf/7bzbA6C8CKoADLsrFSxlst/M7m1eDubdKADlk0tOFQAUhbt/RpLl3Q4A5UdPFQAAQAIIqgAAABJAUAUAAJAAgioAAIAEEFQBAAAkgKAKAAAgAZRUAACUUkOrWtYJrWldExrXfs1olrWwkSOCKgBA6TS0qqNqaEObkqQ1reuoGpJEYIXcDDz8Z2bPMrM/M7P7zezLZnZDEg0DdqnXpelpaWQkuK7X824RgJws68T5gCq0oU0t60ROLQKS6al6QtKvuPsXzOwSSfeY2V3ufn8Cjw0E6nVpfl46eza4vbIS3JakWi2/dgHIxZrWe9oOZGHgnip3/4a7f6H586OSHpDoe0XCFha2AqrQ2bPBdgBDZ0LjPW0HspDo7D8zm5b0U5I+n+TjAjp1qrftACptv2Y0ptFt28Y0qv2ayalFQIJBlZldLOljkt7i7t+NuH/ezI6b2fHTp08n9bQYFpOTvW0HUGmz2qdDmtW4xs5vu4AqQchZIrP/zGxMQUBVd/dbovZx9yVJS5I0NzfnSTwvhsji4vacKknasyfYjoEwLR1l9oTOnf95XRvMAESukpj9Z5LeI+kBd//twZsERKjVpKUlaWpKMguul5ZIUh9QOC09TO4Np6U3tJpzy4DOmAGIokmip+pKSddKapjZvc1tv+buxxJ4bGBLrUYQlbB2H0p800fRMQMQRTNwUOXun5FkCbQFQMb4UEKZTWg88rXKDEDkhaw+YIgxLR1lxgxAFA1BFZCDRl26aVp6x0hw3cipODwfSiizcAZg+CVgQuM6pFmGrpEb1v4DMtaoS7fPSxvNiYxrK8FtSZrNOGUs/PBh9h/Kalb7eL2iMAiqgIzdvbAVUIU2zgbbsw6qJD6UACApDP8BGVuLKQIftx0AUA4EVSiWel2anpZGRoLrek7JRimaiCkCH7cdAFAODP+hOOr17VXTV1aC21Kl6lMdWNyeUyVJY3uC7WmgYjqqjNc3ioSeKhTHwsL2ZWik4PbCQj7tSclsTbp6SZqYkmTB9dVL6eRTUTEdVcbrG0VDTxWK41RMUlHc9hKbrWWTlE7FdFQZr28UDT1VKI7JmKSiuO3oiIrpqDJe3ygagioUx+KitGfP9m179gTbS6QohT0lKqaj2nh9o2gIqlActZq0tCRNTUlmwfXSUqmS1MPCnmsrknyrsCcV04Hk8fpG0ZBThWKp1UoVRO1UxMKeEhXTUU28vlE0BFXIR70ezOo7dSrImVpcLHUwFSpiYU8qpqPKeH2jSBj+Q/bCelQrK5L7Vj2qChT6pLAnAAwvgipkr8L1qA4sBoU8W6VZ2BMYNg2t6rCWdaPu0GEtU5MKhUJQhexVuB5VloU9gWFDsU8UHTlVyN7kZDDkF7W9ArIq7AnkJa+lYSj2WQIVzZftFkEVsre4uH2NP6mU9aiAYRT2FoXBTdhbFEoz2KLYZ8ENyfqt7TD8h+xVoB4VMKzieov+WF9OfWiOYp8FV+F82W4RVCEftZp08qR07lxwTUAFlEJcr9C6NmKH5pJCsc+Cq3C+bLcY/gMAdG1C4z0NtyU5NNep2GdeuV5oqni+bDcIqgAAXduvmW05VVLQW3SBRrSujV37Jz00F1fss12uF4FVRsiXZfgPANC9We3TIc2eD5YmNK5DmtVV+vFch+bazQxERsiXpacKANCbdkvD5DX8xszAgij5+q2DIqgCACQiz3X44nK9mBmYMepUAQCQjqySx+NyvZgZmCHqVJFThQqq16XpaWlkJLiuwELNQBlluaxMXK4XSeoZok4VPVWoGL4pAV1Luxcp62Vl8hx+hKhTJXqqUDV8UwK6kkUvEsnjQyauHtUQ1akiqEK18E0pMQ2t6rCWdaPu0GEtpzJkg/xkUYKAZWWGzOJiUJeq1ZDVqWL4D9VCRd9EUEix+rLoRUojeXyQIUsqrqesVpP+4i+C2lSbm9LoqHTddUOVekFQhWqpSEXfvE/+WefCIH07X1PjGku9AnqnZWV6NUiwzxeFDNTr0pEjQUAlBddHjkhXXjk0gRVBFaolfOOWuE5KEU7+5MJUS9RralSmEZnOyc/vl0YJgiSTxwcJ9vmikIF2Oa0lOgcPgqAK1VPyir5FOPlTSLFaol5Tm3KNa0wX6oLSDIcNEuzzRSED5LQSVAFFU4STP4UUqyXutbOuDV2lHy90INVqkGCfLwoZIKeV2X9A0RRhxhSFFKul3WsnrWKcadivmb4XbR7kd9GlqNl/ZtLBg/m0Jwf0VAEFU5ReIgopVkfUaypUpryiQRLfk06aR4Rw9t/NN0vezNVzH6pkdYIqoGA4+WMQ7WaO3qp7I3+nTHlFgwT7fFHIwLFjWwFVaIiS1QmqgALi5I9+dJo5GgZbO3U7tJx3qQ+UwJAnq5NTBQAV0alK+iB5RVkujowSG/KlauipAoCK6DRzdJCh5SKU+pDoLSu8ihRg7hdBFQBURDdlA/odWi5CqY8iFMZFBxUowDwIgioAqIgrtFfHtTt35Qrt7er32/UCFaHOU1F6y9BByQswD4KcKgBDzczea2bfMrP78m7LoB7U6cjt9+jrulF36LCWY3OgOuVMFaHOUxF6y4B2CKqADDXq0k3T0jtGgutGPe8WQdL7JV2VdyOSEBdceHN9v3bJ5Z2S3ItQELYIhXGBdhj+y0u9PrRjzsOqUZdun5c2mvmbayvBbUma5V+fG3f/tJlN592OJMQN0bXa0KY+ri/qVt27bYivm16gvEt9FKUwLhCHoCoP9fr22RErK8FticCqwu5e2AqoQhtng+0EVUhCu8rprXb2XEntc6aKMuOOwrgoOoKqPCwsbJ9uKg1VxdlhtRZT+y5uO4rDzOYlzUvSZIHr7ewMOkx2PoCKEw7xxfUCXaG9hZpxl3dvGdAOOVV5GPKKs1XRa37URMxncdx2FIe7L7n7nLvP7d3b3Uy6vMxqn27Qfr1dL9XL9RNd/c6a1mNzph7U6ba5VgC20FOVh8nJYMgvajtKoZ/8qAOL239Hksb2BNuBQcUN0cWt99cqDKSieoGKvF5gUYYlgRA9VXlYXAwqzLYaooqzVdAuPyrObE26ekmamJJkwfXVS+RT5c3MPizps5JmzOxhM/ulvNvUq3blEDrNjOuU6F3UGXcsm4MioqcqD0NecbYK+s2Pmq0RRBWNu78m7zYMql05hHbJ69307hR1xh2FQFFEBFV5GeKKs1UwMRkM+UVtB7LWrhzCoDPmijrjjkKgKCKCKqAP5EehSDotITPojLkizrgrwrI5wE7kVAF9ID8KRVKEJWSyNox/M4qPniqgT2XOj2LWVLUUdYguTcP4N6P4CKqAIRPOmipKMUckI8khurIE3XF/c1naj+ohqAKGzB/ry8yaQqyyB939tJ8gDEkhpwoYIg2tal0bkfcxawpS+1IFZdBr+6l3hSTRUwUMkXYfjMyaGh7tembKXqqg1/ZT7wpJoqcKGCLtPhiZNTUcOvXMFLWCerd6bX/Zg0gUSyJBlZm918y+ZWb3JfF4ANIR98EyrjG+lQ+JTsNjZS9V0Gv7yx5EoliS6ql6v6SrEnosACmJ+8C5Sj+eU4uQtU49M7Pap0OaPR9UTGhchzRbmqC71/aXPYhEsSSSU+Xunzaz6SQeC0B6qO2DbiqRF7GCei96aT/vCSQps0R1M5uXNC9Jk5MskLZLvc4CyzgvzSneZf/AxGCKukBynnhPICmZBVXuviRpSZLm5uY8q+cthXpdmp+XzjYXkltZCW5LBFZDqOx1glBs9MwA6aGkQhEsLGwFVKGzZ4PtBFVDhyneSBs9M0A6KKlQBKdO9bYdlcYUbwAop6RKKnxY0mclzZjZw2b2S0k8bqXV69L0tDQyElyikHs2lJjiDRRXQ6s6rGXdqDt0WMtUXsc2Sc3+e00SjzM0duZQbW7u3mfPniBZHUOHRGIgff1MBiHfEZ0w/JeHqBwqSRodlcykqSlpaYl8qiFV9jpBQNH1u95f2ddFRPpIVM9DXK7UuXPBBUOPRGIgPf1OBiHfEZ3QU5WHuFypycntuVbT08FtAEBi+g2OyHdEJwRVeVhcDHKmWu3ZIx08GORaraxI7lv1qgisACAx/QZHLGmDTgiq8lCrBTlTU1Pbc6iOHYuvVwUASES/wRH5juiEnKq81Gq7E9GvvTZ6X+pVAUBiBqkqT74j2iGoKpLJyWDIL2o7ACAxBEdIA8N/RRKXa0W9KgAACo+gqkjicq2oVwUAQOEx/Fc0UblWqIRGXbp7QVo7JU1MSgcWpVn+1QBQGQRVQAaOfWZVx3/2hPxr69Kpca392oxunw/yOQisAKAaGP4DUtbQqo4/tyFNrstGJJtel97d0MY1q7qbahkAUBkEVUDKlnVC2rN9SQx78qb0709ojWoZAFAZDP8BKYtd+mJyXRNdVstoaLWvmjoAgOzQUwWkLG7pC3t4XAe6qJbR0KqOqnE+OFvTuo6qoYZWk2wmAGBA9FQBKduvGR1VQxtqGQI8O6q5UzNdJakv68T235W0oU0t60RsbxU9WwCQPXqqgIQ16tJN09I7RoJr1XevF/aKPbM6+MLugpy44cO47fRsAUA+6KkCEtSoS7fPSxvNdbHXVoLbV2ufbqj111M0ofHIACpuWLGfni0AwODoqQISdPfCVkAV2jirgUon7NeMxjS6bduYRrVfM5H799qzBQBIBj1VQILiSiQMUjoh7F3qNkeq154tAEAyCKqABE1MBkN+UdsHMat9XQ/dRSXGt+vZAgAkg+E/IEEHFqWxPdu3je1RV6UTkjKr3YnxhzRLPhUApIyeKiBBYYmEbhZOTnOB5V56tgAAySCoAhI2W+scHMXNEgx/HwBQPgz/ATlIY5YgACBfBFVADtKYJQgAyBdBFdCHnVXTG/Xefj9uNuCgswQBAPkhqAJ6FOZDra1I8q18qF4CqyLMEgQAJIugCuhREvlQszXp6iVpYkqSBddXL5GkDgBlxuw/oEdJ5UN1M0tQChZI7raaOgAgP/RUAT3KMh+qoVUdVeP8sjNrWtdRNdTQavJPBgAYCEEV0KNe86EGSWpf1olty81I0oY2tawTPbUZAJA+hv+AHvVaNX2QIp9RCyO32w4AyA9BFdCHbvOh2iW1d/P7ExqPDKDCdf0AAMVR3eG/el2anpZGRoLreo+FhIAEDJrUvl8zGtPotm1jGtV+zQzYMgBA0qrZU1WvS/Pz0tlmF8HKSnBbkmrMWUd2Jiab9awitncjnOXH7D8AKL5qBlULC1sBVejs2WA7QRVS1loCYfy+cY38qxmde/9WENRrkc9Z7StUEEWJBwC5qNeDz/FTp6TJSWlxsXCf6dUIqnYe6JWIrgEpuB9IUVgCIZyxt37xukbf3dBFl0jrv7uvbVJ71GMVLXjZ+feFJR4k5d42ABVWkhGo8udUhQd6ZUVyD67NovednNz6HfKtkIKoEgibF2zqov98Qv/unPSWk90HVEWsT0WJBwC5aDcCVSDlD6qiDrT77sBqz56gqzAqCJufJ7BCIpIqgVDU4IUSDxhGgy6gjgTEjTRFbc+x46T8QVXcgXaXRpuzpqampKWloIuwJNEuyimu1EGvJRCKGrwk9fcBZZHEAupIwGTM7J6d23PuOCl/UBV3oCVpc3Orhyocc+0l2gV6lFQJhLyCl4ZWdVjLulF36LCWdw03UuIBwyaJBdSRgMXF4PO8Vfj53iqu4+S66zIJrMofVEUd6FY7e6G6jXaBPsxqnw5p9nzwM6FxHdJsz0nceQQv3eRxJfX3AWWR1ALqGFCtFow4TU0F6T2tI1Ct4jpINjcz6bEqf1DVeqDjtB7kbqNdoE+z2qcbtF+vqL9UPr1ft4zs6zkPI4/gpds8rvDve7teqhu0vxIBlZldZWYnzOwhM3tb3u1BcWS5gDo6qNWkkyelc+eC66hZf+06SDJI9Sl/UCVtHei4wKr1IHcb7QIDSCIPI+vgpah5XGkzs1FJvyfpJZKeLek1ZvbsfFuFouh1AXXkrNPoVcqpPtUIqkLd9kJ1E+0CAyhjHsYQJ6E/T9JD7v41d39c0h9KuibnNqEgZmvS1UvSxJQkC66vXuquNAoyFM74u/ZaaXy8c2mllFQrqOqlF4paVUhRGfMwhjgJfZ+kr7fcfri57Twzmzez42Z2/PTp05k2DvmbrQU15trVmqPsQo52zvg7c0YaGwsurTJI9alWUCV11wtFrSqkLC7fYvxp2bajFyShx3P3JXefc/e5vXv35t0cFAxlFzIS1xkSNePv8celpzwl81Sf6gVV3aBWFVJ2YFEaGdu9/fFHi32irWISehdWJT2r5fblzW1AV8o43F867TpD4vKkHnlkeyeLlPoI1XAGVdSqQspma9JFT9m9ffPxZE60nepJoSd/KekKM/tBM7tQ0qslfSLnNqFEyjjcXzrtOkO6KZWU0QhVtYOqsKvQTLrgguB6elp6WswYDLWq0KeofIr1R6L3HfREW9R1AcvK3Z+Q9MuS/kTSA5L+m7t/Od9WJYdcn/RRdiFl9XoQBEVZWeluklpGI1TVDapao1IpKPwlBbe/+13pwgu370+tKvQpLp8iLn9q0BNtUdcFLDN3P+buP+Luf9/dK3MiINcnG5RdSFH4WR4nnOXXaZJaRiNUFyT6aEUSFZWGNjakSy+VLr44OKCTk9uXsgF6EJdPIQtOrK339XuibWhVyzrRtmZU1etJVUWjHrxm1k4FAfaBxfSm57fL9SlCSYAsj0WawjZX4W8pnHaf5VIwlLew0Lk80uRkdG9XwiNU1Q2qOkWfjzwiffvb2bQFlRY3nLfxPWnuzdKDxwY70YbDfTt7p3YagnpSpRf2HIWBTthzJA3+ARwVoCSd65NkEJTmscjDbK2c7S68uGG/Vt30Ni0uBj1erQFaCiNU1R3+6xR9Rt1P7Sr0od1w3oPHOte36SRquG+nEdkw1JMqvbRmiWUxBJ30UCIz5tCV0dHO+4yMdP68zmg1leoGVe1K1UdFp9SuQp/aDeclMfunm2G9mNrBKJi0ZonFDkEruVyfpIMgZszhvLhJZfX6Vj50O90ulpzBairVDap2LrQcRrujo1sZ/63/AGpXoU+zNWn80uj7bGTwWVfdDOttytsmqlOCIR87Z96lNXkhLhBZfyS5JVaSDoKYMQdJ7SeVveEN3T9OQT6vqxtUSVtRqbt05EjQQ9X6D2uNbKldhQG85PDuHgFJ8k0NPFQStXxMlLgeLUow5CNquOzxR3cXhU1illi7AKWbJVYGfY5+FGHGHOUmCqDTpLJeFODzutpBVatOPVHdFA8DYuxcdNUiYqB+h0p2Lh9jMYN9cT1alGDIR9Rw2ebjQVHYpBfnjQtQrjiYXNCQdBCU90LFlJsoiCQDoQJ8Xld39t9OnXqiMpoZgOpqnf3zjpivK/0Olcxq3/klY6JmA7Zb+DiuB4sSDOlqNyT31oQnHkdN6b/ioPTFI8nNrkujbECeM+aKXm5iaMSVOuhkbGx7T1ZBPq+Hp6cqLoINq6tnNDMAwyHNfJFeFz6O68EaV8TihEhM1jlDO4f5Hjw2eGL5zuExafdQYlmH0EiUL4h2k8rGxnYX6jaT3vxm6X3vK+TndSJBlZldZWYnzOwhM3tbEo+ZuMXF3f8cKaiuHuZVZTAzAMMh7XyRXhY+3q8ZjUYMGX5fT5BXlaK8c4YGDRq6GR4rwxBaXNBHonxBxE0qm5oKAqf3vjco1h162tOkK68s7Of1wEGVmY1K+j1JL5H0bEmvMbNnD/q45yVVO6pWky65ZPf2jY1CzBhAteSdL7KtLdqnCyNG+s91mDGIweT9Gug1aNgZfNx5Q+eerqLXmmoX9OUd9CKCmXT55dKHPrQ9UFpvSVU4c6bQ5Y6SyKl6nqSH3P1rkmRmfyjpGkn3D/zI4VTLMM8pnLEn9ReVPhKzwm0BZgygeqLyRfJalmNd0bNoyKtKV545QwcWt1csl7aChp2vw6j8qzit9xV9CK1d0PeWk1v7sLRMjjp9zrebZFaQ3qlWSQz/7ZP09ZbbDze3Da7TjL1ee7GY4Ycc5TlUEpdXxdI21RXXUybtfh0ev3l38BHLyjOE1inoS6rcBAbQ6XO+ZOWOMktUN7N5MztuZsdPnz7d3S+1O5j9VECPSogryIwBVF+eQyVRta7azRhENUQFDVGvQ3kPD+pbr9miD6EVPeiDOgdN7TpDCri0XBJB1aqkZ7Xcvry5bRt3X3L3OXef27scOMYaAAAa+ElEQVR3b3eP3O5g9lMBnRl+yFhrnkrckEoWQyW9zhhEdSWydNJK8Lq+5VrpgvHmigIJ540lMauw6EEf1HkEKa4z5ODBQi4tl0RQ9ZeSrjCzHzSzCyW9WtInEnjc9j1LcdHtykr7yLWgMwZQPTuH++Jk9a25lxmDqK7Y19uOCaJje4JlluKEr+v1M9IT69IrP5jcEFpSQ+V5TxZAFzqNIMV1hhw7Vsil5cy9l37fmAcxOyjpJkmjkt7r7m2/B8zNzfnx48e7e/B6PThIp04FkeviYnCQp6ejC4aZBVFraM8eeqOQi5um2yf8SsEHV9xJvqFVLeuE1rSuCY1rv2ZKHQiZ2T3uPpd3OwbV0/krR3GTIsKAZWcC++XPl07+ebC0ko1K/+u8dPy/dP98E1Nbyd+9tmmnuPdON8+BEor7nG9nZGT7Z33ILOg0SVi3569Ecqrc/Zi7/4i7//1OAVXP4nqWoqLbnQGVFESuN9xQuHFXVF/bYZYO35pZrw+DaNfTE9V78xPXSQ9/trlWpYLrLx6JXyg8ytpK+56kXnqfij6rEAnrZwSpoBPPyltRPapLMK7X7cyZwo27ovpik2SnOs82Yr0+DKLTpIhuq69L0QuFx2k3RNfLRA0SzNFRQSeelTeoknZHt2FF1k4KMO6K6hskSbbf9foaWtVhLetG3aHDWqZna0j12tPTbp3Cq5eiFwiP0m42ay/PTYI5OiroxLNyB1U7tVtDaKeC1rhAdQySJNtPXSmGDBHqtaen3fbZmvSKI933WMUFT708Nwnm6EoBJ54lUVG9OMID2prw9thjwfDfThT8RAb6rai9XzM6qsa2IcBOdaXaDRmWOcEdvWtXTb2f/cPX8K3XbeVdxYkLnnptU57V6IF+laenqtsiXzsj18OHCznuCrTTT12pfocMUT299vR0s/9sTfIOk6pGxqTHH4uuLUXvE4ZBOXqqBlkDMKr3qpvpmsCABl3nb1b7euphmtB4ZADFUjTDqdeenm72n5iMLxMy9mRp43tB3Sppa3Zf+Nj9tAkom3L0VPVTPb1VAcddUW15rPPHUjRIW1wC+dybo9cOzGoZJqAoyhFUlWxBRQyvcGmNW16b/Tp/LEWDtMUN4T14TLGrBhSttlQSy98Accox/Dc5GV09nWRzFEijLt32Rmnz8fh90v6A6XXIEMOr3+HpqCG8W66N379ItaV2VpOPGqIEBlGOnqqCFvkCWt15Q/uASirWBwyGV9LD0+3WEyxSbaleCpAC/ShHUFXQIl9AqFHfStCNQ/FCFEXSwUVUrpVMmntTsj1Agw7dsfwN0laO4T8pCKAIolBA4bf+diamep/9B6Ql6eAifF0PMtu1kySG7uJmL9KDjKSUJ6gCCirqW3+r8UuD9dWAokgjuEi7XEK73rVun7fXAqRAr8ox/Ndt4U8gB+2+3Y+MSS85nF1bgG6UcW29JHrXKECKtBU/qAoLf66sSO5bhT8JrFAQcd/ubVR6+fs4YaN4kggusi5N0Ou6hVEGLcgLdFL84b+4wp/XXRf8TJ4VchY3pJDEN+CGVrWsE1rTuiY0rv2aoWQCEjHIcF0v+U1JBTKDDt1RTgFZKH5PVVyBz81NeqxQCGkNKTS0qqNqnF96Zk3rOqqGGlodvNHAALqdPZhk6YZB32eUU0AWit9TFVf4U9paqobeKuQsjSTdZZ3Qhja3bdvQppZ1gt4q5Co2v2klGAoMe6OSSC5vNcj7jHIKyELxe6qiCn+2YqkaVFTU4sjttgNZaZfH1NobVaRAJomcLKCT4gdVYeHP0dHo+1mqBhUVruHX7XYgK5HFPluEvVFFCmTKOOMR5VP8oEoKAqsjR1iqBqXRy8youH33a0Zj2v5lYkyj2q+ZlFoNdGdbflOMtVPFCmQop4AsFD+nKhTmTS0sBEN+k5NBQEU+FQqm15lR8fsGeVPM/kMRhflNN03HFxLNotJ6L9IuUAqYu2f+pHNzc378+PHMnxfIQuyHzNTuyuq97Ft2ZnaPu8/l3Y5Bcf7abucXAym5kiJAUXR7/irH8B9QIlFBUtz2IiXyAv1gWA3YUp7hP6AkbFTyzejtO7HAK6qAYTUgQE8VkLCogCpue5ESeQEAg6GnCkjYxFR8npS0Y+mZ2rh+cmpGX3ntvkIk8gIA+kdQBSSs3Rpl4dIzYaX0Na3riy9s6NBJMasPAEqO4T8gYe0Sd+OWnvnjx050XdcKAFBM9FQBKYhL3I1bYmZ9z7rONocM29W1AqqkUS9ODSsgCfRUARmKW2LGT23fHi7zAVRVWN9qbUWSb18zECgrgiogQ1FLz/j3RqVf2730DLWqUGV3L2zPO5T4MoHyY/gPyFCYjN669Mz3f3VG6x/enaROrSpUGYVvUUUEVcCAes0LmdW+bTP9Gj8t3f6e6NmCQFV1U/iWnCuUDcN/wADi8kLuuF5dz+ZjmQ8Mo06Fb0uVc1WvS9PT0shIcF0vYiORBXqqgAHE5YUcv1lSc63ybmbzscwHhk34eo/riWqXc1Wo90q9Ls3PS2ebjV1ZCW5LUq1IDUUWCKqAAcTmf/j2m4N8GGyrwK5x7dcMhUJRCe2+TJQm52phYSugCp09G2wnqBo6DP8BA+glmbyfD4OwAntY32pN6zqqhhpa7f3BgBKJe28VbgLHqZg3dtx2VBpBFTCAqLwQWfS+/XwYxFVgX9aJ3h8MKJHSLDY+GfPGjtuOSiOoAgYQlWQ+96bkPgziKrDHbQeqojQTOBYXpT073vB79gTbMXTIqQIGFJUXMnll/1PBW3OoTCbfmaCl+MrsQJWUYgJHmDe1sBAM+U1OBgEV+VRDiaAKSEG/HwZhDlU45BcVUI1pVPu1uwI7gJzUagRRkMTwH1AoUTlUrSY0rkOaZfZfQszsF8zsy2Z2zszm8m4PgHIjqAIKpF2uVNhDRUCVqPskvVLSp/NuCIDyY/gPKJAJjccGVuGsP4Kq5Lj7A5JkFjNlE4VCzTYUHT1VQIHs14zGNBp7P7P+MKyo2YYyIKgCCmRW+3RIs7KYYlfM+uudmX3SzO6LuFzTw2PMm9lxMzt++vTpNJuLGNRsQxkw/AcUTDic0ToLUGLWX7/c/cUJPMaSpCVJmpub2z0lE6mjZhvKgKAKKKAwsCJ/BAjE5RvSe4siIagCCmpW+wiiUmZmr5D0O5L2SrrDzO5195/LuVmIsF8z9N6i8AiqAAwtd79V0q15twOd0XuLMiCoAgCUAr23KDpm/wEAACSAoAoAACABBFVAnxp16aZp6R0jwXWjnneLAAB5IqgC+tCoS7fPS2srkjy4vuW10jsvI7gCgGFFUAX04e4FaePs7u3rZ4Jgi8AKAIYPQRXQh7VT8fdtnA2CLgDAcCGoQm/qdWl6WhoZCa7rw9klMzHZ/v52QRcAoJoIqtC9el2an5dWViT34Hp+figDqwOL0tie+Ps7BV0AgOohqEL3FhakszsSic6eDbYPmdmadPWSNH7p7vvG9gRBFwBguFBRHd07FTOmFbe9Qhr1IE9q7VTQC3VgMQisZmvR96m2qsMspwEAQ2WgoMrMfkHSr0v6MUnPc/fjSTQKBTU5GQz5RW2vsLB8Qjjbb20luC1tBVaztZb9tbpt4dc1reuoGsH+BFYAUFmDDv/dJ+mVkj6dQFtQdAcPSmbbt+3ZIy1We6wrqnxCuxl+yzpxPqA6v782tawTKbUQAFAEA/VUufsDkmQ7P2hRPfW6dORIkKAeMpOuu06q1eJ/rwLiZvLFbtd6T9sBANWQWaK6mc2b2XEzO3769OmsnhZJiUpSd5eOHdu6XdFyC3Ez+WK3a7yn7QCAaugYVJnZJ83svojLNb08kbsvufucu8/t3bu3/xYjH52S1CtcbiGqfEK7GX77NaMxjW7fX6Par5mUWggAKIKOQZW7v9jdnxNxuS2LBqIg4pLRw+0VLrcQlk+YmJJkwfXVS9uT00ONunT39D49XpuVPTwuedBDdUizJKkDQMVRUgHdWVwMep5aAyezIHldqny5hZ0z/KJsmyW4sk/n/mCfxvZI+2MCMABAtQyUU2VmrzCzhyU9X9IdZvYnyTQLhVOrBUnprZMS3IPk9Xq9c0/WEOh1liAAoFoGCqrc/VZ3v9zdL3L3p7v7zyXVMBTQsWPbZ/9JW0N8i4tBeYVWQ1BuoVWvswQBANXCMjXoXrshvlpNWlqSpqaC3qypqeB2xcsttOp1liAAoFoIqtC9TkN8tZp08qR07lxwPUQBldT7LEEAQLUQVKF7DPHFamhVy7VlPfHYHRr5+rL0mlXZ6FZOVaP8lSUAAB0QVKF7DPFFCtf6W9O6ZJJfvi57d0P+i6uSgrUCb3sjgRUAVB1BFXoz5EN8UaLW+tOTN6V/v7XW3+bj0p03ZNwwAECmCKqAAcWu6Te5ffv6mQwaAwDIDUEVMKDYNf1OsdYfAAwTgipgQFFr/fn3RqVf277W3/ilWbYKAPpUr0vT09LISHBdgTVcs8IyNcCAwjX9lnVCa1rX+GPjWn/zjPzDW2v9jYxJLzmcVwsBoEvXXy/dfPNWoeeVlWCJMokc2i4QVAEJmNW+rQWTL5YaV0l33rmVR3XRU/JrGwB0pV7fHlCFwpUzCKo6YvgPSMkTLXnq62eCxZYpqwCgsBYWdgdUobgVNbANQRWQAhZXBlA67QKnuBU1sA1BFZCCgRZXJkkUQB7iAiczVs7oEkEVkIK+F1eu14Ok0JWVoBs+TBIlsAKQtqilyMykN72JfKouEVShe/SgdK3vxZUXFoKk0FZhkigApClqKbIPflB617vybllpMPsP3Ql7UMIPfKbZtjXbPCR3LwRDfhOTQUA12+lQxeU0kCQKIAu1Guf0AdBThe7Qg9Kz2Zr0lpPSvzsXXHcMqKT4nAaSRAGkiZGIRBBUoTv0oGQjKqdhzx6SRAGkh1zOxBBUoTv0oGQjKqdhaYnueADpYSQiMQRV6A49KNmp1aSTJ6Vz54JrAioAaWIkIjEEVegOPSgAUE2MRCSGoArdowcFAKqHkYjEEFQBADDMGIlIDHWqAAAYdtSnSgQ9VQAAAAkgqAIAYJhQ6DM1DP8BADAsWHIsVQRVQIoaWtWyTmhN65rQuPZrRrPal3ezAAyrdoU+CaoGRlAFpKShVR1VQxvalCStaV1H1ZAkAisA+aDQZ6rIqQJSsqwT5wOq0IY2tawTObUIwNCj0GeqCKqAlKxpvaftAJA6Cn2miqAKSMmExnvaDgCpo9BnqsipAlKyXzPbcqokaUyj2q+ZHFsFYOhR6DM1BFVASsJkdGb/AcBwIKgCUjSrfQRRADAkyKkCAABIAEEVAADDiiVrEsXwHwAAw4glaxJHTxUAAMOo3ZI16AtBFQAAZZLUkB1L1iSOoAoAgLIIh+xWViT3rSG7fgIrlqxJHEEVgKFlZr9lZn9tZl8ys1vN7Kl5twloK8khO5asSRxBFYBhdpek57j7P5D0FUm/mnN7gPaSHLJjyZrEMfsPwNBy9z9tufk5Sa/Kqy1AVyYngyG/qO39YMmaRNFTBQCBN0q6M+oOM5s3s+Nmdvz06dMZNwtowZBdoRFUoVxiZr00tKrDWtaNukOHtayGVnNtJorDzD5pZvdFXK5p2WdB0hOSIrN93X3J3efcfW7v3r1ZNR3YjSG7QmP4D+URU6julP5OR2s/rA1tSpLWtK6jakgS6+5B7v7idveb2eslHZJ0wN09k0YBg2DIrrDoqUJ5xMx6eerCb5wPqEIb2tSyTmTYOJSRmV0l6a2SXubuZzvtDwDtEFShPGJmt1xyKjrHZU3rabYG1fC7ki6RdJeZ3WtmN+fdIADlxfAfyiNm1sujk9E5LhMaT7tFKDl3/+G82wCgOuipQnnEzHr5zuK/1ZhGt20e06j2aybDxgEAhh1BFcojZtbLZO2XdUiz53umJjSuQ5olSR0AkCmG/1AuMbNeZrWvUEFUoy7dvSCtnZImJqUDi9Isk3UAoNIIqoCENerS7fPSRnMu2dpKcFsisAKAKmP4D4gTU2i0k7sXtgKq0MbZYDsAoLroqQKixBQaldSx6N5azLqmcdsBANVATxUQJabQqBY6dzdNxKxrGrcdAFANBFVAlJhCo7HbWxxYlMZ2VH4Y2xNsB3rS5xA0gHwQVAFRJmO6leK2t5itSVcvSRNTkkwav1S6YFy65VrppukgkR3oKByCXlmR3LeGoAmsgMIiqAKixBQa1WJ33U2zNektJ6VXflB6Yl1aPyPJt2YCElihowGGoAHkg6AKiBJTaLTXleGZCYi+DTAEDSAfA83+M7PfknS1pMclfVXSG9z9O0k0DMhdTKHRXjATEH2LWeuymyFoAPkYtKfqLknPcfd/IOkrkn518CYB1cFMQPRtwCFoANkbKKhy9z919yeaNz8n6fLBmwRUBzMB0beEhqABZCfJ4p9vlPSRBB8PKL1wWRrWAURfEhiCBpCdjkGVmX1S0jMi7lpw99ua+yxIekJS7JwmM5uXNC9Jk+QEYIjM1giikJF6PZgdeOpUkHu1uEhQBmSoY1Dl7i9ud7+ZvV7SIUkH3N3bPM6SpCVJmpubi90PANCHAZZWApCMgXKqzOwqSW+V9DJ3P9tpfwBASqhrBeRu0Nl/vyvpEkl3mdm9ZnZzAm0CAPSKulZA7gZKVHf3H06qIQCAHtTr0g03SGfOBLdHRoLlbHYihxXITJKz/wAAWajXpTe8QdrY2Np27tzu/ahrBWSKZWoAoGwWFrYHVK1GWk7r4+PZtAeAJIIqIHONunTTtPSOkeCaxZXRs3Z5Uq09VmfOBDMA67zIgCwQVAEZatSl2+eltRVJHlzfPk9ghR71kifFDEAgMwRVQIbuXpA2dsx63zgbbAe6trgojY11vz8zAIFMEFQBGVqL+WyL2w5EqtWk971PuvTSrW2XXrr9ditmAAKZIKgCMjQR89kWtx2IVatJ3/52UEbBPfj58OFgxl8rZgACmSGoAjJ0YFEa2/GZN7Yn2A4MrFaTlpakqSnJLLheWmKZGiAj1KkCMhQurHz3QjDkNzEZBFQsuIzE1GoEUUBOCKqAjM3WCKIAoIoY/gMAAEgAQRUAAEACCKoAAAASQFAFAACQAIIqAACABBBUAQAAJICgCgAAIAEEVQAAAAkgqAIAAEgAQRUAAEACCKoAAAASYO6e/ZOanZa0IukySd/OvAHdKWrbaFfvitq2orZLSqdtU+6+N+HHzFzL+auTIvx/i9AGiXYUrQ0S7ei1DV2dv3IJqs4/udlxd5/LrQFtFLVttKt3RW1bUdslFbttZVGEY1iENtCO4rWBdqTXBob/AAAAEkBQBQAAkIC8g6qlnJ+/naK2jXb1rqhtK2q7pGK3rSyKcAyL0AaJdrQqQhsk2tEqsTbkmlMFAABQFXn3VAEAAFRCpkGVmf2Wmf21mX3JzG41s6fG7HeVmZ0ws4fM7G0Zte0XzOzLZnbOzGJnAZjZSTNrmNm9Zna8QO3K9JiZ2dPM7C4ze7B5/QMx+202j9W9ZvaJlNvU9hiY2UVm9pHm/Z83s+k029NDu15vZqdbjtM/z6hd7zWzb5nZfTH3m5n952a7v2Rmz82iXWVVhHNIUc4XeZ8finAuKML7vijv8S7a8SIzW2s5Fm9PoQ3PMrM/M7P7m++RGyL2Gfx4uHtmF0n/VNIFzZ9/U9JvRuwzKumrkn5I0oWSvijp2Rm07cckzUj6c0lzbfY7KemyDI9Zx3blccwkvVPS25o/vy3qf9m877GMjlPHYyDpekk3N39+taSPFKRdr5f0u1m9plqe9x9Keq6k+2LuPyjpTkkm6WckfT7rNpbpUoRzSFHOF3meH4pwLijK+74o7/Eu2vEiSUdTPhbPlPTc5s+XSPpKxP9k4OORaU+Vu/+puz/RvPk5SZdH7PY8SQ+5+9fc/XFJfyjpmgza9oC7n0j7eXrVZbvyOGbXSDrS/PmIpJen/HyddHMMWtv8UUkHzMwK0K5cuPunJT3SZpdrJH3AA5+T9FQze2Y2rSufIpxDCnS+yPP8UIRzQSHe90V5j3fRjtS5+zfc/QvNnx+V9ICkfTt2G/h45JlT9UYFEeFO+yR9veX2w9r9h+fJJf2pmd1jZvN5N6Ypj2P2dHf/RvPnv5X09Jj9nmRmx83sc2aW5om1m2Nwfp9mcL8m6dIU29RtuyTp55vdzR81s2el3KZuFf29WFZ5n0Oy+L/meX4owrmgLO/7Ir3Hn29mXzSzO83sx9N8ouZw709J+vyOuwY+HhcM0rAoZvZJSc+IuGvB3W9r7rMg6QlJ9aSff9C2deGF7r5qZn9P0l1m9tfNKDzvdiWuXbtab7i7m1ncNNKp5vH6IUnLZtZw968m3daSu13Sh939+2b2LxV8g96fc5sQoQjnkKKcLzg/DIz3/ZYvKHgtPGZmByV9XNIVaTyRmV0s6WOS3uLu30368RMPqtz9xe3uN7PXSzok6YA3BzF3WJXUGrFf3tyWetu6fIzV5vW3zOxWBd28AwVVCbQrlWPWrl1m9k0ze6a7f6PZPfqtmMcIj9fXzOzPFXw7SOOk2c0xCPd52MwukDQh6UwKbempXe7e2obfV5CPUgSpvRfLqgjnkKKcLwp8fijCuaAs7/tCvMdbgxt3P2Zm7zKzy9w90TUBzWxMQUBVd/dbInYZ+HhkPfvvKklvlfQydz8bs9tfSrrCzH7QzC5UkESY6qyxbpnZk83skvBnBYn3kbMZMpbHMfuEpOuaP18nadc3ZDP7ATO7qPnzZZKulHR/Su3p5hi0tvlVkpZjAvtM27VjzP5lCsb6i+ATkl7XnBHzM5LWWoZ00IeCnEOyOF/keX4owrmgLO/7QrzHzewZYU6bmT1PQWyS6Bfe5uO/R9ID7v7bMbsNfjx6zWwf5CLpIQXjlfc2L+Hsi/9F0rGW/Q4qyMz/qoIu7Sza9goF46ffl/RNSX+ys20KZnJ8sXn5chZt66ZdeRwzBfkHd0t6UNInJT2tuX1O0u83f36BpEbzeDUk/VLKbdp1DCTdqCCIl6QnSfqj5uvwf0j6oYxeW53a9f80X09flPRnkn40o3Z9WNI3JG00X2O/JOlNkt7UvN8k/V6z3Q21mdHGpRjnkKKcL/I+PxThXFCE931R3uNdtOOXW47F5yS9IIU2vFBBPuOXtBWDHEz6eFBRHQAAIAFUVAcAAEgAQRUAAEACCKoAAAASQFAFAACQAIIqAACABBBUAQAAJICgCgAAIAEEVQAAAAn4/wF5qUoYyz5p/QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7cde4b9390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "XPCAplot = gpflow.models.PCA_reduce(Y, 2)\n",
    "GPLVM_X_mean = m.X_mean.read_value()\n",
    "\n",
    "f, ax = plt.subplots(1,2, figsize=(10,6))\n",
    "colors = cm.rainbow(np.linspace(0, 1, len(np.unique(labels))))\n",
    "\n",
    "for i, c in zip(np.unique(labels), colors):\n",
    "    ax[0].scatter(XPCAplot[labels==i, 0], XPCAplot[labels==i, 1], color=c, label=i)\n",
    "    ax[0].set_title('PCA')\n",
    "    ax[1].scatter(GPLVM_X_mean[labels==i, dim1], GPLVM_X_mean[labels==i, dim2], color=c, label=i)\n",
    "    ax[1].set_title('Bayesian GPLVM')\n",
    "    "
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}