{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["# 2A.data - Calcul Matriciel, Optimisation - correction\n", "\n", "Calcul matriciel ([numpy](http://www.numpy.org/))."]}, {"cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["Populating the interactive namespace from numpy and matplotlib\n"]}], "source": ["%matplotlib inline\n", "import matplotlib.pyplot as plt"]}, {"cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [{"data": {"text/html": ["Plan\n", "
run previous cell, wait for 2 seconds
\n", ""], "text/plain": [""]}, "execution_count": 3, "metadata": {}, "output_type": "execute_result"}], "source": ["from jyquickhelper import add_notebook_menu\n", "add_notebook_menu()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Exercice 1: Echiquier et Crible d'Erathosthene"]}, {"cell_type": "markdown", "metadata": {}, "source": ["* Exercice 1-A Echiquier: Cr\u00e9er une matrice \u00e9chiquier (des 1 et des 0 altern\u00e9s) de taille 8x8, de deux fa\u00e7ons diff\u00e9rentes\n", " * en vous servant de slices \n", " * en vous servant de la fonction [tile](http://docs.scipy.org/doc/numpy/reference/generated/numpy.tile.html#numpy.tile)\n"]}, {"cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["[[1 0 1 0 1 0 1 0]\n", " [0 1 0 1 0 1 0 1]\n", " [1 0 1 0 1 0 1 0]\n", " [0 1 0 1 0 1 0 1]\n", " [1 0 1 0 1 0 1 0]\n", " [0 1 0 1 0 1 0 1]\n", " [1 0 1 0 1 0 1 0]\n", " [0 1 0 1 0 1 0 1]]\n", "[[1 0 1 0 1 0 1 0]\n", " [0 1 0 1 0 1 0 1]\n", " [1 0 1 0 1 0 1 0]\n", " [0 1 0 1 0 1 0 1]\n", " [1 0 1 0 1 0 1 0]\n", " [0 1 0 1 0 1 0 1]\n", " [1 0 1 0 1 0 1 0]\n", " [0 1 0 1 0 1 0 1]]\n"]}], "source": ["import numpy as np\n", "#Exo1a-1:\n", "chess = np.zeros((8,8), dtype=int)\n", "chess[::2,::2] = 1\n", "chess[1::2,1::2] = 1\n", "print(chess)\n", "\n", "#Exo1a-2:\n", "chess2 = np.tile([[1,0],[0,1]], (4,4))\n", "print(chess2)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["* Exercice 1-B Pi\u00e8ge lors d'une extraction 2d:\n", " * D\u00e9finir la matrice $M = \\left(\\begin{array}{ccccc} 1 & 5 & 9 & 13 & 17 \\\\ 2 & 6 & 10 & 14 & 18 \\\\ 3 & 7 & 11 & 15 & 19 \\\\ 4 & 8 & 12 & 16 & 20 \\\\ \\end{array}\\right)$\n", " * En **extraire** la matrice $\\left(\\begin{array}{ccc} 6 & 18 & 10 \\\\ 7 & 19 & 11 \\\\ 5 & 17 & 9 \\\\ \\end{array}\\right)$"]}, {"cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["[[ 1 5 9 13 17]\n", " [ 2 6 10 14 18]\n", " [ 3 7 11 15 19]\n", " [ 4 8 12 16 20]]\n", "WRONG: [ 6 19 9]\n", "########\n", "[[ 6 18 10]\n", " [ 7 19 11]\n", " [ 5 17 9]]\n", "(array([[1],\n", " [2],\n", " [0]]), array([[1, 4, 2]]))\n", "[[ 6 18 10]\n", " [ 7 19 11]\n", " [ 5 17 9]]\n"]}], "source": ["#Exo1B:\n", "M = np.arange(1, 21).reshape((4,5), order='F')\n", "print(M)\n", "\n", "idx_row = [1, 2, 0]\n", "idx_col = [1, 4, 2]\n", "#the following line is wrong: it create couples from the two lists\n", "print(\"WRONG:\",M[idx_row, idx_col])\n", "print(\"########\")\n", "# first correct way:\n", "print(M[idx_row][:,idx_col])\n", "# we can also use broadcasted indices to create all the couples we want: \n", "idx = np.ix_(idx_row, idx_col)\n", "print(idx)\n", "print(M[idx])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["* Exercice 1-C Crible d'Erathosthene: On souhaite impl\u00e9menter un [crible d'Erathosth\u00e8ne](http://fr.wikipedia.org/wiki/Crible_d'%C3%89ratosth%C3%A8ne) pour trouver les nombres premiers inf\u00e9rieurs \u00e0 $N=1000$.\n", " * partir d'un array de bool\u00e9ens de taille N+1, tous \u00e9gaux \u00e0 True.\n", " * Mettre 0 et 1 \u00e0 False car ils ne sont pas premiers\n", " * pour chaque entier $k$ entre 2 et $\\sqrt{N}$: \n", " * si $k$ est premier: on passe ses multiples (entre $k^2$ et $N$) \u00e0 False\n", " * on print la liste des entiers premiers"]}, {"cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["[ 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61\n", " 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 151\n", " 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251\n", " 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359\n", " 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463\n", " 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593\n", " 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701\n", " 709 719 727 733 739 743 751 757 761 769 773 787 797 809 811 821 823 827\n", " 829 839 853 857 859 863 877 881 883 887 907 911 919 929 937 941 947 953\n", " 967 971 977 983 991 997]\n"]}], "source": ["#Exo1c\n", "n = 1001\n", "is_prime = np.ones(n, dtype=bool)\n", "is_prime[:2] = False\n", "\n", "for k in range(2, int(np.sqrt(n))+1):\n", " is_prime[k**2::k] = False\n", "print(np.arange(n)[is_prime])"]}, {"cell_type": "markdown", "metadata": {}, "source": ["On remarque que le nombre 6 est barr\u00e9 deux fois car il est multiple de 3 et de 2. Cela signifie que le nombre 6 est barr\u00e9 durant les deux premi\u00e8res it\u00e9rations. En fait chaque nombre ``k*i`` est n\u00e9cessaire barr\u00e9 dans une pr\u00e9c\u00e9dente it\u00e9ration si ``i"]}, "metadata": {}, "output_type": "display_data"}, {"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEPCAYAAABFpK+YAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8U9X7wPHPTdKVNt0tUAqUPSp7CMgeIoiAIJWhgF8U\nBXHwFUVQBFRUEEWU7QIVB6iAosIXhIqAIBX4CWXIKkKBthS6V8b5/VFIrYyW0jZt87xfL140Nyf3\nPHmSPj059+ReTSmlEEII4TR0jg5ACCFE6ZLCL4QQTkYKvxBCOBkp/EII4WSk8AshhJORwi+EEE5G\nCr8QN+nMmTN0794dLy8v9Hq9w+L45Zdf0Ov1nD179obtZsyYQb169UolJp1Ox+eff14qff1bzZo1\nee211xzSd3njVIX/oYceQqfTodPpMBgMVKtWjZEjRxb4i1MUK1asQKcrXHqvxLRu3bqr7rv33nvR\n6XSMGTOmuEMssu+++46OHTsSEBCAl5cXdevW5cEHHyQtLc3RoZWK1157jQsXLvDnn39y7ty5a7b5\n5Zdf7K+rTqcjMDCQbt26sW3btmKL44477uDcuXOEhIQAsH37dnQ6HX///Xe+ds8++yw7d+4stn5L\n06lTp9DpdOj1+nz5/Oe/bt26ARAVFcWECRMcHHH54FSFH6BTp07ExcVx+vRpvvjiC/bu3UtERESx\n96OUQtO0QrevUaMGH3zwQb5t586d48cff6R69erFHV6Rbd68mUGDBtGrVy+2b9/O/v37WbhwId7e\n3mRnZzs6PAAsFkuJ7v/o0aO0adOGWrVqERwcfN12mqaxb98+zp8/z+bNm/Hw8KB3795XFeaiMhgM\n+fq/3nvOaDTi7+9fLH2WturVq3P+/HnOnTvH+fPneeeddzAYDMTFxXH+/HnOnz/Pt99+C0BAQAAe\nHh4OjricUE5k1KhRqmfPnvm2vffee0qn06nU1FT7ttTUVDVmzBgVFBSk3NzcVKtWrdT//ve/fI87\ncuSI6tOnj/Ly8lJeXl7qnnvuUceOHVNKKRUZGak0TVM6nc7+/0MPPXTduDRNUy+//LJyc3NTZ8+e\ntW9/9dVXVY8ePVTXrl3VI488ku8x7777rmrQoIFyd3dX9erVUzNnzlQWi8V+/+eff65uv/125ePj\nowIDA9Xdd9+t/vrrL/v9MTExStM0tXLlStW3b19lNBpVrVq11LJly26Yw6efflq1bt36hm2UUmrz\n5s2qSZMmyt3dXTVt2lRt2bJFaZqmVqxYka//7du353tcnTp11IwZM+y3582bp5o1a6a8vLxU5cqV\n1ZAhQ9S5c+fs91/J9Q8//KA6dOigPDw81OLFi5VSSkVFRak777xTeXl5qaCgIDVw4EB16tSpG8Zd\n0Gtf2Nc1MjJS6XQ6FRsba98WGxurNE1T77//vlJKKbPZrCZNmqSqVq2qXF1dVaNGjdTnn3+ebz/v\nv/++atiwoXJ3d1f+/v6qc+fO9n1eee6xsbH2fF6JTdM01bVrV6WUUtOmTVN16tTJt99ly5apRo0a\nKVdXVxUaGqpefPHFfO+fLl26qIcffli98sorqnLlysrf31+NGDFCpaen3zB/mqapefPmqUGDBilP\nT09VtWpVNW/ePPv9o0aNUnfeeedVj+vatat6+OGHb7jvK3G7uLhc876wsDA1c+bMfLenTp2qxo4d\nq3x9fVVwcLBasGCBys7OVk888YTy8/NTVatWVfPnz8+3n7S0NPXkk0+qqlWrKqPRqFq0aKG+/fbb\nAmMrT5y68MfGxqpOnTopFxcXlZGRYd9+3333qZo1a6qNGzeqw4cPq6eeekq5urqqI0eOKKWUyszM\nVNWrV1c9evRQe/fuVXv27FFdu3ZVdevWVWazWeXk5KgFCxYonU6n4uPjVVxcnEpJSbluXFcKYs+e\nPe1vXJvNpmrWrKm++uor1aVLl3yFf9q0aSosLEytXbtWxcTEqJ9++knVqFFDvfTSS/Y2y5YtU+vW\nrVMnT55U+/btU/3797fHp1Re4a1du7b6+uuv1fHjx9WUKVOUwWBQR48evW6ss2bNUn5+fur333+/\nbpuzZ88qT09PNXr0aHXo0CG1adMm1aRJE6XT6fIVfp1OV2Dhf/fdd9XPP/+sYmJi1M6dO9Udd9yh\nunTpYr//SvFr2LChWrdunYqJiVGxsbHq4MGDysvLS82YMUP99ddf6sCBAyoiIkLVr19fZWdnXzf2\ngl77uLg41b59e/XAAw+o+Pj4676u1yr8iYmJStM0tWDBAqWUUhMnTlSBgYHqm2++UUePHlWvvfaa\n0ul0avPmzUoppf744w9lMBjUZ599pv7++2914MAB9eGHH+Yr/Ff6sFqt6rvvvlM6nU798ccfKi4u\nTl26dEkppdT06dNV3bp17XGsW7dO6fV6NWvWLHX06FG1cuVK5efnl+/906VLF+Xn56f++9//qiNH\njqiNGzcqf3//fG2uRdM0FRAQoBYsWKCOHj2q3n33XWUwGNR3332nlFLqt99+U3q9XsXExNgfc/To\nUaXT6dTu3btvuG+lbr7w+/n5qblz56rjx4+rmTNnKk3TVJ8+fezbXn/9daXT6dShQ4fyPfeuXbuq\nHTt2qJMnT6r3339fubm52V+XisDpCr/BYFBeXl7KaDTaR0jPPfecvc2xY8eUpmlq/fr1+R7bokUL\nNXr0aKWUUh988IHy9PRUFy9etN8fFxenPDw81KeffqqUUuqzzz5TOp2uUHFdKfwrV65UNWvWVEop\n9dNPP6ng4GBlNpvzFf6MjAxlNBrVhg0b8u3jk08+Ub6+vtft40rR2bFjh1Iqr/C/88479jZWq1WZ\nTCa1dOnS6+4nIyND9e/fX+l0OlWlShXVv39/NW/ePJWYmGhv88ILL6iwsDBltVrt29atW1ekEf+/\n7dmzR+l0OvsnoyuF/8p+rxg1apQaOnRovm1ZWVnKaDSqtWvXXnPfhXntlVJX/SG+ln8X/pSUFPXw\nww8rV1dXdfDgQZWRkaHc3Nzsn06uuPfee1X37t2VUkqtXr1a+fr65vs0eqM+tm3bpnQ63VWfav5d\n+Dt27KiGDBmSr828efOU0Wi0Dwy6dOmimjVrlq/N2LFjVfv27W/4vDVNUyNHjsy3bdiwYapTp072\n202aNFFTp061337++eev6ut6brbw33vvvfbbNptNeXt7q379+uXb5ufnZ/9jvGXLFuXh4XHVH/T/\n/Oc/+fZV3jndHH/btm35888/2b17Ny+99BLt2rXjlVdesd9/8OBBNE2jY8eO+R7XqVMnoqOj7W0a\nNWqEn5+f/f7g4GDq169vb1MUAwYMICMjg40bN/L+++8zcuRIDAZDvjbR0dFkZmYyaNAgTCaT/d+j\njz5KamoqiYmJAOzbt4+BAwdSq1YtvL29qVGjBpqmcerUqXz7a9q0qf1nnU5HcHAwcXFx143Rw8OD\nNWvWcPLkSd544w1CQ0N5/fXXqV+/PkeOHAHg0KFDtGnTJt/B7Q4dOhQpJ5GRkdx1111Ur14db29v\n++vyz+ehaRqtW7fO97jdu3ezevXqfDkKDAwkOzubo0ePXrOvwrz2N0MpRf369TGZTPj6+rJx40Y+\n+eQTGjZsyLFjxzCbzVf11blzZ3tfPXv2pGbNmoSFhTF06FDef/99++t7K6Kjo6/Zb1ZWFsePH7dv\n++d7AyAkJOSG740r2rZtm+/2HXfckS9/jz76KB9//DFKKaxWK8uXLy+xxQv/fA6aphEUFESTJk3y\nbQsODiY+Ph7IPUCcnZ1NSEhIvvfOihUrOHbsWInE6AiGgptULB4eHtSsWROA6dOnc+zYMcaPH8/S\npUsdHBm4uLgwatQoZs6cyc6dO9m/f/9VbWw2GwBff/01devWvep+f39/MjMz6dWrFx07dmTZsmVU\nqlQJgEaNGpGTk5Ovvaura77bmqbZ+7iR6tWrM2LECEaMGMHMmTOpW7cus2fP5sMPPyzUc73yR0H9\n6+SwZrPZ/vPp06e5++67GTlyJNOmTSMwMJDTp0/To0ePq56Hp6dnvts2m40HH3yQyZMnX9VHQEBA\noWK8VZqm8b///Y/KlSvj7++Pj49Pvvv/Hde/eXp68scff7B9+3Y2bdrE4sWLee6559i8eTPNmzcv\n9nj/HU9R3xsFefDBB3n++ef54YcfsFgspKSkMHz48Fve77W4uLjku61p2jW3XXleNpsNX19foqKi\nCsxHeeZ0I/5/mz59Oh9//DF79uwBIDw8HICtW7fma7d161YaN25sb3Pw4EEuXrxovz8uLo4jR47Y\n21x5kxT0y/1vjzzyCNu2baNdu3bXLOzh4eG4u7tz/PhxatWqddU/TdM4dOgQFy5cYObMmXTq1In6\n9euTmJh407EUlo+PD5UrV7aPmho1asTvv/+er79/L2MMCgoCyLeUNj4+ntjYWPvt3bt3k5WVxdy5\nc+35OH/+fKFWS7Vq1Yo///yTmjVrXpWjfxfgK2702t92220F9nktNWrUoGbNmlf1WadOHdzc3K7q\nKzIyMl9fmqbRoUMHpk+fzh9//EGVKlWuu07+ynvOarXeMKbw8PBr9ms0Gqldu3ahn9v1/Hvp6Pbt\n22nUqJH9tslkYsiQISxdupQPPviAwYMH4+3tfcv9FodWrVqRlJREZmbmVe+b0NBQR4dXbJxuxP9v\nderU4Z577mHKlCmsX7+eWrVqcd999zFu3DgWL15MjRo1WLhwIdHR0XzxxRcADBs2jJdffpn777+f\n2bNnY7PZmDhxItWqVbMvDb3yqWLt2rV06NABDw+Pq0al11K7dm0uXLiAu7v7Ne/39PRkypQpTJky\nBYAePXpgsVjYv38/e/fu5Y033qBGjRq4ubnx7rvv8swzz3Dy5EkmT55c6O8V3MiMGTNIS0vj7rvv\nJiwsjLS0NJYtW0Z0dLR9DfXYsWOZO3cujzzyCBMnTiQ2NpYXX3wxX8F2d3fnjjvuYPbs2dSvXx+z\n2cyLL76Y73nXrVsXTdOYM2cOw4cPZ9++ffmm5a641h+0KVOmcPvtt/PAAw/w1FNPERQUxMmTJ1m7\ndi1PP/00YWFhVz2mMK/9zbjRH1oPDw+efPJJpk6dSmBgIE2bNmXVqlV8//33bNq0Ccj9vsSJEyfo\n1KkTQUFBREVFcebMGfsfqH/3UaNGDXQ6HT/++CMRERG4ublds6BOnjyZfv36MWvWLAYOHMjevXuZ\nMWMGEydOvGpqsSjWrVvHggUL6NWrFz/99BOrVq3i66+/ztdmzJgxtGvXDk3T+OWXX265z+LSrVs3\nevTowcCBA5k1axZNmjTh0qVL7NixAw8PD0aPHu3oEIuHQ44sOMi1lnMqpdSOHTuUTqdTv/zyi1Iq\nd0nfY489poKDg5W7u7tq3bq12rRpU77H/PXXX+ruu+9WJpNJmUwm1a9fP3X8+PF8bSZMmKAqVapU\n4HLOf652uZZrLef88MMPVfPmzZWHh4fy9/dXbdu2zXeg8JtvvlH16tVTHh4eqkWLFmrr1q3KxcVF\nLV++XCl1/VU1devWveHB1S1btqghQ4aomjVrKg8PDxUUFKQ6dOigvvjii3zt/rmcs3Hjxlct51Qq\ndzVHly5dlJeXl6pXr55avXr1Vf0vXLhQVa9eXRmNRtWxY0e1YcOGfK/VtVbPXHHgwAE1YMAA5e/v\nr4xGo6pbt6569NFH7atdrqUwr/21Xo9/u1FcV5jNZjV58mQVGhqq3NzcVHh4uPryyy/t92/dulV1\n69ZNBQcHKw8PD1WvXj01e/bsG/bx5ptvqtDQUGUwGOzLOf99cFep3MUAjRo1Um5ubio0NFRNnTo1\n38H4az3HV1991b744Hp0Op2aN2+eGjBggDIajSokJCTfAoJ/at68ubrttttuuL9/u9HB3Zo1a+Y7\nuPvv20pd+/3dsGHDfAebs7Ky1OTJk1WtWrWUm5ubqlKliurdu7fasmXLTcValmlK3fjzv9lsZtq0\naVgsFqxWK23btmXw4MGsWrWKn3/+2f4RdujQoTRr1qxU/liJ8kmn0/HZZ58xbNgwR4ciHMxisRAW\nFsbzzz/P+PHjHR2O0ynwc52LiwvTpk3Dzc0Nm83G1KlT7QeW+vbtS9++fW+qw+jo6HwfVZ2Z5CKP\n5CJPRc6FUoqEhASWLFlCRkYGo0aNumH7ipyLm1WcuSjUpK+bmxuQO/r/54GjAj4sXNOtLHesaJwt\nFzc6KOtsubiRipyLv//+m8qVK7NkyRI+/vhjvLy8bti+IufiZhVnLgp1JMdms/H8888TFxdHr169\nqFOnDnv37mX9+vVs3bqV2rVrM2LECIxGY7EFJiqeglabiIqvRo0axbIkVNyaQo34dTods2fPZtGi\nRRw7dowzZ87Qq1cv5s+fz5tvvomvry/Lly8v6ViFEEIUgwIP7v7b119/jbu7e765/YSEBGbNmsWc\nOXOuah8dHZ3vI0pJnAlTCCGcwcqVK+0/h4eHF3nOv8CpnpSUFAwGA0ajkZycHPbv30///v1JSkrC\n19cXgF27dlGtWrVrPv5awZXE+e/LI5PJRGpqqqPDKBMkF3kkF3kkF3lCQkKKbeBcYOFPSkpiwYIF\n2Gw2lFK0b9+eFi1aMH/+fGJiYuznvyhLFwoRQghxfTc91VMcZMSfS0YzeSQXeSQXeSQXea5caa04\nOP25eoQQwtlI4RdCCCcjhV8IIZyMFH4hhHAyUviFEMLJSOEXQggnI4VfCCGcjBR+IYRwMlL4hRDC\nyUjhF0IIJyOFXwghnIwUfiGEcDJS+IUQwslI4RdCCCcjhV8IIZyMFH4hhCjDUlK0Yt+nFH4hhChj\nrlweSyno0iWY06f1xbp/KfxCCFGGPPaYH7t2uQKgv3SRJ+r/gPb0S8XaR4HX3BVCCFFyPvnESGCg\njT59stCfPs29ul8Ie+ZLKmXsRB8fzwv2lkuKrU8p/EIIUYr27HHhzBk9/fqkYThxgub/F03GziME\nvbkBl7/+4qF/tLUZjVgaNsTcoAGexRiDFH4hhChBiYk69u1zoXv3bADUuXgyZ6yl0rR56OPj6fWP\ntjaTCXOjRqS3aod5yCCsNWqAPnd+Xwq/EEKUUVYrHD5sIDzcAkB6usbUp3T0mboK7w+Wcs/Bg3lt\n/fzIadeOc371MfVrh/n2NuDiUuIxSuEXQohblJ0Nbm65P2dladx7byB/bDpC8G/rafLLLxy7tA7D\nf6329rGNu2EZPRyXe3uAwYA7YC7FeAss/GazmWnTpmGxWLBarbRt25bBgweTlpbGO++8Q0JCAsHB\nwUyYMAGj0VgaMQshRJnSq1cQH3xwifocptL333OAHwhrd8h+v03TkVLzNtSoCDIGD0bz9qbkx/XX\npyl1ZcXo9WVnZ+Pm5obNZmPq1Kk89NBD7Ny5E5PJRP/+/VmzZg3p6ekMHz68UJ2ePXv2lgOvCEwm\nE6mpqY4Oo0yQXOSRXOQpq7l45RVvOnfOplPHLFyiojjyzCrCE7bil3La3said8XatjXZHTqQOWgQ\n1qpVb6nPkJCQWw3brlBTPW6XP8OYzWas1tyPK1FRUUyfPh2ALl26MH369EIXfiGEKE82bXLDZoM7\n78w9QNs86zeqvPgZlc9/gy49naDL7WyenmT17k3mwIFkt2sHrq6OC/oGClX4bTYbzz//PHFxcfTq\n1Ys6deqQnJyMr68vAL6+viQnJ5dooEIIUVpiY/UcP66nU6ccyMrC5//2cOnbXfiu2Y/btm2MSUy0\nt7X5+pIeEUFW9+6YW7YEDw8HRl44hSr8Op2O2bNnk5GRwZw5czh9+vRVbTTt2ueTiI6OJjo62n47\nIiICk8lUxHArFldXV8nFZZKLPJKLPKWVC6sVTpzQqFs3d+Y741wq257ZxqA2n2FYt47+ltwVOsTk\n/qd8fbk4+CHcHx2OrV49ANwv/ytJK1eutP8cHh5OeHh4kfZzU6t6jEYjjRo1Yt++ffj6+pKUlGT/\n38fH55qPuVZwZXHOzhHK6vylI0gu8kgu8pRkLqxW+xJ5zp+xMbZHMutmbMb4xy46rVpF55wcWANK\n07DUrMmxBr3wv6M2WvNGmMPDwcWFbIBSeq1MJhMRERHFsq8CC39KSgoGgwGj0UhOTg779++nf//+\ntGzZksjISAYMGEBkZCStWrUqloCEEKKkWa3Qp50L6/67lkrrV1L599/ZlZoM/829X2kasTXbkt2j\nG8ZHB2KrUgVvwOLQqItPgYU/KSmJBQsWYLPZUErRvn17WrRoQb169Zg7dy5btmwhKCiICRMmlEa8\nQghRJK++6k1Ez7M0O/wNnitWsD82Gp7Juz/NVIkL1ZsQ2K4GGQ88gFa3Lu6AzWERl5xCLecsbrKc\nM5d8pM8jucgjuchzK7nYscMVDw9F86bZuG3bRsyLK2l6ch0GW+5XpawGV474tCb0oXZkDh6MNTS0\nOEMvdqW+nFMIIcq61FSN+HgdtWvnLjk/tzuOoHVfEZy6DMPp0wQANjSyOnQgY9gw0nvehY+HG2nF\nf52TMk8KvxCi3DKb805t89tvrny0yIW1w5fhsXYt4yMj0Wy5EzWW0FAyhgzhzxbDqNmpEpoGxXtp\nk/JFCr8Qolw6dUrP8OEB/Lo1Dtc9f3Dvz9/Sa/f/8Pv9HADKxYUjt/XD+GQE+l4dQaejloNjLiuk\n8AshygWLBcaM8WPBgkt4eED1oDTGXVqKV5dV+BzfD4APkBrWEOujI8i6+25MAQGODbqMksIvhCiz\nvv3WQNOmOoKDbRgMkJKssefLM/Tw2IbpvfeYmBQDSWDz8iJj2DAy+/XD3KwZXOcLpSKXFH4hRJmR\nmqqRk6MREJA7N/+//xk4FePOY3dEYfziCzYe+B7PnQn29ub6DUh7YjxZPXqg5NvOhSaFXwhRZixe\n7EVmpsZLz8XjvmEDM2O347EqkmBzrL2N2S8Aa+Nwsvr2JSMiolQuXFLRSOEXQjjM1q2urF5tZO7c\nJHRxcYz2/J59S6OotOpb9Bcv4n+5nTUwkKwePcgYNix3KkfvzGtybp0UfiFEqUlM1PHRR548+2zu\nl7Ja6fZwfvVGfM/swmPXDipbrdwGkAHm8HCsQ4eScvvtWBo0AJ3OobFXJFL4hRAlRin49VdXOnbM\nQdPA29vG9x9nMEG9T+Wta3Hdu5enAXaA0unI6tyZnNatye7UCXOLFpi8vbHIt5iLnRR+IUSxUgps\ntrzZmHemZBP2wFbCL+3Adc8eDqfuxjAv97QJNk9PkgcPQ93Rhpw2bbAFBjowcuchhV8IUayeesqX\nHt0yGFRzN6Y5c9h5cjO8kne/0ulIbNsDbeQgsjt3Rl3nlO6i5EjhF0Lcks2b3UhK0jFwYCa62Fie\nzlpCrYmfEpR5EgCbqxunglsQ2L8F5saNyWnRAtstXn9W3Bop/EKIm5KcrPHXXy60bp0D2dlUitnP\nqUWH8Vv3E+5btlA5JwcAa0AAmQMGkDZ2LG5VqiAz9WWHFH4hRIFyci5fN9xmw7x+J6de3E3vOhtw\nObCfEJuNngBncy9gktm7N+n97yWnV88ye7FxZyeFXwhxQ1lZcFdbN7ZMWknwh+8RcugQTQD+zJ2v\nt4SFcb52G1y7tMJ2V1dsxXjeeFEypPALIa4ye7aJoUMzCMs+QpVZszh0YQP6ibnnubdWrsz+2+4j\nuXVH6v2nBcpoREfFuSyhM5DCL4TgzJnctZehoVa0zEzq79mAz3efE3xqM5rNhlVn4EhgOyo93ZeM\nIUMIdncnGCj1y/eJYiGFXwgnpVTeSSxXr/bAfPgU0ystwLhiBWPT0nLbuLiQPnQoyROewVSlEhkO\njFcUHyn8Qjih335zZdkyT5a+fhyPtWt5cc1KPA//ab8/p3lzDrUeQvATfVD+/jfYkyiPpPAL4QTS\n0jSWL/fk8cdzR/Kt2Y3L+veptH4dOkvut2izXLxI79EL2/j/YG7WjCBkKqeiksIvRAUVG6ujUqXc\nC5gYjYqopYexHPmIkL0bMZw4QXXApunI6tqVjMGDyerVC9zdHR22KAVS+IWooB5+2J8XnjhLz9TV\neC5bxoYLf8I3uffZfHxIuX84WY+Oxla5smMDFaWuwMKfmJjI/PnzSU5ORtM0evToQe/evVm1ahU/\n//wzPpfPszF06FCaNWtW4gELIa5t+XIjfn42+neJx+OHH1ibuoaqj/6Oi+3yN2mNXsR0G4bPf+4i\np2VLMMi4z1kV+Mrr9XpGjhxJWFgYWVlZTJo0iSZNmgDQt29f+vbtW+JBCiGulpioIyZGT8uWZrTM\nTJqf/xX9W99QOf17tKwsfC+3y769LZmD7yNzwAA8PDzIcWjUoiwosPD7+vri65v7FnJ3d6dq1apc\nvHgRAKXk0I8QpemfSzBjjik2PbqBXtXm4rJvH1VsNnu77PbtybjvPrLuvBPl5+egaEVZdVOf9eLj\n4zl16hR169bl8OHDrF+/nq1bt1K7dm1GjBiB0WgsqTiFcHqpqRp9+way8ZsT+K1dSe+lS7kn4Qwk\n5J46wdywIec7D8Bl1ABs1UIdHa4owzRVyGF7VlYW06dPZ9CgQbRu3ZqUlBRMJhOapvHll19y6dIl\nxo4de9XjoqOjiY6Ott+OiIggVa6oA4Crqys5OfLBGyQX//TPXHz8sQv9+lkIVAkY1q1j/0vf0SL1\nF/TW3CWYyZXrEnv/k1SbcA9UwPX28r7IYzKZWLlypf12eHg44eHhRdpXoQq/1WrljTfeoHnz5vTp\n0+eq+xMSEpg1axZz5swpVKdnz569+UgrIJPJJH8EL5Nc5LLZwMXFhNWcjP7UKdaOj6Jn8reEndqG\ndnkqx4oOS7vbSX/4YbLuvLNCX4tW3hd5Qorx5HeFmupZtGgRoaGh+Yp+UlKSfe5/165dVKtWrdiC\nEsJZrZh6jtsPzKHDhdUYYmIYc3m7cnEhq0sXMrp0J/Puu9EqBzk0TlG+FVj4Dx8+zK+//kr16tV5\n7rnn0DSNoUOHsm3bNmJiYtA0jaCgIMaMGVPQroQQ/3LihJ7Iza6MDfsezw8/ZNLWrfb7bCYT2be3\n5Zege2k8pTOaf+5AS3NUsKLCKPQcf3GSqZ5c8jE2jzPl4sIFHYGBNrTUVGwfriTrrU+obTsGgM3D\ngy3+A6kikB8hAAAeSUlEQVQ8qgN+D/UEDw8HR+tYzvS+KEipT/UIIYqHOdPCpI6nWHbXcoJ++BJd\nejoAyb7V4IlRZAwZQptq1aTYiRIlhV+IErZ4sSfdGv5Nq5/nYfzmG35KSYLLizOy27XjaJ8xpHa5\nkxq1HBuncB5S+IUoZtnZkJysIzjIituvvxLxxWc0Ov4TOpW7KiejcnW2e91J04URWMLD8Qcq3kJM\nUZZJ4ReimP285DwBX33MbebVGGJjCQDMGMi48y4ynxyHuVkzwjVNLlUoHEYKvxC3KCFBx9LFRl7u\nuA7jp58yesMGtMtrJqxVqpD+wAP83nQk9Tv5odc7OFghkMIvRJFcuKDD39+GPiuDat+t5MkPVhCw\n+CAAytWVXfUiyBg6jPojGoNORyMHxyvEP0nhF6IIXr3vNG/5zaDy/21Cy87GH7jgWQ3XcUPIGDaM\nasHBjg5RiOuSwi9EIWzY4I6Li6JnzcOY3n6blUe/td+X07Qp8YNGcbrTIGrUlbkcUfZJ4RfiOtLS\nNLw8bRiOHqXx2v/htfEnKmXsA8Dm4sq64JG0XPkwWlgoOqCGQ6MVovCk8AtxDSd+jGHP8+sZ6/8l\nLkePcmXixuLhSU7fPqT+97+0ql7doTEKUVRS+IUg96yYa96I5UE+xbh5EyGHDtEBIBGsfn5k9erF\nofp98b2vPUZ/N0eHK8QtkcIvnJb+1CmMKz7HZf+fGE6fZvzJk/b7bD4+HKzXlz/qDKTn6y3BxYXi\nO1OKEI4lhV84HZe9e/FavBi3H360f5sWwOxq5GevfrSZcydZXbrg7+ZGTwfGKURJkcIvnIPNhnn1\nJjwXLyHw4M7cTQYXvvUYTqd53VAB/mTVD+f4Jn8a35lpv66tEBWRFH5RoRmOHcP45Ze4b9iA4cQJ\nAGwmb9JHjiBt1EPM+M9tvFUtiUaNLOiBe+/NdGzAQpQCKfyiQnL54w88FyzEfcMGdOSePsFStSoL\nXJ7GY3wEfYfmvvW///4CBvktEE5G3vKi4lAKbcMWfBYvwLg7dzonR+fGgRZDqPFEN7K7dCF0t5Hs\nbA3IBpCiL5ySvO1FhaA/dgyf6dNx37IFAJu3N+kjRrCj1aOMndGAX7vHo2nQrl2OgyMVwvGk8Ity\nLWFfHOrtD2gauQTNasXs5cNrajIjtg7ALchEY+CrRhfkYK0Q/yCFX5RLhsOH8Zo/n8rffY/OakFp\nGunDh5P67LP8NbsOJxPTaRCUe8b7qlVtBexNCOcihV+UK/r/28/hBxbR+eJaAJRez/aqA9l/1xP0\ne7keAG++mezIEIUo86TwizIvM1ND2xlFyEfv4L55M5UAs96NnAeHkjZ2LNasMKqlaoDZ0aEKUS5I\n4RdlmuHAAXLGzKbhqZ8BsHl4cG7ASPpunsLaGXoMBqgjFzEU4qYUWPgTExOZP38+ycnJaJpG9+7d\n6dOnD2lpabzzzjskJCQQHBzMhAkTMBqNpRGzqOCSTmdw4c2vaZPwE+5btxIMpGAi6+GHsD31CJq/\nP6syNAwG5ehQhSiXCiz8er2ekSNHEhYWRlZWFpMmTaJp06Zs2bKFxo0b079/f9asWcPq1asZPnx4\nacQsKijDwYMYv/mGSis+p1FqCpB7GcP0kSOZrZvC7Xd70tI/dzrHaJSiL0RR6Qpq4OvrS1hYGADu\n7u5UrVqVxMREoqKi6Ny5MwBdunRh9+7dJRqoqLgMhw9zuNEDBPfsidfixehTUzhauT0fdVpM3Pbt\npEyfzriXXGnZUubwhSgONzXHHx8fz6lTp6hXrx7Jycn4+voCuX8ckpNlJYUoPGVT2PYdInDpu7iv\nW0ewUmTrPbAMGUjG0KGkVG6J8agLtpBsR4cqRIVT6MKflZXF22+/zahRo3B3d7/qfu0635CJjo4m\nOjrafjsiIgKTyVSEUCseV1dX58uFxYLLhx+SPWshPhdyz3+vXF1JHz6a29e+wE+ve+PrC/WAevUA\nXB0ZrUM45fviOiQX+a1cudL+c3h4OOHh4UXaT6EKv9Vq5a233qJTp060bt0ayB3lJyUl2f/38fG5\n5mOvFVxqamqRgq1oTCaT0+TCZoNTy3bR5vMXcDl0CHfgghZIzoB70E0Zh2f9+qycmIZen4qTpOS6\nnOl9URDJRR6TyURERESx7KtQhX/RokWEhobSp08f+7aWLVsSGRnJgAEDiIyMpFWrVsUSkKh4DAcO\nYJz7LqHrfwDAUq0aKS+9xNzj9xFUSRERknsqZF9fOWArRGnQlFI3/G07fPgw06ZNo3r16miahqZp\nDB06lDp16jB37lwuXLhAUFAQEyZMwNPTs1Cdnj17tliCL+8q+mhmw1rovPk1an29AIAcFyMra02k\n2w8jwMMjX9uKnoubIbnII7nIExJSfBf/LLDwlwQp/Lkq6ptad+ECXgsWoPt0FcbMSyidjsxBg4gb\n/xzvranPM8+kXnXStIqai6KQXOSRXOQpzsIv39wVxebEcR1HnvuGh6Ino7v8y/p/hhZkz5pO6JCW\nuAETJ8ovsRCOJoVf3BKbDXQ6MBw5QqtJk+mwexcAWZ07kzppEmdsrahTxwLI/L0QZYUUfnFLpgxN\nYYbLy1TZuhLNaiXDM5D3ar/Jgyt6gqbRXE6cJkSZI4Vf3BSbDS5e1BGcfRrTu+/y8Y4v0dssKL2e\n9BEjuPjMJPq7+oEmI3whyiop/OKmbP/6Euq1hTRO/gAtJwel0/G933CqLR1PcPvq6AFvmdYRokyT\nwi8KdOCAgcaV4zAtXMDgZcvRZWcBkNmvH6n//S8t6tSVSxsKUY5I4Rc3pDv8F+eHfEPXjI9xyU4H\n4FSLu9nUcTK9n6sJgNR8IcoXKfziKn//rSf9ZCJtf3wd44oVjLj8VY/M7t1Je/ZZXBo3preDYxRC\nFJ0UfpGP4cABAmd8RIsd3+KKGWUwkNGvHx+4P06vlxpiMsn8vRDlnRR+QU4OfPNeCo+dfAHP1d8S\nDFjREdOoB8b5U7DUr899gKzFF6JikMLv5HQJCQTOn8/jH36Gu8rC5u5OxvDhRN3xGGlBNWhRX9bh\nC1HRSOF3Ujt/cyHg15/o8MnT6C9dAuBHj4E0/v5p9A1rUw9AvnwlRIUkhd8J6U+f5vZZL1N1948A\nZHfsSPKLL5KR0gKtfo6DoxNClDQp/E4iOxvemunKq24z8FmyCM1qJV1v4tcek2nywYOg09EeKfpC\nOAMp/E7Ca38Uz37xHL4ZR1CaRsbAgRwbNZmqYaGgszk6PCFEKZLCX4Ft3uyGy9Ej9P9tOm6bNqEp\nxVF9fTw+n42uQyuCAJCiL4SzkcJfUZnNtNk4nxqfvIUbOSi9ntSxY4nrP5FajVwcHZ0QwoGk8Fcg\nFgvMnOnNi31/o/LkCYRERwOw87aRhH36NLbgYGo5OEYhhONJ4a9ADNZsum15lZAP3kJvs2CpVo3Y\naW8S0KkTNk/58pUQIpcU/nJuzx4Xjh0zMLzOdnyfeYb7j/4FwIWhozHPmISLpycu8o1bIcQ/SOEv\n5wJ0F4l5fgmB5nloNhuWmjX5a9JcfO9p7ejQhBBllM7RAYib99FHnqSkaLhFRtLmwfaMz56LUpA2\ndizxGzdK0RdC3JCM+Msh3ZZfUQvfJODcVgAyWrYhdsJ0PLs2dXBkQojyoMDCv2jRIvbs2YOPjw9z\n5swBYNWqVfz888/4+PgAMHToUJo1a1aykTqxixd1REW50KvdBbynT2fK5i8BsLgbyXj6SdLGjcNT\nr3dwlEKI8qLAwt+1a1d69+7N/Pnz823v27cvffv2LbHARJ7MTI2Pnozhfp8IPM6cQLm5cWrQY+j/\n+x/0VQIdHZ4QopwpsPA3aNCAhISEq7YrJStFStKxYwb8/a0EpZ2i4UcfsSXtI/SpVswNG3Jp4UJc\n69VzdIhCiHKqyHP869evZ+vWrdSuXZsRI0ZgNBqLMy6n982HOXTcvYDwE3PRsrNRmsZfvR7GtOB5\nlIeHo8MTQpRjRSr8vXr14r777kPTNL788kuWL1/O2LFjr9k2Ojqa6MvfIAWIiIjAZDIVLdoKxtXV\nNV8u4uM1gr0zcVmxgrnfzUCXlASApXdvsidNokqLFo4KtcT9OxfOTHKRR3KR38qVK+0/h4eHEx4e\nXqT9FKnwe3t723/u3r07s2bNum7bawWXmppalG4rHJPJZM9FZqbGpFb7Web9BO5/X/4SVnh7Yh6Z\nTMjgywW/Auftn7lwdpKLPJKLPCaTiYiIiGLZV6EKv1Iq35x+UlISvr6+AOzatYtq1aoVSzDOyGIB\nt9hThLzyCl8n/QRJYA4LI+3ZZ8np358QTXN0iEKICqbAwj9v3jwOHjxIamoqY8eOJSIigujoaGJi\nYtA0jaCgIMaMGVMasVY4679Mw/b6XCJi30XLycFmNLKu8URqLxiFXxU3R4cnhKigNOWA5Tlnz54t\n7S7LFrMZr4UL8Vq8BF1KMgAZgwaRMnkytipVHBycY8hH+jySizySizwhISHFti/55m4p+/TlC4yJ\n/A/eR/4AILHJHbxX5RXGvFvfwZEJIZyFFP7SohTGL7/k6eUzcMtKxVwlBPPiRWS3aoVMlAkhSpMU\n/hJmscD2hccZsG8mHhs2ALCjcn92jZ7Lw10DK/RKHSFE2SSFvwRpmZl4z3qLwe8vxYAV5eZG0qxZ\nVLp7MIPdHR2dEMJZSeEvIeaffqXqK5MwnDqF0jQ+cHucVqv+Q3CLysj3boUQjiSFv7gpRcrkhTT4\n9DUAzA0bkvTmm1QztCGgkdnBwQkhhBT+YqVLTMTn2WcJ2bABGxpf3/YCHdY9DC4uNEaKvhCibJAr\ncBWTXS/9ilf7Hnhs2IDN25vzCz/kl07PYdFcHB2aEELkIyP+W6RlZOD98svc++mnAFxq1p6cpe9A\n1aq80F9W7Aghyh4Z8d+CmG8PEdC9J56ffopydeXX/jN4vP56rFWrOjo0IYS4LhnxF5HHmjW0evoZ\nXK1ZmBs05NK786gdHs4sq4zyhRBlm4z4b5JS4PH11/g++SSu1iy+8RnJ/FGRWC6feloufSuEKOtk\nxH8TLGbFxtvfY3Rc7vUHUp96ikr3TaaqQS5DKYQoP6TwF5ZS+M15g9Fx87GgJ/GlV7A+OpJaWB0d\nmRBC3BQp/AWw2WD36kTu/ngkrnv3ovR63u/1OZWbdeN2chwdnhBC3DQp/AWw/bqbthPG42o9gzUg\ngOQ33qB/nw4gRV8IUU5J4b8B97Vr8Xv6aTRrDlGGNqQu+YT67eTCz0KI8k1W9VzDmdM6fuz8Ef7j\nxqHl5JA+ciTJ339LrVZS9IUQ5Z+M+P/NYiH8vSm0ObYCGxopU6eS8egYGspFz4UQFYSM+C9TCk4f\nzMR/1Cg8V6zA5ubO681XcLjPOJCiL4SoQGTEf1nsrjiMEQ/hbv0/rP7+XFy2jIdatgRZrimEqGCk\n8AOG6GhaPj4CvfU8x/X1yF62HN+W1R0dlhBClAinnuq5eFHHr1N3EjhwIPrz58lu04bd89ZhbCxF\nXwhRcRU44l+0aBF79uzBx8eHOXPmAJCWlsY777xDQkICwcHBTJgwAaPRWOLBFjfPzf9jwEePoCOH\njP79SXr7bTq5y8VwhRAVW4Ej/q5du/LCCy/k27ZmzRoaN27MvHnzCA8PZ/Xq1SUWYElx/+EHwp4Z\njRs5fOg1nr3PLAIp+kIIJ1Bg4W/QoAGenp75tkVFRdG5c2cAunTpwu7du0smuhIQG6vnwzt/xm/s\nWDSLhbSxY+mydzI1a8uJ1oQQzqFIB3eTk5Px9fUFwNfXl+Tk5GINqiTV2rOGFw6OQ1M2ksY9QcaU\nSXjIak0hhBMpllU92g3WuUdHRxMdHW2/HRERgclU+t+AzcoCz6hteEyYgKZsrG46lWqjnqO+t+NG\n+q6urg7JRVkkucgjucgjuchv5cqV9p/Dw8MJv3wdkJtVpMLv6+tLUlKS/X8fH5/rtr1WcKmppXuV\nqpQUjemdo1mRPAgtO5P04cO5fdajoKVQyqHkYzKZSj0XZZXkIo/kIo/kIo/JZCIiIqJY9lWo5ZxK\nKZTKGxm3bNmSyMhIACIjI2nVqlWxBFNSAo/8zqcpA9BnZ3LmzmEkv/66fBtXCOG0Chzxz5s3j4MH\nD5KamsrYsWOJiIhgwIABzJ07ly1bthAUFMSECRNKI9abFheno/qBjfg/8ghadjbH295H7OQ51NHL\ngVwhhPPS1D+H8qXk7NmzJd6HzQYv3bGfpWfuwWAzkz58OMkzZ4KLS4n3XVjyMTaP5CKP5CKP5CJP\nSEhIse2rwn5z123fHpZejMBgM7Oj1aMkz5pVpoq+EEI4SoUr/KdP67Hu3k/AkCEY0lJI7tkXy5vT\nZE5fCCEuq3CF/8f3zuN534Po0tPJ7NeP9PfnU6eezOkLIcQVFersnLoLF5i6bQAGSwI7fXpQ7Z15\naDK9I4QQ+VSIEb/ZDGf+ysZ/5EgMp2LIadwY21eL0dxcHR2aEEKUORVixP/HbgOuw8fjmrMPS7Vq\nXPzkE8KCPRwdlhBClEkVovD33DIDU85qkjUfDr/8OdWCgx0dkhBClFnleqonPl6HccUKTAsXogwG\nzr33AaE9azk6LCGEKNPK7YjfYoG5d+/jo3OTAUh+4w28723v4KiEEKLsK7cjfvcTf/FRyv3olZUf\nG08gY+hQR4ckhBDlQrkr/BYLWGIT8B8xAn1aCum97ybo/WcdHZYQQpQb5a7wr/zMwKVuYzGcPk1O\n8+YkvzePqtXkC1pCCFFY5W6Of8zpGXin/cY5fVXOz1xGJQ9ZtimEEDejXI34XbduxXvxQpReT8ys\nJQTeFujokIQQotwpF4VfKXj5sUw8xz4FQOqECdQY2hy93sGBCSFEOVQupno0ZWPG8ZF4JMVzrFpH\njE8+6eiQhBCi3CoXI37Pjz+m0sHtmP0D2fTwB8hQXwghiq5MF36rFSKXnsH02msApL45i34Pezk4\nKiGEKN/K9FRPUqLitjlPosvKIqH3fZjvusvRIQkhRLlXpkf81X5YRpP0XSSbQng/fI6jwxFCiAqh\nTI74lQJL3EW85+QWe+vclxnVWy6oIoQQxaFMjvijolzY0mURuqQksjt0IEumeIQQotiUyRF/O+8D\nBKYvwYqOFa3e4C65ULoQQhSbWyr8jz/+OEajEU3T0Ov1vP7667cekVJ4T5+OzmYl+YGRdBxXC5Bz\n8QghRHG5pcKvaRrTpk3Dy6t4llieP6/j9ILN9N+6FZuPD5mTJuLpKUVfCCGK0y3N8SulUKr4CnP6\nxRxarHgJgP33PYfN37/Y9i2EECLXLY/4X331VXQ6Hd27d6dHjx63FEyTX97HJ/s4SVXqEXPXQwTJ\nFI8QQhQ7Td3CkP3SpUv4+fmRkpLCK6+8wujRo2nQoEG+NtHR0URHR9tvR0REkJqaevXO4uLxatEc\nLTWVjG+/xXqLf0TKA1dXV3JychwdRpkgucgjucgjuchjMplYuXKl/XZ4eDjh4eFF2tctFf5/WrVq\nFR4eHvTt27fAtmfPns13+/RpPcd7TWNY8hIyunUn6dNPiiOkMs9kMl37j6ATklzkkVzkkVzkCQkJ\nKbZ9FXmOPzs7m6ysLACysrL4888/qVatWpH2Vd0zgYiM5QC84DKrqCEJIYQohCLP8ScnJ/Pmm2+i\naRpWq5WOHTvStGnTIu3L8/PPMZizyOralfELQpDlm0IIUXKKXPiDg4N58803bzmA44dttFu2DID0\n0aPx8JCiL4QQJcnhp2z4eXwkhnPnSPCvS/LtnR0djhBCVHgOL/xTPOcBsLHhWFzdHR6OEEJUeA49\nV4/Ln3/iGhWFzdubrh/3RUndF0KIEuewUnvypJ6EN74CIOP++1Geno4KRQghnIrDCn9qQg7Vtq0B\n4CvjSEeFIYQQTsdhhb/N+R8wWZPIaNiYRvfXcVQYQgjhdBxW+I3ffguAeehgatSwOioMIYRwOg4p\n/F9/CvrIbQBkFuIUD0IIIYqPQwr/bQlbcTFnsk/fkl+PFe00D0IIIYrGIcs5b0/4CYBaT3YmoK2c\neU8IIUqTQ0b8bps2AZDdozt6vSMiEEII5+WQwm84e5YsUwCXahXtpG5CCCGKzmGrev5waUvM3y6O\n6l4IIZyWw07Z0PiR20i7zeKo7oUQwmk5bMSf07Klo7oWQgin5pDCb9N0xAS2cETXQgjh9BxS+OO8\na3MiwccRXQshhNNzSOH3bVebDh1k/b4QQjiCQwq/tWZNR3QrhBACBxX+nXG1HdGtEEIIHFT4Y11l\nxC+EEI7ikMLf6zkp/EII4Si39AWuffv2sWzZMpRSdO3alQEDBhTqcbZKlW6lWyGEELegyCN+m83G\nhx9+yAsvvMBbb73F9u3biY2NLc7YhBBClIAiF/5jx45RpUoVgoKCMBgM3HHHHezevbs4YxNCCFEC\nilz4L168SEBAgP22v78/Fy9eLJaghBBClByHnatHCCGEYxT54K6/vz8XLlyw37548SL+/v5XtYuO\njiY6Otp+OyIigpCQkKJ2W+GYTCZHh1BmSC7ySC7ySC7yrFy50v5zeHg44eHhRduRKiKr1arGjx+v\n4uPjldlsVhMnTlSnT58u8HFfffVVUbuscCQXeSQXeSQXeSQXeYozF0Ue8et0OkaPHs2rr76KUopu\n3boRGhpa1N0JIYQoJbe0jr9Zs2bMmzevuGIRQghRCkr94G6R56QqIMlFHslFHslFHslFnuLMhaaU\nUsW2NyGEEGWeLOcUQggnI4VfCCGczC0d3L0ZRT2hW3mVmJjI/PnzSU5ORtM0unfvTp8+fUhLS+Od\nd94hISGB4OBgJkyYgNFoBGD16tVs2bIFvV7PqFGjaNq0qYOfRfGy2WxMnjwZf39/Jk2a5LS5yMjI\nYPHixZw+fRpN0xg7dixVqlRxylysW7eOLVu2oGka1atXZ9y4cWRlZTlFLhYtWsSePXvw8fFhzpw5\nAEX6nThx4gQLFy7EbDbTvHlzRo0aVXDnxbYw9Aauteb/zJkzpdG1w1y6dEmdPHlSKaVUZmamevLJ\nJ9WZM2fUp59+qtasWaOUUmr16tXqs88+U0opdfr0afXss88qi8Wi4uLi1Pjx45XNZnNU+CXi+++/\nV/PmzVNvvPGGUko5bS7mz5+vNm/erJRSymKxqPT0dKfMRWJionr88ceV2WxWSin19ttvqy1btjhN\nLg4dOqROnjypnnnmGfu2ojz3yZMnq6NHjyqllHrttdfU3r17C+y7VKZ6nPGEbr6+voSFhQHg7u5O\n1apVSUxMJCoqis6dOwPQpUsXex6ioqJo3749er2e4OBgqlSpwrFjxxwVfrFLTExk7969dO/e3b7N\nGXORkZHB4cOH6dq1KwB6vR6j0eiUuYDcT4FZWVlYrVZycnLw9/d3mlw0aNAAT0/PfNtu9rknJSWR\nmZlJnTp1AOjUqVOhamupTPVc64Ru5fkFu1nx8fGcOnWKevXqkZycjK+vL5D7xyE5ORnIzVG9evXs\nj6loJ71bvnw5Dz74IBkZGfZtzpiL+Ph4TCYTCxcu5NSpU9SqVYtRo0Y5ZS78/f3p27cv48aNw83N\njSZNmtCkSROnzMUVN/vc9Xp9vtoaEBBQqJzIwd0SlpWVxdtvv82oUaNwd3e/6n5N0xwQVem6Mo8Z\nFhaGusHqYWfIhc1m4+TJk/Tq1YtZs2bh5ubGmjVrrmrnDLlIT08nKiqKhQsXsmTJErKzs/n111+v\naucMubieknrupTLiL+wJ3Soaq9XKW2+9RadOnWjdujWQ+1c8KSnJ/r+Pjw9wdY4SExMrTI4OHz5M\nVFQUe/fuJScnh8zMTN577z2nzIW/vz8BAQHUrl0bgLZt27JmzRqnzMX+/fsJDg7Gy8sLgDZt2nDk\nyBGnzMUVN/vc/f39SUxMvGp7QUplxF+nTh3Onz9PQkICFouF7du306pVq9Lo2qEWLVpEaGgoffr0\nsW9r2bIlkZGRAERGRtrz0KpVK3bs2IHFYiE+Pp7z58/b5+3Ku2HDhrFo0SLmz5/P008/zW233cYT\nTzzhlLnw9fUlICCAs2fPArnFLzQ01ClzERgYyNGjR8nJyUEp5ZS5UErl+xR8s8/d19cXo9HIsWPH\nUEqxdetW+yDzRkrtm7v79u3j448/tp/QraIv5zx8+DDTpk2jevXqaJqGpmkMHTqUOnXqMHfuXC5c\nuEBQUBATJkywH+BZvXo1mzdvxmAwlPulatdz8OBBvv/+e/tyTmfMRUxMDEuWLMFisVCpUiXGjRuH\nzWZzylysWrWKHTt2oNfrCQsL47HHHiMrK8spcjFv3jwOHjxIamoqPj4+RERE0Lp165t+7idOnGDB\nggX25ZwPPfRQgX3LKRuEEMLJyMFdIYRwMlL4hRDCyUjhF0IIJyOFXwghnIwUfiGEcDJS+IUQwslI\n4RdCCCcjhV8IIZzM/wNEXVBhcjIITAAAAABJRU5ErkJggg==\n", "text/plain": [""]}, "metadata": {}, "output_type": "display_data"}], "source": ["import matplotlib.pyplot as plt\n", "# let's have a quick look at a few random walks\n", "plt.plot(walks[:10,:].transpose())\n", "plt.title('A few random walks')\n", "# Let's see how the root mean square of the position evolves with time/nb of steps\n", "rms_position = np.sqrt( (walks**2).mean(axis=0) )\n", "plt.figure()\n", "t = 1 + np.arange(len(rms_position))\n", "plt.plot(t, np.sqrt(t), ':b', lw=3) #Just to show the fit \n", "plt.plot(t, rms_position, '-r', lw=2)\n", "plt.title('Root Mean Square of Position by Time')\n", "# What are the highest/lowest positions \n", "print('Highest position:{max}\\tLowest position:{min}'.format(max=walks.max(), min=walks.min()))\n", "# How many walks wander further than 50?\n", "bound = 50\n", "hits_the_bound = np.any(np.abs(walks)> bound, axis=1) #for each walk, do we go further than the bound at any time?\n", "print('Number of walks over bound(={}):{}'.format(bound, hits_the_bound.sum()))\n", "# Among the walks that go beyond the bound, what is the mean of the first hits?\n", "# we use argmax on the boolean array to get the first True value\n", "first_hits = (np.abs(walks[hits_the_bound,:])>bound).argmax(axis=1)\n", "print('Mean crossing time:{}'.format(first_hits.mean()))"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Exercice 3 : retrouver la s\u00e9rie al\u00e9atoire \u00e0 partir des marches al\u00e9atoires\n", "\n", "Dans cet exercice, on cherche \u00e0 retrouver la s\u00e9rie initiale \u00e0 partir de la somme cumul\u00e9e de celle-ci. On veut calculer en quelque sort sa d\u00e9riv\u00e9e."]}, {"cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [{"data": {"text/plain": ["array([[ 1, 1, 1, 1, 1, -1, -1, 1, -1, 1],\n", " [-1, 1, 1, -1, -1, 1, 1, -1, 1, -1],\n", " [-1, -1, -1, -1, -1, -1, -1, -1, -1, -1],\n", " [-1, 1, 1, 1, -1, 1, 1, 1, 1, 1],\n", " [-1, 1, 1, -1, 1, -1, 1, -1, -1, -1],\n", " [ 1, 1, 1, 1, 1, -1, 1, -1, -1, 1],\n", " [ 1, 1, -1, 1, 1, 1, 1, 1, -1, 1],\n", " [-1, 1, 1, 1, 1, -1, 1, 1, -1, -1],\n", " [-1, -1, 1, -1, -1, -1, -1, 1, -1, 1],\n", " [-1, -1, 1, 1, 1, -1, 1, 1, 1, 1]], dtype=int32)"]}, "execution_count": 10, "metadata": {}, "output_type": "execute_result"}], "source": ["derivee = walks[:,1:] - walks[:,:-1]\n", "derivee[:10,:10]"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Exercice 4 : simulation, r\u00e9gression, estimation par maximisation de la vraisemblance"]}, {"cell_type": "markdown", "metadata": {}, "source": ["* On commence par simuler la variable $Y = 3 X_1 -2 X_2 +2 + \\epsilon$ o\u00f9 $X_1,X_2,\\epsilon \\sim \\mathcal{N}(0,1)$ \n", "* On souhaite ensuite retrouver les coefficients dans la [r\u00e9gression lin\u00e9aire](http://fr.wikipedia.org/wiki/R%C3%A9gression_lin%C3%A9aire) de $Y$ sur $X_1$ et $X_2$ dans un mod\u00e8le avec constante, par la m\u00e9thode des Moindres Carr\u00e9s Ordinaires. On rappelle que la forme matricielle de l'estimateur des MCO est $\\hat{\\beta} = (X'X)^{-1}X'Y$\n", "* Enfin, $Y$ \u00e9tant normale, on souhaite estimer ses param\u00e8tres par maximisation de vraisemblance:\n", " * La densit\u00e9 s'\u00e9crit: $f(x, \\mu, \\sigma) = \\frac{1}{\\sigma \\sqrt{2\\pi} } e^{ -\\frac{(x-\\mu)^2}{2\\sigma^2} }$\n", " * La log-vraisemblance: $\\ln\\mathcal{L}(\\mu,\\sigma^2) = \\sum_{i=1}^n \\ln f(x_i;\\,\\mu,\\sigma^2) = -\\frac{n}{2}\\ln(2\\pi) - \\frac{n}{2}\\ln\\sigma^2 - \\frac{1}{2\\sigma^2}\\sum_{i=1}^n (x_i-\\mu)^2$ ou encore en divisant par $n$ : $-\\frac{1}{2}\\ln(2\\pi) - \\frac{1}{2}\\ln\\sigma^2 - \\frac{1}{2n\\sigma^2}\\sum_{i=1}^n (x_i-\\mu)^2$\n", " * L'\u00e9criture des conditions au premier ordre donne une formule ferm\u00e9e pour les estimateurs du maximum de vraisemblance: $\\hat{\\mu} = \\overline{x} \\equiv \\frac{1}{n}\\sum_{i=1}^n x_i$, $\\hat{\\sigma}^2 = \\frac{1}{n} \\sum_{i=1}^n (x_i - \\overline{x})^2$.\n", " * V\u00e9rifiez en les impl\u00e9mentant directement que vous trouvez bien la m\u00eame solution que le minimum obtenu en utilisant *scipy.optimize.minimize* pour minimiser l'oppos\u00e9 de la log-vraissemblance.\n", " \n", "**version matricielle**"]}, {"cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["coef X1, coef X2, constante\n"]}, {"data": {"text/plain": ["array([ 2.99395809, -2.00169881, 2.00215411])"]}, "execution_count": 11, "metadata": {}, "output_type": "execute_result"}], "source": ["import numpy as np\n", "import math\n", "from scipy.optimize import minimize\n", "n_samples = 100000\n", "x1_x2_eps = np.random.randn(n_samples,3)\n", "y = 3*x1_x2_eps[:,0] - 2*x1_x2_eps[:,1] + 2 + x1_x2_eps[:,2]\n", "\n", "X = np.hstack( (x1_x2_eps[:,:2], np.ones((n_samples,1))) )\n", "beta_hat = ( np.linalg.inv((X.T).dot(X)) ).dot( (X.T).dot(y) )\n", "print(\"coef X1, coef X2, constante\")\n", "beta_hat"]}, {"cell_type": "markdown", "metadata": {}, "source": ["**version scipy**"]}, {"cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [{"data": {"text/plain": [" status: 0\n", " message: 'Optimization terminated successfully.'\n", " success: True\n", " nit: 30\n", " fun: 1.8165031271016212\n", " nfev: 57\n", " x: array([ 1.99871699, 13.91551151])"]}, "execution_count": 12, "metadata": {}, "output_type": "execute_result"}], "source": ["def log_likelihood(mu,sigma_square, x):\n", " return - 0.5 * math.log(sigma_square) - sum((x - mu)**2)/(2*sigma_square) / len(x)\n", "\n", "def neg_log_likelihood_vectorielle(theta):\n", " return -log_likelihood(theta[0],theta[1],y)\n", "\n", "theta0 = np.array([2., 14])\n", "optim_res = minimize(neg_log_likelihood_vectorielle, theta0, method='Nelder-Mead')\n", "optim_res"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Est-ce bien le r\u00e9sultat attendu :\n", "\n", "* $\\mathbb{E}Y = 3\\mathbb{E}X_1 - 2\\mathbb{E}X_2 + 2 + \\mathbb{E}\\epsilon = 2$\n", "* $\\mathbb{V}Y = 9\\mathbb{V}X_1 + 4\\mathbb{V}X_2 + \\mathbb{V}\\epsilon = 14$\n", "\n", "Toutes les variables sont ind\u00e9pendantes. On v\u00e9rifie que cela correspond aux r\u00e9ponses cherch\u00e9es :"]}, {"cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [{"data": {"text/plain": ["(1.9987233148817718, 13.91554073655678)"]}, "execution_count": 13, "metadata": {}, "output_type": "execute_result"}], "source": ["np.mean(y), np.std(y)**2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Exercice 5 : Optimisation quadratique (sous contraintes) avec cvxopt\n", "\n", "voir [correction](http://www.xavierdupre.fr/app/ensae_teaching_cs/helpsphinx/notebooks/td1a_correction_session9.html#td1acorrectionsession9rst)"]}, {"cell_type": "code", "execution_count": 13, "metadata": {"collapsed": true}, "outputs": [], "source": []}], "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}, "language_info": {"codemirror_mode": {"name": "ipython", "version": 3}, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.1"}}, "nbformat": 4, "nbformat_minor": 2}