LogisticRegression_Circle.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from sklearn.linear_model import LogisticRegression "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dfZxVVb348c93hgFGRAQRH0AuhNzwiUAn0fxl3IxUVETKh67di2WSt8xrlL/kQoqWqelFb9furaiMytLRYEQ0vWShXVN/gsCMz6KiMhoqBj4hDMP398c+Rw/DOfvhnP18vu/Xa15z5uw9+6w9c87+7rXWd60lqooxxhhTSUPSBTDGGJNuFiiMMca4skBhjDHGlQUKY4wxrixQGGOMcWWBwhhjjCsLFMYUiMi/ichPw97Xx7FURPYP41jGREFsHIXJIxE5C/gGMAp4E1gEzFLVjUmWqxwRUWC0qq4ps20ZcATQBSjwDHALcK2qbqn1+Mb4YTUKkzsi8g3gKuBCYADOhfbvgKUi0rvC7/SKr4SBnaeq/YF9cILfGcCdIiLJFsvUCwsUJldEZDfgUuBrqnqXqnap6lrgNJxg8fnCfnNF5FYR+bWIvAmcVXju1yXH+mcReUFENojIt0VkrYh8quT3f114PKLQfDRdRF4UkddFZHbJcQ4XkQdEZKOIvCIi11cKWG5U9R1VXQZMAY4ETvA6vojcV/j11SLytoicLiIDRWSJiLwmIn8rPB4WtDymfligMHnzMaAvsLD0SVV9G/g9MKnk6ZOBW4HdgRtL9xeRA4H/As7EuZMfAAz1eO3/A3wYOAa4WEQOKDzfDXwdGIxzgT8G+ErA8yo9lxeB5cDHvY6vqkcX9vmIqu6qqjfjfO5vwAmcw4HNwPXVlsfknwUKkzeDgddVdVuZba8Uthc9oKptqrpdVTf32PezwO2q+r+quhW4GKePwM2lqrpZVVcDq4GPAKjqClV9UFW3FWo3PwY+EfzUdvAyMKia46vqBlX9naq+q6pvAZeHUB6TY2lulzWmGq8Dg0WkV5lgsU9he9FLLsfZt3S7qr4rIhs8XvuvJY/fBXYFEJG/B+YBLcAuOJ+7FR7H8jIU+Es1xxeRXYBrgeOAgYWn+4tIo6p211guk0NWozB58wCwBZhW+qSI9AOOB+4pedqthvAK8H67vYg0A3tUWab/Bp7EyTzaDfg3oOqOaBHZDzgM+HOVx/8GThPZhML+xeYp6xw3ZVmgMLmiqptwOrP/U0SOE5EmERmBk1K6DviVz0PdCpwkIh8rdAxfSvUX0v44Kbpvi8gY4F+qOYiI7CIinwBuA/4fcKfP468HPtSjPJuBjSIyCLikmvKY+mGBwuSOqn4f5676GpwL6EM4zUjH+B17oKqPAV8DbsKpXbwFvIpTWwnqm8A/Fo4xH7g54O9fLyJv4VzwrwN+Bxynqtt9Hn8usKCQFXVa4RjNOM1wDwJ3BT0hU19swJ0xPojIrsBGnOad55MujzFxshqFMRWIyEmF5p5+OLWTDmBtsqUyJn4WKIyp7GScNNSXgdHAGWpVcFOHrOnJGGOMK6tRGGOMcZXLAXeDBw/WESNGJF0MY4zJjBUrVryuqnuW25bLQDFixAiWL1+edDGMMSYzROSFStus6ckYY4wrCxTGGGNcWaAwxhjjygKFMcYYVxYojDHGuLJAYYwxxpUFCmOMMa4SHUchIj8HTgReVdWDy2yfiDP3fnG2zoWqell8JTTGw5KZsPxntR1DGkC3gzSCdn/wfcB+cMzFMPa0cMpqTJWSHnD3C5xF3X/pss+fVfXEeIpj6lZ7K/z+W7D5jfhfu7isRHEV0uL3TS/BwnOcr1LSCP2GwNuvfPBcQ2+Y+kMLKiYSiQYKVb2vsPqYMfFYMhNW/KJw594AvZqh652kSxWMdu8YJAC2by0fVMrp3Q9OvM6CivEt6RqFH0eKyGqcqZ6/WVh5bCciMgOYATB8+PAYi2dSo70V7rnMuRPvqdicU0q3Zy9IhGHrO5VrKoedBSfOS6RYJr0Sn2a8UKNYUqGPYjdgu6q+LSKTgf9Q1dFex2xpaVGb66mOLJgCz9+bdCnyo6kfdL0LA4ZZH0kdEZEVqtpSbluqaxSq+mbJ4ztF5L9EZLCqvp5kuUwC2lthyQXO3bCJVrGWVamPxPpD6k6qA4WI7A2sV1UVkcNx0nk3JFwsE4f2VrjtPOjeknRJordT1lPDBx3caVSpP2TkJ2D64mTKZCKVdHrsb4GJwGARWQdcAjQBqOqPgM8C/yIi24DN2FKU+fV+/8I6aNole30HTf3gpBA7iL36W3pmPSFAwh+N5++FuQPC/1uYxCXeRxEF66PIgCUzYcUNyd45l8t6ynKH7pKZsPznJB4wikprSln9m9YRtz4KCxQmPu2tcPsFydQWbBCb9+BAaSjEmIiDt9U4UskChUlGkoPYAAaPgfMeSua1syqJ/1nL2VbbSAELFCZe7a2w8Fyg23PXUFgnaviSqP3Vc20vBSxQmHjE0kbeCNN+ZBeTpMQVQKyWETsLFCYa10+A158M/7gNvaHPrrD5bzboK80qZWaVGwVfC6sxxsIChQnfNWN2nm+oVtIAh33B7iTzIIrR8lbLiJQFClOb9lZo+6oz0CpslgGTX25jQWrVPAiOv8reNyGyQGGqF8WdoQWH+hV2c6XNhBsaCxQmuLA/0NbObIoi6RAXaPmiNU3VILOTApqYhTla2kbjmkrGnvZBDSCMFQIBUOc4G9bYDUkErEZhwpuATxrglB9bM4CpTVjNndbEGYjVKMzOwuxotIFSJkzFGkGt43K63vlghlt7b9bEahT1Jsz2Yet3MHEIa1oR6/h2ZTUK42hvhdu+Ct01pLlaLruJW2mfRi03OsUlYF980N7DAVmgqAdhdFLbBHsmDYpBo5axPcXOcwsWvlnTU57V2kltS16atKu1r806vN9n4yjqTa3rS9uoV5NFtWRLWcCwPoq6Uu2HxfoeTNZNX1x9x3fXO7BwhvO4joNFJVajyItq14Cw2oPJo1qaXet0sKg1PeVdNbUIaYRTbF0Hk3O1NMPW2QBSt0DREHdhTIiWzIS5A4IHiYbeFiRMfRh7GvzbyzBtvtMPEYRud9Jpl8yMpmwZYjWKLKq2mck67Ey9q7YPb9d94JsRLNKVIlajyJMFUwrTEgQIEg29nTuq2S9bkDD1bfpi57PQ2DvY7739ilN7r9PaRaKBQkR+LiKvisijFbaLiPxARNaISLuIHBp3GVNlyczgd0ODx8DFr1mAMKZo7Gnw7decTD9pDPa7y38Glw5yavV1JNGmJxE5Gngb+KWqHlxm+2Tga8BkYALwH6o6weu4uWt6qnYqZkt5DUXbyk5m3ryKasa1N0BVv+dm9+Ym5k45iKnjh4Z85DpWzdK+OZvrLNVZTyIyAlhSIVD8GFimqr8t/PwUMFFVXf+juQoU1bSp2nQbgbWt7GT2og7e2Rqw3yeljho1iBvPOTLpYmRLNYt15ShYZHnA3VCgdGz+usJzAUN/BlU7+ZnVIjzNaevg1w++mHQxInX/s28w4qI73v+5X+9GLj/lEKuFuDnvIedzt+jL/udFe/5e52YuJ8GikrTXKO4ArlDV/y38fA/wf1V1RZl9ZwAzAIYPH37YCy+8EGWxo1XNLK85urMJS9vKTr5+86pqVzPIvb369+ah2ZOSLkY6XTEctmzyv38OPn/W9JQ1V430PwWBpbzuZE5bBzc+9CI5zPyOlPV99FBN32CGA0aWA8UJwHl80Jn9A1U93OuYmQwU1cyCmeE3ZVgmzVvGM69+0Dy3W59G3tySj36GNLC+DoL3E2a0jzC1gUJEfgtMBAYD64FLgCYAVf2RiAhwPXAc8C7wBVX1jACZCxTtrXD7+dC12f/v1HmQyEo/QxRZT0kYuEsTl5xUx7WNoM3BGewrTG2giEqmAkXQ6m0G34BhmNPWwW8feonuBN+vTQ1w9anjEr1Ytq3sZO7ix9i4uSuR1xfgzCOG892phyTy+okLUrvIWLOwBYq0ClSlFZj2k8y86cKQRF9DAzDv9GSDQS3aVnZy6e2P8bd34wskn6+3wBHkc9vYG07OxuJfFijSKEhNoqkZTvpBJt5sYenZ9xA2EThzQn1c4M6c/wD3PxtwfYYq1FUWVXtrYSodHwbsB18vO/lEqligSJugb7JjLs59kGhb2cmshe1s7gqnRb9RZIdmqqG7N3PhsR/ObE0hTHPaOrjxwRcjSxuum4AR5GYvA59jCxRpEuTNNTdAHndGRdEpXXdNITWKKjHgugw34fkWtI8xxUkoFijSIsibKsVvqDBEdXGyIFG7KPqGcv9/CfLZTumU5RYo0iDIPDIZzcP2o21lJ9/6XTtbtoXTxNSrQbjm1I/k/841IW0rO7nwllWE0SKY+2lEgqzXncIbQQsUSfvO3tDtc4zEtPmpbsesVpgXnNFD+rF05sTaD2R8a1vZyTdaV9Ed0uUi9zWMuQO890lZqrsFiiQFSaXLYZAIK+9/l6YGvjdtbH7vRjMkzJl2cxsw5u4OftIFUtR6YIEiSX7uLCCVVdFatK3s5IKbV9V8nNxeSHJkzOw7ea+GqkYaBjKGLkifRUpqFhYokhCkvTKlnVvVqnUMRN2P/s2gOW0d/OahF9lew+UkdzcFfhdDkka4JPpxLp7FsEARs/ZWaPsKbPfR3JKiqmetapnWu0+vBq76jDUt5cGEy5ey/q0AU+SXyF1Krd/pylNws2iBIm6X7+tvwaGUVDnDUEstInd3koa2lZ1cffdTdG4MMNFlQaMIn5uwX37eE34zHvsMgFnJTXRpgSJOftsmcxIkaslmam5q4ArroM69aqcQaRT499NyUsPwGywS7Ku0QBEXX0FCoOWLuQgS1dYirAZRf2qZomW3Po20X3pcBKWKmd8MyISChQWKOLS3wsIZuKbESQOc8uNcpMCOveSuqhYIsiBhqqlhCPD8lSdEU6C4+cmETCBV3i1QNMRakjy7/QI886ZzECTOnP8AIy66I3CQGD2kH2uvPMGChOHGc47kqFGDAv2OAiMuuoO2lZ3RFCpOu+7jvc+SC6IvRwAWKMKwYIp353XL2ZkPEpPmLQt8J9jU4GSy2EhqU+rGc47kutPHMXT35kC/d8HNq5jT1hFRqWLyzSedjms3W99xMqZSwpqeauUnVzoHHddBB9DZWssmiKD9XbmYxsVPn0WM6fPW9BSV6yfURZCY09YRKEiMHtLPgoQJZOnMiYGao5559R0mXL40whLFYPpip+Pajd+JRCNmgaIWfv6JGQ8SY2bfGWg68Fzc6ZlEFJuj/Fr/1tZ8BAvEfZ/v7B1LUdxYoKjWgine+zT1i74cEQoyh89Rowax9soTLEiYmkwdP5S1V55A30aPi2fB+re2sv+sOyIuVcRavui+vXuz03qRIAsU1Whv9ZcPfdJ10ZclIpPmLfMdJPbq39uamkyonrx8Mnv17+1r323qZERllp9Wh9efdMZpJcQCRTVuO897n8FjMpvlNOHypb47FutmfWQTu4dmTwrUb5HpZqiWs733WXFD9OWowAJFUO2t0L3FfZ9d98nsRH9+J3T7/BHDWXvlCRYkTKSK/RZ+GqLWv7U1u+Ms/NQqNJxVIauRaKAQkeNE5CkRWSMiF5XZfpaIvCYiqwpfX0qinDv4/be898nolOFjL7nLV5Do2yg2cM7EZur4oTx/5Qm+mqIuuHkVZ85/IIZSRWDafO99/PSNRiCxQCEijcAPgeOBA4HPiciBZXa9WVXHFb5+Gmshe2pv9V5fwk8VMoX8TsnRt1F48vLJMZTImB09NHsSvXxULe5/9o1sBouxp3lfP56/17kOxSzJGsXhwBpVfU5VtwI3AScnWB5vC2d47NCYyXTYM+c/4CtIjB7Sz4KESdSaK/zN93T/s29kcwT3ifO8g8Xt8U/vkWSgGAq8VPLzusJzPX1GRNpF5FYR2a/SwURkhogsF5Hlr732WthlLWQceGQBTftR+K8bsTltHb6m5dirf29LfTWpsNZnM9SvH3wxm30WJ85zJhCtxM9aNyFLMlCUq0T2vBLfDoxQ1bHAH4AFlQ6mqj9R1RZVbdlzzz1DLGaB1/ThzYMyl+XUtrLT12A6y2wyafPQ7Em+Bud9o7X2ddsTcdgX3LdfMyaechQkGSjWAaU1hGHAy6U7qOoGVS2mGM0HDoupbDvy04F0/FXRlyNkfqbl2K1PowUJk0pTxw/1TJ/tVqf/LXO8mrDffiXWQXhJBoqHgdEiMlJEegNnADus1iEipfPxTgGeiLF8H/AzcVfGahN+Pzy5WDDG5NaN5xxJk8dV7M0t3dkcve3VVxHjPFCJBQpV3QacB9yNEwBaVfUxEblMRIq38OeLyGMisho4Hzgr9oJ6Rm3J3JiJOW0dvjqvg8y7Y0xSrj7V+326TcleJtSJ87ynAYqpVpHoOApVvVNV/15VR6nq5YXnLlbVxYXHs1T1IFX9iKr+g6rGP0DBK2pP+0k85QiJ336J607PyVrFJvemjh/K6CHe86pVs2534rymAYqpVmEjs2shjZlrcpq1sN1zHwsSJmuWzpzoK1hkLmXWz/UlhnEVFijceP0DTslWOmzbyk5fi9tbkDBZtHTmRHbr0+i6TyZTZr36KvzMFlEjCxRuvNatzVhtYqaPLCfrlzBZ5if5IsgiXKnglQHlNVtECCxQVNLe6qxbmxNtKzvxqkvs1b+31SZM5vmZcTZzHdteIp6C3AJFJV7VOa8lDFPG6y7KBtWZvPCzNkrmOra9rjcrfhHpy1ugqMS1OtdYWMIwG/zcPVmQMHny+SOGe+6TqY5tr+uNeqe718ICRTUyNqeT191TkMVhjMmC7049xHOm2SBrwaeCV60iwuYnCxTleGU7ZagT20+Ghy1javLIz0yzmeqr8KpVeM1HVwMLFOW4TePbnK27b69xE36q6MZklVdtOXN9FQmxQFGO2zS+GZr8z8+4CVupzuRZ7mrLCSXR9ErkVdPMq50vQ81OV9/9lOv2erhL6OrqYt26dbz33ntJFyV1+vbty7Bhw2hqakq6KJHq06uBLdsq3zCdOf+B7ASU6Yth7oDYX9YCRU8RtvPFrXPjZtft8+pgcN26devo378/I0aMQMTHOpp1QlXZsGED69atY+TIkUkXJ1JXfWasa3q4NT95q4ebyvB4zeSYIn46sethcN17773HHnvsYUGiBxFhjz32qIua1tTxQz2nIs8Ut37SiDKf8vTni57XTI4pMnfxY67b/SwlmRcWJMqrp7+Ln6nIM8OtnzSiFhELFEFkqH9i4+Yu1+02wC4+f/3rXznjjDMYNWoUBx54IJMnT+bpp59m7dq1HHzwwZG85pYtWzj99NPZf//9mTBhAmvXro3kdbLCq/Y8ad6yeAoShgSuQxYoSkU8X0paNNTPjWTiVJVTTjmFiRMn8uyzz/L444/zve99j/Xr10f6uj/72c8YOHAga9as4etf/zrf+lb0M4xm2TOv5mdetyhYoCiVo45sN/84wcZOVNK2spOjrvwjIy+6g6Ou/GPNU1L/6U9/oqmpiXPPPff958aNG8fHP/7xHfZbu3YtH//4xzn00EM59NBD+ctf/gLAK6+8wtFHH824ceM4+OCD+fOf/0x3dzdnnXUWBx98MIcccgjXXnvtTq972223MX36dAA++9nPcs8996CqNZ1L1u3enKfsLpe7vQhueC1Q+DV4TNIl8M1rtKmNnSivbWUnsxZ20LlxM4qTNTZrYUdNweLRRx/lsMMO89xvyJAhLF26lEceeYSbb76Z888/H4Df/OY3HHvssaxatYrVq1czbtw4Vq1aRWdnJ48++igdHR184Qtf2Ol4nZ2d7LfffgD06tWLAQMGsGHDhqrPIw/mTjko6SKEp+WLlbct/3noL2eBwq8MrYtt6X7Vufrup9jctePkapu7uj3Ho4Shq6uLc845h0MOOYRTTz2Vxx9/HICPfvSj3HDDDcydO5eOjg769+/Phz70IZ577jm+9rWvcdddd7HbbrvtdLxytYd66rwux6ufIlPTebiuURF+zdEChTEFL1cYd1LpeT8OOuggVqxY4bnftddey1577cXq1atZvnw5W7duBeDoo4/mvvvuY+jQofzTP/0Tv/zlLxk4cCCrV69m4sSJ/PCHP+RLX/rSTscbNmwYL730EgDbtm1j06ZNDBqUreln4mY3WJVZoKgzzblKKA/Xvrs3B3rej09+8pNs2bKF+fPnv//cww8/zL333rvDfps2bWKfffahoaGBX/3qV3R3OzWbF154gSFDhnDOOedw9tln88gjj/D666+zfft2PvOZz/Cd73yHRx55ZKfXnTJlCgsWLADg1ltv5ZOf/GTd1yjA5jarll016swV08YmXYTUuvDYD9PctOOay81NjVx47IerPqaIsGjRIpYuXcqoUaM46KCDmDt3Lvvuu+8O+33lK19hwYIFHHHEETz99NP06+cM7ly2bBnjxo1j/Pjx/O53v+Nf//Vf6ezsZOLEiYwbN46zzjqLK664YqfXPfvss9mwYQP7778/8+bN48orr6z6HPLE+ueqY1N41Jl6GI1dreLf5uq7n+LljZvZd/dmLjz2wzX/zfbdd19aW8tPXf/oo48CMHr0aNrbP5jpt3jxnz59+vvZS6XK1SJK9e3bl1tuuaXaIhuzg0QDhYgcB/wH0Aj8VFWv7LG9D/BL4DBgA3C6qq6Nu5ymfkwdP9SCqTE9JNb0JCKNwA+B44EDgc+JyIE9djsb+Juq7g9cC2Rnjm9jjMkJz0AhIueJyMAIXvtwYI2qPqeqW4GbgJN77HMysKDw+FbgGLEeOWOMiZWfGsXewMMi0ioix4V4oR4KvFTy87rCc2X3UdVtwCZgj3IHE5EZIrJcRJa/9tprIRXRGGOMZ6BQ1TnAaOBnwFnAMyLyPREZVeNrlws4PUeK+NnHeVL1J6raoqote+65Z41FM8YYU+Srj0KdYZ5/LXxtAwYCt4rI92t47XXAfiU/DwNerrSPiPQCBgA2KsYYY2Lkp4/ifBFZAXwfuB84RFX/BScT6TM1vPbDwGgRGSkivYEzgMU99lkMFHMDPwv8Uet9ZjOTOUlMM37fffdx6KGH0qtXL2699dZIXsPUDz/psYOBaar6QumTqrpdRE6s9oVVdZuInAfcjZMe+3NVfUxELgOWq+pinOauX4nIGpyaxBnVvp4xSShOMz59+nRuuukmAFatWsX69evfn7QvCsOHD+cXv/gF11xzTWSvYeqHnz6Ki3sGiZJtT9Ty4qp6p6r+vaqOUtXLS15vceHxe6p6qqrur6qHq+pztbye8bdEal1rb4VrD4a5uzvf28sPlPMrqWnGR4wYwdixY2losMkXTO1sZHYpaQDdXv75nLj67qdsQFkl7a1w+/nQVZgEcNNLzs9Q9apiQacZ79u3L8888wyf+9znWL58+fvTjM+ePZvu7m7efffdHaYZB9i4cWNVZatHmVrJLkUsUJSSXqBbd34+R90inTXMhJp791z2QZAo6trsPB/x8pNdXV2cd955rFq1isbGRp5++mnAmWb8i1/8Il1dXUydOpVx48btMM34CSecwKc//elIy5YntpJddfJzqxyG7WWCBABaN8uk1rVN64I970NS04yb4DI1kjfm65EFCr8ytEyqTaVcpQHDgj3vQ1LTjJvgrj19XNJF8M9tFTtprLytShYoSjXnY2EXr6mU57R1xFSSjDnmYmjqsfZEU7PzfJWSmmb84YcfZtiwYdxyyy18+ctf5qCDcrQMaESy1Xfn0hx+2Fmhv5rkcVhCS0uLLl++PPgvtrfCwnMqb5+7qfpCxWzERXe4bl975QkxlSRZTzzxBAcccID/X2hvdfokNq1zahLHXBx5/0SSAv99MqxtZScX3Lyq4vZMfSbmDnDZVt11SkRWqGpLuW3WmV1q7GnugcLk39jTch0Y6tnXXYJEptSYsl0Na3oKYsGUpEvg2+7NTUkXwZjUOHP+A26NNRw1KkPNzovO9d4nZBYognj+Xu99UmLuFPc2aeunMPXk/mfdp4i78ZwjYypJCLQ79pe0QJFTXh1zv37wxZhKkrw89sOFoV7+LnU1G0FECTkWKHpqOTvpEoRm4C7uzU/18AHq27cvGzZsqJuLol+qyoYNG+jbt2/SRYncpbc/5rp9l6YcXQaPj2YRUOvM7unEeZkaM+HmkpMOcs3ymNm6KmMpgcENGzaMdevWYYtZ7axv374MG1b9GJGs+Nu7Xa7bvzdtbEwlCcE1Y9y3R5SIYYEiqPbWzGTFTB0/1DVQbFenky9T7bMBNTU1MXLkyKSLYRIy4fKlnvtk6mbp7VcSedkc1blikrH02d36uI/S9OrkMyar5rR1sP6tStPyOPI1i0H4I7KLLFCU49VPkaF5n9ovPc5zH5tR0+SRn4QNr1kMUsUrPX/ajyJ7aQsU5Zw4z337ihviKUdIvO6annn1nbro2Db1w8/7OVNjJ8A7PT/CJnELFNUot2ZFivm5a/LKDDEmSy68xXsUdq765iKYCLCUBYpKdt3HfXsCw+hr4dVX4ZUZYkxWzGnroMvjXm70kH7xFCYsXtebCCYCLGWBopJvPum+/fffiqccIfHTV+E1kaAxWeCnb2LpzInRFyRMC2e4b/dqLq+RBQo30+ZX3rb5jczVKq7zMd/+2EvuiqEkxkTjzPkPeO6TuUynBVNwnVY8BhYo3Iw9zX1I/G1fja8sIZg6fqjnKl5vbum2LCiTSXPaOjzTvY8aNShbmU7g3Yk98hORF8EChRe3IfHdWzM1oyzAmT7upp559R2bNNBkStvKTl9NTrnqwC6avjjyl7BA4cWrVpGhGWXByYDq5WNx4HqaNNBkm9eCREWZS4cFuH6C+/aY5qZLJFCIyCARWSoizxS+D6ywX7eIrCp8RR82K/GaaCtjfRVrrvC3ktf+s6xz26SfnyAxeki/7NUm2lvhdY+kmog7sYuSqlFcBNyjqqOBewo/l7NZVccVvpJr4/EayLLoy/GUI0R+ln3cpta5bdLN7/szc1lO4J3pFONM10kFipOBBYXHC4CpCZUjHLrdu4qYQn6yoN7c0u0rk8SYuE24fClvbvFexCeTTSAOE6sAABL6SURBVE5LZuKZ6RRTbQKSCxR7qeorAIXvQyrs11dElovIgyLiGkxEZEZh3+WRTCndx2Uxc/CuIqbQ1PFD2at/b8/97n/2DcuEMqkyad4yzwn/AHpJRjuwvaYJiiHTqVRkgUJE/iAij5b5OjnAYYaragvwj8B1IjKq0o6q+hNVbVHVlj333LPm8u9klo/O3YxlQAE8NHuSr87tZ159hzGz74y+QMZ4aFvZyTOvvuNrX7/9camyZKb3NEExZDqViixQqOqnVPXgMl+3AetFZB+AwvdXKxzj5cL354BlwPioyuuL2wA8cDKgMtaxDc6HyU+weK9brc/CJKptZSezFvpL3fbTD5c67a3eC6fFXJuA5JqeFgPTC4+nA7f13EFEBopIn8LjwcBRwOOxlbAcP7Mz3n5B9OWIwJorTvCcDwqsz8IkZ9K8ZVxw8yo2d7n3S+zWpzGbQQJgicf1Y9d9Yq9NQHKB4kpgkog8A0wq/IyItIjITwv7HAAsF5HVwJ+AK1U12UAB3pkGXe9kslYBznxQfjr+7n/2DRuQZ2Izp62DERfd4au5aa/+vX3Na5ZKS2bCVo9z9JqDLiKSx0XnW1padPny5dG9wPUT3DuvG/vAt8u2pmXCpHnLfH0o+zYKT14+OYYSmXrl970IzliJTKbBFs3dHddMp+ZB8K3nI3t5EVlR6BPeiY3MrsZ5D0GjS7ZQ95bM1irAf875e93KiIvusEWPTCTmtHX4ChLNTY1cd/q4bAeJ6yfgng7b4D3wN0IWKKp18g/dt2dsGvKe/IyxKLrg5lXWFGVCNWneMt/TyFwx7RCmjh8acYkitGCKd3r9tB9HuoKdFwsU1fL6p21+I5PpskVTxw8NNB3zrx980YKFqVnbyk7f/RHgNDdlOkgsmek9X1zL2YkGCbBAURuvju3n74Xv7B1PWSLw3amHcN3p49ilyd/bxIKFqcWctg5f8zYV9W2UbDc3+UmFhVhHYFdigaIWJ87zzmnu3gyXDo6nPBGYOn4oj3/nePo2+hhogRMsRlx0h43kNoHMaesINGNxLyH7iRQLz/HeJ8b5nNxY1lMYrhrpNDW5GfmJRPKfwzRm9p281x3s/XLUqEHZnELBxKKa99RufRqzmwJbdM0YePsV930Gj3ESZ2JiWU9R85ONkNFR26WevHxy4GUk73/2DcuMMmWNuOiOQEGiUZwki8wHiQVTvINEQ+9Yg4QXq1GEZcEU706pxt5OtlTCHVO1alvZyexFHbyz1XvmzlK9G4Xvf/Yj2e58NDU5c/4DnsuVlpOLWgQ4N4t+mpymzY/9OmE1ijhMXwyNze77dG/NfNosOP0Wj112HJ8/YrjnGtyltnarpdLWsbGX3FVVkDhq1KB8BAnw9/kfPCZ1N5MWKML07b+CNLnvs/mNwlzz2ffdqYdw7enj8NnP/b5fP/giB37799YcVSeKU3D4WTui1C5NDVx3+rj89HEtmendl9lnQKqanIqs6SkKfpqhWs5ORdpbWIJMtVDU1ChcbU1RueV3LetycpUE0d4KbV+F7R7rZ+y6T2JzOYF705MFiqj4yYSKOashDmMvuSvwnSPk7MJQ56q5aSiVq/dCeyssOhfU4zPR0BsujmDBtQCsjyIJx18FDR7NUK8/mcklVN20X3oco4f0C/x7xeyoUbPutD6MjGpb2cnIAKOqe+rVIPlqagJY6CdINMFUjymBEmY1iij5zXCYuyn6ssSsbWUnF96yii6PhboqEeDMI4bz3amHhFouE75qs+CKcpPRVMrvZx8SyXAqx5qekrRkpvcw/QH7wTEXp+LNEra2lZ3828J23q02YgC7Nzcxd8pB1peRItWmufb0+TzeDPj5zBelqK/SAkXS/LxxmprhpB/kMliAEzDmLn6MjZu7aj7WXv1789DsSSGUygQRVnCAnAYICFaTSNlsDRYo0sBrsaOiHNcuisK64DQ3NXDFtLFW04jQnLYObnzwRdeVEoLIdZAPEiRSVJMoskCRFn6DBaTubiNsYd6dAgzcpYlLTrLmqTDMaevgtw+9RHeI14bMrz7nJUhzU0qzHS1QpM21B8Oml7z3y3mwKAo6c6gX69MIJuygXVQ3Nb4gNYk+A2BWeO/1MFmgSJv2Vrj9fOja7L1vnQQLiO6CZXNM7ahtZSdX3/0UnRt9vP+qUHeB2s+YKUh1kAALFOnU3gr3XOavZpHyN1iY2lZ2cuntj/G3d2vv9HbT1ABXnzquLi5mYdfYymkA5p1eH3/P9wX5DKe0uamUBYo0811tFZj2k1x3cvcUxwWuVB7a0dtWdjJrYTuba0hHDmLo7s1ceOyH6ytAQC5bBSxQpN139nRmlvXS0AhTf1RXwQLiDxg9pfFiWNp8JBBaVpIfjQL/flqd1R6KgiSkQGaCBKQwUIjIqcBc4ADgcFUte1UXkeOA/wAagZ+q6pV+jp+5QNHeCgtn4PvjnqE3X9jivmOuVpCR5V7NbcVA0NzUkOh5i8CZE3I6/sFL0M8opDIF1k0aA8UBwHbgx8A3ywUKEWkEngYmAeuAh4HPqerjXsfPXKAA541423nQvcXf/hlo84xarVNHGHf9ejdy+SmH1GfNoVSQ1FfI7Fgot0DRK+7CAKjqEwAirgsZHA6sUdXnCvveBJwMeAaKTBp7mvN1xXDY4mPup9efdLItjr8qc2/IsEwdP3SHi1hcHeF5ZnNs9RAkSOR4doVEAoVPQ4HSdIJ1QL6mWi1n1ov+g8XmN5zq8IsPZqqKG5WegQPCH1mcJzZI0UOQIJHRWoRfkQUKEfkDsHeZTbNV9TY/hyjzXMXPu4jMAGYADB8+3FcZU2vWi4U36c/xbhNV58288le5WI87bN+desgOd8dRjdXIgroZAFcrPwuPlaqDZuBEs55EZBmV+yiOBOaq6rGFn2cBqOoVXsfNZB9FJUHftJC5TrSkRTFlRRrkagGguFwzBt5+xf/+Ca9KF6bUdWa//+LugaIXTmf2MUAnTmf2P6rqY17HzVWggOApeVDXmVFhSUvto1LWk/UnhGTJTFjxC+8FhnrKUZCAFAYKETkF+E9gT2AjsEpVjxWRfXHSYCcX9psMXIeTHvtzVb3cz/FzFyjAyYr6/bf8TRVQqg6qxWlSzTKgdsFPUDU1dsjl5yp1gSJquQwURUFT9QAaesHU/7b+C2OKqvkcIdDyxdw269qa2Xly4jynWSmI7ducaUKuGunUTIypV+2tcOmgKm62ejtT6OQ0SHixQJFF0xc7HdZlE8NcbH7DAoapX9dPcN7/1fRFXPxaXdfILVBk1YnzYO7GQsAIqBgwFkwJv1zGpE17qzOfWtCEEHBq7znqsK6WBYqsO3EeTJsPTf2C/+7z91qwMPnV3gqX7encFPmZdLPUgP2cz5VlDgLWmZ0v7a2w8FygyrmPmvrBSdfVdRXb5ER7Kyw6N3gzUx2PQbKsp3rU3gqLvgxaxWyjNgbDZE3P9HFpCP7er+MgASmcFNDEoFgruO2rwavdz98LcwfU/QfHZES5sRBBgoTVpD1ZoMiz4hu/moF64KQQLv+Z1TBM+gRZhrQSe1/7Zk1P9SLoehflSAMc9gWrZZhk7BAcaljXzwJEWTbgzji1i2+/6mRyNA+q7hi63alhXLanjcMw8Smmty48p6QG4TdIlIw1kganOdWCRGBWo6hX7a3V9V+UsrZdE7VqJsQsNW2+vT99sqwnU1l7K7R9FbbXEDBKWbXehKGquZh6sPdiIBYojLcwA0afAc7iS8b4EUbHdE+WsReYBQrjX62D9sppHlTXa3ubMpbMhBU3VDfOp6ecL0MaFwsUJrhq17/wYkHDBF1FrhJrWgqVBQpTm/ZWWHIBbA22II+3fM/vX/eiaFIqsqal0FmgMOGIqpZR1LsfnGhZVJlW7YpxflktIjIWKEz4og4aOVxqMtei6NsqsoGesbBAYaITdnptKVvCNX1qHdcQhI3TiZVNCmiiM/a0Dz7IYU2xULR9G9x+Prz4IKy+CboKfSR2hxmfqGuOpSzRIbWsRmGiFeWFZuQn4IUHdq7N2AWnektmwopfBF/HoRbW75AK1vRk0iGy7CmfLFMm+f+BNMJhZ9n/IYUsUJj0ibJvI6g8toW/3wy4Dhp71zZrcK2shpcJFihMennm2ofQ1xFEYx//F9U4aig7NAV5/S0K25sHwZa3YHtXtGVzY8Ehc1IXKETkVGAucABwuKqWvaqLyFrgLZycu22VTqInCxQZVnonPGAYjP40LL8BCGGqh9gJ9N7FaeaRRudi3zwItm35oGM+LyzBIPPSmPX0KDAN+LGPff9BVV+PuDwmLUqzqIqGHwG3X7Bj1tOIj8Pz9xFrbSMw/aAvoNg5HEf2UFysz6duJBIoVPUJABHx2tWY8sED0tXPkXc2ar6upX0chQL/IyIK/FhVf1JpRxGZAcwAGD58eEzFM4kqDSBJpHXmkaWqmjIi66MQkT8Ae5fZNFtVbyvsswz4pksfxb6q+rKIDAGWAl9T1fu8Xtv6KMwOop5/KI0ae0PvXWHz33bOemrsAydfb7UDs4NE+ihU9VMhHOPlwvdXRWQRcDjgGSiM2YHXHXKaUklrUsh6svUZTMhS2/QkIv2ABlV9q/D408BlCRfL5FGlPpByAo00DznryQarmYQkEihE5BTgP4E9gTtEZJWqHisi+wI/VdXJwF7AokKHdy/gN6p6VxLlNeZ9QYKKMTmRVNbTImBRmedfBiYXHj8HfCTmohljjOmhIekCGGOMSTcLFMYYY1xZoDDGGOPKAoUxxhhXFiiMMca4skBhjDHGVS7XoxCR14AXajzMYCAPs9baeaSLnUe62Hl84O9Udc9yG3IZKMIgIsv9rn+RZnYe6WLnkS52Hv5Y05MxxhhXFiiMMca4skBRWcW1LzLGziNd7DzSxc7DB+ujMMYY48pqFMYYY1xZoDDGGOPKAkWBiJwqIo+JyHYRqZhmJiJrRaRDRFaJSOrWWw1wHseJyFMiskZELoqzjH6IyCARWSoizxS+D6ywX3fhf7FKRFKz2LPX31dE+ojIzYXtD4nIiPhL6c3HeZwlIq+V/A++lEQ53YjIz0XkVRF5tMJ2EZEfFM6xXUQOjbuMfvg4j4kisqnkf3FxaC+uqvbl9NMcAHwYWAa0uOy3FhicdHlrOQ+gEXgW+BDQG1gNHJh02XuU8fvARYXHFwFXVdjv7aTLWs3fF/gK8KPC4zOAm5Mud5XncRZwfdJl9TiPo4FDgUcrbJ8M/B5nLdkjgIeSLnOV5zERWBLFa1uNokBVn1DVp5IuR618nsfhwBpVfU5VtwI3ASdHX7pATgYWFB4vAKYmWJag/Px9S8/vVuAYKSznmCJZeJ94UtX7ALe1a08GfqmOB4HdRWSfeErnn4/ziIwFiuAU+B8RWSEiM5IuTJWGAi+V/Lyu8Fya7KWqrwAUvg+psF9fEVkuIg+KSFqCiZ+/7/v7qOo2YBOwRyyl88/v++QzhSabW0Vkv3iKFqosfB78OlJEVovI70XkoLAOmshSqEkRkT8Ae5fZNFtVb/N5mKNU9WURGQIsFZEnC5E+NiGcR7k719jzpN3OI8Bhhhf+Hx8C/igiHar6bDglrJqfv28q/gce/JTxduC3qrpFRM7FqSV9MvKShSsL/ws/HsGZr+ltEZkMtAGjwzhwXQUKVf1UCMd4ufD9VRFZhFM9jzVQhHAe64DSO79hwMs1HjMwt/MQkfUiso+qvlJoBni1wjGK/4/nRGQZMB6nXT1Jfv6+xX3WiUgvYAAJNSu48DwPVd1Q8uN84KoYyhW2VHweaqWqb5Y8vlNE/ktEBqtqzZMeWtNTACLST0T6Fx8DnwbKZiCk3MPAaBEZKSK9cTpTU5MxVLAYmF54PB3YqaYkIgNFpE/h8WDgKODx2EpYmZ+/b+n5fRb4oxZ6JFPE8zx6tOVPAZ6IsXxhWQz8cyH76QhgU7HZM0tEZO9iP5eIHI5zfd/g/ls+Jd2Tn5Yv4BScO4stwHrg7sLz+wJ3Fh5/CCfzYzXwGE5TT+JlD3oehZ8nA0/j3H2n8Tz2AO4Bnil8H1R4vgX4aeHxx4COwv+jAzg76XK7/X2By4Aphcd9gVuANcD/Az6UdJmrPI8rCp+F1cCfgDFJl7nMOfwWeAXoKnw2zgbOBc4tbBfgh4Vz7MAl6zHl53Feyf/iQeBjYb22TeFhjDHGlTU9GWOMcWWBwhhjjCsLFMYYY1xZoDDGGOPKAoUxxhhXFiiMMca4skBhjDHGlQUKYyImIh8tTJrXtzC6/zEROTjpchnjlw24MyYGIvJdnNHYzcA6Vb0i4SIZ45sFCmNiUJgr6WHgPZypFboTLpIxvlnTkzHxGATsCvTHqVkYkxlWozAmBoX1vG8CRgL7qOp5CRfJGN/qaj0KY5IgIv8MbFPV34hII/AXEfmkqv4x6bIZ44fVKIwxxriyPgpjjDGuLFAYY4xxZYHCGGOMKwsUxhhjXFmgMMYY48oChTHGGFcWKIwxxrj6/8LhkSOwt53pAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X.shape for the data: (2000, 2)\n",
      "y.shape for the labels: (2000, 1)\n"
     ]
    }
   ],
   "source": [
    "innerRadiusTrue = 1\n",
    "bias = 0.5\n",
    "outerRadiusTrue = innerRadiusTrue+bias\n",
    "Npoints = 1000\n",
    "angles = 2*np.pi*np.random.random(Npoints)\n",
    "xInner = innerRadiusTrue*np.cos(angles).reshape((Npoints,1))\n",
    "yInner = innerRadiusTrue*np.sin(angles).reshape((Npoints,1))\n",
    "plt.scatter(xInner,yInner,label=\"Class 0\")\n",
    "angles = 2*np.pi*np.random.random(Npoints) #Random angles in [0,2pi]\n",
    "xOuter = outerRadiusTrue*np.cos(angles).reshape((Npoints,1))\n",
    "yOuter = outerRadiusTrue*np.sin(angles).reshape((Npoints,1))\n",
    "plt.scatter(xOuter,yOuter,label=\"Class 1\")\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.title(\"Original Data\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "xInnerSq = xInner**2\n",
    "xOuterSq = xOuter**2\n",
    "Xinner = np.concatenate((xInnerSq,yInner),axis=1)\n",
    "Xouter = np.concatenate((xOuter,yOuter),axis=1)\n",
    "X = np.concatenate((Xinner,Xouter),axis=0)\n",
    "print(f\"X.shape for the data: {X.shape}\")\n",
    "y = np.concatenate((np.zeros((Npoints,1)),np.ones((Npoints,1))),axis=0)\n",
    "print(f\"y.shape for the labels: {y.shape}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigma(z):\n",
    "    '''Computes the sigmoid for a numpy array z\n",
    "    '''\n",
    "    return 1/(1+np.exp(-z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def computeGradientNonLinear(x,y,z):\n",
    "    '''Computes the gradient for the cross entropy loss function. x,y are (n,) numPy arrays, signifying the 2 features.\n",
    "       z =  theta0+theta1x+theta2y is an (n,) numPy array. Returns the (3,) gradient at the current iteration.\n",
    "    '''\n",
    "    n = x.shape[0]\n",
    "    sigma_z = sigma(z)\n",
    "#     print(f\"sigma_z: {sigma_z}\")\n",
    "#     print(f\"sigma_z.shape: {sigma_z.shape}\")\n",
    "#     estimation_error = (y/sigma_z)-((1-y)/(1-sigma_z))\n",
    "    estimation_error = sigma_z-y\n",
    "    term1 = sum(sigma_z*(1-sigma_z)*estimation_error)\n",
    "    term2 = sum(x*sigma_z*(1-sigma_z)*estimation_error)\n",
    "    term3 = sum(y*sigma_z*(1-sigma_z)*estimation_error)\n",
    "    return (np.array([term1,term2,term3]))/n\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def GradientDescentNonLinear(X,y,start,step_size,tolerance,Niter):    \n",
    "    ''' Will perform gradient descent on X (nx2), using target y (nx1) \n",
    "    '''    \n",
    "    gradNorm = 10000\n",
    "    i = 0\n",
    "    current_theta_estimate = start\n",
    "    norms = []\n",
    "    n = X.shape[0]            \n",
    "    while i<Niter and gradNorm > tolerance:        \n",
    "        z_estimate = current_theta_estimate[0]+current_theta_estimate[1]*X[:,0] + current_theta_estimate[2]*X[:,1]\n",
    "#         print(f\"z_estimate.shape = {z_estimate.shape}\")\n",
    "#         print(f\"z_estimate = {z_estimate}\")\n",
    "        grad = computeGradientNonLinear(X[:,0],X[:,1],z_estimate)\n",
    "#         print(f\"grad = {grad}\")\n",
    "        gradNorm = np.linalg.norm(grad)\n",
    "        norms.append(gradNorm)        \n",
    "        current_theta_estimate = current_theta_estimate - step_size*grad\n",
    "        i+=1\n",
    "    plt.plot(np.array(list(range(1,i+1))),norms)\n",
    "    plt.xlabel('Iterations')\n",
    "    plt.ylabel('||gradient||')\n",
    "    plt.grid(True)\n",
    "    plt.title(\"Convergence of gradient descent\")\n",
    "    print(f\"Started at location {start}\")\n",
    "    print(f\"Finished with {i} iterations\")\n",
    "    print(f\"Finished with gradient norm = {gradNorm}\")\n",
    "    print(f\"Finished at location {current_theta_estimate}\")\n",
    "    return current_theta_estimate\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def VanillaLogisticRegressionNonLinear(X,y):    \n",
    "    '''Assumes that the input matrix is an nx2 numPy array with n observations and y is an nx1 target vector \n",
    "    '''\n",
    "    n = X.shape[0]\n",
    "#     theta_estimate_start = np.random.rand(3)\n",
    "    theta_estimate_start = np.array([0,1,1])\n",
    "    theta_estimate_start = np.array([0,1.5,2])\n",
    "    learning_rate = 1e-6\n",
    "    err_tolerance = 1e-6\n",
    "    MaxIter = 8000\n",
    "    estimated_coeff = GradientDescentNonLinear(X,y,theta_estimate_start,learning_rate,err_tolerance,MaxIter)\n",
    "    return estimated_coeff\n",
    "\n",
    "   \n",
    "\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Started at location [0.  1.5 2. ]\n",
      "Finished with 8000 iterations\n",
      "Finished with gradient norm = 0.11734195571689643\n",
      "Finished at location [-6.02490666e-04  1.49947988e+00  2.00049799e+00]\n",
      "theta estimated [ theta0,theta1, theta2]: [-6.02490666e-04  1.49947988e+00  2.00049799e+00]\n",
      "Radius estimated : 1.4994800038539466\n",
      "Original Radius : 1.25\n",
      "[0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1]\n",
      "Sum error: 0.0\n",
      "Error probability : 0.0\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEWCAYAAADPZygPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3gU1frA8e+bShEpASO9SEACIgiCCAIqXQVUVNDrBUWxIWpuEX567XjFey+ICir2DoiIEVGKEKpUpSMQioKg9CbS398fc6LruptkQzab8n6eZx9mz5xz5p3dZd+dmZMzoqoYY4wxkRAV6QCMMcYUXZaEjDHGRIwlIWOMMRFjScgYY0zEWBIyxhgTMZaEjDHGRIwlIWMKAPG8KSJ7RWRhhGJQEantll8WkX9FKI4+IjInEts2uc+SkDltInKjiCwWkUMisl1EvhCRVpGOq5BpBbQHqqhqs0gHo6p3quqTp9uPiLQVka25EVN+ICJpInJbpOMoSCwJmdMiIinAc8DTQCJQDRgJdItkXL5EJCbSMeSC6sBmVf3ldDtyR1X2f9/kD6pqD3vk6AGUBg4B12VSJx4vSW1zj+eAeLeuLbAV+BuwA9gO3OLWXQT8BET79HU1sNwtRwEDgQ3AbmAsUM6tqwEo0Bf4AZjlyv8KfO/q/wvYDLQLob/err9dwEM+cUUD/+faHgSWAFXdunOBqcAeYC1wfSavVSUg1dVNB2535X2BI8BJ93o/HqBtNPA/F9smoL+LOcatTwMGA3OBX4HawC3AGhfzRuAOvz7/4d6TbcCtrr/abt1bwFM+da8ElgL7gHlAQ591m4G/A8uB/cAYoBhQ0sVyyu3XIaBSgH1LcK/LAWAh8CQwx2d90NcY6AKsdvv4I/B3n3XdXMwH3HvXyedz/brb9x+Bp3CfQ6APMAf4L7DXvdad3brB7j064vblxUj/Hy0Ij4gHYI+C+wA6AScyvuiC1HkCmA+cBVRwX1BPunVtXfsngFj3hXEYKOvWbwDa+/T1ETDQLd/v+q2Cl+heAT5062q4L8x33BddcSDZfTG0AuLcl8hxfk9C2envVdfX+cBRoJ5b/w9gBVAXELc+wW17C96XfQxwAV6SqB/ktZqJdxRZDGgE7AQud+v6+H7xBmh7p/uyrQKUBabx5yT0A1DfxRILXAGc42Ju4177C3ze25+BBm4/PiBIEnL7tQNojpcMe+MlnowfG5vxkkcloBxe4rvT5zOwNYvP2Wi8HwUlXTw/ZrwWWb3GeInkErdc1mf/muElxPZ4P0AqA+e6dRPc+18S73O7EJeg3ftwHLjd7etdeElafF7n2yL9f7MgPSIeQGF7ANcBq/B+3TXNpN4b7j/uSr/yMXi/zpa6/7xLXXkzn/JlwNW5EGtr4Bu8RNAjB+1vAn7Kos4GoIvP8454p5UyvoB+xSeJudfkIrf8FPCGWy4F/AJUd8/X4L6g3fOK7sshht+TRi2f9Y/gkop7XgI4xu9JKDv9VfFZvxDo6ZbXAt0C7PsNwGy/sleARwPUrYr3K7qUT9m/gbfcch8yT0LT8TmSAdrx5yT0RBbv1QTgPp/P5zM+6+oQPAm9hPth4VN/LdDGLW8G/uKz7lngZZ/PQNAkhPdFfxyXIFzZ0/yehDJ9jfES7x3AmQHqDAuwvUS8HxjFfcp6ATN83od0v8+RAmf7vM6WhEJ42Hnh0+Auqr7lV7wSuAaYlUXzt/B+bf6Bqt6gqo1UtRHwMTDep9+mrrwT8EouXOv4Ae8/1Qc5bL8bKJ9FHJXwToFl+N6V/daHqp7weX4YOMMtfwBcIyLxeK/pN6qa0Vd14BMR2Sci+/CSyEm8L5EMW/zi+O25qh528WfITn8/BYmzKl6y9VcdaJ7Rp+v3JuDsAHUrAXtU9aBP2fd4v9Cz4w/757ccsExEOovIfBHZ42LrApQP0p/ve+ivOvA3v/2syh/f52CvXVYq4P0QCBZLVq/xtXj79b2IzBSRFq48s/csFtju098reEdEf9oX9zkihP0xfgrDBdt8RVXXAIhIVvVmiUiNYOvF6+B64DJX/7DP6mJ4v74y6v4FGIB3mmkBcLeqnsxGrJtd+1NZ1Q3ia7zz392BcUHqbMP7j73KPa/myrKkqqtF5HugM3Ajf0yWW4BbVXWufzuf11V9irfjnS7LqFMc75RZKP0FswXvtNbKAOUzVbV9Fu3Be03KiUgpn0RUDe/UU3ZsxzsVl6FqgDq+n5l4vB85fwU+VdXjIjIB79RcRn++fVTLZNtbgMGqOjibsQaMKYideEfqVYHvAsSS6WusqouAbiISi3edbKzrK+M987cF70iovN+Po+zKan+MHzsSyr8uAX5W1fUZBSLSXERW4V1/uFNVT4hIPbxTEi3dUdJJvF+CYaeq+/FOc40Qke4iUkJEYt0v7GddtQ+Bh0WkgoiUd/XfC2EzH+Al2NZ414QyvAwMFpHqAK7/zEbkjQOuEpGLRSQOeJzfv3Bz0p+v14AnRSTJjTxrKCIJwESgjojc7F6XWBG50L1nf6CqW/Cul/1bRIqJSEO8AQnvZzOGscB9IlJZRMoAD2ZRPw7v2tdO4ISIdAY6+PXXR0SSRaQE8Ggmfb0K3Ok+nyIiJUXkChEplY24fwYSRKR0oJXux9R44DH3+UrGu+aUIehrLCJxInKTiJRW1eN4AxAyfpy9DtwiIpeLSJR73c5V1e3AFOB/InKmW3eOiLTJxr5k7E+tbNY1WBLKERFZICJL8b58uorIUvfomIub6YX3Bf4bVV2gqvWBC4FBIlIMuBxoAixyMV2O+08gIu+IyMoAj7tzK0hVHQqkAA/jfaFtwfvFOcFVeQpYjDcyagXeNainQtjEh3jXDaar6i6f8uF4I6amiMhBvEEFzTOJcxVwL95F7u14o6V24P3qDbk/P0PxvrSn4H3RvY53TeEg3hd7T7wjnZ+AIXhf/oH0wrv+tA34BO+6xtRsxvCq2/5y4FtgEt4RRMAjYhfbABf3XrwjzVSf9V/gjWScjjdSb3qwDavqYrwL9S+6vtLxTvNmSVW/w3uPN7rTX5UCVOuPd7rrJ7zT2G/67Udmr/HNwGYROYA3eOMvrt1CvMEMw/AGKMzEO2IH7+gwDm+gx168HzAVs7M/eJ+jHuL9UfHz2WxTpGWM6DA5ICJtgT6q2ifAujS84aCLM2lfA5ioqg38ymPwTsM0UdWAf8gnIjPwRmW1wBvWOihHO+H19ZaLI9gptUJHRM7AG06cpKqbIh1PbnNHNi+ravUsKxsTQXYklD+1A77zTUAiUtMlJ9wpo7p4o46+wvvldZZbVy7jlJL5IxG5yp3SKYk3RHsF3mtY4IlIcRHpIiIxIlIZ7/TZJ5GOy5isWBLKZSJytXjTkLQAPheRya68kohM8qn3Id6F/boislVE+vp00xO/U3F4f9+yzJ1y+wRv8MEuVV2Ndypsiogsx/ujvWydOnDnzrfiDSt/xV1vKsy68fsfzSbhDbEuLKcCBO86116803Fr8K6/GZOv2ek4Y4wxEWNHQsYYYyLG/k4oROXLl9caNWrkqO0vv/xCyZIlczegXGBxhcbiCo3FFZr8GhecXmxLlizZpaoV/rQi0lM2FLRHkyZNNKdmzJiR47bhZHGFxuIKjcUVmvwal+rpxQYsVpu2xxhjTH5iScgYY0zEWBIyxhgTMZaEjDHGRIwlIWOMMRFjScgYY0zEhD0JiUgnEVkrIukiMjDA+ngRGePWLxCfe7eIyCBXvtZ3hupgfbr51RaIyHrXZ5wrryYiM0TkWxFZLiJdstqGMcaY8AtrEhKRaGAE3k3JkoFe7n4gvvoCe1W1Nt606kNc22S8OdTq491JdKSIRGfR5xC8W/Ym4c2hlTEf28PAWFVt7Pocmdk2cvdV8MzfuJvJm49z8pRNk2SMMRnCfSTUDO9+7BtV9RjevVz8bxTWDXjbLY8DLhcRceWjVfWoelPtp7v+Avbp2lzG73f4fBvvjp/g3e3wTLdcmt/v7BlsG7nu8+Xb+fC7Y/R4eR7rfz6YdQNjjCkCwj1tT2X+eG/4rfz5RmG/1VHvTqH78W67XBnvxmK+bSu75UB9JgD79Pdb8vrWfwxvlul7gZJ4t0rI2HawbfxGRPoB/QASExNJS0sLtr9BXVZaiaqjfLxpH52fm8VV58RyRa1YYqIyvw14Xjh06FCO9incLK7QWFyhsbhCF47Ywp2EAn3D+p+PClYnWHmgo7fM6oN3x8q3VPV/ItICeFdEGmQzPlR1FDAKoGnTptq2bdsAzbImksa917fg8c9W88mybaw5VIxnezSkYZUyOeovt6SlpZHTfQoniys0FldoLK7QhSO2cJ+O2wpU9Xlehd9Phf2pjrtpW2lgTyZtg5XvAspk3PjNb1t98W5jjKp+DRQDymczvlxV/ox4XujVmFf/2pS9h4/RfcRc/j1pDb8eC3gXZmOMKdTCnYQWAUlu1Foc3iCAVL86qUBvt9wDmO4mu0sFerrRczXxbkK2MFifrs0M1weuz0/d8g/A5QAiUg8vCe3MZBth1z45kSkPtOGGC6vyyqyNdB4+i/kbd+fFpo0xJt8IaxJy12f6A5Px7vQ4VlVXicgTItLVVXsdSBCRdCAFGOjarsI7elkNfAnco6ong/Xp+noQSHF9Jbi+Af4G3C4iy/DuWNrHTewacBvhej38lS4ey7+vacgHtzXnlELPUfN56JMVHDxyPK9CMMaYiAr7/YRUdRIwya/sEZ/lI3i3lw7UdjAwODt9uvKNBBjdpt4tsFuGso28dHHt8nx5/yUMnbKON+ZuYvp3Oxh8dQMuOzcxkmEZY0zY2YwJ+USJuBgevjKZj++6mFLFYrj1rcXcP/pb9vxyLNKhGWNM2FgSymcaVyvLxHsv4b7Lk/h8xXbaDZ1J6rJteJe8jDGmcLEklA/FxUTxQPs6fHZvK6qWLc6AD7/l9ncW89P+I5EOzRhjcpUloXzs3LPPZPzdLXmoSz3mpO+i/dCZfLjwBzsqMsYUGpaE8rnoKOH21rX48r7W1K98JoPGr+DGVxfw/e5fIh2aMcacNktCBUSN8iX54LaL+Pc157Hyx/10fG4Wr83eaBOiGmMKNEtCBUhUlNCrWTWmprShVe3yPPX5Gq55aR5rf7IJUY0xBZMloQLo7NLFePWvTXm+V2O27DnMlS/MZtjUdRw7cSrSoRljTEgsCRVQIkLX8ysxLaUNXc6ryPCv1nPlC7NZumVfpEMzxphssyRUwJUrGcfwno15vXdTDvx6gmtGzuWpiattQlRjTIFgSaiQuLxeIlNSWtOzWTVem7OJjs/NYt6GXZEOyxhjMmVJqBA5s1gsT199Hh/efhFRAje+uoBB45dzwCZENcbkU5aECqEW5yTwxX2tuaN1LcYs2kL7oTOZuvrnSIdljDF/YkmokCoeF82gLvWYcE9LypaI4/Z3FtP/g2/YdehopEMzxpjfWBIq5BpWKUNq/1aktK/D5FU/0X7oTCZ8+6NN/WOMyRcsCRUBcTFRDLg8ic8HXEL1hJLcP2Ypfd9ezLZ9v0Y6NGNMEWdJqAipk1iKj++6mH9dmczXG3bTYdgs3pv/PafsqMgYEyFhTUIi0klE1opIuogMDLA+XkTGuPULRKSGz7pBrnytiHTMqk8Rqen6WO/6jHPlw0RkqXusE5F9Pm1O+qxLDdfrkJ9ERwl9W9Vk8v2tOb9qaR6esJIhC4+waZdNiGqMyXthS0IiEg2MADoDyUAvEUn2q9YX2KuqtYFhwBDXNhnoCdQHOgEjRSQ6iz6HAMNUNQnY6/pGVR9Q1Uaq2gh4ARjvs/1fM9apatdcfgnytWoJJXivb3OevbYhPxw8RafnZvHKzA2cOGlT/xhj8k44j4SaAemqulFVjwGjgW5+dboBb7vlccDlIiKufLSqHlXVTUC66y9gn67NZa4PXJ/dA8TUC/gw1/awgBMRrr+wKk+3Kk7rOhX49xffcfXIeazediDSoRljioiYMPZdGdji83wr0DxYHVU9ISL7gQRXPt+vbWW3HKjPBGCfqp4IUB8AEakO1ASm+xQXE5HFwAngGVWdEGhHRKQf0A8gMTGRtLS0wHuchUOHDuW4bTjFnjjMjVWFpLh43lu9n6temE2XWrF0PSeW2CiJWFz59fWyuEJjcYUmv8YF4YktnEko0LeX/xXwYHWClQc6csusvq+ewDhV9Z1UrZqqbhORWsB0EVmhqhv+1JHqKGAUQNOmTbVt27YBNpe1tLQ0cto2nDLiuhS4/ZdjPDlxNeO//ZE1B+MZcm1DmlQvG9G48huLKzQWV2jya1wQntjCeTpuK1DV53kVYFuwOiISA5QG9mTSNlj5LqCM6yPYtnridypOVbe5fzcCaUDj7O5cYVW2ZBxDb2jEm7dcyOGjJ+jx8jwe/2wVh4+dyLqxMcaEKJxJaBGQ5EatxeElAf8RaKlAb7fcA5iu3l9RpgI93ei5mkASsDBYn67NDNcHrs9PMzYiInWBssDXPmVlRSTeLZcHWgKrc23vC7hL657FlJQ23HxRdd6cu5kOw2YxZ71NiGqMyV1hS0Lu+kx/YDKwBhirqqtE5AkRyRiJ9jqQICLpQAow0LVdBYzFSwpfAveo6slgfbq+HgRSXF8Jru8MvfAGOvieoqsHLBaRZXgJ7BlVtSTk44z4GJ7o1oCxd7QgNjqKv7y+gH+OW8b+wzYhqjEmd4TzmhCqOgmY5Ff2iM/yEeC6IG0HA4Oz06cr34g3ei5QX48FKJsHnJfpDhgAmtUsxxf3XcLwr9YzatZGZqzdyZPdGtCpwdmRDs0YU8DZjAkmW4rFRvNgp3OZcHdLyp8Rz53vLeGe979h50GbENUYk3OWhExIzqtSmtT+LflHx7pMXf0z7YbO5OMlW21CVGNMjlgSMiGLjY7inktrM+m+S6h91hn87aNl9HlzET/ahKjGmBBZEjI5VvusM/jojhY8dlUyizbvocPQmbzz9WZOnbKjImNM9lgSMqclKkro09KbEPWC6mV55NNV3DDqazbsPBTp0IwxBYAlIZMrqpYrwTu3NuM/PRqy9qeDdB4+m5Fp6Ry3CVGNMZmwJGRyjYhwXdOqTPtbGy6rexbPfrmW7iPmsvLH/ZEOzRiTT1kSMrnurFLFePnmJrx00wX8fOAo3UbM5T+Tv+PI8ZNZNzbGFCmWhEzYdD6vItNSWnN148qMmLGBLs/PZvHmPZEOyxiTj1gSMmFVpkQc/73ufN65tRlHj5/iule+5rHUVfxy1CZENcZYEjJ5pHWdCkx5oDW9W9Tg7a+9CVFnrtsZ6bCMMRFmScjkmZLxMTzWtT4f3dGC+Ngoer+xkL+NXca+w8ciHZoxJkIsCZk817RGOSYNuIT+l9ZmwtIfaTd0Fl+s2B7psIwxEWBJyEREsdho/t6xLqn9W5J4Zjx3vf8Nd767hB0HjkQ6NGNMHrIkZCKqfqXSfHpPSx7sdC7T1+6g3dCZfLR4i02IakwRYUnIRFxMdBR3tT2HL+67hLpnl+If45bz1zcWsvOwzbZgTGFnScjkG+dUOIMx/VrwZLf6fPP9Xh6e+ytvzt3ESZsQ1ZhCK6xJSEQ6ichaEUkXkYEB1seLyBi3foGI1PBZN8iVrxWRjln1KSI1XR/rXZ9xrnyYiCx1j3Uiss+nTW9Xf72I9A7X62CyLypKuLlFDSY/0Jo6ZaN5/LPVXP/K16TvOBjp0IwxYRC2JCQi0cAIoDOQDPQSkWS/an2BvapaGxgGDHFtk4GeQH2gEzBSRKKz6HMIMExVk4C9rm9U9QFVbaSqjYAXgPFuG+WAR4HmeLcFf1REyub+K2FyokrZEqQ0iWfo9eezYechugyfw4vT19uEqMYUMuE8EmoGpKvqRlU9BowGuvnV6Qa87ZbHAZeLiLjy0ap6VFU3Aemuv4B9ujaXuT5wfXYPEFMv4EO33BGYqqp7VHUvMBUv4Zl8QkS45oIqTH2gDe3rJ/LfKevo+qJNiGpMYRITxr4rA1t8nm/FO+oIWEdVT4jIfiDBlc/3a1vZLQfqMwHYp6onAtQHQESqAzWB6ZnE94c2Pm37Af0AEhMTSUtLC1QtS4cOHcpx23AqCHFdVwlqRcfzzuqDdH1xDp1qxNK9dixx0RLRuPITiys0FlfowhFbOJNQoG8H/yvMweoEKw905JZZfV89gXGqmjGVc3baeIWqo4BRAE2bNtW2bdsGqpaltLQ0cto2nApKXG2B2w4f5+lJaxizeAvfHYzjmWsb0qxmuYjGlV9YXKGxuEIXjtjCeTpuK1DV53kVYFuwOiISA5QG9mTSNlj5LqCM6yPYtnry+6m47MZn8pnSJWIZ0qMh7/VtzrGTp7j+la/514SVHDxyPNKhGWNyIJxJaBGQ5EatxeElgVS/OqlAxqi0HsB09f5KMRXo6UbP1QSSgIXB+nRtZrg+cH1+mrEREakLlAW+9tn2ZKCDiJR1AxI6uDJTALRKKs+UB1pza8uavLfgezoOm8WMtTsiHZYxJkRhS0Lu+kx/vC/2NcBYVV0lIk+ISFdX7XUgQUTSgRRgoGu7ChgLrAa+BO5R1ZPB+nR9PQikuL4SXN8ZeuENdPjtdJuq7gGexEtsi4AnXJkpIErExfDIVcmMu/NiSsbHcMubi0gZs5S9v9iEqMYUFOG8JoSqTgIm+ZU94rN8BLguSNvBwODs9OnKN+KNngvU12NByt8A3gi6A6ZAaFK9LBMHtGLE9HRGpm1g5rqdPN6tPlecVxFv4KQxJr+yGRNMoRAfE01Kh7p8dm8rKpUpTv8PvqXfu0v42SZENSZfsyRkCpV6Fc/kk7svZlDnc5m1bifths5kzKIfbEJUY/IpS0Km0ImJjuKONufw5f2tqVfxTB78eAU3vbaAH3YfjnRoxhg/loRMoVWzfElG334Rg69uwPKt++n43Cxen2MTohqTn1gSMoVaVJRwU/PqTE1pTYtzEnhy4mqufWke6362CVGNyQ8sCZkioWLp4rzeuynDezbi+92/cMXzs3n+q/UcO2ETohoTSZaETJEhInRrVJlpKW3o1KAiQ6euo+uLc1i2ZV/WjY0xYWFJyBQ5CWfE80Kvxrz616bsPXyMq0fO5elJa/j12MmsGxtjcpUlIVNktU9OZGpKG264sCqjZm2k8/BZfL1hd6TDMqZIsSRkirQzi8Xy72sa8sFtzTml0OvV+fzfJys4YBOiGpMnLAkZA1xcuzyT72/N7ZfUZPTCH+gwdBbTv/s50mEZU+hZEjLGKR4XzUNXJDP+7paULh7LrW8t5r7R37L70NFIh2ZMoWVJyBg/jaqW4bN7W3F/uyQmrdhO+2Gz+HTpjzb1jzFhYEnImADiYqK4v10dJt57CVXLleC+0Uu57e3FbN//a6RDM6ZQsSRkTCbqnl2K8XddzMNX1GPuhl10GDqLtC3HOWVT/xiTKywJGZOF6CjhtktqMfn+1jSoXJq3Vh3jxtfms3nXL5EOzZgCz5KQMdlUPaEkH9zenFvqx7HqxwN0Gj6LV2dttAlRjTkNWd5ZVUQeBbLzvyxNVWf5te0EDAeigddU9Rm/9fHAO0ATYDdwg6pudusGAX2Bk8AAVZ2cWZ8iUhMYDZQDvgFuVtVjbt31wGNuP5ap6o2u/CSwwoXzg6pm3HbcmIBEhDZVY7mj60U8PGEFgyetYeLybTzb43zqnl0q0uEZU+Bk50hoM/B9Nh5/mIBLRKKBEUBnIBnoJSLJfn33Bfaqam1gGDDEtU0GegL1gU7ASBGJzqLPIcAwVU0C9rq+EZEkYBDQUlXrA/f7bP9XVW3kHpaATLadXboYr/61KS/0aszWvb9y5QuzGTZ1HUdP2NQ/xoQiyyMhVX07h303A9JVdSOAiIwGugGrfep0wztCARgHvCgi4spHq+pRYJOIpLv+CNSniKwBLgNudHXedv2+BNwOjFDVvW5/duRwf4z5AxHhqvMr0bJ2eZ74bBXDv1rPFyu3M+TahjSuVjbS4RlTIGTndFy1bPa1T1UP+DyvDGzxeb4VaO7X5rc6qnpCRPYDCa58vl/bym45UJ8JbvsnAtSv4/ZjLt4pvMdU9Uu3rpiILAZOAM+o6oRAOyYi/YB+AImJiaSlpQWqlqVDhw7luG04WVyhCRRX97OhxgXxvLP6F64ZOY8O1WO4JimO+BiJaFz5gcUVmvwaF4QntiyTEN5RhQKZ/W9S4C286zsZAtX3v7YUrE6w8kCnDzOrD94+JgFtgSrAbBFpoKr7gGqquk1EagHTRWSFqm74U0eqo4BRAE2bNtW2bdsG2FzW0tLSyGnbcLK4QhMsrrZA3yPHGfLld7w3/wfWHIzjmWsacHHt8hGNK9IsrtDk17ggPLFl53TcpTnseytQ1ed5FWBbkDpbRSQGKA3syaJtoPJdQBkRiXFHQ771twLzVfU43qm9tXhJaZGqbnP7uFFE0oDGwJ+SkDHZVapYLE91P48rG1Zi4MfLufG1BfS8sCqDutSjdPHYSIdnTL6T7SHaIvJudsp8LAKSRKSmiMThDTRI9auTCvR2yz2A6erNjZIK9BSReDfqLQlYGKxP12aG6wPX56dueQJwqYu3PN7puY0iUtaNzssob8kfr1cZk2MX1Urgy/tbc0ebWoxdvIUOw2YydbVNiGqMv1D+Tqi+7xN35NIkWGV3RNIfmAysAcaq6ioReUJEMkaivQ4kuIEHKcBA13YVMBYvKXwJ3KOqJ4P16fp6EEhxfSW4vnF1d4vIarxE9Q9V3Q3UAxaLyDJX/oyqWhIyuaZYbDSDOtdjwj0tKVsijtvfWUz/D75hl02IasxvsjMwYRDwf0BxEckYeCDAMdx1kmBUdRIwya/sEZ/lI8B1QdoOBgZnp09XvpHfR9D5litegkvxK58HnJdZ/MbkhoZVypDavxWvzNzAC9PTmZO+i0evSqZ7o8p4g0GNKbqyPBJS1X+raingP6p6pnuUUtUEVR2UBzEaU+DFxURx7+VJfD6gFTXLl+SBMcu49a1FbNtnE6Kaoi3bp+NUdZCIVBaRi0WkdcYjnMEZU9gkJZZi3J0X88iVyczfuIcOw2bx7vzvbUJUU2RlZ4g2ACLyDN5AgNV4U+mANwx6VtBGxpg/iY4Sbm1Vk/bJiQwav4J/TVjJZ8i8t9wAACAASURBVMu2MeTahtQsXzLS4RmTp7KdhICrgbpuFgNjzGmqWq4E7/ZtxkeLt/Lk56vp9NwsHmhfh9ta1SQm2uYWNkVDKJ/0jYD9oYMxuUhEuP7CqkxLaUObOhV45ovv6D5yLqu3Hci6sTGFQChJ6DCwVEReEZHnMx7hCsyYoiTxzGK8cnMTRt50AT/tP0LXF+fwvylrbUJUU+iFcjoulT//sakxJpeICF3Oq0iLWgk8+flqXpiezhcrf2LItQ1pUt0mRDWFU7aTkKq+LSLF8eZbWxvGmIwp0sqWjGPo9Y3oen4lHvpkJT1enkefi2vw9w51KRkfyu9GY/K/UKbtuQpYijeDASLSSETsyMiYMGlb9ywmP9Camy+qzptzN9PxuVnMXr8z0mEZk6tCuSb0GN6MBPsAVHUpUDMMMRljnDPiY3iiWwPG3tGCuOgobn59If8ct4z9h49HOjRjckUoSeiEqu73K7O/sDMmDzSrWY5J913CXW3P4eNvfqTdsJl8ufKnSIdlzGkLJQmtFJEbgWgRSRKRF4B5YYrLGOOnWGw0D3Y6l0/vaUmFM+K5870l3P3+EnYcPBLp0IzJsVCS0L14M2kfBT4EDgD3hyMoY0xwDSqX5tP+LflHx7pMW7OD9kNn8fGSrXhz9RpTsIQyOu4w8JB7GGMiKDY6insurU3H+mfz4MfL+dtHy2hQPpqkRoepUrZEpMMzJtuyPBISkefcv5+JSKr/I/whGmOCqX3WGXx0Rwse71qf9XtP0nHYLN75erNNiGoKjOwcCWXcPfW/4QzEGJMzUVFC74trUGLfRlK3l+SRT1eRunQbQ3o05JwKZ0Q6PGMylWUSUtUl7t+Z4Q/HGJNTFUpE8c6tzfj4mx95cuJqOg+fzX2XJ9GvdS1ibUJUk09l53TcChFZHuyRRdtOIrJWRNJFZGCA9fEiMsatXyAiNXzWDXLla0WkY1Z9ikhN18d612ecz7rrRWS1iKwSkQ98ynu7+utFpHdWr4Ux+Z2I0KNJFaamtKZdvbP4z+S1dB8xl5U/+v91hTH5Q3Z+Hl0JXIU3U8KXwE3uMQkYF6yRiEQDI4DOQDLQS0SS/ar1Bfaqam1gGDDEtU3Gu3dRfaATMFJEorPocwgwTFWTgL2ub0QkCRgEtFTV+rgRfSJSDngUaI73R7iPiohN0GUKhbNKFWPkTU14+S8X8POBo3QbMZdnv/yOI8dtQlSTv2Tn9t7fq+r3eF/i/1TVFe4xEOiYSdNmQLqqblTVY8BooJtfnW7A2255HHC5iIgrH62qR1V1E5Du+gvYp2tzGb8nxbeB7m75dmCEqu51+7PDlXcEpqrqHrduKl7CM6bQ6NSgIl+ltOGaxpUZmbaBLs/PZvHmPZEOy5jfhDIbYkkRaaWqcwBE5GIgs9tAVga2+DzfinfUEbCOqp4Qkf1Agiuf79e2slsO1GcCsE9VTwSoX8fFOxeIBh5T1S+DxFeZAESkH9APIDExkbS0tGD7nKlDhw7luG04WVyhKYhxXVEBqjeN582Vh7nu5a+5rFoMPerEUTxGIhpXJFlcoQtHbKEkob7AGyJS2j3fB9yaSf1An27/caPB6gQrD3Tklll98PYxCWgLVAFmi0iDbMbnFaqOAkYBNG3aVNu2bRuoWpbS0tLIadtwsrhCU1DjagvcctUJ/jN5LW9/vZnvDsTy9DUNaFOnQkTjihSLK3ThiC3bQ2ZUdYmqng80BM5X1Uaq+k0mTbYCVX2eVwG2BasjIjFAaWBPJm2Dle8Cyrg+/Le1FfhUVY+7U3tr8ZJSduIzplApGR/DY13rM+7OFhSLjaL3GwtJGbuUfYePRTo0U0SFNG5TRK4A7gTuE5FHROSRTKovApLcqLU4vIEG/n/cmgpkjErrAUxXb+6RVKCnGz1XEy9pLAzWp2szw/WB6/NTtzwBuNTFXx7v9NxGYDLQQUTKugEJHVyZMYVek+rl+HzAJfS/tDapS7fRbuhMJq3YHumwTBEUyv2EXgZuwJtDToDrgOrB6rvrM/3xvtjXAGNVdZWIPCEiXV2114EEEUkHUoCBru0qYCywGm9E3j2qejJYn66vB4EU11eC6xtXd7eIrMZLVP9Q1d2qugd4Ei+xLQKecGXGFAnFYqP5e8e6fNq/JWeXLsbd73/Dne8uYccBmxDV5J1QrgldrKoNRWS5qj4uIv8DxmfWQFUn4Q3l9i17xGf5CF4yC9R2MDA4O3268o14o+f8yxUvwaUEWPcG8EZm+2BMYVe/Umkm3N2SV2dvYti0dcwbuouHr0zmuiZV8AaeGhM+oZyOy/h5dFhEKgHHsZvaGVMoxERHcVfbc/jyvks49+wz+ee45fz1jYVs2XM40qGZQi6UJPSZiJQB/gN8A2zGu6WDMaaQqFXhDEb3u4gnuzfgm+/30mHYLN6cu4mTNiGqCZNsJSERiQK+UtV9qvox3rWgc31PrRljCoeoKOHmi6ozJaUNzWuV4/HPVnPdy/NI33Ew0qGZQihbSUhVTwH/83l+NMCtvo0xhUjlMsV5s8+FDLvhfDbu+oUuw+fw4vT1HD95KtKhmUIklNNxU0TkWrErlcYUGSLC1Y2rMC2lDe3rJ/LfKeu46oU5rNhqv0FN7gglCaUAHwFHReSAiBwUkQNhissYk4+UPyOeETdewCs3N2HPL8foPnIuz3xhE6Ka0xfKjAmlVDVKVeNU9Uz3/MxwBmeMyV861j+bqSlt6HFBFV6euYHOw2ezYOPuSIdlCrBQ/lj1ggCPc3ymyjHGFAGli8cypEdD3r+tOSdOneKGUfN5eMIKDh45HunQTAEUyum4kXgzW7/qHvPxbqWwTkQ6hCE2Y0w+1rJ2eSbf35q+rWry/oIf6DhsFjO+25F1Q2N8hJKENgONVbWJqjYBGgErgXbAs2GIzRiTz5WIi+FfVybz8V0XUzI+hlveWsQDY5ay5xebENVkTyhJ6FyfedpQ1dV4SWlj7odljClILqhWlokDWjHg8iQ+W7aN9kNnMnH5NrxZs4wJLpTrOWtF5CW8U3DgTWa6TkTi8abwMcYUYfEx0aS0r0PnBmfz4MfL6f/Bt3yavI0uZ9nfFZngQjkS6oN3m+37gQfwbofQBy8BXZrbgRljCqZ6Fc9k/F0X839dzmXWup3835xfGb3wBzsqMgGFMkT7V1X9n6perardVfW/qnpYVU+p6qFwBmmMKVhioqPo1/ocJt/fmmqlohg4fgU3vbaAH3bbhKjmj7I8HScijxLkttd+0lR11umHZIwpLGqUL8mDzYrxU4lzeHrSGjo8N5O/d6jLLS1rEh1lk6+Y7F0T2pzNvvadRhzGmEIqSoQbm1fj0nMr8NAnK3nq8zVMXL6dZ3s0pE5iqUiHZyIsO6fjZmTzsdm/oYh0EpG1IpIuIgMDrI8XkTFu/QIRqeGzbpArXysiHbPq093ye4GIrHd9xrnyPiKyU0SWusdtPm1O+pT733rcGJOLKpYuzuu9mzK8ZyN+2HOYK56fzfBp6zl2wgYuFGXZORJ6G+90XGbHzgq8BbyTUSAi0cAIoD2wFVgkIqluaHeGvsBeVa0tIj2BIcANIpIM9ATqA5WAaSJSx7UJ1ucQYJiqjna3Iu8LvOTajFHV/gHi/lVVG2XjNTDG5AIRoVujyrSqXZ7HP1vNsGnrmLTCOyo6v2qZSIdnIiDLJKSqOR351gxIz/g7IhEZDXQDfJNQN+AxtzwOeNHN0t0NGK2qR4FNIpLO77fu/lOfIrIGuAy40dV52/WbkYSMMflIwhnxPN+rMV3Pr8TDE1Zy9ci59G1Vk5T2dSkeFx3p8EweCmWIdqgqA1t8nm91ZQHrqOoJYD+QkEnbYOUJwD7XR6BtXSsiy0VknIhU9SkvJiKLRWS+iHTPwT4aY05Du+REpqS0pmezarw6exOdhs/i6w02IWpREs7JRwOdvvMfZResTrDyQEkzs/oAnwEfqupREbkT7yjpMreumqpuE5FawHQRWaGqG/w7EpF+QD+AxMRE0tLSAmwua4cOHcpx23CyuEJjcYUmO3F1KAtVLyzGm6t+pder82lbJYbr68ZRIjZ8I+gK8usVKWGJTVXD8gBaAJN9ng8CBvnVmQy0cMsxwC68hPKHuhn1gvXp2uwCYgJt26d+NLA/SLxvAT2y2q8mTZpoTs2YMSPHbcPJ4gqNxRWaUOI6fPSEDv58tdYcOFGbD56m01b/lC/iykv5NS7V04sNWKwBvlPDeTpuEZDkRq3F4Q008B+Blgr0dss9gOku2FSgpxs9VxNIAhYG69O1meH6wPX5KYCIVPTZXldgjSsv66YcQkTKAy354/UqY0weKx4Xzf91qcf4u1tSungsfd9ezIAPv2X3oaORDs2ESdhOx6nqCRHpj3cUEw28oaqrROQJvIyYCrwOvOsGHuzBSyq4emPxksIJ4B5VPQkQqE+3yQeB0SLyFPCt6xtggIh0df3swZtqCKAe8IqInMI7zfeM/nHknjEmQhpVLcNn97bipbQNvDhjPXPSd/HoVcl0Pb8S3tglU1iE9YZ0qjoJmORX9ojP8hHguiBtBwODs9OnK9/I7yPofMszTtn5l88DzstyJ4wxEREXE8V97ZLofN7Z/HPccu4bvZTUpdt46uoGVCxdPNLhmVwSztNxxhhz2uokluLjuy7m4SvqMXfDLtoPncX7C77n1CmbELUwsCRkjMn3oqOE2y6pxZT729CwSmke+mQlN742n827fol0aOY0WRIyxhQY1RJK8P5tzXnmmvNY9eMBOj43i1GzNnDipE39U1BZEjLGFCgiQs9m1Zia0oZLkirw9KTvuPaleXz304FIh2ZywJKQMaZAOrt0MV79axNevLExW/f+ypXPz2Ho1HUcPXEy0qGZEFgSMsYUWCLClQ0rMS2lDVedX4nnv1rPlc/P4Zsf9kY6NJNNloSMMQVe2ZJxDLuhEW/2uZBDR09w7UvzeHLiag4fO5F1YxNRloSMMYXGpeeexZQHWnNT82q8PmcTHZ+bxdz0XZEOy2TCkpAxplApVSyWp7qfx5h+FxETFcVNry1g4MfL2f/r8UiHZgKwJGSMKZSa10rgi/su4Y42tRi7eAvth85kyqqfIh2W8WNJyBhTaBWLjWZQ53pMuKcl5UrG0e/dJdzzwTfsPGgTouYXloSMMYVewyrehKh/71CHqat+pv2wmczbdiLjNi4mgiwJGWOKhNjoKPpflsSk+1pRq3xJRi0/yi1vLeLHfb9GOrQizZKQMaZIqX1WKT6682JuOjeOBRv30GHoTN6dbxOiRoolIWNMkRMdJbSvEcuUB1rTuFpZ/jVhJT1HzWfjzkORDq3IsSRkjCmyqpYrwbt9m/Fsj4Z899MBOg+fzcszbULUvGRJyBhTpIkI1zetyrSUNrStW4FnvviO7iPnsnqbTYiaF8KahESkk4isFZF0ERkYYH28iIxx6xeISA2fdYNc+VoR6ZhVnyJS0/Wx3vUZ58r7iMhOEVnqHrf5tOnt6q8Xkd7heh2MMfnfWWcW45Wbm/LSTRfw0/6jdH1xDv+dvJYjx21C1HAKWxISkWhgBNAZSAZ6iUiyX7W+wF5VrQ0MA4a4tslAT6A+0AkYKSLRWfQ5BBimqknAXtd3hjGq2sg9XnPbKAc8CjTHuy34oyJSNldfBGNMgdP5vIpMS2lNt0aVeXFGOlc8P5sl3++JdFiFVjiPhJoB6aq6UVWPAaOBbn51ugFvu+VxwOUiIq58tKoeVdVNQLrrL2Cfrs1lrg9cn92ziK8jMFVV96jqXmAqXsIzxhRxZUrE8b/rz+ftW5tx5Pgperz8NY+lruKXozYham6LCWPflYEtPs+34h11BKyjqidEZD+Q4Mrn+7Wt7JYD9ZkA7FPVEwHqA1wrIq2BdcADqrolSHy+bX4jIv2AfgCJiYmkpaUF3uMsHDp0KMdtw8niCo3FFZqCHtfDTYVx62J4a95mJn77PX3qx9OgfHTE44qEcMQWziQkAcr8B+IHqxOsPNCRW2b1AT4DPlTVoyJyJ95R0mXZjM8rVB0FjAJo2rSptm3bNlC1LKWlpZHTtuFkcYXG4gpNYYircztYtHkPD368nP8u/oXrmlTh4SuSKV0iNqJx5bVwxBbO03Fbgao+z6sA24LVEZEYoDSwJ5O2wcp3AWVcH3/YlqruVtWMiaJeBZqEEJ8xxgBwYY1yTBpwCXe3PYfx3/5Iu2Ez+XLl9kiHVeCFMwktApLcqLU4vIEGqX51UoGMUWk9gOnqTeaUCvR0o+dqAknAwmB9ujYzXB+4Pj8FEJGKPtvrCqxxy5OBDiJS1g1I6ODKjDEmoGKx0fyz07l8ek9LKpwRz53vfcNd7y1hx8EjkQ6twArb6Th3jac/3hd7NPCGqq4SkSeAxaqaCrwOvCsi6XhHQD1d21UiMhZYDZwA7lHVkwCB+nSbfBAYLSJPAd+6vgEGiEhX188eoI/bxh4ReRIvsQE8oao2BMYYk6UGlUvzaf+WjJq1keFfrWfeht3868pkrr2gMt44KZNd4bwmhKpOAib5lT3is3wEuC5I28HA4Oz06co34o2e8y8fBAwKso03gDcy3QljjAkgNjqKey6tTcf6ZzPw4+X8/aNlpC7bxtNXN6BK2RKRDq/AsBkTjDHmNNQ+6wzG3tGCJ7rVZ8nmPXQYNou35222CVGzyZKQMcacpqgo4a8tajD5gdY0rVGOR1NXcf0rX5O+wyZEzYolIWOMySVVypbg7Vsu5H/Xnc/6HYfoMnw2I2akc9wmRA3KkpAxxuQiEeHaJlWYltKGdsln8Z/Ja+n24lxW/rg/0qHlS5aEjDEmDCqUimfkTU14+S8XsPPQUbqNmMuQL7+zCVH9WBIyxpgw6tSgItMeaMM1jSvzUtoGugyfzaLN9tcgGSwJGWNMmJUuEct/rjufd/s249jJU1z38tc88ulKDtmEqJaEjDEmr1ySVIHJ97fmlpY1eHf+93QcNou0tTsiHVZEWRIyxpg8VDI+hkevqs+4Oy+meFw0fd5cRMrYpez95VikQ4sIS0LGGBMBTaqX5fMBrbj3stqkLt1G+2EzmbRiO95UmEWHJSFjjImQ+Jho/tahLqn9W1GxdHHufv8bXlx6lB0His6EqJaEjDEmwpIrncknd1/MwM7nsnznSdoNncnYxVuKxFGRJSFjjMkHYqKjuLPNOTzZsjjnVjyTf45bzs2vL2TLnsORDi2sLAkZY0w+cnbJKEbffhFPdW/A0i376DBsFm/M2cTJQjohqiUhY4zJZ6KihL9cVJ0pD7Smea1yPDFxNde9PI/1Px+MdGi5zpKQMcbkU5XKFOfNPhfy3A2N2LTrF654fg4vfLW+UE2IaknIGGPyMRGhe+PKTE1pQ4f6ifxv6jquemEOK7YWjglRw5qERKSTiKwVkXQRGRhgfbyIjHHrF4hIDZ91g1z5WhHpmFWfIlLT9bHe9Rnnt60eIqIi0tQ9ryEiv4rIUvd4ORyvgTHG5IbyZ8Tz4o0XMOrmJuw9fIxuI+bw7y/WFPgJUcOWhEQkGhgBdAaSgV4ikuxXrS+wV1VrA8OAIa5tMtATqA90AkaKSHQWfQ4BhqlqErDX9Z0RSylgALDAb/sbVLWRe9yZS7tujDFh06H+2Ux5oA03XFiVV2ZupPPw2czfuDvSYeVYOI+EmgHpqrpRVY8Bo4FufnW6AW+75XHA5SIirny0qh5V1U1AuusvYJ+uzWWuD1yf3X228yTwLFB0/gLMGFNolS4ey7+vacgHtzXn5Cml56j5PPTJCg4eOR7p0EIWE8a+KwNbfJ5vBZoHq6OqJ0RkP5Dgyuf7ta3slgP1mQDsU9UT/vVFpDFQVVUnisjf/bZfU0S+BQ4AD6vq7EA7IiL9gH4AiYmJpKWlZbLbwR06dCjHbcPJ4gqNxRUaiys0ocb10AUwPj2GDxb8wBfLttC7fhznVwjPV3s4XrNwJiEJUOY/0D1YnWDlgY7cgtYXkSi803x9AqzfDlRT1d0i0gSYICL1VfXAnzpSHQWMAmjatKm2bds2QHdZS0tLI6dtw8niCo3FFRqLKzQ5iatjO/j2h738c9xyhi05RPdGCTxyVX3KlYzLunGYY8tKOE/HbQWq+jyvAmwLVkdEYoDSwJ5M2gYr3wWUcX34lpcCGgBpIrIZuAhIFZGm7lTfbgBVXQJsAOqcxv4aY0zENK5WlokDWnHf5UlMXL6d9kNn8tmybfl+6p9wJqFFQJIbtRaHN9Ag1a9OKtDbLfcApqv3iqUCPd3ouZpAErAwWJ+uzQzXB67PT1V1v6qWV9UaqloD7xRfV1VdLCIV3EAHRKSW28bGcLwQxhiTF+JjonmgfR0mDmhFlbLFuffDb7n9nSX8tD//Xg4PWxJy12f6A5OBNcBYVV0lIk+ISFdX7XUgQUTSgRRgoGu7ChgLrAa+BO5R1ZPB+nR9PQikuL4SXN+ZaQ0sF5FleAMa7lRVu+euMabAO/fsMxl/d0se6lKPOek7aT90Jh8u/CFfHhWF85oQqjoJmORX9ojP8hHguiBtBwODs9OnK9+IN3ous3ja+ix/DHyc6Q4YY0wBFR0l3N66Fu2TExk4fjmDxq8gdek2nrn2PKonlIx0eL+xGROMMaYQq1G+JB/cdhFPX30eK3/cT8fnZvHa7I35ZkJUS0LGGFPIRUUJNzavxpSU1rQ8pzxPfb6Ga16ax9qfIj8hqiUhY4wpIiqWLs5rvZvyfK/GbNlzmCtfmM1z09Zx7ETkJkS1JGSMMUWIiND1/EpMS2lDl/Mq8ty09Vz1whyWbtkXkXgsCRljTBFUrmQcw3s25vXeTdn/63GuGTmXwZ+v5tdjeTshqiUhY4wpwi6vl8iUlNb0bFaNV2dvouNzs5i3YVeebd+SkDHGFHFnFovl6avP48PbL0IEbnx1AYPGr+BAHkyIaknIGGMMAC3OSeDL+1rTr3Utxiz6gfZDZzJt9c9h3aYlIWOMMb8pHhfN/3Wpxyd3t6RsiThue2cxAz78lt2HjoZle5aEjDHG/Mn5VcuQ2r8VKe3r8MXK7bQbOpM1u3N/0IIlIWOMMQHFxUQx4PIkPh9wCQ0qlyaxZKC75pweS0LGGGMyVSexFO/2bU65YrmfMiwJGWOMiRhLQsYYYyLGkpAxxpiIsSRkjDEmYiwJGWOMiRhLQsYYYyLGkpAxxpiIsSRkjDEmYkQ1f9xnvKAQkZ3A9zlsXh7IuznSs8/iCo3FFRqLKzT5NS44vdiqq2oF/0JLQnlIRBaratNIx+HP4gqNxRUaiys0+TUuCE9sdjrOGGNMxFgSMsYYEzGWhPLWqEgHEITFFRqLKzQWV2jya1wQhtjsmpAxxpiIsSMhY4wxEWNJyBhjTMRYEsojItJJRNaKSLqIDMyD7b0hIjtEZKVPWTkRmSoi692/ZV25iMjzLrblInKBT5verv56Eel9mjFVFZEZIrJGRFaJyH35IS7XXzERWSgiy1xsj7vymiKywG1njIjEufJ49zzdra/h09cgV75WRDrmQmzRIvKtiEzMLzG5PjeLyAoRWSoii11Zfngvy4jIOBH5zn3WWkQ6LhGp616njMcBEbk/0nG5/h5wn/mVIvKh+7+Qd58xVbVHmB9ANLABqAXEAcuA5DBvszVwAbDSp+xZYKBbHggMcctdgC8AAS4CFrjycsBG929Zt1z2NGKqCFzglksB64DkSMfl+hTgDLccCyxw2xwL9HTlLwN3ueW7gZfdck9gjFtOdu9vPFDTve/RpxlbCvABMNE9j3hMrt/NQHm/svzwXr4N3OaW44Ay+SEun/iigZ+A6pGOC6gMbAKK+3y2+uTlZyzXvvTskekb3QKY7PN8EDAoD7Zbgz8mobVARbdcEVjrll8BevnXA3oBr/iU/6FeLsT3KdA+H8ZVAvgGaI731+Ex/u8jMBlo4ZZjXD3xf2996+UwlirAV8BlwES3jYjG5NPPZv6chCL6XgJn4n2pSn6Kyy+WDsDc/BAXXhLagpfUYtxnrGNefsbsdFzeyHijM2x1ZXktUVW3A7h/z3LlweILW9zuML4x3hFHvojLnfZaCuwApuL9mtunqicCbOe3GNz6/UBCGGJ7DvgncMo9T8gHMWVQYIqILBGRfq4s0u9lLWAn8KY7hfmaiJTMB3H56gl86JYjGpeq/gj8F/gB2I73mVlCHn7GLAnlDQlQlp/GxgeLLyxxi8gZwMfA/ap6IL/EpaonVbUR3tFHM6BeJtsJe2wiciWwQ1WX+BZHMiY/LVX1AqAzcI+ItM6kbl7FFoN3GvolVW0M/IJ3mivScXkb866tdAU+yqpqXsTlrkF1wzuFVgkoifd+BttGrsdlSShvbAWq+jyvAmyLQBw/i0hFAPfvDlceLL5cj1tEYvES0PuqOj6/xOVLVfcBaXjn4suISEyA7fwWg1tfGtiTy7G1BLqKyGZgNN4pueciHNNvVHWb+3cH8Ale4o70e7kV2KqqC9zzcXhJKdJxZegMfKOqP7vnkY6rHbBJVXeq6nFgPHAxefgZsySUNxYBSW7ESRze4XhqBOJIBTJG0/TGuyaTUf5XNyLnImC/OzUwGeggImXdL6YOrixHRESA14E1qjo0v8TlYqsgImXccnG8/5xrgBlAjyCxZcTcA5iu3snwVKCnG0VUE0gCFuYkJlUdpKpVVLUG3mdmuqreFMmYMohISREplbGM9x6sJMLvpar+BGwRkbqu6HJgdaTj8tGL30/FZWw/knH9AFwkIiXc/8+M1yvvPmO5caHNHtm6ANgFbzTYBuChPNjeh3jneI/j/Urpi3fu9itgvfu3nKsrwAgX2wqgqU8/twLp7nHLacbUiv9v735CrCrDOI5/fxQ4CTGLcj8OKIKhtlBIDWYlKC7UjdDSTQkZCSJiELQb0I3bcBGEupJ0I+pKBUVSzBlUCgQ3LYoo8S9G1uPieY7eLv6bvM57od8HDnfOn/vynnvO8PC+99znySH6NHC5f8anDQAAAqJJREFUlnWt+1XtLQF+qL5dAb6s7eP1z3SdnEKZU9tHav167R/vaeuL6vNPwNoBXc8Jnjwd17xP1YepWq529/SQXMtlwMW6lkfIp8iGoV9zgd+B0Z5tw9Cvr4Af677/lnzCbdbuMaftMTOzZjwdZ2ZmzTgImZlZMw5CZmbWjIOQmZk14yBkZmbNOAiZzTJJd+t1TNJHA257d9/6uUG2bzZoDkJm7YwBMwpCkt54wSH/CkIRsXKGfTKbVQ5CZu1MAh8q68tsrwSqeyRdqBoyHwNImlDWYTpI/nARSUcqcejVLnmopEngrWrvQG3rRl2qtq8oawBt7mn7lJ7U3zlQv5xH0qSka9WXvbP+6dj/wpsvPsTMXpNdwI6IWA9QweRWRCyXNAc4K+lkHbsCeC8ibtT6loj4o1IMXZB0OCJ2Sfo0Mglrv01kJoGlwLv1njO1731gMZnr6yywStI1YCOwKCKiS2lkNmgeCZkNjzVkvrDLZImLd8gcXADf9wQggM8kTQHnycSRC3i+1cChyEzhvwKngeU9bf8cEf+QqZTGgNvAA2C/pE3A/Vc+O7OncBAyGx4CtkXEslrmR0Q3Err3+CBpgkyw+kFELCVz3o28RNvP8mfP33+TxcwekqOvw8AG4PiMzsTsJTkImbVzhyxz3jkBbK1yF0haWBmq+40CNyPivqRFZMmJzl/d+/ucATbX907zyPLvz8xyrKz5NBoRx4DPyak8s4Hzd0Jm7UwDD2ta7RtgHzkVdqkeDviNHIX0Ow58ImmazFh8vmff18C0pEuRZR8635FlmqfITOY7I+KXCmJP8zZwVNIIOYra/t9O0ez5nEXbzMya8XScmZk14yBkZmbNOAiZmVkzDkJmZtaMg5CZmTXjIGRmZs04CJmZWTOPADp3GXhZ7oCjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y=y.reshape(2*Npoints,)\n",
    "theta_est = VanillaLogisticRegressionNonLinear(X,y)\n",
    "print(f\"theta estimated [ theta0,theta1, theta2]: {theta_est}\")\n",
    "radius_true = (innerRadiusTrue+outerRadiusTrue)/2\n",
    "radius_est = np.sqrt(theta_est[0]**2 + theta_est[1]**2)\n",
    "print(f\"Radius estimated : {radius_est}\")\n",
    "print(f\"Original Radius : {radius_true}\")\n",
    "y_pred_algo = list(map(lambda q: 0 if q<radius_est else 1,(np.sqrt(X[:,0]**2+X[:,1]**2))))\n",
    "print(y_pred_algo[Npoints-5:Npoints+6])\n",
    "y_pred_algo = np.array(y_pred_algo)\n",
    "y = y.reshape(2*Npoints)\n",
    "sumErr = sum(abs(y_pred_algo-y))\n",
    "print(f\"Sum error: {sumErr}\")\n",
    "errRate = sumErr/(2*Npoints)\n",
    "print(f\"Error probability : {errRate}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}