5 Neural Networks.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 🧠 Neural Networks Learning\n",
    "\n",
    "In this notebook, you will implement the backpropagation algorithm for neural networks and apply it to the task of hand-written digit recognition. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5000, 401)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# données\n",
    "data = np.genfromtxt('data/neural_nets/data.csv', delimiter=',', dtype=float)\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans ces données (data), les 400 premières colonnes representent les pixels de l'image (20x20), la dernière colonne represente la classe de l'image (chiffres de 0 à 9). (http://yann.lecun.com/exdb/mnist/)\n",
    "\n",
    "Chaque ligne represente un exemple de notre ensemble de données. \n",
    "\n",
    "Mettons ces données dans leus vecteurs correspondants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rajoutons l'ordonnée à l'origine theta 0\n",
    "intercept=np.ones((data.shape[0],1))\n",
    "X=np.column_stack((intercept,data[:,:-1]))\n",
    "y=data[:,-1]\n",
    "# forcer y à avoir une seule colonne\n",
    "y = y.reshape( y.shape[0], 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X (5000, 401)  y  (5000, 1)\n"
     ]
    }
   ],
   "source": [
    "print('X', X.shape ,' y ', y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualisation aléatoire de quelques données"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAABlCAYAAADJa5vzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAB6D0lEQVR4nO39eXBk6ZYfhv1u7vuOfd9r6a5eq1+/mdeP4mhIDxkSR5QsiZQtjSRGUCGGImTJWoZSOIbySlsOhx0hhe1RkMGRKJGyNSGTMmfM0bx5M+/Ne9Ovu193Naq6UAsKKKyJRO57IjOR138Av1Nf3k5UASgsF6j7i0AUCkBm3nvu953vLL9zjqbrOixYsGDBggULFixYsGDBwuXCdtkXYMGCBQsWLFiwYMGCBQsWLOfMggULFixYsGDBggULFkwByzmzYMGCBQsWLFiwYMGCBRPAcs4sWLBgwYIFCxYsWLBgwQSwnDMLFixYsGDBggULFixYMAEs58yCBQsWLFiwYMGCBQsWTADLObNgwYIFCxYsWLBgwYIFE+DCnTNN055rmvbLx/xbXdO02VN+zpGvPbyGuqZp/9Up3/s/0TStevgZjlO8/jrI4O8cvn7zlK+/DjL4K5qmVV7z+q6DHKz9YO0Hax3g9deBcg1vtBwsGVwbGVg6wZLBa8tAuYY3Rg5vcubsn9Z1/V/mfzRN+99omnZf07S2pml/42Uv1HX9NwDcPu8LvAAYZcCFVzn8+r2jXqjr+r8K4M9dxEWeM4wyeFfTtB9rmlbUNG1T07T/1VEv1HX9b+m6HriYyzx3GOXwC5qmfaZpWlnTtEVN07531Auv8X6Y1DTth5qm1TRNe/Syg+E67gdN08YVXcAvXdO0/2WvF1rr4FqtA8A6HwBLBkC3TujXNO3vaZq2fXhG/kTTtO8c9cJrrBN+qGlaWtO0kqZpX2ua9qtHvfA6yuANPhuAC1oLb7JzZsQygP8AwD+67Au5ZPzTuq4HDr/+7GVfzCXgvwHwIwAxAH8KwF/TNO0vXO4lXSw0TYsB+B8A/KcAIgD+TwD+B03Topd5XZeAvwfgKwBxAP8xgP9O07S+y72ki4Ou6+uKLggAeBtAB8BvX/KlXTTe6HVgwJt+PgBvtgwCAD4H8AEOzsjfAvCPNE27LkHK4+LfBjCk63oIwF8F8Hc1TRu65Gu6MFhnQxfOZS1cqnOmadpHmqb9iaZpBU3Tkpqm/WeaprkMf/bnNU1b0TQto2naf6ppmk15/b+uadqSpml5TdP+saZpE6e9Fl3Xf0vX9d8FUD7te5wGZpLBZcFkMpgE8F/rur6v6/ozAH+MC4r4mEgOvwBgR9f1//ehHP4ugDSAf/aU73dsmEUGmqbNA3gfwG/oul7Xdf23AdwH8M+d+uaO/9mmkEEP/CsAfqTr+vMzer8jYRYZXOY6OPx8U8jhMmHJwDwy0HV9Rdf1/4uu68nDs+E3AbgALLzG7R0LZpEBAOi6vqjrepv/BeAEMHba9zsuzCQDAy7sbADMJYfzWguXnTnbB/DvAEgA+C6AfxLAXzP8zV8E8CEODshfBfCvA4B2kDr8j3BgMPYB+DEOIpxXDWaTwX+tHaRof0/TtHde872OCzPJ4P8K4F/RNM2padrC4fX8/mu830lgJjloPf7/1mu833FhFhncBrCi67oarPkaF+Oom0UGAk3TNBwcwL/1uu91TJhFBpe5DgDzyIF4088HwJKBQNO0d3HgnC2fxfu9AqaSgaZp/19N0xoAfgbgDwF88Trvd0yYSgaH73vRZwNgMjmcy1rQdf1CvwA8B/DLR/zufwHgv1f+rwP4FeX/fw3ADw6//10Af0X5nQ1ADcCE8trZU1zD3wXwN45xH5OHn+G4LjIA8IsAvAB8AP46gB0AkZfcxz8BYPM6rQMcZI2WAbQPX/ufHONejvyMqygHHNC3CgD+Mg6iQL+GA8rC//NN2Q8A/mUAnxr+5n8H4O+8SftB+d0nACoAAq+4D2sdvMY6MKscDn9mnQ+WDNTfhXCQRf7rr7iPSVwjnWD4nRMHdYX/7hssg3M/G66IHM50LVw2rXH+0OPc0TStBOB/jwNPWMWG8v0agOHD7ycA/N8O05oFADkcRPdHzvmyzxRmkoGu6z/RD6g7NV3X/w84MNA/Oc17nQRmkYF2UGv1/wPwvwbgwUFq+n+iaZoxInMuMIscdF3P4iDS9O8CSAH4FRxkD0/dge64MIsMcHDYhAw/C+ECaM8mkoGKXwPw27quV17zfY4FE8ng0tYBYCo5vPHnA2DJQLkeLw7qkj89lMO5w2wyAABd11v6QTnMn9UuoDbdjDLABZ8NgDnlcNZr4bJpjf93AI8AzOkHxXT/Eb5Np1K5m+MAtg+/3wDwb+i6HlG+vLqu//Tcr/psYWYZ6D2u5TxgFhlMA9jXdf2/1HW9rev6JoC/D+DPn+K9TgOzyAG6rv+Rrut3dV2P4SB7cAPAZ6d5rxPCLDL4BsC0pmlB5WfvHP78vGEWGQAQQ+yfx8XSVswig8tcB4B55NALb9r50AtvnAw0TXMD+P/gIFj3b5zmPU4J08igBxwAZs7ovV4GU8ngks4GwGRyMOBM1sJlO2dBACUAFU3TbgD4N3v8zb+vaVpU07QxHHRF+W8Pf/7/APDXNU27DQCapoU1TfvnT3sh2kGNkQcHMnFomubRNM1+2vc7AUwhA+2gNeovaprmOrz3fx8HkYifnOb9TghTyADAk4O30P4lTdNsmqYNAvgXASye8v1OCrPIAZqmvXe4J0IA/s8ANnRd/8enfb8TwBQy0HX9CYB7AH7jcD/8RQB3cDHdqEwhAwV/EUAewA9f831OAlPI4JLXAWASOVjngyWDw9c6Afx3AOoAfk3X9c5p3ueUMIsMbmia9uc0TfMenpH/cwDfB/BHp3m/E8IUMlBwGWcDYBI5nOdauGzn7N8D8C/hgCLyX+CF8FT8AwA/x8EB+Y8A/C0A0HX9vwfwfwTw97WDtOYDvN5Mkf8CBwrnL+OgXXIdBxkDaJr2iaZp55WyNYsMgjiIRuQBbOGAyvbnDilu0DTtf6Zp2nlFi00hA13XSzgoEv13cCCHe4fv978FumZ7jJ/m/Y8BU8jhEP8BgAwOokxDOFDCAN6Y/QAAfwkHBcV5AH8TwP9U1/U08GbsBwW/BuC/0g8J84S1Ds59HQDmkcMbfz7AkgFwUJP9TwH4swAK2osZV58Ab4xO0AD8DQC7OOhi/G8D+Bd1Xf8SeGNkQFzG2QCYRw7ntxb0UxTmXfUvAI9x4HX/1ilf/xsAigAaAOyXfT+XJIO/dfj65cu+l0uUwb+Gg5qDBoDpy76fS5SDtR+s/WCtg2uyDiw5WDI4QxlYOsGSwbWQwUXLQTt8gQULFixYsGDBggULFixYuERcNq3RggULFixYsGDBggULFizgNZ0zTdN+RdO0x5qmLWua9utndVFXCZYMDmDJwZIBYMkAsGRAWHKwZABYMgAsGQCWDAhLDpYMjgPHaV+oHXQy/M8B/BkctFT9XNO0f6jr+sOXvEbXtIvoOnsxMFBC3TiGDGw2m26zXa+EpSKHGRxjLWiadp1lkMNBAw1LBsfUCddJDqfRCddNLwKWTgAsnQBYOgE4vU6wZHC9ZACcTi9e47PhjdWLRKfTyei63tfrd6d2zgB8hINC1xUA0DTt7+NgcO3LNhtcLtdrfKS50Ol00G63+X3zODKw2WyIRCIXc4EXhFarhVKphOOuBZvNBp/Pd5GXeO7Y399Hs9nE/v7+qq7rr1wLNpsNXq/3Yi/ynLG/v49Wq3VsGQDXTw6Uga7r2N/fP5YMNE2D2+2+uIu8AHQ6HTSbzRPphEAgcJGXeO5ot9vY29s7kU7w+/0Xe5HnjP39/RPJALB0AmDJALi+dkK9Xj+2XtQ0DR6P5yIv8dyxv7+PdruNTqdzbL143c4GolQqrR31u9dxR0fQPYF7Ez0mbGua9lc1TftC07QvzNR8pEcXlVO9hyGqcaVkcFbodDqvlMN1l8Fx1oIlgwNcZzlYMjhAj3t6I/eDIeL7xsmg0+m8UgbA9ZaDmXTCy7rDdTod+ToL+8j4uZchg1fdp/HrrO+71/UY8FKdcOYXYBK86bbScfA6mbNjQdf13wTwm8ABpe+8P+9VUB90p3MwP1FdKOeRQlZl4HA4Ll0GlwFVBna73ZLBGyoDwJIDYD69eBmw1oElA8KSw/nKgHYPbZ6TGLyapsFms52LbWTEWcmA96f+e9x75n0a/70oWGeDpQ+A13POtgCMKf8fPfyZqaFpmmzS192E6nsd4krIwIheSusksrDZbKaQw0kPnLPEdVkLrwNLBpYMiB77642Tg6ZpYgwf4o2Tgc1mQ6vVUn/0xsnATDqB+1J1tvi9+gV0Z52YVXqdz70oGbyuPWP82x5Zv1PD0osHMMN+eNV6vuxav9dxzj4HMKdp2hQOBPuXcDCx29Sg8rHb7bDb7QAOokn7+/snTmUbHD0XTCwD430ZU/fq96qSPk7EzOFw8HUXvhbU61ejgsb7NWZH1fs6i01os9n4+S6zr4Xzgtll0CtrbtwDhLo2TrJWFBmYXiecJxR5XanzoRdepkteBrvdfmn74VXrW0Uvw1z993WMcrPrhIvAZegE4xlIm8ftdsPhcMDv9yMSicDpdCIajSIej8PhcCAYDMLr9aLVaqFWq6HZbGJjYwMPHz5ErVYTW+mkOG8ZqPfK72nj2Ww2OBwO2O32rjV+1Hsc1kOpNdTfet1pbQZSfK+DXjwtlGD+pemEXnRX4MV5b8RlOGqnds50XW9rmvZvAfjHAOwA/rau69+c2ZWdE7jJHA4HPB4PNE3D3t7et/jIx3kYfJ/DpiBLuAIyUBeikWfN+6YxarPZei5UI/i3nU7nwteCurlUB1uNWPdyzACcKVVDO2zq0Gg05nFF1sJZQzts+LO3t2daGRgjwcbvCdVYdTgcx14rigwAk8rgIqDI6kqdD0YYD27+7LhrwePxoF6vX9h+MAarjhO0UjMnNGDVe3wdB+0yZGA2XIZO4HkIvHAI7HY7AoEAXC4X+vv7MTExAb/fj6mpKUxNTcHr9WJwcBDRaBT1eh3pdBr1eh0//elPsbm5iVarJc7KSbNJ5y0Do3PGtexyuWC32+H1eqUZHYPJvd6Djb14r61WS+wiYwDjNODrdV2/0nrxdaBpGpxOJ5rN5qXqBFU/cq84HI6eTvhZZk+Pi9eqOdN1/XcA/M4ZXcu5gRvXZrPB6XRKV6xwOAybzYZSqYRSqSRKh685zsOw2+3Y399Hp9OZOdeb6IHTHJiqguFBrP7cqOROwtPWdX3+xBd0ChxlRDudzq6fv+x1vSLKr7v5DpX+A13XP3ytNzpDXDTV0+FwYG9vz7QyUA2V42TDTppNBw5kcGjEXLhOMBMOaX0XohPOC0Zj7GVZ1l5wOp2o1+sXth94LXa7XWU0fOv3QO9sQy+GxevCjHpRxevS4I6Di9YJDKzabDa43W64XC54vV709fXB6/UikUigv78fPp8P4XAYPp8PHo8HHo9HOscGAgE4HA74fD64XC44nc4uG+mkOE8Z8Hk5nU64XC7p9Oj1emG32xEMBiUY3yvQpgbq6vU6ms0mGo0GzzO0223pNnkW13rZevE493Gezsgha+3cdUKv+6RtrzrpRirvq97rIhy1c28IYgZQ4C6XC7FYDB6PB5OTk7hz5w6cTifu37+PBw8eoNFooFwuC0f+oopgT4JekVGg24jgv+rvGf232WxdyisUCsHpdMLn88Hn86HT6aBSqaDZbKJarWJ3dxe1Wu3IdO9FQ82UMRpnt9vh8XgQiUTkIPL7/UKloBza7TZ0XZf7qtfr0taV72eGezwLGFP2htoXgUpfpQzMtuZPA3Wf8Hun0wmn09lF33E4HAgEAnA6nUKBAYBmsylfmUwG5XJZbQEM4PI56RbOBi87wFUnp5cuOS71+7yu2Wgw2Gy2LmPU6XTK+nY4HHC73TQOsb+/j/39fdH3zWYTxWIRzWZTxsRcRsT4ImBkj6gwZkmu2v1rmgav1ysO1+TkJPr6+hCLxXDr1i1Eo1F4vV4EAoGuM9LhcEirdwCIRqPQdR0DAwOIxWJot9soFApoNBrShfOyZaPrutAWbTYb+vv7MTMzA5/Ph0Qigb6+PrjdbgwMDCAcDgN4YdcZA7Tc16VSCfV6HcViEcvLyygWi0in01hbWxNHTaV2XrYMXoXTOpS9suxXBcbAE/c417rb7UYoFILb7Ybf70c0GoWmachkMkin03LWMxhBvU8b+qwC+i/DtXfOVAOVaf1AIICJiQl88MEH8Hg8qFQqWFlZgaZpqFarx6Y1XhZUg1ulIvIgNR48Kk2FqX6Hw4FQKISBgQFxbCKRCDqdDnK5HMrlMgqFAnK53ImziecN3jsNETUq6Pf7hU/PrOZhZlMoGYVCAZVKRSJhqrN3XaAqJN5/r8hPL6f+usnBmDl3uVyIRCIIhULweDzo6+uDx+PpMsQrlQrq9XpXFFXTtGttsL6J6JU9Ur+n4QdADDLupf39/TOvXT0JVN0FvMgG2+12yYh4PB4JSHo8HglaURe2Wi3k83lUKhVUKhUxvDm/k59zHde7eo4adSPr0fn/V73Pcf/2IkAqaTAYRCAQwPT0NCYnJzE4OIiPP/4Y/f39ovc7nQ52d3exu7sL4GBt7+3tSYDT5XIhHA4jGAyiUqmgWq2eSfboLKBeB+vL4vE4ZmdnEQ6HMTo6itHRUfh8PoyNjSEej3cFYNUzkt+3223ZD5lMBh6PB5lMBg6HA7u7u2JTcX+Y4Xm/DMfJAr0Mxky72e+3F9Q9TlvI7XZLQiIej2NkZAR2ux26rqNUKgmd1eicAd9mU5wXrr1zRui6DqfTiVgshmg0inA4LNEWRkxUw8sMi7CXsUAlxAPY7XarA5C7HDKXy9XFuWYk1ePxwOVyIRQKIZFIyEJlITA516SAmmFTqrJwu91iZPf19QktY3R0FIFAQAwTOmd8rjRI8vm8KOFyuYxsNitRkqNq1a4K1CyRmiH1eDxdWUG1zmpvb0+KvY374CpBdbL55fF4JEMWDocRiUTgcrmQSCTEOaPxqjpn9XodtVoN1WoV+/v7cDgcqNVqUodxHrI5ilpmxEU9l6v2/F+Fo+hrRzlXNpsNkUhEoqpEu92W4A4d99O0KD8tSMmhjvd4PFJD7fP54HQ6kUgkJHsWi8WEURAIBCTIQKpWsVhEtVpFpVKB1+tFuVxGuVxGOp02dlq80lCNK0bAvV4vvF5vl5PdbrdRr9dFRtSHKo7aGxcRUT/qM+12u7AD+vv75TwcHh5GX18fAoEAarUacrmcrNt2u42dnR2kUilomiZnqd/vFweVWVWVNWAmuN1uJBIJ+Hw+DA8PY3BwEIFAAB6PB81mEwCQyWRkn/KcZ+CWGREjtc1utwsNlBnmcrmMfD6PXC7XlV0xm658mR56mc7r9XpjhtH4GjPde6996vV64Xa7hTXmcDjg9XoxNDSEQCCAWCyGkZER0atutxutVkt0QKvVOlLfnyfeCOeMi8vv9+PGjRsYGxtDX18ffD6fGO2NRkNobpe92Iypdnr8NDb9fj8CgQBu3ryJ0dFRVCoV7OzsoFqtyuIBgNHRUQwPD8Pj8YgxqtJc/H4/QqEQ7HY79vb20Gw2UavVRCm3Wi1xzi4rWmakcdJgSiQSCIfDeP/99zE9PY1AIIChoSH4/f4uyoaqgKlIs9ksvvnmG2SzWayuruLrr78WxVsulwFcXHTkrGCMAodCIQwODsLn82F+fh7z8/Ndxa7tdhu1Wg3tdhtbW1t49OgRarUa8vk88vn8pRgZpwWvVT1o6Zj29/djdnZWCt8nJibgdrsRjUYlMkyHjcYNAOzt7WFvbw/lchlffvklVlZWkEql8MUXX0jGlZmUs5CRmgk/ykkz1g69qq7ydSm6V5XWpaIXvVVd2wxAMQuvNsdwOp2Yn5/HjRs3hPrtcrlE39brdSSTSTx79gyNRkOMNb73eYF1QA6HA4ODg0Lj6u/vR39/P1wul6xvt9uNcDgsgTqv1yvOGTNnzJgVCgU8ffoU+XweT58+xc9+9jNUKpWeDXOuGtSz1OFwSA3VwsICFhYW4HA4JFhbKBTw+PFjFAqFruAdwczay/bXRelPNfvj9XoRjUbh8/nw4Ycf4uOPP4bf75ezv9lsIplMYm1tDel0GhsbG6jVatja2sLW1hacTifGx8fR39+PgYEBfPjhh+jr65OzkfRXs+gDXkckEsGHH36IwcFBjI+P46233oLH40E2m0U2m8X+/j7W19fRarWwt7eHfD6Pvb09sft0XZfMstfrxezsLAYHB+FyufDee+/B6XRiZ2cHCwsLKJfLWFxcxOLiIur1OvL5vGQTzVIScdQ+VfW5WnPN36mvV+2JXk65MQNphjXR675tNhsGBgYwOjoKr9eLkZERxGIx+Hw+sY9isRgGBwdht9uxvb2NZDIpzlmz2USpVMLz589RLpeRTCaxvLws6+Y8s+ZvhHMGvDiIo9Eo+vr6xCmhgWWMkJlpsRmpmYyQDgwMYHx8HOVyGe12Gy6XS4pYNU1DIpHA2NgYPB4PhoeHxQhlxyaPxyNOTKlUQqVSQalUkowcDyszQdM0MTii0SjGxsYwNzcHv9+PoaEhyRL5fL4uWiOdcF3XEQ6HUalUJJIYCAQkWgq8UDZmUTonAdcMa6mCwSBGRkYwPz8Pp9Mp90OjjPUlyWQSmqahUql0UWOvGlQ6mt1uh9/vl/0+Pj6Oubk5WT80thm0cDgcIiMGK4rFYld9IiNw6mF1luvkZRmY43zOUTS900CNJF+1fdALR2XOaGjTUaMB43a7EYlEMDQ0JLqSmSW73S7UV+rJi8oqqKyIcDiMwcFB2ecjIyNy3X6/H06nU9a3y+WS9atmhWq1Gvb29pDL5bC3t4dgMIhcLicBjqvumBFqZp2BmGg0itHRUXHW7HY7dnd3sbOz0xXoPM39X9SeUTNnbrdbKP6Tk5Pwer2SQWWGtFgsIplMYmVlBZVKBZubm9jc3ITL5RLH3WazoV6vS7aAbBo1IHVZOsFIL3O73ejr68Pw8DCGhoakzqxcLkuQjfder9eRSqVQr9ext7eHer0ugXufz4dgMIhEIiFsikQiIcFsXddRLpexs7MDv98PACiVSi/NKpkFdKZ6OWdGJ60X1Rf4Nn3XrPersmfY/C8WiyEQCGB8fBwDAwPw+/0YGBiAz+dDNBpFf3+/BL0CgYA47s1mE7lcTujf1Wq1K8h9nnrx2jpnXGB2ux3RaBQejwdjY2OYmprCzMwM9vf3sbe3JxkzNW1/2QuOm4iGIg1Nu92OgYEB9Pf3IxKJYG5uDvPz86jX6+jr60O9Xkej0UCtVoOu64jFYojFYhIVVOsMNE2TBcg6rHw+j3q9jrW1Nezu7qJYLGJvb+/SCn9Vo4DNS1wuF27evIkbN24gFAphYmICsVgMdrtdjAw1QkSah6qQNE0Tqk+73ZboyLNnz7C8vNxF7TQTzfUoUEZ0vJ1OJ6ampvDuu+8iFAqJwWa327vWONP9IyMjkjV99uwZgAPnrVqtsvWx6Zx0oxNis9kQDAbhcrng8/kwMjICn8+H8fFx3Lp1C8FgUA5u0hrokFFXNBoNiYKqBxRrNDweD1ZXV+FwOFAsFpHNZiW7/DpQ1ymNZ3XNGbPHvbJsBLNAakOH04KBDfUargJ6XauaGaGT4nQ6pVGAy+WS9UPZOxwOTE1NYXJyEk6nU9pxNxoNhEIhNBoNMXqpP9Ua3bPSGdRBvP5IJIL5+XlEIhGMjo7ixo0bUmdLOremaWg0GqLzVTq8zWaT4INqYGmahsHBQUQiETHa8/k8dnd3kc1mr2wGTZUfAPT392NqagqBQAB37tzB7du3uzKoiUQCdrsd+Xwe29vbePr0qVAAyd4Ih8NddEgaq9ybhUIBxWJRzpfzoofSPiAzhPpreXlZgk+kcjFLRobE9vY2Go2GZJJonJKVwnOTupE0LzME79Qznk1MKpUKNjY2JMC4vr6OjY2Nrplte3t7KJVK0tiDwUl2qWRmeWtrC6FQCBsbG6IXuP8nJydRLpdRLBbRarW62DaXCaPesdlsiEajiEQiEqwlxT8YDMqaN2bRdF2XzGK1WkUqlUKtVhMZ81ww1qWqn32R6MUwCYfD0o301q1buHXrlmTLmKRQGQX0BejA03fwer2IxWKYnZ2V4P3W1pbsJzr352EfXVvnjJxiHsCDg4OYn5/HO++8g/n5eWxubgp1oVKpyGGldmK7jIWmGgY0qFnc63K5MDExIQ7J3bt3cePGDSnipcNZq9W6nM9Op4NarSaOGA8apmir1ap0qWm32xIxYC0baY6XufFCoRDm5uYQCoXwC7/wC/j444/h9XoRCoXg8/lQr9clw9FoNFCpVKDruhgspHOyne7o6Cg6nQ4GBgYwNTWFarWKP/qjP0K1WkW1WkU2m+1y9MwK1VhiACIcDuO9997Dn/kzfwbhcFgin6pSZUc3u92OSCSCqakpNJtNfPbZZ9jf30e5XMb29rbUV5kpSmakqJGqFI/HEY1Guyg5w8PDmJubE0eUbaKB7vXMDEKpVJIsNI2TiYkJ6XiWSqUQiUSwtraGQqEgzqvaQOCkULPBpKXx0AReZNO4z1Uqtmowq5z6ZrMpke+TXIcqYzZEMRrlZlkHR0HNkPCLgQuHwyEGSzAYxJ07dzA+Pi61OYyI8/VqTRIdGxrCnU4H8Xhcmgc8efIEuVxODvazvhdmvwYHB/Hd734XY2NjQmv0eDzyt+w4x+ZWRodK13XRc3a7XaLKbrcbU1NTYrA0Gg3kcjksLi6KIXqV1gHwbfnZ7XZMTk7il3/5lxGPx3Hjxg2hfPNva7UaZmdnUa1W8eTJE4TDYZTLZTQaDezt7cHtdmNsbAyJREIaxtCJYfDzyZMnePz4sdBFuW/PA1zjrIWs1Wr48ssvsby83EXt5+9ZR7O3tyfX3G63JeM6OjoqAW3gRY1luVw2FdWZ19HpdFAsFuHxeLC+vo58Pi8Uzp2dna6GWGrQSu3UqVKbnzx5IoyK6elphMNh3LlzB7/0S7+EQCAAAIjFYkin08hkMtja2uraZ2ZwUFg/PTo6ioWFBanJj8ViXSwjUnzVzBr1R61Ww+7uLr788kuxD9fX18VRZyDATHWIlH9/fz/u3r2LWCyGDz74AO+//z7cbrf0XVCfkdqlvVgsyjpnACYSiWB8fFz2+srKClwuF9LptMjgPPbEtXXOCB6wwWAQoVBIGkZwfgWji+pGvazNpaZiOZ/E4XBIjZnb7RaqWjAYlM6E5EwzwsUMia7rYsCpxlqtVpM0LRdYOp2WFqLGWSZmUDaMbJPOGI/HxYCmQVKr1aSjVLFY7FKWVECMfDJTQc6xx+NBOBwWeRaLxQu/55NCNZLooFBGlFMoFJJ0PB3vRqMhWSMqHDaD4f5Q6w3NBnU98kBhS1w2/eAcn3g8jkgkIg0/1GyUMUPWbDYlMthut+U1rEVg9zPW8hidmdPuEwZjPB6PtPdlloPgNdKg2tvbk7lD6ufz8GFQ5TQRe0ajVSdEpTRdBaiZc5vNJlRn7hE6Z/F4HP39/QgEAhgcHOxyzoAXGTeg92gVdrJrNBqiY/j6s9abanfaUCiEWCwmXUcZ/aUTsLe3J4YDn6G63mnEq1Re3i/rmtm9l82EzGCQnxaqzg8EAojH41KLRQo8wa5/Pp8PuVwO0WgUTqdTdCdrVtn9j047HbNWqyUZCuNQ2/OAuuZISaQTpjZ5YsBWbYhhdLTZMIwZfL4v140Z14Fq35RKJWn+kclkkM/nv0VBNwb3gO5gPDvzkuLL2mPqaZ/Ph0gkIl0t1czpZYKfz2fEzE8kEukaLeD3+6Xe6ijnzOv1olqtAgASiYScAblcTgIZPBPU0QKXHcTlZ1NH0h6IRqNSV869znXALCH9AfZdYPMknhnqiBI2YTrP537tnDN1gVLAk5OT0jyDhsvu7i4ePnyIdDotBYBUPpdxzaoT4XQ6EQ6HMT09LYcvuyp6vV7hRzebTXGomKZngXq1WkWhUJBMALODqoIul8vIZDJCaaNDd5kRIMoDeNES2GazYWhoCO+99x7i8TjGxsbg9XrRbrexvLwsSvjJkycoFouo1+sol8vodDqIxWKIx+NSe8FOVMPDw9JKNRQKwev1Yn5+Xigfn3/+eRfFzWzRYuOaYWep9957D6Ojo+jv78fe3h6y2Szu37+P+/fvS1Z1b29PGgj4fD709fVhbGxMahYSiYQUQZsJ6t5mFJw0X7VBDvdOKBSSuhvjc1OjzaTrrK6u4ssvv0S9Xkc8HhcDjRRat9stzq/P5zszQ4XZnNnZWdy9exehUKjLsCMFr91uy94lBYMGFjNnpGuqM6uOkmMv2iQP6f39fWxtbWFzcxONRgNbW1vI5XKvfa/nBXU/0LkNBoMYHx+XwBxrUaLRqMx6Gh8fR19fn9AW2U5Zrf1TqVxqAA04yOjPzMwgGo1iZ2fnXPSD6ijSeOY9MWPWarWwubmJ9fV11Go1rK6uYnd3t2t9qHIiw4AZ52Aw2JVxbjabiMViXTVrKpXWLHrwZVCNJr/fj9nZWUQiEbz11luYmZlBOBxGq9XC0tKS7EEG5xiI0TQNgUBAzleuBZ7TDHjxXxqsaufkizLa1WwQrwnAtzJHRK9gcCAQkHpFBjMv8h6OA6PeqlQqWF5extbWlrBmGKxV/57fG7PIxr9RjXYGr1OpFLLZLHT9oBEUM+rU0WbYD+r5SAdicHAQCwsLogtJ2aWOJKtK3dNcJ+xu/sEHH6BeryOXy2Frawv1el0oo41GA9vb211nw0XqB+O65Hr2er0SgGHgSu3IrcqATKlqtYoHDx7gwYMHcLvdeOeddzAxMSHBPAbyJycnpS6NXTuN7IQzsQte+x1MBFVANG58Ph8mJibw9ttvIxKJwOl0yqZ7+PAhdnZ2kEwmu2pxLtpBU6+bh28sFsNbb72FwcFBxONxDA8Pw+l0olariWfP+1BbOy8tLeGnP/0pCoUCstks8vm8cLLVLnNqbZ0ZFIsRuq7LBnI6nRgeHsadO3cwMDAglItyuYxnz57hm2++QSaTwb1795DNZrtqh9S6u0gkAq/Xi0QigXfeeQfxeByjo6MYGhqC0+mUZhG7u7vY2NjAxsbGt7jVZoGRyuZ2u9Hf3493330X8/PzEgmq1+v4+uuv8bu/+7tdc7tI2QgGg1hYWBDZuN1uxGIxOazNAtVJZrTP5XJJd65wOIwPP/wQMzMzXbV3xvcg1FoNymV1dRV/9Ed/hFKphPHxcaEPj4yMCM04EAiIM2+k+572UCK1bGZmBt///vfFKKYDyqACqYx0KJnlUwdkspmPWktylAx47arRRrm02208fPgQS0tLQv02u3OmZsU5v+6dd97BwMCAzLLhTEcj1bnXTBujXFTnjF/hcBgzMzNIJBJ4+PDhmZ8dqtPJrBabk6jOWbPZxObmJr766isUCgXcu3cP6+vr38oQ8HvW1TocDqE1Tk1NIRaLSbaFwQkGC+iY8D3MeG4QxvsOBAKYm5vDyMgIbt26JR0unz59iuXlZakvJK04EonA7XYjHo9jYWGhS37qbLhGoyFnDoOhdGhUFspFyIpnOwAJ4Bihnvlcq7SVOFKH60p1zvj+ZgL3BimXvXBSx0nd/3t7e8hkMnA6nUin08jlcpJ5eZlzdll7g8+HbA8Goefn56UhViKRkFoynicsg1FBWmQsFsPExAScTidKpZKMkllcXMQ333yDYrHY82y4DBmogQba0GRY8Rxk7SHLCJxOp9DSS6USvv76a/zwhz8UBhFprEwU0DkLBoPIZDJwu93fGjNxVvd9rZwz4MWiIP2DdJZgMCgNINQZRlykgDmcFNUoUIsUy+UyHA6H0PYY6fV4PF0dlWq1GorFomwadplj6lZVPlzIZuKRG7NUbJJAZcOCXXbS4rBstRsTZQFAImms42m1WnA4HMhmswCAYDCIWq0mdSX8v9qpjIeUGeTTC+p8G9ZVsckJO3Cy5pCKBIDU7DETS0qc2lLcbGDUmgYEsyCkLPMe1A56NJB6FTHTuGJkkOMUWPDNLJS6Jrlv+P/XBbMzjUYDxWJRqKb84v7tRTM0ZjN4fTTkezkLxsizCjIO2u22FJCTKm1WqDqfrAN2XYvH4zLbktlPFvazm6uxFpcZSjXLYLPZZH3RQFNl3Ss7e1ZQHWg+f/W584vtn6kHGRHvBUaRGZzQNE3axjN7BEB0hjEybGao10k9QKecgT11v+3u7kr5AACxE3gGqEPIufey2aw4Z7lcTt4rn8+LwdtqtV67Kc9ZyOBVUM8POv8M/tK2MPOzp25XdSBwOsdM/Z5nDbOkPBsByLPtVdN5Wecm9YRqQ9JGAoBcLie2Ezs1q0E+9YxT9w27FJOua7fbhS7Y6XSEFmmmNdJqtVAqlcQW4PXRPqSN5Ha7ZVwIAPEZWMrApjoApIEM6ePnnVG+Vs4ZF6SmafD7/YjH4xgYGMDExASmp6eFzthsNrG6uoqNjQ2k02k0m81LHbasfiYVYSaTweLiojQCoaFFOkowGMQnn3wiUXIamJubm1heXkaxWPzW4U1DTTXYjArMDAY5r5dDptlpjy1u2Skvk8lgaWkJX3/9Ner1unSXZEoeeLFJbTabGL6pVAr5fF6ocFQ2TP2T+kTH14yZM6C7aUEoFOpytJ4/f46f/vSnyOVyePz4MUqlkihZr9eL4eFhfO9738PExASGhoYwMzMDm80mmVazDRxVOfHMoobDYZnh5vV60d/f3xXxBSAHULvdlrk3pLw0m01Uq1VsbW2hUqlgdXUVz58/FyO9UqlgcHAQd+/eFWdIvY6zcmDb7TZ2d3fx1VdfoVKpyD3wkGSmTh0orH42o55Op1PGBtDJ8vl8ALodGBoYvSgeNNRarRYymQxSqRQAyDWZ6QBWjSLKJxgM4u7du5iZmUE8HsetW7dkL3Owrhr82tnZkfqr3d1dCW6lUqmumWV+vx8ffPABpqamhApMQ4UDno2Z2rOAur5YP1wqleSL10Z6Ui6Xk2BVpVIBgG8F4LgOaLQwc7y/v48f/OAHkkUPBoPodDqyJ0j9N8MZ8Sp0Oh0ZPE/2zJ07d6RT8+rqKvb29vDpp5/i008/haYddKpkbQq72amz4UqlkjQOKJVKIhPW89LQbbfbKBQKKBQK4sxddPDzVZ+l2kperxfhcBh9fX3SPI2dD9kYivdgRnA996IovgyqLjPOl9R1HV6vF5OTk9Idlc22kskkkskkdnd3u6jllw2uMTplnU4HDx8+xP7+vgQnfD6fdPXkOVgul791DwzQcmSTz+fD1NSUnLtkqOzu7uLp06dYWVnpCvxf5D0bner9/X1sbGzgBz/4gbAMGHBiwCEQCGBkZET6OcRiMbhcLkxPTwuV+e7du7h9+zaAg7mnqVQKGxsbWFpaQjablTPivHpVXEvnjJSsSCSCWCwmSqdQKGBzcxO5XA67u7vIZDJS4GhsJ3qRUD+Xm6RUKmF9fV0aXqjGZrPZRCKRwPz8PKampiQ6Uq1WkcvlkEqlUCgUJONkzIyZrS16L5DWGAqFpLsanY9cLodMJoPd3V1sbW1hbW2tqx5HbQ/LIm2VzulyuVAsFqWz3dTUFOr1ugzs7HQ6MgvOyMUGzOHAAt3RPb/fL4XcHJT7+PFj7O7uIplMSnTc5/MJtZMjCSKRCAYHB+W+mW0yi3PG+2TGKhqNYn5+HvF4HNPT05ienpbmLjSQ+bzVeSXsNsW5N+zOyAGT7EjFTDoDJSp90Jg5A15/Pezv70tgIZ1OSy0AAKmdIt1oampKGkAw28sDUc32qsOHVRib6LCjLWvz6Gy0221pmsHmI2YE1wafP9fG+++/j3A4LPKiTtD1F50K2+226BF1HaRSKayurqLRaIhjw8Yy0WgU+/v7iMViXZkzRprPWjcYA3ekqDMz5vV6JYvGJiDValWaV9BwVQMKaraP0WR2K242m9Jmur+/HzabDdlsVj7b7HRGoJsKyuYNiUQCk5OTmJmZkfb4hUIBjx49wr1796DrutTese04GyWx9oxNs+iI8bxhl2c1y6QajGZlIajnB4N7zIhwvtPu7q4E7MyM02bJ1H+N553L5ZLB7iwvCQaD2N7eFufbOIz4Mp8zP5sNUuikcBYjGQOqc9ZoNCRwq+4bnm+BQACjo6MIBoOw2Wz44IMPEAqFMDAwIGcF9etl2QtGSikAqSFjEIpNTMgcUs/SyclJYeIMDg7Knp+dncXExERXAJc2ZzabFQrzeTmj18o5Uw9KDpccGBiA1+sVClM6nRbnxYxRQCoJtXGHqjg4x4kRQbfbLVkj8uCBFwdCr5qDXvdstuwZoRoUVMCk7zBFrToSR1276qCSwsKDOxqNIhwOA4BQQpk9Mab7zSQbFbwupt1JrVFpu9wXLJRl18FWq4WdnR20Wi2k02mhQ6o1JpcFo+zZrKJWq8HtdiOfzyOVSkmNEYv02eCmUChgZ2dHmlpsbW3JAG5ShkmJ5UGr0sjUL03TpC6Djj0N/tdR0Pws1gmqwQC1gxYjoKRr0Jni6+mUs7bW7/d3OWfkzbNFNvVLMBiUwme32y33w4Ob+8HYTOCyoK4J0hMTiQSGhoYkGMdAjkpDYnaIGVQ2ztjZ2UG1WkUymZQAFzOrrG9k0M84F4jDmwuFwktphGd136Rf5vN5ZDIZaa3PQNbQ0BC8Xq90rGNNqtr0gdk2novcC9SNzKCyocJ5toE/T7AxACnPvDfaAfl8HuVyWaLfDGCxvs7hcMg5o2manK90hJlJVKmmBHXIZevPl4H6NBwOY2RkRLLBNO4Z6KURelVgzKSoMNo5/FKZA6TG9/f3ywiVaDQq50Yul8P29rbokMvs8n0UeCaQwmi320WfqQ3kOB7COM+ScnG5XN+izzMDVS6XhR5pNtoz96U6RoLUbwbyGYzSdV3KhJgtZDdGBnI583F7e1tYCSq18Tzu/Vo4ZxSMWmd248YN/Ok//acRiUTQ19cnnVW++uorrKysYHNzE81m0zQbSo168F/Savi93W6XIdqxWAxDQ0OIRCIoFotYWVnB2toaksmkKBtG0dUvoHsDqlRHbsjLPlSM3GceqrxWtv7PZDJdyoGv6ZXpUg1dv9+P8fFxibLfuHED4XAYhUIBz549QyqVQiaTETqcmZQOYayrYGajXC5LZjGVSiGdTqNWq4mhMjc3h4WFBYyNjWFoaAjhcBjb29tYXl5GuVzGvXv38OzZM2k8Y6b9oRpXpKORdhMMBvHOO+/A4/GgUqlgfX0d5XIZKysr+PLLLyUzlslkupQ1I9/qzxgF5DwgKni22KYz4PP5hNqkdkU8qcwY0aMDoQZVNE1DoVAQp+D58+ddw9TViCUL1dVaGbVpic1mQygUQjAYBACJ/k9PT+Mv/IW/IF2o3G43Op0OMpkMnj59Ks4HnXXW4FwG1EAVqYoejwfvvPMO7t69i3A4LF07GdntdDool8tIp9NoNBp4+vQpnj17hnK5jEePHmFzc1OiyTRo2EyBRj3wos082/MDQD6fx9LSkuy38zBiuZ64PguFgszo3N/fx8zMDLxeL6ampoS+zcY+6sB1RoSdTmdXFpBOR71eF33BDCJHr5hRB74K3KdTU1MYHh6WZ1koFLC4uIjd3V2sr6/LLD+uEyN1md+zTf7LzlSgu+mG+q9ZoFKcSeX65JNPEIlEEI/Hsb+/j0KhgG+++UbWmZlspZdB1Q+9shpqJpk6lbM+OVJnZGQEsVgMAwMD+M53voP+/n5hM+3t7eH+/fv44z/+Y9kj/LzLzJCqNg/PNe5t0p9V20jt3qlmRSkv9fd05NSf5XI5rK2tYXd3t4sWaRY9obb3J+ic8eccl7S/vy8MrUAgIE1AgIMM3PPnz/EHf/AH0qFyZWVFAlfEeQTur4VzRtBwYLelkZERqddgtCCbzWJnZ0dqssyYDeEiVw1IRvJcLpfQ/NTIODvpMHNmLOA/6iBRf2YGx8wI9aAkGD2m88TNZlRAxvfhv06nU9qj8ovzwNhxi+9t9tbRxpQ+aXzqF2ulaJSxfo80wGazKXV47FpEp8NM980DlZkxrmk6RzSy2AyFka7l5WWhoeTz+SOjqUbar7HxAjMPzFjT+TmLmkRjBO6oQ67X/D1jMIN731hrREpoJBKBruvSJIbvu7e3J93/AEiTg1KpJNSty6zNVcHnQXpnPB6X7F8ikRAHlGg2m0L9TqfT2NzcRLFYxPPnz7GxsdGzwQYDGgC+ZcwZM2e5XE46xJ7nPTPDWigUYLfbUalUoOu60F9Jy+nv75cMKx0yzgF0Op1Sv0snvFwuo1QqodVqwW63S20zI8xmMbqOC/X5qbNNbTYbms0mstmsnJfqTFB1ttfLoJ4zvVgnl70/joKRfseM6/DwMILBYJetxCAf9ar6OrPjqOyZUR+S9sa94fP5MDg4KLV34+Pj6O/vR6FQQCqVksxZMplErVZDvV43hT4kjEF+lVXUa033um6VNcLssPH31KflclnOBjOBdjO/V21qBmQZENV1XZI64XBYSls4H7ZSqSCZTIozyvNQDYqfx/O/8s6ZatCwAJ7zCIaGhoT2V6lUJBVNPqqZNpUKNXIOvKAyut1ujI6OYm5uDoFAADabTQyDXC6HfD4vHXXYpdLv90v0lE1F2FmMtWpsfsDZYOdZ5HiS+6fhrXa9YgpaHcxN2pGRNw28aBOstozmPDDO9+p0OrIBHz9+3JUpMKtRomY5mXHZ29sTg7Wvrw+3b99GoVAQAyscDuPWrVuYn5+XrE86ncba2hoWFxfl0OG9q7TGy4wI8pmycQ8bH3CfT0xMyDqv1WrIZrN49OgRdnZ2sL6+jnw+L93T6FwYP+Nlnw+8CF4wQEJD76yipcb3MDrdx3mtMWOs7mFNOygUZwE7f88IK/+OBzIdXrUr2WWvBSNoSAWDQQwPD2NgYEAyadwXPEjX19fx8OFDqS17/vy5dDY0goabw+HAwMAAZmdnZZwJnRu1AQSdftJizwt8ju12G7lcTpq2sJMvGxu5XC4sLCwgHA6L4clOfKR6+v1+aSQSCASEBh2Px6U5Sl9fH2q1moybUelQvJ6rAO551uRx1iONSjUoc9zgw8v2q5nlwvXp9/sRjUbh8/kwMjKCwcFBOJ1OoeytrKx06U3A3PdFqDYPa7BVe4o/px5nJjwejwutc3h4GLFYTPYHaW0PHz6UgE65XJa98LKgsBlwEr2tBiHtdjv6+/tx69YtDA0NYWhoCK1Wq6suN5fLoVQqfcvpNwPUfc2zmxTnSCSCmzdvdtWvqyMkVDCYVavV4Pf7EQwGu4LgbMDCfaJmbF9HHlfeOQPQZbRzSOzw8DBmZmbQ6XTw4MEDrK+vY2VlBRsbG0ilUl0F/mZaUL3ATjuBQAALCwu4e/eutNVPJpPY3t7G9vY2ksmkRIBcLheGhoYkkjoyMiLDmBlJTaVSeP78OSqVCp4/f47V1VXpYMcD+DJT9WomiOloXddFHq1WS4wx8qdVSgOzRXTk3n77bYyPj2NwcBAfffQRBgYGxBnlIMtPP/0UxWIROzs7QuUwa0E38KIdPOtGyJkfGxvD9773PanL29vbQzgcxscff4yFhQUxunZ3d3H//n388Ic/RDqd7qKFmCV4oToNwAGVzOFwIBgM4s6dO7h9+7YcJsViEVtbW/j000+xvLyMarWKfD7fZVAeZUipgZ6jsq/k4KsOmtGBe537PAt5vyyKSfolnXiHwyFyJUWSDVTU7KvZIueadjBjbHZ2FtFoFHNzc5iZmZHAEzODm5ubKJVKuH//Pv7wD/8QhUJBAlnqHCoVdMzcbjfm5ubwySefyCHe39/f1ZQjn89jZ2dHKNa8tvOQEylZe3t72N7ehsPhwNjYGDKZjNRK9PX1odPpIJFIoNFoCJ2RASpVPuo6YVSZxjiNr3K5jC+++EJmJjLqfBXOTeDFet3f30elUoHX6xUHjXXFdMiMjX7MaHC+DtQMQigUwvz8PCKRCBYWFjAzM4NGo4Gf/OQnWFpawvb2tlDWjIwCs0LTNGnw4PP5kEgkkEgkupomqY2S/H6/1OeSQUPWld/vR71ex+7uLnZ3d7G0tITf//3fRzabxebmJrLZbFcdqtlwmmABZcQsvN/vx9TUFH7xF38R09PT0n6enY3v3buHUqmEXC7XRa03036hDUg78NatWxgbG0NfXx/effdd9PX1IRaLYXh4WM5Do03ApIjT6exqtkKbm/PwmNxQa9R5DafBK50zTdP+NoB/CsCurutvHf4sBuC/BTAJ4DmAf0HX9fypruAMwRoDtWsdZxNUKhUp4jPWKB0HahthUgnpPBwaZf8jzkkO7ERGJyMUCkHTDubSkH5idEY8Hg9CoVDXZHM6Z2wT3G63USqV4HK5JBPBw1+toSFY32Wz2RCJRABAsk6H1LkzkwENCGZw2GWHFD0+Z3YfZJaLRjgNZc6qCAQCiEajQudh63ymravVqjQ/MNaxqVBpgmqbcsWAndM0LXpR+4HpezqvlBfbSNNx5dwqOvksiiWdie231c6lRzmlasTQKINDmZ2LDFS6BbM+6sBpOqJs8sGsIWVzkoOUhvBREfLf/u3fxoMHDwAAoVBIrk/t3nWR66AXeq1fyo4OiHqAqIas+mXMRqtQRwxwaPl56UWVIQFAWp6Hw2HReXQ2uSc4ry6fzwvDgPsd6A4+qVQvzgyk3ohEIjL/iRmkvb091Go1/OQnP+mq9eW1KrVaZ7If1Cwyi/i5vklBpvPlcDhEnzNwwd8DkNfwfvk7h8MhgZBSqSTGLPc5HVAj6LhpmibzwnRdVylxF6oXjaDc1MYdr6JindSoepVeNINO4D53u91C7Wd5BOsLWaekdj8+iQPyMjkc2innaivRsWCXTmNHWw4YZ/ZQDeozeM01rGmaBCRoQzLjqgb3jWvlVXbCecpAxXHXMPU7z33KJxAIIBwOIxaLiRxqtZrMT1VpwcbPohOjaZo0oSID5jztBCPIFmMXWs6/5AxMdmg9qpyBtdzBYLBLh5B5oZbYnCXb6jiZs78D4D8D8F8qP/t1AD/Qdf1vapr264f//w/P5IpeA6FQCNPT01LYx/k1KysrWFxclCYCXDAn8WxZa6DWl1BxHT7EH+CM5MCHS8NhaGgIH330EWKxmEx7Z/3c8+fPUa1WcfPmTYyNjWFgYABTU1MyRI8daNRBmtVqVfj4t27dwv7+PgYGBjA5OYlSqYTFxUWsra2JgcNDjA6RekDX63WJKOzv75+JDPg8SNvhzJVnz54hHo9L+3vSbsbHx4WKoxZ322w2DA8PY2hoCIFAANPT09IemjVIW1tbePLkCcrlMpaWlrCzsyOGVy/nhFEYTpwHDpQQN3ClUimfhQyOIyMqEx4EqVRKqFsqjYnOKdvesovhgwcPsLW1hdXVVWl4cZyDmFF4VQasVXG5XKhWq+ciAzpMzGJEIhGMjo5KFvXRo0dYW1uT+YXqcMnjZj9Vpz4ajSIWi8k8KwCSQXj77bdRKBSwtrYm+5UydLlcNMwvVS8a75fXyQOHTsfAwABCoZDMwCkWi9IBlnVaRp1JUC/S0AfQlY04K72o6iDSsRYWFvDuu+9K8T4/m1mRra0tfPHFF9jc3MTGxgYymYx0YVTprcya8IDmaIlgMIj5+XlMT0/D4/FIfVY2m8Xi4iKy2SyePXsmzITNzU2R+97enmTfSqXSme4HtZ6CxlK9XkcwGBRHizpwb2+vK3NMsF4EePEMPR4P4vG4nB2Tk5NCBex0OigUCrh3715PKiiNH9YwAgd7haUG5YMPu9D9oNaYaJomwalIJIJwOCxGlZoxfx28Si9elk6gkckgr8vlwszMDD7++GPEYjHEYjFks1nk83msrKzg0aNHUnd8Gsrey+RweD1nZisRPA9jsZjQ1cgc0jRN6Nks8VAD+QxmMLgHQJpB0XEdHR3FjRs3kM/nZTAzu+v2on2+zE44NPDPXAangbr2qQMDgQDm5uYwMTEhNXfBYBDJZBLPnj0TaieHr7NcoNfZ4HA4ugL9PBsOg6nnphPoRNlsNgQCAfT39yMWi+Gdd97BW2+91TWrUtd1abpFfaF2ZNU0DQMDAzJOgOt4YmJCKNLJZFIYGY8ePZISIzVxctJgzyudM13Xf6Rp2qThx78K4J84/P63APwhLnGRcXEFg0FMTU0JX5j1ACsrK/jmm28kK8KHdhKB9apV6XQ64pjgjOSgZvRoKA8NDeG73/2uzNwIhUISDV5ZWZHulF6vFxMTE3jrrbfg9XqFksQCaEY5eIDSIXO5XJiYmJAC4FKphEwmI5EztaGEsQMOa5kOldBry0B1mOmc2Ww2rK+v49mzZyiVSpibm8PU1BTa7TaCwSDefvtt6VjJmR2tVgsOh0O6MdKoCwaDyGazePDgAXK5HL755hv8+Mc/lsYHKn+aToq6RtQZVES73ZbGAQCyAP6Z15HBccFIOrOn6XQaDx8+RDqdxtjYGG7evIlAIIDBwUGMjIyI0Ubn7JtvvsGTJ09kds9xnbOjZMDoGM5JBryuSCSCmZkZ9PX1YWRkBOFwGMViEWtra/j000+l2J9zTtROny+DelDRUOVAXhq9bDwyNDTUFY0HIK3uFZy5DF4H6jpmc5DBwUE5fH0+X1f2uFAoSIDrKOes13ybTqcj9E+coV5k5HtwcBCRSARzc3O4c+eO1JLys9ltMJlM4vPPP8fTp08lm6reh5qJYjaBdWY3b96UIedTU1PQtINW6gyI/ehHP8La2lpXAyGV+smB74c4s/3Aa2ZGk+uRz4jZUN4nh6fSwaI+53w3ADKOIBaL4a233pKARF9fn3wWB85ubW1hfX39W2chM24q2u22ZA3OUgYnAR20Xs5ZvV5HNpuVPfy6Dtox9CJwCTKgo65mQaenp/Gd73wH8XgclUoF2WwWu7u7WF1dxePHj8VOOI1R+TI5nKVOALozPpqmiXPW398vgVlNOxi0TkPbmC3nGcpOuTTM7Xa71PaPjIzgxo0b0iGVrKV8Pt+zU+/L7ITDQMml282EykggE2V2dhbvvPMO+vr60NfXJzTPx48fI5PJYH19Hblcrqv8xYheMjgvvdjrnugUseHNwMAA7ty5g48++kiuhTR+9l7gGCJeP58ng38MuDEb6HK50Gw2heq6traGUqnUFSSkvwGcLBt/2pqzAV3Xk4ff7wAYOOoPNU37qwD+6ik/50ioTgxbG3P+F3nEjGyw2JlzCc7yGhRhHykHVQan4SczC+L1eoWiQtoKjapwOCy/V9PwLFzNZrNdrUMZSYvH43C5XHJ4NRoNBINBSfWqRe69FhajzieRwUkWKK+1Xq+jUCjA4XCIYmQGwOfzyWBYGoXMYpC2wSiZOhuKaWjSONU2q2pm6lUwyKAFYOIsZXDEe3UZE6SSsTlAIBCQgySRSIiiUelY6sBa9ZA7DY4rg8NrP5EceG08VNlVKRgMCnWZTVGq1apQrPja04AHFQ1XBmfUaKmaeVGvU8GF68XjglFSUsDV+jka8Wrr5JPI8TR68WXrQH0/NrgIBoMIBAJyQPL5sCU+Ay3MKpGKZNRlXPMqHZz0Z9Il6YAym8hADvVQrwZK6oGMc9QJKs1GpaHSqWbxPoOS/JtsNotcLicypSGVzWYBQM4TRrnZQIT7QaUHvuzazksnvArUj2qzJLJA2BiAQ98ZdDyLzzXCIAPgmDrhrK5F1es8F2gzuN1usZNIa2dm+axrTM/bTiDYqIfUO2aHOe/TmAlkcEYdHcNxJFz3XB9sEBKLxdDf3y8BGe5BtfTlMmVwEqj60OPxCI2RQQzSOylXdQ6skdp53M9T/vZCbCU+I3b5ZqCKjKN6vS7OGXUFADn3Vbon6e42mw3BYFDOS9UmYfM9zix9WVnAy/DaDUF0Xdc1TTvy5NZ1/TcB/CYA2Gy2MyFjqhFrHqYejwczMzN45513ZBZRoVCQiFAqlRIu9EmzZse8piPloMrA4XAcWwZcEKyt4IN3uVwIhUJ47733MDo62mW4VioVfP3112g2m9jY2MDm5qZQW9ihkhuSqfpQKISFhQXMz8/D7/fjxo0baLfbyGQyUjzdK3N4WhnY7fZjyUB9zjs7O/jss88QCoVQLpfRaDSEF036ZiKRkNdwM9CIo7FNuVJmpDfQaePnvubaODMZvAzGa6xUKnj06BGcTic2NzeRTCZluPbY2Bh0XUcymUQmk8HKygq2trZkQPPrOmc9cCyd8Co58FmSNup0OjE8PIzbt2+jr68PPp9PZt5tbW1hc3NTDOaTdlNUo22kNUajUTFKdV2XodfZbFYCBKphZ4jIXqhePAnY2Y90DwZ2SO1gjac6iP2UzsOpdYJqOPCgZBOQ4eFhTExMiAPFw7ZSqeDp06d4/vw51tfXu4bocs+rz5kBm1gshlu3biEajWJmZgbvv/8+gsEgQqGQ1OI8efIEy8vL2NnZkeHVzNIfw3k9U53AZ8G6MwYgqRsfPHiAjY0N5HI5LC0tSQMU3rs6y42OeTQaxdbWlsjgww8/FON0eHhYgl2kv9K4Ua/nNDJ4HTm8DKR6dzod2O12JJNJOSu+853voFwuw263I5fLCV1ZbY5z1jiJTnhdGajnPPcPu/UODQ1hbm4OPp8PmqZhfX0dX3zxhXRqVKl4xud61OWfRDecpZ2gOlm0E7744gupFQ2Hw9B1XVq+G1/DLwYqNU2TYK7f78fAwICMYxgdHYXNZsPg4CDu3LmDcrmMr7/+GmtrayiXy9jY2OhqCnQWMjivs8EY2LXb7RgdHcXNmzcRj8fxwQcf4Pbt2wAO2FG7u7t4/vw5vv76a6TTaeRyuZ42wynOiHOxlVSbuNFoIJvNot1u48svv5SxMaxJp95ktoszzMhCUO0I0redTifGxsYwOjoqHT4nJyeh6zpGRkYkEMBmKac5O0/rnKU0TRvSdT2padoQgN1Tvs+podYVcbDq4OAgJicnEY/HJTJYLpelvuisOzSq2YvzkIMa8XK73eKRs3B7enoao6OjEi1utVpYWVnB+vo6isUiHj16hCdPnqDRaEgUWX1vGuakN7GBxOjoKKrVKpxOJ5aWlrrqdnpdoxJlOxMZGBUuABQKBSwvL8Pr9UrkPBgMYnJyEuFwuKsLk/p8SeMzXr/RQeNGVA+1k1yvEil3noUMTgOm51kfUi6XEQ6HcePGDVn7+XweyWQSOzs7MjKAkaXX2RPnLQMWKXO4LjsusTA5n89L0we1AchJ70ltlqGOa2DNZrVaRaFQkP2mGvqqDA7X0KWsg1eB654Zd2aguAd6ZWJUQ+9VMj1LvahmTckW4Pwh1gywtqPVaqFarSKZTOLp06cyk6Zer8s+57XRYGc3w0gkgsnJSQwODmJubg5vv/22UHlID+fMPHZwy+fzcl3GLBKzbeexHwwBAMkc08Eol8tYXV3FN998g3Q6jXv37iGTyUggga9VdSAd30ajIQGdubk5MUii0SiazaZkXLTDLPzL9MZl6UXeF+lnwAEVOp/Py4Duubk5NBoNPH/+HKFQCJVKRWRIuZ6VjXBZOkE1nrnfx8bGZJ3zOWYyGTx69EhKJUgDVWldRmNevb/j4DzsBBXUUaVSCcvLy1311p3OwXBxZsaMr2EdPp2zQCAAr9crnR+DwSBmZmYwMzMDv9+PoaEhKQ1gc4t0Oo1UKvVSltF5y+C4UJke/N5ms6Gvrw/z8/NIJBKYnp7GxMQEqtUqtre3US6XZexOOp3ump92eD/H+mzD3jpXncBrIotM13Wsrq4KmyKVSkmWlVRErgMAXWwS1TljMP/GjRtoNBqIxWIYHR1Ff38/KpUKEomEMDcMOvBEOK1z9g8B/BqAv3n47z845fucGFxMjKRys8TjcfT390vbYLUNKmlOL9s4p4HNZlNrsM5cDtw85Mbb7XYZnkeFQ8Mhn89jb28Pm5ub2NzcRKVSkfQzO8gYFwgHTXY6B7M8WPjP+qxisSiU0aMWl8vlUiNtZy4Dbmamn6kINzY2pBsl08tUGKS3kbLI5h4shq5Wq/D5fIjH4xgaGsL09LR0KmR2tRcF6ig4HA60221m5uIA/puzlMFp4Ha7u6gJXKscoql24yJeZ1+clwyM0V9GrRigYFc57oXTZHrUQ4rtxgOBAGKxmHRGZfMAduxS28sTbBjEDC0uUC++CuohzEPG5/MhFArJ7Cvg4CBjp0s6oErnxWN91lnpRfW5MOPRi3LT6XSkpTNb2zO6qw5SV88OHrIDAwNydrD+LhwOizwajYbMMWMzENa1qo6qUTZcm4f1FWeuEyibcrmMzc1N6bpKtsP6+jp2d3dRKBS6mrUYax9U2bAZzP7+PnK5nBhhgUBAghQMjLFGR31vI85bBi8D74u0xkKhgI2NDQAHtemkwHNOYrlcxvr6elfW6Cxg0IvABeoEBnZZkzk2NobBwUH09fXB4XDI802n0xLMU+unetGxaKjy98cF5XCIc7MZm82mBGRUClu1Wu3ZgZrnPPW7atC3Wi1puc/ux6SFskNwIpHA+Pg4XC4XNjY2xH6gzlRldFEyOA54Xex47fF4MDg4iKGhIUQiEbGhS6UStra2kMvlsLu7K9TN0wYveEYedo09N52gOoHM8uu6jkwmI8wQZszZ9Ij2otrxm2ud/7Ijttvt/tZ4nl5JhdfBcVrp/z0cNP9IaJq2CeA3cOCU/b80TfsrANYA/AuvfSXHgKosSEMbHBzE9773PUxNTWFiYkK6+e3u7uKzzz5DLpfDzs7OtzrqnHRhqV0L2YmLXcoOF/ov44zlwE2wu7uLxcVFRCIR6SDXarVkqHalUpHUerFYRC6Xk9lWarcYt9vdpSyq1SqePHkCv9+PhYUFJJNJmfnBLjZUQsBBhy/ebz6fF4OpXC7TGDtTGajPio1N6KBubW0hHo9Li1MO3wYgAyPpiOTzebhcLunwSIfe6/VicHAQo6OjKJfL+OM//uOu6DMjbWoGhlQeUiFYl8SDAEAIB/vjUsDrjUajWFhYkOGadrsd9XodyWQSS0tLSKVSchirWYXjgGtKP+Tq01kitQpnLAMqSda9MCLqdrtl7dMYV6/tOPdDx1/XdXHMSFNYWFhAKBSCw+FAoVAQ6uT6+jp++tOfolAoSPBEXQdKC/VLWwdGUB7spBYIBDAwMCAZKBZpM+uSyWSQTCZliDPQW2eq9SmNRkPe/yz0Ip+Nw+FAIpFAX18fpqenpYMiI/+NRkMyRblcDp999hkeP34shjmzFmpra44k+e53v4uFhQWZ90QDjIPtSQFOp9NYWlrC/fv3JcPCTKMa/CuXyzLGRXFeznw/AAeGwubmJn7wgx/A5XLJ5/JsyOVyQnuks0n0oqs1Gg1sbGxIlrqvrw+JREIaLITDYVkzbLder9cBdOvFSqUi18OAxlnL4GVQjSR2G2232/i93/s9RKNRvP322/ilX/olmXsUCASQy+Xwwx/+UAKUahvws9KLF6kTeN2BQADvvvsu5ubmMDQ0hI8//hiDg4NYX1/HV199hWKxiC+//BJPnjwRQ1XNyFL38nwgywR4YZ+oOpT78Sg5HOqKM7eVCAbp1Iwwr+9lBrOxfor3+vz5c7EzHz9+jFAohI8++ggffvghvF4v3n77bczPz2NlZQXZbBYOhwOlUgnr6+tivBvthPOWwVFQM2aUBWdEhsNhfPTRR/jFX/xFofAzCP6Hf/iH2NzcxPr6utD/jA57L3APAS+6erN74yGl+lx1Ah0q1lEyS0ymg1qHrK7lXvWWqu+haZrod3X0EP+Oa63XDM2T4DjdGv/yEb/6J0/9qa+Bl9FcEokEPB6PKIXd3d0ug83o4Z4EdFCMYLeWTqfzy697b0D3LBvgxayYTCbTNeuLDgqjXk+ePJH2t8buaqqi4nsCkE5fe3t7UgzMBgEsDlWL7dkRjddJhMNhFAoFtNvtM5GBUR4AvjX8lE08mAHi3LZ2uy2d2orFotQbqo0DAMDn8yEajYrjVS6X8eTJE/h8PnG8em0sQ+ctAR3DSqXyRNf13FnL4bigvDiom23H+Qw5lLlYLIpyOul+eJUMqtXqmcuAjgUzZ1SKXAN0NE+jEFVqIjnlnIHDrCODHdVqFeVyGRMTE9IqV91blMGhwX5p60CFqkuoN5k5Y3MANXPGjo081NROnkYoGYFv/fywAdOpdYJ6UHo8HqkBY+aM2N/fR6lUknbGqVQK6XQaQDcNU82OUr8NDg5K7RrbJZPOwgJyZs44J43geysdCbvAzF6pVDrz/cB7IR1PNaBJa2ZwibPPep1/qtGhRphZFmCz2US/cu5bIBAQh4+vVbrVdoGyKZfLF6oXeY4ySFAqlbCxsYF8Po+hoSGhtUajUQAQGhsdj9PgVXrxInRCrww5szsDAwPo7+9HIpHA9vY2MpkM0um0NAtTs8x8LyOdWW0lThnTEOVrXiaHQxmcuZ1AUCerOu9VMO4LlTWj2koApOaVwTz2N6jX6zJDtNFowOv1yh5RAyNer5f22bnJ4CRwu90y72tgYAADAwOw2WxSElQsFrG9vY21tTVks1lJUhynq/NRZwMDgfV6/dx0grqOmckEDpx3Y0ZTzQb38g9U543gmj+qrv1CMmdmgqooOJQ5Go3KQDmHw4FsNotms4nt7W2he/CQOis643nDeJ3lchkrKytSb+X3+6VhB++PA7bplB2n6YnqtKkbjRuQDpCavr1oGGltdrsd/f39GB8fRzwex9jYGPr7+9HpdJBKpbC/v48nT57gwYMHYmRy0Han08H29jb6+vqwt7cnM31Y9D46Ooq5uTkUi0XhH6vXYdb1Q2XD5jFsiz0yMoJEIgGn0ykZxHQ6jXQ6LdSs09RlXRTUSBZpJXSYGPkiFz6Xy0lQAji6oF89dHVdl3bJdrsd8Xgco6OjCAaDmJ2dlZl4bB28sbGBjY0NJJNJqdW7KnpFzZwx4+33+xEKheD1eiXDwmwz5fk6kb+zuGaVKaFmTNmVlfOLarWadJwjvZU6g1lX1qQODQ1heHgY4XAY09PTGBoaEmqP3W6XLPje3h6ePXuGb775RuqWjcbqZYNBK+59GkSkq/K5qk1rjM9UpQAR7XZb5iWyK53NZkM8HsfExAScTidWVla+lZEzE9RnxPoz2gZPnz5FLBaDz+dDLBaDpmno6+tDPB4XlgadFTPfowrqSpvNJgEM1g9NTU3B5XJhc3MTOzs7uH//vmSad3d35YxXaeMsEWEgg3W4ZNVQTtVqFVtbW13jGi76vrmme3WWNe7To/7Pf9XmZ/yXHa8bjQZSqRSSyaSwD0hzjEQisn7MoBtUUEa6rnd1oZyYmMCtW7cQi8UQDodlyPbS0hI2Nzexvb0tMyJZt2UW3fcy9Dq3jA6Z+r16T0Y9qeqAUCiE/v5+aZw3NzcHv98vDDcGO/L5vPR4OK2NdaWcM+BFZxlOLSdPdnh4WHjllUoFT548waNHj0SBnGbOgFnA7JjqdDF9ykyBOpnceKAcdc9qJEx9XzplbLNOOiFpPBcFY8SfBsjMzAzu3r2LeDyOmzdvYmJiAul0Gl9++SXS6TQWFxfx4x//WIo82bnv6dOnEknc2tpCX18fbt26hU8++QRerxcLCwtoNpvIZDIoFApdBb68DjOuHx5ITqcToVAIbrcbIyMjUty7v7+PZDIpNAXWVjDifhUcNKfTib6+Pmlcwyg/HSfWTNAo59o2RslUh6/T6QiN1+fzYW5uDh9//DEikQhu3ryJyclJVKtVLC0tYWlpCclkEg8fPsTOzk6XIWRW2RGqDDjHhoZbX1+fODp0cLa3t6WZhhmcM+AgCxWJRBAKhRAIBODz+brmOJZKJezs7HTVzjCYo2ma0HyDwSAWFhawsLCAQCCAyclJccKp33K5HFZXV1EqlfDll1/iJz/5CWq1mjTVAMxzjjAqzPVOxzUSicDtdqNarSKdTn9r2HIvA1V91mRrsF6Ho0k4L9Hj8WBxcVEoQqojbSaoIxY4JiAUCiGRSCAej+P999/H/Pw8QqEQxsfHsbOzI82F1Iyx2e6rF2gTAEB/fz9u3bqFRCKBd999F3fu3EE6ncYXX3yBdDqNBw8e4Cc/+Yl0rms2m8JEYnZ0eHgYwWBQBtVzXUUiEcnOsgHPj3/8YwmYnXXjtZfdL/9V12CvAEov4/tl1DxjEI+jZ9xuN549eyYlNBy5wcwT6917nT2XBaOcGKAKBoN466238L3vfQ/RaBS6rgtD4Ec/+hG++uqrri6UahD4KoDZK/X+ez1zY+ZMhcrUs9kO5t198MEHiMViQm8FgN3dXQngbm9vY2dnR7J0x52zasSVc84ACOeTURxS8biJSqWSzDtiu+CrDNY4EGr04lUP/SjlqC46GjE0ZNSItFHhXTTUz+W1srsau0s6HA7oui5Rfw6G5Vw3vodKQctkMgAgSkd932azKV0fjddgNhgVD+lHarcqzmPiV6PROHGTh8uGzWaTWX80CJk14D1xvb4qW6wqZZfLJQYJHZZoNCq1ljabTYaNsmmOuo6uivzU/U5KKOvDSE9jsTObDl10BPyo6wZe7H1SDo3ZjF40FTaQYY0dnTvWYZLWSZoNjZe9vT2h9DAb96rmF5cJlZrMeybDwm63S5DK2H3zZWcD14I6+5HBMWYZ6cwaHTuzgVRPNnpiHRyDEqyPo2NLGuhVBAOzPp9PHClmdgqFQlfnZnW0EDvQsW08m6YEg0FEo1GpS+V7knXhdDqlzvK0Rujr3Cv/PevnZTTi1dok7g1mG2k7MeOojuS5bPTKIqqzIsPhMMLhsIwo4hcbILG1PJuqqZRWs0Nt5tGLsnqS9yEjiY3C4vG41CdTT5ZKJWkWRnvkdXClnDMeKg6HA0NDQ5idncX4+LjUmWWzWSwuLiKTyWBzc7NrKOBVSMUeBdUJU43wvr4+hEIhtNttVCoV6TbDgaNHcWF5QHMWSCAQkHlHnU4H6XRa6vXUTMRlgBuJzpPX65WMqd/vR7VaxdraGtbX17G4uIi1tbVvzWpRD4xOp4NarYbnz58jk8kgFArhrbfeQjQahdvtxuzsrBhwnOejFoebEdwXLpdLMkvDw8OIx+MIhULY3d3F6uqqzAPj+rgKUTCufY/Hg76+PgwPDyMUCkn3TnYlZQ0d0J0dU/c+a9WYYfT7/V3zrcbHx3H79m2pwWJ3vufPn2N5eVlaJ191fUIYqTtAN8WjV4TZDFAdEdYhDg4O4tatWxKQYXBieHgYfr8fY2NjuH37NgKBAIaGhjA0NASXyyX1QOze1Wg08PjxY/zsZz9DoVDAysoKSqWSOClme/Zc55xBNjk5Cb/fj9nZWQwMDCCfz+Phw4dCVU0mk13ZDWMdBb+Ab9N9AEhQgo692+1Gq9Xq2ntmkc1RYFDHGIBQmSlXCXxe6tzPhYUF3L17F6FQCNFoVO7T5XLB4/Fgenoa4XAYrVZLgrIOh0OcOTbN8vl84rCpLBuOrWAToIvcF1xjNJhpEzBgwP4AXJc8K0jvpV3Is/2oumuVFh0KhSSgc/PmTYyNjQlLhR1t2fHxZQ2ULhpqtsvn88nMwg8++ADDw8OYnp6WWsuNjQ08ePAAuVxOGpxQVmrTCzPcVy+oSQS32y33Ozg4iIGBg3nfpKuznpi2M0tZ6HhzP9jtdoRCIczNzSESiWB+fh4fffSR7KtisYhyuYzFxUU8fvwYu7u7kjVTRzidZm9cKeeMcDqdGBoawsLCAgYHB6XoOJPJYHFxETs7O9ja2uq5sK4i1AerGuLDw8MYHR2VsQGM8HNhqEqTIMWD6W124err60MsFpMOj6lUqmvmyWVBdc644aLRKIaGhuB0OqXBxfLyMu7fv49nz56JgqRSYrSHhketVsPa2hrsdjv6+vqwubmJZrOJRCKBoaEhBAIBmaGkNmExG6WJUNcEi3rpnHFWE5UtnTMaWGZ1zozUFI/Hg3g8LlQbGsuVSgWFQkEaI/DgNj4rHubMtE9OTqKvrw9DQ0P47ne/i8HBQcRiMQwPD8PhcCCZTCKZTCKVSolzxqYgV9GAI4y1GUcVLqtyM0MEGPh2pJyRe13Xu5wzdlKMRqO4ffs2otEoxsbGcPPmTfh8PvkCumtKdnZ2UCqVsLS0hD/5kz+ROjPSU45TBH/RUA0SOmfRaBQff/wx5ubmkEql4PV6kUqlsLq62jXv05gZNXYYM54dqtypj9nZjfr2KuwLXi/rynoFJ/h3Zoe6hxnACgQCWFhYkI6Cbre7yznzer1IJBJ4//33JWBFQ5Sdmt1utzRYUz+Lc6DYKEd97UXoRdXZIEskGo1ifn5eqLx+vx+adtDki7X4tIvoVLbbbemyCvR+7lzjdrsdw8PDuH37NsLhsDhnZGzRCaRzdsyh9OcONUhps9mkznh4eBjvv/8+ZmZmEA6HEQwG0W63sbm5KeUgbABiZGiZfX+rZRDhcBg+nw+3b9/G22+/DeDAR2ATOQZbU6mUNJEhu0BlXiQSCbz33nsYGRnB7Ows3n//ffh8PjkfWErz2WefoVqtIpVKdXX6Pu2ZcaWcM3WhUMnwcGCamdQtlZJ3HcFMgMfjgaZpUjdDLjjRSwakAapzg6iEdV2XYcZUXpcZKTYah6rCUZuf8HDgv+R89zIu1dfw9zzgjIeMmZWRSvkCDoIWrBHweDxCh61Wq9LR0KzUrFeBz4dRTj47dvFj50Y1s8zoF79nLQVHKCQSCSQSCYRCIYmydTodtFqtLoqs2rXwuuoTs4N7lBFvlabHhjFerxftdhvRaBSVSkWoz+rzZTc+OiF8rsViEdlsVg5cDirtNa/IrKDeolFBYywWi0kgg0En3p+aLWOAR9f1rrEVpHirdEfK/zQ0oYuE6mjy2XP8SjAYhMvlElqS6rDxNVcJKl1NLfngfdBA5++ZMSEtzxjI5GwwPvNOp4NKpSLOWSaT+dbMzLPoUndc0P5hOQKzZ3TOOMOLs/YCgYA4ZxyzUalUvsWKUZ0z0qLZ6ZIdY9n9mQOOK5WKUADVc+iyGUe8BofDIRRVlepqs9mkCQhtBP6/V4DzqsBut8vaCAQCiMfjXR0W+axYS6uOiOA+YPMxjhRhAyFd14XGSDaSWvJwVuUAV8o5U5tCRCIR9Pf3IxAIyCwbHqzkVgNXa0EdBypdi4vO7XYjFArB5XJhaWkJ+/v7KBaLqFQqEvVVqV2RSAQulwuzs7P48MMPJWPAFO/a2pp0O2R9zWXKkZxvFi2TE82akXA4jHK5jOHhYTSbTemqZjxsuek8Hg8SiYTQIdxutzisfA0PI+MwarOsJ7VuhNG9gYEB3LlzB9PT0wgGg9KFa2lpCV9//bWk73tlVM0G1alm8KVcLqNQKIiT5XA4cPv2bVm3m5ubyGQyQnWx2WxyIDGSxsjwwMCARFpJad3b28Pa2hqazSbu37+PBw8eoFAoYG1tTWoTr1KHxqsOVc6kILFmhsOVD9v1w+PxYGhoCPF4HIFAALdu3ZI20TTSSGVhF8pWqyWdS3d3d6VZws7ODlKpVFdNDq/HbFCvqdlsolariXHKWWV3795Fo9FAMpnE/Py8dNcju4TZBTWbPjU1hXfeeQfhcFhYCnx/Bi1okBppkGYB9yudTuqBubk5/MIv/ILUjrAJQiqVwvb2tpQIXHQN1WlBmZN2FYvFMDY2hmg0CofDIc2ffD4f5ufn5UynI0bna39/H/l8Ho1GQ7oT0vk4HJXTZdRypli1WsXOzk5Xt9zzWgt8prRj+vr6pEnD0NCQdJYk9ZIOI2XQarVQq9WEhshmN+q18tynrWC325FIJDA6OioB7VAohGazKZ0M19bW8OjRo66z4rIYW2owibVw4XAY77//Pm7duoWBgQHpRry1tYXFxUUUi0V8/fXXeP78ucx6Vde+mfa1EaqdwO/9fj9GRkYQiURw69YtvP/++wAgTaN2d3fFRmbQgs4+Hfm33npLuoK/9dZbwizb2toSG+H+/fsolUp4+PChDHU/qy7YV8I5UzMEjAL4fD5xSNj2nBEMMxdvvy7UBejxeBAKhRAOhzEzMyPp6cePH0PXdTEuVOeMBd1+vx+Dg4OYn59HNBpFNBoVB4hd/dg44jKNExrpNNBVpappB13YKINYLCYd5orFIjRNE4MMeOHcU1kxysgiXjUDq3Z/MiNU54z0DnYcm5ubQ71el25aW1tbeP78eVfnTcDcCpdQnz07YbXbbXg8Hng8HoyOjmJ/fx/VahU+nw/hcFhaBbM9/sDAALxer3TXcrlcUuzO+hM2CiDtc2VlBd988w3K5TIymYy0R77qFOmrCtaKsAEM1wEzW06nUzqOscaGgTy1XpfBF77P5uYm0uk0tre38dVXX2FnZweNRqOrbtfsBjrXI/U36yKpF9gqfmhoCP39/ahWq3jy5Ak8Hg8ajQZyuRyq1ao0knA6nZicnMTk5KTUVnB4LGua1U6+vAYz7guV1sWgzsjICBYWFhCLxWQcB5u/5PN5cfhPWyty0eD18Uzv7+8XSrvNZpPAhtvtxtDQEABIfTpHLfBZVioV5PN5FAoFPH36VILdOzs7MqSac2PVzKnxrDwvmRmD04lEAoODg5ienpb+AxyNotLx1HOEzhlZJcZsH+9LpTWyDl2d3VUqlZDNZvHs2TOsr69ja2sLOzs7kpU3g96gvez1ejE1NYX33nsPkUgEg4ODCAaDWF1dxdOnT5FOp6Uu3RiYMPv67wWPxyPBl9HRUUxPTwM4yJrn83k5G7jPGaAn2yIajeLdd9/F7du3EQwGMTExgUAggM3NTRmls7S0hJ/97GfSybdUKnXVKb4uroRzRhj57j6fT9LyhBk2xHmCh0Wn05EZT1Q6nG8yMzMj3WTYPIEHjbFInk4KueSk9jB6ZoY6K5VS02q1UCgUsLm5KXQlUjeGh4fhcrmQyWTE8GAHHQBCUQiFQpiYmEAoFMLY2BgikYjMBFJnPBlrzcwI0jdUChJnNaVSKTE81O6bVxGcu5TNZtHf3y91RTQ+SV9hRo0RT9K6SGmjMZ/P58XAZ70agxKVSgXb29soFAoS6LkKRtpJ8KrIXq9mIRcN1biq1+soFotIp9NYWVmRiLdxhAgACa4wE2QEnby9vT2sr68jl8tJZzK1O+FVAddmrVaTdvHJZBKbm5twu90y34363+FwYHBwUOQQj8dRrVZl3zgcDgwPDyMWi8Hv90vDLY6vYACUhrnZ9oUauGJ310gkgqmpKUQiEYyMjCAQCMDtdiOfzyOVSgk1iRmWq6gn9/b2UCqV4HA4sLm5iWg0Ck3ThAGjgsEJBi/ZTj+VSsmZubu7K+wb6kHuD5VmDFycnqBO0HVdApA+n0+CjwzaMwNM/cFzUaXx8kxUgzA8T1Wnjg5suVyGzWaTAEipVMLm5qYEeFSasHq9FwVVD/J+1Hl3fX19ErzM5/OoVCrCEshmszKXy2z7+TRg0MXhcAhdk04TO5mOjo7C5XKhUqkI7ZFMLDaAYZA3l8uhWCxiY2MDKysrKBaL2N3dRa1WEyoj9fBZyc/0zpm66Zn18Hg8QmvkoaE2/lAFdB0WmhGapqHVauH58+col8uYnJzEwsIC7HY7ZmdnMTw8jFarJYpDpWIxCsBF2NfXBwC4d+8eFhcXkc/n8eTJE+Tz+a4o1WXdJ8F0cavVwuPHj1GtViXbRyPsk08+gd1uRzabxdbWlkSFc7mDIfTkiYfDYczNzUkb2b6+PmiahuXlZSwvLyOdTktRp1kNEODFzL94PI7+/n4MDg6K0ZHL5fD5558jm81ieXm5i75hxnvpBV4rD+LV1VXk83l4vV7cvHkTAITaous6bt68KQ1seNCqNYqM+LdaLTE4arUakskkarUaUqmUNJQh1Ylr7iplG18G4/Ubg1m9KCKXZahyX+u6jnQ6LdS83d1dGfxp7LanXrNxrfP3NMpUY5TGllovchWypGp2J5PJoFwuS5F/Lpfr6kZKp4xZtHfffVdYJ41GQ2psSAtmoIPPoNVqIZvNYnNzE7u7u10dcc0APntmddggiQ1hvv/972NkZES6vgLAV199hR/96EfI5/N49uxZ16Dxy6bzHwfq9ZVKJTx58gR+vx/pdBoPHz6UnxudM44cUqm+DFTRaWO9mRoEUeuzL0MvULd3Oh2kUikUi0VpXsYZiPF4XNYtg3IMwPL85xgFtVEKW+GzJotMCspnd3dXulmn02kUi0V8/vnncmaQBqpe60VBzV6yg63L5cLCwgJu374t1LyJiQkUi0U8evQIxWIRDx48wOeffy7jh1S9afa1T/A6Vd1fLpfx/PlzZLNZzMzMYG1tDX6/X2w+1tyxeyNtPf6cTAy/349isYhvvvlGssk/+9nPUCqVkMlkZBYka9fOUmamd85UqAcmF5/aHOJNAo1N0nhI4wwGgxgcHISu68KfVTNnwWAQU1NTEl1yu91SS7CxsSE0Bh68l22gqJQdfnEgd61Wk0nsoVBIUvVsiNFoNCRNDbxwziKRCGZnZ0VJs5FAu91GLpdDPp/vaipjViVF6gSpHH6/X2Zz7e3tycR6du+8yp1LmTljxpiUW7/fLxFi9UChomaEmLRIZkwYCSuXy9jc3JQhxjxoSeEB0BWBfRNglqyBMUrOltitVkvoRb0aEBjXt3EPqw0gcrlc14BV/u1Vet68N9I9W60Wkskk/H4/Wq0WxsfHpTkIM0nhcFgyCWxiwOAn1ztlwCYsNORVA96MUJ0rr9eLaDSKRCKByclJjI+PS1MgOqbMkrOxBalgV01PNptNmcnHBjCkfLMbHXAgHw6cZzaUr+FaULMvwIsgjlEmFy0jY8COGT9mvEKhEIrFojhaTqcTLpcL9XodoVAIXq9XasmM+50BONoEaoZNPUeSySS2t7clc7a1tdVFj78sx0a9HzpozBBFo1HEYjFxSJgx3tnZkbpbtRHOVVv7wLfrb9lvgec8AKk559pmTwIGL8hAogxtNhuKxSJyuRx2dnawsbGB1dVVCe4aS2bOEqZ3zozRUEZ3FhcXu5o47O/vSwE/5xeoHfuuC1Tjk7TDnZ0dfPXVV0gmk4hEIhgYGJDOc5zvoHbyopFKfnmtVsPi4qJQulhfYyaoxre68ZaWllCv1xGNRpHP52V2i9o9igYH56RxcKaxS9eDBw/w6NEjGTxrNFLNJBMqA6/Xi4mJCdy8eROJRAKapqFSqaBUKknWsF6vAzDX9Z8ENpsN+/v7ktlYXV3FF198gUgkgvHxcYyNjXVlzbk+2JmLLXPpeDebTaE1qtHRYrHYc5j1VZWbCjXIoXY6pJ4kBz8QCGBwcBB2u11ozq87TPMsrptUs7Oo52Awj1lRIy3yqj5vGpn7+/tIpVIAgFwuJwE81pqQ6sigFTOJqlOiNonIZrNCd3/8+DG2t7dN0/mVzjlrqdlZb2hoCF6vFyMjI4jFYojFYvB6vWi1Wtje3sajR49Qr9fx4MEDJJNJ6Ux4FbJlRvB6GXRyOBzSuIV2Ap8V5cVuhdwDDN6pAXDVIbtsp6wXGJhXnVKyIdiBj134UqkUfD6fDNJmcId7X82cxeNxBINB1Go1ydirM9M4/1J1cNWmOJchB362ph107h4bG0MgEMCNGzewsLAAj8eDvb09rK6uIplMYmlpCclkUnoLqM//KoMyYCBvf38fq6ur+OyzzxAIBDAxMYFEIiHOq8qI0XUdpVIJ+Xxe7A2egV9//TUymQy2t7clgKHK6zzkZnrnDHhx4zyYC4UC/uAP/gBfffUVgBdZFXbdYre9q5aePQ7Ue1Gnt1cqla5hqxycNzs7KxEAGu6bm5uoVCp49uwZfv7zn6NUKkkURe02Y/y8y4R6WFAJF4tFcSzVmUakKnIGEusoSG3odDrI5XLodDrIZrPY3t5GtVrFw4cPsbS0JIa68fPNALULEwuV3377bXz/+9+XtZ7P55FOp5FMJoULD1zNiBifeavVQiaTEecrl8shGAzi3XffxbvvviuDR202G8rlsnTNymQySCaTUn/I7JtaX8QGAMbunGot61WG+szVVvSsIeFBxe6GMzMziEQiKJfL2N7elujiRQ+kVz9HzWQa9+ZZft5V2x8qVEd2ZWUF6+vr8Pl8ePr0qTSLmJ2dhd/vx9TUFCYnJ7vWuKpjaXg2Gg0sLS3h0aNH0ignlUpJQ4nLlJnajZFrNxAI4E/9qT+F733ve/D7/aLzyQ6p1+t48uQJ/uRP/gT5fB6rq6tYXV3tagJBWVw18Fw8Do7KjmuaJrrU7DKgQ6Q2/TJmvbk+aQMwa8bAvkphZvaMIwjUemSuDzaOYlMU2prqZ10kjGU/mqZhYGAAH3zwAeLxOO7evYuPPvoIzWYTS0tLePz4MdbX1/HTn/5UHDM2P+I9qP9eNfDZk2ljs9lw7949bG5uIhAIYHZ2Vjo1s4meOn4hm80ilUqhWq1KIKpYLOLZs2cSHOazP2+GxZVwzlTQuGBXJeDFAmVGyOxNHF4X9PJ5OHGgXq1Wg8/nE8eDPFp14zGbwHb0bCnK5ge9ajXMAtVJVx1yZhAzmYxEgDgDT515RmOUFLd2u41CoSAd+vL5vGxAM0SFXwaVikHlotJSGo2GFHrTqL7Ke4LPD4DMICPPP5/PyywzANJhkQ1E1A58rL2jXNTBu5SpWdf/WUDtWsb9oDq2LpcLfr8fjUZDmqrY7fZLzZ6p160a5MDJjAgji+I6Pm/ucxqr7XYbTqcTtVpNAlSke7MOjfeuZo3YPr3RaIjBwjolnhNmiLT3onIFAgH09fXB5/MhEAjA4/GgXq9L7ZVap8R7vCqdOV8GtXui8dw3rntj4PqqBrKNOqFXraymadJAqlcHV+DFfFOWCTidzq5B1WrgjtR4vjdxmbLjs1NHTbGTtdfrFRmwbKVYLKJUKnU1PbuKWWMVRpYdnxvb5TcaDaFz7+/vC9NIHaXAcT1strezsyO2sjpOBzj/cocr5ZxRCatGKH8OQCLCqsK+rlAPVHWOx/b2trQaX11dxeeff961gJh14BBJZhUYBeV7m1l2alaPz5vFmqQuhkIhiYKqFEDSPblWarUaisUims2mUBWoiK9CFKnZbGJnZwdPnz7terbPnj2T+UXXYTaXMYuSzWbhcDhw7949pFKproj33t6ezOthhpU1AWz4YHTI+Nqr8MxPA97j3t4eisUi7HY7Njc3sby8LAM2Wa85MTGBcDiMtbU1mSPJBhqXeYirOu+s3uuq74ujQEOFEWSel+VyGW63G8vLy4jFYt+SJf+v1tikUimkUimhwpup1kyl67J2JJfLYWNjoytgUywWsby8jGKxiM3NTayurkqwRnXSrzJ4H2o2SP1dr78/6ndXBapOeFkDI03TupyrXjRNfu3t7Uk9pkpZpF2pGvOXKTv1fjkGg01/ZmdnEY/HAUC6C/785z/H/fv3pdxBbX50XdArCM16S9p2yWQSXq8XT58+FUecNFd2dmw2m5K4UEeGAEfXX541rpRzBrygMlSr1Z41Qdf1sO0F3qsaKa3X69jZ2QFwtGevRpmMbV+vgvx48FBRkmJDutNRRrb6varYjlLqZpWDel3tdhs7OztYXl6WzoPVahUbGxsol8syXPY6UPS41jlrR9M0pNNpLC4ufutvjQdqr/ciroNh9jKoeoKyA4DNzU2Ew2E0m03cvHlT5v+Nj48jGo3i4cOHCAQCACAH22XJyQzG0FUEh7cDQD6fx8bGRleU/WVQux9eVnvwV0F1zmhUZ7NZrK+vw2azoVKpyOzOe/fuyciEQqHQ1SjrKpx7r8J1uIfT4Dj3rWbVXhVcUM+Lo5xas+gj3hOdM4/HI0Om4/E4MpkMVlZWkMlk8POf/xw///nPZfC82tn4KmeMjTBm0NQzj6NGjAEZVY+oeu8om/AinvuVc86Ilxld1y0acFwY0/rHxVVW6uom7NW57SQwk9I9DnjQVCoVZLNZmfvC2oPXlYcZwf19mrV+1Dq/Cs/6LMDoMamclUoFxWIR4XAY5XK5qwufSpECzNPB0cLxYdwrxGmfZS8nplek+qKhZgmZHcxkMrDZbNIqm7QkdvhTA5RX+fyzcLbope+OWhtmWjOq40n6os1mQyaTwe7urnSlJZ3/svfsecPooL3M4XrZexDGOsaLwJV0ztT0vfHnZtowFwWjQ3Ea50z996rgrI2Eo7JsZoIxwlOpVHD//n2srq5KBpHtkRlRM+u9vA5UZfmq534VDteLAnVnq9XC2tqadCSLRqPSEIj0MFJDrws19k3EWet24/tcppGnZv84u7BarWJxcRGrq6uw2WxSU8MBzeyydp0yZhZejeM+4+P8nZnWi2rzcRjy06dP8Tu/8zvw+Xwyk61er8ugbu7Z684YIc7CrrsMOV1J5wx4cx2xo3AVHIvzwJt638wUttttrK+vCw+eCpcDfK8TXYHoVS9g4XigrPb392UGXqvVwrNnz4QW5nA40Gw2Jeugjhaw9O7Vw3V+Xup6rtVq0HVdWmGra5b6UF3D11kuFo7GdXnuKkuMNZf7+/sSZCPFN5vNSoCi1Wodm9Z81XHVbcPr/XQsWHiDcN2pChbOF0fRRK11ZeGqwUjjPA3d34KFqwhLh18PaBf50DRNSwOoAshc0EcmLuizJnRd7zvOH16CDICLkcNJZbCGi3s+uMDPOpYcLBkcwNIJlgwAaz8AlgwIaz9YMgAsWwmwdALw5srgQp0zANA07Qtd1z+8bp91Elz0dVlysGRw0Z91ElgysGRAWHKwZABYMgAsGQCWrURYa+HNk4FFa7RgwYIFCxYsWLBgwYIFE8ByzixYsGDBggULFixYsGDBBLgM5+w3r+lnnQQXfV2WHCwZXPRnnQSWDCwZEJYcLBkAlgwASwaAZSsR1lp4w2Rw4TVnFixYsGDBggULFixYsGDh27BojRYsWLBgwYIFCxYsWLBgAljOmQULFixYsGDBggULFiyYABfqnGma9iuapj3WNG1Z07RfP+P3/tuapu1qmvZA+VlM07T/UdO0p4f/Rs/yM08DSwbnK4PD9ze9HCwZyDVZ+8GSgSUDWDrh8HreeBkcXpO1HywZWPsBb64MLsw50zTNDuA/B/DnANwC8Jc1Tbt1hh/xdwD8iuFnvw7gB7quzwH4weH/Lw2WDC5EBoDJ5WDJ4ADWfrBkAFgyACydAFgyIKz9YMkAsPYD8GbL4CIzZx8BWNZ1fUXX9SaAvw/gV8/qzXVd/xGAnOHHvwrgtw6//y0A/8xZfd4pYcngnGUAXAk5WDI4gLUfLBkAlgwASycAlgwIaz9YMgCs/QC8wTK4SOdsBMCG8v/Nw5+dJwZ0XU8efr8DYOCcP+9VsGRwOTIAzCUHSwYHsPaDJQPAkgFg6QTAkgFh7QdLBoC1H4A3WAZvTEMQ/WBmwBs9N8CSwQEsOVgyACwZAJYMAEsGhCUHSwaAJQPAkgFhyeHyZHCRztkWgDHl/6OHPztPpDRNGwKAw393z/nzXgVLBpcjA8BccrBkcABrP1gyACwZAJZOACwZENZ+sGQAWPsBeINlcJHO2ecA5jRNm9I0zQXgLwH4h+f8mf8QwK8dfv9rAP7BOX/eq2DJ4HJkAJhLDpYMDmDtB0sGgCUDwNIJgCUDwtoPlgwAaz8Ab7IMdF2/sC8Afx7AEwDPAPzHZ/zefw9AEkALB7zUvwIgjoNOK08B/D6A2EXeryWDi5fBVZGDJYPzl4MlA0sGlgyulhwsGZy/HCwZWDK4SnJ4U2WgHV6cBQsWLFiwYMGCBQsWLFi4RLwxDUEsWLBgwYIFCxYsWLBgwcywnDMLFixYsGDBggULFixYMAEs58yCBQsWLFiwYMGCBQsWTADLObNgwYIFCxYsWLBgwYIFE8ByzixYsGDBggULFixYsGDBBLCcMwsWLFiwYMGCBQsWLFgwASznzIIFCxYsWLBgwYIFCxZMgP8/V4+HIxi00vYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1080x576 with 13 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(15,8))\n",
    "for i in range(13):\n",
    "    c = random.randint(X.shape[0])\n",
    "    a = X[c,1:].reshape((20, 20))\n",
    "    a=np.transpose(a)\n",
    "    plt.subplot(1,13,i+1)\n",
    "    plt.title('label '+ str(y[c]))\n",
    "    plt.imshow(a,cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Partie 1, régression logistique\n",
    "Dans cette partie, nous aimerions utiliser la régression logistique pour classifier nos images\n",
    "\n",
    "Rappelons que la régression logistique nous donne la probabilité d'appartenance (oui ou non) à la classe 1 (elle permet une classification binaire).\n",
    "\n",
    "Pour étendre la régression logistique à une multi-classification, nous allons utiliser une stratégie 1 contre tous.\n",
    "\n",
    "Nous paramètres $\\theta$ seront donc une matrice avec un nombre de lignes égale au nombre de classes, et avec un nombre de colones égale au nombre decaractéristiques (chaque ligne de la matrice $\\theta$ correspond aux paramètres d'un classifieur."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Sigmoid(z):\n",
    "    return 1 / (1 + np.exp(-z));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lrCostFunction (X, y, initial_theta, alpha, MaxIter,lambda_):\n",
    "    # dans cette fonction vous devez appliquer la regression logistique avec tout ce que ceci implique\n",
    "    # calcul du coût, minimisation du coût avec descente du gradient, et retour des paramètres theta pour une classe\n",
    "    theta = np.array(initial_theta, copy=True)\n",
    "    for i in range(MaxIter):\n",
    "        y_hat = Sigmoid(X @ theta)\n",
    "        error = y_hat - y\n",
    "        theta_grad = np.expand_dims((X * error).mean(axis=0), axis=1)\n",
    "        theta[0] = theta[0] - alpha * theta_grad[0]\n",
    "        theta[1:] = theta[1:] * (1 - lambda_ * alpha / X.shape[0]) - alpha * theta_grad[1:]\n",
    "    return theta.T\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predictOnveVsAll (all_theta, X):\n",
    "    # ici en utilisant les paramètres calculés par la régression logisitique, \n",
    "    # nous aiemrions retourner les etiquettes predites\n",
    "    \n",
    "    # Ici chaque classifieur retournera de probabilité, il faudra choisir \n",
    "    # la probabilité maximale \n",
    "    # de tous les classifieurs d'une exemple donné\n",
    "    # répeter pour tous les exmemples\n",
    "    y_pred = np.argmax(Sigmoid(X @ all_theta.T), axis=1)\n",
    "    \n",
    "    return classes[y_pred[..., np.newaxis]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 401)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "classes= np.unique(y)\n",
    "number_classes=classes.shape[0]\n",
    "all_theta = np.zeros((number_classes, X.shape[1]));\n",
    "all_theta.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MaxIter= 10000\n",
    "lambda_= 0.1\n",
    "alpha = 0.01\n",
    "# initial_theta pour chaque classifieur\n",
    "initial_theta=np.zeros((X.shape[1], 1));\n",
    "for i in range (number_classes):\n",
    "     # appel pour chaque classifieur\n",
    "     theta = lrCostFunction(X,(y==classes[i]).astype(int),initial_theta,alpha, MaxIter,lambda_);\n",
    "     all_theta[i,:]=theta;               "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = predictOnveVsAll(all_theta, X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Qualité du classifieur RL\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "precision = np.mean(y==y_pred)*100\n",
    "precision\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Partie 2: Réseaux de neurones\n",
    "Pour cette partie, nous choisissons une réseau simple:\n",
    "\n",
    "- une couche d'entrée avec 400 noeuds (20 x 20 pixels) + le biais\n",
    "- une couche cachée avec 25 noeuds\n",
    "- une couche de sortie avec 10 noeuds (nombre de classes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# poids de la couche 1\n",
    "W1 = np.genfromtxt('W1.csv', delimiter=',', dtype=float)\n",
    "W1.shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# poids de la couche 2\n",
    "W2 = np.genfromtxt('W2.csv', delimiter=',', dtype=float)\n",
    "W2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_layer_size  = 400;  \n",
    "hidden_layer_size = 25;   \n",
    "num_labels = 10;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Prédiction**\n",
    "\n",
    "Appliquer une propagation en avant en utilisant les paramètres données pour prédir les classes de l'ensemble d'apprentissage.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict (W1, W2, X):\n",
    "    # appliquer une propagation en avant\n",
    "    # !--- n'oubliez pas d'appliquer la sigmoid à chaque couche afin d'avoir les probabilités ---!\n",
    "    \n",
    "    # prédire la classe en choisissant la probabilité maximale parmi les 10 noeuds de sortie\n",
    "    \n",
    "    h1_output = Sigmoid(X @ W1.T)\n",
    "    \n",
    "    intercept = np.ones((X.shape[0],1))\n",
    "    h2_input = np.column_stack((intercept, h1_output))\n",
    "    \n",
    "    h2_output = Sigmoid(h2_input @ W2.T)\n",
    "    \n",
    "    y_pred = np.argmax(h2_output, axis=1)[..., np.newaxis]\n",
    "    \n",
    "    y_pred = y_pred + 1\n",
    "    y_pred[y_pred==10] = 0\n",
    "    return y_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred=predict(X)\n",
    "precision = np.mean(y==y_pred)*100\n",
    "precision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Vérification de l'implementation\n",
    "Comparer vos algorithmes à ceux de scikitlearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import log_loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = LogisticRegression(penalty=\"l2\", max_iter=10000)\n",
    "model.fit(X[:, 1:], y.squeeze())\n",
    "y_prob = model.predict_proba(X[:, 1:])\n",
    "y_pred = np.argmax(y_prob, axis=1)[..., np.newaxis]\n",
    "\n",
    "np.mean(y_pred==y)*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TP 6 Réseaux de neurones - Multi classification et Propagation en arrière"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rajoutons l'ordonnée à l'origine theta 0\n",
    "intercept=np.ones((data.shape[0],1))\n",
    "X=np.column_stack((intercept,data[:,:-1]))\n",
    "y=data[:,-1]\n",
    "# forcer y à avoir une seule colonne\n",
    "y = y.reshape( y.shape[0], 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('X', X.shape ,' y ', y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transformer y de tel sorte à avoir un vecteur pour chaque exemple\n",
    "\n",
    "Equivalent de tocategorical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "YY=np.zeros((y.shape[0], int(np.max(y))+1))\n",
    "YY.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range (YY.shape[0]):\n",
    "    YY[i,int(y[i])]=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15,8))\n",
    "for i in range(13):\n",
    "    c = random.randint(X.shape[0])\n",
    "    a = X[c,1:].reshape((20, 20))\n",
    "    a=np.transpose(a)\n",
    "    plt.subplot(1,13,i+1)\n",
    "    plt.title('label '+ str(y[c]))\n",
    "    plt.imshow(a,cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Réseaux de neurones\n",
    "Pour cette partie, nous choisissons une réseau simple:\n",
    "\n",
    "- une couche d'entrée avec 400 noeuds (20 x 20 pixels) + le biais\n",
    "- une couche cachée avec 25 noeuds\n",
    "- une couche de sortie avec 10 noeuds (nombre de classes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# poids de la couche 1\n",
    "W1 = np.genfromtxt('data/neural_nets/W1.csv', delimiter=',', dtype=float)\n",
    "W1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# poids de la couche 2\n",
    "W2 = np.genfromtxt('data/neural_nets/W2.csv', delimiter=',', dtype=float)\n",
    "W2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_layer_size  = 400;  \n",
    "hidden_layer_size = 25;   \n",
    "num_labels = 10;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO  1 calcul du coût\n",
    "**Tâche 1**: Modifier la fonction computeCost afin d'obtenir un coût avec régularisation\n",
    "\n",
    "Rappelons que le coût avec régularisation est calculé comme suit:\n",
    "\n",
    "$ J(\\theta)= \\frac{1}{m} \\sum_{i=1}^m\\sum_{k=1}^K \\left[ -y_k^{(i)} \\log (h_\\theta(x^{(i)}) )_k - (1-y_k^{(i)}) \\log (1-h_\\theta(x^{(i)}) )_k \\right] + \\frac{\\lambda}{2m} \\left[ \\sum_{j=1}^{25}\\sum_{k=1}^{400} \\left( \\theta_{j,k}^{(1)} \\right)^2 +   \\sum_{j=1}^{10}\\sum_{k=1}^{25} \\left( \\theta_{j,k}^{(2)} \\right)^2 \\right] $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Sigmoid(z):\n",
    "    # pour une valeur donnée, cette fonction calculera sa sigmoid\n",
    "    return 1 / (1 + np.exp(-z));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def computeCost(X, YY, theta1, theta2, hidden_layer_size, lambda_):\n",
    "    a1 = X\n",
    "\n",
    "    z2 = np.zeros((X.shape[0], hidden_layer_size))\n",
    "    a2 = np.zeros((X.shape[0], hidden_layer_size))\n",
    "\n",
    "    intercept = np.ones((X.shape[0], 1))\n",
    "    z2 = np.column_stack((intercept, np.matmul(a1, np.transpose(theta1))))\n",
    "    a2 = Sigmoid(z2)\n",
    "    a2[:, 0] = 1\n",
    "\n",
    "    z3 = np.zeros((YY.shape[0], YY.shape[1]))\n",
    "    a3 = np.zeros((YY.shape[0], YY.shape[1]))\n",
    "\n",
    "    z3 = np.matmul(a2, np.transpose(theta2))\n",
    "    a3 = Sigmoid(z3)\n",
    "\n",
    "    J = np.zeros((YY.shape[0],1))\n",
    "    J = (1 / y.shape[0]) * (np.sum(np.sum((-YY * np.log(a3)) - ((1 - YY) * np.log(1 - a3))))) + \\\n",
    "        (lambda_ / (2 * X.shape[0])) * (np.power(theta1[:, 1:], 2).sum() + np.power(theta2[:, 1:], 2).sum())\n",
    "\n",
    "    return J"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO  2 calcul du gradient\n",
    "**Tâche 2**: Implémenter la fonction NNCostFunction afin de retourner:\n",
    "\n",
    "- Le coût avec régularisation\n",
    "- Le gardient du coût par rapport à chaqu'un des paramètres du réseau\n",
    "\n",
    "Rappelons le gradient de la sigmoid:\n",
    "\n",
    "$ sigmoid(z) =\\sigma(z)= \\frac{1}{1-e^{-z}}$\n",
    "\n",
    "$ \\sigma(z)^\\prime = \\sigma(z)(1-\\sigma(z)) =a(1-a)$ tel que $a$ est l'activation d'une couche donnée\n",
    "\n",
    "L'algorithme à implementer est comme suit:\n",
    "\n",
    "Pour chaque exemple de l'ensemble d'apprentissage faire\n",
    "\n",
    "Pour chaque noeud de la couche de sortie, calculer la dérivée (gradient):\n",
    "\n",
    "$\\delta_k^{(3)} = \\left( a_k^{(3)}- y_k \\right)  $\n",
    "\n",
    "Pour chaque noeud de la couche caché calculer la dérivée (gradient):\n",
    "\n",
    "$\\delta^{(2)} =  \\left(\\Theta^{(2)}\\right)^T \\delta^{(3)} \\quad .* \\quad \\sigma^\\prime\\left(z^{(2)} \\right) $\n",
    "\n",
    "Notons que l'opérateur $.*$ represente la multiplication élement par élement et non pas la multiplication matricielle\n",
    "\n",
    "Le gradient de chaque noeud et de chaque couche sera finalement:\n",
    "\n",
    "$\\Delta^{(l)} =  \\Delta^{(l)} + \\delta^{(l)} * \\left( a^{(l)} \\right)^T $\n",
    "\n",
    "Fin pour\n",
    "\n",
    "Fin pour\n",
    "\n",
    "Diviser le gradient cumulé par le nombre d'exemples:\n",
    "\n",
    "$ \\frac{\\partial J(\\Theta)} { \\partial \\Theta_{i,j)}^{(l)}} = \\frac{1}{m} \\Delta_{i,j)}^{(l)}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def NNCostFunction(X, YY, theta1, theta2, hidden_layer_size, lambda_ = 0):\n",
    "    a1 = X # 5000x401\n",
    "    z2 = X @ theta1.T # 5000x401 @ 401x25 = 5000x25\n",
    "    a2 = np.column_stack((np.ones((z2.shape[0], 1)), Sigmoid(z2))) # 5000x26\n",
    "\n",
    "    z3 = a2 @ theta2.T # 5000x26 @ 26x10 = 5000x10\n",
    "    a3 = Sigmoid(z3) # 5000x10\n",
    "\n",
    "    delta3 = a3 - YY # 5000x10\n",
    "\n",
    "    delta2 = ((delta3 @ theta2) * a2 * (1 - a2))[:, 1:] \n",
    "\n",
    "    theta2_grad = (np.broadcast_to(delta3[:, :, np.newaxis], (*delta3.shape, a2.shape[1])) \\\n",
    "        * np.broadcast_to(a2[:, np.newaxis, :], (a2.shape[0], delta3.shape[1], a2.shape[1]))).mean(axis=0) \\\n",
    "        + (lambda_ / X.shape[0]) * np.column_stack((np.zeros((theta2.shape[0], 1)), theta2[:, 1:])) \n",
    "\n",
    "    theta1_grad = (np.broadcast_to(delta2[:, :, np.newaxis], (*delta2.shape, a1.shape[1])) \\\n",
    "        * np.broadcast_to(a1[:, np.newaxis, :], (a1.shape[0], delta2.shape[1], a1.shape[1]))).mean(axis=0) \\\n",
    "        + (lambda_ / X.shape[0]) * np.column_stack((np.zeros((theta1.shape[0], 1)), theta1[:, 1:])) \n",
    "\n",
    "    J = computeCost(X, YY, theta1, theta2, hidden_layer_size, lambda_)\n",
    "    \n",
    "    return J, [theta1_grad, theta2_grad]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO  3 entrainement du réseau\n",
    "**Tâche 3**: Entrainer le réseau en utilisant la descente du gradient:\n",
    "Ici il faudra reprendre l'algorithme de la descente du gradient afin de mettre à jour les paramètres du réseau\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradientDescent(X, y, theta, alpha, iterations, lambda_ = 0):\n",
    "    losses, gradients = [], []\n",
    "    for _ in range(iterations):\n",
    "        loss, grad = NNCostFunction(X, y, theta[0], theta[1], hidden_layer_size, lambda_)\n",
    "        \n",
    "        theta[0] = theta[0] - alpha * grad[0]\n",
    "        theta[1] = theta[1] - alpha * grad[1]\n",
    "        \n",
    "        losses.append(loss)\n",
    "        gradients.append((grad[0].mean(), grad[1].mean()))\n",
    "        \n",
    "    return theta, np.asarray(losses), np.asarray(gradients)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_cost_gradient(costs, gradients):\n",
    "    plt.figure(figsize=(15, 5))\n",
    "    plt.subplot(1, 2, 1)\n",
    "    plt.title(\"costs\")\n",
    "    plt.plot(range(costs.shape[0]), costs)\n",
    "    plt.xlabel('iterations')\n",
    "    plt.ylabel('cost')\n",
    "    plt.subplot(1, 2, 2)\n",
    "    plt.title(\"gradients\")\n",
    "    for i in range(gradients.shape[1]):\n",
    "        plt.plot(range(gradients.shape[0]), gradients[:, i])\n",
    "    plt.xlabel('iterations')\n",
    "    plt.ylabel('gardient')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Entrainement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iterations = 500\n",
    "alpha = 1.5\n",
    "theta = [np.random.normal(0, np.sqrt(1/212), W1.shape), np.random.normal(0, np.sqrt(1/17), W2.shape)]\n",
    "\n",
    "theta, losses, gradients = gradientDescent(X, YY, theta, alpha, iterations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_cost_gradient(losses, gradients)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparaison avec les poids W1, W2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "computeCost(X, YY, theta[0], theta[1], hidden_layer_size, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO 4 : Régularisation\n",
    "Mettre à jour le calcul du coût en ajoutant le terme de régularisation et comparer les résultats (avec et sans régularisation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "iterations = 500\n",
    "alpha = 1.5\n",
    "lambda_ = 500 \n",
    "theta_reg = [np.random.normal(0, np.sqrt(1/212), W1.shape), np.random.normal(0, np.sqrt(1/17), W2.shape)]\n",
    "\n",
    "theta_reg, losses_reg, gradients_reg = gradientDescent(X, YY, theta_reg, alpha, iterations, lambda_=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plot_cost_gradient(losses_reg, gradients_reg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "computeCost(X, YY, theta_reg[0], theta_reg[1], hidden_layer_size, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "l'erreur est un peu plus élevée que celle sans régularisation. La précision est plus basse aussi."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TODO 5 : Prédiction\n",
    "\n",
    "Appliquer une propagation en avant en utilisant les paramètres données pour prédir les classes de l'ensemble d'apprentissage.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict(W1, W2, X):\n",
    "    \n",
    "    # appliquer une propagation en avant\n",
    "    # !--- n'oubliez pas d'appliquer la sigmoid à chaque couche afin d'avoir les probabilités ---!\n",
    "    \n",
    "    # prédire la classe en choisissant la probabilité maximale parmi les 10 noeuds de sortie\n",
    "    h1_output = Sigmoid(X @ W1.T)\n",
    "    \n",
    "    intercept = np.ones((X.shape[0],1))\n",
    "    h2_input = np.column_stack((intercept, h1_output))\n",
    "    \n",
    "    h2_output = Sigmoid(h2_input @ W2.T)\n",
    "    \n",
    "    y_pred = np.argmax(h2_output, axis=1)[..., np.newaxis]\n",
    "    \n",
    "    \n",
    "    return y_pred"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La précision du modèle :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = predict(theta[0], theta[1], X)\n",
    "precision = np.mean(y==y_pred)*100\n",
    "precision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La précision du modèle avec régularisation (qui est plus bas que celui du modèle sans):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = predict(theta_reg[0], theta_reg[1], X)\n",
    "precision = np.mean(y==y_pred)*100\n",
    "precision"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = predict(theta_reg[0], theta_reg[1], X)\n",
    "precision = np.mean(y==y_pred)*100\n",
    "precision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# optionnel 1 : Vérification de l'implementation\n",
    "Comparer vos algorithmes à ceux de scikitlearn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import log_loss\n",
    "from sklearn.neural_network import MLPClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = MLPClassifier(solver='sgd', alpha=0.01, activation='logistic', max_iter=5000,\n",
    "        hidden_layer_sizes=(25,), random_state=42)\n",
    "clf.fit(X[:, 1:], YY)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = np.argmax(clf.predict(X[:, 1:]), axis=1)[..., np.newaxis]\n",
    "precision = np.mean(y==y_pred)*100\n",
    "precision"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "clf.score(X[:, 1:], YY) * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optionnel 2 : Visualisation des poids\n",
    "    \n",
    "La visualisation de poids permet de voir quel partie du réseau est activé et pour quelle classe.\n",
    "Il est possible de visulaiser les paramètres theta1 . \n",
    "Ceci peut se faire en utilisant un reshape de ces paramètres afin d'avoir 25 images de taille 20x20 (n'oubliez pas d'ignorer la premuère colonne, celle du biais)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def visualize_weights(w1, w2):\n",
    "    plt.figure(figsize=(50, 20))\n",
    "    for i in range(w1.shape[0]):\n",
    "        plt.subplot(3, 13, i + 1)\n",
    "        plt.imshow(np.reshape(w1[i][1:], (20, 20)))\n",
    "    for i in range(w2.shape[0]):\n",
    "        plt.subplot(3, 13, w1.shape[0] + i + 2)\n",
    "        plt.imshow(np.reshape(w2[i][1:], (-1, 1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "visualize_weights(theta[0], theta[1])"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "4d48602314cfe5a9900ec6b938bb325d180d9d20e0954f86c002443dd920ba8a"
  },
  "kernelspec": {
   "display_name": "Python 3.8.8 ('base')",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}