Chapter_4_Homework.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# HW3_Yue Shuai_3036226"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4_3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAHVCAYAAADLvzPyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XlUlPXiP/D3ZwZwGVxywauoqWl53S6Iim1Xh30ELuJSYZqae5qnMk1TQ6HytpfX61dts8w0N1SQkXWwrDQXyO1qWW5ouVuBIjDz/P5Q/CEiDDDD55mZ9+uczrnAMLzxHn0/z/PZhKIoICIiIvXQyA5AREREt2M5ExERqQzLmYiISGVYzkRERCrDciYiIlIZljMREZHKsJyJiIhUhuVMRESkMixnIiIilXGT9YObNWumtGvXTtaPJyIiqlV79uy5oChKc2teK62c27Vrh927d8v68URERLVKCHHC2tfysTYREZHKsJyJiIhUhuVMRESkMixnIiIilWE5ExERqQzLmYiISGVYzkRERCrDciYiIlIZljMREZHKsJyJiIhUhuVMRESkMixnIiIilWE5ExERqQzLmYiISGVYzkRERCrDciYiIlIZljMREZHKsJyJiIhUhuVMRESkMixnIiIilWE5ExERqQzLmYiISGVYzkRERCrDciYiIlIZljMREZHKsJyJiIhUhuVMRESkMiznmwoKCrB9+3bZMYiIiFjOAFBUVIRBgwYhICAAmzdvlh2HiIhcHMsZgJubG7p06YKioiIMGTIESUlJsiMREZELYzkDEELgrbfewnPPPYeioiIMHjwYycnJsmMREZGLqrSchRCfCCHOCSEO3OXrQgixUAhxVAixTwjR0/Yx7U8IgXfffRdTp05FYWEhBg0ahK1bt8qORUQEs9mMpKQkxMfHIykpCWazWXYksjNr7pyXAwir4OsGAJ1u/jcewP/VPJYcQgi8//77mDJlCq5fv46BAwciNTVVdiwicmFmsxnRoeGIjZmIq7EJiI2ZiOjQcBa0k6u0nBVF+RrApQpeEgXgc+WGHQAaCyFa2ipgbRNCYOHChZg0aRKuX7+OqKgopKeny45FRC7KaDTi9M4D2JFnwAKlN3bkGZC78wCMRqPsaGRHthhz9gZwqtTHuTc/dwchxHghxG4hxO7z58/b4EfbhxACixYtwoQJE1BQUIDIyEhkZmbKjkVELig7Oxsh+V5wv/nPtTs0CM33Qk5OjuRkZE+2KGdRzueU8l6oKMoyRVF6KYrSq3nz5jb40faj0WiwePFijB07FgUFBYiIiIDJZJIdi4hcjK+vL1J151AECwCgCBak6M7Bx8dHcjKyJ1uUcy6ANqU+bg3gjA3eVzqNRoOlS5fi6aefxrVr1xAREYFt27bJjkVETuhuk74MBgO8/bvB39OIWWIX/D2NaO3fDQaDQXJisiehKOXe5N7+IiHaAUhSFKVbOV8LBzAFwAAA/gAWKorSp7L37NWrl7J79+6q5pXCYrFgzJgxWL58OXQ6HYxGIx599FHZsYjISZRM+jq98wBC8r2QqjsHb/9uSEjZAq1WC7PZDKPRiJycHPj4+MBgMECr1cqOTVUkhNijKEovq15bWTkLIVYB6A+gGYCzAGIBuAOAoihLhBACwCLcmNF9FcBoRVEqbV1HKmfgxl+e0aNHY8WKFdDpdEhJScHDDz8sOxYROYGkpCTExkzEjjwD3KFBESzw9zQibtUSREREyI5HNlKVcnar7AWKosRU8nUFwGQrszksrVaLTz/9FBaLBStXrkRYWBhSU1Px4IMPyo5GRA6uoklfLGfXxB3CqkCr1eKzzz5DTEwM8vLyEBoaih07dsiORUQOjpO+qCyWcxVptVp8/vnnePzxx/HXX38hNDQUP/zwg+xYROTAOOmLyrJqQpg9ONqYc1nFxcUYNmwY1q5di0aNGiE9PR29elk1lEBEdAdO+nJ+Np0QZi+OXs7AjaMmY2JisH79ejRu3Bjp6enw8/OTHYuIiFSoKuXMx9o14O7ujlWrViE6OhpXrlxBcHAwsrOzZcciIiIHx3KuIXd3d6xevRpRUVG4fPkygoKC8OOPP8qORUREDozlbAMeHh5Ys2YNIiMjcenSJQQGBmLfvn2yYxERkYNiOduIh4cH1q5di/DwcFy8eBGBgYE4cKDcI7CJiIgqxHK2oTp16mDdunUICwvDhQsXEBAQgIMHD8qORUREDoblbGN169ZFQkICQkJCcP78eQQEBOB///uf7FhERORAWM52ULduXWzcuBHBwcE4d+4c9Ho9Dh8+LDsWERE5CJazndSrVw8bN25EYGAgzp49C71ejyNHjsiORUREDoDlbEf169fH5s2bodfr8fvvvyMgIAA///yz7FhERKRyLGc7q1+/PhITE9GvXz+cOXMGer0eR48elR2LiIhUjOVcC3Q6HbZs2YJHH30Up0+fhl6vxy+//CI7FhERqRTLuZbodDokJyfjkUceQW5uLvR6PY4dOyY7FhERqRDLuRZ5enoiOTkZDz30EE6dOgW9Xo/jx4/LjkVERCrDcq5lDRo0gNFoRN++fXHixAno9XqcOHFCdiwiIlIRlrMEDRs2xNatW9GnTx8cP34cer0ep06dkh2LiIhUguUsSaNGjZCSkoLevXvj2LFj0Ov1yM3NlR2LiIhUgOUsUePGjZGSkgI/Pz/88ssv0Ov1OH36tOxYREQkGctZsnvuuQepqanw9fXF0aNHERAQgN9++012LCIikojlrAJNmjRBeno6fHx88NNPP0Gv17OgiYhcGMtZJUoKukePHjhy5AgCAgJw9uxZ2bGIiEgClrOKNG3aFOnp6ejWrRsOHz6MgIAAnDt3TnYsIiKqZSxnlWnevDkyMjLQtWtXHDp0CAEBATh//rzsWEREVItYzirk5eWFjIwM/P3vf8fBgwcREBCACxcuyI5FRES1hOWsUi1atEBmZiY6d+6MAwcOIDAwEBcvXpQdi4iIagHLWcX+9re/ITMzE/fffz/27duHoKAgXLp0SXYsIiKyM5azyrVs2RImkwmdOnVCTk4OgoKCcPnyZdmxiIjIjljODqBVq1YwmUzo2LEjsrOzERwcjCtXrsiORUREdsJydhDe3t4wmUzo0KED9uzZg5CQEBY0EZGTYjk7kNatW8NkMqF9+/bYtWsXwsLC8Oeff8qORURENsZydjBt27aFyWRCu3btsHPnToSGhrKgiYicDMvZAd17770wmUxo27YtduzYAYPBgL/++kt2LCIishGWs4Nq164dsrKy0KZNG3z33XcYMGAA8vLyZMciIiIbYDk7sPbt28NkMqF169bYvn07wsPDkZ+fLzsWERHVEMvZwd13330wmUzw9vbG119/jYiICFy9elV2LCIiqgGWsxPo2LEjTCYTWrZsiaysLERGRrKgiYgcGMvZSXTq1Akmk+nWlp9RUVG4du2a7FhERFQNLGcn8sADD8BkMqFFixZIT0/HwIEDUVBQIDsWERFVEcvZyXTu3BmZmZnw8vJCamoqC5qIyAGxnJ1Qly5dkJmZiebNmyMlJQWDBw/G9evXZcciIiIrsZydVNeuXZGRkYFmzZohOTmZBU1E5EBYzk6se/fuSE9PR9OmTbFlyxYMHToUhYWFsmMREVElWM5O7h//+AfS09PRpEkTJCYm4rHHHmNBExGpHMvZBfj4+CAtLQ333HMPNm3ahJiYGBQVFcmORUREd8FydhE9e/ZEWloaGjdujA0bNmDYsGEsaCIilWI5uxA/Pz+kpqaiUaNGWLduHYYPH47i4mLZsYiIqAyWs4vp3bs3UlNT0bBhQ6xZswYjRoxgQRPZiNlsRlJSEuLj45GUlASz2Sw7EjkoN9kBqPb16dMHKSkpCAkJwerVq6HRaPD5559Dq9XKjkbksMxmM6JDw3F65wGE5HshVrcUy/y7ISFli0v83TKbzTAajcjOzoavry8MBoNL/N72wjtnF9W3b19s3boVnp6e+PLLLzF69Ghe5RPVgNFoxOmdB7Ajz4AFSm/syDMgd+cBGI1G2dHsruTCJDZmIq7GJiA2ZiKiQ8P5b0oNsJxd2EMPPQSj0QidTocVK1ZgzJgx/MtEVE3Z2dkIyfeC+81/Vt2hQWi+F3JyciQnsz9XvjCxF5azi3vkkUeQnJyM+vXr47PPPsO4ceNgsVhkxyJyOL6+vkjVnUMRbvz9KYIFKbpz8PHxkZzM/lz5wsReWM6Ef/7zn0hOTka9evXw6aefYsKECSxooioyGAzw9u8Gf08jZold8Pc0orV/NxgMBtnR7M6VL0zsRSiKIuUH9+rVS9m9e7eUn03ly8zMREREBK5du4Zx48ZhyZIl0Gh4/UZkrZJJUTk5OfDx8XGZSVElY865Ow8gNN8LKbpzaO1Ck+GsJYTYoyhKL6tey3Km0jIyMhAREYGCggJMnDgRixcvhhBCdiwiUjm1XJioedY4y5lqJC0tDZGRkbh+/TqeeeYZLFq0iAVNRKpXdjlbqu4cvFV0B1+VcuYzS7pDcHAwNm3ahDp16mDx4sWYOnUqZF3EERFZy5lmjbOcqVyhoaFISEiAh4cHFi1ahOeff54FTUSq5kyzxlnOdFcGgwEbNmyAu7s7PvjgA0ybNo0FTWRj3PLTdpxp1jjHnKlSiYmJGDx4MIqKivDiiy/izTff5Bg0kQ2ofYzU0ah91jjHnMmmIiMjsXbtWri5ueHtt9/GrFmzeAdNZAPONEaqBlqtFgkpWxC3agl0cYMQt2qJaoq5qqwqZyFEmBDiiBDiqBBiZjlfbyuEMAkhsoUQ+4QQA2wflWSKiorCmjVr4ObmhjfeeAOzZ89mQRPVUHZ2NoLymiMFpxCPvUjBKQTnNXfIMVK10Gq1iIiIwJw5cxAREeGQxQxYUc5CCC2A/wIwAOgCIEYI0aXMy+YAWKMoii+AJwAstnVQki86OhqrV6+GVqvFggULMHfuXBY0UQ306NEDn2l/Qiz24CqKEYs9+Ez7E7p37y47GklmzZ1zHwBHFUX5VVGUQgCrAUSVeY0CoOHN/90IwBnbRSQ1GTx4MFatWgWtVovXXnsN8+fPlx2JyCGZzWbs2rULZouCueiJV9ELOzAQTVBXdjRSAWvK2RvAqVIf5978XGnzAAwXQuQCSAbwrE3SkSoNHToUK1euhEajwfz58xEXFyc7ErkANcxqtlWGkolLm95ahqeV+xGPvYhGKjQA/mW5F/v377dtcHI4bla8prxpuWWfZcYAWK4oyjtCiAcBrBBCdFMU5bbTE4QQ4wGMB4C2bdtWJy+pxOOPPw6LxYLhw4cjNjYWWq0Ws2fPlh2LnFTZWc2xuqVYVsuzcG2ZoWQi2N7CKLhDgyJY4I8EJOIkUnXnENu9O5KSkqzaglLN21VS9VlTzrkA2pT6uDXufGw9BkAYACiK8r0Qoi6AZgDOlX6RoijLACwDbiylqmZmUomYmBiYzWY89dRTmDNnDjQaDWbNmiU7Fjmh0rOa3aFBXJ4F/juNMBqNiIiIcLgM5W2WEYzWeMbje/j1eRAfLVyMMz8crPQiQA0XLWQf1jzW3gWgkxCivRDCAzcmfG0u85qTAAIBQAjxdwB1AZy3ZVBSp+HDh2P58uUQQuDll1/GG2+8ITsSOSE17PxkywzlbZaxxeMMxr70HMZOfQZnfjh4+/KqHfvLXV7FpVjOq9JyVhSlGMAUACkA/ocbs7IPCiHihBD/uvmyaQDGCSF+BLAKwCiF03hdxlNPPYVPPvkEQgjMnDkTb7/9tuxI5GTUsPOTLTOUd/Zzh0d7IjY2Fjk5OQjKa3b7XXV+c+zdu/eO91HDRQvZhzWPtaEoSjJuTPQq/blXSv3vQwAetm00ciSjRo2CxWLBmDFjMH36dGg0GrzwwguyY5GTMBgMWObfDf47jbft/GQwGBwyQ8lmGSVHLMaVOmLRbDZjM07gVfS+NR69CSfwWHHxHe/j6+uLWN1SxOVZbr02RXcOcQ64XSXdzqpyJrLG008/DbPZjPHjx2PatGnQaDR47rnnZMciJ1BRmTlqhpLNMsqOV2u1WggI+CMBoWiDFJyCBgJubnf+c62GixayD+6tTTa3dOlSTJw4EQCwcOFCPPssV9a5OkebUSwzb1JSEl55YgJeye+GA7iEbmiC+boDiF+9tNyJZyVZc3Jy4CPhooWsV5W9tVnOZBf/93//h2eeeQYAsGjRIkyePFlyIpLF0Q53kJ1X7Yc3UPVVpZz5WJvsYtKkSbBYLJgyZQqmTJkCrVZ7626aXIsalkFVhey8aniET/LxVCqym8mTJ2PhwoUAbpT1smXLJCciGRxtRrEa8jrL4Q1UfSxnsqtnn30W7733HgBgwoQJ+OijjyQnotpW28ugarrFphqWbRFxzJlqxbvvvotp06ZBCIGPP/4Yo0ePlh2JakltjqHaYrzYkcZ8HW2inavjhDBSpbfeegszZsyAEAKffvopRo4cKTsS1ZLamlGclJSE2JiJt8aLi2CBv6cRcauWVGm82BFmQMueuEZVxwlhpErTp0+HxWLBzJkzMXr0aGg0GowYMUJ2LKoFd1vTa2sVjRdX5WfXVt6akD1xjeyLY85Uq1566SW89tprUBQFo0aNwsqVK2VHIifiSuPFapi4RvbDO2eqdS+//DIsFgvmzp2Lp556ClqtFk888YTsWOQEamvHrJLH3nv27IHZbIZWq4Wfn5/Vj79tMVbMrTude8ydY84kzfz58zFv3jxoNBqsWrUKjz32mOxI5ATsPV5ceqw3KK8ZNuMEBATq6uqhdd/ulY752mqs2JEmrlVXReXriGPuVRlzhqIoUv7z8/NTiObOnasAULRarbJ27VrZcYgqlZiYqPT09FYKMVZRMF4pxFjFF02VBIQovp7eSmJiYtW/34rvK09xcbGSmJioxMfHK4mJiUpxcXF1fy3VKS4uViIDQ5Went7KTOGr9PT0ViIDQ2/9jrb8c6wtAHYrVnYkx5xJqvnz52P27Nkwm8144oknsH79etmRyMHUdF1zVZU71os2OIBLVo352nKs2Jk3K6nsrGpnH3NnOZNUQgjEx8dj5syZtwo6ISFBdixyECWPNmNjJuJqbAJiYyYiOjTcrgVd7qQznEI3NLFq8pmMSWu1fQFTXaVzrl69GkF5ze9avk4/+c/aW2xb/8fH2lSaxWJRpk+frgBQ3NzclI0bN8qORA5AxqPNksetvrpWygz8Q3kAjZW/4x7FR9fqtseulX7/zce1vmUe19or790eD6tF2Zw96rZQvDT1lPUIUuLQS0lAiOKja3Xr/9va/nO0BVThsTYnhJFqKIqC6dOn45133oG7uzvWr1+PyMhI2bFIxeLj43E1NgELlN63PjdL7IIubhDmzJljt59bMlFp7969KC4uhpubG3r27HlrVnhlM4hrc5MTW23MYm9lcxagGG3wJZqiLqJwLzbjBK7e44Gffz8FDw8PAI6xWUxp3ISEHJIQAm+99RbMZjPef/99DBkyBAkJCRgwYIDsaFRLqro0RtZyorttUlJ2BnGsbimWlTODuDY3ObHVxiz2VjZnOk6jBeohG4PhDg1eRW/4FyYjNTX1Vm5H2CymujjmTKoihMC7776LqVOnorCwENHR0bcmgJBzq874scFggLd/N/h7GjFL7IK/p9Eu65qtYTabMX/+fBz7Zm+5k5hkjfs6yths2Zy7cQED0Pb2i4qrLZxmwldl+FibVElRFEydOhWLFi1CnTp1sGnTJoSGhsqORXZU3cevani0WXJhsfub7zCisAPegP+tr80Su1Bv3kDs/vq7Ctfk2mtDDUdZD10255o6J1CnCPjRPEjVj+Orgo+1yeEJIbBw4UJYLBYsXrwYAwcORGJiIoKCgmRHIzup7uNXNTzaLFn2s7jwQcRjL4pw+2P2iOLiCvfBtvZxeHVotVokpGy5dQETp9Kx2bI53+3eHR8tXAz/H+y725tasZxJtYQQ+M9//gOz2YylS5ciMjISSUlJCAwMlB2Nqqmiu0NH3o6y5MIiEm3xCQ7DHwkIRmts8TiDDv49odVqK7zwsPchFmq4gLFG2ZwlfzZqvqiwF445k6ppNBosXrwY48aNQ0FBASIjI2EymWTHomqobExZTePHVVUyXmoBkIAQvAI/rPD4FYNemoiElC3w8/OrcNzX2TfUqC5n3mSlMhxzJodgsVgwfvx4fPzxx6hfvz6Sk5PRr18/2bGoCqwZU1bD+HF1VDauW9nXrR1vr2xcWk0HQagpi1pUZcyZ5UwOw2KxYMyYMVi+fDnq16+PrVu34tFHH5Udi6wka01ybanswqKir1szaauygx7UdBCEmrKoCQ++IKdVXFysPPXUUwoARafTKdu3b5cdiazkiAcVlFZyyERcXJxdDpmo7BCLyv781PTnq6YsagIefEHOSqvV4pNPPsHw4cORn5+PsLAwfPfdd7JjkRUceUy5Nvbwrmx8tbJxaTWNW6spi6NiOZPD0Wq1WL58OYYNG4a8vDyEhYVhx44dsmNRJUqWysStWgJd3CDErVpy2yNZNRzMcLcclZ2QVBsq20xETZuNqCmLo+JSKnJIWq0Wn332GSwWC1avXo3Q0FCkpaWhT58+sqNRBcpb0nO3Nb7rkjcjNTXVphOKKpqkVNFaY9lbYJrNZpjNZhR66eBXnIgB11shtcy6X4PBgGX+3eC/U/66YDVlcVScEEYOrbi4GMOGDcPatWvRqFEjpKWloXfv3pV/I6lGeTOV++iSUa9jS1z/5XebTSiqbJJSRTOmAUg7PKJ07qC85lhf9yTcWjXBG+++fcfjbzXNdldTFrXgDmHkMtzc3LBy5UpYLBasX78ewcHByMjIgJ+fn+xoZKXy7ko75tfF4f+dxN7CqDs25TAYDNVaolPZRh8V3R3PmjVL2p1g2dyvFvSC/zkjtFrtHb+3mjYbUVMWR8QxZ3J47u7uWLVqFaKjo/HHH38gODgY2dnZsmORlcobn/zW4zzCi7zvKMq9e/dWe2JWZZOUKhonrWi8vDR7jJ1zcpVrYjmTU3B3d8fq1asRFRWFy5cvIygoiP94SVSVkipvFnebv3dEmu78HUVZXGqP6qpOzKpsklJls8krm01trxndnFzlmjjmTE6lsLAQQ4YMQWJiIpo2bYrMzEz06NFDdiyXUp0NKMqOT4aEhGDIgH/dsSmH36MPomD+pmptZGLtRh/VHSet7qlalXGkU6W4I1jFuEMYubTr169j8ODB2LJlC5o1awaTyYRu3brJjuUybFVS5RWl0Wis0Xvbc5KSPXdAU/vkKu4IZh2WM7m8goICDBo0CEajEc2bN4fJZELXrl1lx3IJ9i4ptd5F2uvO2RG48u9eFVUpZ445k1OqW7cuNmzYgNDQUJw/fx4BAQE4dOiQ7FguwZ5jpNZOzJLBkXdAqylOWrM93jmTU7t27RqioqKQlpaGFi1aICsrC507d5Ydy6mp+e7W3mry+LmqY7ZqGuPlnbN1+FibqJRr164hMjISGRkZaNmyJbKysnD//ffLjuXU1D5GWhO2KsXS79OjRw98tHAxzvxw0KoxW7WN8bryBVlVsJyJyrh69SoiIiJgMpnQqlUrZGVloVOnTrJjkUpYW7i2KsWy75Nc5wzOFv2F4+YnUBduld55qvFO1ZkvyGyFY85EZdSvXx+JiYno378/zpw5A71ej6NHj8qORSpQlfXJNT0Ao2T998iRI3Hi25xb77O7IBLNzHWQjtMAKh+zrckYr70OGalsHThVDcuZXIZOp0NSUhL++c9/4vTp09Dr9fj1119lxyLJqlK45ZVicF5zfPXVV7dK7m7lV/oi4OTKTIQV/O229wlHW+zBBQCVT6Kr7qS72jj6kmyDe2uTS9HpdNiyZQsMBgO2b98OvV6PrKwstG/fXnY0kqQqJ075+voiVrcUcXmWW4+Tt+AECtadRPRv4ViXvBlDBvzrtpOtXr/vbwiNjoTZbEbujv3YmT8AKTiFWOxBEf7/+yRqT+G6O1Bw3VLp3t3VPfWpsv3FST1YzuRyPD09kZycjLCwMHz33Xe3Crpdu3ayo1E11GSCltlsRlFREZLccxFX6HerKFN05xBXzl1oSSn6fZsIQ8HfkIZcdEADrCkIwkM7U/Dqq6/eUX7df1yHAz9+jsMef8FcWAwNAAPaYBn+B1+sRzjuRZrnedzXpxfGTn0G+/fvR2z37gCA119/vdzfqWRJWckYb5yVY7yyj76siJpmn6sBJ4SRy/rzzz8RFhaG77//Hu3atcO2bdvQtm1b2bGoCmoyQevWDOMd+1GQfw0WKIjCvUjTnUfrvt0rnCk9YsQI5K4yYQb+AQPaQAsNZold+EF/D/qYLt+2ActM7IQn3PESfNAd6/Bv9MFAtEMBitG17kY8NGQAHn/88VtlZM+Z2GqcSAaob/a5vXBCGJEVGjZsCKPRiD59+uD48ePo378/Tp06JTsWVUFNJmiVfO/O/AHYj8H4N/pghceviHhhTKWl8MuhwziPAoTeLOYiWJBS/ywefvjhO8aCU5ELHzSFOzSIwr14xuN7zBK78JBnCro+3BvLly+/bQJVTSedVUStG6XY83d2VCxncmmNGjVCSkoKevXqhWPHjqF///7Izc2VHYusVJNZy6W/VwsNBqIdRhZ1hIeHR4XFbDQaUXT0LO5DQ/gjATOxE92xDnU7tsScOXNuld9M8QO6Yx1aQQcD2qAIFqTpzmPsS89VuLuZNZPOqqsmO6zZa5Y3wB3GysNyJpfXuHFjpKamws/PD7/++iv0ej1Onz4tOxZZoSZbhVb3e7OzsxF6tQU2IQRx6AVPuKM7miBs0L/g4eFxq/zqz4tGo3+0x2ldMeaIPTfuUvt2R2xs7G3LjcqWXo8ePe7ItQUn8P26ZJvMrK7Okid7z/LmsZjlUBRFyn9+fn4KkZpcvHhR8fX1VQAonTp1Uk6fPi07ElWiuLhYiQwMVXw9vZWZwlfx9fRWIgNDleLiYrt9b2JiotLT01spxFhFwXilEGMVX09vJTExsdyfkZiYqMTHxyuJiYl3vHdJhp43M/T09FYiAkKUiIAQpXvdFsoM/EPxRVMlEm2Va3j6rj/Hmt81MTFRiYuLKzdHZaryO1dHTf5/dCQAditWdiQZ+7IzAAAgAElEQVQnhBGVcunSJQQGBiInJwcPPPAATCYTWrZsKTsWVcAW+1mXPkc6NTW1whnDttyq8m4TtGK/+C/Wrl2LkyszMQM9bpt0VtXTvWwx2cqeJ42VzunsO4xx+06iGrh48SICAgKwb98+dO7cGVlZWWjRooXsWGRnVSkxa4uksuVBFZWej4/PXYtbq9VaveTIFjO01TrL29FUpZy5zpmojKZNmyIjIwMBAQHYv38/AgICYDKZ4OXlJTsa2VFVNugoGbetqJjKln2sbimWlSn78jY1KVljXd5GI959ut52QEbZ9yzvYqA6a5vLvk9ISEi1Nj2h6mM5E5WjWbNmyMjIgF6vx8GDBxEYGIjMzEw0b95cdjSqQE02srD1Bh3WlH1FO32Vt9GI2WxG3PDJ5b6nwWAo92JgzLOTEHeXC4C7/RmW9z7rkjcjNTW1SpueUPWxnInuonnz5sjMzIRer8eBAwduFXSzZs1kR6NyWHOnWpGK7mKrw5qyr2ynr7J36PHx8RUuOSrvYgDP4sbyLivveu92UZGamlrp0wKyHS6lIqqAl5cXMjMz0blzZ+zfvx9BQUG4ePGi7FhUjppuZGHrDTqsXR5UlaVNFb3n3S4G9u/fX6W1zVxzrA4sZ6JKtGjRApmZmXjggQfw448/IigoCJcuXZIdi8qoaanUZIOO8thjN66K3rOi4rbVBQDVHs7WJrLSmTNn0L9/f/z888/o2bMn0tPTcc8998iORTepcUaxPZYH3e09bbXEy5ZLxeh2XEpFZCenT59G//79cfToUfj5+SE9PR2NGzeWHYvAUgH+f3Hv3bsXxcXF0Gq18PPzq/JFgSusOZaB5UxkR7m5uejXrx9+/fVX9O7dG2lpaWjUqJHsWASWCuA6Jzw5IpYzkZ2dPHkS/fv3x7Fjx+Dv74/U1FQ0bNhQdiyyUm2dHSzjjOKKHu8bDAaemSwRy5moFpw4cQL9+/fH8ePH8eCDD2Lr1q0saAdgNpsxMGQATn73I8IKWmJr3d/Q9qF/YGNqcrUe/d6t6GTdwd5t17F68wZi99ff8Y5aIp7nTFQL7r33XphMJrRt2xbff/89BgwYgL/++kt2LKpEUlISftm2G7sLIvEG+mB3QSR+2bYbSUlJVr9H2VOaXnliAh7y64358+ffOk7R3mcU3+0Ix7vNti4uLuaZyQ6E5UxUA+3atYPJZEKbNm3w7bffYsCAAcjLy5Mdiyqwdu1aRJrb3LbkKsLcBuvWrbP6PUoX76uKH1rlu+HPH4/h2vyNt45T3LNnj93WC1d0cWA2m9GqT9c7lltptVquX3YgLGeiGurQoQNMJhNat26N7du3Izw8HPn5+bJjUQW24ORtd5bJOFml7y+9ptqIU/gNV7EPQ/Bvpc+tO1Kz2Wy39cIVXRzEDZ8MAIj94r+3rdf28/Pj+mUHYlU5CyHChBBHhBBHhRAz7/Kax4QQh4QQB4UQX9o2JpG63XfffTCZTGjVqhW+/vprREREsKBVaujQobigvQ5/JGAWfoA/EnBBex1Dhgyx+j1KPzrOxkUEw/uOO1I3Nzebb0JSorKLg9M/HIRWq71t0xF7bIpC9lNpOQshtAD+C8AAoAuAGCFElzKv6QRgFoCHFUXpCuA5O2QlUrWOHTsiKysLLVu2RFZWFiIjI3H16lXZsaiMiIgI9O73CP6oq+BbnMUfdRX07vcIIiIi7jqOW1bpotuHS9iEE3fckfbs2dOmO46VZs3FQdnH1bbeAY3sq9LZ2kKIBwHMUxQl9ObHswBAUZQFpV7zJoCfFEX5yNofzNna5KyOHDmC/v374/fff0dgYCASExNRr1492bGolPLWQwOo0uzq0ht+bN2wGQVHf0Po1Ra1svlJ6Q1X7surgwO4hH0Yopqd0ah8Nl1KJYQYAiBMUZSxNz8eAcBfUZQppV6zEcBPAB4GoMWNMt9aznuNBzAeANq2bet34sQJ634jIgdz+PBh9O/fH2fPnkVwcDA2bdrEgla5mmz/KWPzE5kXB1Q9ti7noQBCy5RzH0VRni31miQARQAeA9AawDcAuimKcuVu78s7Z3J2hw4dgl6vx7lz5xAaGoqNGzeibt26smNJI2NDjqq42/pgXdwgzJkzR2KyynFnNMdg63XOuQDalPq4NYAz5bxmk6IoRYqiHANwBEAnawIQOasuXbogMzMTzZs3R0pKCgYNGoTr16/LjiVF2aU/JcuN7jamK4Mjn8ZUlVOnyDFYU867AHQSQrQXQngAeALA5jKv2QhADwBCiGYA7gfwqy2DEjmirl27IiMjA82aNYPRaMTgwYNdsqDtvSGHLXA2M6lJpeWsKEoxgCkAUgD8D8AaRVEOCiHihBD/uvmyFAAXhRCHAJgATFcUhSfSEwHo3r07MjIy0LRpU2zZsgVDhw5FYWGh7Fi1qqZnLdeGimYzWzuLm8hWuLc2US3JyclBYGAgLl26hIEDB2LNmjVwd3eXHatWqPGsZWvxlCeyFe6tTaRCPj4+SEtLQ+PGjbFx40Y88cQTKCoqkh2rVjjiI+OSu+WRI0fi0Le78G1eqGofyZPzYTkT1aKePXsiPT0djRs3xoYNGzBs2DCXKGhH2wCj9AQ275UHUbdAwWNIhxkWVT6SJ+fjJjsAkavx8/NDamoqgoODsW7dOmg0GqxcuRJubs7917FkRnFtPMa25jjHir5eegKbOzR4Fb3gjwQYcQqhaIMU3TnEOcAsbnJczv2vAZFK9e7dGykpKQgJCcGaNWug0WiwYsUKpy/o2lB2jDhWtxTLSo0RV/Z1oPwJbEFojTexD694HlD9I3lyfHysTSSJv78/UlJS0KBBA6xevRojR47kLGAbqGzZljXLuspb87y17u+4d3ig6h/Jk3NgORNJ1LdvX2zduhWenp748ssvMWrUKBZ0DVW2bMuaZV3lTWBr97APli9fzk0+qFawnIkke+ihh2A0GqHT6fDFF1/g6aefZkHXQGU7fVmzE5ijTWAj58N1zkQq8c0338BgMCA/Px+jRo3Cxx9/DI2G189VVfrEptB8rzsOgqjs6zJzq3nvcao5mx58YS8sZ6I7bdu2DQMGDMDVq1cxZswYLFu2jAVdDZUdBKG2gyK40YlrYDkTOTCTyYTw8HBcu3YN48aNw5IlS1jQdqSGO1ZH3kGNrMcdwogcmF6vR2JiIurWrYsPP/wQkydPhqyLaGenltOyHGHvcapdLGciFQoMDLxV0EuWLMGUKVNY0HagltOyHPm4SrIPljORSgUFBWHTpk2oU6cOFi9ejKlTp7KgbUwtd6yOuPc42RfLmUjFQkJCsHHjRnh4eGDRokV4/vnnWdA2pJY7Vi7dorI4IYzIASQnJyM6OhqFhYV4/vnn8c4770AIITuWw1PrsipyTpytTeSEkpKSMGjQIBQVFeHFF1/Em2++yYK2AbUtqyLnxXImclKbNm3CkCFDUFxcjBkzZuDf//43C5rIQXApFZGTioqKwtq1a+Hm5oY333wTs2fP5hg0kRNiORM5mIEDB2L16tXQarVYsGABXnnlFRY0kZNhORM5oMGDB2PVqlXQarV49dVXMX/+fNmRiMiGWM5EDmro0KFYuXIlNBoN5s+fj7i4ONmRiMhGWM5EDuzxxx/HF198AY1Gg9jYWLz22muyIxGRDbjJDkBENRMTEwOLxYKnnnoKc+bMgUajwaxZs2THUg01HGxBVFUsZyIn8OSTT8JisWDkyJF4+eWXodVqMWPGDNmxpCt7FGOsbimWcZMRcgB8rE3kJEaMGIFPP/0UQgi89NJLePvtt2VHkk4tB1sQVRXLmciJjBw5Eh999BEAYPr06XjvvfckJ5JLLQdbEFUVy5nIyTz99NP48MMPAQAvvPACPvjgA8mJ5FHLwRaOymw2IykpCfHx8UhKSqr1c65dGbfvJHJSy5Ytw4QJEwAACxcuxLPPPis5Ue3jwRbVV3a8PlV3Dt78s6uRqmzfyQlhRE5q/PjxsFgsmDRpEqZOnQqNRoPJkyfLjlWrSo5iLDnYIo4HW1it9Hi9OzSIy7PAf6cRRqMRERERsuM5PZYzkRObOHEizGYzpkyZgilTpkCr1WLixImyY9UqrVaLiIiIOwqFS6wqVtF4PcvZ/jjmTOTkJk+ejIULFwIAJk2ahGXLlklOJF/JI9vYmIm4GpuA2JiJiA4N55hqKRyvl4vlTOQCnn322VsztydMmICPP/5YciK5uMSqcgaDAd7+3eDvacQssQv+nka09u8Gg8EgO5pL4GNtIhfx3HPPwWw248UXX8S4ceOg0WgwevRo2bGk4CPbynG8Xi6WM5ELmTZtGiwWC2bMmIExY8ZAo9Fg5MiRsmPVOl9fX8TqliIuzwJ3aG49so3jI9vb3G28nuyPj7WJXMz06dOxYMECKIqC0aNHY8WKFbIj1To+siW1450zkQuaOXMmLBYLZs+ejVGjRkGj0eDJJ5+UHavW8JEtqR03ISFyYfHx8XjllVeg0WjwxRdfICYmRnYkIqdVlU1I+FibyIXNnTsX8+bNg8ViwfDhw7FmzRrZkYgILGcilxcbG4u5c+fCYrFg2LBhWLdunexIRC6P5UxEmD9/PmbPng2z2YwnnngCGzZskB2JyKWxnIkIQgjEx8dj5syZMJvNePzxx7Fx40bZsYhcFsuZiADcKOjXX38d06dPR3FxMR577DFs3rxZdiwil8RyJqJbhBB44403MG3aNBQVFWHIkCFISkqSHYvI5bCcieg2Qgi89dZbeP7551FUVITBgwcjOTlZdiwil8JyJqI7CCHwzjvvYOrUqSgsLER0dDS2bt0qOxaRy2A5E1G5hBB4//33MWXKFBQWFmLgwIFITU2VHYvIJbCcieiuhBBYuHAhJk2ahOvXryMqKgrp6emyYxE5PZYzkWRmsxlJSUmIj49HUlISzGaz7Ei3EUJg0aJFmDBhAgoKChAZGYnMzEzZsYicGsuZSCKz2YzQ8GjETIxFbMJVxEyMRWh4tOoKWqPRYPHixRg7diwKCgoQERGBrKws2bGInBbLmUgio9GInQdOI8+wA0rvBcgz7MDOA7kwGo2yo91Bo9Fg6dKlePrpp3Ht2jWEh4dj27ZtsmMROSWWM5FE2dnZyPcKATTuNz6hcUe+VyhycnLkBrsLjUaDDz/8EKNGjcLVq1cRHh6Ob775RnYsIqfDciaSyNfXF7pzqYCl6MYnLEXQnUuBj4+P3GAV0Gg0+OijjzBixAjk5+fDYDBg+/btsmMRORWWM5FEBoMB/t284Wn0h9g1C55Gf/h3aw2DwSA7WoW0Wi0+/fRTPPnkk7cK+rvvvpMdi8hpCEVRpPzgXr16Kbt375bys4nUxGw2w2g0IicnBz4+PjAYDNBqtbJjWcVsNmPEiBFYtWoVGjRogNTUVPTt21d2LCJVEkLsURSll1WvZTkTUU0UFxdj+PDh+Oqrr9CwYUOkpaWhT58+smMRqU5VypmPtYmoRtzc3PDFF19g6NCh+PPPPxESEgJeeBPVDMuZaoXaN9qgmnFzc8PKlSsxePBg/PHHHwgODsbevXtlxyJyWCxnsjtH2WiDasbd3R2rVq1CdHQ0rly5gqCgIGRnZ8uOReSQWM5kd4600QbVjLu7O1avXo2oqChcvnwZQUFB+PHHH2XHInI4LGeyO0fbaINqxsPDA2vWrEFkZCQuXbqEwMBA7Nu3T3YsIofCcia7c8SNNqhmPDw8sHbtWoSHh+PixYsIDAzEgQMHZMcichgsZ7I7R91og2qmTp06WLduHQwGAy5cuICAgAAcPHhQdiwih8B1zlQrHHmjDXsp+TPJzs6Gr6+v0/6ZFBQUYODAgUhJSYGXlxeysrLw97//XXYsolpn801IhBBhAD4AoAXwkaIo/77L64YAWAugt6IoFTYvy5lcWckM9p0HTiPfKwS6c6nw7+aNlC0JTlnQ165dQ1RUFNLS0tCiRQtkZWWhc+fOsmMR1SqbbkIihNAC+C8AA4AuAGKEEF3KeV0DAFMB7KxaXCLX42oz2OvVq4dNmzYhMDAQZ8+ehV6vx5EjR2THIlIta8ac+wA4qijKr4qiFAJYDSCqnNfFA3gTQIEN8xE5JVecwV6vXj1s3rwZer0ev//+O/R6PX7++WfZsYhUyZpy9gZwqtTHuTc/d4sQwhdAG0VRkmyYjchpueoM9vr16yMxMRH9+vXDb7/9Br1ej6NHj8qORaQ61pSzKOdztwaqhRAaAO8BmFbpGwkxXgixWwix+/z589anJHIyrjyDXafTYcuWLXj00Udx+vRp6PV6/PLLL7JjEalKpRPChBAPApinKErozY9nAYCiKAtuftwIwC8A8m5+y98AXALwr4omhXFCGLk6V5/BnpeXh7CwMHz77bdo06YNsrKy0KFDB9mxiOzGprO1hRBuAH4CEAjgNIBdAIYpilLugkUhRBaAFzlbm4gq89dffyEsLAzfffcd2rZti23btqFdu3ayYxHZhU1nayuKUgxgCoAUAP8DsEZRlINCiDghxL9qFpWIXFmDBg1gNBrRt29fnDx5Enq9HidOnJAdi0g6bkJCRNL98ccfCAkJwQ8//ID27dtj27ZtaNOmjexYRDZl0ztnIiJ7a9SoEVJSUtC7d28cO3YM/fv3R25uruxYRNKwnIlIFRo3bozU1FT4+fnh119/hV6vx+nTp2XHIpKC5UxEqtG4cWOkpaWhZ8+eOHr0KPR6Pc6cOSM7FlGtYzkTkarcc889SEtLg4+PD37++WcEBATgt99+kx2LqFaxnIlIdZo0aYL09HT06NEDR44cQUBAAM6ePSs7FlGtYTkTkSo1bdoUGRkZ6N69Ow4fPoyAgACcO3dOdiyiWsFyJqohs9mMpKQkxMfHIykpCWazWXYkp9GsWTNkZGSga9euOHToEAICAsCtf8kVsJyJaqDkXOaYibGITbiKmImxCA2PZkHbUPPmzZGZmYkuXbrg4MGDCAwMxIULF2THIrIrljNRDbjaucyyeHl5ISMjA507d8b+/fsRGBiIixcvyo5FZDcsZ6IaUOu5zM74qP1vf/sbMjMzcf/992Pfvn0ICgrCpUuXZMcisguWMzml2ionNZ7L7MyP2lu2bAmTyYROnTohJyeHBU1Oi+VMTqc2y0mN5zI7+6P2Vq1awWQyoWPHjsjOzkZISAguX74sOxaRTbGcyenUZjlptVqkbEnAqiVxiBukw6olcUjZkiD1XGa1Pmq3JW9vb5hMJtx3333Ys2cPQkNDceXKFdmxiGyG5UxOp7bLSavVIiIiAnPmzEFERITUYgbU+ajdHlq3bg2TyYT27dtj165dCAsLw59//ik7FpFNsJzJ6bhKOd2NGh+120ubNm1gMpnQrl077Ny5kwVNToPnOZPTKRlz3nkgF/leodCdS4F/t9bSHzfXJrPZDKPRiJycHPj4+MBgMDj17378+HH0798fJ06cwMMPPwyj0YgGDRrIjkV0m6qc58xyJqfkauVEwLFjx9CvXz+cOnUKjzzyCIxGIzw9PWXHIrqF5UxELunXX39Fv379kJubi3/+859ITk6GTqeTHYsIQNXKmWPOROQ0OnToAJPJBG9vb3z99deIiIjA1atXZcciqjKWMxE5lY4dO8JkMqFVq1bIyspCZGQkC5ocDsuZnHKrR3JtnTp1QmZm5q0tP6OionDt2jXZsYisxnJ2cc681SO5tgceeAAmkwktWrRAeno6Bg4ciIKCAtmxiKzCcnZxzr7VI7m2zp07IzMzE15eXkhNTUV0dDQLmhwCy9nFucJWj+TaunTpgszMTDRr1gxbt27F4MGDcf36ddmxiCrEcnZxrr6bFrmGrl27IjMzE02bNkVycjILmlSP5eziXGmrR3Jt3bt3R0ZGBpo0aYItW7Zg6NChKCwslB2LqFzchIS4mxa5lJycHAQEBODy5cuIiorCmjVr4OHhITsWuQDuEEZEVIG9e/ciMDAQV65cQXR0NL766iu4u7vLjkVOjjuEERFVoGfPnkhPT0fjxo2RkJCAYcOGoaioSHYsoltYzkTkkvz8/JCamopGjRph3bp1GD58OIqLi2XHIgLAciYiF9a7d2+kpqaiYcOGWLNmDUaMGMGCJlVgORORS+vTpw9SUlLQoEEDrF69GiNHjuQOeSQdy5mIXF7fvn2xdetWeHp64ssvv8To0aNZ0CQVy5mICMBDDz0Eo9EInU6HFStWYMyYMSxokoblTER00yOPPAKj0Yj69evjs88+w7hx42CxWGTHIhfEciYiKuXRRx9FcnIy6tevj08//RTjx49nQVOtYzkTEZXRr18/JCUloV69evj4448xceJEFjTVKpYzEVE59Ho9EhMTUbduXXz44YeYPHkyZO2oSK6H5UxEdBeBgYHYvHkz6tSpgyVLlmDKlCksaKoVLGciogoEBwdj06ZNqFOnDhYvXoypU6eyoMnuWM5ERJUIDQ1FQkICPDw8sGjRIjz//PMsaLIrljMRkRUMBgM2bNgAd3d3fPDBB5g2bRoLmuyG5UxEZKXw8HCsX78e7u7ueO+99zBjxgwWNNkFy5mIqAoiIyOxdu1auLm54e2338bMmTNZ0GRzLGcioiqKiorCmjVr4ObmhjfffBOzZ89mQZNNsZyJiKohOjoaX331FbRaLRYsWIC5c+eyoMlmWM5ERNU0aNAgrFq1ClqtFq+99hrmzZsnOxI5CZYzEVENDB06FCtXroRGo0FcXBzi4uJkRyInwHImIqqhxx9/HF988QU0Gg1iY2Px6quvyo5EDo7lTERkAzExMfj888+h0Wgwd+5cLFiwQHYkcmAsZyIiG3nyySexfPlyCCHw8ssv480335QdiRwUy5mIyIZGjBiBTz75BEIIvPTSS3j77bdlRyIHxHImIrKxUaNG4aOPPgIATJ8+He+++67kRORoWM5ERHbw9NNP48MPPwQATJs2DR988IHkRORIWM5ERHYyduxYLF26FADw3HPP4T//+Y/kROQo3GQHILIls9kMo9GI7Oxs+Pr6wmAwQKvVyo5FLmz8+PEwm8145plnMHXqVGg0GkyePFl2LFI5ljM5DbPZjNDwaOw8cBr5XiHQLY2Ff7dlSNmSwIImqSZNmgSLxYIpU6ZgypQp0Gg0mDRpkuxYpGJ8rE1Ow2g0YueB08gz7IDSewHyDDuw80AujEaj7GhEmDx58q1x52eeeQbLli2TnIjUjOVMTiM7Oxv5XiGAxv3GJzTuyPcKRU5OjtxgRDdNnToV7733HgBgwoQJt2Z0E5XFcian4evrC925VMBSdOMTliLozqXAx8dHbjCiUp577jm88847AG6MR3/yySeSE5EasZzJaRgMBvh384an0R9i1yx4Gv3h3601DAaD7GhEt3nhhRfw1ltvQVEUjB07Fp999pnsSKQynBBGTkOr1SJlSwKMRiNycnLg4xPH2dqkWi+++CLMZjNmzpyJ0aNHQ6PRYMSIEbJjkUoIWYeD9+rVS9m9e7eUn01EpBavv/46Zs+eDY1Gg88//xxPPvmk7EhkJ0KIPYqi9LLmtbxzJqfE9c7kKF5++WVYLBbMnTsXTz31FDQaDWJiYmTHIsmsGnMWQoQJIY4IIY4KIWaW8/UXhBCHhBD7hBAZQoh7bR+VyDol651jJsYiNuEqYibGIjQ8GmazWXY0onLNmTMH8+bNg8ViwfDhw/HVV1/JjkSSVVrOQggtgP8CMADoAiBGCNGlzMuyAfRSFKUHgHUAeE4aScP1zuSIYmNj8corr8BiseDJJ5/E2rVrZUciiay5c+4D4KiiKL8qilIIYDWAqNIvUBTFpCjK1Zsf7gDQ2rYxiazH9c7kqObNm4fZs2fDbDYjJiYGGzZskB2JJLGmnL0BnCr1ce7Nz93NGADl3qIIIcYLIXYLIXafP3/e+pREVcD1zuSohBCIj4/HzJkzYTab8fjjj2Pjxo2yY5EE1pSzKOdz5U7xFkIMB9ALwFvlfV1RlGWKovRSFKVX8+bNrU9JVAVc70yOTAiB119/HTNmzEBxcTGGDh2KzZs3y45Ftcya2dq5ANqU+rg1gDNlXySECAIwG0A/RVGu2yYeUdVxvTM5OiEE/v3vf8NsNuOdd97BkCFDsGHDBkRERMiORrWk0nXOQgg3AD8BCARwGsAuAMMURTlY6jW+uDERLExRlJ+t+cFc50xEVDFFUfDCCy/g/fffh4eHBzZs2IDw8HDZsaiaqrLOudLH2oqiFAOYAiAFwP8ArFEU5aAQIk4I8a+bL3sLgCeAtUKIHCEEn8EQEdWQEALvvvsupk6disLCQgwaNIirDlwEdwgjIlI5RVEwdepULFq0CHXq1MGmTZsQGhoqOxZVkU3vnImISC4hBBYuXIhJkybh+vXrGDhwINLT02XHIjtiORMROQAhBBYtWoSJEyeioKAAkZGRyMzMlB2L7ITlTETkIDQaDf773/9i3LhxKCgoQEREBLKysmTHIjtgORMRORCNRoMlS5ZgzJgxuHbtGsLDw7Ft2zbZscjGWM5ERA5Go9Fg2bJlGD16NK5evYrw8HB88803smORDbGciYgckEajwYcffoinnnoK+fn5GDBgAL799lvZschGWM5ERA5Kq9Xik08+wfDhw5GXl4ewsDB8//33smORDbCciYgcmFarxfLlyzFs2DDk5eUhNDQUO3fulB2LaojlTETk4LRaLT777DM88cQT+OuvvxASEoIffvhBdiyqAZYzEZETcHNzw4oVKzB06FD8+eefCAkJAXdhdFwsZyIiJ+Hm5oaVK1di8ODB+OOPPxAcHIy9e/fKjkXVwHImInIi7u7uWLVqFaKjo3HlyhUEBQUhOztbdiyqIpYzEZGTcXd3x+rVqxEVFYXLly8jKCgIP/74o+xYVAUsZyIiJ+Th4YE1a9YgIiICly5dQmBgIPbt2yc7FlmJ5UxEtcJsNiMpKQnx8fFISkqC2WyWHcnpeXh4YN26dRgwYAAuXryIwFIuTiMAAA3FSURBVMBAHDhwQHYssgLLmYjszmw2IzQ8GjETYxGbcBUxE2MRGh7Ngq4FderUwfr16xEWFoYLFy4gICAABw8elB2LKsFyJiK7MxqN2HngNPIMO6D0XoA8ww7sPJALo9EoO5pLqFu3LhISEhASEoLz588jICAAhw4dkh2LKsByJiK7y87ORr5XCKBxv/EJjTvyvUKRk5MjN5gLqVu3LjZu3Ijg4GCcO3cOAQEBOHz4sOxYdBcsZyKyO19fX+jOpQKWohufsBRBdy4FPj4+coO5mHr16mHTpk0IDAzE2bNnERAQgJ9++kl2LCoHy5mI7M5gMMC/mzc8jf4Qu2bB0+gP/26tYTAYZEdzOfXq1cPmzZuh1+vx22+/Qa/X4+eff5Ydi8oQiqJI+cG9evVSuLUckeswm80wGo3IycmBj48PDAYDtFqt7FguKz8/H+Hh4di2bRu8vb2RlZWFjh07yo7l1IQQexRF6WXVa1nORESuKT8/HwaDAd988w1at26Nbdu2oUOHDrJjOa2qlDMfaxMRuSidTofk5GQ88sgjyM3NhV6vx7Fjx2THIrCciYhcmqenJ5KTk/HQQw/h5MmT0Ov1OH78uOxYLo/lTETk4ho0aACj0Yi+ffvixIkT0Ov1OHnypOxYLo3lTEREaNiwIbZu3Yo+ffrg+PHj6N+/P06dOiU7lstiORMREQCgUaNGSElJQe/evXHs2DHo9Xrk5ubKjuWSWM5ERHRL48aNkZKSgp49e+KXX35BQEAATp8+LTuWy2E5ExHRbe655x6kpaXB19cXP//8MwICAvDbb7/JjuVSWM5ERHSHJk2aID09HT4+Pvjpp5+g1+vx+++/y47lMljORERUriZNmiAtLQ09evTAkSNHoNfrcfbsWdmxXALLmYiI7qpZs2ZIT09Ht27dcPjwYQQEBODcuXOyYzk9ljMREVWoefPmyMjIQNeuXXHo0CEWdC1gORMRUaW8vLyQmZmJLl264ODBgwgKCsKFCxdkx3JaLGciIrJKSUF37twZ+/fvR2BgIC5evCg7llNiORMRkdVatGiBzMxM3H///di3bx+CgoJw6dIl2bGcDsuZiIiqpGXLljCZTOjUqRNycnIQFBSEy5cvy47lVFjORERUZa1atYLJZELHjh2RnZ2N4OBgXLlyRXYsp8FyJiKiavH29obJZEKHDh2wZ88ehISE4I8//pAdyymwnImIqNpat24Nk8mE9u3bY9euXQgNDWVB2wDLmYiIaqRt27YwmUxo164ddu7cCYPBgL/++kt2LIfGciYiohq79957YTKZ0LZtW3z//fcs6BpiORMRkU20a9cOWVlZaNOmDb799lsMGDAAeXl5smM5JJYzERHZTPv27WEymdC6dWts374d4eHhyM/Plx3L4bCciYjIpu677z6YTCa0atUKX3/9NSIiInD16lXZsRwKy5mIiGyuY8eOyPp/7d1tjJX1mcfx7wXDoCiiAj6kurVNUYtVYx0Ja3RdQ62U6BBQEKIWqtbYSk3sWlOjcRf0jbtqY2MNldBYi4oIIoM6DjgDabNZTWmaNrWNCbVPpJuUbbu+qYoDV1+cky47BblVOPd/5nw/yST3Pecfzo8r58yP++EMW7dy4oknsnXrVi6//HIL+n2wnCVJh8SUKVPYsmULJ5xwAgMDA3R3d/PWW2/VHWtYsJwlSYfMaaedxpYtWzj++OPp7+9n9uzZFnQFlrMk6ZA6/fTTGRgY4LjjjmPz5s3MmTOHt99+u+5YRbOcJUmH3NSpUxkYGGDy5Mn09fUxd+5c3nnnnbpjFctyliS1xBlnnEF/fz+TJk2it7eXK664woLeD8tZktQyZ555Jv39/UycOJEXXniBefPmsWvXrrpjFcdyliS11FlnnUV/fz/HHnssGzduZP78+Rb0EJazJKnlzj77bF5++WWOOeYYNmzYwIIFC3j33XfrjlUMy1mSVItzzjmHzZs3c/TRR7N+/XoWLlxoQTdZzpKk2px77rls2rSJCRMmsG7dOq6++moGBwfrjlU7y1mSVKvzzjuPTZs2cdRRR/HMM89w7bXXtn1BW86SpNpNmzaNvr4+xo8fz+rVq1m0aBG7d++uO1ZtLGdJUhGmT5/OSy+9xJFHHsmTTz7J4sWL27agLWdJUjHOP/98ent7OeKII1i1ahXXXXddWxa05SxJKsoFF1xAb28v48aN4/HHH+eGG25gz549dcdqKctZklScCy+8kBdffJFx48bx2GOPceONN7ZVQVvOkqQiXXTRRTz//PMcfvjhrFy5kptuuqltCrpSOUfEzIh4PSK2R8TX9/H42Ih4uvn4qxFxysEOKklqPxdffDEbN27ksMMOY8WKFdx8881kZt2xDrkDlnNEjAa+BXwOmAosjIipQ5ZdD/w5Mz8BfAO472AHlSS1pxkzZtDT08PYsWNZvnw5S5YsGfEFXeXIeRqwPTPfyMxdwGpg9pA1s4HvNrfXAjMiIg5eTElSO7vkkkvYsGEDY8eO5ZFHHuGWW24Z0QVdpZw/Avxur/0dze/tc01mDgJvAhMPRkBJkgAuvfRSnnvuOTo7O3n44Ye59dZbR2xBVynnfR0BD51GlTVExI0RsS0itu3cubNKPkmS/mbmzJmsX7+ezs5OHnroIW677bYRWdBVynkHcPJe+ycBv9/fmojoACYAfxr6B2Xmo5nZlZldkydP/mCJJUltbdasWaxbt44xY8bw4IMPcvvtt4+4gq5Szj8EpkTExyKiE1gA9AxZ0wMsam5fCQzkSJuUJKkYl112GWvXrmXMmDHcf//93HHHHSOqoA9Yzs1ryEuAPuAXwJrMfC0ilkVEd3PZSmBiRGwHvgr83cetJEk6mLq7u1mzZg0dHR3cd9993HnnnSOmoKOuv0hXV1du27atlueWJI0czz77LPPnz2f37t3cddddLFu2jBI/MBQRP8rMripr/Q1hkqRhbe7cuTz11FOMHj2ae++9l6VLl9Yd6UOznCVJw968efN44oknGDVqFEuXLuWee+6pO9KHYjlLkkaEq666ilWrVjFq1CjuvvtuVqxYUXekD6yj7gCSJB0sCxcuZM+ePTzwwAPMmTOn7jgfmDeESZJGnMHBQTo6yjr+9IYwSVJbK62Y3y/LWZKkwljOkiQVxnKWJKkwlrMkSYWxnCVJKozlLElSYSxnSZIKYzlLklQYy1mSpMJYzpIkFcZyliSpMJazJEmFsZwlSSqM5SxJUmEsZ0mSCmM5S5JUGMtZkqTCWM6SJBXGcpYkqTCWsyRJhbGcJUkqjOUsSVJhLGdJkgpjOUuSVBjLWZKkwkRm1vPEETuB39Ty5Ps3CfifukMME86qGudUjXOqzllVU+KcPpqZk6ssrK2cSxQR2zKzq+4cw4GzqsY5VeOcqnNW1Qz3OXlaW5KkwljOkiQVxnL+/x6tO8Aw4qyqcU7VOKfqnFU1w3pOXnOWJKkwHjlLklQYy1mSpMK0ZTlHxMyIeD0itkfE1/fx+NiIeLr5+KsRcUrrU9avwpy+GhE/j4ifRkR/RHy0jpwlONCs9lp3ZURkRAzbj3h8GFXmFBHzm6+r1yLiyVZnLEWF998/RMSWiPhx8z04q46cdYqI70TEHyLiZ/t5PCLim80Z/jQiPt3qjB9YZrbVFzAa+CXwcaAT+AkwdciaLwPLm9sLgKfrzl3onC4GxjW3v9SOc6o6q+a68cD3gVeArrpzlzgnYArwY+CY5v5xdecueFaPAl9qbk8Ffl137hrm9E/Ap4Gf7efxWUAvEMB04NW6M1f9ascj52nA9sx8IzN3AauB2UPWzAa+29xeC8yIiGhhxhIccE6ZuSUz/9LcfQU4qcUZS1HlNQVwD/DvwNutDFeQKnP6IvCtzPwzQGb+ocUZS1FlVgkc1dyeAPy+hfmKkJnfB/70HktmA49nwyvA0RFxYmvSfTjtWM4fAX631/6O5vf2uSYzB4E3gYktSVeOKnPa2/U0/oXajg44q4g4Bzg5M59vZbDCVHlNnQqcGhH/GRGvRMTMlqUrS5VZ/RtwTUTsAF4EvtKaaMPK+/05VoyOugPUYF9HwEM/T1ZlzUhXeQYRcQ3QBVx0SBOV6z1nFRGjgG8Ai1sVqFBVXlMdNE5t/zONMzE/iIhPZeb/HuJspakyq4XAY5n5QET8I/C95qz2HPp4w8aw/VnejkfOO4CT99o/ib8/HfS3NRHRQeOU0XudOhmJqsyJiPgMcCfQnZnvtChbaQ40q/HAp4CtEfFrGte+etrwprCq770NmfluZv4KeJ1GWbebKrO6HlgDkJn/BRxG4z970P+p9HOsRO1Yzj8EpkTExyKik8YNXz1D1vQAi5rbVwID2by7oI0ccE7NU7XfplHM7XptEA4wq8x8MzMnZeYpmXkKjevz3Zm5rZ64tany3nuOxo2GRMQkGqe532hpyjJUmdVvgRkAEfFJGuW8s6Upy9cDfL551/Z04M3M/O+6Q1XRdqe1M3MwIpYAfTTuiPxOZr4WEcuAbZnZA6ykcYpoO40j5gX1Ja5HxTn9B3Ak8EzzfrnfZmZ3baFrUnFWba/inPqAz0bEz4HdwNcy84/1pa5HxVn9C7AiIm6lcap2cbsdRETEUzQugUxqXnv/V2AMQGYup3EtfhawHfgL8IV6kr5//vpOSZIK046ntSVJKprlLElSYSxnSZIKYzlLklQYy1mSpMJYzpIkFcZyliSpMH8Fs2uBPuFDFugAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x4ea20278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This file is associated with the book\n",
    "# \"Machine Learning Refined\", Cambridge University Press, 2016.\n",
    "# by Jeremy Watt, Reza Borhani, and Aggelos Katsaggelos.\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import csv\n",
    "\n",
    "# sigmoid for softmax/logistic regression minimization\n",
    "def sigmoid(z): \n",
    "    y = 1/(1+np.exp(-z))\n",
    "    return y\n",
    "    \n",
    "# import training data \n",
    "def load_data(csvname):\n",
    "    # load in data\n",
    "    reader = csv.reader(open(csvname, \"r\"), delimiter=\",\")\n",
    "    d = list(reader)\n",
    "\n",
    "    # import data and reshape appropriately\n",
    "    data = np.array(d).astype(\"float\")\n",
    "    X = data[:,0:2]\n",
    "    y = data[:,2]\n",
    "    y.shape = (len(y),1)\n",
    "    \n",
    "    # pad data with ones for more compact gradient computation\n",
    "    o = np.ones((np.shape(X)[0],1))\n",
    "    X = np.concatenate((o,X),axis = 1)\n",
    "    X = X.T\n",
    "    \n",
    "    return X,y\n",
    "\n",
    "# YOUR CODE GOES HERE - create a gradient descent function for softmax \n",
    "# cost/logistic regression \n",
    "def softmax_grad(X,y):\n",
    "    w1 = w0\n",
    "    alpha = alpha0\n",
    "    w_path = [] \n",
    "    max_its = 2000\n",
    "    for k in range(max_its):  \n",
    "        t1 = y*np.dot(X.T,w1)\n",
    "        grad = -np.dot(X,sigmoid(-t1)*y)\n",
    "        w1 = w1 - alpha*grad\n",
    "        \n",
    "        w_path.append(w1)\n",
    "    w_path = np.asarray(w_path)\n",
    "    w_path.shape = (np.shape(w_path)[0],3)\n",
    "    return w_path\n",
    "\n",
    "# plots everything \n",
    "def plot_all(X,y,w):\n",
    "    # custom colors for plotting points\n",
    "    red = [1,0,0.4]  \n",
    "    blue = [0,0.4,1]\n",
    "    \n",
    "    # scatter plot points\n",
    "    fig = plt.figure(figsize = (8,8))\n",
    "    ind = np.argwhere(y==1)\n",
    "    ind = [s[0] for s in ind]\n",
    "    plt.scatter(X[1,ind],X[2,ind],color = red,edgecolor = 'k',s = 25)\n",
    "    ind = np.argwhere(y==-1)\n",
    "    ind = [s[0] for s in ind]\n",
    "    plt.scatter(X[1,ind],X[2,ind],color = blue,edgecolor = 'k',s = 25)\n",
    "    plt.grid('off')\n",
    "    \n",
    "    # plot separator\n",
    "    s = np.linspace(0,1,100) \n",
    "    plt.plot(s,(-w[0]-w[1]*s)/w[2],color = 'k',linewidth = 2)\n",
    "    \n",
    "    # clean up plot\n",
    "    plt.xlim([-0.1,1.1])\n",
    "    plt.ylim([-0.1,1.1])\n",
    "    plt.show()\n",
    "    \n",
    "# load in data\n",
    "X,y = load_data('/Users/Administrator/Desktop/python/imbalanced_2class.csv')\n",
    "\n",
    "w0 = np.array([0,2,2])\n",
    "w0.shape = (3,1)\n",
    "alpha0 = 10**(-1)\n",
    "\n",
    "# run gradient descent\n",
    "w_final = softmax_grad(X,y)\n",
    "w=w_final[-1]\n",
    "# plot points and separator\n",
    "plot_all(X,y,w)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4_9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAGKCAYAAAARyQg4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xl8FdX5+PHPk30hhC0JOwnbJUEUJCKKIEtAEkSwamvrT0WlWJeqra3la6vSarVUqa1LVRRF/bpb+eIGCgiyyCIKKhpk3yEkIJsJkOX5/TFDmoQkzA25WZ/36zWve2fmnJnn3qDPnTPnnBFVxRhjjDENR1BtB2CMMcaY6mXJ3RhjjGlgLLkbY4wxDYwld2OMMaaBseRujDHGNDCW3I0xxpgGxpK7afBEZJKIaIklV0S+EZEJtR1bWSLyUxEZ50f5piLyFxH5TkTyROSwiCx0j9Mg/vsu87fLE5FtIvKOiIyu4vH6icikag7Tn/PHu/8mE2srBtPwNYj/+I3x4CBwnruMBuYBz4jIL2o1qpP9FBjnpaCIxAPLgPHA80AGcIW77Tmcz9lQTMH5240AJgLHgZki8nwVjtUPuK8aY/NXvHv+xFqMwTRwIbUdgDE1pEBVl5VYnyci5wNjgVcrqiQikaqaF/DoquYpoDmQqqo7S2yfLSJPALG1E5b/PHzPW8r8/V4VkY+BaSLyqaq+GOAQjalX7MrdNGaHgdATKyIy2G36vUhE3hWRI8AT7r4gEZkoIhtE5JiIrBORa0seTERGicgcEdkrIodEZJmIjChTpr2IvOmWyRORjSJyv7tvOnAZcGGJZuhJ5QUuIp2AS4EHyyR2AFR1m6p+45btISKvi8h295bEtyJyR8lm+xKffbCIvCUiR0Rkk4jcXM65B4nIfLfMQRFZICJ9Suzv6J5vv3u+j0TEV2J/onuuq0TkJRE5ALxX4V+pAqr6PLAcuKnEsc9z/3a7RORHEVktIleV2D8OeNx9f+I7XuDH9xQqIo+4twaOueeZISJhXj6/2xT/jVt0/okY/P3sxpyKXbmbRkNETvx7jwIuAS4Eri+n6DTgBeCfwFF32+PAtcBfgC+B4cDzIrJPVd93yyThJKlHgCIgHZglIoNUdYlb5iUgEpgAHAA6Az3cffcDHYFmwImkuqOCjzMIEGC2h4/eDvgeeAXnB01v4M9uHA+VKfss8CIwFfg58KSIrFTVFeD8CADmAPPd7+NHYIB7jlUi0gJYDOwDfgXk4jSjzxWR7mWuzh8B3sG5lVDo4XOUZw4wUURCVTUf6AQsAZ7G+dsNAF4QkSJVfQ34AKeJ/06cZn6AQ358T/8DXOV+ps1Aa5zbIcHu91Pp5wd2u/VfAW7B+bdkTPVTVVtsadALMAnQcpZ/lSk32N3+aJntXXGS9bVltr8EfF7BOYNwfjx/BDxfYvsRYHQlsb4NLPDwmSa6sYb7+V2IG9fdwKZyPvtfSmwLBbKBv5XYthRYCUgFx78fJ7G1KLGtOU6fh1vc9UT3XDM8xqzArRXsu9Hdn1DJZ30G+KTE9lud//VV6Xt6H5hSST0vn/8MN+bBtf3fhi0Nd7Erd9NYHATS3PfhQF/gLyKyX1X/XKbsB2XWh+Ek9xklrv7B6ZT3cxEJVtVCEWkP/NU9TxucBAHOleQJq4GHRKQlTsLZdpqf65RNuiISwX+vODtS+lZEiKoWlCj+cfGBVfNFZD3Q3i0bDZwL3K6qFZ03Dedq+lCJ7+ow8AWQWqZs2e+5KqTUikhznKvtMThX4sHurpNuXZx0IG/f02rgJhHJwmk1+abMd+HP5zcmYOyeu2ksClR1pbssUdXHcK6y7nabUkvKKrPeCidJHATySyzTca7u2rj3Zd8FzgfuBYYA5wCzgIgSx/oZzpXvo8BW957wsCp8nhPJqqOHspOB3+E0tWe4cT3g7osoU/ZAmfXjJco0x0mmuys5Vyucz5hfZhkCdChTtuz3XBXt3OPvd9enu+d/GKdn/Tk4IwnKfs7yePmeHgCexLlt8hWwXURuL3EMfz6/MQFjV+6mMfsOCAO68N/kACdfDe8HCnDu3xaVc5y9OE33fYB0VS2+Dy4ikSULqtP5bZz7Y6Afzi2Dd0Wko6ru8yP2hW6cFwEbTlH2CuBxVf17ibhG+XGuE37A+fxtKimzH+dHzv3l7DtcZr06OpKNAL5wWxkigFE4TfhPnygg3sf7n/J7UtWjOD/e7hWRbjj31f8pIt+7f3d/Pr8xAWPJ3TRmZ7iv209R7hOcK/dYVZ1TXoESSfxYiW2dcH4QfF22vKoWActE5M/AZzgdwfZR+kq5Qqq6VURm4LQ8vKOqpa6mRaQD0EydHvORZeIKBq481TnKOeePIrIcuEZEnqigaX4ezlj9bzXAQwhF5HqcH0gnRi2E4/ydSn7WGJzOkyVjPe7ui3CT9Ql+fU+qul5EfofTMS4Fp5ney+c/7r56aU0wpkosuZvGIkRE+rvvw3Duuf8JmKmqeyqrqKrfi8jTwOsi8necZvUIoCfQXVXHA2txerZPEZF7gBice7/F93pFJBang91LwDqcZHQnsAfIdIutBcaIyFj3eLtUdVcFod0EfAqsFJF/4NzXDccZBXALcA3OsKs5wC0isgHnyvIWt1xVTATm4owCmIrTW/48YKU6owb+Afw/4BMRedz9/AluTIvV6bFeFYnu3y8Upw/AGJwk+ryqvgSgqgdF5HOcq+pDOK0ME3FupzQtcay17uvtIvIJcEhVv8fD9+T+oPoCWAXkAZfj/H90oVvEy+ff5ta9VkQOAvmqurKK34sx5avtHn222BLohZN7yx8H1uPcY40pUW6wu/+Mco4hwB3AtzhXd9k4ifWaEmXOAVbg/I97Pc5Mc9NxEh84ieJZnOFWuUAOTu/rXiWO0QqYgZNcFJh0is8Wi9MEvBZn6NdhnERzPRDslklwj3kI5z7334FfusdvUtlnBxYAb5fZdqF7jlyce/Tzgd4l9rfFGUqY5X5XW4D/BXq6+xPdc13s8e9X8m93FKel5R3KGXWAc3vkE5wfHduAu9y/f06Zv+XfgV04PwAW+PE9/R7nx91B97teDowpE0Oln98tcxXOD7zjnKLnvi22VGURVZs/wRhjjGlIrLe8McYY08BYcjfGGGMaGEvuxhhjTANjyd0YY4xpYCy5G2OMMQ1MvR3n3qpVK01MTKztMIwxxpga8cUXX+SoapyXsvU2uScmJrJypc37YIwxpnEQka1ey1qzvDHGGNPAWHI3xhhjGhhL7sYYY0wDY8ndGGOMaWAsuRtjjDENTL3tLW+MMVVx6NAh9u7dS35+fm2HYkyx0NBQ4uPjadq06akLe2DJ3RjTaBw6dIisrCzatWtHZGQkIlLbIRmDqpKXl8fOnTsBqiXBW7O8MabR2Lt3L+3atSMqKsoSu6kzRISoqCjatWvH3r17q+WYltyNMY1Gfn4+kZGRtR2GMeWKjIystttFNZrcRSRCRFaIyFci8q2I/Nnd/oqIfC8ia0TkeREJrcm4jDGNh12xm7qqOv9t1vSV+zFgqKqeBfQGRopIf+AVoAfQC4gExtdwXMYYY0yDUaPJXR1H3NVQd1FV/dDdp8AKoH1NxnVw7zZWzX29Jk9pjDHGBEyN33MXkWARWQ3sBeao6vIS+0KBq4HZFdSdICIrRWRldnZ2tcW05j+T6bnoZrKqqSODMcYE0vTp0+nbty8xMTE0b96cPn368Nvf/tbv4/z4449ceeWVtGzZEhFh+vTpfPzxx/zzn/8MQNSmJtV4clfVQlXtjXN13k9Eziix+9/AQlVdVEHdqaqaqqqpcXGennrnSdeBVxAmhXz20RvVdkxjjAmEhx56iPHjx3PRRRfxzjvv8NJLLzFmzBjeffddv4/11FNP8d577zF16lSWLl3KqFGjLLk3ELU2zl1VD4jIAmAksEZE7gPigBtrOpaElIEcDo4ldMNsDh/9FTER1p/PGFM3PfHEE9x44408+OCDxdtGjx7Nfffd5/ex1q5di8/n47LLLqvOEE0dUNO95eNEpJn7PhJIA9aKyHjgIuDnqlpUkzEBEBRMfucRXMAq3li2ucZPb4wxXh04cIDWrVuftL1sT+ucnByuvfZaWrZsSVRUFIMHD2blypXF+xMTE5k2bRqrVq1CRBARJk2axJQpU9i6dWvxtnHjxgEwbtw4UlNT+eCDD0hJSSEqKopRo0axf/9+NmzYwJAhQ4iOjiY1NZWvv/66VCxTpkzhnHPOITY2loSEBEaPHs2GDRuK92/evJmYmBjuvffe4m35+fmcffbZDBo0iKKimk8L9V1NN8u3AeaLyNfA5zj33N8HngYSgKUislpE7q3sIIHQ4uwxNJMf+XLJbPIL7R+SMaZuOvvss3n88cd58cUX2bdvX4Xlxo4dy0cffcQjjzzCG2+8QVFREUOGDClOqjNmzCAjI4MePXqwdOlSli5dyvjx4/nFL35B69ati7fdc889xcfctm0b9957Lw888ABTp07ls88+Y8KECVx55ZVceeWVvP322xQUFHDllVfi9I927Nixg1tvvZWZM2fy7LPPUlhYyIABAzh48CAASUlJPPLIIzz00EPFP0AeeOAB1q1bxwsvvEBQkE3J4q8abZZX1a+BPuVsr/1pcDsPoSgolD55S3nvq5/wk7NrtMO+MaaW/Pm9b/lu16FaOXdK26bcN7qnX3WefPJJxo4dy7hx4xARkpOTueyyy/jd735XPG3p7NmzWbJkCQsWLODCCy8EYOjQoSQmJvLwww/zzDPP0KdPH+Li4sjKyqJ///7Fx2/Tpg3h4eGltp2wf/9+li5dSpcuXQD4+uuvefjhh3nxxRe55pprAGcq1VGjRrF27VqSk5MBePTRR4uPUVhYyPDhw4mPj2fmzJnF9W688UZmzJjBtddey7PPPsuDDz7IY489Vnwu4x/7OXRCeBOk82AyQr9k6qcbS/3qNMaYuuLMM88kMzOTd999l5tvvhlV5f777yc1NZUjR5yRxitWrCAuLq44sQNER0dz8cUXs3jx4iqfOzExsVSy7dq1K+D8cCi77cQ86QDLli1j+PDhtGzZkpCQEKKiojhy5Ajr1q0rdfxp06axa9cuhg4dypAhQ7jpppuqHGtjV/tXzHWI+NJpt2EO+Xu/Z9H6FAZ1r74e+caYusnfK+e6IDw8nNGjRzN69GjASYrjx49n2rRp3H777ezevZuEhIST6iUkJLB///4qn7dZs2al1sPCwk7afmLb0aNHAacpf8SIEfTr149nnnmGtm3bEhYWxqhRo4rLnNCuXTsGDhzIe++9Z4n9NNmVe0ndRwJwaeRXTF24qZaDMcYYb2644QZatGjB2rVrAadpvbwHkGRlZdGiRYsajW327Nnk5uYyc+ZMLr/8cs4//3x69+5d7o+MGTNm8P7773PWWWdx1113kZubW6OxNiSW3EuKbQdtenNZk29YvCGHNTsP1nZExhhTSnlJOzs7m4MHDxZfrZ977rns3buXhQsXFpfJzc3lgw8+4IILLqj0+GFhYSddUZ+OvLw8goKCCAn5b0Pxm2++SUFBQalyOTk5/OpXv+Lmm29m9uzZ7N+/nz/84Q/VFkdjY8m9LF8GrQ99TYewIzy7yK7ejTF1S69evZgwYQJvv/02Cxcu5OWXXyYtLY2oqCiuvfZaAC666CIGDBjAz372M1588UXef/99MjIyyMvL4/e//32lx+/RowdZWVlMnz6dlStXsmXLltOKd+jQoRQWFnLdddcxb948HnvsMSZOnHhSE/9NN91ETEwMkydPpnXr1jz55JM8+eSTzJ8//7TO31hZci/Ll46g/L7zVt7/ejc7frBmIWNM3XHvvfeyZcsWbrvtNkaMGME999xDz549WbFiBUlJScXlZsyYwfDhw7njjju44oorUFU++eST4g5vFfnpT3/KuHHjuOuuuzjnnHOYNGnSacXbq1cvXnjhBZYvX87FF1/Mq6++yltvvUVsbGxxmddee4133nmH6dOnEx0dDcCVV17J5ZdfznXXXcfhw4dPK4bGSOprr/DU1FQtOSFDtVGFR88gL+5Mzsi8mnHnJ3LPxSnVfx5jTI3LzMwsHp5lTF1U2b9REflCVVO9HMeu3MsSAd9IIrct4NIzWvD6im0czMuv7aiMMcYYzyy5l8eXDvm5/Lrzbn48Xsgry7fWdkTGGGOMZ5bcy5M4EMKa0Cn7Uy7o2ooXlmzhWEFhbUdljDHGeGLJvTwh4dB1GKybzYSBiWQfPsbMVbtqOypjjDHGE0vuFfFlwOHdDGyyk+Q2TZm6aBNFRfWz86ExxpjGxZJ7RbqNAAlCvp/FhEFJbNh7hAXrTp48whhjjKlrLLlXJKoFdOgP38/i4jPb0iY2gmc+tUltjDHG1H2W3CvjS4esbwg9vJPrBySxfPN+vtp+oLajMsYYYyplyb0yvgzndd1sruzXgZjwEKbalLTGGGPqOEvulWnVFVp2g+8/JCYilF/078isb3azbZ9NSWuMqT3Tp0+nb9++xMTE0Lx5c/r06cNvf/vb2g6r0RIRnnjiidoOoxRL7qfiS4fNi+DoIa4fkERwkDBtsV29G2Nqx0MPPcT48eO56KKLeOedd3jppZcYM2YM7777bm2H1mgtXbqUK664orbDKMWS+6n4MqAoHzbOI6FpBGN6t+PNlTv44cfjtR2ZMaYReuKJJ7jxxht58MEHGT58OKNHj2bSpEmsX7++tkMrJT8/n8LC+jn5V15enl/l+/fvX/y43brCkvupdOgHkS3g+1kA/HJgZ/LyC3l5mU1Ja4ypeQcOHKB169YnbReRUuvbt28nIyODyMhIEhMTee6557j88ssZPHhwcZlx48aRmlr6OSRbtmxBRHj//feLt02ZMoVzzjmH2NhYEhISGD16NBs2bChVb/DgwVx++eVMnTqVLl26EBERwa5dzuRfa9asYdSoUcTExBATE8MVV1zBnj17PH/mE8d+4YUXSEpKokmTJlx99dUcO3aMFStW0K9fP5o0acLgwYPZtm1bqboTJ06kV69eNGnShPbt23PVVVeddO7ExETuvPNO7r//ftq3b0/Tpk2L9z3xxBN06NCB6Ohoxo4dy7x58xARFixYUOq7L9ksfyLeV199la5du9K0aVPS09PZsWOH5898ukJq7Ez1VVAwdL/ISe6FBfhaxzDYF8eLn21hwqDORIQG13aExphG5Oyzz+bxxx+nY8eOXHzxxbRs2fKkMqrKmDFjyMnJYdq0aURERHDfffexf/9+unXr5vc5d+zYwa233kqnTp04dOgQTz/9NAMGDGDdunWlHt26ZMkSNm7cyOTJk4mKiiI2NpYNGzYwYMAAUlNTefnllyksLOSee+5h9OjRrFix4qQfJRVZtmwZOTk5PP7442zbto3f/OY3REZGsnz5cu666y6io6O57bbbmDBhArNnzy6ut3fvXu6++27atm1LdnY2U6ZMYejQoXzzzTcEB//3/9+vvvoqPXv25N///jcFBQWA89jcX//619x8882MGTOGxYsXc8MNN3iKd/ny5ezatYspU6aQl5fH7bffzoQJE/jwww891T9dlty98KXDV6/B9uWQOIAJgzrzi2eX886XO/nFuR1rOzpjzOmYNRH2fFM7527dC9L/5leVJ598krFjxzJu3DhEhOTkZC677DJ+97vfFV9xzpo1i1WrVrFs2TLOPfdcAPr27UuXLl2qlNwfffTR4veFhYUMHz6c+Ph4Zs6cyTXXXFO878CBA6xatapUy8Itt9xC69atmTVrFmFhYQCceeaZ9OjRgw8//JBRo0Z5iuHIkSPMnDmz+MfEggULePbZZ/n0008ZNGgQALt27eKWW24hNzeXqKgoAJ5//vlSsZ933nm0b9+eJUuWFNc74f333yciIqJ4/cEHHyQjI4Mnn3wSgBEjRpCTk8NTTz11yngPHTrEBx98QPPmzQHYs2cPv/nNb8jLyyMyMtLTZz4d1izvRZehEBwG3zu/uM7r3JJe7WJ5zqakNcbUsDPPPJPMzEzeffddbr75ZlSV+++/n9TUVI4cOQLAihUrSEhIKE7sAJ06daJv375VOueyZcsYPnw4LVu2JCQkhKioKI4cOcK6detKlevbt+9Jtwzmzp3LpZdeSlBQEAUFBRQUFJCUlERiYiIrV670HENqamqpVoKuXbsSFhbGBRdcUGobUHw7AJwfOueffz6xsbGEhITQvn17gJNiHzZsWKnEXlhYyOrVq7nkkktKlSu7XpFzzjmnOLEDpKSkALBz505P9U+XXbl7ER4DSYOc5D7iAUSECYM68+vXVjEnM4uLep58/8sYU0/4eeVcF4SHhzN69GhGjx4NwLRp0xg/fjzTpk3j9ttvZ8+ePcTHx59ULz4+nsOHD/t1rm3btjFixAj69evHM888Q9u2bQkLC2PUqFEcPXq0VNnyOpXl5OQwefJkJk+efNK+7du3e46jWbNmpdbDwsKIiYkhKCio1DagOK7PP/+cSy65hEsvvZSJEycSHx+PiNC/f/9Txp6dnU1BQQFxcXGltpdd9yfekrEFmiV3r3zp8MGdkLMe4rqTfkZr2jePZOrCTZbcjTG16oYbbuCuu+5i7dq1ALRu3Zq9e09+FsbevXtLNQlHRERw/HjpkT/79+8vtT579mxyc3OZOXMm0dHRABQUFJxUDk7u1AfQokULLr30UsaPH3/SvlatWnn4dFU3Y8YM4uLieOONN4pj27q1/M7QZWOPi4sjJCSE7OzsUtvLrtdV1izvVfd059Vtmg8JDuKGC5L4YusPfLH15H/kxhgTCOUl7ezsbA4ePFh89XnOOeeQlZXF8uXLi8ts27aNL7/8slS99u3bs2XLllJXk3PmzClVJi8vj6CgIEJC/nst+OabbxZ3OjuVYcOGsWbNGvr27UtqamqpJTEx0dMxqiovL4/Q0NBSifuVV17xVDc4OJjevXszc+bMUtvry3wClty9im0Hbc4qHhIH8NPUDsRGhjJ1oU1qY4ypGb169WLChAm8/fbbLFy4kJdffpm0tDSioqK49tprAcjIyOCss87iiiuu4LXXXmPGjBmMGjXqpKb6sWPHcuTIEcaPH8/cuXN5+OGHeeGFF0qVGTp0KIWFhVx33XXMmzePxx57jIkTJ57U7FyRSZMmFQ+Fe/vtt1mwYAGvvPIK48aNKzWcLBCGDx/Otm3buOOOO5g3bx73338/L774ouf6d999N7NmzeLWW2/l448/5r777uODDz4AKHU7oC6q29HVNd3TnR7zP+YAEB0ewtX9O/Hxd1lsyj5Sy8EZYxqDe++9ly1btnDbbbcxYsQI7rnnHnr27MmKFStISkoCnCbmd999l5SUFK6//nruuOMObr31Vs4777xSxzrjjDN4/vnnWbp0KZdccgmffvppqd7l4PyYeOGFF1i+fDkXX3wxr776Km+99Vapzm2V6d69O8uWLSMqKooJEyaQnp7OfffdR3h4eHEHuEDJyMhg8uTJ/Oc//yn+fCXH75/KpZdeymOPPcb//d//MXbsWD7//HMeeeQRgFJj4esiUa2fvb1TU1PVn56W1WLXaph6IYz5N/S5CoDsw8cY8LdPuDy1PQ9e2qtm4zHG+CUzM5Pk5OTaDqPWXH755eTk5AT8irkhe+CBB/jrX//K/v37AzKkrbJ/oyLyhaqmlruzDOtQ5482Z0FMW1g3qzi5x8WE85Oz2/H2Fzv47fDutGoSXstBGmOMqQ7Z2dk89NBDDBkyhKioKBYtWsTkyZO54YYbamSs+umw5O4PEXdCm9ch/yiEOmMixw/szOufb+elpVv57fDutRykMcbUH4WFhVTWglyyI19NCwsLY+3atbz00kscPHiQNm3acPvtt3P//ffXWkxeWXL3ly8DVk6DLYug23AAusY3IS05gZeXbuGmC7sQGWZT0hpj6p633367tkM4ybBhw/j0008r3F+bt45jY2NrbLrY6mbJ3V9JAyGsiTMkzk3uADde2Jkrns7irS+2c815ibUXnzHG1CPPPPOM3xPrmFOz5O6vkHBnOtrvZ8GofzhN9UBqp+b06diM5xZt5qpzOxEc5O1hCMYY05j5fL7aDqFBsqFwVeFLh8O7Yffq4k0iwoSBndm2P5fZa7w/ytAYY4ypbpbcq6LbCJCgUhPaAIzo2ZrEllFMXbixVu8TGWMqZv9tmrqqOv9tWnKviuhW0OHck5J7cJBww8DOfLXjICs225S0xtQ1oaGh5OXl1XYYxpTrxHS51cFTcheRZBHpX2I9UkQeFJH/E5FfV0sk9Y0vHfZ8DQd3lNp8Rd/2tIgOsylpjamD4uPj2blzJ7m5uXYFb+oMVSU3N5edO3eW+zS/qvDaoe7fwGfAMnf9EWAcsAiYLCIRqvpwtURUX/gyYM69ztV7v18Wb44IDeaa8zrxz7nrWZ91mG4JMbUYpDGmpBNThu7atYv8/PxajsaY/woNDSUhIaHaprX1mtzPAKYAiEgo8P+AO1T1WRG5A7gRaFzJvVU3aNn1pOQOcHX/Tjy1YCPPLdrM5MvPrKUAjTHladq0aZ2fF9yY0+X1nns0cMh9399df8dd/xLo5OUgIhIhIitE5CsR+VZE/uxuv1VENoiIikhgH/BbnXzpsHkhHD1UanPLJuFckdqeGat2svfQ0QoqG2OMMYHhNblvwknqAJcCq1R1n7veCvA6A8ExYKiqngX0Bka69/KXAGnAVo/HqRt8GVCUDxs/OWnX+As6k19UxPTPttR8XMYYYxo1r8n9UeABEfkcuA14rMS+wcDXXg6ijhPPRg11F1XVVaq6xWMsdUf7fhDZ/KRe8wCJraIZ2bM1/7tsK0eOFdRCcMYYYxorT8ldVafhXFm/Dlykqi+X2L0f+KfXE4pIsIisBvYCc1R1uR/x1i3BIdDtIlj/ERSenMAnDOrMoaMFvPH59loIzhhjTGPleZy7qi5U1SmqOq/M9kmq+oEfxylU1d5Ae6CfiJzhta6ITBCRlSKyMjs722u1wPKlQ94PsGPFSbv6dGzOOYnNeX7xZgoKi2ohOGOMMY2RX5PYiEh3ERkqIhllF39PrKoHgAXASD/qTFXVVFVNjYuL8/eUgdF1GASHOQ+SKceEQV3YeSDv81qlAAAgAElEQVSPD77ZXcOBGWOMaay8TmKTIiLfAJnAXOD9Mst7Ho8TJyLN3PeROE39a6sQd90RHgOJA8u97w4wrEc8XeKimbpwk02aYYwxpkZ4vXJ/BggDfgL4gKQyS2ePx2kDzBeRr4HPce65vy8it4nIDpym+q9F5Dk/PkPt86XDvg2Qs/6kXUFBwi8HdubbXYf4bOO+ciobY4wx1ctrcu8D3KmqM1V1vapuLbt4OYiqfq2qfVT1TFU9Q1X/4m5/TFXbq2qIqrZV1fFV/UC1wpfuvFbQND+2TztaNQnnGZuS1hhjTA3wmtw3AhGBDKRei20PrXtV2DQfERrMuPM7sXBdNpm7D5VbxhhjjKkuXpP7ncDdIuK1+b3x8WXA9uXwY065u/9f/05EhQXz7CK7ejfGGBNYXpP7Q0A7YK2IrHOnkC21BDDG+sGXDloE6z8ud3ezqDB+mtqBd1fvYvdBe+SkMcaYwPGa3NcAHwKv4EwV+205S+PWpjfEtKmwaR7ghguSUOCFJVtqLCxjjDGNj6enwqnqdYEOpN4Tca7ev3oD8o9C6MldFDq0iCKjVxteXb6NW4d2pWlEaC0EaowxpqHzaxIbABFpJSLdRKRlIAKq13wZkP8jbFlcYZEJAztz5FgBr6/YVoOBGWOMaUw8J3cR+ZmIZAJZOBPP7BWRTBG5ImDR1TeJAyE0usIhcQC92sdyXueWPL94C8cLbEpaY4wx1c/rDHU/B17DefTrdUCG+7oJeF1ErgxYhPVJaAR0GeLcd69kNroJF3Zmz6GjvPfVrhoMzhhjTGPh9cr9j8BUVR2lqi+p6kfu6yjgWeBPgQuxnvFlwOFdsPurCosM7h6HLyGGZxfZlLTGGGOqn9fk3hX4TwX7/uPuNwDdLwKk0l7zIsL4gUms3XOYT9fVkafbGWOMaTC8JvcsILWCfanufgMQ3Qo6nAvrKk7uAGN6tyOhabhNamOMMabaeU3uLwCTRORPItJDRJqLiE9E/gTcBzwfuBDrIV+60yx/cGeFRcJCgrhuQBJLNuxjzc6DNRicMcaYhs5rcv8L8AgwEWfCmhzgO3f9EXe/OcHnPt7+FFfvvzi3I03CQ5hqD5QxxhhTjTwld1UtUtU/Ah2AwcDP3dcOqvontV5hpbXqBi26VHrfHaBpRCg/79eBD77ZzY4fcmsoOGOMMQ2dX5PYqOoPqrpIVd90X38IVGD12onZ6jYvhGOHKy163YAkBHh+8ZYaCc0YY0zDV2FyF5EMEWla4n2lS82FXE/40qHwOGz8pNJibZtFMvqstrz++TYO5ubXUHDGGGMassrmln8f6A+scN8rIBWUVSC4ekOr5zr0h4hmTtN8yphKi/5yYGdmrNrJ/y7fyi1DbFShMcaY01NZck8Cdpd4b/wRHOKMeV/3ERQVQlDFv31S2jZlYLdWTP9sC+MHJhEeYr+TjDHGVF2FzfKqulVVj59YBXa520otwE53vynLlw55+2H7qR93f+OgLmQfPsb/rap4+JwxxhjjhdcOdZuBPhXsO8vdb8rqMgyCQit9kMwJA7q2JKVNU55dtJmiIvutZIwxpuq8JveK7rUDRADHqiGWhieiKSQNPOWQOHCmpJ0wqDMb9h5h/vd7ayA4Y4wxDVWF99xF5Eygd4lNGSLSo0yxCOCnwLoAxNYwdE+HWb+HnPXO+PdKjDqzDX+fvZZnFm5iWHJCDQVojDGmoamsQ92lOFPLgnNP/d4Kym0GbqzOoBoU30gnuX8/65TJPTQ4iOsvSOKBDzJZvf0AvTs0q6EgjTHGNCSVNcs/CMQATXGa5Ye66yWXcFXtoqpzAx1ovdWsIyT08tQ0D3Blv47ERIQwdeHGAAdmjDGmoaqst3y+qv6oqkdUNUhVF7jrJRebdcULXzpsXwa5+09ZtEl4CFed24nZa/awdd+PNRCcMcaYhsav6WdFpL2IDLUZ6vzkSwctgvUfeyp+3YBEgoOEaYttEIIxxhj/VXbPvZiIxABvAiNObHJfS47ZsplXKtKmN8S0cYbEnXXlKYsnNI1gbO92vLlyO3ekdadFdFgNBGmMMaah8Hrl/hDQERiIk9gvxXkq3DScDnX9AxFcgxEUBN1HwoZ5UOBt1OCEQZ05ml/Ey0u3Bjg4Y4wxDY3X5J4B/BVY7q7vUtWFqjoBmAn8PhDBNSi+dDh+BLYs8lS8W0IMQ3vE89LSLRzNLwxsbMYYYxoUr8k9AdiuqoXAj0CLEvs+5L/N9aYiSYMgNMpzr3lwHiiz78fj/OfLHQEMzBhjTEPjNblvB1q579cDF5fYdy5wtDqDapBCI6HLUCe5q7fpZft3bsGZ7WN5btFmCm1KWmOMMR55Te5zgDT3/aPALSLymYjMB+4HXgpEcA2OLx0O7YQ9X3sqfmJK2s05PzLnu6wAB2eMMaah8Jrc/4A7W52qvgxchtOR7gfgVmBiQKJraLpdBIhfTfMje7amQ4tIm9TGGGOMZ56Su6rmqmpOifUZqnqVqv5EVZ9S1aLAhdiANImDDv38Su4hwUGMv6AzX247wMotp54ExxhjjPGU3EVkmIiMq2DfOBEZUq1RNWS+dNi9Gg56f277FantaRYVytSFmwIYmDHGmIbCa7P8X3F6zJenFc489MaL7unO67rZnqtEhYVwdf9OzMnMYlP2kQAFZowxpqHwmtx7Aisr2LcKSKmecBqBOB80T/KraR7gmvMSCQ0O4tlFNiWtMcaYynlN7gWUHtteUstqiqVxEAFfBmz+FI55vwqPiwnnsrPb858vd5B92Nssd8YYYxonr8l9MfB7ESk1ybm7fifgbdo14/ClQ+Fx2PiJX9XGD0wiv7CIl5duCUhYxhhjGgavyf2PQA9gg4g8LCK/FZGHcSa08WFD4fzTsT9ENPO7ab5LXBPSkhN4adlWco8XBCg4Y4wx9Z3XoXBfA+cAS4Crgcnu62Kgn6quCViEDVFwKHQbAes/giL/5o2/cVBnDuTm89ZKm5LWGGNM+Tw/z11Vv1fVn6tqa1UNdV+vUtV1Xo8hIhEiskJEvhKRb0Xkz+72JBFZLiLrReSNss3/DZIvHXL3wY7P/aqWmtiCszs247nFmygotOkFjDHGnMxzcq8mx4ChqnoW0BsYKSL9cVoCHlXVbjiz3t1Qw3HVvK7DICjUeca7nyYM6sL2/XnM/nZPAAIzxhhT34VUtENE3gT+R1U3uu8ro6r6s1OdTFUVONFFPNRdFBgK/MLd/iIwCXjqVMer1yJiIXGAc999+F/8qjo8JYGkVtFMXbiJUb3aICIBCtIYY0x9VNmVexxO8gWId9crWuK9nlBEgkVkNbAX54E0G4EDqnqih9gOoF0FdSeIyEoRWZmdne31lHWXLwNy1kHOBr+qBQcJ1w1I5OsdB1m753CAgjPGGFNfVZbcr8VJvKjqYFUdUtni9YSqWqiqvYH2QD8gubxiFdSdqqqpqpoaFxfn9ZR1V/eRzus6/3rNg/NAGYB5mfa0OGOMMaVVltw3A30AROQTEelRnSdW1QPAAqA/0ExETtwiaA/sqs5z1VnNO0HCGX4PiQOIbxrBWR2aMSdzbwACM8YYU59VltzzgCj3/WCg6emeTETiRKSZ+z4S5xnxmcB84HK32LXAzNM9V73hS4dtSyHX/ye+DU+O56vtB9h76GgAAjPGGFNfVdihDmfO+H+JyBx3/dcisruCsqqqf/BwvjbAiyISjPPD4k1VfV9EvgNeF5EH3PNO8xh//edLh4UPw/o5cNYp+ySWkpaSwCMfr2Pe2r38vF/HAAVojDGmvqksuf8SeBgYg3MPfBjOULbyKHDK5O5OhtOnnO2bcO6/Nz5t+kCT1s6QOD+Tuy8hhvbNI5mXmWXJ3RhjTLEKk7uqrgVGA4hIETBWVVfUVGCNRlAQdL8I1rwDBccgJNxzVREhLTmB11ZsI+94IZFhwQEM1BhjTH3hdRKbJGB1IANp1HwZcPwwbFnsd9W05ASOFRSxeENOAAIzxhhTH1WY3EUkqsRqNhAiIlEVLYEPtQHrfCGERFap13y/pBbEhIcw9zsbEmeMMcZR2ZX7YRE5cR/8CHD4FIupqtBI6DLUSe5a7hD/CoWFBHGhL455a7MoKvKvrjHGmIapsg511+NOYuO+t8wRSL50+P4D2PMNtDnTr6rDUxJ4/+vdrN5xgLM7Ng9QgMYYY+qLyjrUvVji/fQaiaYx634RIM7Vu5/JfXD3eIKDhLnfZVlyN8YY461DnYiEiEh4mW0jROQOETk7MKE1Mk3iof05VZqKNjYqlH6JLZhns9UZY4zBe2/5NyjxlDYRuQ2YDTwELBORiwMQW+PjGwm7VsEh/2ffHZYcz/dZh9m2LzcAgRljjKlPvCb3/kDJB4//HpiiqpHAc8AfqzuwRsmX4byum+131eEpCQDMtQfJGGNMo+c1ubcE9gCISC+gLfC0u+8tIKX6Q2uE4npA88QqDYnr1DKabvFNLLkbY4zxnNyzgET3/Uhgq6qe6EkfCRRVc1yNk4hz9b7pUzh2xO/qaSkJLN+8n4O5+QEIzhhjTH3hNbm/BUwWkYdx5pB/qcS+PsD66g6s0fKlQ+Ex2DTf76ppyQkUFikL1lnHOmOMacy8JveJwDNAD5yOdQ+W2NcXp8OdqQ4dz4OI2Co1zffu0IxWTcKYa73mjTGmUatsEptiqloA/KWCfT+p1ogau+BQ6DbC6VRXVAhB3h8GExwkDO0Rz6w1e8gvLCI02OtvN2OMMQ2J13Hu8SKSVGJdRGSCiPxTREYHLrxGqvtIyN0HO1b6XTUtOYHDRwv4fPP+AARmjDGmPvB6aTcd+E2J9T8D/8bpXDdDRMZVb1iNXNc0CApxnvHupwu6tSIsJIg51mveGGMaLa/J/WzgEwARCQJuAu5W1R7AX4E7AhNeIxXZDDoNqNJ996iwEC7o2oq5mVmonw+hMcYY0zB4Te6xwD73fV+gBfCKu/4J0LWa4zK+DMj5HvZtPHXZMtKSE9i+P491Wf4PpzPGGFP/eU3uO/jvRDWjgLWqutNdjwWOVndgjZ5vpPNahav3YcnxgM1WZ4wxjZXX5P488HcReQu4C5haYl9/ILO6A2v0midCfM8qJfeEphGc1T6WOd9ZcjfGmMbIU3JX1YeAX+NMQftr4LESu1vgzC9vqpsvHbYthVz/e76nJSewevsB9h62RhVjjGlsPA+EVtWXVPXXqjpNS/TUUtVflXz2u6lGvnTQQtgw1++qae6DZOavtQltjDGmsfFrlhP3ue6dRSSl7BKoABu1tmdDdHyVhsT1aB1Du2aRzPnOkrsxxjQ2nmaoE5FQnKb4a4HwCop5n0rNeBMU5HSsWzMDCo5DSJjnqiJCWnI8b6zcTt7xQiLD7M9jjDGNhdcr93uBi4EbAAFuBa4D5gFbAJulLlB8GXD8MGxd7HfVtJQEjuYXsWRDTgACM8YYU1d5Te4/BSYBb7rrK9x78COAxcCYAMRmAJIuhJDIKvWaPzepJU3CQ2xInDHGNDJek3sHYJ2qFuKMaW9eYt8rwGXVHZhxhUVBlyFOcvdzxrmwkCAu9MUxN3MvRUU2W50xxjQWXpP7bqCZ+34zMKjEvi7VGpE5mS8dDm6HrDV+Vx2enEDOkWN8teNAAAIzxhhTF3lN7guAge77Z4G7ReRVEXkBmALMDEBs5oRuFzmvVWiaH+yLIzhImGfPeDfGmEbDa3L/I/ASgKr+E/g90Ak4C3gcuC0g0RlHTAK0S61Scm8WFcY5ic3tvrsxxjQiXmeo26Oqa0qsP6qqA1T1bFX9g6r+GLgQDeA0ze/6Eg7t9rtqWnICa/ccZvv+3AAEZowxpq7xaxIbU4t8Gc7rutl+Vx2W7MxWZ1fvxhjTOFQ4iY2IfA547mKtqv2qJSJTvvhkaNbJaZpPvc6vqkmtouka34S5mVlcNyApQAEaY4ypKyqboe5b/EjuJsBEnKv3lc/D8R8hLNqv6mnJCTy3aBMH8/KJjQwNUJDGGGPqggqTu6qOq8E4jBe+dFj+FGycD8kX+1V1eEo8T3+6kU/XZXPJWW0DFKAxxpi6wNM9dxGJEZE2FexrIyJNqjcsU65O50N4bJV6zffu0JyW0WHMs/vuxhjT4HntUDcN+EsF+yZhz3OvGcGh0C3N6VRXVOhf1SBhaI945q/dS35hUYACNMYYUxd4Te6DgA8q2PchpWesM4Hky4DcHNj5hd9VhyUncOhoAZ9v2R+AwIwxxtQVXpN7LFDRIOmyc82bQOo6DIJCqvSM94HdWhEWEsRce8a7McY0aF6T+3pgVAX7MoCN1ROOOaXI5s699yrcd48OD2FAl5bMydyD+vkQGmOMMfWH1+T+OHCriDwsIj1FpIX7+nfgFuBfXg4iIh1EZL6IZIrItyJyu7v9LBFZKiLfiMh7ItK0ah+nkfBlQPZa2Of/b6q0lAS2789j/d4jAQjMGGNMXeB1+tlngfuAm4GvgWz39RbgT+5+LwqAO1U1GegP3CIiKTgd8iaqai9gBs7c9aYi3Uc6r1WZra6HM1vdnO+s17wxxjRUnqefVdUHgLY4zfPXuK9tVfVvfhxjt6p+6b4/DGQC7QAfsNAtNgd7PnzlWiRBfEqVmuZbx0ZwZvtYGxJnjDENmF9zy6vqQVWdraqvuK8Hq3piEUkE+gDLgTXAJe6uK4AOVT1uo9F9JGz9DHL97/melpzAqu0HyD58LACBGWOMqW1eJ7G5TERuKLGeJCKficgBEfmPiDTz56TupDf/Ae5Q1UPA9ThN9F8AMcDxCupNEJGVIrIyOzvbn1M2PL4M0ELYMM/vqsOS41GF+Wut17wxxjREXq/c/wSU7OT2ONAK+BtwNvBXrycUkVCcxP6Kqr4DoKprVXWEqvYFXqOC3veqOlVVU1U1NS4uzuspG6Z2fSE6rkpD4lLaNKVtbARzrGneGGMaJK/JvTPwDYCIxAIjgN+499v/CIz2chAREZzZ7jJV9R8ltse7r0E4PySe9voBGq2gIKdpfsNcKCi3oaNCIkJaSgKL1mdzNN+/me6MMcbUff7ccz8xMPpCoBCY667vALxeRg8ArgaGishqd8kAfi4i64C1wC7gBT/iarx8GXDsEGxd4nfVtOQEjuYXsWRDTgACM8YYU5sqe+RrSV8BV4nIMmA8MF9VT/TG6gh4unmrqosBqWC3p7HypoTOgyEkwuk132WIX1XP7dyCJuEhzM3MYlhyQkDCM8YYUzu8XrnfDVwKHMK5cv9ziX1jcXq8m5oWFgWdhzjJ3c8Z58JDgrmwexxzM/dSVGSz1RljTEPidRKbxThX6P2ATqpaMpk/j3Of3NQG30g4uA2yvvW7alpKPNmHj/HNziqPaDTGGFMH+TOJzWFV/UJVD5TZ/qGqrqv+0IwnxbPV+T+hzRBfPMFBwlzrNW+MMQ1KhffcReRm4C1VzXbfV0ZV9anqDc14EtPaGRb3/SwY5N+svc2iwujbqTlzvsvizhG+AAVojDGmplXWoe4JYCXOPPJPnOI4Clhyry2+dPjkATi8x0n2fhienMBfP8xk+/5cOrSIClCAxhhjalKFzfKqGqSqK0q8r2wJrrmQzUl8Gc5rFR4kk5bi9JS3ueaNMabh8GtueVNHxadAs45VepBMUqtousRFMzfTpqI1xpiGorJ77h39OZCqbjv9cEyViDhX719Mh+M/Qli0X9XTUhKYtmgzh47m0zQiNDAxGmOMqTGVXblv9rhscV9Nbeo+EgqOwqYFflcdnpxAQZGycF0jfxiPMcY0EJV1qBPgCDATeBfIrZGITNV0GgDhTZ0HyfQY5VfVPh2b0yI6jLnfZXHxmW0DFKAxxpiaUllyHwRcCVyOMwvde8DrwCxVza+B2Iw/QsKgaxqs+wiKipwHy3gUHCQM8cUz57s95BcWERpsXTGMMaY+q6y3/GJVvRVoC/wEyAOmA1ki8ryIjHCf4mbqCl8G/JgNO7/wu+rwlHgOHS1g5ZYfAhCYMcaYmnTK5KyqRao6R1VvABKAa4FI4APglQDHZ/zRLQ0kuErPeB/YLY6w4CCbrc4YYxoAf6+8z8Rprh+A89jX76s9IlN1kc2h0/lVGhIXHR7C+V1bMjczC/XzITTGGGPqllMmdxFJEZG/iMh6YCnQC7gXiFfVSQGOz/jLlwHZmbB/k99V05IT2Lovlw17jwQgMGOMMTWlwuQuIv8jIl8DX+NcrU8B2qrqSFWdrqqHaipI4wef+yCZ7/2frW5YcjyATWhjjDH1XGW95f8KHAbeAHYCnYG7RKS8sqqqf6j+8IzfWnSGuB7OfffzTvW8n9LaxEbSq10sczOzuGlwlwAFaIwxJtAqS+7bcB4Ic76H4yhgyb2u8KXDkscg7wfnPrwfhiXH869568k5coxWTcIDFKAxxphAqmwoXKKqJnlcOtdk0OYUfBmghbBhnt9V05ITUIVP1lrTvDHG1Fc2Tr0hatcXouOqNCSuZ9umtImNYO53NiTOGGPqK0vuDVFQMHS/CNbPhYLjflUVEdKSE1i0Poej+YUBCtAYY0wgWXJvqHwZcOwgbPvM76ppKQnk5Rfy2cacAARmjDEm0Cy5N1SdB0NIRJUmtOnfuQXRYcE2JM4YY+qpysa5dxQRe7h3fRUWDUkXOvfd/ZxxLjwkmAt9cczLzKKoyGarM8aY+uZUz3PvAyAin4hIj5oJyVQbXzoc2AZ7M/2umpacQNahY6zZdTAAgRljjAmkypJ7HhDlvh8MNA14NKZ6dT8xW53/veaH+OIJEqzXvDHG1EOVTWKzCviXiMxx138tIrsrKGsz1NVFTdtA27Od++6DfudX1ebRYaR2asGczL38doQvQAEaY4wJhMqS+y+Bh4ExODPQDQOOVVDWZqirq3wZMP8BOJwFMQl+VU1LiefBD9ey44dc2jePOnUFY4wxdUJlM9StVdXRqtoNEGCszVBXD/nSndd1/j9IJi3Z+TEwz3rNG2NMveJ1KFwSsDqQgZgASegJsR2rNCSuc1wTOsdFMzfT7rsbY0x94im5q+pWoEhEfiYij4vIK+7rT0WksqZ9U9tEnKv3TfPheK7f1YcnJ7Bs0z4OH80PQHDGGGMCwVNyF5F4YCXwGjAK5/Gvo4DXgc9FJC5gEZrT5xsJBUdh86d+V01LSSC/UFm4zmarM8aY+sJrs/w/gJbAuaraWVXPc++zn+tu/0egAjTVoNMFEBZTpSFxZ3dsTvOoUGuaN8aYesRrcs8A/qCqn5fc6K7/D85VvKmrQsKgWxp8PxuKivyqGhwkDOkRzydr91JQ6F9dY4wxtcNrcg8HDlew7zAQVj3hmIDxZcCPe2HXl35XHZ6cwMG8fFZu/SEAgRljjKluXpP7MuAPIhJdcqO7/gd3v6nLuqaBBFepaX5g9zjCgoNstjpjjKknvCb3O4GewHYReV1E/iUirwHbgRR3v6nLolpAp/OrNCSuSXgI53VpydzMLNTPh9AYY4ypeV6Hwq0GugFTgThgOBAPPA10U9WvAhahqT6+dNj7Hezf7HfVtJQEtuzLZWP2jwEIzBhjTHXy/Dx3Vc1R1YmqOkxVU9zXu1XVxkjVFyceJFOl2eriAazXvDHG1AOek7tpAFp2gVa+KjXNt4mNpGfbpnbf3Rhj6gFL7o2NLx22LoG8A35XTUtO4IttP7DvSEXPDzLGGFMX1GhyF5EOIjJfRDJF5FsRud3d3ltElonIahFZKSL9ajKuRsWXAUUFsGGu31WHpySgCp+stQfJGGNMXVbTV+4FwJ2qmgz0B24RkRTg78CfVbU3cK+7bgKhfSpEtapS03zPtk1p3TTC7rsbY0wdd8rkLiLhIvJHETnrdE+mqrtV9Uv3/WEgE2iH8zz4pm6xWGDX6Z7LVCAo2OlYt34OFPr3MBgRIS0lnoXrcjiaXxigAI0xxpyuUyZ3VT0G/BFoVp0nFpFEoA+wHLgDeFhEtgOP4ExpW16dCW6z/crs7OzqDKdx8aXDsYOw9TO/q6YlJ5CXX8jSTfsCEJgxxpjq4LVZfjnQt7pOKiJNgP8Ad6jqIeAm4Deq2gH4DTCtvHqqOlVVU1U1NS7OHkRXZV2GQHB4lZrmz+vSkuiwYOs1b4wxdZjX5H4XcJOI3CoinUUkWkSiSi5eTygioTiJ/RVVfcfdfC1w4v1bgHWoC6SwaOh8oTMVrZ8zzoWHBDOoe5zNVmeMMXWYP1fuXYDHgPXAIZwHxpRcTklEBOeqPFNVSz4mdhdwoft+qHsOE0i+dDiwFbLX+l11WHICWYeOsWbnoQAEZowx5nSFeCx3PU6nt9M1ALga+EZEVrvb7gZ+CfxLREKAo8CEajiXqUz3kcBvnKv3+GS/qg7xxREkMCczi17tYwMTnzHGmCrzlNxVdXp1nExVFwNSwe5qu6dvPGjaFtr2ce67D/TvuT8tm4TTt1Nz5n6XxW+Hdw9QgMYYY6rKr3HuIpIiIleLyN0i0trd1lVEYgITngkoXwbsWAmH/e8cl5acwHe7D7HzQF4AAjPGGHM6PCV3EWkiIm8Ca4DngPuBtu7uB4H7AhOeCShfOqCw/iO/q6alJADwiU1oY4wxdY7XK/d/AOcDw4AYSjetfwiMrOa4TE1IOANiO1RpSFyXuCZ0bhXNnEybitYYY+oar8n9J8AfVHU+UHZqsq1Ap2qNytQMEadj3cb5kO9/83paSgJLN+Zw+Kh/M90ZY4wJLK/JPRKoaEqyGE5O+Ka+8KVDQR5s+tTvqsN6xJNfqCxanxOAwIwxxlSV1+T+OXBNBfsuB/yfx9TUDYkXQFiMMyTOT307NadZVKjNVmeMMXWM13HufwLmishcnBnkFMgQ+f/t3XuYHXV9x/H3Zy+5XyAJWUIuJCTZQyIoSLiJAc3N7Gq9oyLaqrWhj1YB8bEVe1Frn/qoxdoKRW5Vn1KwSlQUDBAIhaAECCCYewgJBMImIYQkQG673/4xZ3Xd5nLO5syZ3dnP63nOc/bMzG/ON/Pk2c/OzNgaSF0AABMsSURBVG9+P11KEu7nplSfpa2uL0yaCasXQFsb1JT+AEVdbQ0zCiO5Z9Vm9re2UVdb7UkGzczsQEr6bVx8Pn0m0Bf4LkmHuq8AJwCzIuLh1Cq09BWaYVcLPP9Y2U1nTW1g+6v7WLrhpRQKMzOzrij5VCsiHoiI6SRTs44BBkfEORHxQGrVWXVMng2q7dKl+XMbj6FPbY3neDcz60a6ch11N7AP8OgleTFgGIw7u0uPxA3qW8dZE4dztx+JMzPrNkoOd0nNkn5NEu4vALsl/VrS21Orzqqn0ASbl8FLG8puOnvKSNZtfYWntuxKoTAzMytXqSPUXQT8AtgFXAycX3zfBdxaXG89WaEpeV+9oOymM6Yko9W517yZWfdQ6pn75cA1ETEnIq6OiPnF9znAtcCX0ivRqmL4RBjR2KX77qOP6s/UUUN8393MrJsoNdyHA/MPsu4WYFhlyrFMFZpg/WLY/XLZTWdNbWDphpd4cdeeFAozM7NylBrui4DzDrLuPOC+ypRjmSo0Q9t+WLuw7KazpzTQFrBo1ZYUCjMzs3IcdBAbSVM7fPw34DpJw4GfAZuBkcB7gCbgk2kWaVUy5nQYMDzpNX/S+8pqetLoITQM6cvC5S28/7QxKRVoZmalONQIdb8jGYmunYCLiq/gj2eGWwDUVrw6q66a2mQimZW/hNZ9UFtfclNJzJrSwE8fe47d+1rpV+//DmZmWTlUuL+1alVY91FogsdvhGd+AxPKG1V41tQGblzyDA+ue5G3FEamVKCZmR3OQcM9IsqfJsx6vhPeCrV9YNWCssP97BOGM6BPLQtXtDjczcwyVPYIdZLqJA3o/EqjOMtA30Ew4bzkkbiIw2/fQb/6WqZPHsHC5ZuJMtuamVnllDqIzVBJV0naRDJC3c4DvCwvCk3w0tOwZVXZTWdNaeCFHbtZ9vyOFAozM7NSlDrl6/dJHnm7FlgL7E2rIOsGGufCbZ9Lzt5HnlhW0xknjkSCu5a3cNLooSkVaGZmh1JquM8ELoqIm9IsxrqJoaNh1CnJI3HTP1dW0+GD+nLauKNZuKKFS2c3plSgmZkdSqn33J8BXk2zEOtmCs2w8WHYVf5sb7OmNrDs+R1setkTB5qZZaHUcP8C8LeSxqVZjHUjhSYgYPUdZTed1T6RjKeBNTPLREnhHhG3A/cDayWtlvRQ51e6ZVrVHXsyDBnTpVniJh4zkAkjBnqWODOzjJR0z13St4BLgIdxh7reQYLCXHj8v2Hfa1Dfv4ymYuaJI/nhbzawa89+BvUttWuHmZlVQqmX5T8JfCkizoyICyPi451faRZpGSk0wb5X4eny5wWaNbWBva1t3L/aE8mYmVVbqeH+KrA0zUKsGxo/HfoM6tIc79OOP5qh/eu5y3O8m5lVXanh/h1gniQddkvLj7q+MGlmMhRtW1t5TWtrmHHiSBat3Exrm0erMzOrplJvho4AzgRWSboX2N5pfUTEX1eyMOsmCs2w/Oew6TEYfVpZTdtniXv0mZc4ffywlAo0M7POSg339wP7gXpg9gHWB+Bwz6PJc0A1yYA2ZYb7uY0jqK8VC5e3ONzNzKqo1EfhJhzmdULahVpGBgyDcWcnl+bLNLhfPWedMNz33c3MqqzsWeGsF2qcCy1PwvZnym46e2oD67a8wlNbdqVQmJmZHUipz7l/6nDbRMRVR16OdUuFZrjr75Kz9zPnldV0xokj+fufL+PuFS1MPGZQSgWamVlHpd5z/+4h1rV3hXa459WISTB8cvJIXJnhPuboAUwZNYSFyzcz79yJKRVoZmYdlXrPvabzCxgGXAD8FpiaZpHWDRSaYP1i2P1y2U1nTxnJIxu2se0VD2xoZlYNXb7nHhHbI+JHwNXA9ypXknVLhWZo2wdr7y676aypDbQFLFrpiWTMzKqhEh3qngamVWA/1p2NPQP6D0seiSvTSccNpWFIX+5e6V7zZmbVcEThLmkUcBlJwFue1dQmvebX3Amt+8trWiNmTmngf1dtYc/+1pQKNDOzdiWFu6QtkjZ3em0HNgLTgc+XuJ+xkhZJWiFpmaSLi8t/JOnx4mu9pMe7/C+y9BSaYPd2ePbBspvOntLAK3tbeXDdthQKMzOzjkrtLX8lf+gV3243SbgviIgXS9zPfuCyiHhU0mBgqaS7IuKD7RtI+heg/F5blr6JM6C2T3Jpfvyby2p69sTh9K+vZeHyFs5rPCalAs3MDEoM94j4ciW+LCI2AZuKP++UtAIYDSwHKE5M8wFgRiW+zyqs7yCYcC6svA3mfC2Z871E/eprmT55BAtXtPDB08emWKSZWffQr76GSSMHZ/LdpZ65V5yk8cCpwJIOi6cDLRGxJouarASFJrjtMti6Go4plNV0zuuO5c7lLbzj3xenVJyZWfcxddQQbr94eibffdBwl3RPGfuJiJhZ6saSBgG3AJdExI4Oqy4AbjpEu3nAPIBx48aVUZ5VTGMx3FfdXna4v/uU4zhmcF/27i9v+lgzs55oUN/Mzp9RxIHn2pb04xLajwLeRBLutSV9oVQP/BK4IyKu6LC8DngOOC0iNh5uP9OmTYtHHnmklK+0SvveuVDXD/78zqwrMTPrNSQtjYiSHj0/6J8VEXH+Ib5gHMkUr+8AtgLfLrEwAdcDKzoGe9EsYGUpwW4ZKzTDvV+HXVtgkDvHmZl1N2U95y5pkqTrgTXAO4EvAsdHxD+XuItzgI8CMzo8+tZcXPchDnFJ3rqRQhMQyTPvZmbW7ZQ6K9zrgC8B5wPPAhcDN0REWYOFR8Ri4IBdrCPiY+XsyzJ07OthyOjkvvupF2ZdjZmZdXLIM3dJp0maDzxB0rP9k8DkiLi63GC3HJGS0eqeugf27c66GjMz6+Sg4S7pV8BDwATgQxExJSJ+EBEeP9SS++77XoWn78u6EjMz6+RQZ+5vI7mEPha48gDDz/7RqzrlWrcxYTr0GZRcmjczs27lUPfcv1K1KqznqeubDEe7egG0tUFNJSYYNDOzSjjUo3AOdzu0QjOsuBU2PQ6j35h1NWZmVuTTLeu6yXNANcnZu5mZdRsOd+u6gcNh7Fm+725m1s043O3IFObCC0/C9mezrsTMzIoc7nZkCsUBBn1p3sys23C425EZMRmGT/KleTOzbsThbkeu0ARP3w+7dxx+WzMzS53D3Y5coRna9sFTd2ddiZmZ4XC3ShhzBvQfBqt8393MrDtwuNuRq62DxrfBmjugdX/W1ZiZ9XoOd6uMxrnw2kvw7JKsKzEz6/Uc7lYZk2ZCbR/3mjcz6wYc7lYZfQfD+OlJuEdkXY2ZWa/mcLfKKTTBtnWwdU3WlZiZ9WoOd6ucQlPy7kvzZmaZcrhb5QwdA8e+3kPRmpllzOFulVVoTnrMv7I160rMzHoth7tVVqEJog3W3Jl1JWZmvZbD3Spr1Btg8HG+725mliGHu1WWlMzxvvYe2Lc762rMzHolh7tVXqEZ9r0C6+/PuhIzs17J4W6VN3461A/0pXkzs4w43K3y6vvBpBnJLHEerc7MrOoc7paOQjPsfB42/TbrSszMeh2Hu6Vj8hxQDaz6VdaVmJn1Og53S8fAETDmDN93NzPLgMPd0lNogheegJc3Zl2JmVmv4nC39BSak3dfmjczqyqHu6VnxGQYNtHhbmZWZQ53S4+UXJpffz/s2Zl1NWZmvYbD3dJVaIbWvfDUPVlXYmbWazjcLV1jz4T+R/vSvJlZFTncLV21dTD5bbD6Dmjdn3U1Zma9gsPd0leYC69tg40PZV2JmVmv4HC39E2cCTX1HtDGzKxKHO6Wvn5DYMJ033c3M6sSh7tVR6EZXlwLW9dkXYmZWe5VNdwljZW0SNIKScskXdxh3WckrSou/0Y167IqaJybvPvs3cwsdXVV/r79wGUR8aikwcBSSXcBDcC7gNdHxB5JI6tcl6XtqLFw7MlJuJ/z2ayrMTPLtaqGe0RsAjYVf94paQUwGvgL4OsRsae4bnM167IqKTTDfd+E62ZlXYmZWfqGT4L3XJ3JV1f7zP33JI0HTgWWAN8Epkv6J2A38PmIePgAbeYB8wDGjRtXtVqtQk65EDY9Aa17sq7EzCx99QMy++pMwl3SIOAW4JKI2CGpDjgaOAs4HfgfSSdERHRsFxHXANcATJs2LbCe5ejj4cM3Z12FmVnuVb23vKR6kmC/MSLmFxdvBOZH4iGgDRhR7drMzMzyoNq95QVcD6yIiCs6rPoZMKO4TSPQB9hazdrMzMzyotqX5c8BPgo8Kenx4rLLgRuAGyT9DtgL/FnnS/JmZmZWmmr3ll8M6CCrP1LNWszMzPLKI9SZmZnljMPdzMwsZxzuZmZmOeNwNzMzyxmHu5mZWc443M3MzHLG4W5mZpYzDnczM7OccbibmZnljHrqKK+StgAbKrjLEXg8+2rwca4eH+vq8HGuDh9nOD4ijillwx4b7pUm6ZGImJZ1HXnn41w9PtbV4eNcHT7O5fFleTMzs5xxuJuZmeWMw/0Prsm6gF7Cx7l6fKyrw8e5Onycy+B77mZmZjnjM3czM7OccbgDkuZKWiVpraS/ybqePJI0VtIiSSskLZN0cdY15ZmkWkmPSfpl1rXklaSjJP1E0sri/+uzs64pjyRdWvyd8TtJN0nql3VNPUGvD3dJtcCVQBMwFbhA0tRsq8ql/cBlETEFOAv4tI9zqi4GVmRdRM59B1gQEScCb8DHu+IkjQY+C0yLiJOAWuBD2VbVM/T6cAfOANZGxLqI2AvcDLwr45pyJyI2RcSjxZ93kvwiHJ1tVfkkaQzwduC6rGvJK0lDgHOB6wEiYm9EbM+2qtyqA/pLqgMGAM9nXE+P4HBPAubZDp834tBJlaTxwKnAkmwrya1/Bb4AtGVdSI6dAGwB/rN4++M6SQOzLipvIuI54FvAM8Am4OWIuDPbqnoGhzvoAMv8CEFKJA0CbgEuiYgdWdeTN5LeAWyOiKVZ15JzdcAbgf+IiFOBVwD316kwSUeTXEmdABwHDJT0kWyr6hkc7smZ+tgOn8fgyz6pkFRPEuw3RsT8rOvJqXOAd0paT3KLaYak/8q2pFzaCGyMiParTz8hCXurrFnA0xGxJSL2AfOBN2VcU4/gcIeHgcmSJkjqQ9JZ49aMa8odSSK5P7kiIq7Iup68iogvRsSYiBhP8n/5nojwmU6FRcQLwLOSCsVFM4HlGZaUV88AZ0kaUPwdMhN3XCxJXdYFZC0i9kv6K+AOkp6YN0TEsozLyqNzgI8CT0p6vLjs8oi4PcOazI7EZ4AbiycF64CPZ1xP7kTEEkk/AR4leeLmMTxSXUk8Qp2ZmVnO+LK8mZlZzjjczczMcsbhbmZmljMOdzMzs5xxuJuZmeWMw92sh5L0ZUlbiz83Fj8flUEdH5D0sQMsv7f4GJOZVZnD3SwfGoF/AKoe7sAHgI8dYPmngC9WtxQzAw9iY2YHIKl/RLx2JPuICI/YZpYRn7mb9XCS3gL8ovjxaUlRHFu+ff04STdL2ibpVUl3dBg2FUnji20ulPRDSdvb9yfpTyUtLrZ9SdIiSdM6tP0+8D7gvOI+QtKXi+v+32V5STMkLZG0W1KLpKuKkwn9/t9S3MdbJP1Y0i5J6yR9qsKHzSzXfOZu1vM9CnyeZGrM95JMjbkHQNIwYDHwIvCXwKsks5ctlNTY6ez8WyQTc5wPtBaXjQd+CDwF9AE+DNwn6aSIWAf8IzCO5HZAewBvPFCRkqYCC4C7SP4gGAt8nWT61LmdNr8W+AHJUKMXAFdKeiQiHirjuJj1Wg53sx4uInZIWlX8+FhErO+w+lJgIHBKRGwDkPQAsB74BHBlh20fjIhPd9r3V9t/llRDEsynAx8BvhoRT0naBtRExIOHKfXvgQ3AOyOitbjPbcCPJJ0dEb/psO1NEfG14jb3An9C8oeLw92sBL4sb5Zvs0gCeYekOkl1wE5gKTCt07a3dW4saYqkn0pqITmb3wcUSDrwlesM4KftwV50C8mEIG/utO2d7T8Up/pcQzIds5mVwGfuZvk2AjgL+OAB1t3d6XNLxw+SBpOEbAvwOZKz7t3AdUC/LtQyqvN3RESrpBeBYZ223d7p894ufqdZr+RwN8u3bcCtJPfGO9vZ6XPnKSLPJjlbnh0RK9sXShraxVo2ASM7LpBUCwwv1mlmFeJwN8uHvcX3zme3d5M8h76sC4+29S++72lfIOlNJJ3slnb67lLOqpcA75F0eYdL8+8l+T20uMzazOwQfM/dLB/aO9RdJOlMSScXP19B0sv9HkkflnRecUS5KyVdcJh9PgjsAq6VNEfSJ4Cbgec6bbcSOFnSuyVNk3TcQfb3NZI/DH4mqVnSPJLe8Hd06kxnZkfI4W6WAxGxgeRxuPcCD1B8Tj0itpLcc18JfJvkHvo3gKHAE4fZZwvJY3HHAj8HLiF5nG5tp02vKu73BuBhYN5B9rcMaCK5ND+fJOxvAt5fzr/VzA5PEZ1vs5mZmVlP5jN3MzOznHG4m5mZ5YzD3czMLGcc7mZmZjnjcDczM8sZh7uZmVnOONzNzMxyxuFuZmaWMw53MzOznPk/lKaP3msHFMwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x524d47f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This file is associated with the book\n",
    "# \"Machine Learning Refined\", Cambridge University Press, 2016.\n",
    "# by Jeremy Watt, Reza Borhani, and Aggelos Katsaggelos.\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import csv\n",
    "\n",
    "# sigmoid for softmax/logistic regression minimization\n",
    "def sigmoid(z): \n",
    "    y = 1/(1+np.exp(-z))\n",
    "    return y\n",
    "    \n",
    "# import training data \n",
    "def load_data(csvname):\n",
    "    # load in data\n",
    "    reader = csv.reader(open(csvname, \"r\"), delimiter=\",\")\n",
    "    d = list(reader)\n",
    "\n",
    "    # import data and reshape appropriately\n",
    "    data = np.array(d).astype(\"float\")\n",
    "    X = data[:,0:8]\n",
    "    y = data[:,8]\n",
    "    y.shape = (len(y),1)\n",
    "    \n",
    "    # pad data with ones for more compact gradient computation\n",
    "    o = np.ones((np.shape(X)[0],1))\n",
    "    X = np.concatenate((o,X),axis = 1)\n",
    "    X = X.T\n",
    "    return X,y\n",
    "\n",
    "# YOUR CODE GOES HERE - create a gradient descent function for softmax \n",
    "# cost/logistic regression \n",
    "def softmax_grad(X,y, w0):\n",
    "    w1 = w0\n",
    "    alpha = alpha0\n",
    "    w_path = [] \n",
    "    max_its = 10\n",
    "    mis = []\n",
    "    steps = []\n",
    "    for k in range(max_its):  \n",
    "        t_temp = y*np.dot(X.T,w1)\n",
    "        grad =  -np.dot(X,sigmoid(-t_temp)*y)\n",
    "        \n",
    "        r_temp = sigmoid(-t_temp)*(1-sigmoid(-t_temp))\n",
    "        r_temp = r_temp.T\n",
    "        r_ = np.diag(r_temp[0])\n",
    "        hess = np.dot(np.dot(X,r_),X.T)\n",
    "        \n",
    "        w1 = w1 - np.dot(np.linalg.pinv(hess),grad)\n",
    "        \n",
    "        pred_label = np.dot(X.T[:,:], w1)\n",
    "        pred_label[pred_label>0] = 1\n",
    "        pred_label[pred_label<=0]  = -1\n",
    "        mis.append(sum(pred_label!= y))\n",
    "        \n",
    "    return mis\n",
    "\n",
    "\n",
    "def square_margin(X,y,w0):\n",
    "    w1 = w0\n",
    "    alpha = alpha0\n",
    "    w_path = [] \n",
    "    max_its = 10\n",
    "    mis = []\n",
    "    steps = []\n",
    "    for k in range(max_its):  \n",
    "        t_temp = 1 - y*np.dot(X.T,w1)\n",
    "        t_temp[t_temp<0]=0\n",
    "        grad = np.dot(X, -2 * t_temp *y)\n",
    "        thre =  1-y*np.dot(X.T, w1)\n",
    "        keep = []\n",
    "        for x in range(thre.shape[0]):\n",
    "            if (thre[x]>0): keep.append(x)\n",
    "                \n",
    "        filtered_X = X[:, keep]\n",
    "        hess = 2 * np.dot(filtered_X, filtered_X.T)\n",
    "        w1 = w1 - np.dot(np.linalg.pinv(hess),grad)\n",
    "\n",
    "        pred_label = np.dot(X.T[:,:], w1)\n",
    "        pred_label[pred_label>0] = 1\n",
    "        pred_label[pred_label<=0]  = -1\n",
    "        mis.append(sum(pred_label!= y))\n",
    "        \n",
    "    return mis\n",
    "\n",
    "\n",
    "# load in data\n",
    "X,y = load_data('/Users/Administrator/Desktop/python/breast_cancer_data.csv')\n",
    "\n",
    "fig, ax= plt.subplots(1,1,figsize=(8,6))\n",
    "\n",
    "ini = 0.000001\n",
    "w0 = np.array([ini,ini,ini,ini,ini,ini,ini,ini,ini])\n",
    "w0.shape = (9,1)\n",
    "alpha0 = 750\n",
    "\n",
    "# run gradient descent\n",
    "mis_square = square_margin(X,y,w0)\n",
    "#print (mis_square)\n",
    "mis_softmax = softmax_grad(X,y,w0)\n",
    "# print (mis_softmax)\n",
    "\n",
    "ax.plot(range(10), mis_softmax, label='Softmax')\n",
    "ax.plot(range(10), mis_square, label='Square_margin')\n",
    "ax.set_xlabel('Iteration',fontsize=15)\n",
    "ax.set_ylabel('Number of Misclassifications',fontsize=15)\n",
    "ax.set_title('Breast Cancer Dataset',fontsize=15)\n",
    "\n",
    "ax.legend(fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 4_16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training for class 0\n",
      "Iter 1 \n",
      "Iter 1001 \n",
      "Iter 2001 \n",
      "Iter 3001 \n",
      "Iter 4001 \n",
      "Training for class 1\n",
      "Iter 1 \n",
      "Iter 1001 \n",
      "Iter 2001 \n",
      "Iter 3001 \n",
      "Iter 4001 \n",
      "Training for class 2\n",
      "Iter 1 \n",
      "Iter 1001 \n",
      "Iter 2001 \n",
      "Iter 3001 \n",
      "Iter 4001 \n",
      "Training for class 3\n",
      "Iter 1 \n",
      "Iter 1001 \n",
      "Iter 2001 \n",
      "Iter 3001 \n",
      "Iter 4001 \n",
      "Training for class 4\n",
      "Iter 1 \n",
      "Iter 1001 \n",
      "Iter 2001 \n",
      "Iter 3001 \n",
      "Iter 4001 \n",
      "Training for class 5\n",
      "Iter 1 \n",
      "Iter 1001 \n",
      "Iter 2001 \n",
      "Iter 3001 \n",
      "Iter 4001 \n",
      "Training for class 6\n",
      "Iter 1 \n",
      "Iter 1001 \n",
      "Iter 2001 \n",
      "Iter 3001 \n",
      "Iter 4001 \n",
      "Training for class 7\n",
      "Iter 1 \n",
      "Iter 1001 \n",
      "Iter 2001 \n",
      "Iter 3001 \n",
      "Iter 4001 \n",
      "Training for class 8\n",
      "Iter 1 \n",
      "Iter 1001 \n",
      "Iter 2001 \n",
      "Iter 3001 \n",
      "Iter 4001 \n",
      "Training for class 9\n",
      "Iter 1 \n",
      "Iter 1001 \n",
      "Iter 2001 \n",
      "Iter 3001 \n",
      "Iter 4001 \n",
      "Accuracy for training set is: 0.839917\n",
      "Accuracy for test set is: 0.836300\n"
     ]
    }
   ],
   "source": [
    "import csv\n",
    "import numpy as np\n",
    "\n",
    "def sigmoid(t):\n",
    "    return 1/(1 + np.exp(-t))\n",
    "\n",
    "def checkSize(w, X, y):\n",
    "    # w and y are column vector, shape [N, 1] not [N,]\n",
    "    # X is a matrix where rows are data sample\n",
    "    assert X.shape[0] == y.shape[0]\n",
    "    assert X.shape[1] == w.shape[0]\n",
    "    assert len(y.shape) == 2\n",
    "    assert len(w.shape) == 2\n",
    "    assert w.shape[1] == 1\n",
    "    assert y.shape[1] == 1\n",
    "\n",
    "def compactNotation(X):\n",
    "    return np.hstack([np.ones([X.shape[0], 1]), X])\n",
    "\n",
    "def readData(path):\n",
    "    \"\"\"\n",
    "    Read data from path (either path to MNIST train or test)\n",
    "    return X in compact notation (has one appended)\n",
    "    return Y in with shape [10000,1] and starts from 0 instead of 1\n",
    "    \"\"\"\n",
    "    reader = csv.reader(open(path, \"r\"), delimiter=\",\")\n",
    "    d = list(reader)\n",
    "    data = np.array(d).astype(\"float\")\n",
    "    X = data[:,0:-1]\n",
    "    y = data[:,-1]\n",
    "    y.shape = (len(y),1)\n",
    "\n",
    "    # pad data with ones for more compact gradient computation\n",
    "    o = np.ones((np.shape(X)[0],1))\n",
    "    X = np.concatenate((o,X),axis = 1)\n",
    "#     X = X.T\n",
    "    \n",
    "    return X,y\n",
    "\n",
    "def softmaxGrad(w, X, y):\n",
    "    checkSize(w, X, y)\n",
    "    ### RETURN GRADIENT\n",
    "    t1 = y*np.dot(X,w)\n",
    "    grad = -np.dot(X.T,sigmoid(-t1)*y)\n",
    "    return grad\n",
    "\n",
    "def accuracy(OVA, X, y):\n",
    "    \"\"\"\n",
    "    Calculate accuracy using matrix operations!\n",
    "    \"\"\"\n",
    "    \n",
    "    pred_vector = np.dot(X[:,:], OVA)\n",
    "    pred_label  = np.argmax(pred_vector, axis =1)\n",
    "    pred_label = pred_label.reshape(X.shape[0],1)\n",
    "    acc_ = sum(pred_label==y)\n",
    "    acc = float(acc_) / X.shape[0] \n",
    "    \n",
    "    return acc\n",
    "    \n",
    "\n",
    "def gradientDescent(grad, w0, *args, **kwargs):\n",
    "    max_iter = 5000\n",
    "    alpha = 0.001\n",
    "    eps = 10^(-5)\n",
    "\n",
    "    w = w0\n",
    "    iter = 0\n",
    "    while True:\n",
    "        gradient = grad(w, *args)\n",
    "        w = w - alpha * gradient\n",
    "\n",
    "        if iter > max_iter or np.linalg.norm(gradient) < eps:\n",
    "            break\n",
    "\n",
    "        if iter  % 1000 == 1:\n",
    "            print(\"Iter %d \" % iter)\n",
    "        iter += 1\n",
    "\n",
    "    return w\n",
    "\n",
    "def oneVersusAll(Y, value):\n",
    "    \"\"\"\n",
    "    generate label Yout, \n",
    "    where Y == value then Yout would be 1\n",
    "    otherwise Yout would be -1\n",
    "    \"\"\"\n",
    "    y = np.copy(Y)\n",
    "    yes = y==value\n",
    "    no  = y!=value\n",
    "    y[yes] = 1\n",
    "    y[no] = -1\n",
    "    return y\n",
    "\n",
    "\n",
    "trainX, trainY = readData('/Users/Administrator/Desktop/python/MNIST_train_data.csv')\n",
    "\n",
    "#trainX = trainX[:50]\n",
    "#trainY = trainY[:50]\n",
    "\n",
    "# # training individual classifier\n",
    "Nfeature = trainX.shape[1]\n",
    "Nclass = 10\n",
    "OVA = np.zeros((Nfeature, Nclass))\n",
    "for i in range(Nclass):\n",
    "    print(\"Training for class \" + str(i))\n",
    "    w0 = np.random.rand(Nfeature, 1)\n",
    "    OVA[:, i:i+1] = gradientDescent(softmaxGrad, w0, trainX, oneVersusAll(trainY, i))\n",
    "\n",
    "\n",
    "print(\"Accuracy for training set is: %f\" % accuracy(OVA, trainX, trainY))\n",
    "\n",
    "testX, testY = readData('/Users/Administrator/Desktop/python/MNIST_test_data.csv')\n",
    "print(\"Accuracy for test set is: %f\" % accuracy(OVA, testX, testY))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}