experiment_SVM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SVM experiment\n",
    "\n",
    "Here we supply code to reproduce the results of the SVM experiment in the "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import packages\n",
    "import scipy as sp\n",
    "import numpy as np\n",
    "import matplotlib.pylab as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import ot\n",
    "import time\n",
    "from ot.bregman import sinkhorn\n",
    "from ot.utils import dist\n",
    "from ot.optim import cg\n",
    "import gromovWassersteinAveraging as gwa\n",
    "\n",
    "import os \n",
    "from os.path import dirname, join as pjoin\n",
    "import scipy.io as sio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importing the Shape Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 100, 1300)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_dir = os.getcwd() # Get the current working directory name.\n",
    "mat_fname = pjoin(data_dir, 'planarShapes.mat') \n",
    "# Add the file name to the current working directory.\n",
    "\n",
    "mat_contents = sio.loadmat(mat_fname) # Read the file\n",
    "\n",
    "planarShapes = mat_contents['planarShapes']\n",
    "planarShapes.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classification\n",
    "\n",
    "Let's try some simple classification task. We will classify shapes by class using the Riemannian algorithm. This is an easy classification task if we are given the 'labeled' points -- i.e. points in the matrix are given in order along the curves. To make things more interesting, I'll permute the order of the points in each curve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first import a few classes of shapes (we will want to increase the number for a real experiment)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "num_classes =  20\n",
      "num_shapes =  400\n"
     ]
    }
   ],
   "source": [
    "indices = [0,20,40,60,80,100,120,140,160,180,200,220,240,260,280,300,320,340,360,380]\n",
    "\n",
    "num_classes = len(indices)\n",
    "\n",
    "num_shapes = 20*num_classes\n",
    "\n",
    "labels = []\n",
    "\n",
    "for j in range(num_classes):\n",
    "    labels = labels + 20*[j]\n",
    "\n",
    "shapes = planarShapes[:,:,indices[0]:indices[0]+20]\n",
    "\n",
    "for j in range(1,len(indices)):\n",
    "    shapes = np.concatenate((shapes,planarShapes[:,:,indices[j]:indices[j]+20]),axis=2)\n",
    "    \n",
    "print('num_classes = ',num_classes)\n",
    "print('num_shapes = ', num_shapes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABZ4AAAEyCAYAAAB6e8UYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3c+O3MbZKO637QMbsZCchZWt4gvyFQTGwVEAA0GOAe91AdoH8GcYCBABHwwDZ+8LkrTNaJPAyk/6Reqz6B6rh1PF5p8iWSSfBxiMxJmeZrPZVcWXb711OB6PAQAAAAAApXy09A4AAAAAALAtAs8AAAAAABQl8AwAAAAAQFECzwAAAAAAFCXwDAAAAABAUQLPAAAAAAAUJfAMAAAAAEBRAs8AAAAAABQl8AwAAAAAQFECzwAAAAAAFCXwDAAAAABAUQLPAAAAAAAUJfAMAAAAAEBRAs8AAAAAABQl8AwAAAAAQFECzwAAAAAAFCXwDAAAAABAUQLPAAAAAAAUJfAMAAAAAEBRAs8AAAAAABQl8AwAAAAAQFECzwAAAAAAFCXwDAAAAABAUQLPAAAAAAAUJfAMAAAAAEBRAs8AAAAAABQl8AwAAAAAQFECzwAAAAAAFCXwDAAAAABAUQLPAAAAAAAUJfAMAAAAAEBRAs8AAAAAABQl8AwAAAAAQFECzwAAAAAAFCXwDAAAAABAUQLPAAAAAAAUJfAMAAAAAEBRAs8AAAAAABQl8AwAAAAAQFECzwAAAAAAFCXwDAAAAABAUQLPAAAAAAAUJfAMAAAAAEBRAs8AADCDwyH+eDjE88Mh3p+//3HOnwMAwJwOx+Nx6X0AAIBNOweB/xYRn11sfh0RXx+P8dPUPy//igAAoJ3AMwAAFHAO/j6NiEcR8TIintwGfQ+HeB4Rf0g87MXxGF/M8PPsvgEAwBQEngEAYKQOGcnvI+KQeOjxeIyPpvx5RHzVtm9XXxwAAAygxjMAAHTUUkf5adwN7Mb5/0/P/36Z+ZMvZ/j5tX1THxoAgOIEngEAoIOLrOY/xCm7+A8R8bfz9keZh91ufxKnLONLr8/bp/55675deV0AADCIwDMAAFwYmNXcmpF8LmnxdUS8iFP5ixdxUepi4p9fy5ZuzYiWDQ0AwBBqPAMAwFlbreaI+DFWWEd5TP3pqPh1AQBQN4FnAAA4OxzieZxKTTS9OH9P/ux4jC/OAd6ncSph8TIintQSnG3btzGvufyeAgCwFUptAACwOy3lI9rqIbfWWT4e46fjMb44HuOj8/cqgs4RV/dtVH1oZTgAAEgReAYAYFeuLKaXrYd8rc7yWg2tD21RQgAA2ii1MVJq2uL5R1VOswQA2LsrpSWehJrGv7pS8/ppKMMBAECGjOeOUtMIM1kef4+IZyHzAwBgUUPKaWw1q3moK8dDGQ4AALJkPHfQkunxOiIedvwzMj8AAGYiU3d6MscBAGgj4/nsSkbG07g7aI7z/7sGnSMiHsn6AACYTW789jSuLBJIZ23Hse34AwCwAzKeoz0j5niMnw6HeB+nshlj3Jz/vqwPAICJtYzfjsdjfJRap8OYrL/ccbx2/OfcRwAAliHwHO3TBI/H+KLl56lg8ps4DbI/udjWVpbDlE4AgIFaAp/PQzmNxXQowyHoDwCwcbvNNrgsexHpQXHEhwVTctMIv437i638KSIeN7Z9HRGfX3kOAAB6yCz0fLuos3Iay8od/58j/54BALAhu8x4zpTWSPk1I2bsdExZNwAAZXWYtaacxoJSxz8s7AgAsBt7DTw/j3yW862i9Zdb6kg/i4gvwwURAEAv6givj/cMAGA/9jq4aytv8Wt5jJIB4PPfapbleBanshymGgIA9Pey53aW5z0DANiJvQaecwPbF8djfHQ8xhdTZB0fj/HT+W9/dJ5K+GXcL/fxWZymIAIAcHa5Psf5uzrO69T6nmXeZwAAVmivgedaLlJymdcWHAQAOMstInj+cXNGWdFZa5SVmQX49fEYP11ZLBIAgJXZZY3niPGLBRbah+dhcRUAgFbGTPvgfQYA2JZdZjzXEHQ+S2VeHyPikamFAAC/MktsH7zPG6WECgDs0+4CzzVN4UtMNTye98nUQgCADyxItw/e5w2q6foLAJjX7kpt1DqFr9b9AgCYU2pm2vlHf4u7izK/DvWcN+UiQHnvfT7/u4YZi/TkOgcA9muPgef3cbrT3nQ8HpfLAK91vwAA5iLwiBsP2+M6BwD2a4+B5+dR4R33WvcLAGAuxkOkOC/WzfsHAPu1xzvMqQX9XseHbIql1LpfAABzsbgcKc6LdXOdAwA7tbvAc2JBvxdRwTS98/M/i4h3503vIuLZ0vsFADAji8uR4rxYsVqvvwCA6e0u8JyqG1fDoOe8X48j4uPzpo8j4rHVngGAHZEZSYrzYgUOh/jj4RDPD4d4f/7+63XM8Rg/HY/xxfEYH52/L379BQBMb1eB54sFa/4QpwUu/hARf6skuPs07i6YEuf/P11gXwAAJtcMVJ03y4zkjlzGbMSpfnAq0Mm8Kr/O2pW2GwAAMLddLS5Y88IWVnsGAPbkIlB1eeP9dQg004Hzpy7XrrNqnXW6NT4XANRmb4HnaoO7NQfFAQBKM/ZhDOdPXdqusyLiqxAMnYXPBQC12Xwm7eVUo4h4n/m1GhYmUbsO4Irc9NG+24EqPOq5HS45f+rStgCkkoLz8bkAoCqbDjwnao19nPi1RYO7t0GRiPjxvC83oaYhsHOpgHFL/cjvem4XrIY6tAWq4BrnT13akmgEQ0foOS65+rkwzgFgTpsutdEy1ehdnILui9YXU4ML2LtUzcfzj1Jt4+uIeJj4M+8ifWMxt/3F+XlSz/EsIh4ntn99/rf6lFCIcRBjOH/qk6vjrPzDcH3P82u/73MDwNw2FXhODHYeRaU1nSPU4AL2reXiJxdgLuUYpz4id2MyFay+idN+ulCDgiw4xhjOn3UQ7BxuyPVi2+fC9ScAc9tM4DkzoDlGOvBcRcda82KHACVlMpufRvrip68hGc+5G5N93WZPC3zAFYKEzMF5VifvyzClrxddfwIwty11LqlFKw5xCj5fqmnBPrXpgM3L1WaO/rUdbyJdP/KHntufRL6dfddznx5Fpo60GorwQUuNdp8LinGe1et4jJ+Ox/jieIyPzt93F3QeOC7odb3Y4TlcfwIwq9UGnpudarQHMF5EnQv2tS3CAbA6mQue3Gr27zN/Jhdg/jZOtZabbfo3Pbf/FPn2Nxesvsns6/vMa/trCH7ApVw78HSBfWG7nGdUacRNkc7Xix2fw/UnALNaZamNNZbVyDHtDNiKlhqOv4nMtM6I+Hfi92dZyK9lEaQ+Cx7mXlvOTUT80nzOMa8D1sD0bubgPFvemGubLV8Xjamt3PW4dH2OLR9nAOqz1sDz80h3qs3gc7WLVujwgbVLtGMPIr0oYFut5dXUR56oTnW1/RSUZEEr5uA8W9aYRQS3vgDhHDdF3HgBoEZr7YDWWFbjV+rPAWuXacdSQeeIU1+TnNa5ppqPmX3NTVnNleZo+iwinqoHzQ6Y3s0cnGcT6dhPjSl1sqoyKQP67TlqK6vfDEB1VhN4vuzcI18X9OVKAhirGlgB+9ajbnPOy8jXWl6182u499riVI+6GfzIyS5QWHyHYUaNsdvTiHgWG2wHqEemTX4Wpxt8buwN1CNpJpcc1GUx4TGPndXAJKI5boq48QJAdVZRaiMz9appNVOxTIMC1mJA3eam1bTNpZUoRWJqOGu19WnzrIPzsIwetYM7/d6Y56jB0H2do9Sico4A1GYtgefnke7c38Upa3tVneqaBlbAvl1pf1PBUovnZQxcfPFlOJaskLEONXAeltE1aWbuGs9LBVnnTCISSAZg7daSXZubYvXRCspqpDyJiDeNbW/CNChgYc2yGtHS/kZ6Oue3Kyl5NLuWshxttReV32CtVjNtnk1zHpbRqXZwrp/rMhbo+9iF18yZpZaydYEA2IK1ZzyvMlvhPFh4FhGfXGx+GxGPBWmApWSyjY6Rzup5EaebZbJwRup73NfY77E/Wxu7sU7OwzJqLFmy5Hs71/Fw/gKwBWvJeN7aQglP427QOc7/t7ggsKTUgoGHOAVBL72Oc5BZdvN4mSyvnEeZxR6hNlsbu7FOzsMCKl20cZJs9i597JjM7p5k7AOweqsIPM/Yuc/FIAJYXI+yGhHbaX+r1AziR3667qsw7ZYV2ODYjRVyHpZz2U/FKXD/OK70RRPfKC1e7qJPaYshN98HHI9ZSnqM4WY4ANesIvAcMaxzr1j1gwhg2zIXVzkvN9T+rkUuSy/iflb6Z2HGDBVI3MwKbQdLS9zYC4Gy0VIzpO70RUPqE/cMYk6RzX71dQ01sF5z1esCqUENQBerCTxv7G5qaqD0JiIebOT1AfXrVVZjlj3iVy0LEX6eecijDfWRrJAABGvgPC2my+zNXkHcvu/NkGz2DteTU85KHRrUbq73kFr/oYgB19uDXtPGrusBuGItiwtWt6DFWOfXdLso16uI+F3crfu86tcH1O1wiPeRvng5xmn2hQUDK9Sy0FBzMUJ9CLOyCBZr4Dwto8txbBtnnMt19P6bY3S5nuyzD41ruavjpb7Ho+/+jDXkenvga9rcdT0A7dYSeH4eGx4kbv31ActLXCA9iIiHiV/V7lQsc8HWDDrf8l4ymyEBCJib87SM0kHc8+9P+t50DJZ3Cor2DZ6ef/+/I+LjtudPPG7wMRkQGH8ePa9H53oMAOu2lgHW1hfj2/rrAxaUmb7624h42/hVZTUql5lanKP8BnOyfgVr4DwtoGOZi741mAe9Nz3KNly93upRvqNziYmLMVgq6Hxt3PUqs/3qMYn+JWWGXI8OqbPtuhdgZ9YSeN76IDH3Ol4JGgAFpC6QPo2If0aP2ojUIbFQVltfqI7pxiWCLt8tNHaYYqEvKM15OkAquHtt4fcBNZh7vzc9A6ydric7LmjfJ3iaGoNFRLyL9jIWf4xTKcamLosLDqm93Pt6e0id7SHPA8C6raXUxqZrQWVe35s4DaDUfQZGMbV425Tf2K/Me98029ih79RuWILztJ85r8OmLA9R8nX0fN5BY7CW57g5HuP3V/av2trLY57HZxdgnVYReI7Yfkej/ipQivZkfxLv+aNws2HzWgITTZN91rc+PmO7nLvd1FyTd8AChkXe8z7B06HHb2R956HPOctnYsjzbD0RDWDLVhN4jtjXALFtsBERX8VOjgPQjxkURLRnSkXEL6H/2ISWsULTJDccBAJYK+dudzXPmloyKN71unTouTbmtW3x/K75BggA7VaT9TRwkYQ1y9Z9jn0dB6Af9ZyJSNfKfBOnepH6j+3oWhNzqtqZQ+qIQg2cu91VU5O3WWs6In6OhWp2d6wFPbQOcsSIeuQjnrNmvRcl7LHwJAATWk3G897ucrbcqX4dpswDGTVnJjEvJVe2b+kaz9ob1sq5210t2bMt+/EsIr6MDc7k2dNs32v6xgJqOW8BWFfgeXcDxNRgIyJ+jJ0dByBPcJGu9tiP7kGiDfg5ZgrC7C0pgO1w7vZzLQA6R4DUe7ZvfQPJzheAeqzpQrOaaV5zyUzhyh4H04lgXzIliH4bEW8bvzrLtFOqly3hpO9Yl8v+Pk7BnicXY4VvLscO59+f6v0dPBUcFubc7aGtrMSM5RB7l1pgOwaUDxl8vrimBihrTRnPqbuctzu/m6lHV6aZPU5sN50INsoCcvRh4clt6JP1NcdUY1PBWSvnbhlzZZbKYKWPoeeLEh0A5a0m8Bxxb4AYcXfK8G46hEwJjqdhMAa7onQCfSnNsn59LqanDNQI2rEVzuVx5hqLCAjSx9DzxQ0OgPJWFXi+pUO4TwAK9kdbyFj6jvXp855N9f4KALEVzuXx5hyLuElAH0POl6H9pnMTIG+tF5VqfN2ndifsQKO264NQz5lxdrd+wgb0ec+men+fxt1AXZz//3Tk34W5OZfHm61edlutaWgaeL707jdnrHMO1UnVRB+zbezzzPna6U7G87j9SN7ZXOKOp9qdsH0tn/N/RcTnIcOCnpxT61NDjWeZ8myFc3m4xvXOq/Nm/QarNqTfrCU2ASUlYlo/R8SXcbfca0S3GFTXbbdrlw19nuTj9UfLW2vgefGFBmtc5E/tTtg2A1umkAge/C7csKxanxvcU9wM1xaxFc7lYUre1FKigNr0PSeH3MBy3lOzTBvf9Pr8lYo3jXGM+2u59Xme1ONzwezmumn3tvlclrHKwHPEvAsN9lzM711EfJzYPvsAVhYHbIvPNFMThNm2Uhe66uKyFc7lYUr1FY4/W9D38+C8pzY9EhjXqhmMLpGBLUDdw2oDz7emvkhu6Rh+E+kAUM4xIr6KGU/QtmMTH4LnPixQscZA4H1UcmOLbXJzY7tKX+jK1mIrnMv9leor3OxkC/r2r857atIxu3nPugStlfi4YguB50kvkls6hlxmc277TZw+zLPd2ayxHAjQXY9pTj67FHGlz/soDKQWNzRI5kIXKKVgxrObnUxi7pJUPZ/Pec+iOiY2dZGKcY2p8dwM8vZ9ntzj59CpxMder6G20LBlV57tu8pl5vcfZX79o0iv4PxDZnvEzKtmn0/qr+OU4Xw8f/86Tie/Fbyhfk8jHXR+Fxef6b12YEziSdzvwyJOA1IrtS/s4mbUH6L/+5Ebz+S2t+6HVcTZGud1L6m+4nV8mILcVfY6rvcewVmfvnJkv/qr4zF+Oh7ji+MxPjp/bxubO+9ZTOKcHxp0fh0R38b9eNOf4pTkOGTb95HuW7o+T+rxc2XaNgPen0XEX+Ju2/L3wyH+sctxxvF4XPVXxPGPEcdfIo7Hi69fIo7fZbb/seff+Udj2+3X8/Njnkcc39/+/+Jv3dl+/nfq77xf4Jhl9yX1enz58rXMV03thq/9fDX6sP/k+sCl93OPX+f3ZdD7MeaxifOj8/jKl681fDmvBx+z1uuGa7/juPua4qtPfzeyX736GWh5nPPe1yxfifM0F99qfv2j8bjv5ogVDf1ctTw+FRf8/yKObzpsy12Ll/q6jVs2Y4ejjkGNX6svtRHRe/G/ZH3jlt8vUiKjpimuLftyjAkXaQT6qandYJ9MB63LmPejVI1n7RJb5Lwur2ubo8Y2pfXpK4f2q2P7VOc9cxhRv3lTcaBMvDA6bPs57peoLV3io/nYTS5yuInAc0pbJxIR/47uiwUeo8CigFfqLc9a9yWzL7kPiwE3zKjRMb6KiN/F/Y5nMwMB6iYYU5ex70eifYmI+Dz61Yp2M4LNcV6Xp/9gKX3OvaHn6ZDHCTYzhxH1m63nktAxaJ0KUI8JRqesOkC95cDz8yizKGCxwVHipE2doLMElRL78igMuGFRmZtCbyLiX9EzOAQlOCfrUjBrefDfEUxii5zX5Qnms5SWJKvvj8f4psPvXu0P+57fpfpv0vYa1O8YX+rCuThSwfdirMEB6inf/y13+rlFL3KvObdYYN9FMrKOjUUHYsFF/hL7UmyRRmCw1GKCn0bEL8dui5VAUcf7i9TexGnw8jAsNji7xPsxdIHRVFvTdfxRalExqInzujwLqLGIc5/4LOLOomKHiHjcHK+M6Ff7nt9j+l0uJGIT30V6gcjvUjGMVGyj9LaZnif1uv8S3QKdNzF+LMmFZnztfJOreQ31tvmwCXaleUPs07gbdI5IL3w46fXcZjOeI8rUfp7yA3ilHMjo8h4996WtFMgiWdmwN7KDqJ2swDqMzewZ29bsNbOIbXNelyXDkyVNPV7pe3537Xe1Q+16lgxNrV+Vim2kMkLHbJvreYaWctAOL6RjVnTpGtJ9THY9t+nAc0pNg6CWDrHIgoYD9qdXoF6QAcoS1KN2bo4sr8Q4RlsDzEEQjaXMMV7pc3536XdrilPUquU4dpUrr1raXM/TlfrNlRpRQ3qKAPVk13O7u0gsOE21hNy0vogFpuIkpgf8FKeTPSW3HejhctpURDyI+1NwTPWlJqZOL6/EdN3BZQWU32LLnN/jNI9fRETi2gLmMPl4JXPtnNOl31WO47qxMYi5gsFLBp2bmaWvI+J/aYfrlGpHOpTteBERf4pTMPpy2/dxv515E91LfBymGvvsLvAc0buTmHQ/IhEEj9OCTSmPFhgM5zrnVwbmMM5FZsNtfaWH8aEG1NI3xiAldeF0jA/9k75geqNvCA+9CZ9os9T4ZjOc3+M4flSmqrrtHftdCV8Ztze1Ip/N2Qyk5QJr74rtVLu5nicVZP4+6kiypKAZAtS3Jum7d1dqYw1appCk6hTNUX6jOeUnV8tIowY9mO7OGjWmhEXM3C/t3Ux1K5NTh7VZbJnzexzHj9qsrdSLz1BaJh5x6bam8pfRXppgazWeU6+76nOcely0j7nSNUXbHYHnCvUsmj95R5TotB/EKTNz9n2BLVEvlzVzgbSMKWtAXvvb2iy2zPk9juMH46jxnHalrvOLyARbczceOtbUHbVtzufJHBfoZK6+W+C5UomG6lFUMphrOzkj4qvQIEInAnesmSDDcjKrYo/OeLnWJmmz2DLn9ziOH4y3tiztORhvwnRa+u53caoNXqT98UGtVLNeS7QskLDAQijZus+hthtcdVGn7FGka3NZTJA1sNDgQi7HCHFqLx5Hmb43W1/y/PceJH6mzWIrqqoJu0KOH4xUy1pUlTHehOmk+u6I0wKZxWJ5As/rkRvM/RzzB3tz+xJhJV5olVh85xCn4LMFIFgbQYY6PI1yfe+1G8vNMls3oc1iI4YuusmJ48faLJC8xTA/h0QdmMRF351aELNYLE+pjRXJ1AnKFQSfdFpbZl9+DNNgoJWpqGxJoy94dd78eZgeOpuSU1Bb6ku+Dms7ALAR6imvQ8vaV98fj/HNMnsF2zN1SRuB55WrqeZRS0DtJiJ+CbWqoKrPLJTiAm45pW9mubHMnqmvCvsgEWQdvE8wj6k/ay4W1i87LXaBqUOpaddvIuJ3oe4z3FKnjC0qWe6BfoqWPMnUl9RusXmJUljGrLBd2TUNZt0LrvE+wTyexCl2d+lNFCppI/C8fqkLzrcR8duYeeCcqe32r4j4pPGrghHsmbq4bJELg4V0qataoI6ldos9cAMN9sMN1XXIvR/v3RSE4pqzG1OzHYf9YaU21i8xLfBBVFKLUVkBOFELly0zFbJepcqgKEHA1hmzwn4oEbYOmffplvcLClFqg6ua02LjFMxKebTA6r3Zu8lWEmYvEtN3H8ZpAPXVxVR2WLNURuwxTv2O9n1ZRbI4MyU4YEtkQMJOdJktxPIu3qd3iR+bkQLlTDp7VcbzBl1Z5O+zmPHObsvd5GcR8XjOfYGlyAZlDxoZsRF3Mwe17wuRxQndyIAEqFPLWCbidOPATCwYQcYzQ+RqMUbMXLsudzc5Ir6ce19gQerfsnm3GbFxGvw3Lw6078uRxQkdyIAczixGYGJtY5bLNa2+a7ZF2idod/5MPEj8qNh6LjKeNypVizEifoxKsp5kYLEnMp7ZE+17XTJZnG8j4p+hzjz8Sh3zYWSKA1O7Uuv50jHujkHfnrd9erHtdvb1l6G9Z+daPls3EfFtqc+FC8CNytRirKnecm5fXrkjyZbMcQcRKiPDtiKJLM6b8/eHcTdLSH/LbiXWYvC56K5IHXmAnMRYJqeZ+PBJ3A06R5zap7/E3fb+2eEQ/5ApzQ6l+vCIiF9K3oyR8bwjNdVbzuzLmzg1/p/MtR8wpbnuIEJNnPd16zsDQxYoe2Bm0nBmuQBza2mzS5EpzS7M1YcbDOxITfWWM/vyr7gbdJ58P2Bis9xBhJpctO83jR89DBmENehcc14WKDtiLYbhzHJhU2S6rkJqTauSGZVdM6WTdaUL7gdMbZY+XOB5ZzIlOLKD7Skb0ea+xKnWZG4/DABYIxey7NK5b/kl8SM3E5eXG0i+T/SxptCzF4Knw+UWNVdOjNVxw3UdMkls38f9tuhtnGZV33n4iKdOLZ7dKRidimeIcVCBJ3H/M/ImCvfhSm3QNlWlWZx/6vIbuf24iVOjbtESVsXUXfbM9Ou6NEpmRKTfm1uv43RBl1uUOOI0RjDNlE2wQN44SvKwFcbu65Zqi84/utz2c9wvM1pal0UOlfNgcefPzLO4W3ngbUQ8VuOZojKD7WZjeWuyTrdl0P86TlO0Z9sXGONKgMeFLLvg4q0eLf18RMT7iPg48bAX5+/Xaihq09gEwVPATfN9aLT3ryLit3E9KJyLj5Q2a/If+zbX9ZrGk9xUlZzJSgS01KDOluCYal9gqMQUvUOczudfz2kDB3bC9Ot6pEpmHOIUXMuNBR9F+j1sUnqDTciUowP2RdmdHWi097+PiD/F3RjE48S2VCmPKbI4U+U8/ls5DiYyS2lQgWciIllvOdvpTlmLKDPoz+3LKzWRqFA2wONClj1J3Ex8FxG/iYin2uvZtQ0qs/194j3s+/dhNdTaBMJN811KxSAS276JbnWlpwhGfxwWNGQas9xsU2qDpJayF8/ifk2kqWs/p/blTZwa38taNKahsDhT9OAutVOXd2UNhW+j4/ujfApbpZ0Cbim7Qx+J8yVVQ7prjeeu5TyU46CI8/n797h7Hr6JiD+VPJ8EQUhqKXvxZcy8wn1mX/4Vd4POk+8HdGSKHtyVmgWgvZ7Xkzhd4DT99vz9Xn+fGWzKBGOrtFNARCi7Qz8dM6NTpTu6lvNIUY6DkprnU/Fa5jKe6aUtmzNOgbVZ7gzLKqVWsqbgLu31vHKZWodD/CMKLNQrE4wt0k4BUIPGOCu3APQ1rj3pxOKC1Kota/N2MbXb2kNT3mlbpAY1tLkYKPwmTvVsLSgIZgHMJrG46WV/XGShXplgbJR2aiTjb4DxLsdZEfG/YlgNaTN26MriglQpNc02VYto6sYuN93358hfdMNkEgGfjyPi3yEbEFLtdUTEA21zcW3lAkYH1gSW2DBlZEa4ctMLgAEyJUe7luN4ZMxGBxYXpE6JabaPYoHpianpvuf/W/iI2Vl0C/LO7fVf436pB1MBC7pSDuurGFEGSBkhtk4ZmeGMgQDm07EchwUIuepwiO8i4i8x8bki8MxobYPN+BAMVvuZTXPuQTuBiem1HON3cZpg0znqAAAbJ0lEQVSuGTGwT/b+ATnGQKyRm01sQSYxIDUjPeI0HvwonO9E67nz/XmRzGIMBCihprIXuSkBr0w1YWLqQ0K7WWqI7UFLyYtcWZOP49Qfx4j6zN4/Nk0pmVGMgVgV5WHYikw5jpyPw/nOB6kSfYeI+LL0E8l4pohayl5k7tq8idMH6JOLbaaaUJRp6NBOxmwZ19qa88//O9LTLgcfa+8fW6YPH8fxY230aWxZy/nd5HzfsTlnKwk8M5mlpt0lguAP4n5d0QgNLYUk6myZwgQNAhNltFxM3ByP8fvz7xTvf71/bJkg1HjKFrAmysOwZZkxW84xtNm7M1WiSo5GlSktUvbieIyfLqcSR8TnmV99ZFolYyWm6n0cEf8OnTfckZgK+C4ifhMRT7W97S77qshnsDy8OI7Fp71npnIKOrMVSsmM1Bx/axuonPIwbFZmzJ2j9MbOXMQvUkHn13GqXlD2OWU8M5XMnba3cWr8Pr3YNmnGVFt22HnfZG8xmCwp6EfmbD89s1ZeHI/xhWMM/ejLYV/0k+xJj7GkxQd34Npi5FO87wLPTKqGshctA4vXc+8L22OqHvQjwNNPjzp9ERftjmnv0J0gFOyPfpI9aZzvqWvXJn3gRi0RvxAUYVI1lL3ITQ9u25cSz8tumKoH/ZjS3k+f4/JruzPFtHflqdgqpWSmoc2gZsrDsCeX53uc+rhrPouIv2rDN2n2+IXAM3PL1n2Ou3Vyi9YZygwsFqlBzXacz40HiR9NUhsJNsLNmn5yx6U5ZW3SdidRz149QDZFEKosbQZAtZ7Eadx4zcO424Y/OxziH+Ikq/dzzHwdIfDM3FKN3O3/mzWHPovTdJA59+VNRPwuDJK54uKCqlmu5SZkSUGbXD/gZk1a7nh9H/NmZz6N+ftpmIXM3EloMwAq1HPxwUufxOna9zJO8l2z/9Sn1uv8XjyOu6U2jhHxbMrrCDWemV2qnlZE/BgL1MmtoQY166ROLQzXaHvfh4VMWtVQh1I9e7ZKfef+urRJ2gyAdei5kHXTMe629W/P2z692PY6Ip5FxJehpvqiloph6PSZXc+yFy+nvGNWQw1qVkudWhjo3O4/iYh/R8THYYZJq0pKACiRwlbJzO2hRwmNTm2G8TVTcF5Bd5l1Dm46Prx5g/GTuBt0jjj1qX+J+/3GvWzpoa+BzhaJYQg8U4vcVOKfY976cIvUoGaVBGFgHMGedVEiha1yI7mfrm331TZDHWimsNbzSrCcJSUS8r6NbnWgu2oGqDsHo1OfDZ+XwRaJYSi1QTUyJTiexoxTAVqmW74OJTi4YGoujGMa9vrUUPIDSlM6q58+bfe1NsOxZwo1nlcdPguuK6hO47x9FRG/jbvZzM0yGyV0Kd2hnMdAS7U1As9UrW1wG6fGpHijUlMNaup1Pk/+Gh9uSNxExLc6N+imxgtDYH8EfPop2Xa7AckUajuvurQxxkSsQSJO8nOcFqq7PLdTQeEpAtQpqaD1P+NUTlUg+uxwiO8i4s9xKnf4LiJ+OB7jmymfU4dO7dpS/ieZPtWzBvUrUzz252IAeZkFP2QxBtgzpRtWyNRGtiZT21LQOa9k261sGVO4el4N7csGPq5LeRolf6heIk7yTdzvPx9HxJ8a276P+/3GFBmwqXrTD0Nd6V+dX+/jOAWd4/z98dTHQcYzVcvcIc7dMZvsjnBmP96c9+OTi20yZHZAVgKUYebAusgMBSLKld3RplBKh5IAv55XQ8+7EY+7moHt2oKt65gtPUdmdPM5XscpeB6x4XJyF8c/1c5ETNzWyHimapkslJzJ7ghn9uNfcTfoHGFhrL2QlQDlXA44H8YKFgDaMQtCsjmy+PvLzA4c9HdCtjkd5T6ricUEbzMcbyJ9Xg3ty4Y+rktmv1lgbFrHbOlUZvTbOCX8XdvWNaM2tcjhX+P+gqTPDof45+EQx/PXP7oufBiRbq9Kb8u+wPu/+93F68uZNI4h45nVabsjHB8WJJz8TlVt9cOYj6wEKMNnaV30e2yNjNsyLDzK1No+q9FzMfqhfdmIx3VqZ3yOILveVnTYlsqgnlrbIoddal+P2ZbL1B6aST7ptZfAM6vT0nmnPuCTXTzUEgBnfi5UoQyBzHVxo4CtcU6PZ0zEHK5cdz2KHmOJoZ/7Me1FLqgs2AzldCi5M9cih+/iQw3lKd3Eqe8dW65k8j7bRR2rk5uWFxFfxrxTgHNTon6O+1M1TB3fiIsO7Tdx6lRMDYXhLCy1LqYCszVKZ403eQke5VCI9s9q37HE0L5scB+YKk+TKBHimhFGaHzOfh/dFjl8HacAbklzBJ0jTmWFmv1v36DzLHEMGc9sRlvmXJwGHsXvJGemgvSa7sV6yOqBsnym1kd2Flsi43m8qWeu6CeI6DTTtNc5MrQvK9kHan9gfi2lPJptyBhzZTz3kVxUca5+VMYzW9KWITfJneTM4iqLZs/ICpmUhbWgoIsZLJeZBs1MBCpy2+9FxFfnTT/qa1gxWfwdXBlbTj1zxdiLiJbP6pBFKocukNn2uAHXYGZcwMxSn+FEG3ITEf9/hz+XWuTwdUT8EMMXSey6rS1Tu5ld/DpO2d6LLeQr8MyWpAYkqRo3Uw9WcwPtV1MHhFumbN1bgbX0c++EASJM4zKo8DBMNa2a6cFsxZCA1d50+LxPHbw39uLqZ3VoILmUgf2icmNQiUSZjv8ddwO7NxHxX3G3DXoc98t5fH08xjdxv71K/e6YbV9HxLeR7n9TQeZvFm0jldpgSxJTJ3otNlFwH5pTNd6c9+OTi22Dpze0LFDxPNJTthadWrEVpsRBeT5X6+M9YwuUjemmy+d9ymOpvWENhpynysgAY61lLCPwzKZ1qAc21SC52QA8iFMW37396DtobhukRMSP0b2gvAF7TwaIUN7U9UEpz3vG2unPu1v6897nvVrLBTj163suDf2cOGeBPXBxwNblpv/9HBNOE25O94qIzzO/+ihXD6ylTlhbrbs+U7NMURzm8ny6CRepMJappuvjPWPt1A3ubtHPe9dyKEoAUcqcZTOWLhECMAeBZzYtN1iNiC9j3guObN3nyNRkzmz/Y7TXusvVuU55r+ZzdxeD0MvM9VIr38KeWdxrfbxnrJ26wd0t/nnvGJxzM4F7Bi66PuRcWvxzAlArgWc2LzNYzV5wDBygXJMbjESkBzZ/zmxvy2p+mQm0f5947oiIj0NGSB8uaGACF+3W5QIeqTaLSliQjQ2Qtd/Rij7vbiZwx4gs+N7n0oo+JwCzE3hmr/pmII8KyrZkXudKcHyc2Z7Lav71jnoi0N5cVfVd4u8KoF7nggamdXlj52G4IVa1274mIr46b/rRDBpWRHZiD2PLAQxN6uj5uEE3E8YknAx57NzPt3NDk0aUzQAoSOCZveqbgTw6KJsZjOQGMKngcEQ+q7n1jvrlc0f+cy+A2k52FIsaUA9+TcwoWCE1VVkr2YnzGdpODHhc75sJY9qwIY+d+/luH7eBMcJQQ5NG3JgCKOhwPObKv8K2pVYRjogfY8aVu1tW6n4WEY8T20dfFB0O8TxOg9WmF+fsNRL6rKoOY2TapoiebcX536tYKX3oavAsS3/C2qTa11rbxa0Y2k4MeVzf93dMGzZw/+Z+vs2NXfu8xyOPt7YCoBCBZ7jQMkC5iYhfYoLBR25gM9WAJzMIfRMR/4pT6Q+Dq4TzcftrfFhc8CYivnWcKKnlIvF13F3Y8ta7SJfmuTn/jVVcbApgrpMbBqzJFoNwS+o6Th3aTszRvox5jiGPXeD5nseG+ta+n2GfeYA6uCiAu1JTq95GxG9joqnEuXpgU9UJS0wvvYnT63oYpkonXQxcLwN/zbIA0FtzCmycbm6kSk6kgs4R+XrwDzN/p9bSFaa1rpMSRKyJkj6F9Cz7MLSdmKN9GfMcQx7b6TGZ8hhDnm9r65P0+gwrqQNQB4FnuJAZoPwzIj5t/OqqL1QaNZ9/iYhPGr+y6tc3ARerFJe5cM8FmHNy9eBzHtVY69HF4Wq5YcCabC0It6Q+46Kh7cSgx/WsaTymDWt9bGY/rj5fLqgfET8P2NfewerKa0L3/gxb8A9geQLP0NAcoMSp/ETKVi5UXIhd5xgxSuZCLnXhnnMT6QvOHzLbb1r+VpULwZ0vBp/E6YL4UUQ8rWXfSHPDgJWRoV9O53HR0HZiyOP6LsA3pg1re2xL8Dg6PF8uqP/lgH3tFVife8HYAUFun2GAFVLjGa5You7znLZW/20KWzxGFk2ZT0uNwd9Eul5jU+tigT0WIzxmnu/F+TGLng9qMa6fdoWaaWOGS3y2H0R6hs6i46JaxmsjF7UrWtu61sX4hnwefYYB1kngGa7IDHLeximIc1mCY5UDH4O469Z8jHoEJVuDm/Ps7fpljvfTSF/ItS0MOPqmVmJfHkXmYjYi/h0Ln9+1BAwYZs3tJPvh5kh/axoHT7kgYc8A7phFBJ/HQn3hiIUghwSRn8eA1+kzDLA+As/QwVoyPYZqvL5X582fhwHdrw6H+C4i/hynQOG7iPjheIxvlt2rdi0XAq8jff7enH9X4GiAAZnNswZ8Wy7ycgHwWTOhpwwYMD03DqiVQNU4Y2b+zX3sp2qH+gZWC2QOL3ITb0QwuPfj9PkA+6FRhw761H2ufFGOpNvXFxFfxWmg+zAqrAG7lPPrfxwfgnMfR8Tjmo5LjxrCt+9vysPM7/91bef0HJrHPCL+Gunj9z7zJ17GvPVxc7Uec2OBR5Gu9fjdROeD2o3rphY+1Zm7Zu1G5T7Dn7ct2rbQsZ9qwdO+i0wP3o9rtaMnHo8N3e8h7b8+H2AnZDzDAFeyP1abMSpjLa3241KghnBfuy/LkTnmOVWUsogoVgqkWSu6yGtRqmHdam8n2Sfn5XhzZsFm/s6Q2sFFxyZDsnNL78dcfeSQ/R6Y8azPB9gJgWcYYEAJg1Vc4Jj2llbTcSlYQzh1gyR3/qas+iZLHz2PeUoVi/flFLpxUeQ1mhK/XoII1Kim/nuthn62Sxz7WtqVGm5g1LAPOSPOEX0+wA4YcMEAuWlw0VKCY6ZdG8u0t7TZj0tqOmVu2mrkz6+PIj1l8ttIn7/fJn4/J1eW4+kay83kDDjmTa/jfCHVNiV5SS3tWZ/zO1eWY7XvPf20TQ9fdMfYO+OakUZ8tksc+74lLqYyVQmPPqotZzT0HKl5bARAOTKeoaCWbISIc0ZgzYOqWjJLajP3cRmQUV9scbgeC2nmZMtKnP9dbWZLwWzy1sWW1iJzHjbLbNzKnoNdM7G0P9sjm42laVeWU+LY15Sx3qc9m6jcx/OoNOMZANoIPENBHeq+Vn+xI1CQNuVxmTLYO2Hd3b5B8KpLcwwoN1FN3eYpJc7Nn+O00Gaf4/RVdPjsuKjeFgE/amFcs5yxx36N/cJUbZ82FYC1EniGwi4G2dnM58oHyy7OJpLJqI3ovkhdzqQ1hHvs95AFDRepf9wj2F8sm3wremaGd77hUFNmG+OtMWDEdhjPbMMag61Ttn3OawDWSOAZJrK2IMoaB/dzKpC10zdzOKWqrOECi+5FLJA53GFmwqVdZDaPNWTB1bgbuH8Vpxr5qTZToHKF1tYHsh3GM+vQdVy1tmDr2LZvba8XAK4x8IfprG1Bm1oWcKlOboG53MJpmcX1cse3a9A5uyjgUhckmUVhcgvw3GT+zPuYcJHCHu9Fzsuo6JjXauCCq5efqYeRL9Ux5+JNlLO2PpDtMJ6pXJ9x1QoXoBvc9vUdb7b9na0s8gzA+sl4homsLeNGdlpen2mTA2oF56x2kbpCpTlymcbPIuLLaByXicqBVPt5XYuWz06uhEnT7tuftVpbH8h2GM/Ub8uleMa0fSWOi7YXgNoIPENhjQDYq/Pmz6Py4OGWLwLGaruIjdP72qXUxCoX3SupZ2mO3PE6xt334jYYnVrwbsgCiKsM9teqwI2Y3bc/a5ZZnPLejaPl9pAtMp6p39ZvDgwtl1HiuDj/AajN6jt2qEliitzDOAVcvlrB9MBciQTT3NunR96ZDhmni4yUjyJ9fKsqnzGlnqU5cv1T84Lss4j4c/QrY5J9L1Y2nbd6LSU4upRb0P6s3OVnPk7v5eMYOYUcOjCeWUDP8g6bLsUzojxIr+OSOea5cWhuOwBMSuAZylplXcGLzIzfxCkbdNMB0AGeRMR/EttTQdD3mb+RrRW8wvqFxYwMTN7qUrLhkrrNM+pxw+FtnLLOvSfbtMr+kfo1g2/nzdr4GQ2oTezmQFrn45I75vFhtmXTJoL6AKyPUhtQ0BqnDqoFd935GP0Y3UoD5GoTO549ZM7LZpmNW7svY7JGQ6cis05r7B+pnzFMHYaUd9AHpHU9Li3H3NgHgKoIPENBa6yrtsZ9nlvLMUp5ER9qF7uYGiFTHzZVyzlX4/nr87+9F7AwfQ1TcF7VwY2l+V1Zf+SrMPYBoBICz1DQGjNvXCxc13KMmqp+r7cglwkkcwrqtsb+kfoZw9TBDYD5yTIHYC0MyKCA2/qCcSrH8Doqr1N6WQ8x2msSc9J2LKp+r7cmVw97z3WyYQ1y9dxvbxz1WJQMLm16kboVUbN5fr2O+YA63ABQhMAzjJQYyD2MU0bXV1MGwFIX6rmL98b2f0TE3y/2N1Ub18XCXanB/TEi/ut4jN8LdgJcl7pBJBjCSAKeFWi7sbTojm3YgGNugVcAFqHUBoy0xPTCzJTlt3EaeH56sS1X/zblXZxuRpl6l2B6IkB5pujTR6ovPv9I/1wp46c6KEsDwFIEnmGkJQZyPRe7exfprOYmA08AZiUYQlfqhK+P96webvIBsBQDehhvifqCj3r8bpegc4R6iADMT41eulIqYH28Z/VQlgaARQg8wwjnTI4HiR9NPZDrc0H+rsPvGHgCsATBELrK3XTvczOeeXnPKqEONwBLEXiGgS6mDz5s/Ogmph/IpS7U30bEm8a21xHxQ+Z3b8LAE4AF5YIhEaep4c3Fctk12fHr4z2rSGqB16X3CYDtE3iG4VLTByMifpl6IHf++8/iQzbzuzgFwf8U9zMZvon7F/WPj8f4vYEnAEtrBkPOm/8Wp3qkh/P3vwk+757s+EodDvHd4RD/ORzieP7+3flH3jMA2DmLC8JASy6IZLEWALbKIlicxzlP41SS4WVEPDke46fc9uX2lHOQ+f8kfvRfx2N84z0DgH0TeIaBlrwwdlEOwFYteWOX5bm5vi6HQ/wn0gtZvzse43/MvT8AQF0M3mG41PTBiIgHM0wHtlgLAFulLuy+pUqZfXbeTn1SQee27QDAjgg8w0AXCyLdNH70MKavRemiHICtUhd239xcX5d3PbcDADsi8AwjnIPPvyR+NHVmjotyADbp4sZuc7Hcnw6H+OPhEM8Ph3h//m7Bwe1xc31dfui5HQDYETWeYaSlalFarAWAPVH7d3tSY5nzj7zPK3JeYPDPcSqv8S4ifjge45tl9woAqIHAM4xkoT8AmJ7+dlvabiSc/+3mOgDAygk8w0hLZWDJeAZgT5aaYcQ03EgAANg+g3QYaYlalBfB7j/E6SL8DzH9goYAsKRs7V+1n1fJIoIAABsn4xkmMmUmtCwhAPampV99FhGPE9vVBK5AboaWsQwAwPbJeIbpPI27F8Fx/v/TPn8kk8UlSwiAXcnNMIqIL6NAf0t5V2ZoPYnTDYJLr+PDAoMAAKycjGeYSN9alD1Xdn8dEQ8Tf1uWEAC7ovZzva5lNVuvAgBg2wSeYSItF1s35+8PL/7/fyM9TTgXYL45/65pxQDsmuBmvdwUAADYNwM+mE5qCunbiPifcTeY/DAi/k+kpwmngs4REZ9HZkHDkfsMAGuTLdlgMd75ZEqDZReEnHPfAABYhoxnmFAiy+pB5IPJfSipAQBnFrBbloUfAQBIEXiGGbVMOc1RUgMABrpW6kEZjjLaAvxxykh3jAEAdkipDZhX29TS5l2g1xHxbSipAQBDZUs9KMPRX6acRsQpqJzy6HiMn47H+OJ4jI/O341hAAB2QsYzzOh8gfb3iPi08aO3cbr4/TJkBAFAES0lIL6OUxauMhwdOZYAAPQl4xlmdA4k/ylOJTRu3UTE4+MxvpERBADlnPvS3MyhbJZuRGt272Zdec1PI70Q8tNoWeBxur0FAKB2Mp4BANidDnWJk9m9W70x3JbRfF6oUb1sAAB6EXgGAGB3xpSO2GKQtS0Qf37NrT+fcNcAAFgppTYAANidoWU41r4o4ZAFAs/fldMAAKAXgWcAAHbpeIyfMusrvMw85GW01zr+1ZI1onPPfSVo3vaarwXqAQDgHoFnAAC4qy2791pmcGuAtxEU/k8qONwWsO7y89xzx8gFAlsC9QAAcI8azwAA0JCr49yl1nHL79zEKdDbDP5GnIK8zyLiceQX+GtdAPDKc784vxYLBAIAMAuBZwAA6Khj8Pd9pAO817yLiI8T2zsv8Nfy3Mc4BZMtEAgAwCyU2gAAgI461jrO1Uu+JhV0jvhQxuNqmY+W534ZFggEAGBGAs8AANBDh1rHuQDvzZU//S6z/WXje+7nbc/9xAKBAADMSeAZAAAKygV4I+LbuB8UvvU6In5I/PwyI7nTAoCp574NLlsgEACAuajxDAAAM2ks4Pc+Tokgl4sXti7wZwFAAADWQuAZAAAAAICilNoAAAAAAKAogWcAAAAAAIoSeAYAAAAAoCiBZwAAAAAAihJ4BgAAAACgKIFnAAAAAACKEngGAAAAAKAogWcAAAAAAIoSeAYAAAAAoCiBZwAAAAAAihJ4BgAAAACgKIFnAAAAAACKEngGAAAAAKAogWcAAAAAAIoSeAYAAAAAoCiBZwAAAAAAihJ4BgAAAACgKIFnAAAAAACKEngGAAAAAKAogWcAAAAAAIoSeAYAAAAAoCiBZwAAAAAAihJ4BgAAAACgKIFnAAAAAACKEngGAAAAAKAogWcAAAAAAIoSeAYAAAAAoCiBZwAAAAAAihJ4BgAAAACgKIFnAAAAAACKEngGAAAAAKAogWcAAAAAAIoSeAYAAAAAoCiBZwAAAAAAihJ4BgAAAACgKIFnAAAAAACKEngGAAAAAKAogWcAAAAAAIoSeAYAAAAAoCiBZwAAAAAAihJ4BgAAAACgKIFnAAAAAACKEngGAAAAAKAogWcAAAAAAIoSeAYAAAAAoCiBZwAAAAAAihJ4BgAAAACgKIFnAAAAAACKEngGAAAAAKAogWcAAAAAAIoSeAYAAAAAoCiBZwAAAAAAihJ4BgAAAACgKIFnAAAAAACKEngGAAAAAKAogWcAAAAAAIr6f0lJg0LBHQCdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1800x360 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize = (25,5))\n",
    "\n",
    "for j in range(5):\n",
    "    fig.add_subplot(1,5,j+1)\n",
    "    plt.plot(shapes[0, :,50*j], shapes[1, :,50*j], 'ob');\n",
    "    plt.axis('equal');\n",
    "    plt.axis('off')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert the shapes to cost matrices. We include a permutation step to mix up the points along the curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "CList = []\n",
    "pList = []\n",
    "p = ot.unif(100)\n",
    "\n",
    "for j in range(num_shapes):\n",
    "    shape = shapes[:,:,j]\n",
    "    C = sp.spatial.distance.cdist(np.transpose(shape), np.transpose(shape))\n",
    "    perm = np.random.permutation(np.eye(100))\n",
    "    C = np.matmul(np.matmul(perm,C),perm.T)\n",
    "    CList.append(C)\n",
    "    pList.append(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparing the Training Data\n",
    "\n",
    "There are a few preliminary steps to prepare the training data for the Riemannian algorithms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train/Test Split\n",
    "\n",
    "We split the cost matrices into a training set and a testing set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import some machine learning functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.svm import SVC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reshape `CList` into a big matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "CListMatrix = np.zeros([num_shapes,100*100])\n",
    "\n",
    "for j in range(num_shapes):\n",
    "    C = CList[j].reshape(1,100*100)\n",
    "    CListMatrix[j,:] = C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Split the feature vectors and labels into training and testing sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = CListMatrix\n",
    "y = labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "300 100\n"
     ]
    }
   ],
   "source": [
    "train_size = X_train.shape[0]\n",
    "test_size = X_test.shape[0]\n",
    "\n",
    "print(train_size,test_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we reshape the training and testing vectors into lists of cost matrices once again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "CList_train = [X_train[j,:].reshape(100,100) for j in range(train_size)]\n",
    "CList_test = [X_test[j,:].reshape(100,100) for j in range(test_size)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute the Average\n",
    "\n",
    "To train any classifier in the Riemannian setting, we first find the average cost matrix from the testing set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter Frechet_Loss\n",
      "0 49847302.89841316\n",
      "current size  (100, 100)\n",
      "1 19090568.235538132  full steps\n",
      "current size  (100, 100)\n",
      "2 17506092.506424382  full steps\n",
      "current size  (100, 100)\n",
      "3 17317609.272966314  full steps\n",
      "current size  (100, 100)\n",
      "4 17235061.913078688  full steps\n",
      "current size  (100, 100)\n",
      "5 17173868.357179705  full steps\n",
      "current size  (100, 100)\n"
     ]
    }
   ],
   "source": [
    "CnBase, pnBase, Frechet_Loss = gwa.network_karcher_mean_armijo_sched(CList_train[0],pList[0],CList_train,pList,20,50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Log Matrices to the Average\n",
    "\n",
    "Next we log all matrices (training and testing) back to the average and reshape. The result is our training and testing set for the Riemannian-based classification algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "CAvg = CnBase\n",
    "\n",
    "vectors_train = np.zeros((train_size,100*100))\n",
    "\n",
    "for j in range(train_size):\n",
    "    C0, C1, p, q, vec, opt_coup = gwa.log_map(CAvg,CList_train[j],pList[0],pList[j])\n",
    "    vec = vec.reshape(1,100*100)\n",
    "    vectors_train[j,:] = vec\n",
    "    \n",
    "vectors_test = np.zeros((test_size,100*100))\n",
    "\n",
    "for j in range(test_size):\n",
    "    C0, C1, p, q, vec, opt_coup = gwa.log_map(CAvg,CList_test[j],pList[0],pList[j])\n",
    "    vec = vec.reshape(1,100*100)\n",
    "    vectors_test[j,:] = vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create Naive Baseline\n",
    "\n",
    "We should also test everything without doing any alignment or anything interesting. We could just treat each cost matrix as a vector and do classification in the space of unaligned cost matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "vectorsBaseline_train = np.zeros((train_size,100*100))\n",
    "\n",
    "for j in range(train_size):\n",
    "    vectorsBaseline_train[j,:] = CList_train[j].reshape(1,100*100)\n",
    "    \n",
    "vectorsBaseline_test = np.zeros((test_size,100*100))\n",
    "\n",
    "for j in range(test_size):\n",
    "    vectorsBaseline_test[j,:] = CList_test[j].reshape(1,100*100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create Less Naive Baseline\n",
    "\n",
    "We could also do *some* alignment by picking one of the shapes in the training set, aligning everything to it and pulling back to tangent vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "CAvgLNB = CList_train[0]\n",
    "\n",
    "vectorsLNB_train = np.zeros((train_size,100*100))\n",
    "\n",
    "for j in range(train_size):\n",
    "    C0, C1, p, q, vec, opt_coup = gwa.log_map(CAvgLNB,CList_train[j],pList[0],pList[j])\n",
    "    vec = vec.reshape(1,100*100)\n",
    "    vectorsLNB_train[j,:] = vec\n",
    "    \n",
    "vectorsLNB_test = np.zeros((test_size,100*100))\n",
    "\n",
    "for j in range(test_size):\n",
    "    C0, C1, p, q, vec, opt_coup = gwa.log_map(CAvgLNB,CList_test[j],pList[0],pList[j])\n",
    "    vec = vec.reshape(1,100*100)\n",
    "    vectorsLNB_test[j,:] = vec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Classification\n",
    "\n",
    "Now that we have suitable training and testing vectors for both the Riemannian algorithms and the naive algorithm, we can run whatever classification algorithms we like."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Support Vector Machine\n",
    "\n",
    "Classify the shapes using SVM classification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "modelSVM = SVC(kernel=\"linear\",C=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Riemannian Algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "0.94\n"
     ]
    }
   ],
   "source": [
    "modelSVM.fit(vectors_train, y_train)\n",
    "print(modelSVM.score(vectors_train,y_train))\n",
    "print(modelSVM.score(vectors_test,y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Naive Baseline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "0.19\n"
     ]
    }
   ],
   "source": [
    "modelSVM.fit(vectorsBaseline_train, y_train)\n",
    "print(modelSVM.score(vectorsBaseline_train, y_train))\n",
    "print(modelSVM.score(vectorsBaseline_test,y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Less Naive Baseline\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "0.84\n"
     ]
    }
   ],
   "source": [
    "modelSVM.fit(vectorsLNB_train, y_train)\n",
    "\n",
    "print(modelSVM.score(vectorsLNB_train,y_train))\n",
    "print(modelSVM.score(vectorsLNB_test,y_test))"
   ]
  }
 ],
 "metadata": {
  "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
}