{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Assignment 9: Improving performance\n", "\n", "This homework provides practice in making Python code faster. Note that we start with functions that already use idiomatic `numpy` (which are about two orders of magnitude faster than the pure Python versions)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "%load_ext Cython" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import math\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "from sklearn.datasets import make_blobs \n", "from numba import jit, vectorize, float64, int64" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "sns.set_context('notebook', font_scale=1.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Functions to optimize**" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def logistic(x):\n", " \"\"\"Logistic function.\"\"\"\n", " return np.exp(x)/(1 + np.exp(x))\n", "\n", "def gd(X, y, beta, alpha, niter):\n", " \"\"\"Gradient descent algorihtm.\"\"\"\n", " n, p = X.shape\n", " Xt = X.T\n", " for i in range(niter):\n", " y_pred = logistic(X @ beta)\n", " epsilon = y - y_pred\n", " grad = Xt @ epsilon / n\n", " beta += alpha * grad\n", " return beta" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEECAYAAADK0VhyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4XGXd//H3N5nse5u26Z62dEuhZSkgshQolH2xKPgosoj6Ax8BBcQNpQX1EZAWBRQVEFGwIJuClm7QFpC1lEL3fQ1d0ibNnsnM3L8/ZooxpM1kPTOZz+u65prkPudMvqedzCfn3Oe+jznnEBGRxJbkdQEiIuI9hYGIiCgMREREYSAiIigMREQEhYGIiKAwEBERFAYiIoLCQEREAJ/XBUSrsLDQFRcXe12GiEjcWLJkSZlzrk8068ZNGBQXF/Pee+95XYaISNwwsy3RrqvTRCIiojAQEZEow8DMDjOz35nZh2YWNLOFUW6XZ2Z/NLNyM9tvZk+YWe8OVSwiIp0u2j6DccC5wFtAShte/2lgFPA1IATcBbwAnNyG1xARkS4WbRi86Jz7O4CZPQMUtraBmZ0ATAEmOecWR9p2AG+b2RnOufntrFlERDpZVKeJnHOhdrz2OcCuA0EQeZ13gE2RZSIiEiO68tLSMcDqFtpXRZaJiEhEMORoDIYiD0cgGMIfDBEIOob2zsTMuvTnd2UYFAAVLbSXA8OjeQEzmwbcDtC/f/9OK0xEpC2cc9T4g1TWNVLdEKCqPkB1Q4Dq+gA1/gC1DQFq/EFq/QFq/UHq/EHqGoPUNwapbwxR3xikIRB+9gdCNHzyCNIYDOEPhAgd4g7Ea356Nmm+5C7dx5gedOacmwZMA5g4caJu1iwinaIxGGJPVQO7qxooq2pgb00DZdV+9lb7qaj1s6/WT3mNn/11jVTUNVJZ13jID+topPqSSPMlkZ6STGpyErnpPlJ9qaT6kkhNTsKXbKT6wstSkg1fchIpSYYvuWuPCA7oyjAoB1oaBl0QWSYi0iX21zayZV8N2/bVsb28lh0VdZRW1LOzso6d++spq/a3+hqpyUnkZ6bQJzuNw/pkk5uRQk66L/JIITvNR1ZqMllpPrLTfGREvs5ISSYjNTn8nJJMekoyab4kkpK650O9vboyDFbT8iWkYwhfXioi0m6hkGNbeS1rdlaxfk81G3bXsGFPNZv31lBR29jiNukpSQzIy+Cwvtn0zUmnb04afXLSKMxOo3d2Kr2z0ijISqFXVioZKcldfp4+lnRlGMwGfmxmJznnXgcws4mE+wtmd+HPFZEepjEYYs3OKj7asZ/lO/azvLSStTurqGsM/td6KcnG4F6ZHDU4n6G9sxjcK5NBBRkMKshgYH4GeRkpCfUB3xZRhYGZZRIedAYwEMg1s89Hvv+Xc67WzNYDi5xz1wA45940s7nA42Z2C/8ZdPa6xhiIyKHsr2vkvc37eGfTPt7fWs6H2/fTEPjPFe4pycaIPtmMKcphVFEOo/rmMKJvNoMLMvAla5ad9oj2yKAv8LdmbQe+HwZsjrxW8+7uy4CZwKOExzS8BNzQnkJFpOdqCARZsrmcxevKeH39HlaUVuIiHbZJBqOLcjlycD4TBuVx+MA8RvbL7vKraxJNVGHgnNsMHPLYyjlX3EJbBXB15CEi8om91Q0sWL2beSt38fq6sk9O+aQmJ3FscS8+M6wXxw3rzZFD8slOi+kLH3sE/QuLSLepqPXz8vKd/GNZKW9t3PvJ5ZrD+2Rx6qi+nDKqkOOH9SYjVX/1dzeFgYh0qUAwxMI1e3j6vW28sno3gUgCHD0kn7PGFXFmST+G98n2uEpRGIhIl9i5v56/vLWFp97bxp6qBgDGFOVw8VEDOe+I/gzulelxhdKUwkBEOtXSreU8/Pom5izfSSDkyE33ceUJQ/nCxMEcPjDP6/LkIBQGItJhzjneWL+XB19dz5sb9wLho4CrTyzmwgkD1QcQBxQGItIh/95Qxj1z1rB0a3heylNG9eHaScM5YXhvDfCKIwoDEWmX5Tv2c9fLq3ltXRkAZ43rx7dOG8kRg3QqKB4pDESkTcqqG/jlnDU89d42nIOTDivku2eNZsLgfK9Lkw5QGIhIVIIhx+NvbmbGvLVU1QcY1S+bH59fwskjW5qcWOKNwkBEWrV2VxW3PvMhH2yrIDfdx/QLx/Hl44doHqAeRGEgIgcVCIb4zcIN3P/KOhqDjgsnDOAnF5RQmJ3mdWnSyRQGItKirXtrufGppSzdWkFRbjo/+9zhTB7bz+uypIsoDETkU55fup0fv7CC6oYAF04YwJ0XH05eRorXZUkXUhiIyCcaAkGm/WMlf31nK9lpPmZeNoGLjxyo8QIJQGEgIgCUVtRx3RPvs2xbBWP75/LQ5UcztHeW12VJN1EYiAhLtuzjG48vYW+Nn6lHDeRnnztCU0gkGIWBSIJ76cNSbnp6GcGQY/qF47jihKE6LZSAFAYiCco5x28XbeDul9eQnebjwSuOZtIoDSBLVAoDkQQUCjnueGklj/17M/3z0nn0qmMZ2z/X67LEQwoDkQQTCIb4/nMf8cyS7Yzql82frzmefrnpXpclHlMYiCQQfyDEjbOWMnv5TsYPyuNPVx9HQVaq12VJDFAYiCSIxmCIbz35PnNX7uL4Yb14+MqJ5KRrIJmEKQxEEkAgGOLbT33A3JW7+OyI3jxy5bG6dFT+i6YcFOnhgiHHLX9bxj8//JjjisNHBAoCaU5hINKDOef4yd+X88IHpRw9JJ9Hrz6WzFSdEJBPUxiI9GD3v7KeJ97eytj+uTz21ePITlMQSMsUBiI91Kx3tjJj3loG5mfwp6uPJVedxXIIUYWBmZWY2QIzqzWzUjO7w8xaPeloZhPNbK6Z7Ys85pvZ8R0vW0QO5ZXVu/jh8x9RkJnC49ccR1+NI5BWtBoGZlYAzAcccBFwB3AzML2V7QZHtvMBX4k8fMA8MxvasbJF5GDW7Kzi+ieXkpKcxKNXHcuIPtlelyRxIJoTiNcCGcBU51wl4Q/zXGCamd0daWvJeUAO8Dnn3H4AM/s3UAacC/y2w9WLyH/ZV+Pna4+/S40/yANfOoqjhhR4XZLEiWhOE50DzGn2oT+LcEBMOsR2KUAAqGnSVh1p05SIIp3MHwhx7V+WsG1fHTdOHsn54wd4XZLEkWjCYAywummDc24rUBtZdjDPRta518z6mllfYCZQDvytfeWKyMFMf3EF72zax3lH9OfGySO9LkfiTDRhUABUtNBeHlnWIudcKXAacAmwK/KYCpzlnNsTTXFmNs3MnJm50tLSaDYRSUjPLtn+ySWkv/zCBJKSdPAtbdNll5aaWX/CRwBLCJ9qOify9T/NbEg0r+Gcm+acM+ecDRigQ16Rlqz6uJIfvfAROek+fvvlozW6WNolmg7kciCvhfaCyLKD+S7hfoPPO+caAczsFWAdcAtwQ9tKFZHm9tc1ct1fllDfGOLXXzyK4kLds1jaJ5ojg9U06xuIXDaaSbO+hGbGACsOBAGAc84PrABGtL1UEWnKOcf3n/2QzXtrue7UEUwZV+R1SRLHogmD2cBZZpbTpO0yoA5YdIjttgCHm9knk6WbWRpwOLC57aWKSFOz3t3G7OU7Oa64FzefOcrrciTORRMGDwENwHNmdoaZfQOYBsxoermpma03s0eabPcwMAB43szOM7PzgReA/sDvO2sHRBLR+t1VTH9xBbnpPmZ+8Uh8yZpZRjqm1XeQc64cmAwkAy8SHnk8E7i92aq+yDoHtlsCnE144NmfgccJn1o60zm3rDOKF0lE9Y1Brv/rB9Q3hrjrkvEMzM/wuiTpAaKawtA5txI4vZV1iltoWwAsaFdlItKie+asYdXHlfzPcYM554j+XpcjPYSOLUXiyFsb9/LI65sYXpjFj88v8boc6UEUBiJxorohwHefWUaSwb2XTtBNaqRTKQxE4sTP/7WKbfvquO7UEZqATjqdwkAkDixau4cn397KmKIcbtC8Q9IFFAYiMa66IcAPnv0QX5Jx76UTSPNpugnpfAoDkRh3z8urKd1fzzdPO4xxA1qaGUak4xQGIjFsyZZ9PP7WFkb0yeJ/T9MsLtJ1FAYiMaohEOR7z36Ec3DXJeN1eki6lMJAJEb9duEG1u+u5iufGcrE4l5elyM9nMJAJAZtKqvhN69uoCg3nVvPHu11OZIAFAYiMcY5x0/+vhx/MMTtF5SQk57idUmSABQGIjHmXx/t5LV1ZZwyqg9nH657FEj3UBiIxJDqhgB3vrSSVF8Sd1w4DjPdy1i6h8JAJIb8av5adlbWc+2kEbqFpXQrhYFIjFi/u4o/vrGZIb0y+eapGlMg3UthIBIDnHNMf3ElgZDjx+eXkJ6iMQXSvRQGIjFgwardvLaujJNHFnLG2L5elyMJSGEg4rGGQJA7/7mS5CTjJ+eXqNNYPKEwEPHYo69vZsveWq44YSgj++V4XY4kKIWBiIf2VDXwwCvr6JWVyrfPGOV1OZLAFAYiHpoxby01/iDfOXMUeRkaaSzeURiIeGTNziqeencrh/XN5n+OHex1OZLgFAYiHvnZv1YRcvDDc8fgS9avonhL70ARDyxau4fFa/dw4mG9OW20LiUV7ykMRLpZMOT4+T9XYQY/OleXkkpsUBiIdLNn39/Oml1VfP7oQZQMyPW6HBFAYSDSreobg8yct5Y0XxI3TdGlpBI7ogoDMysxswVmVmtmpWZ2h5lFNXmKmU01s3fNrM7M9prZy2am6RglIf3xjc18vL+eq08cRv+8DK/LEflEq2FgZgXAfMABFwF3ADcD06PY9mvAk8Bs4Bzga8A6wNf+kkXiU0Wtn98sXE9+ZgrXaVZSiTHRfChfC2QAU51zlcA8M8sFppnZ3ZG2TzGzQmAmcL1z7g9NFj3f0aJF4tGDr66nqj7AbeeN1QAziTnRnCY6B5jT7EN/FuGAmHSI7S6NPP+pnbWJ9Bg7Kur405tbGJifweWfGep1OSKfEk0YjAFWN21wzm0FaiPLDuZ4YA1wjZltN7NGM3vbzD7b7mpF4tSv5q/FHwhx05mjdK8CiUnRhEEBUNFCe3lk2cEUAaOB24DvARcANcDLZtYvmuLMbJqZOTNzpaWl0WwiEnPW767mmSXbGdUvm4uPGuh1OSIt6spLSw3IBq5xzj3hnHsZuBgIAt+K5gWcc9Occ+acswEDBnRhqSJd5965awg5uHnKaJKTNMBMYlM0YVAO5LXQXhBZdqjtHLDwQEOk32EJUBJ9iSLx68PtFcxevpMjB+czpSSqA2IRT0QTBqtp1jdgZoOBTJr1JTSzivDRQfM/hQwItaFGkbh1z5w1ANx69mhNOyExLZowmA2cZWZNb8F0GVAHLDrEdi9Fnk870GBmecAxwLI21ikSd/69oeyT+xp/dkSh1+WIHFI0YfAQ0AA8Z2ZnmNk3gGnAjKaXm5rZejN75MD3zrn3gL8Dj5jZlWZ2HvAPoBF4sBP3QSTmOOe4d+5aAG6ZMtrjakRa12oYOOfKgclAMvAi4ZHHM4Hbm63qi6zT1OXAC8AM4BnCQXB65DVFeqyFa/awZEs5U0r6MWFwvtfliLQqqmkhnHMrgdNbWae4hbZq4LrIQyQhhEKOX85dgxmajE7ihmYtFelkc1bsZEVpJReMH8CYIk1RLfFBYSDSiYIhx73z1pKcZHznTB0VSPxQGIh0on8s28H63dV8/uhBDCvUTO0SPxQGIp2kMRjivvnrSEk2rp98mNfliLSJwkCkkzz3/na27K3li8cOYVBBptfliLSJwkCkE/gDIX69YD2pviT+9zQdFUj8URiIdIKn3tvGjoo6Lj9+KEV56V6XI9JmCgORDqpvDPLAK+vISEnW7SwlbikMRDroibe3squygSs+O5Q+OWlelyPSLgoDkQ6o9Qf47cL1ZKUmc+0pOiqQ+KUwEOmAP7+5hbJqP189aRgFWalelyPSbgoDkXaqbgjw0KIN5KT7+NpJw70uR6RDFAYi7fTYG5sor23k6ycPJy8zxetyRDpEYSDSDpX1jfx+8UbyM1O4+sRir8sR6TCFgUg7PPLaJirrA/y/U0aQk66jAol/CgORNiqv8fPI65sozE7lys8O9bockU6hMBBpo9+/tpHqhgDXnXoYmalR3R9KJOYpDETaoKy6gcfe2Ey/3DS+fPwQr8sR6TQKA5E2+O3CDdQ1BvnWaYeRntL8lt8i8UthIBKlXZX1/OWtLQzMz+DSYwd7XY5Ip1IYiETpgVfW0xAIcf3ph5Hm01GB9CwKA5EobNtXy6x3t1LcO5NLjhnkdTkinU5hIBKFXy9YR2PQ8e0zRpGSrF8b6Xn0rhZpxYY91Tz7/nZG9cvmggkDvC5HpEsoDERacd/8dYQc3HTmaJKTzOtyRLqEwkDkEFZ9XMmLy0o5YmAeZ43r53U5Il0mqjAwsxIzW2BmtWZWamZ3mFnUl1OYWZKZvWdmzszOb3+5It3r3rlrALh5yijMdFQgPVerY+nNrACYD6wELgJGAPcSDpLbovw5XwN0CYbElSVb9jF/1W6OK+7FpFF9vC5HpEtFc2RwLZABTHXOzXPOPQRMB24ys9zWNo6Eyc+AH3WoUpFu5JzjrpfDRwW3nj1aRwXS40UTBucAc5xzlU3aZhEOiElRbH8n8AawoO3liXhj0do9vLNpH5PH9GVicS+vyxHpctGEwRhgddMG59xWoDay7KDMbDzwVeCW9hYo0t1CIcc9c9ZgBrecNdrrckS6RTRhUABUtNBeHll2KPcDDzjn1re1MAAzmxbpdHalpaXteQmRNvvX8o9ZUVrJhRMGMLZ/q2dCRXqELru01My+CIwGftre13DOTXPOmXPOBgzQYB/peo3BEL+cswZfknHTmaO8Lkek20QTBuVAXgvtBZFln2JmKcA9wF1AkpnlAwf+xMoys5x21CrS5Wa9s5XNe2v50vFDGNo7y+tyRLpNNGGwmmZ9A2Y2GMikWV9CE1mELyWdQTgwyoFlkWWzgKXtKVakK9U0BPjVgnVkpSZzw+SRXpcj0q2iuWffbOC7ZpbjnKuKtF0G1AGLDrJNNXBas7Yi4K/AD4FX2lGrSJf6w2sbKav28+0zRlKYneZ1OSLdKpoweAi4AXjOzO4ChgPTgBlNLzc1s/XAIufcNc65ALCw6YuYWXHky4+cc293uHKRTrSnqoE/LN5IYXYqXz95uNfliHS7Vk8TOefKgclAMvAi4QFnM4Hbm63qi6wjEnfuf2UdNf4gN04eSVaabnIviSeqd71zbiVweivrFLeyfDOgYZwSc9bvruaJt7cyrDCLLx6nm9xLYtKspZLwfjF7FcGQ4/vnjNGNayRh6Z0vCe3fG8rCk9EN68WUEk1RLYlLYSAJKxhy/PSlVQDcdt5YTUYnCU1hIAnrufe3s/LjSqYeNZDxg/K9LkfEUwoDSUg1DQF+OXcNab4kTUYngsJAEtSDr65nV2UD/2/SCAbkZ3hdjojnFAaScLbsreHh1zYxIC+d6yaN8LockZigMJCE89N/rsIfDPGDc8eSkapxkiKgMJAE89q6PcxbuYvjintx/vj+XpcjEjMUBpIw/IEQ019ciRn85IISXUoq0oTCQBLGo29sYv3uar503BAOH9jSLTpEEpfCQBJCaUUdv5q/jt5Zqdx61iFv3S2SkBQGkhDueHEldY1Bvn/OGPIyU7wuRyTmKAykx3t1zW5eXrGTY4sLuOToQV6XIxKTFAbSo9X5g9z+9xUkJxl3Xnw4SUnqNBZpicJAerT7Fqxl675arjlpGGOKcr0uRyRmKQykx1q+Yz8Pv7aJIb0y+c4Zo7wuRySmKQykRwoEQ3zv2Q8Jhhw//9wRGmks0gqFgfRIj7y+iRWllXz+mEGcNLLQ63JEYp7CQHqcjXuqmTFvLYXZqfzo3LFelyMSF3xeFyDSmYIhx81/W0ZDIMTMiw6nICvV65JE4oKODKRH+f3ijSzdWsGFEwZw7hGaiE4kWgoD6TFW76xk5ry19MlJ446LxnldjkhcURhIj+APhLj56WX4gyHuuuQI8jN1ekikLRQG0iPcO3cNK0oruXTiIE4f08/rckTijsJA4t7itXv43eKNDCvM4vYLdHpIpD0UBhLXyqobuOnpZaQkG7/64pFkpekCOZH2iCoMzKzEzBaYWa2ZlZrZHWZ2yCGdZnasmf3RzNZHtltjZrebWXrnlC6JzjnHd/+2jLLqBm6ZMprxg/K9LkkkbrX6Z5SZFQDzgZXARcAI4F7CQXLbITa9LLLuXcA6YDxwZ+T5kg5VLQL8bvFGXl2zh5NHFvL1k4d7XY5IXIvmmPpaIAOY6pyrBOaZWS4wzczujrS15BfOubIm3y80s3rgd2Y21Dm3pWOlSyJ7c8Ne7n55Nf1y05hx6ZGamlqkg6I5TXQOMKfZh/4swgEx6WAbNQuCA5ZGngdEXaFIM7sq67n+r0tJMuPBLx1Nn5w0r0sSiXvRhMEYYHXTBufcVqA2sqwtTgBCwIY2bicCQGMwxLeefJ+y6gZ+eO5YJhb38rokkR4hmjAoACpaaC+PLIuKmRUR7mP4s3Nud5TbTDMzZ2autLQ02h8lPZRzjmn/WMG7m8s5b3x/rj6x2OuSRHqMbrm01MxSgaeBauA70W7nnJvmnDPnnA0YoDNLie7xN7fwxNtbGds/l7svGY+Z+glEOks0HcjlQF4L7QWRZYdk4d/Yx4FxwInOuVa3EWnutXV7uOOllRRmp/LwlRM1nkCkk0XzG7WaZn0DZjYYyKRZX8JB3Ef4ktQznXPRrC/yX9bvruKbT7xPshm/+8pEBuZneF2SSI8TzWmi2cBZZpbTpO0yoA5YdKgNzewHwLeAy51zr7e7SklYO/fXc+Wj71JVH+AXlxzBMUOj7qYSkTaIJgweAhqA58zsDDP7BjANmNH0ctPISONHmnz/JeDnhE8R7TCzzzR59OnUvZAeaX9dI1f98R12VNTx3bNGM/XoQV6XJNJjtXqayDlXbmaTgQeAFwlfWTSTcCA0f62mU1RMiTxfFXk0dTXwWFuLlcRR3xjkG4+/x+qdVVx5wlC+eeoIr0sS6dGi6oVzzq0ETm9lneJm31/Fp0NApFX+QIj/feJ93t60j3OPKOInF4zTlUMiXUyzlkpMaQyGuP6v77Ng9W5OHlnIjEuPJFlTTYh0OYWBxIxAMMS3n/qAOSt2ccLw3vzhiomkpxxyclwR6SQKA4kJ/kCIG2d9wD8//JjjinvxyFUKApHupJE74rk6f5DrnljCwjV7OK64F49efSyZqXprinQn/caJp6rqG7nmsfd4Z/M+Th3dh99++RgyUnVEINLdFAbimdKKOr762Lus3lnFeeP7M/PSI0n16cyliBcUBuKJ5Tv289XH3mV3VQNf+cxQpl04TlcNiXhIYSDdbs6KnXx71gfUB4Lcdt5YrjlpmMYRiHhMYSDdJhhyzJi3hgdf3UB6ShIPXX4MZ40r8rosEUFhIN2kvMbPDbOW8tq6Mob0yuShy4+hZECu12WJSITCQLrcvzeUcdNTy9hZWc/pY/oy89IjyctM8bosEWlCYSBdxh8IMWPeWn63eANJZtwyZRTfPPUwktRRLBJzFAbSJZbv2M/3nv2QFaWVDO2dyX2XHclRQ3QvApFYpTCQTlXfGORXC9bx+8UbCYYcl04cxE8uGEe2blMpEtP0GyqdwjnHglW7ufOfK9myt5ZBBRn8Yup4ThpZ6HVpIhIFhYF02Prd1dz50koWrd1DcpJxzUnDuHnKKM0vJBJH9Nsq7fbx/jp+vWAdT7+3nWDIcdJhhdx+QQkj++W0vrGIxBSFgbTZ7sp6fr94I4+/tQV/IMSIPlncevYYppT000hikTilMJCobd1by0OLN/DMe9vxB0MMyEvn22eOYupRA/Ela4I5kXimMJBDcs7x5oa9/PHfm1mwahchB0N6ZXLdqSOYevRA0nyablqkJ1AYSIv2VjfwwgelPPXuVtbuqgZg/KA8rjlpGOcd0V9HAiI9jMJAPlHfGGThmt28sLSUBat30Rh0pCQbF04YwFUnFnPU4Hz1CYj0UAqDBFfrD7B4bRlzV+xk7spdVDcEABhTlMOlEwdz8VED6ZWV6nGVItLVFAYJaHNZDa+t28PCNXt4fX0ZDYEQAAPzM7j8M0O5YEJ/Svrn6ihAJIEoDBLA9vJa3t64j3c27ePNjXvZuq/2k2Wj+mVzZkk/ziwpYsKgPAWASIJSGPQwVfWNrCyt5MPt+1m6rZz3t1Sws7L+k+U56T7OHlfEyaMKOWVkHwb3yvSwWhGJFQqDONUYDLF1Xy3rdlWxemcVa3dVserjKjaV1fzXeoXZaZxZ0o/jh/XiM8N7M7Z/ru41LCKfElUYmFkJcD9wAlABPAxMd84FW9kuD7gPuBhIAl4CbnDO7e1I0Ymizh9kR0Ut28rr2Lavli17a9myt4ZNZTVs2VtLIOT+a/3cdB+fHdGbwwfmcfjAPI4anM+gggyd+hGRVrUaBmZWAMwHVgIXASOAewl/uN/WyuZPA6OArwEh4C7gBeDk9pcc35xz1PqD7KvxU1bdQFm1nz1VDeyuqmd3VQO79tfz8f56dlbWs6/G3+Jr5Kb7GD8ojxF9shnRN5vRRTmMKcqhKDddH/wi0i7RHBlcC2QAU51zlcA8M8sFppnZ3ZG2TzGzE4ApwCTn3OJI2w7gbTM7wzk3v3N2ofs456hvDFHrD1DXGKSmIUiNP0BNQ/hRVR+gOvJcWddIVX2Aijo/++saqagNP/bV+vFHrt45mMzUZIry0inpn8ugggwG5mcwuFcmQ3tnUtw7i/zMFH3oi0iniiYMzgHmNPvQn0X4r/xJwIuH2G7XgSAAcM69Y2abIsu6PAxq/QH+/OYWAiGHPxAiEAoRCDr8wRCNwRD+QIjGYHhZQyCEPxiioTFIQyBEfWMQf+S5rjFIfWOI+kAQ51r/uS3JSfNRkJXK2P65FGSm0DsrjcLsVHplpdInJ42+Oen0zU2jKC+dnDSfPuxFpFtFEwZjgFeaNjjntppZbWTZwcJgDLC6hfZVkWVdrs4f5P9mt1TCoaUmJ5GWkkSaL5n0lCQKs9NITwl/nZnqIzM1mYzUZLJSfWSmhZ+z03xkp/uQKG8ZAAAHQElEQVTISfORk55CTrqPnHQf+Zmp5Kb7NH2DiMS0aMKggHCncXPlkWXt2W54FD8XM5sG3A7Qv3//aDb5LznpKTxy5UR8yUmkJBspyUmRR/jr1OQkUnzh5wMBkJqcpBu2i0jCielLS51z04BpABMnTmzzCZpUXxKTx/br5KpERHqeaM5dlAN5LbQXRJZ19nYiItLNogmD1TQ7x29mg4FMWu4TOOh2EQfrSxAREY9EEwazgbPMrOmNbS8D6oBFrWxXZGYnHWgws4mE+wtmt6NWERHpItGEwUNAA/CcmZ1hZt8gfB5/RtPLTc1svZk9cuB759ybwFzgcTObamYXA08Ar8fjGAMRkZ6s1TBwzpUDk4FkwpeRTgdmErnKpwlfZJ2mLiN89PAo8DiwBPhcx0oWEZHOFtXVRM65lcDpraxT3EJbBXB15CEiIjFKI6FERARz7Z1foZuZ2R5gSzs3HwCUdmI5Xuop+9JT9gO0L7Gop+wHdGxfhjrn+kSzYtyEQUeYmXPO9YhhxT1lX3rKfoD2JRb1lP2A7tsXnSYSERGFgYiIJE4YTPe6gE7UU/alp+wHaF9iUU/ZD+imfUmIPgMRETm0RDkyEBGRQ1AYiIiIwkBERBQGIiKCwkBERFAYiIgICRgGZjbUzP5qZvvMrNbMlpnZ2V7X1RFmdqOZOTN7xuta2srMcs1supm9Y2b7zWynmT1vZqO8ru1QzKzEzBZE3kOlZnaHmTWfwj3mmdkXzOwfZrbDzKrNbImZ/Y/XdXWUmQ2M7I8zs2yv62krM/OZ2ffNbJ2ZNZjZdjOb2ZU/M6oprHuKyO063wSWEZ5WuwY4Esjwsq6OMLO+hG82tMfjUtprCPB14BHgR4Rvp/oD4G0zG++c2+ZlcS0xswJgPrASuAgYAdxL+I+r2zwsrT1uAjYB3wHKgHOBJ82s0Dl3v6eVdcw9QDWQ5XUh7fQY4dsGTCd8m+DBQElX/sCEGnRmZrOAgcAk51zI63o6Q+TucqmE3yxlzrnPe1xSm5hZFhByztU1aesFbAXucc7F3EhSM/sBcCvhGSErI223Eg7loqZ3AIx1kQ/9smZtTwInOOeGeVRWh5jZKcALwM8Jh0KOc67a26qiFzlT8SIwIXIvmW6RMKeJzCwPmAr8pgcFwXHApcD3va6lvZxzNU2DINK2j/B05QO8qapV5wBzmn3ozyJ8hDnJm5Lap3kQRCwldv/tDylyqu5+4A7CRzrx6KvAK90ZBJBAYQAcDaQAzszeMLPGyHm4H5hZ3E11G6n5fuBu59wOr+vpTGbWBzgMWOt1LQcxhvCh+yecc1uB2siyeHcCsftv35prgTTgQa8L6YDjgbVm9oCZVUb6pZ4zsy4N6EQKg6LI8++A14AphO/N/FPgOq+K6oCrgX7AL70upAvcS/h872Me13EwBUBFC+3lkWVxy8wmAxcT/j+IK2bWG7gTuMk51+h1PR1QBFxFuD/zi4R/148Bnu/KP1zjugM5cuqnf2vrOedWAwf+EWc75w6cVnnVzAYR7rD8TddUGZ227Etk3f8Drm9+iiUWtPH/pfm21wGXA5c45/Z2QXlyEGZWDDwJ/N0595inxbTPz4C3nHP/8rqQDrLI46IDvwNm9jGwiHCn8oKu+KFxHQbAF4A/RLGeEf6rDeDVZsteAa42s1yPO/7asi8/JNzBOtfM8iPtPiAl8n2Vcy7YNWVGpS378p9vzC4kfOrre86557uisE5SDuS10F7Af95ncSXSaT+bcF/Nlz0up83MbBzhc+2nNPmdyIw855lZMBb/cDqIcmBjsz+GXgf8hK8o6pIwiOvTRM65h51z1tojsvqqyHPzw6wD33vaqdzGfRkNTCT8pjnwOBG4MPL1CZ7sREQb9wUAMzuRcCfsQ865e7ypPGqradY3ELlsOZNmfQnxwMwygZcIX5V2vnOu1uOS2mMk4T7BN/nP78SBfoPthP/IiBer+PTnFJG2Lvucivcjg6g55zab2QrCh1kPNVk0GdgQT5eeEb6W/b5mbfcB+4HbgY+6vaIOiPxV9yLwMnCDx+VEYzbwXTPLcc5VRdouA+oIH8rHDTPzAX8j/GH6Wefcbo9Laq/XgdOatZ0NfI/w2ImN3V5R+70ETG922e8phMNuWVf90EQbZ/A54FnCnWNzgVMJX5Z5hXPuCQ9L6zAzW0h8jjPoCywBHHAFUN9kcWV3X14Xjcigs5XAcuAuYDgwA7jPORdXg87M7PeEB/3dCLzTbPFS51xD91fVOczsKuCPxN84g1zC760dhMdK5BB+n612zp3ZVT83YY4MAJxzz5vZFYRHut5I+Lz7/8Z7EMS5EmBQ5Ovm/TmLCAd2THHOlUeuunmA8BFNBTCT8KCzeDMl8vyrFpYNAzZ3XykC4JyrNLPTgV8TPnXqB/5OeJR4l0moIwMREWlZXHcgi4hI51AYiIiIwkBERBQGIiKCwkBERFAYiIgICgMREUFhICIiwP8HEwaKZ4z9imIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x = np.linspace(-6, 6, 100)\n", "plt.plot(x, logistic(x))\n", "pass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Data set for classification**" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "n = 10000\n", "p = 2\n", "X, y = make_blobs(n_samples=n, n_features=p, centers=2, cluster_std=1.05, random_state=23)\n", "X = np.c_[np.ones(len(X)), X]\n", "y = y.astype('float')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Using gradient descent for classification by logistic regression**" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEECAYAAADOJIhPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VFX6+PHPSSGFBAhNOqEIEVBaRCw0wYJK+YoI6q4oawELlv2p2FkU664F3UUEu2IDpCiKgNJsdKVK78UAgQAJqef3x0lIZuZO6sy9M5nn/XrlJXPPzcxDDPe595TnKK01QgghQlOY0wEIIYRwjiQBIYQIYZIEhBAihEkSEEKIECZJQAghQpgkASGECGGSBIQQIoRJEhBCiBAmSUAIIUJYhNMBlKR27do6MTHR6TCEECKorFy58rDWuk5J5wV8EkhMTGTFihVOhyGEEEFFKbWrNOdJd5AQQoQwSQJCCBHCJAkIIUQIkyQghBAhLOAHhoU4fhxeeglWr4boaLjiCrjjDlDK6ciECH6SBERAS02Fvn3ht98Kj82YYV6/+65zcQlRWUh3kAhoL73kmgAAtIYpU2DJEmdiEqIykSQgAtrKldbHMzPh66/tjUWIykiSgAhoUVHe26pUsS8OISorSQIioF16qfXxmjXhlltsDUWISkkGhkVAOHoU3nnHdPNcdx0kJZnjo0bB8uUwdSpkZ5tjCQnw+OPQooVz8QpRWSittdMxFCs5OVlL7aDK7Z13YMwY2LvXvK5eHW69FV55xUwD1RrmzoV580z30LBh0Lq1oyELEfCUUiu11sklnidJQDhp925IToaUFNfjYWFmCuiwYc7EJUSwK20SkDEB4ahJkzwTAEBeHjzyCBw5Yn9MQoQSSQLCUadOeW87dAjuvtu+WIQIRZIEhKO6dYPwcO/tCxaYVcNCCP+QJCAcNXAg9Orlvf3oUfMlhPAPSQLCp7SGTZtg61bz55IoBTNnQh0vm+CddRY0aWL+/OOPcO210LYtXHIJjBtnxg6sLFkC991nupNmzixdLEKEIlknIHxmxgx44QVYscJ08XTtCs88A927F/99sbFw5ZXw0UeebcePmxlEO3fCkCGuA8U//QRbtsD777t+zyOPwPjxcPq0eT1hAnTqZM5r164Cf0EhKiGZIip8Yu1auOwyM5hbVPPm8Ouv3u/08/IgK8us/v38c+tzHnkEpk0zTxfuwsPhl19g1SqoWhUaNoSrr4aMDM9zo6PNQrR334XIyDL99YQIOqWdIipPAsIn3nrLMwEAbN8Ob7wBY8e6Hs/KMhf3774zd/cFq4GtLF1qnQAAcnOhR4/Ci35CgnUCAPNk8PHHUL++qU7qD/v2mS6uBg388/5C+JqMCQifsEoABQ4edH2dmgqdO8Nrr5nxg5QUOHbM+/evXl38Zxe96JdmJtF335V8Tln9+KMZ4G7Vynz17i2lrkVwkCcB4RMNG3pva9y48M9Hj5qL5bp1pX/v9PTyx2Xl2DEzUOyrncl27TLdWbt3Fx774QfYscOMW9Sv75vPEcIf5ElAeNizB557zsy+2bmzdN9zzz2Fs3iKSkoyReAKvPgi/PGHT8IsN2/jE+X15puuCaDAjh2mK0yIQCZJQLh47jlo08ZU6XziCVOo7fbbC9sLdvUaNMjM8X/2WbPqt0YNMwgckf9sqZRJAFOmmIJwBVat8l2sYeX87f39d+jf3zxhpKWZaaTt2kHLljB4sJndVBYFhe/K2iZEIJDuIHHGkiWmmmfRQdqsLJg82VzQH3jATPtcvrywfeZMMx0zL891+mbBeoFRo8w00blzTd/+zz/7Jta4OO8DwCXJzTW7kvXpY7qlTpwobNu2zexm9s03cM455u+xdKl5Ourb1ww8uyuuK6y4NiECgUwRFWcMHw7vvWfd1qSJmQL6zjv2xuSUTp3MjKcHHzRTXHNyoFEjuPlm8/RTdDxh504zQ8m9Sygx0YwJyEwh4QSZIirKbPNm722HD5uNXULFqlXmSSEtrfDY3r1mamnDhjByJEycCLNmmXPatIG6dc2ThVLmienppyUBiMAnSUCccd555s7VStWq1iWfK7OiCaBATg589ZVZtzB+vOlaKtCsGUyfbn6O0g0kgoUMDIsz/vMfiI+3brv4YntjCWT798OHH7omADCzgT7+WBKACC7yJCCYMwc++cQstOrVCxYvLly8FRZmEsD69c7GGEgiIrxvdvP77/bGIkRFSRIIca++Ck8+6bq5S7Nmps87Otr0bU+eLKtfC4SFmUFjb2sdqla1Nx4hKkqSQAg7dcosZnLf3WvHDjNIPHWqafvb35yJLxDl5ZlpsbGx1iuZi9sbQYhAJGMCIWzGDHPBt7JypUkEAweG3oBwSVJTrRNAw4amZpAQwUSSQAjzNggMZhbMFVfA/Pn2xRPs9u0zT01z5zodiRClJ0kgBOXlmUVhr71m+v2t7N1b+rpBotBff5mpo0IEizInAaVUS6XURKXUH0qpXKXUQotzlFLqMaXUHqVUhlJqsVKqg08iFhWSl2fq/gwfbsofF+y+JXxn0SLzJCVEMCjPk0Bb4CrgT8DbGtPRwJPAi0A/4CQwXylVrzxBCt/55BMzFmClWrXyF2UThU6dgqeeKnx98qTZQKd7d/P10ENm20whAkF5ZgfN1lrPBFBKTQVqF21USkVjksDzWus384/9AuwE7gGeqEjAomKK21Dl5EnvG7eLsnntNfOzHD0a/u//YOHCwrYlS0w9ou+/h5gYx0IUAijHk4DWuqTLxEVANeCLIt9zCpgN9C3r5wnfKq6bIi+vsBS0qJiMDLN3QteurgmgwNKl8L//2R6WEB788fCfBOQCW9yOb8xvEw4aONB7W6NGcmfqa3/+6b3Nl3srCFFe/kgCCcBJrbVbZRVSgVilVBU/fKYopaFD4aKLPI+Hh5v9AopuACP8KzbW6QiECNApokqpMUoprZTS+/fvdzqcSkUp0xXx0ENmj4CaNaFjRxgxwlQQlYVhvmc12B4TYxKyEE7zRxJIBeKUUuFuxxOAdK11VklvoLUeo7VWWmvVQAqyl9s775gyBi1aQLduZi9cMIngpZfMBulHjsB115lzp0+HzExnY66M3Afba9Y0s4VkdbEIBP4YBtwEhAMtMdNICyTltwkbjB8PDz9ceFHfvt3MSDl2zOwdXCAtzWyOIusF7NG9O3zwgdl1TIhA4I8ngZ+BNGBwwQGlVCxmvcC3fvg84SY316wIdr+rz8kxdfCLXvBnzPDcFlH4z6ZN8O678PjjvttvWYiKKPOTQP4F/ar8lw2Bakqp6/Jfz9FapyulXgCeVEqlYu7+H8QknDd8ELMowaFD3melbNliCsOdd555XauW6R4K8K2mK42//oJnnjF/fuUVMy7wzjuySE84pzzdQXWBL92OFbxuhlkU9gLmov8oUAtYAVymtT5UvjBFWVSvbi7ue/d6tiUkQP36ha+vvBLat4c1a+yLTxinT8P775uB+VGjnI5GhKryLBbbWTBoa/G1M/8crbUep7VupLWO0Vp301qv9nn0wlLVqmaTdCu9e0OdOoWvw8Ndk4Kwn1QdFU6S9aGVxB9/mJWprVqZEtBvvGEGgefNM7VsYmKgZ094663C78nKgltuMeULhHMyMpyOQIQySQJBLjMThg2Db74xtX/Cw81isHfega++gtWr4ZdfoHNnuOAC1+8dPRo+/dSZuEORt7GXjh3tj0WIApIEgtzo0fD554Wvc3NNgbKRI82GMB07er/IyBOAvawSQIcOZuGeEE6RJBDEtPZeFXTpUrNFZOfO1u15eaa7SNjrkksgLs4MCnfoYNZy1JMC68JBkgSCWF6e97r0mZmwbZv3JBAWBklJZktEYZ/LL4cnn3Q6CiEKyezkIBYebi7kVho08D5DqMDIkVCjhu/jEtaaNYO773Y6CiFcSRIIciNGeFb+DAuDG24wNWoKZGebmUKLFxfWshk0CP7zH6hdG+Fn9eubmVlF/58IEQikOyjIXX89REbCpEmmPlDt2mbPgH/+s/CcDz6Al1+G9evNDJVOnWDsWLjqKvjiCzh82Ln4Q8WwYaYrSIhAo3SA1wtITk7WK1ascDqMoLVsGfTtC0ePuh5v3NjcmQ4aJMXj7LB0KezZAwsWmCe1/v1NElbK6chEZaWUWqm1Ti7pPHkSqOQmT/ZMAGAuSBMmSAKwQ4MGpoz3558XThN97z24/XZzXBKBcJKMCVRyxXX1xMdLH7W/RUWZi/1nn7muE8jONgl63jznYhMBKjMT7rsPHnzQlg0+JAlUcsXVre/UCQYMsC2UkHTvvd4rumZlwcyZ9sYjAtiUKeaxMDrabAjy6quwdavfP1aSQCV3333QvLnn8fbtzRTRhg3tjymUbN1afG0g913HRIhZv97MHVYKbrrJte3tt6FtW7+HIEmgkmva1NQH6t/fTFNs1MjMKJo61dxwTJ/udISV24wZ3jePCQ83g/YixKSlmTncSkG7drBzZ2Hb0KFmBajWph/RBjIwHAK6dDHdDtnZ5vcuIv//+tGjsH+/s7GFgpQUqFbN/Nsv6sYboV8/Z2ISNtPadPHcf79nW2IizJ5tEoIDJAmEkMhI19fVq5uZK1JDyP/cEwAUJmVRif30k6ntfuqUZ9snn5g7AYdJd1AICw836wSEM7791uw5LCqZQ4dMzRalTMXAoglg1CgzL1vrgEgAIE8CIW/MGLMB/Suv2DIbrVLq0MFUAk1IgOXLSz+h4/hxWLTIe/0nEURycsw/pnHjPNu6dDFL85s2tT2s0pAkEOLCwuC550zX0BNPmN9lUXoJCeZCXq2aeT17tqnnVHSsJTraelFedLRJICKIzZ5tZl1YmTs3KGqFSBIQgKlrHxZmuik3bjRz2EXJ7ryzMAGAGeg95xyzGvvQITPmFxYGzz7rualMjx6eu72JILBtG1x7rdnT1d0zz8Cjj5q+1iAhSUAApvvyoYfgttvMJAWZNVSyiAjT5euuZUtTnbWA1qar7bPPYPdus6lMr14wcaJ9sYoKysgwM3veftuz7eqr4f33g7YcryQB4WLpUkkA7tq3h82bPRd9XXaZKQJXEqXgxRfh8cfNmEFiIrRo4ZdQhS9pbS7uw4d7ttWqBXPmmP7+ICdJQLhITISYmOJXuVZWVhvBd+liZvktWwb//jf8/jtUrWq6cl58sWxTPKtVg969fRuz8IPVq012P3jQs+2tt+COOyrV3F5JAsLFueeaLo5QK2wWFmYGeQFSU80Fu1s3M2sqIgIuukhWV1dqqanmjn/GDM+2W24x5V6rVrU9LDvIOgHhYcIE6N69cGWx+yKzyigvD44cMV8dO5qSLrNmmf59UUnl5RU+ztWs6ZoAkpLMIg6tTd3vSpoAQJKAsNCiBSxcaGa/vfqqqX0zdqyZ9RIZae6aK7OVK+Gbb5yOQvjNDz+YO5zwcBg92rVt6lRz4d+4EVq3diY+m0l3kLCkFFx5pfkCSE42/14OHjSb08+aZfrK27aFe+5xNlZ/WL/e6QiET+3bZ4qzLV3q2fbww2YObyg88lqQJCBKLTLSbEsJpurtTTeZJ+rnnvM+o8hqsDUYBOlsP1FUVhY89pjrfN0C3bub8roNGtgfV4Cp5A/2wt/Cwsz6Am/atLEvltKIiyt5HU9iItx1ly3hCH+YOtXcfURFuSaAyEj48UdzV7JokSSAfJIERIXdf7+pheU+a65dO1My5eqrPccRwsJKf7cdF+e9LSEBunaFs88u3Xu1bg1//7t1m1JmUHjCBNl2M+hs2mT+5yoFgwe7tr30EuTmmieDnj0dCS+QSXeQ8ImPPzbjBzNnwokTJgH885/mZmv2bPj+ezOGsGsXdO5s6uvk5cGkSabu1tGj3t+7c2dz42alZUto1crU7N+505RnLk6XLvDGG9CkianieeyYeY+OHc2U0Msvr/wD35XGyZNw993w4YeebYMGmV+ugnm/wiulA7zDNjk5Wa9YscLpMIQftW0LGzZYtyUkmDpcPXtCerpne1hY2bZofPJJM9NJBCmtzYItq/66Bg3MtC6pygeAUmql1jq5pPP8cs+jlLpFKaUtvkb44/NEcGvSxHvbu+/C+eeb/ZCrVPFsL+sevfXqle18ESB++83cEYSFeSaA9983yWHfPkkA5eDv7qBLgaIFCLb7+fNEELr1VliyxHPzpauvhoEDzZ///W9TcmH6dFOMbf58OHCgbJ/Tpo11GRgRoFJSYNgw02/nbsQIs5w7Jsb+uCoZfyeB5Vrrk37+DBHkrr/edO9OnGjW6NSoYS74r7/uel7fvuYrPb1sd/RKmXGFV14xNfxFAMvKMtsxLlzo2dahA0ybBs2b2x5WZSYDwyIgDB9ungiOHDEr9Iu7wQsLM+ecOFH8ezZqZBa4JSaa5CEDvgFs3Dizq5GV2bPhmmvsjSeE+DsJbFNK1QK2Aa9oraWCuvBKqdJNG42OhgsvhK++8n5OtWpm/cLdd/suPuFjS5eau36rEf9u3cyUMnl08zt/JYEDwJPAMiAcGAq8pZSK1Vq/WtI3K6XGAE8D1K9f308himA2bpyp8V+0vEN0tOn2adMGbr7ZesOXkuTlmbUNixebtUY33WRKZggfOXbMrNZdu9azrVYt0w3Urp3tYYUy26aIKqU+B/oAdbTWpZ7TIVNEhTdpaabC77ZtULeuqWHUsGH53y8724xPzJxZWOoiPt50KT32mG9iDlkjR5qpnVYmTjQ1+oVPlXaKqJ1JYDDwBdBCa13qWUKSBIRdXnjBbA/rrkYNM0OxVSv7Ywpq06bBkCFmta67666Dzz+XgRo/cnSdgBfa7b9CBBRvq5KPHYOPPrI3lqC1e7dZ+KGUudAXTQCJibB3r3nM+vJLSQABws7/C9cBh4FdNn6mEKVmdcNaICfHvjiCTl6eWdChFDRtCnv2FLZFRBT2r+3YUbH+OuEX/loxPE0p9YhSqq9S6hql1EfAEGBsWcYDhLCTtwHg2FhTika4+e9/zYU/PNxc6Iu6915z4c/Ohv79nYlPlIq/Zgf9CQwHGgMK2ADcrLWWh2oRsB591KxcLrrvSHi4WcMgM4Ty/fEH9Ohh+sjcdexoZvdUq2Z7WKL8/JIEtNaPATKfQgSV+HhTrO7112H5cjPltF8/syFVSDt92mTCTz/1bIuLg3nzTD1vEZRkxbAQRcTGWs8QCkkTJnjfXeeFF+CRR+yNR/iFJAEhRKEVK2DAAOv9Qm+/3YwDhOhevJWVJAEhQt3x42Y3rnnzPNs6dYIZMwo3lxaVjkzUFSIU5eXBU0+Z2T01argmgKpVYc4cM7tn5UpJAJWcPAkIEUrmzjWLuE5aVHh//HGz7Zos4gopkgSEqOy2bTMLHX7/3bOtTx9TMU/24g1ZkgSEqIxyckwN/rlzPdvq1TP9/BdcYH9cIuDIc58QlcnLL5t+/shIzwQwfrzp5z9wQBKAOEOeBIQIdsuWmf04rfr5u3aFBQvMAgghLEgSECIYnTxpyjesWuXZVqMG/Pij2ZNXiBJId5AQweSBB0x3T3y8ZwIo6O5JTZUEIEpNngSECHRffw3XXmsqcrrr3x+mTzeV7oQoB0kCQgSiAwfg4otNDX53jRqZUqdNm9ofl6h0pDtIiECRl2e2Y1QKGjRwTQDh4WY+v9Zm0xZJAMJHJAkI4bTJkws3Z/niC9e2O+4wF/6cHFPfRwgfk+4gIZywYQN07w5Hjni2tWtndrepUcP+uETIkScBIeySlWWmdSoFbdu6JoDYWLPTvdawdq0kAGEbSQJC+NvYsebCHxUFixe7tv3rX+bCf+qUeTIQwmbSHSSEPyxZAldcARkZnm09esD330OVKvbHJYQbSQJC+MqxY3DJJbB+vWdb7dqmu6dNG/vjEqIY0h0kREXdeafp7klIcE0ASsGkSaa7JyVFEoAISPIkIER5TJ0KQ4dCbq5n25AhMGWKbM4igoIkASFKa9cu092zd69nW7NmZhVvgwb2xyVEBcitihDFyc019XmUgsRE1wQQGQmzZ5vunu3bJQGIoCRJQAgr48ebC39EhLnQF3X//ebCn5Vldu8SIohJd5AQBVavhl694Phxz7bOnWHhQoiLsz0sIfxJngREaEtPN7tvKQWdOrkmgPh4s2uX1rBihSQAUSlJEhCh6dFHzYW/alX47TfXtpdeMhf+tDQ4/3xn4hPCJtIdJELH999Dv36mL9/dlVeavv8I+SchQov8xovKLSXFbM6yZYtnW716Zlpnixb2xyVEgJDuIFH55OXBzTeb7p66dV0TQFgYfPyx6e45cEASgAh5kgRE5fHhh+YiHx4OH33k2nbLLWbOf24u3HSTI+EJEYj8lgSUUm2UUguUUulKqf1KqbFKKdkNW/jWli1w1lnmrn/YMHOHX6B1a9MdpDW8956UcRDCgl/GBJRSCcB8YAMwAGgB/AeTdJ7wx2eKEJKdDVddBfPne7ZFRcE330Dv3vbHJUQQ8tet0QggBrhWaz1Pa/0W8C/gQaVUNT99pqjsXnjB3PFXqeKZAB591Nzxnz4tCUCIMvDX7KC+wFytdVqRY58BLwI9gNmW3yWEu19/hT59zM5b7i66CBYsgOho++MSopLwVxJIAn4oekBrvVsplZ7fJklAeJeWZnbfWrPGsy0hwZRvOO8828MKJBrNfLbzNZvJJJd+tOJqWjkdlghC/koCCcAxi+Op+W3FUkqNAZ4GqF+/vk8DEwFs1Ch44w3rtjffhLvv9svH/s5BZvEnsURyG52ojn+eLA5xkk9ZRyRh3Ex74okq1/toNCP4mndZRU7+sYmspBk1mMFQzuMs3wUtKr2AXCymtR4DjAFITk7WxZ4sgtusWTBoEOTkeLYNHAjTpvltVo9GM5JvmMJaTmBWEb/Ob4zjUv5Oe59+1rMs4k2WcwjTrfUyP/Mk3fkHncr8XlPZwCRW4f4PYwfH6M8U1nE3cZR//+JdHON5lrKag0QRTk8SeYLuVEEm91VG/koCqUB1i+MJ+W0ilO3bZ1bx7trl2dakidmkvUkTv4cxgeW8zUqXi+ke0hjNAq6kJXWo6pPPmcMWxrGU0xQmul0c5xHm042mtKJWmd7vG7Z4JIDC903jDZbxKJeUK9b9nKAfn7KWv84cW8Ju/uAQXzEEhXI5fx9pjGURy9lPGIoLacRYepFATLk+X9jPX7ODNmH6/s9QSjUGYvPbRKjJy4PrrjOzexo1ck0A4eFmu0atzXEbEgDAHLZaXkz3c4K3Wemzz/mcdS4JoMARMpjMqjK/Xw55xbbvrMB91r/52SUBFPiazcxms8ux45ymP5/yNqtYzUFWcoA3WU4/PiXT4u8rApO/ksC3wBVKqfgix4YAGcAiP32mCEQTJ5oLf3i46dopasQIc+HPyTFdQjZLJ9tr20ksisyVUxqZXttOFNPmTTeKT5Lud+tlsc4iAQDkolmM65PbK/zCKg56nPsTe3i7HMlNOMNfSeAtIBOYrpTqo5S6A9PH/4rbtFFRGa1bB7VqmYv/iBGubeedB6mp5uI/YYIz8eVrSx3L41UIpzfNffY5bbx8DkAnPCc+nCaHz1nPdDaSjedG9rfSkSrF/NPtRL3yBQrFjiW4t23ksNdzf7dIDiIw+SUJaK1Tgd5AOGY66L+AV8mf8SMqodOnoVs3c+E/91w4erSwLTbW9PNrDb//DjVqOBdnEQ9xsWUiGEhretPMZ5/zIBfS3mLGTg+acgsdzrz+k8M8xgLOYwJDmcogvqAjE/mS9S7fV4VwkrHez7g6VejL2eWO9RpaEW7xJFGPqtxJZ5djxc1uqlaBgWlhL7/NDtJabwAu9df7iwDx9NMwdqx12zPPwBOBWyWkCdWZzQ28xE+s4RAxRHApzRjNJRXqUnFXi1hmcQPPspgV7CecMC6mMWPpSSThrGQ/DzGPJewix22UYj0p3Mu3dOAszqb2meO30pHVHCTDre/9WtrQ2HJORuncSgfWcoj3WMPx/K6qxlTj75zLKL5lE0eoThRX0oIhtOEz1nl0q4UB20jlACeoT7zFp4hAorQO7BmYycnJesWKFU6HIYpauBD69jV3/+4uvRS+/daUdhBnaDRfsZGl7CGGCG6nE4kkkEE2XZjstS++QDyRDKINT9CNnRynDXWYzWbeZw1bOUptYrmKs3mBPkS4PeAvZCfvsIo9pNGQeIbTscTurj85zDQ2EEMkbajDcGaxnxMu5wynA62pzb/5iRQyPN6jKw1ZxK0ytdQhSqmVWuvkEs+TJCBKJS0NhgyB777zbKtTBxYvhqQkzzZBFrkM5gu+YQu5+Xf6dYjldjrxHVstB1e9CQPyMF1CF9GIaVzPG/zGInaxjxM0oToT6UdzEsgkhxF8zRTWklVkRlEC0YynL3+jdKuuB/MFU9nocbw6UfzCP3iNX70OBP+PqxlJidch4QelTQIBuVhMBAitTZfO0xZDOUrB5MkwfLj9cQWZcSxmltv0yhTSeZ6lXuf7e1NwKc8il4XsogmvcapId8xmjtKK8bSnHls5SprFLKdUTvM6v3Ej5xJWim4vbwPAx8lkDls9nhCKWl/CE45wniQB4Wn+fDNlM81iItfo0TBunNTmL4PF7LY87otn8FMW01xzocSnizUcYAeptKBmiZ9RvZgB4AbEU5NYr+21iiwaK6h3tJa/6Ex9epBY4mcL/5MkIIyDB2HAAFi2zLOtVy+zmKtmyReMULaMfbzHao5ymlbU4gEuoCaxZFlM83RaDJFULeUMnitpyc/s9TjegXoMpg11iOUrNp4pvVGgCdW5my4AHOAEf2M6S9hNNnlUIZxLacYUrpXVxQ6TJBDKcnNN0bb//c+zrW5dmDEDLrzQ/riC0ERWMJoFHKNwsHwGm5jOEKICcBfXGCK5ik9I4RQac8H+Fz25DNc9lzPJ4SyqcjY12ckxsvM7pDpQj//SlwjC6ENznuNSXmcZWzFTgxOpwUNcRN380ht3MYcf2HnmfbPI5Tu2ci/f8jHX2vJ3FtZkYDgUffaZ2Wc3z6L8wCuvwAMP2B9TEMsgm3OZwDaLcg3daMwqDlp22wSiAbTiKXrwGr8yi82kkenRbVWXWG6lIw9y4ZmLPMBJMrmWz1nOfo6RSVUi6UUiz9ObS3jODQR7AAAUUElEQVTvzJTTohoQz5/cU6GCd8KaDAwLV3/+Cf36mT153V13ndmkPUYey8tjPMssEwDASg6QHkR1dGayme/Z7rH+oKi/SOdFfmI6G/mUQXTOX7j2PEuZx44z550im6/ZwnEyLRMAQCoZpJEpScBBgfecKnzn9GkYPNjM5ElKck0ALVvCxo1mBtCXX0oCKEYmOZwmh6ms9ygGl0kOE1ju9XtLKvYWiIpLAEVt4Sj/yi8FptEeBeYKrGQ/zbFeJX4OdahHXPkCFT4hTwKV0euvw/33ex4PC4OPPoIbb7Q/piD0JRv4H8tYxUEyyD7TH96KmjzERdxGZyazil0c9/oe1YkmhXS7QrbdEnZzkiyiifD690wnh4tozD5OkFlkkLwqkdxOp1JNUxX+I0mgskhJMat1163zbLvrLhg/3lTyFKXyDZu5k9mk4rkqejNHeYj5nMtZXruBClTmBACmnPQK9tGTZjQngYOc9DinLrG8SB960YwprOUAJ2hMdYbRnhs414GoRVGSBILZqVNmy8UPPoCICNfdubp0ga++ggbWhcaEtdPk8Aa/8Rq/WiaAAsc4zb18i29m+wcvDTzBjywhkVtpz2oOeHQnXUMr6hNPU6pzKx3oR2uqlXNrTeF7kgSCjdamRv/Ika7H8/Lg9tvNGMBllzkTW5A7RRZXM4VFWOx4ZmE5+/0cUXD4iT3cx3e8zpVo4D3WsJnDnMzvQnuPNbzP7+TlJ8wmVOcfdOApelq+3zQ28iXrSSWDJGrzABeS6GVMQVScJIFgsXw5XHmla4nmAu+9B8OGmQFgUW4vsLTUCUC4eoNlHCWD1tQigxxSyTxz0QczcFxgN8cZxxISSeBmt72cn2ER41hyZuzge7bzHVuZxhDaUdeev0yIkdlBgezwYbj6anNx79LFNQHceSekp5sng1tukQRQBhNYTk/epwWvk8zbvJNf/Ow39jkcWXD7hLU8xULWcNAlAVjJIo8v2eBy7AjpTGCFy+AxmDGYF1ji83iFIU8CgSY3F557Dp56yrPtvPNg+nRo0cKzTZTKsyxiLIvPzPSBY9zGbD5jvcxSsdlRt/LTn7OeAxYDywArZacyv5EkECjmzjXdPVZmzTILvUSFZJDNB/xRJAEUms92etHUgahCl/vageIWjEXLpcpvpDvISbt2mW4epTwTwBNPQHa26e6RBOAT6/jrTG0bK3+RztW0tDGi0FWXWO7JLy5XYAhtSaKW5fmX0MSOsEKSJAG7ZWbCvfeaC39iohnwLdCnj6nmWVDHP0LufnypAfHEFHNHeYR0pnI9d9KJphXYolGUrAeJnEMdXuZnnuAHfmAHUUQwjt40otqZ8xRwKYk8JzvV+o1cZewyZYop2uYuLs7s1nXxxfbHFGIaUo1LaMI8tlu2N6UGQ5nGLP4M8dn//neIU7TnLXZyDICX+Zn+tGIKg+hGEyawglQySKYBN5Ry8xtRPpIE/GndOtOVs3OnZ9vrrxc+EQjbvM8AOjPJY2VrFcKoTxwz+NOhyELLWg65LMbLIpepbKQNi/kXvXiKHg5GF1qkO8jX0tLghhvMxf3cc10TwI03wvHjprtn1ChJAA5oQDV2ch//RxK1iSGOSDpTn39zuUthOOE/9ajqdTX2D0WqkAp7yJOAL2ht7uyt6vA3a2Zm97RrZ39cwlIUEUxnCDnkcZIsqhFFGMprFUxRcRGEUYMoutOU0+Qwh62W5xW3X7HwD3kSqIilS6FqVVOd0z0BTJliksP27ZIAApS5MEWf6W/uQkOv5yYQbVdYlVIOeVQnmnfoTxW8FzIMI4w1HGRbMbO4hG9JEiirgwehd2/TldOtm1m1W+C++0wNf61Nl5AIKg9zEdW8zFVPD5KdwQLZNlJ5k+XFFo/7i5N0ZCLtmEAfPmCV1GfyO0kCpZGTA48/bi789evDDz8Utl1wgZnvrzW89hpESXXEYJWLRnmZheJeykCUzxoOsK+YLp+0/M3qT5PDAnYyjBlkSAL2KxkTKM6sWTBggHXb3Llw+eX2xiOKdYosNpJCE2q47H1bWtFEEEcVr1shioqbzWayyrDb2jpSeJtV3McFfowqtMmTgLutW82sHqU8E8Czz5qnAq0lAQQQjeZR5tOOCZzPZJJ4kyF8yTG32jQliSGSbrIy1W8UWCaAVtTkHC8rhQH2FrNzm6g4SQJg+vXvuMNc+M8+23V3rmuuMdU8tTZdQrI7V8B5gaW8yE9nFh6lcpov2MCtzCz2+zaSwiRWspoDZ469whXFDlyK8lF4336nGlF0oZHX7z27mAQhKi50u4O0hnffhdtu82yrXRvmzIHzz7c/LlFmU9loeYGZzw42kEIb6rgcTyebW5jBd2zlBFnEEEFvmvEB/0d94ulGYxaw05bYQ0VxK7BPkMUIkvmWLfzlth3n+TTgVjr4N7gQF3pPAqtWQb16ZlqnewKYONHs0JWSIgkgSOShLfe1BThJlstdfoFRfMuXbOBE/iBkBjl8zRZG8nV++wUh+A/DOfWJoyuNmEQ/etCU6kRxFlUZRBKfcx2R8mTmV6HxJJCaCsOHw4wZnm233gpvvGHm+4ugE4aiCdUtFxklEM2Fbt0M6WR7rR00j+2kcIr+JHEbnXg7f7MZ4T8KGI+poNufJPqTxGHSiSKceNmH2BZ+ueFRSi1USmmLL/tW3OTlwYsvmn7+mjVdE0CbNvDnn4VdQpIAgtrfOJcqFr/KV3E2zanpcuw4pzni1uVQIJXT7CUNgIn04zY6Sh17P6tFDG05y+VYbWIlAdjIn0+9PwIXun3ZN/fuqqtg9GjXY9OmmQv/+vXQqpVtoQj/upsuPE8fOnAWcVShGTUYQWcm09/j3LpU9TrQ2JIEkqh95vXLXM6v/IOLihm0FBVzmAxe5ienwwhp/rzNOaq1/tWP71+8AQPMXP6HHzZTOyMjHQtF+N+DXMh9XMBh0qlOtNc7+HDCGEZ7NpLisgAsHMVNnEsmOdzDHKaxkZNkkYumlpSM8Kuv2cIjXOJ0GCGr8j7rjhxpvkTICCeMs4gr8bz76UoMEXzEH2zhCLWIZTgduY8L6MNHLGaXy/lHvFS8FL5x1Ev3nLCHP5PA5Uqpgv+7S4CHtNZ/+PHzhCi1CMI4QgZ/kU4K6XzKOvaS5pEAhP+d7TZuI+zlrzGBRcB9wBXAHUATYIlSKrE036yUGlMwmLx/vxSQEr61lF38P+axicOAmcO+igO8w2pnA6vEYrxcasJQ3O2217CwV6meBJRS1YH6JZ2ntd6U/9+nixxeopSaD2wC7s//Kul9xgBjAJKTk2WnP+FT77OGYxZdPCfz1w0I34ohgmtoxZds8GgLR8lMIIeVtjtoMDCpFOdZlmDUWh9USv0EdCptYEL4y+Fi+qBjiSBddhjzqSRqk+2laFw2eXzJerrKDCzHlKo7SGs9WWutSvoq6W0ofvW4ELZIJMFr22Da0sqtjzq2Es+fsMPFNEIX80+/uDbhf7asjldK1QMuAVba8XlCFOc+LqC5RSLoRH0mcDV/MJIPGcirXMEa7mAkUkKkIk6Ty8U0tmyLIpyBJNkckSjK50lAKXWeUuobpdQtSqleSqlhwEIgD3jN158nRFk1I4HPGMQAWtOQeJpSnRtox3SuJ4ZIoojg77TnfrrSnvpcQQvLFcmidH5jL0Noy+W0cDkejuIfdKQ7ic4EJgD/TBE9ghkbeB6oBZzAJIGBWuvdfvg8IcrsfBoyg6HkkEcY6sw+w1b60JyzqcV6UmyMsPJYSwpJ/Jd6xNGdJtQnjmgiuYazGUQbp8MLeT5PAlrrfcBVvn5fIfwhooQ7fI3mTmZLAqigDHLYwTF2cIxhtOd9Bjodksgnz7hCFGMKa5ks6wd86ms2y25hAUSSgBDF+IYtMnfFx46QwTJkEWigkCQgRDGyihSZE75RnSjau5WPFs6RJCBEMc6nodMhVDp9aE4LqRcUMCQJCFGM+7iAnjQt8/fFU4WEEC9BXZdYRtKZRKoDZgOZobTlXYt9HoRzZCmkEMWIJoJvuIl/8xOL2MWv7C2xrEQM4czn79zKTFJDuAz143RjFF05RRYbSaEpNaiD7OIXaCQJCFGCWCJ5ip4ALGMfY1nIMvaThyaDbI+kkEEuV/NpsTWKQsGXbOBSmtOOuiRLt1rAku4gIcqgCw35mpvYxih2cj/XeVnsFOoJAGApe7id2VIbKMBJEhCiHOKJIo4qbCfV6VAcpYCbaOd1O89l7GU2m+0NSpSJJAEhKmB3iC96iqcKkxlAk/zBX3d5wI4QT5SBTpKAEBXgrU5+qEgji0msJIlalu3ViKIPzW2OSpSFJAEhKiCWSKdDcNx2UhlBMrWI8WjrT2vaUteBqERpSRIQogI6UK/C71HX4uLpD9FEEFfOpFXcxjpNqUFfzuZ9BnIlLUikBh2oxyNcLGsCgoBMERWiAh7iIlZygJ0c82irlj94vJ8Tlt8bgeJBLmQ7qUxlo1/jDANOl3PbzASiuYtkXuc3TpLt0nYudbmTzgBcQyuuoVVFQxU2kyQgRAVcQCNmM5Tx/MYWUomnCjWJoTW1GEASSdRmLIt4k2WkkkEe5uLfjrrMYChNqcE9fOPXGBsQx35Olvv7UznNOJZSj6rUJJb9pBFFBJfQhJe4jBjpEgtqkgSEqKB2nMXbxXR7jKEnY+jJVo7yC3voRH2XfvLi9jwurTiqcJIsy7bMUjwB1CSKY2QWO8x9kFPEEckk+nE5LWlAfDmjFYFExgSEsElLavJ32nsMlI4kmTbUKdd7xubv0NWrmC0ao0pxp36aXJpSgxjCiz3vJNksYIckgEpEkoAQDlrCLt5iBXfRmT40I5YIIlB0oWGxBehqE8N0rmcfDzKbG/knF1qefw61GUjrEuNIz9/5qyfNeILudC2mzMMRMkr3lxNBQbqDhHDAcU5zI9NYwA4yySUCxcU04UeGUZMYWlCT51nK4/zg8b2xRLKbB1z64nuQyOtcySv8yhoOEkU4F9KIF7mMNtRhDlstB6/dreQAn3AtjanGr+yzPKe5D7qvROBQWgd2XY/k5GS9YsUKp8MQwqeG8RUf8ofH8X60YhY3AGZ/44eZxyRWcZxMwoA21GE2N3gdR8gljz84RBxVOLvIAq5scmnNm+woRSLYwF20pCY9+YCf2ePS1pwafMffOVv2Awh4SqmVWuvkEs+TJCCEvTLIpjVvsoc0j7Z4qrCBu2hUpAzDCTL5mT00pQZJ1C735+aheZ81fMzvLGM/p9yme4IZt1jLSKKJ4BAnGc18lrKbLPJIpj6juUQ22gkSpU0C0h0khM1Okc0xL/sMnCCLA5x0SQLxRHEFLSv8uWEohtOR4XTkeZbyFD+QU6TCZxgwlLZnisGdRRzvMfBMFVCFqnAMIvBIEhDCZrWIIYnaLLfYbL0VNTnXhv13R3MxMUTwKevYy3HqE89g2vAwF3ucKxf/yk2SgBA2UyhG0JmNpLiswK1CGLfQ0WtZZl/HcD9duZ+u5JJHuEwUDFmSBIRwwHA6EU8U77GG3RynHnEMpR230cn2WCQBhDZJAkI4ZDBtGUxbp8MQIU5uAYQQIoRJEhBCiBAmSUAIIUKYJAEhhAhhkgSEECKEBXzZCKVUCrDL6Ths1AAsVhGFNvmZWJOfiyf5mRRqqrUusUZ5wCeBUKOU0lprWaJZhPxMrMnPxZP8TMpOuoOEECKESRIQQogQJkkg8PzL6QACkPxMrMnPxZP8TMpIxgSEECKEyZOAEEKEMEkCQggRwiQJCCFECJMkIIQQIUySgBBChDBJAkIIEcIkCQQwpdRCpZS2+Ip2OjY7KKXaKKUWKKXSlVL7lVJjlVLhTsflFKXULV5+H0Y4HZtdlFItlVITlVJ/KKVylVILLc5RSqnHlFJ7lFIZSqnFSqkODoQbFGR7ycD3I/CY27FMJwKxk1IqAZgPbAAGAC2A/2BuXJ5wMLRAcCmQUeT1dqcCcUBb4CrgVyDSyzmjgSeBh4BNwIPAfKVUO631QVuiDCKSBALfUa31r04H4YARQAxwrdY6DZinlKoGjFFKvZR/LFQt11qfdDoIh8zWWs8EUEpNBWoXbcx/Sh4NPK+1fjP/2C/ATuAe5AbCg3QHiUDVF5jrdrH/DJMYejgTknCa1jqvhFMuAqoBXxT5nlPAbMzvlHAjSSDwXZ7fJ56ulJqrlDrP6YBskoR5lD9Da70bSM9vC2XblFI5Sqk/lVJ3Oh1MgEkCcoEtbsc3Ir83lqQ7KLAtAj4AtgJNgceBJUqp9lrrnU4GZoME4JjF8dT8tlB0ANPXvQwIB4YCbymlYrXWrzoaWeBIAE5qrXPdjqcCsUqpKlrrLAfiCliSBGyklKoO1C/pPK31pvz/Pl3k8BKl1HzM3fH9+V8ihGit5wJzixz6Nr8P/Aml1Oul6CoRwoMkAXsNBiaV4jzLnZG01geVUj8BnXwaVWBKBapbHE/IbxPGVOB6IJHQmiXkTSoQp5QKd3saSADS5SnAk4wJ2EhrPVlrrUr6Kult8r8qu0249eEqpRoDsbiNFYQ47fbfULcJ01XW0u24xxiTMCQJBBGlVD3gEmCl07HY4FvgCqVUfJFjQzDz4xc5E1JAug44DOxyOpAA8TOQhnnqBkApFQv0w/xOCTfSHRSg8mcBPQ98ifkH3gR4FMgDXnMwNLu8BYwCpiulXgSaA2OAV0J1jYBSahpmUPgPzN3ukPyvUaEyHpB/Qb8q/2VDoJpS6rr813O01ulKqReAJ5VSqRQuFgsD3rA94CAgSSBwHcGMDTwP1AJOAAuBgflTJSs1rXWqUqo38CZmjvcx4FVMIghVfwLDgcaY340NwM1a648cjcpedTE3RkUVvG6GWRT2Auai/yjm384K4DKt9SGbYgwqsr2kEEKEMBkTEEKIECZJQAghQpgkASGECGGSBIQQIoRJEhBCiBAmSUAIIUKYJAEhhAhhkgSEECKE/X8BAqnGkD2IMQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# initial parameters\n", "niter = 1000\n", "α = 0.01\n", "β = np.zeros(p+1)\n", "\n", "# call gradient descent\n", "β = gd(X, y, β, α, niter)\n", "\n", "# assign labels to points based on prediction\n", "y_pred = logistic(X @ β)\n", "labels = y_pred > 0.5\n", "\n", "# calculate separating plane\n", "sep = (-β[0] - β[1] * X)/β[2]\n", "\n", "plt.scatter(X[:, 1], X[:, 2], c=labels, cmap='winter')\n", "plt.plot(X, sep, 'r-')\n", "pass" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**1**. Rewrite the `logistic` function so it only makes one `np.exp` call. Compare the time of both versions with the input x given below using the `@timeit` magic. (10 points)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "np.random.seed(123)\n", "n = int(1e7)\n", "x = np.random.normal(0, 1, n)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**2**. (20 points) Use `numba` to compile the gradient descent function. \n", "\n", "- Use the `@vectorize` decorator to create a ufunc version of the logistic function and call this `logistic_numba_cpu` with function signatures of `float64(float64)`. Create another function called `logistic_numba_parallel` by giving an extra argument to the decorator of `target=parallel` (5 points)\n", "- For each function, check that the answers are the same as with the original logistic function using `np.testing`. Use `%timeit` to compare the three logistic functions (5 points)\n", "- Now use `@jit` to create a JIT_compiled version of the `logistic` and `gd` functions, calling them `logistic_numba` and `gd_numba`. Provide appropriate function signatures to the decorator in each case. (5 points)\n", "- Compare the two gradient descent functions `gd` and `gd_numba` for correctness and performance. (5 points)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**3**. (30 points) Use `cython` to compile the gradient descent function. \n", "\n", "- Cythonize the logistic function as `logistic_cython`. Use the `--annotate (-a)` argument to the `cython` magic function to find slow regions. Compare accuracy and performance. The final performance should be comparable to the `numba` cpu version. (10 points)\n", "- Now cythonize the gd function as `gd_cython`. This function should use of the cythonized `logistic_cython` as a C function call. Compare accuracy and performance. The final performance should be comparable to the `numba` cpu version. (20 points)\n", "\n", "Hints: \n", "\n", "- Give static types to all variables\n", "- Know how to use `def`, `cdef` and `cpdef`\n", "- Use Typed MemoryViews\n", "- Find out how to transpose a Typed MemoryView to store the transpose of X\n", "- Typed MemoryVeiws are not `numpy` arrays - you often have to write explicit loops to operate on them\n", "- Use the cython boundscheck, wraparound, and cdivision operators" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**4**. (40 points) Wrapping modules in C++.\n", "\n", "Rewrite the `logistic` and `gd` functions in C++, using `pybind11` to create Python wrappers. Compare accuracy and performance as usual. Replicate the plotted example using the C++ wrapped functions for `logistic` and `gd`\n", "\n", "- Writing a vectorized `logistic` function callable from both C++ and Python (10 points)\n", "- Writing the `gd` function callable from Python (25 points)\n", "- Checking accuracy, benchmarking and creating diagnostic plots (5 points)\n", "\n", "Hints:\n", "\n", "- Use the C++ `Eigen` library to do vector and matrix operations (include path is `../notebooks/eigen3`)\n", "- When calling the exponential function, you have to use `exp(m.array())` instead of `exp(m)` if you use an Eigen dynamic template.\n", "- Use `cppimport` to simplify the wrapping for Python\n", "- See [`pybind11` docs](http://pybind11.readthedocs.io/en/latest/index.html)\n", "- See my [examples](http://people.duke.edu/~ccc14/cspy/18G_C++_Python_pybind11.html#) for help" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "\n", "\n", "\n" ] } ], "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.5" }, "latex_envs": { "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 0 } }, "nbformat": 4, "nbformat_minor": 2 }