{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learning on Tangent Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lead author: Nicolas Guigui.\n",
    "\n",
    "In this notebook, we demonstrate how any standard machine learning algorithm can be used on data that live on a manifold yet respecting its geometry. In the previous notebooks we saw that linear operations (mean, linear weighting) don't work on manifold. However, to each point on a manifold, is associated a tangent space, which is a vector space, where all our off-the-shelf ML operations are well defined! \n",
    "\n",
    "We will use the [logarithm map](02_from_vector_spaces_to_manifolds.ipynb#From-substraction-to-logarithm-map) to go from points of the manifolds to vectors in the tangent space at a reference point. This will enable to use a simple logistic regression to classify our data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We import the backend that will be used for geomstats computations and set a seed for reproducibility of the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Using numpy backend\n"
     ]
    }
   ],
   "source": [
    "import geomstats.backend as gs\n",
    "\n",
    "gs.random.seed(2020)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We import the visualization tools."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use data from the [MSLP 2014 Schizophrenia Challenge](https://www.kaggle.com/c/mlsp-2014-mri/data). The dataset correponds to the Functional Connectivity Networks (FCN) extracted from resting-state fMRIs of 86 patients at 28 Regions Of Interest (ROIs). Roughly, an FCN corresponds to a correlation matrix and can be seen as a point on the manifold of Symmetric Positive-Definite (SPD) matrices. Patients are separated in two classes: schizophrenic and control. The goal will be to classify them.\n",
    "\n",
    "First we load the data (reshaped as matrices):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geomstats.datasets.utils as data_utils\n",
    "\n",
    "data, patient_ids, labels = data_utils.load_connectomes()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We plot the first two connectomes from the MSLP dataset with their corresponding labels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApsAAAFbCAYAAACakkVNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHB0lEQVR4nO3deXjU9bk+/nv27JOEhCwQYhIMiAi0CJRaFo8s0lMr4FpbBaviArZKrRZrRdQ21fZgK18U/Z1WtNalVHFrjxtrewqKVASqB1kCiUASCGTPTGZ5//6wpI2JkzuSD5lM7td15bpk8vjZ55knk8nnthljDERERERELGDv6Q0QERERkdilYVNERERELKNhU0REREQso2FTRERERCyjYVNERERELKNhU0REREQso2FTRERERCyjYVNERERELKNhU0REREQso2FTuuy0007DN77xjU7r1q9fD5vNhvXr11u/URGsXLkSNpsN7733Xo9ux2fdc889sNlsPb0ZIiJtdKU3nag9evSoxVslvZmGzT5kx44duPjii5Gfn4+4uDgMGDAAU6dOxbJly3p600RE5J86+wF58uTJGD58+Cndpp/97Gd46aWXTuk6JXZo2Owj/va3v+Hss8/GBx98gOuuuw7/7//9P1x77bWw2+349a9/bck6J06ciObmZkycONGS5fd2d911F5qbm3t6M0REOqVhU06Gs6c3QE6Nn/70p/B6vdiyZQtSU1PbfK+qqsqSddrtdsTFxVmy7FPFGAOfz4f4+PhuX7bT6YTTqaegiIjENr2z2Ufs3bsXZ555ZrtBEwD69+/f5t9PP/00xo4di4SEBKSlpWHixIl488032/1/f/3rXzF27FjExcWhsLAQTz31VJvvf/Yzmyd+NdTR1+TJk1v/v2AwiPvuuw9FRUXweDw47bTTcOedd8Lv97dZ/onPjr755psYNWoU4uLiMGzYMLz44osdHgO/34+FCxciMzMTiYmJmDVrFo4cOdLhMt944w2cffbZiI+Px2OPPQYAqKmpwS233IK8vDx4PB4MHjwYDzzwAMLhcOv/v3//fthsNvzyl7/E448/3roPY8aMwZYtW9qs6/M+F8UefxGRE55++mmMHj0a8fHxSE9Px+WXX47y8vI2NX/5y19wySWXYNCgQfB4PMjLy8Ott97a6W9YbDYbGhsb8eSTT7b27Llz57apqampwdy5c5Gamgqv14urr74aTU1Nrd+fNGkSRo4c2eHyhwwZgunTp3+xHZdeQcNmH5Gfn4+tW7di586dEeuWLFmCK6+8Ei6XC/feey+WLFmCvLw8rF27tk3dnj17cPHFF2Pq1Kn4r//6L6SlpWHu3Ln4xz/+8bnLnjhxIn73u9+1+br//vsBtB14r732Wtx999348pe/jIceegiTJk1CSUkJLr/88nbL3L17Ny677DLMmDEDJSUlcDqduOSSS/DWW2+1q7355pvxwQcfYPHixbjxxhvx6quvYsGCBe3qdu3ahW9961uYOnUqfv3rX2PUqFFoamrCpEmT8PTTT+Oqq67Cww8/jHPOOQeLFi3CwoUL2y3jmWeewS9+8Qtcf/31uP/++7F//37Mnj0bgUDg8w8++OMvIrGvtrYWR48ebff12T7y05/+FFdddRVOP/10LF26FLfccgvWrFmDiRMnoqamprVu1apVaGpqwo033ohly5Zh+vTpWLZsGa666qqI2/G73/0OHo8HEyZMaO3d119/fZuaSy+9FPX19SgpKcGll16KlStXYsmSJa3fv/LKK7F9+/Z2r0FbtmzBxx9/jO985ztf8ChJr2CkT3jzzTeNw+EwDofDjB8/3tx+++3mjTfeMC0tLa01u3fvNna73cyaNcuEQqE2/384HG797/z8fAPAbNy4sfWxqqoq4/F4zA9+8IPWx9atW2cAmHXr1nW4Tc3NzWb06NEmNzfXHD582BhjzLZt2wwAc+2117apve222wwAs3bt2nbb8cILL7Q+Vltba3JycsyXvvSl1seeeOIJA8BMmTKlzX7ceuutxuFwmJqamnbLfP3119us/7777jOJiYnm448/bvP4j370I+NwOExZWZkxxpjS0lIDwPTr188cO3aste7ll182AMyrr77a+tjixYvNvz8F2eMvIrHtRM+K9HXmmWcaY4zZv3+/cTgc5qc//WmbZezYscM4nc42jzc1NbVbV0lJibHZbObAgQOtj322NxljTGJiopkzZ067//9E7Xe/+902j8+aNcv069ev9d81NTUmLi7O3HHHHW3qvve975nExETT0NDQyVGR3kzvbPYRU6dOxaZNm/DNb34TH3zwAR588EFMnz4dAwYMwCuvvAIAeOmllxAOh3H33XfDbm97aXz2173Dhg3DhAkTWv+dmZmJIUOGYN++ffQ23XTTTdixYwdeeOEFZGdnAwD+/Oc/A0C7dwt/8IMfAAD+9Kc/tXk8NzcXs2bNav13SkoKrrrqKrz//vuoqKhoUztv3rw2+zFhwgSEQiEcOHCgTV1BQUG7X+msWrUKEyZMQFpaWpt3GKZMmYJQKISNGze2qb/sssuQlpbWZl0AIh6frhx/EYl9y5cvx1tvvdXua8SIEa01L774IsLhMC699NI2vSk7Oxunn3461q1b11r77589b2xsxNGjR/HVr34Vxhi8//77J7WtN9xwQ5t/T5gwAdXV1airqwMAeL1eXHjhhXj22WdhjAEAhEIhPP/885g5cyYSExNPav0S3fTXCX3ImDFj8OKLL6KlpQUffPABVq9ejYceeggXX3wxtm3bhr1798Jut2PYsGGdLmvQoEHtHktLS8Px48epbXnsscfwxBNP4LHHHsNXvvKV1scPHDgAu92OwYMHt6nPzs5Gampqu8Fw8ODB7Qax4uJiAJ9+fvLEENvRNp8YBj+7zQUFBe22d/fu3di+fTsyMzM73J/P/pEVu65/15XjLyKxb+zYsTj77LPbPX7ih17g095kjMHpp5/e4TJcLlfrf5eVleHuu+/GK6+80q4X1dbWntS2Rup5KSkpAICrrroKzz//PP7yl79g4sSJePvtt1FZWYkrr7zypNYt0U/DZh/kdrsxZswYjBkzBsXFxbj66quxatWqLi3D4XB0+PiJn1gjeffdd/H9738f1157LebNm9dhjRXv5LHb3NFfnofDYUydOhW33357h8s4MeB2dV0iIicjHA7DZrPhf/7nfzrsO0lJSQA+fRdx6tSpOHbsGO644w4MHToUiYmJOHjwIObOndvmDx2/CKbnTZ8+HVlZWXj66acxceJEPP3008jOzsaUKVNOat0S/TRs9nEnfmo+fPgwBg8ejHA4jA8//BCjRo2yZH1HjhzBxRdfjFGjRmH58uXtvp+fn49wOIzdu3fjjDPOaH28srISNTU1yM/Pb1O/Z88eGGPaDKcff/wxgE//sry7FBUVoaGhwdKmWFRUZPnxF5HYUlRUBGMMCgoK2v3Q++927NiBjz/+GE8++WSbPwjq6I8pO9IdbwA4HA5cccUVWLlyJR544AG89NJLuO666z53UJXYoc9s9hHr1q3r8F21E5+RHDJkCGbOnAm73Y5777233U+53fGOXCgUwuWXX46Wlha88MILcLvd7Wq+/vWvAwB+9atftXl86dKlAID//M//bPP4oUOHsHr16tZ/19XV4amnnsKoUaPa/Ar9ZF166aXYtGkT3njjjXbfq6mpQTAYPOl1WH38RST2zJ49Gw6HA0uWLGnXJ4wxqK6uBvCvdx7/vcYYQ4d6JCYmtvnL9i/qyiuvxPHjx3H99dejoaFBf4XeR+idzT7i5ptvRlNTE2bNmoWhQ4eipaUFf/vb3/D888/jtNNOw9VXX43U1FT8+Mc/xn333YcJEyZg9uzZ8Hg82LJlC3Jzc1FSUnJS27BixQqsXbsWN9xwQ5sPrQNAVlYWpk6dipEjR2LOnDl4/PHHUVNTg0mTJuHdd9/Fk08+iZkzZ+Lcc89t8/8VFxfjmmuuwZYtW5CVlYXf/va3qKysxBNPPHFS2/pZP/zhD/HKK6/gG9/4BubOnYvRo0ejsbERO3bswB//+Efs378fGRkZJ7WOwYMHW3r8RST2FBUV4f7778eiRYuwf/9+zJw5E8nJySgtLcXq1asxb9483HbbbRg6dCiKiopw22234eDBg0hJScELL7xAf85+9OjRePvtt7F06VLk5uaioKAA48aN6/L2fulLX8Lw4cOxatUqnHHGGfjyl7/c5WVI76Nhs4/45S9/iVWrVuHPf/4zHn/8cbS0tGDQoEG46aabcNddd7Xe7P3ee+9FQUEBli1bhh//+MdISEjAiBEjuuUD3CduoL5ixQqsWLGizfcmTZqEqVOnAgD++7//G4WFhVi5ciVWr16N7OxsLFq0CIsXL263zNNPPx3Lli3DD3/4Q+zatQsFBQV4/vnnu/0GwQkJCdiwYQN+9rOfYdWqVXjqqaeQkpKC4uJiLFmyBF6vt1vWY+XxF5HY9KMf/QjFxcV46KGHWu9tmZeXh2nTpuGb3/wmgE//UOjVV1/F9773PZSUlCAuLg6zZs3CggULPvdm6/9u6dKlmDdvXmvM7pw5c77QsAl8+odCt99+u/paH2Iz+v2c9FKnnXYahg8fjtdee62nN0VEREi//vWvceutt2L//v0d3tlEYo8+sykiIiKnhDEGv/nNbzBp0iQNmn2Ifo0uIiIilmpsbMQrr7yCdevWYceOHXj55Zd7epPkFNKwKSIiIpY6cuQIrrjiCqSmpuLOO+9s/Syp9A36zKaIiIiIWEaf2RQRERERy0Tdr9HD4TAOHTqE5ORkSyILRUSMMaivr0dubi7s9tj8mVu9VESs1JU+GnXD5qFDh5CXl9fTmyEifUB5eTkGDhzY05thCfVSETkVmD4adcNmcnIyAODA309DSlLkSXniz6/hllnGRQl6qpupuob8RKou7ODeTQgmcHUJVXwkoqsuQNX5MjxUXW0hl12bXBai6mzkJ4UDidy7TsllfqouGMftR1wVdy3UDEum6pLLfVRdIJ57SrZ4uf1w+sKdF/1TYxa3TPbcuRq5wuQD3LHxp7aPN+1I/JHOz10w5Mdfti1t7TexqCu9dNLPuF4ad4y7nlK2HabqTGICVVc3hAtNYK/hjO1NVB0ABFJcVJ09yF3vgURuG1PeKaPqas7Jp+oSDnPPM+Pkeq5r68dUXfXss6g6R4A7fs5msgGRZe467nU1bm8Vt0AA1RO4H2BtYW4jmzO4cxJ3jFteQkULVdeZYNCPd/7yc6qPWjZsLl++HL/4xS9QUVGBkSNHYtmyZRg7dmyn/9+JX/ekJNmRkhz5ADvccdS2OF3cxeR0cI3U6eLWyw6bxs3VsfsBAE4nOYy4uGHT4WGX173DZtjNPcmcTvLXhC5yP8hrgb4GyWeacXGFYXY/Qvyw6XB377BJv3iQxybk4oZNp4P/m8do//XyF+2jgFW9lOyRdq6vGAdXx/Zcuk85+eeFcZHDJjndGPa5ayevd/LY0D2IHDadNm772GvLQTYWJznUs8Om00nOB+Q1DfD7zA6bDg93Thxutud270eHmD5qyYeVnn/+eSxcuBCLFy/G3//+d4wcORLTp09HVRX/k4GISF+mPioiscKSYXPp0qW47rrrcPXVV2PYsGFYsWIFEhIS8Nvf/rZdrd/vR11dXZsvEZG+rit9FFAvFZHo1e3DZktLC7Zu3YopU6b8ayV2O6ZMmYJNmza1qy8pKYHX62390gfaRaSv62ofBdRLRSR6dfuwefToUYRCIWRlZbV5PCsrCxUVFe3qFy1ahNra2tav8vLy7t4kEZFepat9FFAvFZHo1eN/je7xeODx8B+8FRGR9tRLRSRadfs7mxkZGXA4HKisrGzzeGVlJbKzs7t7dSIiMUd9VERiSbcPm263G6NHj8aaNWtaHwuHw1izZg3Gjx/f3asTEYk56qMiEkss+TX6woULMWfOHJx99tkYO3YsfvWrX6GxsRFXX301vYyJP7+m03tVbV38KLWskQ/cRNUN2MrdUsSVEU/V1RZy92drziTvs+nj7s8GAPEHG6k6h5fbRge57sRPuBsHt5A36A6T98907zxA1XmSk6g6GO5+ZSkHuP2wb9xG1cWfOYSqSzhA3icwwN+btSWpP1Xn8JM3Iu7H/Sxrb+ICCNx27lpwHKnttMaEuRCAntQdfRT49IbtnfXS95ZwvbR45Y1U3fFi7o+TbNxteQHydqiJB7nnRelM7j6IAABD9qA6ri6QzD1//CmFVF1jLrfehlzuBvphcirwZpA3a2/h9jfpIHejcc+eys6LADSOyOWWV9lA1e27ZhBVBwDp/+Au7ONDuNdV7z4yTOG5LVSdM5+76Xz9yKyI3w8G+JnEkmHzsssuw5EjR3D33XejoqICo0aNwuuvv97uw+4iItIx9VERiRWW/YHQggULsGDBAqsWLyIS89RHRSQWWHJTdxERERERQMOmiIiIiFhIw6aIiIiIWEbDpoiIiIhYRsOmiIiIiFhGw6aIiIiIWEbDpoiIiIhYxrL7bJ6s5LIgnK7I6SdsMtAHdzxC1X39fy6i6pw+Lh0ge81xqu7g17mbNMcf5dNgQglcMlDCHm4bPUe41KSGfC6lwru9mqrD4DSuLkgeGzJRx6QkUnUuMrnIxHPHz3aQS8cIHefOmyMlhaoDAHd9BldXyx3DkJtLV2oo5FKdkg5wqVjhqqOd1xgurSQWxB0Lw+mKnEDCJgN9PJdLGjrjMa43s+xcyBSasrn3T9ydh0y1ctWTdY1cUo6vhUv8CcaTyTvlXF2IDE0KxnHbF3Jzda5mbvsac7h+EYrjkoH8qVy6jbOJ6/WWIJOxDPm2oD2eO8mhDO51oSUp8opDLfz7lXpnU0REREQso2FTRERERCyjYVNERERELKNhU0REREQso2FTRERERCyjYVNERERELKNhU0REREQso2FTRERERCyjYVNERERELBO1CUKe6mY4HZFTLwZsraKWxSYD/Xn9C1TdWUu5dIysrVw8QO46Lg0mkM6l0ADA8SFcko+rmVtmKnmsw0VcGkMwg0uNcTZyaU0H55xJ1aXu5dJv7IHI194J5jQvVdcwgHuquRq4tI2449xxOTqCS5ICgLzXOk/eAYDjo9KpuvT3j1F1DcWpVN2xM5OpOt+kUZ3WhPw+YNlz1PJ6u5Rth+G0eyLWHC/Oo5bFJgN9dD2X2lb49nepuviPuGSUuKPc86d6FFcHAJ7j3Hsyvn5cv3dyQVioLea2seiPPqru+BCu1we41ozGgdz+xldxx8/ZxO2vu4Fbb3MGt96ESm55viw+wS+wn0svMlwZLdzIXVw1Z3An2Zce+diE/EoQEhEREZEooGFTRERERCyjYVNERERELKNhU0REREQso2FTRERERCyjYVNERERELKNhU0REREQso2FTRERERCyjYVNERERELBO1CUIN+YlwuiKnRrgyuEQEp49LW2GTgXYs5NIxJsy/nqqrHxg53eMENtkBANx1XBqDw8/VBXK5pJyQm0tjODKSSziKP8Yl+aTu49IdEjfvpepMdiZVF0zlkk0SHWRKRSoXKdGSTKZjHOaTUg6fm0HV2YPcMj85n1ueu5ZbXkM+VQZ3LVfXV5jEBBhH5B5j41okjU0G2jflt1TdiG1cbw4kcs8zgH9eGHKR9gC9SHLFXFnI070xNGE+dIxi41o4QnHcgQ4kcL0vwIXZIewkTzB/ycCfyi0z7OAW6vdy++xISaHqgmQYYbCTl+lQFy49vbMpIiIiIpbRsCkiIiIiltGwKSIiIiKW0bApIiIiIpbRsCkiIiIiltGwKSIiIiKW0bApIiIiIpbRsCkiIiIiltGwKSIiIiKWidoEobDDhnAnqSu1hVzUQfaa41Rd1lburv9sMtBflj9G1U2+9jqqzhbmb9efUtpM1fnTufSixhyuLvONfVRd9XkFVJ33Iy4OJpTgpuqQmU6V2WsbqLpj56RRdf12NFF1fi8X7dCQy10LrgY+9iLpMBcjYw9wy6wt4NpLxt/ruPWGuHQMh7/zyJJQgIw1iQF1Q7ydprGhm1Ny4j/ikrXYZKDtt3GpbWcvvpGqSy7l32dhE3DcNdzzgk3o6f8eVxeK5/bFkC8fmdu5NDZjY1PRuO1z+rjjl1TOvba5mrjXLEMmCCWU8+NSxnY/VXdsGLeNKWXcObElczGD7OtRKC7yPgeDfuymlqR3NkVERETEQt0+bN5zzz2w2WxtvoYOHdrdqxERiVnqoyISSyz5NfqZZ56Jt99++18rcUbtb+tFRKKS+qiIxApLupfT6UR2drYVixYR6RPUR0UkVljymc3du3cjNzcXhYWF+Pa3v42ysrLPrfX7/airq2vzJSLS13WljwLqpSISvbp92Bw3bhxWrlyJ119/HY8++ihKS0sxYcIE1NfXd1hfUlICr9fb+pWXl9fdmyQi0qt0tY8C6qUiEr26fdicMWMGLrnkEowYMQLTp0/Hn//8Z9TU1OAPf/hDh/WLFi1CbW1t61d5eXl3b5KISK/S1T4KqJeKSPSy/BPnqampKC4uxp49ezr8vsfjgcfD3WtKRKQv6qyPAuqlIhK9LL/PZkNDA/bu3YucnByrVyUiEpPUR0WkN+v2dzZvu+02XHDBBcjPz8ehQ4ewePFiOBwOfOtb3+rScoIJNhh35Dv7N2dyd/4/+PUsqi53HZc0VD+Qe/eATQZa/9//H1X3pfu5tA0AMC7u5whXI5dM4Evj9jlQxL0YBhO4c1c6m0voKXjwA6ou1Oyj6mwu7qmRtZZLLgqXRv7jjhP6lXHXqvFySRH+bK4OAI4N5c5x3DEuUiW5nEskah6YSNWxyUWpr+zotCZoWqhl9ZTu6qMA0JjlgMMTOT4m8SB3Tpuyub4Sd5Q7V4FErg+wyUDvLXmUqht7J7c8AEguJ68VMqyrYSDXMxIqubimkJs7J9593XvNGzJ1KpDIRSbVFXD7YTNcylqok/nhhEACVYbUvXzqWGMut89sL60p4l6PQm7uc9ohD3dsOuu5wQA/Qnb7sPnJJ5/gW9/6Fqqrq5GZmYmvfe1r2Lx5MzIzM7t7VSIiMUl9VERiSbcPm88991x3L1JEpE9RHxWRWKJsdBERERGxjIZNEREREbGMhk0RERERsYyGTRERERGxjIZNEREREbGMhk0RERERsYyGTRERERGxjOXZ6F9UQlUQTlfkdBunL3IqxgnxR7mUnEA6l0wQIENZbGFu+9hkoPfveoRbMYAxPyYTMsgUCBsZnuDaV0HV+b9aSNUlH+BiOfznnEHV2QPcjvhTuQSI+jzuHCeXczfjbiQTWtz13HGxcSE+AIC03VzCiKOZW6g/nTuGiftqqbqGc9K5umnDO60JBnzAq9Tier2M7U1wOiNf96Uz46hlublThepRZJwOGbuTXMo9L9hkoHd/xiUNAcAZj3H9mX2utXi5fW7K4hK9nM3ceoNkUg67HyFu85BQwe1v2i5uxUkHGqm6miHcC3X22iNU3Ue39qPqACBtO/e64MvgXoC9e8jUtj9xSXrBMdzrZVM2l3bF0DubIiIiImIZDZsiIiIiYhkNmyIiIiJiGQ2bIiIiImIZDZsiIiIiYhkNmyIiIiJiGQ2bIiIiImIZDZsiIiIiYhkNmyIiIiJiGQ2bIiIiImKZqI2rdNUF4HRGjnyKP8jFVoUSuNi840O4PC93HRe/lVLK5YgZFzfz0xGUALb8lItjG/Jbbpmpu7j1hgZysYye49wxdPi5upoi7hwbBxcPZsgYT1b1MC6+zOnjluep4Y5Lw0D+58nECi4SrTGXy6lL2VNP1QVTuKjEUBx3UlqSO9/nUEvf+Tk7kOKCcXXy/CAveBd3SuE5zh1f9nnGxuUml3ORq2wEJQB8dD0XE1z8JNdLbWSSp4uMpG3uzx3EUBy3PM8xMkJxH3dSEg/7qTrnce710lZ5jKpLjuPGG+PmXjviKvhxKXUvdx3WGG7dnjouytPm5uIl3eXVVJ0vMyfi9+1diEPuOx1XRERERE45DZsiIiIiYhkNmyIiIiJiGQ2bIiIiImIZDZsiIiIiYhkNmyIiIiJiGQ2bIiIiImIZDZsiIiIiYhkNmyIiIiJimahNEPJleOB0RU4qcXi5u+8n7DlO1bma46k6NtXGn84lrbgag1QdupBqwyYD7foulzR0xuNc4kbah1yqhNPHHcOQm9tpNlGnOZNbnr2FWx6bSOThLkHYg9x6wy5uve5aMq4EgC1Anjs/GedCsge55bEpMszz0wT449Lb2YMGdkTeX3cdmSDUyB03Xz/yeRagyuAmn9+d7GYrWxeST9hkoI/ncL106H9zy/NlcMcwoYLb6dpiMvmujKsLerjtqyniEsK8pWQvreWSAx0+7nU1mMa97rPXFgA05HKziZtMibKzPTcvcuLPCSEyXcnVEPmJYgvyTyS9sykiIiIiltGwKSIiIiKW0bApIiIiIpbRsCkiIiIiltGwKSIiIiKW0bApIiIiIpbRsCkiIiIiltGwKSIiIiKW0bApIiIiIpaJ2gSh2kIHHB5HxBqHL/L3T/Ac4RICUrdWUXWBXC9V15jDJQj50rg6NkEFAFJ3cXVsMtBH8x6h6v5z5TepusZs7hiyqUk5/9tE1QWSEqg6V0P3JqV493NRKb407pr2e7mfE+OP8hdNfT6X9JF4mNuXUJKbqnMdqqXqkg5y5y75w+pOa4IhP7WsWBBIdMC4Il9XgWTyem/hrncnF/JCC3OBLGgYyF1zLd4uJGuRpWwy0P9dyyUNfeX2G6i62kLyPSN2lw3ZM2zctZBYwSX5OJq5OtPI9Xp7HPe6Gkjh6rry1pydTCgzdu4YOnxkUk+AO4b2yqNUnSfYP+L3HV3oo3pnU0REREQs0+Vhc+PGjbjggguQm5sLm82Gl156qc33jTG4++67kZOTg/j4eEyZMgW7d+/uru0VEen11EdFpC/p8rDZ2NiIkSNHYvny5R1+/8EHH8TDDz+MFStW4J133kFiYiKmT58On8930hsrIhIL1EdFpC/p8mc2Z8yYgRkzZnT4PWMMfvWrX+Guu+7ChRdeCAB46qmnkJWVhZdeegmXX355u//H7/fD7//X7/3r6uq6ukkiIr1Kd/dRQL1URKJXt35ms7S0FBUVFZgyZUrrY16vF+PGjcOmTZs6/H9KSkrg9Xpbv/Ly8rpzk0REepUv0kcB9VIRiV7dOmxWVFQAALKysto8npWV1fq9z1q0aBFqa2tbv8rLy7tzk0REepUv0kcB9VIRiV49fusjj8cDj4e89YCIiHRIvVREolW3vrOZnZ0NAKisrGzzeGVlZev3RETk86mPikis6dZhs6CgANnZ2VizZk3rY3V1dXjnnXcwfvz47lyViEhMUh8VkVjT5V+jNzQ0YM+ePa3/Li0txbZt25Ceno5Bgwbhlltuwf3334/TTz8dBQUF+MlPfoLc3FzMnDmzS+tJLgvB6Yp81/zET7jbgDTkc8kj4aJEqi7k5u76n/nGPqouUJRD1bn2ff7ntT4rNDCTqkv7kEuLYJOB/vS3V6i6GYVfoeoC44dRdfYWLmEh890aqs4WJFM0qo5RZcHiAVRd4tYyqs40NVN1tn5pVB0ANBVz10zV2dyvalP3cOekLj+r8yIAiYdaqLqmws73ORjwAT1428pT1UcBIOWdMjjtkZN1/CmF1LKC8VwySm0xGVdDlvV/j6tLqOTSrZqy+I8buOrJdKUM7nWBTQba/OAKqm7axXOoutoi7nWwMYdLMQuT04Oxc4UhD/fely2riKqry+fW238LF3cV8nDpVADQnMHtS4gLbUMwnktB7P8+19SqLxlB1TX3j3xNh/w+YAe1qK4Pm++99x7OPffc1n8vXLgQADBnzhysXLkSt99+OxobGzFv3jzU1NTga1/7Gl5//XXExZFHVUQkxqmPikhf0uVhc/LkyTDm83/Ss9lsuPfee3Hvvfee1IaJiMQq9VER6UuUjS4iIiIiltGwKSIiIiKW0bApIiIiIpbRsCkiIiIiltGwKSIiIiKW0bApIiIiIpbRsCkiIiIilunyfTZPFZv59CuSllTujv7e7dVUXTAjiao7MpJLYqg+r4BbbwKXPOH/KpfyAQCe41zqhdPH1TVme6k6Nhnof/ZtpuqGPfJlqi75AJcIknAkSNWF4rhz0vTVdHK9XJrO0RHcOW7O4rbPwYVsAQDcNdy1kLWFW2hLCtde0v5RR9UdGZ1C1SVWdn6swzbu+MWCmnPy4XRFvhl8Yy53PJLKuWuk6I/cNRLycGk1oXgykcXN1Tm5AC4AnaeonJBQQaYrFXLbyCYDvfnHJ6m6Mx67iarrv5Xrke56Lq3J0cQtz3GYS2ODi+sryeR+NI7I5Za3nyoDAHhquAS6uGruGLoauLrmSVziXubbXFJdOCPy634w5Mcuakl6Z1NERERELKRhU0REREQso2FTRERERCyjYVNERERELKNhU0REREQso2FTRERERCyjYVNERERELKNhU0REREQso2FTRERERCwTtQlCgUQ7wp2kQYSdZArI4DSqzNnIpbzEH+PSAbwf1VJ1pbO57Us+wCVUAIDDz9WG3OQxJMsC47kEAzYZ6MObHqHqpl0yl6oLu9iEEe74eY6TCRBHm6i6hINcKlYghas7XszVAUCYvBbq87i0pn5buUSQcAK3jXHHuedd0j+OdFoTDPmpZcWChMM+ODvp9A25XCpaKHIQUavjQ+K5QpLhgobg3ddC1QW53QUAhOLIZKBisj+TZbVF3EayyUAfXc/10ilXfJeqqynkLobkcu6cxDVxCX5m7wGuLsT1C+MYQNU1kaltAJB0iJsl2GQg5/5Kqi7s4tKQjI9L+LIfiXyx2sN8H9U7myIiIiJiGQ2bIiIiImIZDZsiIiIiYhkNmyIiIiJiGQ2bIiIiImIZDZsiIiIiYhkNmyIiIiJiGQ2bIiIiImIZDZsiIiIiYpmoTRBKLvPD2UlCkHsnlySAYJAqOzjnTKoudR+3vBCZjFLw4AdUnf+cM6g6AKgpclF1nhouziLnf7kEHHsLl5yQfIBLoWGTgd5ctZKqG/KbG6k6dy2XFmHnLgW0pHKJKsbOnY8k8tJv8fKpF7kbGqk6fz/u3Bkn97Osva6Zqos/yl3TZRfndFoT8vuApdTiej3jtHd6LsLkK0EwjrueAlwYDMLcKUXmdvKJRrJxbQoA4DnG7XNKGZsgxCXbNOZwsUn9t3LHhk0GevuZ31J1Zz3EJRc5AtxJbsrmkvTsw7m62kKu/+S/Uk2utx9VB/DphqUXJlJ1nmNFVN3AR7hZ4shlI6i62tMjfz/s8wH3UIvSO5siIiIiYh0NmyIiIiJiGQ2bIiIiImIZDZsiIiIiYhkNmyIiIiJiGQ2bIiIiImIZDZsiIiIiYhkNmyIiIiJiGQ2bIiIiImKZqE0QCsY5AFfkBAVPMhlTEeASFlL3cnWJm/dy681Mp8pCzT6qzh7gkicAwDi4BIPmTDYRJIGqy3y3hqpLOMId67CL+3mITQbadc2jVN3gZ2+g6lz13PbZyAAUe4g7H04/dy3EV3HrBQBXaQVV52jkEjxsh45QdeH6BqrOjTyqLnVP54klwUAXImR6OdfWj+G0RU4z82acRS0r5Oauz8aBfHIVw9i45RlytSEuBAsA4N3HPdeCHnLl5L6wqU7u+gBVV1MYR9WxyUA7bn2Eqhv5AJk01MIlMCVUcecjZT9XV3dGKlUXTCQTogCEyGvBu5tbnreUmxFsg3KpOvZ5krIn8vdDLdxyAL2zKSIiIiIW6vKwuXHjRlxwwQXIzc2FzWbDSy+91Ob7c+fOhc1ma/N1/vnnd9f2ioj0euqjItKXdHnYbGxsxMiRI7F8+fLPrTn//PNx+PDh1q9nn332pDZSRCSWqI+KSF/S5c9szpgxAzNmzIhY4/F4kJ2dTS3P7/fD7/e3/ruurq6rmyQi0qt0dx8F1EtFJHpZ8pnN9evXo3///hgyZAhuvPFGVFdXf25tSUkJvF5v61deHvdHACIisawrfRRQLxWR6NXtw+b555+Pp556CmvWrMEDDzyADRs2YMaMGQiFOv7rz0WLFqG2trb1q7y8vLs3SUSkV+lqHwXUS0UkenX7rY8uv/zy1v8+66yzMGLECBQVFWH9+vU477zz2tV7PB54PF24D4WISIzrah8F1EtFJHpZfuujwsJCZGRkYM+eTm7YJCIiHVIfFZHezPJh85NPPkF1dTVycnKsXpWISExSHxWR3qzLv0ZvaGho89N1aWkptm3bhvT0dKSnp2PJkiW46KKLkJ2djb179+L222/H4MGDMX369C6tJ66qGU5HJwkAhrujv0lJpOrYhB6Tncktr5ZLRrG5uNPgT+08GeUENiHATqY2uBq4OluQO4ahODIpp5lbr7uWWx6bDLTnWyuoumGPcOkYNjKwxl3P1nHH2R7gUy/g4q6vUDKXROLsxyUN2cmkrXAi9yvikKvzayGM7k246apT1UcBoHr2WXC4I58zNr3FRT4f46vIZC0yFM2Xyi0vkMhdwwkV/PMi8bC/8yIANUXc8yKxgosTM3budcHRxC0vuZyLe3EEuGPIJgN9cAeXNFS8kkuBY98js5OpbZ5DXHP27uHfm6s7jatlX6eDCdy1lUTOCL4MbsVxR7vw+tGJLg+b7733Hs4999zWfy9cuBAAMGfOHDz66KPYvn07nnzySdTU1CA3NxfTpk3Dfffdp88SiYj8k/qoiPQlXR42J0+eDBPhHcU33njjpDZIRCTWqY+KSF+ibHQRERERsYyGTRERERGxjIZNEREREbGMhk0RERERsYyGTRERERGxjIZNEREREbGMhk0RERERsUyX77N5qtQMS+409SLlgJtalmvnAarOnOal6oKp3N38j53DJahkreX2oz7PQdV1hXFwSQK+fmTUQdUxqqzpq1xqjOd4gKpj0yJc9dzPV2wy0Ic3cekY7PKasrjEhkAS99R1NlJlAIC4Su76b/FyKRXGlkTVuXYdpOqaiguoupC782s11MMJQqeSI2DgsEW+rpIOcukyjTlcr3I2cdcxnSTm45ZXV8A9v9N2kZFeAJzHm6k6bym3L45mrlmFPNy+OA5zPTeuiXs+NmVzr1ts6hSbDPTx3EepujMe797UNnYMiq8m464ABBK412obmYJIJ/4c4+r8/chEwFDk5YX8fB/VO5siIiIiYhkNmyIiIiJiGQ2bIiIiImIZDZsiIiIiYhkNmyIiIiJiGQ2bIiIiImIZDZsiIiIiYhkNmyIiIiJiGQ2bIiIiImKZqE0QSi73wdnJ1tk3bqOWZeLjqbqGAdzhSCRTd/rtaKLqwqVlVF1yeSZVBwDVw7gEA89xbnne/VyST7B4AFWXcISLd3Ad5Y5hSyp3jm1k0hCbPtHdSUNnrCDTMcgwCz5FAwglkOkwzdxC3QeOUnX1X+WSgZoyuWs6kEwkCHUh+aK3czYbOIORE0M8eyqpZYXicqk6dwN3fAMJ3PsdSeVcio/NcH0g6QAfrWWr5BJ6PLXcMk0j19NsWUVUHVzc65bZyyXp2YdzCUIJVWyiDneO2WSgj+ZxvfT0p7jkInc9l6aTtK+eqgOA48Wp3LprueXZuZdfBOK5513Yye1zMLGTBKEuhBrqnU0RERERsYyGTRERERGxjIZNEREREbGMhk0RERERsYyGTRERERGxjIZNEREREbGMhk0RERERsYyGTRERERGxjIZNEREREbFM1CYIBeKdMJ0kI8SfOYRalu0gl47hauDuqu9L5W6b7/dyaRb9yrKousZs/mcDp4+rs3eSLHKCL43b58StXBrS0RGFVF3CQS7Vxti5/bCHuIQFNxkW0ZTFrZdNBvroBi4dg00uQgJXBgAhD3d9xR1uoOoaRnJpM0kfcklD/vHk88TbeU2IfH7EBPPPrwgaR3Dnyk/2vuYM7loKJFJlcDV5qLqQm3t+1wxJ4lYMIDmOe5l0+Lh4Mnscty91+dx6k7dy6zUhLvGntpA7dyn7ueXZuzm1jU0G2n3Vo1Td1NfmUnX1RclUHQDYW7i6hjyuLhTPHetBL1Rx683neqkvJ/LJCzeTJxd6Z1NERERELKRhU0REREQso2FTRERERCyjYVNERERELKNhU0REREQso2FTRERERCyjYVNERERELKNhU0REREQso2FTRERERCwTtQlCLV4Hwq7IaRUJB7i76oeOH6fq4o5zEQYtydyM3pDLpW0YL5dm4a7n0moAwFPD1YZdXOKG38vts2lqpuqas7j1BlK4BKGkA1QZnH7umnHXc3WBJO4pZOMWRycDfXgTmTS0nEwaAuDr56Lq4srYSBCurGVAKlXX1J9bYFJZ59d+qIV/LvV27rognM7I58xTyaVCOZu4yJ+ESu5chZ1cnSHrAmRiVvbaI1whAOPmnhfBNC4xLpDCJQj139JI1bHpT8YxgKrLf6Waqqs7I5Wq8xwio4HIcYR9HWSTgd76w0qq7rzvXEPVAUBmDdcj4/Zy12FgQDpVVzWZSwZK3cUdw/i/Ra4LBgw+oZakdzZFRERExEJdGjZLSkowZswYJCcno3///pg5cyZ27drVpsbn82H+/Pno168fkpKScNFFF6GykssmFxGJdeqjItLXdGnY3LBhA+bPn4/NmzfjrbfeQiAQwLRp09DY+K+3+2+99Va8+uqrWLVqFTZs2IBDhw5h9uzZ3b7hIiK9kfqoiPQ1XfrM5uuvv97m3ytXrkT//v2xdetWTJw4EbW1tfjNb36DZ555Bv/xH/8BAHjiiSdwxhlnYPPmzfjKV77SfVsuItILqY+KSF9zUp/ZrK2tBQCkp3/64dWtW7ciEAhgypQprTVDhw7FoEGDsGnTpg6X4ff7UVdX1+ZLRKSv6I4+CqiXikj0+sLDZjgcxi233IJzzjkHw4cPBwBUVFTA7XYjNTW1TW1WVhYqKio6XE5JSQm8Xm/rV15e3hfdJBGRXqW7+iigXioi0esLD5vz58/Hzp078dxzz53UBixatAi1tbWtX+Xl5Se1PBGR3qK7+iigXioi0esL3WdzwYIFeO2117Bx40YMHDiw9fHs7Gy0tLSgpqamzU/llZWVyM7O7nBZHo8HHg933zERkVjRnX0UUC8VkejVpXc2jTFYsGABVq9ejbVr16KgoKDN90ePHg2Xy4U1a9a0PrZr1y6UlZVh/Pjx3bPFIiK9mPqoiPQ1XXpnc/78+XjmmWfw8ssvIzk5ufXzQ16vF/Hx8fB6vbjmmmuwcOFCpKenIyUlBTfffDPGjx/f5b+gdPrCcIY6iV0JcHfpd6SkUHVHR3BJEQmHubvvuxq4On82lyBkY4MYADQM5H6OcNeSSQJHuQgcW780qs7ho8pwvJhLEGrxcgkj8VXceu0B7rg4uZAP/tyRCShsMtCH87mkIQCYdP08qq65kEuziD/EHRx7HZc65RrKpWMEkjq/FkJ+Mt7IAqeyjwJA3N4qOO2R3/Hcd82gL7Qvn8eXRaZMkUFOCeXcS1XqXq5PfXRrP27FAOIqyJdJNpSKfIsn5OF6X/J+bnlNZGqbfTh3bIKJ3A5793A7HF/NnbukffVUXX1RMlXHJgOtefo3VB0AjLvjRqquavTAzosApO7mXkAyf/8BVeebeCZV15wROQUx1MKlJAJdHDYfffRRAMDkyZPbPP7EE09g7ty5AICHHnoIdrsdF110Efx+P6ZPn45HHuFf8EREYpn6qIj0NV0aNo3p/CeZuLg4LF++HMuXL//CGyUiEqvUR0Wkr1E2uoiIiIhYRsOmiIiIiFhGw6aIiIiIWEbDpoiIiIhYRsOmiIiIiFhGw6aIiIiIWEbDpoiIiIhYRsOmiIiIiFimSzd1P5UasxxwuCNHIbUk9aeW5a7PoOryXjtK1R0+l1te0mEuYurY0MhRciek7W6h6gAgsYKMlwxwdfX5cVRdU3EmVeeu4aLOwm4uYi13AxeN6CqtoOrg4qJL4yq9VF0ogYueC3m4n/98/bjtYyMoAWDDY49TdWPu4qLYGnK4mFjjJOsc3LXQf2tTpzXBIJmXGgOqJwyEwx35+Zv+jy5k4RIC+7kYO38qd04ztvupusZc7nmRtp2P2Uvdy/XdBnLdbBRucwbXCzw1ZMzjIe4ch53cOQl5uLq607j9CCRw5+R4cSpVZydfLjNruGhVNoISAN554FGqbvjDXOzwsaHcsWlJGknVNWVz566zWOmQn3+/Uu9sioiIiIhlNGyKiIiIiGU0bIqIiIiIZTRsioiIiIhlNGyKiIiIiGU0bIqIiIiIZTRsioiIiIhlNGyKiIiIiGU0bIqIiIiIZaI2QchmPv2KxOHnkhjctVxCwPFR6VSdPcitl02KiDvGJUA4mvmUj8ZcLpXI6efWnXg4QNVVnc2tN2sLl+BSn8ctz9+Pq3M0plF1oWQuManFy6WGOMlzF3e4gasr467p5kLumgb4ZKAt93PpGF9deANVZwtxzxNj51IvWlI7PyfBQPcm5kQzW9jAFo58jI8PIRN1uFMAQy4u7ODO/bFh3POb7aW+DHJHANQY7jnuru/e6zjEtSDEVXO92dXA1ZVemEjVeXdTZTDkobYZ9vWcW15DHlcXt/cIVVc1eiC3QPDJQDu/9whVN+Q3XG9Of+EDqi58OZc0VJ8f+eSFuhDEpnc2RURERMQyGjZFRERExDIaNkVERETEMho2RURERMQyGjZFRERExDIaNkVERETEMho2RURERMQyGjZFRERExDIaNkVERETEMlGbIORqNHB0ksDT3I+blUNuN1WX/v4xqu6T8zOoutoC7vAml3NpJv50LskCAFL21NO1jFASdwxT93D70pLCHZt+W7lzYpzctWA7xKVFOPtxSUPGlkTVuQ8cpeoaRuZSdWySS/yhRq4QQENOClXHJgP9bekKqu5r37ueqvOlcjttI+JrggEy4iYGNGfY4fBEfn5493HJO6ab357we7kFppCJWTVFXF/x7uH2FwA8dVxPs5NpbA4ft7xgfDxVxyYDOfdXUnWeY0VUnbeUi48JJnBRSGyqk53bXYTiufMRGMClrKXu5lPHjg3l+gubDLTrGi617bx111B1tYOpMoTdkY9h2ME/j/TOpoiIiIhYRsOmiIiIiFhGw6aIiIiIWEbDpoiIiIhYRsOmiIiIiFhGw6aIiIiIWEbDpoiIiIhYRsOmiIiIiFhGw6aIiIiIWCZqE4SSD/jg7GTr7E1clEBDIZfy0lCcStW5ayMnG52Q8fc6qq55YCJVl7ivlqoDgGAKl9pgD3IJAK5D3Lrr8rOourR/cMcmnMAlF9nrmrnl1Tdwy8vkUiVcuw5SdfVfLaDqkj7kkoZaBqRSdexxAQDj5BKEbCHu+meTgf768GNU3ZnLbqLqUvZ33hfsQS6RJhbEHTNwuCOfs5TntlDLssdzfSXcyCVXOVLIay6Z6+Ehdx5Vl/ynD6g6ALCRCXTIy+HqAty11//93VRd86RhVF3YxaWTDXyEOza2QdzyklK55Lu4Y1yCUCCeqxv0QhVVVzWZe83K/D1/zbQkjaTq0l/glskmA615+jdU3Zd+xvXS7I3HI34/GPLjALWkLr6zWVJSgjFjxiA5ORn9+/fHzJkzsWvXrjY1kydPhs1ma/N1ww1cvJ2ISKxTHxWRvqZLw+aGDRswf/58bN68GW+99RYCgQCmTZuGxs/8FHvdddfh8OHDrV8PPvhgt260iEhvpT4qIn1Nl36N/vrrr7f598qVK9G/f39s3boVEydObH08ISEB2dnZ3bOFIiIxRH1URPqak/oDodraTz/Hl57e9vNtv//975GRkYHhw4dj0aJFaGpq+txl+P1+1NXVtfkSEekruqOPAuqlIhK9vvAfCIXDYdxyyy0455xzMHz48NbHr7jiCuTn5yM3Nxfbt2/HHXfcgV27duHFF1/scDklJSVYsmTJF90MEZFeq7v6KKBeKiLR6wsPm/Pnz8fOnTvx17/+tc3j8+bNa/3vs846Czk5OTjvvPOwd+9eFBUVtVvOokWLsHDhwtZ/19XVIS+P+4tCEZHerLv6KKBeKiLR6wsNmwsWLMBrr72GjRs3YuDAgRFrx40bBwDYs2dPh03S4/HA4/F8kc0QEem1urOPAuqlIhK9ujRsGmNw8803Y/Xq1Vi/fj0KCjq/d+C2bdsAADk55D3IRERimPqoiPQ1XRo258+fj2eeeQYvv/wykpOTUVFRAQDwer2Ij4/H3r178cwzz+DrX/86+vXrh+3bt+PWW2/FxIkTMWLECEt2QESkN1EfFZG+xmaM4eJAANhsHd+5/4knnsDcuXNRXl6O73znO9i5cycaGxuRl5eHWbNm4a677kIKmRRRV1cHr9eLr5x/L5yuyGkV7louQcjRzNUdOzOZqqvr+LdY7aTs4+rsAe4UBBK55AQACMVxtTYuQAhJB0NUnaueS8eoLeRSOeKOcxsYf5Q7x+5DZHJRMpeU0jQggavLdFB1jhaqDE39ufPraqCf3jAObpnxR7lz4k/lludL5+r+cfMjVN3Zd9/YaU2oxYftT/4YtbW1dG/qLqeijwL/6qXnnLsYTmfk6zl+L5dcFcrg1n/8DC7xJxhPlaHfjsh/iX9C/WncApPK/dyKAbjLq6m6UD/u9cN+oIKqO/qNYqou8+0yqs74fFTdkW9y6zXky5Evgyv09+N6VdjJ1bkauPWm7uKWF1/NvQYCwNGzyNSkam7dtYO59SYe4vb5/Tu5Xlr8ZOReGvb5UHoP10e7/Gv0SPLy8rBhw4auLFJEpE9RHxWRvuak7rMpIiIiIhKJhk0RERERsYyGTRERERGxjIZNEREREbGMhk0RERERsYyGTRERERGxjIZNEREREbGMhk0RERERsUyXbup+KsUfaYbTEfnmx44jtdSywlVcOoZv0iiqzs2tFg4/l7SS+soOqq5h2nBuxQBakrmfIxx+LsEg+UMuRaOpMI2qS6zk0hiS/nGEqiu7mMuMTt3DJTuEXFwSQ8jN1QWSubpGL1WGpDIydSqJT53qv5VLaWlJ5Y6hzXCpSSn7ufQnJhkIAN6799FOa+rqw0h7klpcn1A/Mouqa0ni+gqbChXkArgQiuNeqtg0tqZsLsEMAHyZXG9xNXA9zRPsT9U1kylh4QyuadiPkGk1p1NlSNnD1cUd5dZrC5HXDJmk58vh0uzi/8ZtX3MG188AwMGFNaE+nzzHbm6WyN54nKorHsD10o/nRO6ldfVhpN1DLUrvbIqIiIiIdTRsioiIiIhlNGyKiIiIiGU0bIqIiIiIZTRsioiIiIhlNGyKiIiIiGU0bIqIiIiIZTRsioiIiIhlou6m7sZ8eoPVYMjfeW248xoACJsWqi7kJ+/ESgoFuBuxBsntCwb47Qu1cD9HGPImyMz5APhtDNvIG/iS62XPXTDA3Xg5DPKm7mydn6wjT3GohTtv7HoBIBjs3mMYDHA3QbYHuZsvh1q49dbVd/68q2v4tOZEv4lFrb002PlziD1XbF8J+ck68j7ZzD4AQDDQ/S9pdu6ygy3IFTq6u6eRy7Ozr5c+br0h7mWLRvdI8poJN3N9JUi+BoZa+Ju609c/2e/DDnKWIK8F9hx31ku70kdtJsq67SeffIK8vLye3gwR6QPKy8sxcODAnt4MS6iXisipwPTRqBs2w+EwDh06hOTkZNj++e5XXV0d8vLyUF5ejpSUlB7ewpMTK/ui/Yg+sbIvp2I/jDGor69Hbm4u7PbY/DRRLPdS7Uf0iZV90X7wutJHo+7X6Ha7/XMn5JSUlF598v9drOyL9iP6xMq+WL0fXi8ZRN9L9YVeqv2IPrGyL9oPDttHY/NHehERERGJCho2RURERMQyvWLY9Hg8WLx4MTweT09vykmLlX3RfkSfWNmXWNmPaBQrx1b7EX1iZV+0H9aIuj8QEhEREZHY0Sve2RQRERGR3knDpoiIiIhYRsOmiIiIiFhGw6aIiIiIWEbDpoiIiIhYplcMm8uXL8dpp52GuLg4jBs3Du+++25Pb1KX3HPPPbDZbG2+hg4d2tObRdm4cSMuuOAC5Obmwmaz4aWXXmrzfWMM7r77buTk5CA+Ph5TpkzB7t27e2ZjI+hsP+bOndvuHJ1//vk9s7ERlJSUYMyYMUhOTkb//v0xc+ZM7Nq1q02Nz+fD/Pnz0a9fPyQlJeGiiy5CZWVlD21xx5j9mDx5crtzcsMNN/TQFvd+vb2PAr23l6qPRhf10VPfR6N+2Hz++eexcOFCLF68GH//+98xcuRITJ8+HVVVVT29aV1y5pln4vDhw61ff/3rX3t6kyiNjY0YOXIkli9f3uH3H3zwQTz88MNYsWIF3nnnHSQmJmL69Onw+XyneEsj62w/AOD8889vc46effbZU7iFnA0bNmD+/PnYvHkz3nrrLQQCAUybNg2NjY2tNbfeeiteffVVrFq1Chs2bMChQ4cwe/bsHtzq9pj9AIDrrruuzTl58MEHe2iLe7dY6aNA7+yl6qPRRX20B/qoiXJjx4418+fPb/13KBQyubm5pqSkpAe3qmsWL15sRo4c2dObcdIAmNWrV7f+OxwOm+zsbPOLX/yi9bGamhrj8XjMs88+2wNbyPnsfhhjzJw5c8yFF17YI9tzMqqqqgwAs2HDBmPMp8ff5XKZVatWtdZ89NFHBoDZtGlTT21mpz67H8YYM2nSJPP973+/5zYqhsRCHzUmNnqp+mj0UR+1XlS/s9nS0oKtW7diypQprY/Z7XZMmTIFmzZt6sEt67rdu3cjNzcXhYWF+Pa3v42ysrKe3qSTVlpaioqKijbnx+v1Yty4cb3u/ADA+vXr0b9/fwwZMgQ33ngjqqure3qTOlVbWwsASE9PBwBs3boVgUCgzTkZOnQoBg0aFNXn5LP7ccLvf/97ZGRkYPjw4Vi0aBGampp6YvN6tVjqo0Ds9VL10Z6nPmo95ylfYxccPXoUoVAIWVlZbR7PysrC//3f//XQVnXduHHjsHLlSgwZMgSHDx/GkiVLMGHCBOzcuRPJyck9vXlfWEVFBQB0eH5OfK+3OP/88zF79mwUFBRg7969uPPOOzFjxgxs2rQJDoejpzevQ+FwGLfccgvOOeccDB8+HMCn58TtdiM1NbVNbTSfk472AwCuuOIK5OfnIzc3F9u3b8cdd9yBXbt24cUXX+zBre19YqWPArHZS9VHe5b66KkR1cNmrJgxY0brf48YMQLjxo1Dfn4+/vCHP+Caa67pwS2TEy6//PLW/z7rrLMwYsQIFBUVYf369TjvvPN6cMs+3/z587Fz585e8Zm1SD5vP+bNm9f632eddRZycnJw3nnnYe/evSgqKjrVmylRQL00uqmP9pxo76NR/Wv0jIwMOByOdn8BVllZiezs7B7aqpOXmpqK4uJi7Nmzp6c35aScOAexdn4AoLCwEBkZGVF7jhYsWIDXXnsN69atw8CBA1sfz87ORktLC2pqatrUR+s5+bz96Mi4ceMAIGrPSbSK1T4KxEYvVR/tOeqjp05UD5tutxujR4/GmjVrWh8Lh8NYs2YNxo8f34NbdnIaGhqwd+9e5OTk9PSmnJSCggJkZ2e3OT91dXV45513evX5AYBPPvkE1dXVUXeOjDFYsGABVq9ejbVr16KgoKDN90ePHg2Xy9XmnOzatQtlZWVRdU4624+ObNu2DQCi7pxEu1jto0Bs9FL10VNPfbQH+mjP/n1S55577jnj8XjMypUrzYcffmjmzZtnUlNTTUVFRU9vGu0HP/iBWb9+vSktLTX/+7//a6ZMmWIyMjJMVVVVT29ap+rr6837779v3n//fQPALF261Lz//vvmwIEDxhhjfv7zn5vU1FTz8ssvm+3bt5sLL7zQFBQUmObm5h7e8rYi7Ud9fb257bbbzKZNm0xpaal5++23zZe//GVz+umnG5/P19Ob3saNN95ovF6vWb9+vTl8+HDrV1NTU2vNDTfcYAYNGmTWrl1r3nvvPTN+/Hgzfvz4Htzq9jrbjz179ph7773XvPfee6a0tNS8/PLLprCw0EycOLGHt7x3ioU+akzv7aXqo+qjVuhNfTTqh01jjFm2bJkZNGiQcbvdZuzYsWbz5s09vUldctlll5mcnBzjdrvNgAEDzGWXXWb27NnT05tFWbdunQHQ7mvOnDnGmE9v2/GTn/zEZGVlGY/HY8477zyza9eunt3oDkTaj6amJjNt2jSTmZlpXC6Xyc/PN9ddd11UvhB3tA8AzBNPPNFa09zcbG666SaTlpZmEhISzKxZs8zhw4d7bqM70Nl+lJWVmYkTJ5r09HTj8XjM4MGDzQ9/+ENTW1vbsxvei/X2PmpM7+2l6qPRRX301PdR2z83WERERESk20X1ZzZFREREpHfTsCkiIiIiltGwKSIiIiKW0bApIiIiIpbRsCkiIiIiltGwKSIiIiKW0bApIiIiIpbRsCkiIiIiltGwKSIiIiKW0bApIiIiIpbRsCkiIiIilvn/AQDF2SRtu52xAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 800x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "labels_str = [\"Healthy\", \"Schizophrenic\"]\n",
    "\n",
    "fig = plt.figure(figsize=(8, 4))\n",
    "\n",
    "ax = fig.add_subplot(121)\n",
    "imgplot = ax.imshow(data[0])\n",
    "ax.set_title(labels_str[labels[0]])\n",
    "\n",
    "ax = fig.add_subplot(122)\n",
    "imgplot = ax.imshow(data[1])\n",
    "ax.set_title(labels_str[labels[1]])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to compare with a standard Euclidean method, we also flatten the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(86, 378)\n"
     ]
    }
   ],
   "source": [
    "flat_data, _, _ = data_utils.load_connectomes(as_vectors=True)\n",
    "print(flat_data.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Manifold"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentionned above, correlation matrices are SPD matrices. Because multiple metrics could be used on SPD matrices, we also import two of the most commonly used ones: the Log-Euclidean metric and the Affine-Invariant metric [[PFA2006]](#References). We can use the SPD module from `geomstats` to handle all the geometry, and check that our data indeed belongs to the manifold of SPD matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "from geomstats.geometry.spd_matrices import (\n",
    "    SPDMatrices,\n",
    "    SPDAffineMetric,\n",
    "    SPDLogEuclideanMetric,\n",
    ")\n",
    "\n",
    "manifold = SPDMatrices(28, equip=False)\n",
    "\n",
    "spd_ai = SPDMatrices(28, equip=False)\n",
    "spd_ai.equip_with_metric(SPDAffineMetric)\n",
    "\n",
    "spd_le = SPDMatrices(28, equip=False)\n",
    "spd_le.equip_with_metric(SPDLogEuclideanMetric)\n",
    "\n",
    "print(gs.all(manifold.belongs(data)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Transformer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great! Now, although the sum of two SPD matrices is an SPD matrix, their difference or their linear combination with non-positive weights are not necessarily! Therefore we need to work in a tangent space to perform simple machine learning. But worry not, all the geometry is handled by geomstats, thanks to the preprocessing module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geomstats.learning.preprocessing import ToTangentSpace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What `ToTangentSpace` does is simple: it computes the Frechet Mean of the data set (covered in the previous tutorial), then takes the log of each data point from the mean. This results in a set of tangent vectors, and in the case of the SPD manifold, these are simply symmetric matrices. It then squeezes them to a 1d-vector of size `dim = 28 * (28 + 1) / 2`, and thus outputs an array of shape `[n_patients, dim]`, which can be fed to your favorite scikit-learn algorithm.\n",
    "\n",
    "Because the mean of the input data is computed, `ToTangentSpace` should be used in a pipeline (as e.g. scikit-learn's `StandardScaler`) not to leak information from the test set at train time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.model_selection import cross_validate\n",
    "\n",
    "pipeline = Pipeline(\n",
    "    steps=[\n",
    "        (\"feature_ext\", ToTangentSpace(space=spd_ai)),\n",
    "        (\"classifier\", LogisticRegression(C=2)),\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have all the material to classify connectomes, and we evaluate the model with cross validation. With the affine-invariant metric we obtain:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7098039215686274\n"
     ]
    }
   ],
   "source": [
    "result = cross_validate(pipeline, data, labels)\n",
    "print(result[\"test_score\"].mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And with the log-Euclidean metric:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.6862745098039216\n"
     ]
    }
   ],
   "source": [
    "pipeline = Pipeline(\n",
    "    steps=[\n",
    "        (\"feature_ext\", ToTangentSpace(space=spd_le)),\n",
    "        (\"classifier\", LogisticRegression(C=2)),\n",
    "    ]\n",
    ")\n",
    "\n",
    "result = cross_validate(pipeline, data, labels)\n",
    "print(result[\"test_score\"].mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But wait, why do the results depend on the metric used? You may remember from the previous notebooks that the Riemannian metric defines the notion of geodesics and distance on the manifold. Both notions are used to compute the Frechet Mean and the logarithms, so changing the metric changes the results, and some metrics may be more suitable than others for different applications.\n",
    "\n",
    "We can finally compare to a standard Euclidean logistic regression on the flattened data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7333333333333334\n"
     ]
    }
   ],
   "source": [
    "flat_result = cross_validate(LogisticRegression(), flat_data, labels)\n",
    "print(flat_result[\"test_score\"].mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example using Riemannian geometry does not make a big difference compared to applying logistic regression in the ambiant Euclidean space, but there are published results that show how useful geometry can be with this type of data (e.g [[NDV2014]](#References), [[WAZ2918]](#References)). We saw how to use the representation of points on the manifold as tangent vectors at a reference point to fit any machine learning algorithm, and compared the effect of different metrics on the space of symmetric positive-definite matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ".. [PFA2006] Pennec, X., Fillard, P. & Ayache, N. A Riemannian Framework for Tensor Computing. Int J Comput Vision 66, 41–66 (2006). https://doi.org/10.1007/s11263-005-3222-z\n",
    "\n",
    ".. [NDV2014] Bernard Ng, Martin Dressler, Gaël Varoquaux, Jean-Baptiste Poline, Michael Greicius, et al.. Transport on Riemannian Manifold for Functional Connectivity-based Classification. MICCAI - 17th International Conference on Medical Image Computing and Computer Assisted Intervention, Polina Golland, Sep 2014, Boston, United States. hal-01058521\n",
    "\n",
    ".. [WAZ2918] Wong E., Anderson J.S., Zielinski B.A., Fletcher P.T. (2018) Riemannian Regression and Classification Models of Brain Networks Applied to Autism. In: Wu G., Rekik I., Schirmer M., Chung A., Munsell B. (eds) Connectomics in NeuroImaging. CNI 2018. Lecture Notes in Computer Science, vol 11083. Springer, Cham"
   ]
  }
 ],
 "metadata": {
  "backends": [
   "numpy",
   "autograd",
   "pytorch"
  ],
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}