{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Discriminative Classification" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Preliminaries\n", "\n", "- Goal \n", " - Introduction to discriminative classification models\n", "- Materials \n", " - Mandatory\n", " - These lecture notes\n", " - Optional\n", " - Bishop pp. 213 - 217 (Laplace approximation) \n", " - Bishop pp. 217 - 220 (Bayesian logistic regression) \n", " - [T. Minka (2005), Discriminative models, not discriminative training](https://github.com/bertdv/BMLIP/blob/master/lessons/notebooks/files/Minka-2005-Discriminative-models-not-discriminative-training.pdf)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Challenge: difficult class-conditional data distributions\n", "\n", "Our task will be the same as in the preceding class on (generative) classification. But this time, the class-conditional data distributions look very non-Gaussian, yet the linear discriminative boundary looks easy enough:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAG2CAYAAACQ++e6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOUUlEQVR4nO3de5QU5Z3/8c/QOFwijIuIDHI1awR1BQET/Cnh4i7GeOeoMRE18coPNBD25yoGHBgQNBiDxqCgLjHe4jkGjG6ikRMvkGU9QUajUdFVUVgYlyGaQVGHMFO/Pyo1U11T1+7qS3W9X+f06enq6qqna0bq4/d56qkqwzAMAQAApESXUjcAAACgmAg/AAAgVQg/AAAgVQg/AAAgVQg/AAAgVQg/AAAgVQg/AAAgVQg/AAAgVQg/AAAgVQg/AAAgVRIVfvbv36958+Zp2LBh6tGjhw4//HDV19erra2t1E0DAAAJ0bXUDYjilltu0d133637779fRx99tF566SV973vfU01NjWbNmlXq5gEAgARIVPj5r//6L5111lk67bTTJElDhw7VI488opdeeqnELQMAAEmRqPBz0kkn6e6779bbb7+tr3zlK/rTn/6kP/zhD1q+fLnr+i0tLWppaWl/3dbWpo8++kgHH3ywqqqqitRqAACQD8Mw9Mknn2jAgAHq0iWGETtGgrS1tRnXX3+9UVVVZXTt2tWoqqoylixZ4rl+XV2dIYkHDx48ePDgUQGP7du3x5InqgzDMJQQv/zlL3Xttddq2bJlOvroo/XKK69o9uzZuu2223TJJZd0Wt9Z+WlubtbgwYO1fft29e7du5hNBwAAOdqzZ48GDRqkv/71r6qpqcl7e4kKP4MGDdL111+vmTNnti9bvHixHnzwQW3ZsiXw83v27FFNTY2am5sJPwAAJETc5+9EXer+2Wefderry2QyXOoOAABCS9SA5zPOOEM33XSTBg8erKOPPlovv/yybrvtNl166aWlbhoAAEiIRHV7ffLJJ5o/f77Wrl2rXbt2acCAAfr2t7+tG2+8UdXV1YGfp9sLAIDkifv8najwk68wB88wDO3fv1+tra1Fbl3lyWQy6tq1K9MKAADyEnf4SVS3V6Ht27dPjY2N+uyzz0rdlIrRs2dP1dbWhqrMAQBQDISfv2tra9PWrVuVyWQ0YMAAVVdXU7HIg2EY2rdvn5qamrR161YdccQR8UxMBQBAngg/f7dv3z61tbVp0KBB6tmzZ6mbUxF69OihAw44QB988IH27dun7t27l7pJAAAk61L3YqA6ES+OJwCg3HBmAgAAqUL4AQAAqUL4icGCBdKiReHWXbTIXB8AAJQG4ScGmYx0443BAWjRInO9TKY47QIAAJ0RfmIwf75UX+8fgKzgU19vrp9UK1as0LBhw9S9e3eNGTNGGzZsKHWTAACIhEvdY2IFmhtvzH4tVU7wefTRRzV79mytWLFCJ554olauXKlTTz1Vb7zxhgYPHlzq5gEAEAqVnxi5VYCKFXwGDhyoFStWZC3buHGjevbsqQ8++CCWfdx222267LLLdPnll2vEiBFavny5Bg0apLvuuiuW7QMAUAyEn5jZA1C3bsWr+IwbN06bNm1qf20YhmbPnq3Zs2dryJAhWesuWbJEBx54oO/D2Z21b98+bd68WVOmTMlaPmXKFG3cuLFwXwwAgJjR7VUA8+dLixdL+/ZJ1dXF6eoaN26cfv7zn7e/fuCBB7Rt2zbNnTu307rTp0/X+eef77u9ww47LOv17t271draqkMPPTRr+aGHHqoPP/ww94YDAFBkhJ8CWLSoI/js22e+Lkbl57rrrtOnn36qLl266IYbbtDixYvVq1evTuv26dNHffr0yWk/zvudGYbBPdAAAIlCt1fM7GN8WlqCrwKLy9ixY5XJZNTQ0KCbb75ZBx98sC699FLXdXPp9urbt68ymUynKs+uXbs6VYMAAChnVH5i5Da42e8qsDh1795dI0eO1Jo1a7Rq1So9+eSTnvfVyqXbq7q6WmPGjNG6det0zjnntC9ft26dzjrrrPy/AAAARUL4iYnfVV3FCkDjxo3THXfcodNPP10nn3yy53q5dnvNmTNHF110kcaOHasTTjhBq1at0rZt2zR9+vR8mg0AQFERfmIQ5nL2YgSgUaNGqWvXrlq2bFn8G5f0rW99S3/5y19UX1+vxsZGHXPMMfrtb3/b6WoyAADKGeEnBq2t4S5nt95vbS1MOx566CHNmDFDRx55ZGF2IGnGjBmaMWNGwbYPAEChEX5iEOVGpXFXfNra2tTU1KT77rtPb731ltauXRvvDgAAqDCEn4Rbv369Jk+erOHDh2vNmjWqqakpdZMAAChrhJ+Emzhxotra2krdDAAAEoN5fgAAQKoQfgAAQKoQfgAAQKoQfgAAQKoQfgAAQKoQfgAAQKoQfgAAQKoQfgAAQKoQfgqlsdG870VjY6lbAgAAbAg/hdLYKC1cSPgBAKDMEH4Q2vr163XGGWdowIABqqqq0uOPP17qJgEAEBnhB6Ht3btXI0eO1J133lnqpgAAkDPCT4UYOHCgVqxYkbVs48aN6tmzpz744INY9nHqqadq8eLFmjp1aizbA7IUcpxcpYzBq5TvAZQY4ScOjY1SQ0Pnh+S+vAD/cI0bN06bNm1qf20YhmbPnq3Zs2dryJAhWesuWbJEBx54oO9jw4YNsbcR8FXIcXKVMgavUr4HUGJdS92AirBypfkPkpsrrui8rK7O/L+3GI0bN04///nP218/8MAD2rZtm+bOndtp3enTp+v888/33d5hhx0Wa/uAgmpsNP87vOoqqba21K0BUOYIP3G46irpzDOzlzU0mMHnnnuk0aOz3yvAP87jxo3Tddddp08//VRdunTRDTfcoMWLF6tXr16d1u3Tp4/69Onjua2dO6WPP5Z69Aje786d5vOAAbm2HIiBVRE588zg/76SGpSsdiOcpP6eURR0e8WhttYMOM6H5L68AP8hjh07VplMRg0NDbr55pt18MEH69JLL3VdN6jb64gjDtRvfrOhPdh4+egjBa4DlJ1y7DoKM5ansVFatapoTfJtRxLGHZXj77nYkvK7KgEqPxWie/fuGjlypNasWaNVq1bpySefVJcu7tk2TLdXVdVhgVWdjz4y36Pqg8gaGzv/g2wfJ2fX1CQ984z0//5fuP9xCNq2VdJ8803zuba29JUBt8qV83tY7ZU6HyOpeN8jSpXNbxvlVpUpxzblK47fVYUi/FSQcePG6Y477tDpp5+uk08+2XO9oG4vy86d2d1an376qd555x3t2mUu+/TTrdq16xXt399HgwcPjuMrIC2ijpOTpClTwp2cwm572jTzuQBj8GJRBmMJC6YcT8rl0qZKDGFliPBTQUaNGqWuXbtq2bJlsWzPquhYAejtt1/SpEmT2t9fuHCOFi6ULrnkkqzB1kAgv3Fyy5ZJ//3f0tSp0iGHmBWPadOk3bvNk9MJJ/hXdk44QXrwwY73+vaVtm/vGIPXo4e5vQcflEaM8D7BlPok5DxG1nGQpFmzpM8/7zhGEifKfDU1ZT+XSrmEsApH+CmU2lrz/8SK+Mf70EMPacaMGTryyCNj26Y9APXqNVGbNhl0dcFb2MDg10VTWytde625DefFApK0Zo332Be3isiVV0pf/3rHZ/fvN3/+/POONlthyt6uXE9CcYQmt647q72StG+feQy+/nUz/OTS5VXqcBenKN2okvvx2r07+xkVjfBTKLW1RSlBt7W1qampSffdd5/eeustrV27NvZ9DBhg/rtiGFJVFcEniRYskDIZaf784HUXLZJaW3P88y3G/7VOnWqesO28rq5cudIMCVZYeuqpjvfcgtKcOVKvXtLZZ+d+ZVXQMQhzorba7eWuu8znfLruwvyu4ggVxVDJXYRhJeV3VSYIPwm3fv16TZ48WcOHD9eaNWtUU1MT+z527uwIPoZhviYAJUsmI914o/mzXwBatMhcr76+OO1ytXWr+WwN8LX+4d6yxXzevTu7q8f+D7j9SkvJPMFddVV2l5HFbRqKpibpG9+QjjiiI3w0NXVsJ46TRdQT9ZVXmvu2Ap4k/d//awYgZ9ddlGpOmG6eOEJFMU7KUacb6dKl876tv68tWzq/l4SgQACMhPCTcBMnTlRbW1vBtm8Nera6upyDoJEMVuDxC0D24BOmQlQw1s6dYWXx4s7LrWDgdQL3O2k5g5LkfjK2xhrFVc1ynqibmqSf/Ux68snsE3VTk7nvvn3N1/aJt/bty96mFS6iVN7CdPMEhYow46eKcVKO+ntesMC7TYsXd/ytOdsUd1ehXzBcsECaObMj6Ev+37MM5ptLEsIPPDmDj9R5EDQBKDn8AlDJgo/9H/+mJunCC82T/e23h/u81aV15ZX+69XWmut4dSVZ7bCqTVb1SeqoCNgvNc+nEuD8bEODGXyk7BO13wn6vvvMZ2cQtMY2xSUoVFhGjHAfnyWV50n57LPN6p7dM89Iv/iFdPHF5pWFdkcfbT7H3a3rFwyffLLj78LiFwyjBsCUI/zAlVvwsRCAksstAOUcfOLozvD7x9+LV3eVPdg0Nkqvvppd1ejXL/szv/1tR6B56inpoYc63rMfCLeKk70SUKguHSs0BI3/sdjHNjm5tbNY3TxxnZTjrLo8/rj3390vfmE+7OrqpFGj8tunG79gKHVU1Cwpr9bEykiR5uZmQ5LR3Nzc6b3PP//ceOONN4zPP/+8BC3Lz44d5iOudXfsMIxNm+JZL8nHtZLV1xuGZBjV1eZzfX0OG6mrMz8c9lFX13kbO3caxubN2Y977jHXv+ce8/WDD5qvFy0ynzdvdt9OXZ35nEvbJMOYNMkwLrvM/Pm00zqWn3WW+fxv/2a25emnc9+P2zHYvLnjfa/vZj8OkmHMm5d9jKyHtY5zO3G00/L00+Y6Tz/d0Xa3dvvJ5XNRP+O3vtvfnXVM583r/J71+871+0YR9PcQxP7fQjHaWyR+5+9cUPmpEGGqMPZqTpAwl7NT8Umu+fPNgsa+fVJ1dY5dXXF0Z0SpCgwb5t0W59WVV11lzvdjr/xs2dJRxRk/XjrwwOyrv557znxI0m9+07H81782n3/0I/P5yivNcRiNjWb3SaG7dKxjZO+6Gz7cfM95Az7rcnhnJecrX5F++lPpH/6hY1mYbh47q3q0aZP5etMmqX9/9/3Z210qftONuLXNqgIOH57sLiJ71xw8EX4qQJhuKL9uLK/tRdk3kmXRoo7gs2+f+TpyACrnMQZubWto6Ag/y5eb7zc2dgwslsxgsX69/7bt3Ut+YzDcjkFjY+egYL0+4wyzLfb37d+jttYMdfauLeegcIvXLNluwnbzOLso7X8w5Xg1Ub7TjRTr0vFC3sakBPPNJUXiws+OHTt03XXX6amnntLnn3+ur3zlK7rvvvs0ZsyYUjetpPwCUJTgg8rnHONjvZZKfJWXF+sfcOuqpzi36zwp9O1rhh97Gpw3zwxNcVRy4hzgKnUeE+JVdbKuhrNfObR2rfm95s2Tzjkne7tu38uq9FnTBjz4oFlpyqXKVY4nZevvy3ou1qXjhdxPkeabS6JEhZ+PP/5YJ554oiZNmqSnnnpK/fr107vvvquDDjqo1E0rC24BiOADO7fBzWEugy8p6x/wxsbCnzCtcGDvYrO6mOKoZuXbVegMgl5XWYVpa9RuHmdYtIeuqMcm6KRcign7rN+99Vysq9TCTn2Q736QJVHh55ZbbtGgQYO0evXq9mVDhw71XL+lpUUtLS3tr/fs2VPI5pUFewCyZmUm+EDyv6qrrAKQV1WgFP8Xe+WV8Vac8u0qLGYQlIJDiDXmyApScY3zKcWEfc6/u2J167rt55BDzPATtB9mdc5ZosLPE088oVNOOUXnnXeeXnjhBR122GGaMWOGrvDo3166dKkWRr2MtgJwOwo4hbmcPZYAFEd3RqFCTm2tefsK6+cwnLfRCLufQgeTOI6Rs5vHTdgQks9tNtyUYm6gJHYRMatzzhIVft577z3dddddmjNnjm644Qb98Y9/1Pe//31169ZNF198caf1586dqznWP3YyKz+DBg0qeDtbW1u1YcMGNTY2qra2VuPHj1cmkyn4fi3cjgJOra3h5vGx3m9tzXFH5XwCqa2Vfvzj4HXq6syrnawA4/w/6zD7KddjYOfs5nETxwzPuSjnwfTlpBwnkEyIRIWftrY2jR07VkuWLJEkHXfccXr99dd11113uYafbt26qVu3bkVt45o1azRr1iz9z//8T/uygQMH6vbbb9fUqVMLvv9C3o7i5ptv1gMPPKCtW7dq4MCBWrBggb7zne/k32gUXJRzcdmN+Skme3CxX+1UboNzvUSpOoVZNyiEWNuYPDkZx6fSEBJz1qXUDYiitrZWRx11VNayESNGaNu2bSVqUbY1a9bo3HPPzQo+knmF2rnnnqs1a9YUdP9et6NwhqBcbdiwQT/5yU/05z//WdOmTdPFF1+s9957L/+GA+XMCkSFOrnH2U0Wpa1xfK9CH5tyU6yr1MrxargKk6jwc+KJJ+qtt97KWvb2229ryJAhJWpRh9bWVs2aNUuGYXR6z1o2e/Zstebcn+Cvtnagbr99RVbw2bhxo3r27Km//e2DWALQb37zG02ZMkWHH364rr76arW2tmpnvokKSLu0BYgkK9bvir+JgktU+PnBD36gF198UUuWLNE777yjhx9+WKtWrdLMmTNL3TRt2LChU8XHzjAMbd++XRs2bIh93zt3SkcdNU7vv7+pPfgYhqHZs2dr9uzZGjJkSFYF6Prrl+jAAw/0ffi10zAM/eu//quOOeYYffWrX439+wBAJ1RDEKNEjfk5/vjjtXbtWs2dO1f19fUaNmyYli9frgsvvLDUTVNjyEGRYdeL6sQTx2nNmp+3v37ggQe0bds2zZ07t32ZFYwuumi6Lr/8fN/tHXbYYZ7vXX755dq4caOeffZZVVdX59VuAAlTqhCSlIHkpUZIDKXKcOunqVB79uxRTU2Nmpub1bt376z3vvjiC23dulXDhg1T9+7dI2/7+eef16RJkwLXe+655zRx4sTI2w/yhz/8QRMmTFBzc7O6dOmir3zlK1qwYIEuv/zyWPfz6quvauTIkdqyZYuOPPLIwPXzPa4AAPidv3ORqG6vcjZ+/HgNHDhQVVVVru9XVVVp0KBBGj9+fEH2P3bsWGUyGTU0NOjmm2/WwQcfrEsvvdR13SVLvLu9vvQl8+HV7bV161ZJ0pFHHhnLIGoAAIotUd1e5SyTyej222/Xueeeq6qqqqyBz1YgWr58ecHm++nevbtGjhypNWvWaNWqVXryySfVpYt7tp0+fbrOP9+92+t//1fatUsaMMC922vChAnatGlTpDvEAwBQTgg/MZo6daoee+wx13l+li9fXvB5fsaNG6c77rhDp59+uk4++WTP9fr06aM+ffq4vveP/9hxVZhbuHnuued07bVz9ctfbuG2GTFasEDKZMLNsbNokTkJIcMfACA3hJ+YTZ06VWeddVZJZngeNWqUunbtqmXLluW1Hb87xL//frPeffctgk/MMplwt5Ww36YCAJAbwk8BZDKZggxqDvLQQw9pxowZoQYiB/G6Q/z48d/Vjh3fJfjELMx9tcLcnwsAEIzwk3BtbW1qamrSfffdp7feektr166NbdvFuEM83T0d/AIQwQcA4kP4Sbj169dr8uTJGj58uNasWaOamppYt2+/Q3wUVsXIY2hRO7p7srkFIIIPAMSL8JNwEydOVFtbW8G273aHeMm/+hPlSjC6ezqzH5PFi6V9+9Lz3QGgGJjnB57sIWbMGKlXr+zlfp+x1g1j/nzz5H7jjWbQsUtb8LHMny9VV5vBp7o6Xd8dAAqN8JNSQRMUOqs3zkDj9nl78Pnkk2jtcQtAaQ0+kvndreCzb1/nUAgAyB3hJ8W8ApBb8JE6D3a2f94ZfJzr3nln8AncHoC6dk138LG+e0uLd1UMAJAbxvw4pOVWZ15z+bgFH3uQcX7OfiWYW/CxjmeXLuEGNltaW9PZ3eNW7QozLgoAEB7h5+8OOOAASdJnn32mHj16lLg1xeEVgLyCj9fnrLzoVvH57LPPJEmzZh2gvXvDDWyWsrt70nKy9+vmIwABQHwIP3+XyWR00EEHadeuXZKknj17et6ktJL06SPt328Gmf37pX79zHt77dpl/tynj/TFF+6fa26W9u7tWPalL3WsbxiGPvvsM7377i69885BGjEiE2oeG6nj5G9fVuiTfannGwozvokABADxIPzY9O/fX5LaA1Ca/O1v0htvmA9Jqqkxg83WrdJf/2ouO+igjvX/+lcz/Njt3i3t2GGO1zn0UHPZO+8cpGnT+uvdd82Ttf0E/vzz0vjx2XP9TJ5cmu6eQsw3FCVQPfusNHFi8LrW+62twdsEALgj/NhUVVWptrZW/fr109/+9rdSN6fojj224wqjV1/tWL5ihXTHHdL3vy/NmNHxetw46cUXO7qoBg6UrPu5fv/70t69B0jKaOHC7GAxf74ZfJ59VnrhhY4T+dCh0u9/n92mYgWgQsw3FCVQPf98+AkcqfgAQJ6MFGlubjYkGc3NzaVuStmprzcMyTCqq83n+nr39ydPzn621rNeDxtmPjvXsT5v/zmT6VjXbZ9u+/dbJw5e+8l1/0GfK9b3AoAki/v8TfhBpxOw1wnZGXCc6zuDkVcAcq5jvQ7Tzrq6WL+6537CHA/DMNsTFFzsx8fe/rQEnzDHyFKs3zGAZCH85CHN4cfrBOR2Aq6vN4yJEzsvr6vrCD5WWPEKCpMnd2zDLRQ5g0+5hYCgSphzvaC2249BlM9VgrDfNU3HBEA0hJ88pDX81NW5Bwy3k41XV5V9fWudoC4yty4ue7dY2OpKXMcgavXB+n7V1cHrhwlIQcetktH9ByAfhJ88pDX8OE++bsHGMDpXJ+yftQKBtX5QMLB3X1jrVlV1rvoEhbG4eG3bGYr8gopfl0zYsUJhA1Ulins8FYD0IPzkIa3hxzA6n9QnTgwOPvbPuo1VCVPB8BrcXF3t3eVWzADkVqVyVsnsY5WiDMr2ep3Gyo8lyngqALAQfvKQ5vBjGN4nd7/g47UN58nLGaac61rjf5xVH69AYlVk4h786heAvMYfhQ0/9nWdAYeTfgdCIICoCD95SHv4MYzOAciqyEQNPvbuIrdwYF/XPtD5oIOyn/0CULEqQPb2h6kMBbXJOqlnMv7fxf67SNsVTmnu/gMQHeEnD4Qfk/Nkb52kw3wmqJLhrJ64XfZuDXq2XzKfS8hwE3Zgs7M7zt5Ov+pM2IG7zgHeXutHqbpVCio/AKIi/OSB8NMhjgkGgwKQX3eS25xB9fWGMXRo+BOiW7dY2OBkb1t1dUdocp6Yw4yBspY5pwewhzy/NpTjpf6FQvcfgFwQfvJA+DG5dXn5nYD8qinOk9fEiR0Bxm8cjb0dbvMAhQ0vbuuFrc54VR/s3VZRg5R9Xbd5kbzaWIjxTeUmbIgGACfCTx6SHn7imCnXa7K9fKoPfhWgoADh7CYK0+0V5mRpreMciO03Psk+MaO98hM2SIXt2krjyT6O3ymA9CL85CHp4Sdql45fpcVt/TgCkL2SEnZQq9t6cVQJnMHEL6j4jfnxOi5hgo9z+2GrSZUk379bACD85CHp4ccwcv8/6KCBtV4n+ijVJuvEbp/Dx2/sjH2/9vWsilUc40PsAcUrqDjXcasUec3943aJvxf78UkT7u0FIF+EnzxUQvgxjOhVEfsJPMx2/QKIF2dlw/6ZoIqTc72hQzuHDL9baQSdLK3P+rXBr8vNGYByuUop6HsAALwRfvJQKeHHMKJVRfL9P++gAGQPN15dQUHjXpwBw/6eV/dZlLE/Xbp4d2u5hR2vY5BL9cZ5FZhfu619Uy0BgA6EnzxUUvgxjOJWE8J0p9nXcVvf7fJ2t23bZ4P2qra4VamC2mwPO9b2onRb5VL58QtYfhWmsMHOGawAoBIRfvJQaeHHMIo7U65XN5XXydptmXPuG7/PObuqnNUhvwHaQV2AUY9ZLuOPgsKg2/v2eY6CKkTOUOVXLXK7gatfaKSSBKCcEH7yUGnhpxTjSJzdP1Yw8TrxOk/g9rlv/NptBRv7bTDsy/1mTg4bGqJWb8KOsfJ6z23sldeYI2s9t8qOc52o7YnabgAoNcJPHiop/ORSiYiL895VQdwqCX4VK68Kj/OWEV6zL1vdWV7vOUNG2OpNlPe9wqDX2Cd7e/0CkFfwcW4/bLUoyncFgFIh/OShUsJPLpWIuPcdZb4aZ/jxq1h5hbqqquxnr5O81xgYr0qM35iZMMezri78/EjWcQia88et28/+Oij45TIvEcEHQDkj/OShEsJPrpWIuPeda+Uk6LVbxcR58rePF3ILVVFP7FGrN26fHTbMf5yMPZjU1XUEH2f1zOsY26cQcKu4RTnOpawaAkAuCD95SHr4CXuSKsTJzG8MS9jKST5jZ6yTv73yYxidb0jqt2237rcwbQgSdQJJt3urebXBfuWb/eH2e3Be+Ramwsa8QwCSgPCTh6SHH7dKhN9AY7e5anK5iidMOAnbfeV1abpfuHJeYm4f8xMUwILCkbMNuV7l5DZBo9t3sHdJ1dd3XN3l1QXmrHo550EK6uryG1tVzCsFASAfhJ88JD38uCl0NSjM57xOyG7BJ2z3k/Mz1nw89m4mZ3eQfX9husPiFhRMnMHH+jmoC8xtvI/ba6+qF5UfAElH+MlDJYYfw8h9TEsYYWcatk7AVreU/TNRBwVbXT1eYcH62dpXJuMdsIo9nsUrmAR9F68uMLcKlzP0helKZMwPgCQj/OShUsOPYeQ2niZuYQbxhuEVluzbcQtbQbfQKBZnMAkKPlb7nHMY2b+L/ftYD+u7B3U5un2+lH8nABAV4ScPlRx+DKO0/0cfVMGI0oawY4ycVRX7z6XqznHrqnJ22fmNQfIKPha3K9+sfdj3H2V8lvN9AhCAckP4yUOlhx/DKP5YDrcqjfNqJvvJP2hAcX29YUyYENyFZw8R9iui6uuDB/IW6qahzu46ewXI/vuwvoPbRIzO7+L23e3btY6t39gf+2f9gpdzPQAoF4SfPKQh/BhGca/i8brM294FFlRxsEQZu+QW8uz78Qt/YU/wUYKAs+LjDIJWm4KOhdfVac6AY78xq/N9N9zbC0CSEX7ykIbwU8zKT9AgY+ftKOzruJ2M/cKKfSZm67POkOdWFfKqsEQJWlGPg9uxsB72G5eG2a71Xf2Ocym6OQGgmAg/eaj08FPMk2FQZcJabr+Lu1fXS5huGKt7x77MHvKcg4rtoSFqBSif4GP/PvZlznFJUfhdLef2e6BqA6DSEH7yUMnhJ44TeRR+t6FwDna234HdLfSECSf2/Tk/57Zft5uBhglAUY+XVzDxOhZhuv+C2hj1fQBIOsJPHio1/JTDydFrTIpbqHB2zQWFNnvAcKuw2Je7dbEFHYN8uwq9uvDcjoUzmIU9rkHrEoAAVDLCz98tWbLEkGTMmjUr9GcqMfyUw8nRuW2vAdf29ezrBFVfnBUk68amhpEdJJzBxwpO9kHEXl1CVnsymfy6jXI5Fn4KdWVasbYPAHGI+/zdVQm0adMmrVq1Sscee2ypm1Jyra1Sfb00f77/etb7ra3573PBAumFF6TJk83t2tuwaJG0b59UXW0+n3yy1NYmTZhgfk6Snn02ex3J/PyNN0qLF5vL7N9p/nzp+efNzx10kLR1q7nd3/9eeu45c5833pjdRmv9+nrz9Y03mj9bbbCz2pzJmN9lwwb/757JeB9v+7E4+eTs77loUfZ3stb349ZeL0F/A24ymY5j5/d56xhbxxMAEi2WCFVEn3zyiXHEEUcY69atMyZMmJD6yk8peI3VcVYz/K5+cnvtd4m+vevLfmNT536qq4Nnh3b7Lm7dan7fPahaUi6zTYdRDt2mAOAn9d1eF198sTF79mzDMIzA8PPFF18Yzc3N7Y/t27cTfmLiDEBeocZvHbdtBY27sdY96CDz2X6LB78xRPbPOtsY9X5gQe97zX1UziEi6HcTps10oQEolFSHn0ceecQ45phjjM8//9wwjODwU1dXZ0jq9CD8xMPrCiavYBQUSqzAEGZwsvMRZj/2zwdVenINQF7BJ+x2LaUIEkFjr6J+Pt/1AMCS2vCzbds2o1+/fsYrr7zSvozKT+k5u7bcAohfKLGfCO3Bwe0EaV/mvHmo9dqvwuRss/0zbqIGoKDgE3a7YdeJsl5Y1vZyvfKNLjQAhZDa8LN27VpDkpHJZNofkoyqqiojk8kY+/fvD9wGY37iZx+L41YFstZxu9rK7UToFYDcQpLzYb+bvN9Jtq6uI/gE3QLEr6pi/+72yRbDbCdMtaZUQSLf26PE0YUGAHapDT979uwxXnvttazH2LFjjWnTphmvvfZaqG0QfgrHPt4mzEnT70ToFoDcgo/9Z2vsjz18BJ2E870FiLWdoCCVz0m/EEHCr0vNeWwmT86tSy3fLjQAsEtt+HHD1V7lwW0cTr5jWtzuDG8tGzrUPQQ5rwKz2uZXbcr3pBzUhRbHST/uNocNVGG78YL2U4z7zAGobIQfG8JP6UUZ3Bx1EK9Vzclk/KtB9nY4A5DbbTGcbbAmSgxT5bCHqaDB03FWO+IOEkGByvndct1fvl1oAGAYhJ+8EH7i5RV0gpaHHcRrr6j4BR/n59zuJeZ3Erdv06/K4TUGyf5e2PmCchF3kPA6NnFVmqj8AIgL4ScPhJ/4+IUQv/fDDuK1rzN0aOcg5Pd5+x3d7eOEgi4/D/N97NvzCl9WO4O6i6JUwrzumZavsJMxRg1AjPkBECfCTx4IP/EICj5B6wWdYN3eswJFVVW0E7B1b6+gwGIfQzRxovs6QUHKfhVZmMvow97gtJCzRbtdsRYULIPEFaAAwEL4yQPhp7NcJtOrqzNP3GE+Z53knSdNr8pAHBUVv7aE6d7xWscelNw4Bz8PG9bxvb3GHlk/2wOXPWQUa7bouLrUwlb2CEAAoiD85IHw01nUcThxnrScY0LCVAuCKi9hv0fQ1VnOSkjYk7rXlWdB42vcwpDb1Wte+7SHy6iB1qpAxXXpfyn+lgBUNsJPHgg/7vxORn5dRm7biTonjHXC7dIlfDdJrreksNi7ptw4A1DUaob12usGrMOGdf6cMwRZnz3oIP/v4hYGo4aQuMbmcG8vAIVC+MkD4cdbUNCI45YNXp9xqzjk2n0S9Lkwt7awh4KgsU1B7bBCjNezWwXIXgkLW20K6lb0+46MzQFQ7gg/eSD8+POqYARVfvIJPn7jYHLtPgkb5NzWsy+Lepd4r/et+Yqsio89AAXtM9eQF+b4MDYHQFIQfvJA+AnmVY3J5yTstQ+vbYUdTG19xq37JOxYIa/Q5TwObu316+ax3rM+bwUg67UVgOwzWDsrTl7HPsox9/osY3MAJAnhJw+En3C8rvzJ5yTstY2o70dhbSvoKjFnIMgnMLht0zqeznuAOSdwDNOGXAYlOz8bR7gEgGIi/OSB8JPNrXLhV/mxTx6Yz0m4mBWHoMHN9v35dfFFDUDO9Z1jjbzuTO/VBWcY+V2OXqjbTDDIGUAxEH7yQPjJFlTJ8Rrzk+uJtNgnSmflJ6iLKuzYJuckiF7r2Ss49uNodXk5Q5BbQHP+DuKo/MTZjVWKQAsgfQg/eSD8dOYVcCxeg4TL/X5NUbvo4jqJB3WZeVV8/Lrm8pnhOY6uyqDQ6gx7Ua88A4AghJ88EH7cBQ0GLsZdy+MUNKDarb1xVKX8goZzMLNkzuFjP+5+c/Y4L2cPc2sMr89G/b2FWd+rW7Bc/0YAJAvhJw+En868Kj/Ok1axbrOQr3wrPPmEIK9bWNhfW1Uzr0BpX+7X1qAxSEEVvbgDkFt7yu1vA0ByEX7yQPjpzO2E7XWZu1VBcAaEoJN0sQa5xtF1FVf3l/MYOW8h4Ty+9mPkvPIrqA1e1Za452cKU1FLSrcogGQh/OSB8BPMOZjZ7YSXz7IoolZhJkzIv+vKeq9QoSGoYmMJe8d3a3v2QdiFvCVJmDFEhbqyDEB6EX7yQPjx5/Z/7V4BxCvsWCfSOLo84qrCxLnvfKsl1vEMO44mTDAJE0ji5FfdofIDoBAIP3mo1PCTyziVoK4rZ7UizADfoOW5iLsKk8++CzVOJs7jVKzQ4VbdKXYIA5AehJ88xHXwym1it1wqJF4/u60f5QRezlUYL36/T2eocLuJaJQ2O4Xt3rK26bfvYnU3uQWtYgRhAOlF+MlDXAevlN0xuezLawyIW7hx67oKO/C2kFWHQoasoG05b0/hdYzty4MClfMYR5k40e/9Qld+vH4PparOAUgHwk8e4jx4peyOidqmoMvUg6pBYb5LMaoOhTzBB1Uuwt4fLJdAnO/fkv39Yg52ti8LG5AJQAByQfjJQ9wHrxxL/V7/Z+52QnQOvvULO34nymKONylkyAo6dvleOu73fq5/S1F+31G2G3b9QgzgBgAnwk8eCjHguZDdMfm2yWu+Hq8TbC4Bppjfv5jda7nOap1PII56LIP2VawJDnPdLgCERfjJQ6Gu9ipm5SMsrwqJ3wk2alWlmJWvYoYs5x3Y46rGhGlr2L+lsG2I45Yk5TbAH0D6EH7yUMhL3ctpYregE6jb+9ayoDugO7fhNVA3znBSipBlHQev36dz315TB3jNhZPvVVthAok9AFnrl0s4B4AoCD95SEPlJ2zVwX6Cta9j7wILqioUY5BrPuNoorDCgX3AsF94tK/v1Y4wc+F4fZ+4/pbK6W8TAHJF+MlDpY/5CVshcZ4Q3SoTfuEmzrlpcvk+ua4XZhtB3UR++3ILkWHmwvH6HnH9LZVTVRIAckH4yUMlX+2VyxgQv5Dj9V4xv1uxx5rYBzm73fDVbwCx39xIUappcf8tUfkBUAkIP3mo1Hl+wu7LfnJ3q1IEBaBKHjPiDDjOoBM0V5Jf0AkKP4X6WyqnqiQA5IPwk4dKneE5l8GvQff2ci7v0qVyT5xeIcEZhKJ0iYXt9irU31I5VSUBIF+EnzxU6r29wrQhlxOs9T3DjBkph++ZK7ffp3UsnDM7W9/THpDs69u34xxU7nbFVSH+lsqpKgkAcSD85KFS7+oeJNcTrL27xm/MSKWeTL3u6eU2BijoUnfn5wsVFMutKgkAcYj7/N1VqHgLFoRfd/589+Xz5pnPN96Yvd6iReay+nrvzybRokXSvn1SdbX5PHmy+T0XLzZfW9/X7fvbl0kdx8xSyOPU2hrud2G939pauLYAQLki/MCV2wm8vt582E/mlRp87N/Lep3JdAQi6/taz27HxO1n+2cKIY6gCwAVL5b6UUKktdsrqqCBuWEu3U4qr+6gsLe8iGNwMwAgG91eKCivbix7haO6ujRti9uCBWY1J6gLb9Ei6dlnpaFDpfff7+gCk9yrJ86uMYuzSuT8rLM9fhYtMrusolR6AAAmwg+y+I0ZmT+/48ReXW2OA0rymJFMJrgLzwpEknTppR3rOQOQddzsx8frGErux83eHr8A5OySBABEQ/hBFquS4FaFcA4CttaLuwpRrAqIvRIzcaJ/8HG+19pqfsZ63zoO9uOzaJF/APJrj9d6lTrAHACKKpbOs4RgzE94QbMDu13uHVbQpffO8UVel4XHNYbGb3xO0Pa9jkM+bWOCQgDIxjw/eSD8RBMUcLxu+RB2u2FCRbHmFnJuL8rNW/3u+RVXACL4AEgzwk8eCD/RBd3TKpfKj/3zQbMQF7MC4nZLirCfKcRsyrm0BwAqEeEnD4SfaMJ26eR6gg/TvVPsCkiYW3m4tTOO9fJtDwBUKi51R9HYr/yyBtq6Xcad62zBbgN8vQb0eu07TmEHLFsKPZty1PYAAEKKJUIlBJWf/BSqChGme6fQFZByG2NTbu0BgFKi2ysPhJ/cFXr8iV+4KfS+y+3qqnJrDwCUGuEnD4Sf3ORahQh7N3lre263jih0BaSQA5YroT0AUA4IP3kg/ESXTxUil0va/QY7R9luGMUYsJzk9gBAuWDAM4rGbzbhMLMRB63jN4Oy1/Kw+w6j0AOWk94eAKhUhB+4CnMbhXwCkF/wCSOOABTllhjFuMqq3NoDAJUqp/Dz+eef66OPPtJhhx2Wtfz111/X0UcfHUvDUFpxViGcQcX+s9s+rH3b14t6k1AAALxUGYZhRPnAY489ph/84Afq06ePDMPQPffco6997WuSpNGjR6uhoaEgDY3Dnj17VFNTo+bmZvXu3bvUzUkdq9rTpYvU1hYuXMV901QAQPLEff7uEvUDixcvVkNDg/70pz/p3//933XppZfq4YcfliRFzFGRLV26VMcff7x69eqlfv366eyzz9Zbb71V0H0iPvPnmxP2tbWZz2G6bubPJ/gAAOIVOfz87W9/0yGHHCJJGjt2rNavX6+VK1eqvr5eVVVVsTfQ7oUXXtDMmTP14osvat26ddq/f7+mTJmivXv3FnS/iIfbjMVps2BB+O+9aBHBL/UaG80/gsbGUrcEqCiRw0+/fv306quvtr8++OCDtW7dOr355ptZywvh6aef1ne/+10dffTRGjlypFavXq1t27Zp8+bNBd0v8mcfQN3SYj7feGP6AlAmE+57W8crkylOu1CmGhulhQsJP0DMQg94/uSTT9SrVy898MAD6to1+2PV1dV65JFHdPXVV8feQD/Nzc2SpD59+ri+39LSopaWlvbXe/bsKUq7kM3tyrG4LlcvlAULzOARpl1RxiWF+d5hrrQDAOQh7IRAI0eONBobG2OZXCgObW1txhlnnGGcdNJJnuvU1dUZkjo9mOSweJI6Y3GhJxzkFhYIZfNm8w9i8+ZStwQoqZLN8HzZZZcZgwcPNt58882s5Q0NDcapp54aS2OimDFjhjFkyBBj+/btnut88cUXRnNzc/tj+/bthJ8iSvqMxYUObty8FIEIP4BhGCWc4fnee+/VwoULddJJJ+nxxx9Xv379NG/ePP3qV7/SmWeeWYiilKdrrrlGTzzxhNavX6+BAwd6rtetWzd169atiC2DXdJnLPbrooqja8q+/cWLzUHgdHWlWGNj57E91tQhblOI1NaaDwDRRU1LS5YsMbp3724ccMABxumnn25sLuL/kbS1tRkzZ840BgwYYLz99tuRP8+9vZCLQldo/O5ojxSpq+u40V2YR11dqVsMFE3JKj+NjY1aunSp7r33Xh111FHasmWLLrjgAo0ePbpwycxh5syZevjhh/XrX/9avXr10ocffihJqqmpUY8ePYrWDqRLISs0bpf/U/lJqauukpxV9IYG6YorpHvukZz/1lL1AXIXNiV1797dGDVqlPEf//EfhmEYxtNPP2307t3buPnmm2NJYWHIZfCyJGP16tWhPk/lB/mIu0LDmB8EYswPYBhGCSs/q1ev1gUXXND++pRTTtFzzz2n008/XR988IFWrFgRayhzYxR4BmnAS9wVmiRe/g8AlSL0JIf24GMZPXq0Nm7cqOeffz7ONgFlJe4JGv0GS8+fn94JIAGgWHK6q7vd0KFD9Z//+Z9xtAUoO3FXaMJcJUYFCAAKK+/wI0n/8A//EMdmgILIdbbmoAqNFD2gJP3yfxRZba1UV8fgZiBmsYQfoJxZ99OS/EOHPewUqkIT5UalVHyg2lrubgsUAOEHFS+X+2ktWECFBgAqFeEHqRB1tmYqNABQuQg/SA23AMQd1AEgfQg/SBXupwUAqDJSNHPgnj17VFNTo+bmZvXu3bvUzUEJdevWMWlhS0upWwMA8BP3+Tv0JIdApXCbrRkAkB6EH6RK3LM1AwCShzE/SA3up4WK0tgorVxp3g2eSRCBSAg/SIVCzNYMlFRjo7RwoXTmmYQfICLCDyoe99MCANgRflDxuJ8WAMCO8IOKx2zNAAA7wg8AlLvGRvNh19CQ/WxXWxs8DogB00gxwg8AlLuVK83BzW6uuKLzsrq64JInA6aRYoQfAGWrtbVVGzZsUGNjo2prazV+/HhlMplSN6v4rrrKDCl2a9ea92iZN08655zs9wgzgC/CD4CytGbNGs2aNUv/8z//075s4MCBuv322zV16tQStqwEnN1YjY3SG2+YPw8fLo0eXZp2AQnFDM8Ays6aNWt07rnnZgUfSdqxY4fOPfdcrVmzpkQtKxONjVKxjkFjo9mF5hxzBCQYlR8AZaW1tVWzZs2S2z2XDcNQVVWVZs+erbPOOiudXWC5yGfANGODUIEIPwDKyoYNGzpVfOwMw9D27du1YcMGTZw4sXgNKyVneHnzzY6ft2zpHGCc3WSFGDANJBjhB0BZaQzZvRJ2vYrgF14WLzYfds7w4jZguqHBDD7z5pmff/BBacQI8z0qPKhwhB8AZaU25Ik37HqJ19gonXCCGU4sW7Z0BJ5588xBz5a+faVjj83eht+8P9ZnR4zoPHDamgsIqDCEHwBlZfz48Ro4cKB27NjhOu6nqqpKAwcO1Pjx40vQuhLwq/pI7lWfU07JbV9u3WurVpk/5zqZIlCGCD8Aykomk9Htt9+uc889V1VVVVkBqKqqSpK0fPny9Ax29uuykrK7q6T8wghjg5AShB8AZWfq1Kl67LHHXOf5Wb58ebrm+Qmqrrh1V3mxV3asQdNbtpjPDQ2du9c+/li65hrz53vu6bwfqj5IKMIPgLI0depUnXXWWczwHCe3yo7VbeZW2bnyyo6fR49mMkVUDMIPgLKVyWTSczl7Lq68Mlr1xa8LbdmyztvyurxeYrwPEo3wAwBJFfWO7H6B5b//W7r2Wu/PTpuW/bqcxvtwh3pERPgBgKSprTXDR5wn+qlTzfBgV6iB1XFjFmpERPgBgKSprc296uJVJTnkEO8xPVdeKU2eHD1YUJFBmeLGpgCQJlaVJMoM2bmGl1z2BRQB4QcA0qwQXWhAmaPbC0Dy0b0SXlNT9nM+XWilkM8d6oG/I/wASD4GvIa3e3f2s59yrAoxCzViQPgBgErlViWxZnTesqVzpaRLF+nxxzsqaFGqQsWqyPjNVcQs1AiJ8AMAlcqvSrJ4ceebol55pXkj01wqaMWqyPiFJmahRkiEHwCoVG5VkrVrzdAzb550zjnZ7zU1ddzFPY59UZFBmSL8AEgWBryWJyoySBDCD4BkYcBreLl0e0XBVXZIKMIPgGSheyU8v26vWbOk44/Pfs+qqIWtoJXLVXbleFUayhrhB0Cy0L0Sntuxsu7Ovm2bdPvt7p9LWgUtaXMVoeQIPwBQCOXaJdS3r/l87rnmoGc7twpaU5O0Zo109tnR90VFBmWK21sAQCEE3deqsdGsVhT7vleHHGI+Dx/eUSmzP6Ts14ccYl4B1tYWfV9WRYbwgzJD5QcASqFcxsuE1dTUeSwQV9khoQg/AJKP7pXwcj1Wa9Z4zwGUtDFCSD3CD4DkY8BreLkeq6lTzfFLdlxlh4Qi/ABAvnKZeLEcdelizvVj7+Kynrdv7xgvZBk0yHzmKjskTOLCz4oVK7Rs2TI1Njbq6KOP1vLlyzV+/PhSNwtAmkWdePHKKzvm2IkyXqbQV5A9/rjZteXWveX1PYAEqjIMwyh1I8J69NFHddFFF2nFihU68cQTtXLlSt1777164403NHjw4MDP79mzRzU1NWpublbv3r2L0GIAFckZQrwqP1aX0KZN0e6Z5TVepqFBGjNG2ry5MJWWoO/h3GdTk/SNbxSuPcDfxX3+TlTl57bbbtNll12myy+/XJK0fPly/e53v9Ndd92lpUuXlrh1AFLDeaWW35VNgwaZEws+/XR2t1E5jpcJM4GkPfg1NZnvWc9AQiQm/Ozbt0+bN2/W9ddfn7V8ypQp2rhxo+tnWlpa1NLS0v56z549BW0jAHSye7d0223SlCnSE0907rJK2ngZe/ADEioxkxzu3r1bra2tOvTQQ7OWH3roofrwww9dP7N06VLV1NS0PwZZg/MAoNh273af9HDlyuxlxZr8MI79WJUs50BooMwlJvxYqqqqsl4bhtFpmWXu3Llqbm5uf2zfvr0YTQSA8Fat6hx+Fi6UXn3V7BpzPiT35VFDTNAM1Pkq1QzWQAiJ6fbq27evMplMpyrPrl27OlWDLN26dVO3bt2K0TwAlSqXy9itR11dx720oiqXSQVznRQxaTNYI1USE36qq6s1ZswYrVu3Tuecc0778nXr1umss84qYcsAJEKul4nnchm7NRngmWd2BKQtW8xn667qTU3ShRdKDz3kvu1ymFTQCn5nnmlWop591rxyTZJ+8QupVy/z57Vrze/Vt29HFxiDoFHGEhN+JGnOnDm66KKLNHbsWJ1wwglatWqVtm3bpunTp5e6aQDKXa6ViKuu6jy41yuErFzpPU/O4sXm87RpHcumTjWf3YLCIYd4D4Qu1iBpv+B3++0dP1vfzY45gFDGEhV+vvWtb+kvf/mL6uvr1djYqGOOOUa//e1vNWTIkFI3DUClCnP5t2XBAjMsWSEoyJo15vNNN5kDoqWOCtGzz3ZeP9dqSq5dd/bg19RktvG556T77pMuu8xcZ/Fiad488y7xzspPlLmNgCJKVPiRpBkzZmjGjBmlbgYAdGaFBisESWZ30LRp0te/Lq1f7/65DRvMh92113ZeL9dqStSuO2v8kFfwu+8+adIkacQIM/ycdFLnK76sC0y44zvKUOLCDwCUPbeT+3nnmeHnwQfN0CCZ1Z1rr+2onEhm5WfxYmnZMmny5Oxt5FpNidJ1Z7U/inIZnA2ERPgBUHly7eYpBGdX1fr1Zqix7++cczoCSEODGX4mT+4cShobc7vyKkrXnXN/zuNodctZz5L0T/9khjqr26upSfrZz6QnnyyvGayBvyP8AKg8uXbzhBXl8m9rLM9f/2o+r1oV7Yoz51VqxayY+B1H+yDna64xn63j2NBgBh8peTNYIxUIPwAqT6G7eaKEEGuen4MOirYPSynny7GOY1OT2bU1dar02mtmV92yZWbFx1p+yCFUdJAYhB8AlSfXbp5CsAYCH3GEOWA5SVdAWcexoaGjYmWNQ7K65U45JXg7uc6xBBRI4m5vAQCJdMghnSctdOs+y3VG5ais/XTpUvjbUBT6VhpARFR+ACAu9gHC1rw41sDgtWs71rPPiHzmmR2fsV8qX2jWfhoa8utW8xtcfsYZHcdBMr8zl7mjDBB+ACAuYQcIu82IXFdnVobK5Sq1sPy+85NPdgx8lsz5jrjMHWWA8AMgHYrRneQ2I7I1b8+8eeZyrxmRa2sLf5VaGEHTBAwaZI5damoyX599dvDgcmuixwcf7Dx3EVAChB8A6VCM7iS3Sow1b491Q2brZ7dB14W+Si2MsAHMGrjtDGCNjVKPHubP1vPnn3c8O8NVOVSvkDqEHwCIqlBXLxX6KrUwkz+ecIJZoZE6KlNRApg9PNlv4iox2zPKBuEHAKIq5dw7+ci3Wy1MALvqKvOyfquba8SI4levgACEHwBIi2J0q9XWdty7bMSI7G0y2zPKBOEHAArJOdC6kIOug7rjijX5Y7HmKgJyRPgBgEJyDrSOOr4lSpAol+64Yt+DDIiI8AMAfvwGCT/7rFlpse5tJcV/9VI5BAkqOagwhB8A8OM3SPjaa81n+/26KvHqpXwDGOEJZYbwAwB+/AYJz5tnzttjXdUkJe8EX4xgUg7VK8CG8AMAfvy6sYYPN5+dVzUVQ1zdcQQTpBDhBwCSiO44IGeEHwAIw15pefNN89m6Y3spbjpa6d1xQAERfgDAjzV3ziefSLfdlv2edXf2Uty2oZy64wp1uw+gQAg/AODHmjvn6aelCy80lzU1SWvWmLdxuPZabttQLvMLASERfgAgjEMOyQ44p5zS0d1Vyts2lFt3HJAAhB8ASDK3gc+l7I4DEoDwAwBJZh/4THccEArhBwAsfnPnuHUhNTUVvk1BnN1Yhe6Oi3qM6GZDGSL8AIDFb+4cty6kOXPSd9uGqMeIbjaUIcIPAFj85s7x6kJKU/CRcjtGQJkh/ACAxS/MlPKKrmIKmrOHY4QK0KXUDQAAxCyfm5Vac/Y4x/UAFYTKDwBUGm5WCvii8gMAAFKF8AMAfvLpQkoLjhEShm4vAPBTyV1Icc3ZU8nHCBWJ8AMAacWcPUgpwg8ApBVz9iClCD8AkFblOGdP0DxDQAwY8AwAKB/MM4QiIPwAAIBUIfwAAIBUIfwAADowZw9SgAHPAIAOxZyzJ655hoCICD8AgNJgniGUCOEHAFAazDOEEiH8AABKoxznGUIqMOAZAACkCuEHAACkCuEHAACkCuEHAFA+mGcIRcCAZwBA+SjmPENIrURUft5//31ddtllGjZsmHr06KEvf/nLqqur0759+0rdNAAAkDCJqPxs2bJFbW1tWrlypf7xH/9Rf/7zn3XFFVdo7969uvXWW0vdPAAAkCBVhmEYpW5ELpYtW6a77rpL7733XujP7NmzRzU1NWpublbv3r0L2DoAABCXuM/fiaj8uGlublafPn1812lpaVFLS0v76z179hS6WQAAoMwlYsyP07vvvquf/vSnmj59uu96S5cuVU1NTftj0KBBRWohAAAoVyUNPwsWLFBVVZXv46WXXsr6zM6dO/WNb3xD5513ni6//HLf7c+dO1fNzc3tj+3btxfy6wAAgAQo6Zif3bt3a/fu3b7rDB06VN27d5dkBp9Jkybpa1/7mn7+85+rS5do2Y0xPwAAJE9Fjfnp27ev+vbtG2rdHTt2aNKkSRozZoxWr14dOfgAAABICRnwvHPnTk2cOFGDBw/Wrbfeqqampvb3+vfvX8KWAQCApElE+HnmmWf0zjvv6J133tHAgQOz3kvolfoAAKBEEtF39N3vfleGYbg+AAAAokhE+AEAAIgL4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKQK4QcAAKRK4sJPS0uLRo0apaqqKr3yyiulbg4AAEiYxIWff/u3f9OAAQNK3QwAAJBQXUvdgCieeuopPfPMM/rVr36lp556KnD9lpYWtbS0tL9ubm6WJO3Zs6dgbQQAAPGyztuGYcSyvcSEn//93//VFVdcoccff1w9e/YM9ZmlS5dq4cKFnZYPGjQo7uYBAIAC+8tf/qKampq8t1NlxBWjCsgwDH3zm9/UiSeeqHnz5un999/XsGHD9PLLL2vUqFGen3NWfv76179qyJAh2rZtWywHL8327NmjQYMGafv27erdu3epm5NoHMt4cBzjw7GMD8cyHs3NzRo8eLA+/vhjHXTQQXlvr6SVnwULFrhWZuw2bdqkjRs3as+ePZo7d26k7Xfr1k3dunXrtLympoY/wpj07t2bYxkTjmU8OI7x4VjGh2MZjy5d4hmqXNLwc/XVV+uCCy7wXWfo0KFavHixXnzxxU5BZuzYsbrwwgt1//33F7KZAACggpQ0/PTt21d9+/YNXO+OO+7Q4sWL21/v3LlTp5xyih599FF97WtfK2QTAQBAhUnEgOfBgwdnvT7wwAMlSV/+8pc1cODA0Nvp1q2b6urqXLvCEA3HMj4cy3hwHOPDsYwPxzIecR/HRAx4dgo74BkAAMApkeEHAAAgV4mb4RkAACAfhB8AAJAqhB8AAJAqhB8AAJAqqQw/77//vi677DINGzZMPXr00Je//GXV1dVp3759pW5aIqxYsULDhg1T9+7dNWbMGG3YsKHUTUqcpUuX6vjjj1evXr3Ur18/nX322XrrrbdK3azEW7p0qaqqqjR79uxSNyWRduzYoWnTpunggw9Wz549NWrUKG3evLnUzUqc/fv3a968ee3nmMMPP1z19fVqa2srddPK3vr163XGGWdowIABqqqq0uOPP571vmEYWrBggQYMGKAePXpo4sSJev311yPvJ5XhZ8uWLWpra9PKlSv1+uuv6yc/+Ynuvvtu3XDDDaVuWtl79NFHNXv2bP3whz/Uyy+/rPHjx+vUU0/Vtm3bSt20RHnhhRc0c+ZMvfjii1q3bp3279+vKVOmaO/evaVuWmJt2rRJq1at0rHHHlvqpiTSxx9/rBNPPFEHHHCAnnrqKb3xxhv68Y9/HMt9lNLmlltu0d13360777xTb775pn70ox9p2bJl+ulPf1rqppW9vXv3auTIkbrzzjtd3//Rj36k2267TXfeeac2bdqk/v3761/+5V/0ySefRNuRAcMwDONHP/qRMWzYsFI3o+x99atfNaZPn561bPjw4cb1119fohZVhl27dhmSjBdeeKHUTUmkTz75xDjiiCOMdevWGRMmTDBmzZpV6iYlznXXXWecdNJJpW5GRTjttNOMSy+9NGvZ1KlTjWnTppWoRckkyVi7dm3767a2NqN///7GzTff3L7siy++MGpqaoy777470rZTWflx09zcrD59+pS6GWVt37592rx5s6ZMmZK1fMqUKdq4cWOJWlUZmpubJYm/wRzNnDlTp512mv75n/+51E1JrCeeeEJjx47Veeedp379+um4447TPffcU+pmJdJJJ52k3//+93r77bclSX/605/0hz/8Qd/85jdL3LJk27p1qz788MOsc1C3bt00YcKEyOegRNzeotDeffdd/fSnP9WPf/zjUjelrO3evVutra069NBDs5Yfeuih+vDDD0vUquQzDENz5szRSSedpGOOOabUzUmcX/7yl2poaNCmTZtK3ZREe++993TXXXdpzpw5uuGGG/THP/5R3//+99WtWzddfPHFpW5eolx33XVqbm7W8OHDlclk1Nraqptuuknf/va3S920RLPOM27noA8++CDStiqq8rNgwQJVVVX5Pl566aWsz+zcuVPf+MY3dN555+nyyy8vUcuTpaqqKuu1YRidliG8q6++Wq+++qoeeeSRUjclcbZv365Zs2bpwQcfVPfu3UvdnERra2vT6NGjtWTJEh133HG66qqrdMUVV+iuu+4qddMS59FHH9WDDz6ohx9+WA0NDbr//vt166236v777y910ypCHOegiqr8XH311brgggt81xk6dGj7zzt37tSkSZN0wgknaNWqVQVuXfL17dtXmUymU5Vn165dnZI4wrnmmmv0xBNPaP369ZFu0gvT5s2btWvXLo0ZM6Z9WWtrq9avX68777xTLS0tymQyJWxhctTW1uqoo47KWjZixAj96le/KlGLkuvaa6/V9ddf334++qd/+id98MEHWrp0qS655JISty65+vfvL8msANXW1rYvz+UcVFHhp2/fvurbt2+odXfs2KFJkyZpzJgxWr16tbp0qagiWEFUV1drzJgxWrdunc4555z25evWrdNZZ51VwpYlj2EYuuaaa7R27Vo9//zzGjZsWKmblEgnn3yyXnvttaxl3/ve9zR8+HBdd911BJ8ITjzxxE7TLbz99tsaMmRIiVqUXJ999lmnc0omk+FS9zwNGzZM/fv317p163TcccdJMseivvDCC7rlllsibauiwk9YO3fu1MSJEzV48GDdeuutampqan/PSpZwN2fOHF100UUaO3Zse8Vs27Ztmj59eqmbligzZ87Uww8/rF//+tfq1atXezWtpqZGPXr0KHHrkqNXr16dxkl96Utf0sEHH8z4qYh+8IMf6P/8n/+jJUuW6Pzzz9cf//hHrVq1iqp4Ds444wzddNNNGjx4sI4++mi9/PLLuu2223TppZeWumll79NPP9U777zT/nrr1q165ZVX1KdPHw0ePFizZ8/WkiVLdMQRR+iII47QkiVL1LNnT33nO9+JtqM4LkdLmtWrVxuSXB8I9rOf/cwYMmSIUV1dbYwePZrLs3Pg9fe3evXqUjct8bjUPXdPPvmkccwxxxjdunUzhg8fbqxatarUTUqkPXv2GLNmzTIGDx5sdO/e3Tj88MONH/7wh0ZLS0upm1b2nnvuOdd/Gy+55BLDMMzL3evq6oz+/fsb3bp1M77+9a8br732WuT9VBmGYeSX0wAAAJKDgS4AACBVCD8AACBVCD8AACBVCD8AACBVCD8AACBVCD8AACBVCD8AACBVCD8AACBVCD8AACBVCD8AEueRRx5R9+7dtWPHjvZll19+uY499lg1NzeXsGUAkoDbWwBIHMMwNGrUKI0fP1533nmnFi5cqHvvvVcvvviiDjvssFI3D0CZS+Vd3QEkW1VVlW666Sade+65GjBggG6//XZt2LChPficc845ev7553XyySfrscceK3FrAZQbKj8AEmv06NF6/fXX9cwzz2jChAnty5977jl9+umnuv/++wk/ADphzA+ARPrd736nLVu2qLW1VYceemjWe5MmTVKvXr1K1DIA5Y7wAyBxGhoadN5552nlypU65ZRTNH/+/FI3CUCCMOYHQKK8//77Ou2003T99dfroosu0lFHHaXjjz9emzdv1pgxY0rdPAAJQOUHQGJ89NFHOvXUU3XmmWfqhhtukCSNGTNGZ5xxhn74wx+WuHUAkoLKD4DE6NOnj958881Oy3/961+XoDUAkoqrvQBUnFNOOUUNDQ3au3ev+vTpo7Vr1+r4448vdbMAlAnCDwAASBXG/AAAgFQh/AAAgFQh/AAAgFQh/AAAgFQh/AAAgFQh/AAAgFQh/AAAgFQh/AAAgFQh/AAAgFQh/AAAgFQh/AAAgFT5/wtIrRvL+4fBAAAAAElFTkSuQmCC", "text/plain": [ "Figure(PyObject
)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "using Random; Random.seed!(1234);\n", "# Generate dataset {(x1,y1),...,(xN,yN)}\n", "# x is a 2-d feature vector [x_1;x_2]\n", "# y ∈ {false,true} is a binary class label\n", "# p(x|y) is multi-modal (mixture of uniform and Gaussian distributions)\n", "using PyPlot\n", "include(\"./scripts/lesson8_helpers.jl\")\n", "N = 200\n", "X, y = genDataset(N) # Generate data set, collect in matrix X and vector y\n", "X_c1 = X[:,findall(.!y)]'; X_c2 = X[:,findall(y)]' # Split X based on class label\n", "X_test = [3.75; 1.0] # Features of 'new' data point\n", "function plotDataSet()\n", " plot(X_c1[:,1], X_c1[:,2], \"bx\", markersize=8)\n", " plot(X_c2[:,1], X_c2[:,2], \"r+\", markersize=8, fillstyle=\"none\")\n", " plot(X_test[1], X_test[2], \"ko\") \n", " xlabel(L\"x_1\"); ylabel(L\"x_2\"); \n", " legend([L\"y=0\", L\"y=1\",L\"y=?\"], loc=2)\n", " xlim([-2;10]); ylim([-4, 8])\n", "end\n", "plotDataSet();" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Main Idea of Discriminative Classification \n", "\n", "- Again, a data set is given by $D = \\{(x_1,y_1),\\dotsc,(x_N,y_N)\\}$ with $x_n \\in \\mathbb{R}^M$ and $y_n \\in \\mathcal{C}_k$, with $k=1,\\ldots,K$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Sometimes, the precise assumptions of the (Gaussian-Categorical) generative model \n", "$$p(x_n,y_n\\in\\mathcal{C}_k|\\theta) = \\pi_k \\cdot \\mathcal{N}(x_n|\\mu_k,\\Sigma)$$ \n", "clearly do not match the data distribution." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Here's an **IDEA**! Let's model the posterior $$p(y_n\\in\\mathcal{C}_k|x_n)$$ *directly*, without any assumptions on the class densities." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Model Specification for Bayesian Logistic Regression\n", "\n", "- We will work this idea out for a 2-class problem. Assume a data set is given by $D = \\{(x_1,y_1),\\dotsc,(x_N,y_N)\\}$ with $x_n \\in \\mathbb{R}^M$ and $y_n \\in \\{0,1\\}$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- What model should we use for the posterior distribution $p(y_n \\in \\mathcal{C}_k|x_n)$?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "##### Likelihood\n", "\n", "- In Logistic Regression, we take inspiration from the generative approach, where the **softmax** function \"emerged\" as the posterior. Here, we **choose** the 2-class softmax function (which is called the [**logistic** function](https://en.wikipedia.org/wiki/Logistic_function)) with linear discrimination bounderies for the posterior class probability:\n", "$$\n", "p(y_n =1 \\,|\\, x_n, w) = \\sigma(w^T x_n) \\,.\n", "$$\n", "where $$\\sigma(a) = \\frac{1}{1+e^{-a}}$$ is the _logistic_ function.\n", "\n", "- Clearly, it follows from this assumption that $p(y_n =0 \\,|\\, x_n, w) = 1- \\sigma(w^T x_n)$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "

\n", "\n", "- (Bishop fig.4.9). The logistic function $\\sigma(a) = 1/(1+e^{-a})$ (red), together with the scaled probit function $\\Phi(\\lambda a)$, for $\\lambda^2=\\pi/8$ (in blue). We will use this approximation later in the [Laplace approximation](#gaussian-cdf).\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "- Adding the other class ($y_n=0$) leads to the following posterior class distribution:\n", "$$\\begin{align*}\n", "p(y_n \\,|\\, x_n, w) &= \\mathrm{Bernoulli}\\left(y_n \\,|\\, \\sigma(w^T x_n) \\right) \\\\\n", "&= \\sigma(w^T x_n)^{y_n} \\left(1 - \\sigma(w^T x_n)\\right)^{(1-y_n)} \\tag{B-4.89} \\\\\n", " &= \\sigma\\left( (2y_n-1) w^T x_n\\right)\n", "\\end{align*}$$\n", " - Note that for the 3rd equality, we have made use of the fact that $\\sigma(-a) = 1-\\sigma(a)$.\n", " - Each of these three models in B-4.89 are **equivalent**. We mention all three notational options since they all appear in the literature. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- For the data set $D = \\{(x_1,y_1),\\dotsc,(x_N,y_N)\\}$, the **likelihood function** for the parameters $w$ is then given by\n", "$$\n", "p(D|w) = \\prod_{n=1}^N \\sigma\\left( (2y_n-1) w^T x_n\\right)\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "- This choice for the class posterior is called **logistic regression**, in analogy to **linear regression**:\n", "$$\\begin{align*}\n", "p(y_n|x_n,w) &= \\mathcal{N}(y_n|w^T x_n,\\beta^{-1}) \\quad &&\\text{for linear regression} \\\\\n", "p(y_n|x_n,w) &= \\sigma\\left( (2y_n-1) w^T x_n\\right) &&\\text{for logistic regression}\n", "\\end{align*}$$\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ " \n", "- In the discriminative approach, the parameters $w$ are **not** structured into $\\{\\mu,\\Sigma,\\pi \\}$. In principle they are \"free\" parameters for which we can choose any value that seems appropriate. This provides discriminative approach with more flexibility than the generative approach. \n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "##### Prior \n", "\n", "- In *Bayesian* logistic regression, we often add a **Gaussian prior on the weights**: \n", "$$\\begin{align*}\n", "p(w) = \\mathcal{N}(w \\,|\\, m_0, S_0) \\tag{B-4.140}\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "### Some Notes on the Model\n", "\n", "- Note that for generative classification, for the sake of simplicity, we used maximum likelihood estimation for the model parameters. In this lesson on discriminative classification, we specify both a prior and likelihood function for the parameters $w$, which allows us to compute a Bayesian posterior for the weights. In principle, we could have used Bayesian parameter estimation for the generative classification model as well (but the math is not suited for a introductory lesson). \n", "\n", "- In the optional paper by [T. Minka (2005)](https://github.com/bertdv/BMLIP/blob/master/lessons/notebooks/files/Minka-2005-Discriminative-models-not-discriminative-training.pdf), you can read how the model assumptions for discriminative classification can be re-interpreted as a special generative model (this paper not for exam). \n", "- As an exercise, please check that for logistic regression with $p(y_n =1 \\,|\\, x_n, w) = \\sigma(w^T x_n)$, the **discrimination boundary**, which can be computed by\n", " $$\\frac{p(y_n\\in\\mathcal{C}_1|x_n)}{p(y_n\\in\\mathcal{C}_0|x_n)} \\overset{!}{=} 1$$\n", "is a straight line, see [Exercises](https://nbviewer.org/github/bertdv/BMLIP/blob/master/lessons/exercises/Exercises-Classification.ipynb). \n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Parameter Inference\n", "\n", "- After model specification, the rest follows by application of probability theory.\n", "\n", "- The posterior for the weights follows by Bayes rule\n", "$$\\begin{align*}\n", "\\underbrace{p(w \\,|\\, D)}_{\\text{posterior}} &\\propto p(w) p(D|w) \\\\ &= \\underbrace{\\mathcal{N}(w \\,|\\, m_0, S_0)}_{\\text{prior}} \\cdot \\underbrace{\\prod_{n=1}^N \\sigma\\left( (2y_n-1) w^T x_n\\right)}_{\\text{likelihood}} \\tag{B-4.142}\n", "\\end{align*}$$\n", "\n", "- In principle, Bayesian inference is done now! \n", "\n", "- Unfortunately, the posterior $p(w \\,|\\, D)$ is not Gaussian and the evidence $p(D)$ is also not analytically computable. (We will deal with this later)." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Application: the predictive distribution\n", "\n", "- For a new data point $x_\\bullet$, the predictive distribution for $y_\\bullet$ is given by \n", "$$\\begin{align*}\n", "p(y_\\bullet = 1 \\mid x_\\bullet, D) &= \\int p(y_\\bullet = 1 \\,|\\, x_\\bullet, w) \\, p(w\\,|\\, D) \\,\\mathrm{d}w \\\\\n", " &= \\int \\sigma(w^T x_\\bullet) \\, p(w\\,|\\, D) \\,\\mathrm{d}w \\tag{B-4.145}\n", "\\end{align*}$$" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- After substitution of $p(w | D)$ from B-4.142, we have closed-form expressions for both the posterior $p(w|D)$ and the predictive distribution $p(y_\\bullet = 1 \\mid x_\\bullet, D)$. Unfortunately, these expressions contain integrals that are not analytically computable. " ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Many methods have been developed to approximate the integrals in order to get analytical or numerical solutions. Here, we present the **Laplace approximation**, which is one of the simplest methods with broad applicability to Bayesian calculations." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### The Laplace Approximation\n", "\n", "- The central idea of the Laplace approximation is to approximate a (possibly unnormalized) distribution $f(z)$ by a Gaussian distribution $q(z)$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "- The Laplace approximation usually serves one or both of the following two purposes: \n", " 1. To approximate a posterior distribution without closed-form expression by a Gaussian distribution.\n", " 2. To approximate (part of) the integrand in an integral with purpose to get an analytical solution for the integral." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Example \n", "\n", "

\n", "\n", "- (Bishop fig.4.14a). Laplace approximation (in red) to the distribution $p(z)\\propto \\exp(-z^2/2)\\sigma(20z+4)$, where $\\sigma(a)=1/(1+e^{-a})$. The Laplace approximation is centered on the mode of $p(z)$. " ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Working out the Laplace Approximation \n", "\n", "- Assume that we want to approximate a distribution $f(z)$ by a Gaussian distribution $q(z)$.\n", "\n", "- Note that, if $q(z)$ is a Gaussian distribution, then $\\log q(z)$ is a second-order polynomial in $z$, so we will find the Gaussian by fitting a parabola to $\\log f(z)$. \n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "\n", "##### estimation of mean \n", "\n", "- The mean ($z_0$) of $q(z)$ is placed on the mode of $\\log f(z)$, i.e., \n", "\n", "$$z_0 = \\arg\\max_z \\left( \\log f(z)\\right) \\tag{B-4.126}$$ " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "##### estimation of precision matrix\n", "\n", "- Note that since $\\nabla \\log f(z) = \\frac{1}{f(z)}\\nabla f(z)$ and the gradient $\\nabla \\left. f(z) \\right|_{z=z_0}$ vanishes at the mode $z=z_0$, we can (Taylor) expand $\\log f(z)$ around $z=z_0$ as \n", "$$\\begin{align*}\n", "\\log f(z) &\\approx \\log f(z_0) + \\overbrace{\\left(\\nabla \\log f(z_0)\\right)^T (z-z_0)}^{=0 \\text{ at }z=z_0} + \\ldots \\\\\n", "&\\qquad + \\frac{1}{2} (z-z_0)^T \\left(\\nabla \\nabla \\log f(z_0)\\right) (z-z_0) \\\\\n", " &= \\log f(z_0) - \\frac{1}{2} (z-z_0)^T A (z-z_0) \\tag{B-4.131}\n", "\\end{align*}$$\n", "where the [Hessian matrix](https://en.wikipedia.org/wiki/Hessian_matrix) $A$ is defined by\n", "$$\n", "A = - \\nabla \\nabla \\left. \\log f(z) \\right|_{z=z_0} \\tag{B-4.132}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "##### Laplace approximation construction\n", "\n", "- After taking exponentials in eq. B-4.131, we obtain\n", "\n", "$$\n", "f(z) \\approx f(z_0) \\exp\\left( - \\frac{1}{2} (z-z_0)^T A (z-z_0)\\right) \n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- We can now identify $q(z)$ as\n", "$$\n", "q(z) = \\mathcal{N}\\left( z\\,|\\,z_0, A^{-1}\\right) \\tag{B-4.134}\n", "$$\n", "with $z_0$ and $A$ defined by eqs. B-4.126 and B-4.132. \n", "\n", "- All we have done up to now is approximate a function $f(z)$ by a Gaussian $q(z)$. This procedure is called the **Laplace Approximation**. Often, the required integrals (for Bayesian marginalization) can be approximately computed if we replace $f(z)$ by $q(z)$. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Bayesian Logistic Regression with the Laplace Approximation\n", "\n", "- Let's get back to the challenge of computing the predictive class distribution (B-4.145) for Bayesian logistic regression. We first work out the Gaussian Laplace approximation $q(w)$ to the [posterior weight distribution](#logistic-regression-posterior) \n", "$$\\begin{align*}\n", "\\underbrace{p(w | D)}_{\\text{posterior}} \\propto \\underbrace{\\mathcal{N}(w \\,|\\, m_0, S_0)}_{\\text{prior}} \\cdot \\underbrace{\\prod_{n=1}^N \\sigma\\left( (2y_n-1) w^T x_n\\right)}_{\\text{likelihood}} \\tag{B-4.142}\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "##### A Gausian Laplace approximation to the weights posterior \n", "\n", "- Since we have a differentiable expression for $\\log p(w | D)$, it is straightforward to compute the gradient and Hessian (for [proof, see optional slide](#gradient-hessian)):\n", "$$\\begin{align*}\n", "\\nabla_w \\log p(w | D) &= S_0^{-1}\\cdot \\left(m_0-w\\right) + \\sum_n (2y_n-1) (1-\\sigma_n) x_n \\\\\n", "\\nabla\\nabla_w \\log p(w | D) &= -S_0^{-1} - \\sum_n \\sigma_n (1-\\sigma_n) x_n x_n^T \\tag{B-4.143}\n", "\\end{align*}$$\n", "where we used shorthand $\\sigma_n$ for $\\sigma\\left( (2y_n-1) w^T x_n\\right)$. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- We can now use the gradient $\\nabla_w \\log p(w | D)$ to find the **mode** $w_{N}$ of $\\log p(w|D)$ (eg by some gradient-based optimization procedure) and then use the Hessian $\\nabla\\nabla_w \\log p(w | D)$ to get the variance of $q(w)$, leading to a **Gaussian approximate weights posterior**:\n", "$$\n", "q(w) = \\mathcal{N}\\left(w\\,|\\, w_{N}, S_N\\right) \\tag{B-4.144}\n", "$$\n", "with\n", "$$\n", "S_N^{-1} = S_0^{-1} + \\sum_n \\sigma_n (1-\\sigma_n) x_n x_n^T \\tag{B-4.143}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Using the Laplace-approximated parameter posterior to evaluate the predictive distribution \n", "\n", "- In the analytically unsolveable expressions for evidence and the predictive distribution (estimating the class of a new observation), we proceed with using the Laplace approximation to the weights posterior. For a new observation $x_\\bullet$, the class probability is now\n", "$$\\begin{align*}\n", "p(y_\\bullet = 1 \\mid x_\\bullet, D) &= \\int p(y_\\bullet = 1 \\,|\\, x_\\bullet, w) \\cdot p(w\\,|\\, D) \\,\\mathrm{d}w \\\\\n", " &\\approx \\int p(y_\\bullet = 1 \\,|\\, x_\\bullet, w) \\cdot \\underbrace{q(w)}_{\\text{Gaussian}} \\,\\mathrm{d}w \\\\\n", " &= \\int \\sigma(w^T x_\\bullet) \\cdot \\mathcal{N}\\left(w \\,|\\, w_N, S_N\\right) \\,\\mathrm{d}w \\tag{B-4.145}\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- This looks better but we need two more clever tricks to evaluate this expression. \n", " 1. First, note that $w$ appears in $\\sigma(w^T x_\\bullet)$ as an inner product, so through substitution of $a:=w^T x_\\bullet$, the expression simplifies to an integral over the scalar $a$ (see Bishop for derivation):\n", "$$\\begin{align*}\n", "p(y_\\bullet = 1 \\mid x_\\bullet, D) &\\approx \\int \\sigma(a) \\, \\mathcal{N}\\left(a\\,|\\, \\mu_a, \\sigma_a^2\\right) \\,\\mathrm{d}a \\qquad &&\\text{(B-4.151)}\\\\\n", "\\mu_a &= w^T_{N} x_\\bullet \\qquad &&\\text{(B-4.149)}\\\\\n", "\\sigma_a^2 &= x^T_\\bullet S_N x_\\bullet \\qquad &&\\text{(B-4.150)}\n", "\\end{align*}$$\n", " 2. Secondly, while the integral of the product of a logistic function with a Gaussian is not analytically solvable, the integral of the product of a Gaussian cumulative distribution function (CDF, also known as the [probit function](#scaled-probit)) with a Gaussian _does_ have a closed-form solution. Fortunately, \n", "$$\\Phi(\\lambda a) \\approx \\sigma(a)$$\n", "with the Gaussian CDF $\\Phi(x)= \\frac{1}{\\sqrt(2\\pi)}\\int_{-\\infty}^{x}e^{-t^2/2}\\mathrm{d}t$, $ \\lambda^2= \\pi / 8 $ and $\\sigma(a) = 1/(1+e^{-a})$. \n", " Thus, substituting $\\Phi(\\lambda a)$ with $ \\lambda^2= \\pi / 8 $ for $\\sigma(a)$ leads to \n", "\n", "$$\\begin{align*}\n", "p(y_\\bullet = 1 \\mid x_\\bullet, D) &= \\int \\sigma(w^T x_\\bullet) \\cdot p(w|D) \\,\\mathrm{d}w \\\\ \n", "&\\approx \\int \\underbrace{\\Phi(\\lambda a)}_{\\text{probit function}} \\cdot \\underbrace{\\mathcal{N}\\left(a\\,|\\, \\mu_a, \\sigma_a^2\\right)}_{\\text{Gaussian}} \\,\\mathrm{d}a \\\\ \n", "&= \\Phi\\left( \\frac{\\mu_a}{\\sqrt(\\lambda^{-2} +\\sigma_a^2)}\\right) \\tag{B-4.152}\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- We now have an approximate but **closed-form expression for the predictive class distribution for a new observation** with a Bayesian logistic regression model. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that, by [Eq.B-4.143](#Laplace-posterior-logistic-regression), the variance $S_N$ (and consequently $\\sigma_a^2$) for the weight vector depends on the distribution of the training set. Large uncertainty about the weights (in areas with little training data and uninformative prior variance $S_0$) increases $\\sigma_a^2$ and takes the posterior class probability eq. B-4.152 closer to $0.5$. Does that make sense?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Apparently, the Laplace approximation leads to a closed-form solutions for Bayesian logistic regression (although admittedly, the derivation is no walk in the park). " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Exam guide: The derivation of closed-form expression eq. B-4.152 for the predictive class distribution requires clever tricks and is therefore not something that you should be able to reproduce at the exam without assistance. You should understand the Laplace Approximation though and be able to work out simpler examples. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### ML Estimation for Discriminative Classification \n", " \n", "- Rather than the computationally involved Laplace approximation for Bayesian inference, in practice, discriminative classification is often executed through maximum likelihood estimation. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- With the usual 1-of-K encoding scheme for classes ($y_{nk}=1$ if $x_n \\in \\mathcal{C}_k$, otherwise $y_{nk}=0$), the log-likelihood for a $K$-dimensional discriminative classifier is \n", "\n", " $$\\begin{align*}\n", " \\mathrm{L}(\\theta) &= \\log \\prod_n \\prod_k {p(\\mathcal{C}_k|x_n,\\theta)}^{y_{nk}} \\\\\n", " &= \\log \\prod_n \\prod_k \\Bigg(\\underbrace{\\frac{e^{\\theta_k^T x_n}}{ \\sum_j e^{\\theta_j^T x_n}}}_{\\text{softmax function}}\\Bigg)^{y_{nk}} \\\\\n", " &= \\sum_n \\sum_k y_{kn} \\log \\big( \\frac{e^{\\theta_k^T x_n}}{ \\sum_j e^{\\theta_j^T x_n}} \\big)\n", " \\end{align*}$$\n", "\n", " " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ " \n", "- Computing the gradient $\\nabla_{\\theta_k} \\mathrm{L}(\\theta)$ leads to (for [proof, see optional slide below](#ML-for-LG)) \n", "\n", "$$\n", "\\nabla_{\\theta_k} \\mathrm{L}(\\theta) = \\sum_n \\underbrace{\\big( \\underbrace{y_{nk}}_{\\text{target}} - \\underbrace{\\frac{e^{\\theta_k^T x_n}}{ \\sum_j e^{\\theta_j^T x_n}}}_{\\text{prediction}} \\big)}_{\\text{prediction error}}\\cdot x_n \n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ " \n", "- Compare this to the [gradient for _linear_ regression](https://nbviewer.jupyter.org/github/bertdv/BMLIP/blob/master/lessons/notebooks/Regression.ipynb#regression-gradient):\n", "\n", "$$\n", "\\nabla_\\theta \\mathrm{L}(\\theta) = \\sum_n \\left(y_n - \\theta^T x_n \\right) x_n\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- In both cases\n", "\n", "$$\n", "\\nabla_\\theta \\mathrm{L} = \\sum_n \\left( \\text{target}_n - \\text{prediction}_n \\right) \\cdot \\text{input}_n \n", "$$" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- The parameter vector $\\theta$ for logistic regression can be estimated through iterative gradient-based adaptation. E.g. (with iteration index $i$),\n", "$$\n", "\\hat{\\theta}^{(i+1)} = \\hat{\\theta}^{(i)} + \\eta \\cdot \\left. \\nabla_\\theta \\mathrm{L}(\\theta) \\right|_{\\theta = \\hat{\\theta}^{(i)}}\n", "$$\n", " - Note that, while in the Bayesian approach we get to update $\\theta$ with [**Kalman-gain-weighted** prediction errors](https://nbviewer.jupyter.org/github/bertdv/BMLIP/blob/master/lessons/notebooks/The-Gaussian-Distribution.ipynb#precision-weighted-update) (which is optimal), in the maximum likelihood approach, we weigh the prediction errors with **input** values (which is less precise)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Code Example: ML Estimation for Discriminative Classification\n", "\n", "- Let us perform ML estimation of $w$ on the data set from the introduction. To allow an offset in the discrimination boundary, we add a constant 1 to the feature vector $x$. We only have to specify the (negative) log-likelihood and the gradient w.r.t. $w$. Then, we use an off-the-shelf optimisation library to minimize the negative log-likelihood.\n", "\n", "- We plot the resulting maximum likelihood discrimination boundary. For comparison we also plot the ML discrimination boundary obtained from the [code example in the generative Gaussian classifier lesson](https://nbviewer.jupyter.org/github/bertdv/BMLIP/blob/master/lessons/notebooks/Generative-Classification.ipynb#code-generative-classification-example)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "P(C1|x•,θ) = 0.4016963547369798\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAG2CAYAAACQ++e6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB3sUlEQVR4nO3dd3wT9f8H8Ffa0sEqo4xCKQ2yQTYKAlIQxQGCgOuLCoIMAb8Mv8gQLEtQcaDiAFQcOBAZTlQUEX6gX4GCCwpfsZVCi5TVUkah7ef3R7z0kt4ll+SSyyWv5+ORR2lyuXvnWrgXn3UWIYQAERERUZiIMLoAIiIiokBi+CEiIqKwwvBDREREYYXhh4iIiMIKww8RERGFFYYfIiIiCisMP0RERBRWGH6IiIgorDD8EBERUVhh+CEiIqKwYqrwU1xcjFmzZsFqtSIuLg6NGjXCvHnzUFpaanRpREREZBJRRhfgiSeffBKvvvoq3nrrLbRq1Qq7du3C/fffj/j4eEycONHo8oiIiMgETBV+fvjhBwwYMAC33HILACAlJQXvv/8+du3aZXBlREREZBamCj/du3fHq6++ioMHD6Jp06b4+eef8X//939YsmSJ4vZFRUUoKiqyf19aWopTp06hZs2asFgsAaqaiIiIfCGEwNmzZ1GvXj1EROgwYkeYSGlpqZg+fbqwWCwiKipKWCwWsXDhQtXt09LSBAA++OCDDz744CMEHtnZ2brkCYsQQsAkPvjgA0ydOhWLFy9Gq1atsHfvXkyaNAnPPvsshg0bVm5755af/Px8JCcnIzs7G1WrVg1k6UQh5ZtvvsHgwYPRokUL/Pjjj0aXQ0QhrqCgAA0aNMCZM2cQHx/v8/5MFX4aNGiA6dOnY/z48fbnFixYgFWrViEjI8Pt+wsKChAfH4/8/HyGHyIfHDx4EM2aNUOlSpVw9uxZdiMTkV/pff021VT38+fPl+vri4yM5FR3ogBr2LAhLBYLzp07hxMnTui67xEjRiA1NRU7duzQdb9ERBJTDXju378/Hn/8cSQnJ6NVq1bYs2cPnn32WYwYMcLo0ojCSkxMDOrVq4ejR48iMzMTtWrV0m3fO3fuxG+//Ybz58/rtk8iIjlThZ8XX3wRs2fPxrhx43D8+HHUq1cPY8aMwWOPPWZ0aURhx2q12sPPVVddpdt+X3rpJeTm5qJNmza67ZOISM5U4adKlSpYsmSJ6tR2Igocq9WK//u//0NmZqau+7322mt13R8RkTNTjfkhouBhtVoBQPfwQ0Tkbww/ROQVf4SfCxcuYOPGjfjvf/+r2z6JiJyZqtuLiIKHFH6ysrJ02+dff/2Fm2++GdWqVcPp06d12y8RkRzDDxF5JSUlBYAtsJSWluqy5HxpaSnatWuHypUr+7wvIiI1DD9E5JWkpCRERUXh0qVLyMnJQVJSks/7bNmyJfbs2aNDdURE6jjmh4i8EhkZieTkZAAc9ExE5sLwQ0Re44wvIjIjhh8i8pre4ee1115DamoqXn75ZV32R0SkhOGHiLymd/jZv38/vv/+e11nkBEROeOAZyLymt7h5/7778dVV12FJk2a6LI/IiIlDD9E5DW9w0/r1q3RunVrXfZFRKSG3V5E5DUp/Bw9ehSXLl0yuBoiIm0YfojIa7Vr10ZcXBxKS0uRnZ3t8/62b9+OH3/8EYWFhTpUR0SkjOGHiLxmsVjsKz3r0fU1dOhQdO3aFb///rvP+yIiUsPwQ0Q+0XPcT3JyMlJSUlCzZk2f90VEpIYDnonIJ3qGn61bt/q8DyIid9jyQ0Q+4SrPRGQ2DD9E5BOGHyIyG4YfIvKJXuHn999/R2pqKkaPHq1HWUREqjjmh4h8IoWf48eP49y5c6hUqZJX+zl69Ci+//57nD59Ws/yiIjKYfghIp9Uq1YN1apVw5kzZ5CVlYVWrVp5tZ82bdrggw8+QGxsrM4VEhE5YrcXEflMWuvHlxuS1q1bF3feeScGDBigT1FERCoYfojIZxz0TERmwm4vIvKZHuHn4MGDOH36NKxWK2rXrq1XaURE5bDlh4h8pkf4efrpp9GlSxcsW7ZMr7KIiBQx/BCRz/QIP/Hx8WjYsCESExP1KouISJFFCCGMLiJQCgoKEB8fj/z8fFStWtXocohCxr59+9CqVSvEx8fjzJkzRpdDRCFG7+s3W36IyGfSbK/8/Hyu00NEQY/hh4h8VrFiRdSpUwcAZ3wRUfBj+CEiXfg67qd///647bbbkJubq2dZRETlMPwQkS58WeiwtLQUn3/+OTZs2KBrTURESrjODxHpwpeWHyEE3n//fZw+fRo1a9bUuzQiIgcMP0SkC1/CT2RkJO688069SyIiUsRuLyLSBW9xQURmwfBDRLqQwk9WVhY8XT7s5MmT+PGzz5A5cSLgjwHPubnAnDn+2XcghcrnIDIYww8R6SI5ORkRERG4cOEC/v77b4/eu2XLFnTt3x/3vvCC/8LP3LnmDw2h8jmIDMbwQ0S6qFChApKSkgB43vUVERGBlHr10MDbg7NFhIg8wPBD5cyZA8yfr23b+fNt2xMB3o/7ue2225D56ad439sDe9IiYtaglJsL8Kav2pn150wBwfBD5URGAo895j4AzZ9v2y4yMjB1UfAzxaDnYOw60nKhzs0Fli8PWEku6zBDqAjGn3OgmeVnZQBOdadyZs+2fX3sMcfv5aTgM2+e8usUnjSHn9zc8v8gp6c7fpXk5QFffw385z+Alju+u9t3XJztz/v3274mJmrbrz9JF+pbby2rxflzSPUC5c8RELjPoVSrN/tYtgwYM8b4cy8Jxpp8pcfPKlSJMJKfny8AiPz8fKNLMYV584QAbF+1PE/01ltvCQDiuuuuc71hWprtl+ifxzxADADEV7Lnyj2+/NL2vpyc8vvbvdu2ze7d5fbt9pGWpv+J8JS8fkmwfg6lWo3Yh97HCERNWuTkqP+eeypYPpMO9L5+s+WHVCm1ALHFh1zR3PIzZoztf6P/2D5hAr764QfcBgCLFwP/+x8waBBQq5atxeOee4ATJ2z/i+3a1XXLTteuwKpVZa8lJADZ2cCoUcCKFbaWn3vusW3TooX6/4iNbglwOkf28wAAEycCFy6UnSOA/7P3VV6e41ejsLUmIBh+yCV5AFqwALh0icGH1FkrVgQAHD58GCUlJYhUGxDm1EXzn3nzMHDzZnRdtMj2/NSptot/hw7l37tunfrYl1Gjyj83ejRw7bVl7y0utv35wgXbV3n3krwuby9CeoQmpa47qV7A9hdx+XLb56pVy7suL6PDnZ486UYFlM/XiROOXymkMfyQW7NnlwWf6GgGHzOaM8c2MF3Lz27+fKCkxLtZfPWEQDSAS8XFOHLkCBo2bKjpfX369EGfGjWARYvcbzxokO2CLZeeXtayIw9My5bZQoIUljZuLHtNKShNmQJUqQIMHOj9zCp3oUnLhVqqW80rr9i+Si1BaWme/8C0hDs9QkUgLFtm+yxKlH7O3pyvYGeWn1WQYPght+bPLws+ly7ZvmcAMhdpBh/g+mcn79b0RkREBBoC+B9sXV9aw48DqctMGuAr/cOdkWH7euKEY1eP/B/wDh0cw8+cObagJO8ykjgHJcDW5XHjjUCTJmXhIy+vbD96XCw8vVCPHm07thTwAODBB20ByLnrzpPWHC3dPHqEikBclJ27CKV9KwViAIiIKH9s6fcrI6P8a2YICgyAHmH4IZecx/hI3wMMQGYSyBl8VpSFn9TUVLfbCyHw3//+F9VzctAYQKR0cOewsmBB+eelYKB2AXd10XIOSoDyxVgaa6TXGAznC3VeHvDSS8CnnzpeqPPybMdOSLB9L81SA2z/C5GTwoUnXXVaunnchQot46cCcVH29Oc8Z456TQsWlP2uOdekd1ehq2A4Zw4wfnxZ0Adcf05PA2Cwhzk/Y/ghVUoXQy0XUQpOrn52eg5kt/7zVdNaP7m5OHfoELr26AEAOHvHHaicmAg8/7y2g0ldWqNHu94uMdG2jVpXknQRklqb5LVLLQLyqea+tAQ4vzc93RZ8AMcLtasL9Ouv2746B0FpbJNe3IUKSYsWyuOzgOC8KA8caGvdk/v6a+Dtt4H77gNuuMHxtVatbF/1HozsKhh++mnZ74XEVTD0NACGOYYfUuTqYsgAZF66zuBT+V+rPfzs3u3YkqL0j/OyZTg3dy4aAigEUOnDD90fV627Sh5scnOBX35xbNWoXdvxPV98URZoNm4E3n237DX5iVBqcZK3BPirS0cKDe7G/0jkY5ucKdUZqG4evS7Kera6bNigHjreftv2kEtLA9q18+2YSlwFQ6CsRU0S5q01utJlwrxJhOo6P2lp2tfcmTfP/XIgWtfx4Xo/5iX97KKjffgZqqxDsxoQAEQ3LevQ5OTY1iCRP1assG2/YoXt+1WrbN/Pn6++Zonz2iierpEDCNGrlxAjR9r+fMstZc8PGGD7+sgjtlq+/NL74yidA2ktFlefTX4eACFmzXI8R9JD2sZ5P3quGfTll2XrLnm7jow379Nz3R6l3zvpnM6aVf416ecdyPWJvD2O/O8C1/lRxZafEKD3YNaSEm2tANLrJSXaa6XgoMsMPpX/tab887/WrNq1HWdXKf2v1ZNWAatVeTtpP/LugDFjbOv9yFt+MjLKWnF69AAqV3as77vvbA8A+Pzzsuc//tj29amnbF9Hj7aNw8jNtXWf+LtLRzpH8q675s1tr8nHAQFl0+GdW3KaNgVefBGoXr3sOS3dPHJS69HOnbbvd+4E6tZVPp68bqMkJtpabLT+3kmtgM2bm7uLSN41R6oYfkKA3oNZPRlryC4vc9JlBp/KxU2KKDl5eShq1QoxMTE+1+sxpdrS08vCz5Ilttdzc8sGFgO2YLF1q+t9y7uXXI3BUOrSyc0tHxSk7/v3t9Wi1lWYmGgLdfKuLedB4RKlgcRqtHbzOI9Pkf/CBONsIudA7KlATR33521MXAXAcKdL+1EAHTlyRAwdOlTUqFFDxMXFibZt24pdu3Zpem+odntJeDsK0sL590HX34/du0UpICrFxQkA4sCBA27fsnHjRjFgwADx7LPP2vfh0FQvNeNL3S3eNuFr6U6QjiF1sbnqYpJ3hygdR+kYvnY/SfuWurdWrXLdZSg9vvyyrJvKk24eOefut1Wr1I+ntg/5vjy9hYO/u3Dk3XlCBO72InofJ4S6uuTCutvr9OnT6NatG3r16oWNGzeidu3aOHToEKpVq2Z0aUGBt6MgdwIxg88CwFqvHn47dAiZmZlo2rSpy+3379+Pjz/+GBX/WR26HOl/8Lm5/v9frDStWN7FJnUx6TFjxteZT9L/5KXp72qzrLTU6mk3j3OLg3wgrqfnxl2rjBEL9kk/e+lroGapaV36wNfjkANThZ8nn3wSDRo0wMqVK+3PpaSkqG5fVFSEoqIi+/cFBQX+LC8o8HYUpCaQM/is9evbw487119/PV555RX7fcFUm+p97cbwxujRZUFDD77OfApkEATchxBpzJEUpPQa52PEgn3Ov3eBmjqudJxatWzhx91xuKqz93RpPwqQFi1aiEmTJokhQ4aIWrVqiXbt2only5erbp+Wlibwz8wT+SNUu73kpFk80dFGV0LBIGAz+P7pzvj3yJECgHjkkUe83JEf5OQIMWWK7aHW3eLcrSR14XjSjeBpl46e3RSe7Mu5m0eJUXeW1zILUGsXm54C1aWk9ThG/XwMENbdXn/++SdeeeUVTJkyBTNnzsRPP/2Ef//734iJicF9991XbvsZM2ZgypQp9u8LCgrQoEGDQJZsCN6OgpwFbAbfPy0T1iVLAGhc6DBQEhOBZ55xv01amm22k9QK4Pw/ay3HMcNtA5y7eZToscKzN7hgnzbBuICkSZgq/JSWlqJTp05YuHAhAKB9+/b4/fff8corryiGn5iYGGNmmhiIt6MgJYGewSd1YWkJPwcOHAAANGjQQH3cT6DIg4t8tpNZZsx4MrtHy7buQoi0j969zXF+Qg1DotcijC7AE4mJiWjZsqXDcy1atMDhw4cNqii4qA1mnTfP9vz8+cbWR+HDk/AzYsQING/eHF9++aW/y/KOFIj8dXHXczqyJ7Xq8bn8fW6CTaCmjnOKut+ZquWnW7du9v8lSg4ePOjdnaNDDG9HQcFEmohw8uRJnD17FlWqVFHdtmLFiqhatSpq1KgRoOqCjFm6yShwPyv+TvidqcLP5MmTcc0112DhwoW444478NNPP2H58uVYruW+NyFMy3R2BiAKJCnMnDp1CllZWbjyyitVt920aVMAKyPTYmsI6chU4adz585Yv349ZsyYgXnz5sFqtWLJkiUYOnSo0aUZirejoGBktVpx6tQpZGZmugw/ZDJGhRC2hmjDkKiJRQghjC4iUAoKChAfH4/8/HxUrVrV6HKIQtrtt9+Ojz76CEuWLMHEiRONLoeITEzv67epBjyT/82Zo31g9Pz5/I8YqdMy6PnEiRO47bbbMMqTe1EREfmI4YccSHeIdxeApHFGkZGBqYvMR0v4OX78ODZs2IC1a9cGqiwiInON+SH/0/sO8aTNnDm2IKnlfM6fbxu3FeytblrCT+3atfHKK6+gtLQ0UGURETH8UHmuAhCDj39ILW6A6/MqP//BTh5+hBCwWCzltklISMDYsWMDXRoRhTmGH1LEO8QHVii2uEnrbxUWFuLkyZNI0PMGoUREPmD4IVWBuEN8KHb3eCvUWtxiY2ORmJiI3NxcZGZmKoafv//+G2fOnEGdOnVQrVq1wBdJRGGJA57Jpdmzy26QGuHBb4vWmWAcYO1I6XYkZgw+EqnrKysrS/H1ZcuWoXnz5pg2bVoAqyKicMeWH3LJ+Q7xeo9LCcXuHl8FosUtUKxWK3bs2KE66NlisSA+Pj58b21BRIZgyw+pkoeOoiIgNdX2vKuWGuk9vXtrX0na1c1Xwy34SOQtbtHR5v3s7mZ8zZ49G2fOnMHChQsDWRYRhTmGnzDlbjFD59abOXNsgUbiKqj07g1s3uxZF1Wodff4yrnFTevCk8FG693dlWaCERH5jQgj+fn5AoDIz883uhTDzZsnBGD76uo15+2k76WH8/O9eyvvNy1N+Vhqx46MVK8v1KmdczOei82bNwsAomnTpkaXQkQmpvf1m+EnjCldVF0FH+dtpEd0tOvgo3YsVzVJ+w037s652QJQZmamACCio6NFSUlJudcnTpwohg8fLg4cOGBAdURkFgw/PmD4Kc/5oiq10Li72EqvWyzaW2q07lMeqMx2sfeF1vNjpnNy+fJlERkZKQCII0eOlHs9KSlJABA7d+40oDoiMguGHx8w/CjztptFaumRHr17l9/GubtLS2uSEd09WrvlhLBtl5am7/E9bRkzUwCyWq0CgNi2bVu519544w2xaNEikZeXZ0BlRGQWel+/OdWd3E6tVlqIcP78skHNJSW2r5s3A1YrkJICfPedbTvn2zbIj7VlC9Cjh+M2vXuXbaNlGrxe/HF7CU8WcNy82Tabzt220utaZ9IFA6vViszMTGRmZqJ79+4Or91///0GVUVEYU2XCGUSbPlxTepqch5ro9YyJLX8SO+zWh1bb1x1oUnvlbrLACFSUpTrClRrh97dTqHcmuOJESNGCABi7ty5RpdCRCbFbi8fMPyoky7AamNtnAOP8+Bm6Xt5AJJvozSQWh58fB0vpBe9BxyH4jgeT82fP18AEPfff7/D8+fPnxcZGRni+PHjBlVGRGbB8OMDhh9lWsf8OAcctZYg+VggpQDkvI3aeCGlOvUea6N2HK1jj7SMFZKfH3n94RB8hBDitttWCQAiNTXV4fkffvhBABDJycn25wL1MyYic2H48UE4hx+1i7TadPfU1PLPp6WVBR8prLjqEpP2oRSKnINPsIUAdy1hzttpHRyudt5C2QMPbBcARLVqDR2e/+6770R8fLxo3769ECK8zgkReYbhxwfhGn7S0pQDhquxOEpjddTG+rjqInLu4pJ3iwVyZpc3s7nUxkApba8lILk7b6EqJydHABBAhEhLu1Tu9ZKSEgYfInKJ4ccH4Rp+nC++7gYhy7ug5NvJA4S7YCDvvpC2ldYEUhvn488LoNq+1abiKwUVV10yWscKaQ1UoaS0tFTExsb+E4AO6TaeiojCB8OPD8I1/AhR/qKemuo++MjfqzRWRUsLhtrg5uho9S63QAYgpVYq51Yy+VglTwZlq30fbi0/QgjRvHlzAUAMG/ZNQFv8iCg0MPz4IJzDjxDqF3dXwUdtH84XL+cw5bytNP7HudVHLZBILTKBWExQaSC22utaBzc7BxxPBlGHoptuukkAECtWrJAF4jcFMFz8618fG10eEQU5hh8fhHv4EaJ8AJJaZDwNPvLuIqVwIN9WPtC5WjXHr64CUKBagOT1a2kZcleTFHwiI11/FvnPItRnOI0bN04AEDNnzhRCSOdopAAgFixYYHB1RBTsGH58wPBj43yxly7SWt7jriXDufVEadq7NOhZPmXem5ChxNu7x8vrdNU6o3Vws/MAb7XtPWl1M7OnnnpKABB33323/RxFRX0lgIXigQd2GF0eEQU5hh8fMPyU0WOBQXcByFV3ktKaQfPm2VZ51hp8lLrFtAYneW3R0eVXo5bfqV7rcZ2XB5CHPFc1BONUf72tWbNGABANGnQN6+4/IvIOw48PGH5slLq8XF2AXLWmOF+8UlPLAoyrcTTyOpTWAdIaXpS209o6ozQ2RwjHbitPg5R8W6V1kdRq9Mf4pmCya9cuYZvtVddtiCYicsbw4wOzhx897jyuttieL60PrlqA3AUI524iLd1eWi6WagOxXY1Pki/MKG/50RqktHZthePFfvr0k/+EH4jz588LIYQ4dOiQ+Pvvv0VxcXFYnhMi0o7hxwdmDz+edum4amlR2l6PACRvSdG6po3Sdlq72rTU5ByqXJ0btVDkqhYt9Ti3toXTRd52nkpFTExVAUDs27dPCCFEfHy8ACAyMjJk24XXuSEibRh+fGD28COE960i7gbWql3oPWltki7s8jV8XI2dkR9Xvp3UYuVukLUW8oCiFlSct1FqKVJb+0dpir8a+fkJJ9LvUNu2bQUA8fnnn4uSkhJRvXp1AUD8/fff9m1DvfuPiLzD8OODUAg/QnjeKiK/gGvZr6sAosa5ZUP+HnctTs7bpaSUDxmubqXh7mIpvddVDa663JwDkDeLFLr7HOFg4MCBAoBYunSp/bni4mJRWlpqYFVEZAYMPz4IlfAjhGetIr6OFXIXgOThRq0ryN24F+eAIX9NrfvMk7E/ERHq3VpKYUftHHjTeuM8C8xV3dKx9RjfFWwmT54sAIiHH37Y6FKIyGQYfnwQSuFHiMC2JmjpTpNvo7S90vR2pX3LV4NWa21RaqVyV7M87Ej786TbypuWH1cBy1ULk9Zg5xysgtkLL7wgAIhBgwYZXQoRmQzDjw9CLfwIEdgbZap1U6ldrJWec177xtX7nLuqnFuHXA3QdtcF6Ok582b8kbswqPS6fJ0jdy1EzqHKVWuR0g1cXYVGf7QkffLJJwKAaN++vdi7d68YPny4eOqpp/Q/EBGFHIYfH4Ra+DFiHIlz948UTNQuvM4XcPnaN67qloKN/DYY8uddrZysNTR42nqjdYyV2mtKY6/UxhxJ2ym17Dhv42k9ntatl19//VUAENWqVROrV68WAET37t31PxARhRyGHx+EUvjxpiVCL873rnJHqSXBVYuVWguP8y0j1FZflrqz1F5zDhlaW288eV0tDKqNfZLX6yoAqQUf5/1rbS3y5LP66uzZs0Ja6+eHH34QCxcuFCtXrvTPwYgopDD8+CBUwo83LRF6H9uT9Wqcw4+rFiu1UGexOH5Vu8irjYFRa4lxNWZGy/lMS9O+PpJ0Htyt+aPU7Sf/3l3w82ZdokCF51q1agkAYs+ePf49EBGFFIYfH4RC+PG2JULvY3vbcuLue6UWE+eLv3y8kFKo8vTC7mnrjdJ7rVbX42TkwSQtrSz4OLeeqZ1j+RICSi1unpxnI1sNr7rqKgFArFu3zv8HI6KQwfDjA7OHH60XKX9czFyNYdHacuLL2Bnp4i9v+RGi/A1JXe1bqftNSw3ueLqApNK91dRqkM98kz+Ufg7OM9+0tLAFet2hO++8UwAQc+bMEX///be4dOlSYA5MRKbG8OMDs4cfpZYIVwONldaq8WYWj5ZworX7Sm1quqtw5TzFXD7mx10AcxeOnGvwdpaT0gKNSp9B3iU1b17Z7C61LjDnVi/ndZDcdXW5GlsVyJmCkunTpwsAonHjxgKAeP755wN3cCIyLYYfH5g9/Cjxd2uQlvepXZCVgo/W7ifn90jr8ci7mZy7g+TH09Idpjd3wcQ5+Eh/dtcF5vxQC0RqrV7B1PKzbNkyAUDUqVNHABBvv/12YA5MRKbG8OODUAw/Qng/pkULrSsNSxdgqVtK/h5PBwVLXT1qYUH6s3SsyEj1gBXI8SxCqAcTd59FrQtMqYXLOfRp6UoMljE/X3/9tQAgWrZsKYqLi8Xly5f9f1AiMj2GHx+EavgRwrvxNHrTMohXC7WwJN+PUthydwuNQHEOJu6Cj1Sf8xpG8s8i/zzSQ/rsrs6TnKsuwkCdq4MHDwoAomLFirynFxFpxvDjg1AOP0IY+z96dy0YntSgdYyRc6uK/M9G3UBUqavKucvO1RgkteAjUZr5Jh1DfnxPxmc5v+7Pc3bx4kVhsVgEAHHs2DH/HYiIQgrDjw9CPfwIEfixHEqtNM6zmeQXf3cDiufNE6JnT/ddePIQIZ8RNW+e+4G8/rppqHN3nbwFSP7zkD6D0kKMzp9F6bPL9yudW6UuP6X3ugpeztv5S/369QUAcdddd4mzZ8/670BEFDIYfnwQDuFHiMDO4lGb5i3vAnPX4iDxZOySUsiTH8dV+NN6gfckCDi3+DgHQakmd+dCbXaac8CR35jV+XUlwXBvL0nXrl2FtNLz+fPn/XcgIgoZDD8+CIfwE8iWH3eDjJ1vRyHfRuli7CqsyFdilt7rHPKUWoXUWlg8CVqengelcyE95Dcu1bJf6bO6Os9GdHP64u677xYARM+ePY0uhYhMguHHB6EefgJ5MXTXMiE9L7+Lu1rXi5ZuGKl7R/6cPOQ5DyqWhwZPW4B8CT7yzyN/znlckidczZZT+jn4s9VGD2lpaQKAGDVqlNGlEJFJMPz4IJTDjx4Xck+4ug2F82Bn+R3YlUKPlnAiP57z+5SOq3QzUC0ByNPzpRZM1M6Flu4/dzV6+nqwefPNNwUA0adPH6NLISKTYPjxQaiGn2C4OKqNSVEKFc5dc+5CmzxgKLWwyJ9X6mJzdw587SpU68JTOhfOwUzreXW3rZkC0MaNGwUA0ahRI6NLISKTYPj5x8KFCwUAMXHiRM3vCcXwEwwXR+d9qw24lm8n38Zd64tzC5J0Y1MhHIOEc/CRgpN8ELFal5BUT2Skb91G3pwLV/w1My1Q+1fy1FNPCQDCYrGI4uJi33dIRCFP7+t3FExo586dWL58Odq0aWN0KYYrKQHmzQNmz3a9nfR6SYnvx5wzB/j+e6B3b9t+5TXMnw9cugRER9u+XncdUFoK9Oxpex8AbN7suA1ge/9jjwELFtiek3+m2bOBLVts76tWDcjMtO3322+B776zHfOxxxxrlLafN8/2/WOP2f4s1SAn1RwZafss27a5/uyRkernW34urrvO8XPOn+/4maTtXVGqV4273wElkZFl587V+6VzLJ1PX1SoUAEAIITA0aNHkZyc7PtOiYg8oUuECqCzZ8+KJk2aiE2bNomePXuGfcuPEdTG6ji3Zria/aT0vasp+vKuL/mNTZ2PEx3tfnVopc+i1K3m6rO7ay0JltWmtTCi27RRo0YCgNiyZYt+OyWikBX23V733XefmDRpkhBCuA0/Fy9eFPn5+fZHdnY2w49OnAOQWqhxtY3SvtyNu5G2rVbN9lV+iwdXY4jk73Wu0dP7gbl7XW3tIzMGIE9q9qQLrVGjPgKAWLlypaelElEYCuvw8/7774vWrVuLCxcuCCHchx9pSq3zg+FHH2ozmNSCkbtQIgUGLYOTnR9ajiN/v7uWHm8DkFrw0bpfiRFjcdyNvfL0/a63GyUAiLRgn5dPREEhbMPP4cOHRe3atcXevXvtz7Hlx3jOXVtKAcRVKJFfMOXBQelCKn/O+eah0veuWpica5a/R4mnAchd8NG6X63beLKdVtL+vJ35pu2cPSuSk7sJAOK+++7zvWgiCnlhG37Wr18vAIjIyEj7A//MGImMjNQ0a4RjfvQnH4uj1AokbaM020rpQqkWgJRCkvNDfjd5VxfhtLSy4OPuFiCuWlXkn12+2KKW/WhprTFiLI4Qvt8exV0XWlJSF3srbI8ePXwvmIhCXtiGn4KCAvHrr786PDp16iTuuece8euvv2raB8OP/8jH22i5aLq6cCsFIKXgI/+zNPZHHj7cXYR9vQWItB93QcqXkKLHWBxnrrrUnM9N797edam56kJbvny5uOeeewQAkZSU5PnOiSjshG34UcLZXsFBaRyOr2NalO4MLz2XkqIcgpxngUm1uWpt8rX1xF0Xmh6tM3rXrDVQae3Gc3ccpZD5999/21tuL1686N0BiChsMPzIMPwYz5PBzZ4O4pVacyIjXbcGyetwDkBKt8VwrkFaKFFLK4c8TLkbPK1nt5RerVXO+1Or1fmzeXs8tS600tJSUbFiRQFAHDx40PsPQkRhgeHHBww/+lILOu6e1zqIV96i4ir4OL9P6V5iri7i8n26auVQG4Mkf03rekHe8HUsjjO1c6NXS5NSYCsuLhZ///23KCoqEq1atRIAxFdffaXPByKikMXw4wOGH/24CiGuXtc6iFe+TUpK+SDk6v3yO7rLxwm5m36u5fPI96cWvqQ63XUXedISpnbPNF9pXYzR0wCkFqD+858jAoCIiooSN998swAgXn31VX0+DBGFLIYfHzD86MNd8HG3nbsLrNJrUqCwWDy7AEv39nIXWORjiFJTlbdxF6Tks8i0TKPXeoNTf64WrTRjzV2wdMf1z/d3AUAkJCSICRMmCABi2rRpvn4MIgpxDD8+YPgpz5vF9NLSbBduLe+TLvLOF021lgE9WlRc1aKle0dtG3lQUuI8+NlqLfvcamOPpD/LA5c8ZARqtWi9utS0tewVi5kzz4hnnnlGABB33nmnbwclopDH8OMDhp/yPB2Ho+cYFucxIVq6W9y1vGj9HO5mZzm3hGjtrlObeeZufI1SGFKavaZ2THm49DTQSi1Qek391/ozuOuudQKAuOqqq7w7IBGFDYYfHzD8KHN10XLVZaS0H0/XhJEuuBER2seZeHtLCom8a0qJcwDSGnycA43aDVit1vLvcw5B0nurVXP9WZTCoKchxFVrmCc8DV1jxuwRAEStWrU8PxgRhRWGHx8w/KhzFzT0uGWD2nuUWhw8DRxa36fl1hbyUOBubJO7OqQQo/ZVqQVI3hKmtbXJXbeiq8/obYj0xtdffy0mT54s1q5dK86cOSOklZ7Pnj2r/8GIKGQw/PiA4cc1tRYMdy0/vgQfV+NgvO2K0xrklLaTP+fpXeLVXpfWK5JafOQByN0xvQ15Ws6P3j9TLR577DEBQIwdO1YIIUT16tUFAM2rtBNReGL48QHDj3tqrTG+XITVjqG2L62DqaX3KHW1aR0rpBa6nM+DUr2uunmk16T3SwFI+l4KQPIVrJ1bnNTOvSfnXO29voRLX2zatEk88sgjYt26dUIIITp06CAAiE8++US/gxBRyGH48QHDjzZqM398uQir7cPT1z0h7cvdLDHnQOBLYFDap3Q+ne8B5ryAo5YavBmU7PxePcKlXgYPHiwAiOeff95/ByEi02P48QHDjyOllgtXLT/yxQN9uQgHssXB3eBm+fFcdfF5GoCct3cea6R2Z3q1LjghfJuOrvfq0BJvlkqQe/jhhwUAMWnSJH0LI6KQwvDjA4YfR+5actTG/Hh7IfX1Qukp55Yfd11UWsc2OS+CqLadvAVHfh6lLi/nEKQU0Jx/Bnq0/OjZjeVpoJ0x45QoKiqyP7906VIBQAwYMEC/oogo5DD8+IDhpzy1gCNRGyTsjwupnjztotOrVcpdl5lai4+rrjlfVnjWo6vSXWh1DnuuZp5J9/P69ttvhRBCfP755wKAaNu2rfaCiCjsMPz4gOFHmbvBwIG4a7me3A2oVqpXj1YpV0HDeTAzYFvDR37eXa3Z4zydXcutMdTe6+nPTcv2at2Czu+tX7++ACB2794thBBi3759AoCoWrWqKC0t1VYQEYUdhh8fMPyUp9by43zRCtRtFnzlawuPLyFI7RYW8u+lVjO1QCl/3lWt7sYguWvR0zsAKdWj9J7i4mJx8uRJcenSJSGEEOfPnxfSWj8nT57UVgwRhR2GHx8w/JSndMFWm+YutSA4BwR3F2l/zhZyPpavXVd6dX85nyPnW0g4n1/5OXKe+eWuBrXWFr3XZ9LSouZNt2jdunUFALFr1y5thRBR2GH48QHDj3vOg5mVLni+POcJT1thevb0vetKes1focFdi41E6x3fpf3JB2H785YkWsYQeTogvmvXrgKAWLNmjfZCiCisMPz4gOHHNaX/tasFELWwI11I9egO06sVRs9j+9paIp1PreNotAQTLYFET65ad1y9dvToUTF58mTx9NNPO+zvX//6lwAgnnrqKf8UTESmx/Djg1ANP96MU3HXdeXcWqFlgK+7572hdyuML8f21zgZPc9ToGbhKbXuuDtfO3bsEACEVbqvxz8effRRAUA8+OCD/i2aiEyL4ccHep28QK9Xo+UYnraQqP1ZaXtPLuDB3AqjxtXP0zlUKN1E1JOanWnt3pL26erY/lrIUKkO56Cl5Wd06NAh8cgjj4j58+c7bLNixQoBQNx4443+LZyITIvhxwd6nTwju2O8OZbaGBClcKPUdaV14K0/Wx38GbLc7cv59hRq51j+vLtA5XyOPVk40dXr/m75Ufs5+NLC9c033wgAolmzZv4pmohMj+HHB3qePCO7Yzytyd00dXetQVo+SyBaHfx5gXfXcqH1/mDeBGJff5fkrwdysLP8Oa0BWWmbQ4cOCQAiJiZGlJSUaC+IiMIGw48P9D55gRjz4mtN0vdKF0Tnwbeuwo6rC2Ugx5v4M2S5O3e+Th139bq3v0ue/Lw92a/W7T0ZwD1r1jkxa1ZRuecvXbokIiIiBACRk5OjrSAiCisMPz7wx4Bnf3bH+FqT2no9ahdYbwJMID9/ILvXvF3V2pdA7Om5dHesQC1wqHW7yZMnCwBi5syZ5V5r2LChACC2b9+urSgiCisMPz7w12yvQLZ8aKXWQuLqAutpq0ogW74CGbKc78CuV2uMllq1/i5prUGPW5LoNcD//vvvFwDEwoULy72WmpoqAIhVq1ZpL4yIwgbDjw/8OdU9UDNttHB3AVV6XXrO3R3QnfehNlBXz3BiRMiSzoPaz9P52GpLB6ithePrrC0tgUQegKTtjQznxcXF4tSpU6KgoKDca1Iwcp4JRkQkBMOPT8Kh5Udrq4P8AivfRt4F5q5VwZdBrt5+Hn8cQ4iycCAfMOwqPMq3V6tDy1o4ap9Hr9+lYPrddGXevHkCgBgxYoTRpRBREGL48UGoj/nR2kLifEFUaplwFW70XJvGm8/j7XZa9uGum8jVsZRCpJa1cNQ+h16/S8HUKqnmnXfeEQBEr169jC6FiIIQw48PQnm2lzdjQFyFHLXXAvnZAr2YpHyQs9INX10NIHa1NpInrWl6/y4FU8vP3Llzxdy5c8WxY8fKvbZt2zYBQKSkpBhQGREFO4YfH4TqOj9ajyW/uCu1UrgLQEaPGfEn54DjHHTcrZXkKui4Cz/++l0KplZJIYSoUqWKACAOHjxY7rUjR44IACIyMlJcvnzZgOqIKJgx/PggVFd49mbwq7t7ezk/HxER+sFHraVH67R3b7q9/PW7FEytkkIIUVpaKqZNmyZGjx4tzpw5U+71kpISERMTIwCIP//8M7DFEVHQY/jxQaje20tLDd5cYKXPqWXMSDB8Tm8p/Tylc+G8srP0OeUBSb69fD/Og8qVZlz543cpmFolPdG0aVMBQGzevNnoUogoyDD8+CBU7+rujrcXWHl3jasxI8F6MfWV2j29lMYAuZvq7vx+fwXFYGuV9ETfvn0FAPH6668bXQoRBRm9r99RoJA3Z472bWfPVn5+1izb18cec9xu/nzbc/Pmqb/XjObPBy5dAqKjbV9797Z9zgULbN9Ln1fp88ufA8rOmcSf56mkRNvPQnq9pMR/tcgVFxejtLQU0dHRqttYrVYAQGZmZmCKIqKwxfBDipQu4PPm2R7yi3moBh/555K+j4wsC0TS55W+Kp0TpT/L3+MPegRdf9i8eTP69u2Lq6++Gj/++KPiNgw/RBQoDD9UjlprjlIACvXgA9i+btkCbN5cFoDmz1cPQM7Bx/kcyt8TLk6fPg0AiI2NVd2G4YeIAoXhhxyoBR/5Bd5Fz4WpzJljCzPuuvDmz7cFn5QUICurrAsMUA4xzl1jEueQ5Pxe53pcmT/f1mXlSUuPkYYMGYKTJ0/i8uXLqtsw/BBRoDD8kANXY0Zmzy67sEdH28YBBWrMiD9ERrrvwpMCEQCMGFG2nXMAks6b/PyonUNA+bzJ63EVgJy7JM0gMjISNWrUcLmNFH5yc3Nx4cIFxMXFBaI0IgpHugybNolwne3lDVfTv51nfuk9eymQSwlInyk11f0ij87HVFvnx5fVlM06TV0PpaWlonLlygKA2L9/v9HlEFEQ4VR3HzD8aOdudWCl6d5auQs3zosFqoUbvYKAq8UH3e1f7Tz4UluwLVCoh9WrV4t58+aJn376yeV2V155pQAgvvjiiwBVRkRmwPDjA4Yfz7gLOGq3fNC6Xy2hIlBrCznvz5Obt7q655deAcjMwUcIIQYPHiwAiKVLl7rc7tZbbxUAxEsvvRSgyojIDBh+fMDw4zl397TypuVH/n533TuBbAHxptvKn91UenSjBYsVK1aIUaNGia1bt7rcbuLEiQKAmDp1aoAqIyIzYPjxAcOPZ7R26Xh7gdfSvRPoFhAtt/JQqlOP7XytJxQsWbJEABBDhgwxuhQiCiJc4ZkCRj7zS5phpDSN29vVgpWmfrtaY0htCrlenFd1lq/lo8Tfqyl7Wk8o4HR3IgoIXSKUSbDlxzf+aoXQ0r3j7xaQYBtjE2z1+KqoqEjTdr/88osAIGrUqOHniojITNjt5QOGH+/5e/yJq3Dj72MH2+yqYKvHVyUlJSIiIkJUrFhR/P333y63PXv2rADAv6dE5IDhxwcMP97xthVC63o90v6c757uy7G1CrZ1dYKtHj2cPn3aHmguXrzodvuEhAQBQOzduzcA1RGRGTD8+IDhx3O+tEJ4M6Xd1WBnT/arRSAGLJu5Hr2UlpaKU6dOiUOHDmnavlOnTgKAWL9+vX8LIyLT4IBnChi1wceA+/tUadlGfusItRuAentsLfw9YNns9ejFYrGgevXqqF69uqbtrVYrdu3axUHPROQ3DD+kyFXwkfgSgFwFHy30CECe3BQ0ELOsgq0eo3DGFxH5m1fh58KFCzh16hTq16/v8Pzvv/+OVq1a6VIYGUvPVgjnoCL/s9IxpGPLt/P0JqEUPH7++Wd8/PHHaNWqFQYPHux2eyn8ZGVl+bkyIgpXHoefjz76CJMnT0aNGjUghMCKFStw9dVXAwDuvfdepKen614kBZ7erRDyABQRYfuzWrhyPraWYEXBa+fOnUhLS0O/fv08Cj9s+SEif4nw9A0LFixAeno6fv75Z7zxxhsYMWIE3nvvPQCAEEL3AuUWLVqEzp07o0qVKqhduzYGDhyIAwcO+PWYpJ/Zs20L9pWW2r5qDU2eBDEKPk2aNMGoUaPQt29fTdvLw4+//00hovDkccvP5cuXUatWLQBAp06dsHXrVgwaNAh//PEHLBaL7gXKff/99xg/fjw6d+6M4uJiPProo7jhhhuwb98+VKpUya/HJt+F44rFzubMASIjtX3u+fNtrV5mD389e/ZEz549NW/fsGFDWCwWnDt3DidOnLD/exOWcnOBZcuAMWOAxESjqyEKGR63/NSuXRu//PKL/fuaNWti06ZN2L9/v8Pz/vDll19i+PDhaNWqFdq2bYuVK1fi8OHD2L17t1+PS76TD6AuKrJ9fewx2/PhJDJS2+eWzldkZGDqCiYxMTGoV68eAHZ9ITcXmDvX9pWIdKO55efs2bOoUqUK3nnnHURFOb4tOjoa77//PiZMmKB7ga7k5+cDAGrUqKH4elFREYqKiuzfFxQUBKQucqQ0c0yv6er+4q8WGi2fW8tMOzMpLS1FRIRn/8+yWq04evQoMjMzcdVVV/mpMiIKV5r/RerRoweOHTuGpKQk1K1bV3Gbbt266VaYO0IITJkyBd27d0fr1q0Vt1m0aBHi4+PtjwYNGgSsPrJxt1ZQsLYA+bOFxtXnDrXgAwA33XQTKlWqhDVr1mh+T0pKCgC2/BCRn2hdDXHkyJEiOTlZ7N+/3+H59PR0cdNNN+my4qInxo0bJxo2bCiys7NVt7l48aLIz8+3P7Kzs7nCcwCZfcVif99qItRuXqqmc+fOAoD45JNPNL9n9uzZAoAYPXq0Hyszgd27bb8Uu3cbXQmRoQxb4fm1117D3Llz0b17d2zYsAG1a9fGrFmzsHbtWtx6663+ymaKHnroIXzyySfYunUrkpKSVLeLiYlBTExMACsjObOvWOyqi0qPFhr5/hcssA0CD6UWH8nXX3+NU6dOoXbt2prfE5bT3XNzy4/tkZYOUVpCJDGRg6CJvOVpWlq4cKGIjY0VFSpUEP369RO7A/g/ktLSUjF+/HhRr149cfDgQY/fz3t7kTf83ULj6o724eq7774TAETjxo2NLiVw0tLKbnSn5ZGWZnTFRAFjWMtPbm4uFi1ahNdeew0tW7ZERkYG7rrrLnTo0MFvwczZ+PHj8d577+Hjjz9GlSpVcOzYMQBAfHw84uLiAlYHhRd/ttBw+r8yqeXnr7/+8mrAtCmNGQM4t6KnpwOjRgErVgDO/9ay1YfIe1pTUmxsrGjXrp347LPPhBBCfPnll6Jq1ariiSee0CWFaQFA8bFy5UpN72fLD/lC7xaacBjzc/78eTF37lzxwgsviOLiYs3vKy4uFlFRUQKAy3F9IY9jfoiEEAa2/KxcuRJ33XWX/fu+ffviu+++Q79+/fDXX3/h5Zdf1jOTKRJc7ZUMoncLjRmn/3vjxIkTSEtLQ4UKFTxaCiMyMhLJycn4888/kZmZ6XJsHxGRpzS3JcuDj6RDhw7YsWMHtmzZomdNREFF7wUazTr93xsVKlTA6NGjcd9993m8AnxYDnomooDw6q7ucikpKdi+fbsetRAFHb1baLTMEgulFqC6deti2bJlXr2X4YeI/MXn8AMA1atX12M3RH7h7WrN7lpoAM8Ditmn/wcSFzqEbVBzWhoHNxPpTJfwQxTMpNWaAdehQx52/NVC48lNSs3c4iMRQnh9w2O2/MAWesx+Z1uiIBQG80cp3GkZR+McdjxpoZk3L7xbaFx57rnnULFiRYwfP97j9zL8EJG/sOWHwoKnqzWHWwuNv5w5cwYXLlzw6r1S+Dly5AguXbqE6OhoPUsjojDG8ENhQykAheKNRIPJI488guHDhyM2Ntbj99apUwdxcXG4cOECsrOzccUVV/ihQiIKRww/FFbC5X5awaJy5cqoXLmyV++1WCxISUnB/v37kZmZyfBDRLrhmB8KO7Nnly1WGB3N4BPMOO6HiPyB4YfCjtJqzeQfK1euxIsvvojs7Gyv3s/wQ0T+wPBDYUXv1ZrJtaeeegr//ve/cejQIa/ez/BDRP7AMT8UNsLlflrBpH///mjdujUaNGjg1fu50KELubnAsmW2u8FzEUQijzD8UFjwx2rN5N5TTz3l0/vZ8uNCbi4wdy5w660MP0QeYvihkBdu99MKJVL4OX78OM6dO4dKlSoZXBERhQKGHwp5vJ+WeVWvXh3x8fHIz89HVlYWWrVqZXRJRBQCGH4o5HG1ZmPs27cPnTt3xhVXXIFffvnF6/1YrVbs3bsXmZmZDD9EpAuGHyLyi9OnT+P8+fM4f/68T/uRwk9WVpY+hZlRbq7tIZee7vhVLjHR/TggDpimMMbwQ0R+0alTJxw6dAhFRUU+7YeDnmELKXPnKr82alT559LS3Dd5csA0hTGGHyLyi5iYGDRq1MinfZSUlKC4uBgAsHPnTpSUlCAyMlKP8sxlzBhbSJFbv952j5ZZs4DbbnN8jWGGyCWGHyIKSuvWrcPEiRNx5MgRAMC2bduQkpKC559/HoMGDTK4ugBz7sbKzQX27bP9uXlzoEMHY+oiMimu8ExEfrF161a8+OKL+O9//+vxe9etW4chQ4bYg4/k6NGjGDJkCNatW6dXmeaUmwsE6hzk5tq60JzHHBGZGFt+iMgv1q1bh+effx7Tpk3D1Vdfrfl9JSUlmDhxIoQQ5V4TQsBisWDSpEkYMGBAeHaBecOXAdMcG0QhiOGHiPyiXbt2uP3229G+fXuP3rdt27ZyLT5yQghkZ2dj27ZtSE1N9bFKk3AOL/v3l/05I6N8gHHuJvPHgGkiE2P4ISK/GD58OIYPH+7x+3I1dq9o3S4kuAovCxbYHnLO4UVpwHR6ui34zJple/+qVUCLFrbX2MJDIY7hh4iCSqLGC6/W7UwvNxfo2tUWTiQZGWWBZ9Ys26BnSUIC0KaN4z5crfsjvbdFi/IDp6W1gIhCDMMPEQWVHj16ICkpCUePHlUc92OxWJCUlIQePXoYUJ0BXLX6AMqtPn37encspe615cttf/Z2MUWiIMTwQ0R+0axZM5w7dw5fffWVR7eliIyMxPPPP48hQ4bAYrEoBqAlS5aEz2BnV11WgGN3FeBbGOHYIAoTDD9E5BdHjhzB+fPnERcX5/F7Bw0ahI8++shhnR8AiIqKwurVq8NrnR93rStK3VVq5C070qDpjAzb1/T08t1rp08DDz1k+/OKFeWPw1YfMimLUPpvVYgqKCiw3yG6atWqRpdDFNL+/PNPnD59GldeeSWio6O92kdJSQm2bduGvXv3YvLkyYiLi8O5c+dgsVh0rtZk0tOBjh1tf969W3v4mTPHdReas9Gjy7q9PDkOkc70vn6z5YeI/MLXW1sAti6w1NRUXHPNNXj44Ydx4cIF/P3336hbt64OFYaA0aM9a31x1YW2eHH5falNrwc43odMjeGHiIJedHQ06tevj+zsbGRmZjL8SDy9I7urwPK//wFTp6q/9557HL8PpvE+vEM9eYjhh4h0d+TIEaxfvx7JyckYMGCALvu0Wq328NO1a1dd9mlaiYm28KHnhX7QIFt4kPPXwGq9cRVq8hDDDxHp7rfffsO///1vtGvXTtfws3XrVmRmZuqyP1NLTPS+1UWtlaRWLfUxPaNHA717ex4s2CJDQYo3NiUi3dWsWRO33347+vTpo9s+rVYrADD8+EpqJfFkhWxvw4s3xyIKALb8EJHuOnfujA8//FDXfTL8+Ik/utCIghzDDxGZghR+srKyyr/I7hXt8vIcv/rShWYEX+5QT/QPhh8iMgUp/Bw+fBglJSWOKzxzwKt2J044fnUlGFuFuAo16YDhh4h0N2HCBKxfvx5paWkYPXq0LvusV68eoqOjcenSJRw5cgQNGzbUZb8hTamVRFrROSOjfEtJRASwYUNZC5onrUKBapFxtVYRV6EmjRh+iEh3OTk5yMnJQXFxsW77jIiIQMOGDfG///0PmZmZDD9auGolWbCg/E1RpRWdvWlBC1SLjKvQ1KEDV6EmTRh+iEh3S5cuxaxZs1C/fn1d95uSkmIPP6mpqbruOyQptZKsX28LPbNmAbfd5vhaXl7Z7Sz0OBZbZChIMfwQke7q1auHevXq6b5f+4yvHTuAtm3LXuCAV+OxRYZMhOGHiEzDHn5eew147bXyG3DAqyNvur08wVl2ZFIMP0Sku1dffRWVKlXCoEGDUKlSJd32aw8/bdsCb7xR9gK7V5S56vaaOBHo3NnxNWnAstYWtGCZZReMs9IoqDH8EJGuSkpK8OCDDwIAbrzxRv+En7w85W4Udq84Ugos0t3ZDx8Gnn9e+X1ma0Ez21pFZDiGHyLSVVFREYYMGYLTp0+jWrVquu5bCj85OTkoKipCTEyMrvvXVbB2CSUk2L4OGWIb9Cyn1IKWlwesWwcMHOj5sdgiQ0GK4YeIdFWxYkWsWbPGL/tOSEhApUqVcO7cOfz1119o2rSpX46jC3ddQkaFo1q1bF+bN1dvJZO3oKWn22aAOd/xXQu2yFCQYvghItOwWCywWq347bffkJmZGdzhx51gGS+jVV5e+bFAnGVHJsXwQ0SmIg8/duxe0c7bc7VunfoaQGYbI0Rhj+GHiHS1fv16PPTQQ+jTpw/efPNN3fefkpICAOXDDy+02nh7rgYNKt/1xVl2ZFIMP0Skq7y8PBw9ehSnT5/2y/7tM77k4cdo3tzXKhhFRNjW+pF3cUlfs7PLxgtJGjSwfeUsOzIZ04Wfl19+GYsXL0Zubi5atWqFJUuWoEePHkaXRUT/uP3229GxY0fExcX5Zf9BGX48va/V6NFla+x4Ml7G34OkN2ywdW0pdW+pfQ4iE7IIIYTRRWi1evVq3HvvvXj55ZfRrVs3LFu2DK+99hr27duH5ORkt+8vKChAfHw88vPzUbVqVdXtSkpKcPnyZT1LJ1JUoUIFREZGGl2Gqfz8889o164datasiRMnThhThHMIUWv5kbqEdu707J5ZauNl0tOBjh2B3bv909Li7nM4HzMvD7jxRv/VQ/QPrddvrUzV8vPss89i5MiReOCBBwAAS5YswVdffYVXXnkFixYt8nn/QggcO3YMZ86c8XlfRFpVq1YNdevWhcViMboUU5Bafk6ePImzZ8+iSpUqgS/CeaaWq5lNDRrYFhb88kvHbqNgHC+j5f5c8uCXl2d7TfpKZBKmCT+XLl3C7t27MX36dIfnb7jhBuzYsUPxPUVFRSgqKrJ/X1BQ4PIYUvCpXbs2KlasyIsR+ZUQAufPn8fx48cBAInBOg7EQ19++SVOnDiBbt262YOKnqpWrYoaNWrg1KlTyMzMRJs2bXQ/hq5OnACefRa44Qbgk0/Kd1mZbbyMPPgRmZRpws+JEydQUlKCOnXqODxfp04dHDt2TPE9ixYtwly1fngnJSUl9uBTs2ZNn+sl0kIaF3P8+HHUrl07JLrAnnnmGXzzzTd4++23/RJ+AFvrz6lTp5CVlRX84Udy4oTyuj7Lltm6uKTnArX4oR7HkVqynAdCEwW5CKML8JRza4wQQrWFZsaMGcjPz7c/srOzVfcrjfGpWLGifsUSaSD9zoXKOLNOnTqhT58+fgs+QJAOevbW8uWO42yklpVffrF1jTk/AOXnncfquCMdx9P3ebL/OXP8t38iH5im5SchIQGRkZHlWnmOHz9erjVIEhMT4/G9f9jVRYEWar9zeoy/cyeg4cebaezSIy2t7F5angqWRQW9XRTRbCtYU1gxTfiJjo5Gx44dsWnTJtx222325zdt2oQBAwYYWBkRBZpX4cfbbh5vprFLiwHeemtZQMrIsH2V7qqelwcMHQq8+67yvoNhUUEp+N16q60lavNm28w1AHj7bUAabL5+ve1zJSSUdYFxEDQFMdOEHwCYMmUK7r33XnTq1Aldu3bF8uXLcfjwYYwdO9bo0ogogBRXeXbH25aIMWPKD+5VCyHLlqmvk7Ngge3rPfeUPTdokO2rUlCoVUvbjUf9yVXwe/75sj9Ln02OawBREDNV+Lnzzjtx8uRJzJs3D7m5uWjdujW++OILNGzY0NC65swBIiOB2bPdbzt/PlBSwpX4KTQVFBSgRYsWqF69Ovbs2YMKFSr45Tjylh9X4/50oWX6t2TOHFtYkkKQO+vW2b4+/rhtQDRQ1kK0eXP57b1tTfG2604e/PLybDV+9x3w+uvAyJG2bRYsAGbNst0l3rnlx5O1jYgCyFThBwDGjRuHcePGGV2Gg8hI4LHHbH92FYDmz7dtN29eYOoiCrRTp04hJycHp06d8lvwAcpafgoLC3Hy5EkkeDuuRm9SaJBCEGDrDrrnHuDaa4GtW5Xft22b7SE3dWr57bxtTfG0604aP6QW/F5/HejVC2jRwhZ+uncvP+NLmmDCO75TEDJd+AlGUuBxFYDkwUdLC1Gw4u1FyJXExETs3r0b586d8+txYmNjkZiYiNzcXGRmZgZP+JEoXdxvv90WflatsoUGwNa6M3VqWcsJYGv5WbAAWLwY6N3bcR/etqZ40nUn1e+JYBmcTaQRw49OXAWgUAk+q1evxqRJkxxuL3LTTTdpvr0Ihb6YmBh0CNCCfVar1R5+Okv3yZJ4283jD85dVVu32kKN/Hi33VYWQNLTbeGnd+/yoSQ317uZV5503Tkfz/k8St1y0lcAuPJKW6iTur3y8oCXXgI+/TS4VrAmkogwkp+fLwCI/Pz8cq9duHBB7Nu3T1y4cMGnY8ybJwRg+6r0vb/Ur19fvPTSSw7Pbd++XcTFxYmsrCxdjnHVVVeJsWPHOjzXvHlzMX36dF32H670+t0LN0OHDhUAxJNPPln+xbQ02188rY+0NM8OnpNje09OjvttV62yHWP+/LLj7d5te233bsfvlZ7z5FieUDq2M2/Po7Rvd/sn0sjV9dsbbPnRmbwFaMEC4NKlwLT4dOnSBTulKaiwLf44adIkTJo0qdyA8IULF2LhwoUu97dx40aH7ixvbi9C4Wffvn1IT09Hs2bNyrfG6MzldHd/d/NI43q0kLrkqlXz7BgSI9fLkc5jXp6ta2vQIODXX21ddYsX21p8pOdr1WKLDpkGw48fzJ5dFnyiowPT1dWlSxe8+eab9u/feecdHD58GDNmzCi37dixY3HHHXe43F/9+vUdvvfm9iIUfr744gtMnToVQ4cOxapVq/x6LJfhx9tuHn+QBgI3aWIbsGymGVDSeUxPt9U9ZkzZOCSpW65vX/f7CdQtO4g0Mt3tLcxg/vyy4HPpku17f+vSpQv279+PwsJCnD9/HjNnzsSCBQsU73hdo0YNNG7c2OVDuueUM09uL0Lhp379+ujTpw/atm3r92OZ7hYXtWqVX7RQafVkb1dU9pR0nIgI/9+Gwt+30iDyEFt+dOY8uFn6HvBvC1CnTp0QGRmJ9PR0fPPNN6hZsyZGjBihuK033V7e3F6Ews/dd9+Nu+++OyDHkqa7Z2VlobS0FBERQfB/OfkAYWldHGlg8Pr1ZdvJV0S+9day98inyvubdJz0dN+61VwNLu/fv+w8ALbPzGnuFAQYfnSkNKtLyzR4PcTGxqJt27ZYt24dli9fjk8//VT1YuBNtxdvL0LBpkGDBoiMjMSlS5eQm5tb7nfWEK7W05Gvgqy0InJamq1lKFhmqWnl6jN/+qntIbnnHk5zp6DA8KMTV9PZAxWAunTpghdeeAH9+vXDddddp7pdjRo1UKNGDY/3z9uLUDCJiopCgwYNkJWVhczMTPfhJxDdSUorIkvr9syaZXtebUXkxETvFyPUk7tlAho0sI1dysuzfT9woPvB5dJCj6tWlV+7iMgADD860LKOTyACULt27RAVFYXFixfrv3ME7+1FKHgMGDAAhw4dwiuvvBKQxS+tVqs9/HTv3t31xoHoTlJqiZHW7ZFaTKU/Kw269vcsNS20BjBp4LZzAMvNBaQxg9LXCxfKvjqHq2BovaKww/Cjg5ISbdPZpddLSvxTx7vvvotx48ahWbNm/jkAgvP2IhQ8MjIycPDgQZSWlgbkeFarFd99913gBz37a/aSv2epaVn8sWtXWwsNUNYy5UkAk4cn+U1cAa72TEGD4UcHnvy91bvFp7S0FHl5eXj99ddx4MABrJcPqiQKsDVr1iAvLy8gs72AshlfWVlZATmenZFr7/jC1241LQFszBjbtH6pm6tFi8C3XhG5wfBjclu3bkXv3r3RvHlzrFu3DvHx8UaXRGGsTZs2AT2e6aa7Gy0Q3WqJiWX3LmvRwnGfgV5jiUgFw4/JpaamBqyLgSjYmCL8OA+09uega3fdcYFa/DFQaxUReYnhh4h0cfbsWWzYsAEJCQm46aabAnJMKfxkZ2fj8uXLqFChQkCO6xHngdaejm/xJEgES3dcoNYqIvISww8R6SIrKwv33XcfatWqhePHjwfkmHXq1EFMTAyKioqQnZ2NRo0a6X8QV4OEN2+2tbRI97YC9J+9FAxBgi05FGIYfohIF1FRUejTpw+qVq0asGNGREQgJSUFBw4cQGZmpn/Cj6tBwlOn2r7K79cVirOXfA1gDE8UZBh+iEgXLVq0wKZNmwJ+XKvVag8/fuFqkPCsWbZ1e6RZTYD5LvCBCCbB0HpFJMPwQ0Sm5vdBz666sZo3t311ntUUCHp1xzGYUBhi+CEiUzPFjC9/YHcckdcYfohIF88++yzeeOMNDB8+HP/5z38CdtyAhR95S8v+/bav0h3bjbjpaKh3xxH5EcMPEekiKysLv//+O06dOhXQ4/p9lWdp7ZyzZ4Fnn3V8Tbo7uxG3bQim7jh/3e6DyE8YfvyF/xhQmJk8eTJuvfVWNGjQIKDHlcLPsWPHcOHCBcRJN9PUi7R2zpdfAkOH2p7LywPWrbPdxmHqVN62IVjWFyLSiOHHX/iPAYUZq9VqDyKBVL16dVStWhUFBQXIyspCC6mbR2+1ajkGnL59y7q7jLxtQ7B1xxGZQITRBZB5bN26Ff3790e9evVgsViwYcMGo0sigsViQUpKCoAwHPQM2FqYO3a0PaS7qMu746TXpMeyZcbVShQk2PJDmp07dw5t27bF/fffj8GDBxtdDgWZTz/9FBEREejWrRuqVasW0GNbrVb88ssv4Rl+5AOf2R1HpAnDT4hISkrCzJkzMW7cOPtzO3bsQJ8+fbB//340bNjQ52PcdNNNAbtnE5nPmDFjkJubi/T0dLRv3z6gx9ZtxpertXOUupDy8nw7nh6cu7H83R3n6TliNxsFIYYfPQTBPwZdunTBzp077d8LITBp0iRMmjSpXPBZuHAhFi5c6HJ/GzduRI8ePXStkUJbx44dkZubi4SEhIAfW7fw42rtHKUZXVOmhN9tGzw9R1xfiIIQw48eguAfgy5duuDNN9+0f//OO+/g8OHDmDFjRrltx44dizvuuMPl/urXr69rfRT6Pv30U8OOrVv4cbV2jloXUjgFH8C7c0QUZBh+9BAE/xh06dIF06ZNQ2FhISIiIjBz5kwsWLAAVapUKbdtjRo1UKNGDd1rIDKKbuHHVZgxckZXILlbpoPniEIAw48eguAfg06dOiEyMhLp6en45ptvULNmTYwYMUJxW3Z7UaiRZnudOXMGZ86cCfiA66Djy81KuUwHhQGGnxARGxuLtm3bYt26dVi+fLl95o0SdnuR3nbv3o1hw4ahdevW+OCDDwJ+/MqVK6NWrVrIy8tDVlYW2rVrF/AaggpvVkrkEsNPCOnSpQteeOEF9OvXD9ddd53qdt52exUWFuKPP/6wf5+ZmYm9e/eiRo0aSE5O9qpmCg3Hjh3D77//jtjYWMNqsFqtyMvLQ2ZmJsMPEbnE8BNC2rVrh6ioKCxevNgv+9+1axd69epl/37KlCkAgGHDhjkMtqbw06VLF3zzzTeoUKGCYTWkpKTgp59+0n+tH1+6kMIFzxGZDMOPvxjwj8G7776LcePGoVmzZn7Zf2pqKoQQftk3mVvNmjVdtjYGgt/u7h7KXUh6LdMRyueIQhLDj78E6B+D0tJS5OXl4fXXX8eBAwewfv16vx+TKBj5LfyEsiBYpoPICAw/Jrd161b07t0bzZs3x7p16xAfH290SRSGdu3ahWPHjqFVq1aG3NwUYPjxShAs00FkBIYfk0tNTUVpaanRZVCYe+mll/Dmm29i4cKFigtrBoIUfrKysiCEgMViMaQOUwmCZTrKcbfOEJEOeFd3IvJZcnIyOnXqZF9vx6gaLBYLzp8/j+PHjxtWB/lIWmfIeSwSkY4YfojIZ3PnzsXOnTtx9913G1ZDTEyMfX0qdn0RkSsMP0QUMuRdX0REahh+iChkcNCzDrhmD4UBDngmIp9dffXViI6OxocffohEAy+aDD86COSaPXqtM0TkIYYfIvLJ5cuX8dNPPwEAoqOjDa1FGnDN8GMSXGeIDMLwQ0Q+sVgs+Oabb3D69GnD76bOlh+T4TpDZBCGHyLySVRUlOG3tpBI4efw4cMoKSlBZGSkwRWRS8G4zhCFBQ54JqKQUb9+fVSoUAGXL1/G0aNHjS6HiIIUw48flJSUYMuWLXj//fexZcsWlJSUGF0Skd/k5OTgs88+w969e40uBZGRkUhOTgbAri8iUsfwo7N169YhJSUFvXr1wr/+9S/06tULKSkpWLdundGl+eyJJ55Aq1atULFiRTRt2hTvvfee0SVRENi6dSv69++PyZMnG10KAI77ISL3GH50tG7dOgwZMgRHjhxxeP7o0aMYMmSI6QPQtm3b8Nxzz+G3337DPffcg/vuuw9//vmn0WWRwSpVqoROnTqhZcuWRpcCgAsdmh7XGaIAYPjRSUlJCSZOnAghRLnXpOcmTZrkty6wpKQkvPzyyw7P7dixAxUrVsRff/2lyzE+//xz3HDDDWjUqBEmTJiAkpIS5OTk6LJvMq/+/ftj586deOmll4wuBQBbfkxPWmeI4Yf8yBThJysrCyNHjoTVakVcXByuuOIKpKWl4dKlS0aXZrdt27ZyLT5yQghkZ2dj27Ztfjl+ly5dsHPnTofjTZo0CZMmTULDhg0dtl24cCEqV67s8uGqTiEEHn74YbRu3RpXXXWVXz4PkbcYfojIHVNMdc/IyEBpaSmWLVuGxo0b47fffsOoUaNw7tw5PP3000aXBwDI1XgHYq3beapLly5488037d+/8847OHz4MGbMmFFu27Fjx+KOO+5wuT/pBpFKHnjgAezYsQObN282fFE7Imdc6JCI3DFF+Lnxxhtx44032r9v1KgRDhw4gFdeeSVowo/WJf39tfR/ly5dMG3aNBQWFiIiIgIzZ87EggULUKVKlXLb1qhRAzVq1PDqOL/88gveeOMNZGRkuAxIFD6mTJmCn376CdOnT0e/fv2MLsfe8nP06FEUFRUhJibG4IqIKNiYottLSX5+vtsLeFFREQoKChwe/tKjRw8kJSXBYrEovm6xWNCgQQP06NHDL8fv1KkTIiMjkZ6ejieeeAI1a9bEiBEjFLf1pdtL+t90s2bN/PI5yHx+/vlnbN++Hfn5+UaXAgCoXbs2KlasCCEEDh8+bHQ5RBSETNHy4+zQoUN48cUX8cwzz7jcbtGiRZirdt8YnUVGRuL555/HkCFDYLFYHAY+S4FoyZIlfltxNjY2Fm3btsW6deuwfPlyfPrpp4iIUM62vnR79ezZ02FsEdGiRYtw+PDhoBn/ZbFYkJKSgn379iEzMxNNmjQxuiQiCjbCQGlpaQKAy8fOnTsd3nP06FHRuHFjMXLkSLf7v3jxosjPz7c/srOzBQCRn59fbtsLFy6Iffv2iQsXLvj0mdauXSuSkpIcPkODBg3E2rVrfdqvFhMmTBAWi0X079/fb8dYt26daNasmd/2H470+t2jMrfccosAIF599VWjSyEiHeTn56tev71haMvPhAkTcNddd7ncRhq8CNhWku3Vqxe6du2K5cuXu91/TExMwPv7Bw0ahAEDBmDbtm3Izc1FYmIievToEZB7DLVr1w5RUVFYvHix346Rn5+PAwcO+G3/RHrgjC8icsXQ8JOQkICEhARN2x49ehS9evVCx44dsXLlStUunWAQGRmJ1NTUgB/33Xffxbhx4/w6Hmf48OEYPny43/ZP5lJaWoovvvgC1apVQ5cuXRAVFRw96Qw/RORKcPxL5UZOTg5SU1ORnJyMp59+Gnl5efbX6tata2BlxistLUVeXh5ef/11HDhwAOvXrze6JAoj+fn56N+/PwDbBINgwVWeicgVU4Sfr7/+Gn/88Qf++OMPJCUlObwmFFZUDidbt25F79690bx5c6xbtw7x8fFGl0RhpKioCB07dsTFixeDas0ntvwQkSsWEUbpoaCgAPHx8cjPz0fVqlUdXrt48SIyMzNhtVoRGxtrUIUUjvi7p78zZ86gevXqAICzZ8+icuXKBldERL5wdf32RvAOnCEi8lK1atVQrVo1AOz6IqLyGH6IKCSx64uI1DD8EJHX1q5dix49emDBggVGl1IOww8RqTHFgGciCk6HDh3C//3f/6FRo0ZGl1IOww8RqWH4ISKvDRo0CI0aNQrKm9wy/BCRGoYfIvJa48aN0bhxY6PLUMTwQ0RqOOaHiEKSfKHDMFrRg4g0YPghIq/997//xfbt23H69GmjSylHui9gQUFBUNZHRMZh+AkzFosFGzZsMOTYWVlZsFgs2Lt3ryHHd2fLli2wWCw4c+aM0aWYxoQJE9C9e3ds27bN6FLKiYuLQ506dQCw64uIHDH8hIDhw4fDYrHAYrGgQoUKqFOnDq6//nq88cYbKC0tddg2NzcXN910k0GVUqhp0KABrrjiCtSuXdvoUhRx3A8RKWH4CRE33ngjcnNzkZWVhY0bN6JXr16YOHEi+vXrh+LiYvt2devWRUxMjO7HLykpKRe0wtGlS5eMLiGg1q1bhz/++ANdunQxuhRFDD9EpIThxwUhBM6dOxfwhzeDM2NiYlC3bl3Ur18fHTp0wMyZM/Hxxx9j48aNePPNN+3bybu9Ll26hAkTJiAxMRGxsbFISUnBokWL7NueOXMGo0ePRp06dRAbG4vWrVvjs88+AwC8+eabqFatGj777DO0bNkSMTEx+OuvvzTVmpGRgWuuuQaxsbFo1aoVtmzZ4vD6999/j6uuugoxMTFITEzE9OnTHQJcSkoKlixZ4vCedu3aYc6cOQ6f87XXXsNtt92GihUrokmTJvjkk08c3vPFF1+gadOmiIuLQ69evcrdBuHkyZO4++67kZSUhIoVK+LKK6/E+++/77BNamoqJkyYgClTpiAhIQHXX389RowYgX79+jlsV1xcjLp16+KNN97QdI5IHww/RKSEU91dOH/+vCE3RCwsLESlSpV83k/v3r3Rtm1brFu3Dg888EC511944QV88skn+PDDD5GcnIzs7GxkZ2cDAEpLS3HTTTfh7NmzWLVqFa644grs27cPkZGR9vefP38eixYtwmuvvYaaNWtq7vqYOnUqlixZgpYtW+LZZ5/FrbfeiszMTNSsWRNHjx7FzTffjOHDh+Ptt99GRkYGRo0ahdjYWIdwo8XcuXPx1FNPYfHixXjxxRcxdOhQ/PXXX6hRoways7MxaNAgjB07Fg8++CB27dqFhx9+2OH9Fy9eRMeOHTFt2jRUrVoVn3/+Oe699140atQIV199tX27t956Cw8++CC2b98OIQROnTqFa6+9Frm5uUhMTARgC1qFhYW44447PPoM5BuGHyJSJMJIfn6+ACDy8/PLvXbhwgWxb98+ceHCBftzhYWFAkDAH4WFhR59rmHDhokBAwYovnbnnXeKFi1a2L8HINavXy+EEOKhhx4SvXv3FqWlpeXe99VXX4mIiAhx4MABxf2uXLlSABB79+7VXGdmZqYAIJ544gn7c5cvXxZJSUniySefFEIIMXPmTNGsWTOHml566SVRuXJlUVJSIoQQomHDhuK5555z2Hfbtm1FWlqaw+ecNWuW/fvCwkJhsVjExo0bhRBCzJgxQ7Ro0cLhONOmTRMAxOnTp1U/w8033ywefvhh+/c9e/YU7dq1K7ddy5Yt7Z9JCCEGDhwohg8frrhPpd89M8jJyRHdunUTgwcPNroUVd98840AIJo1a2Z0KUTkA1fXb2+w5ceFihUrorCw0JDj6kUIAYvFovja8OHDcf3116NZs2a48cYb0a9fP9xwww0AgL179yIpKQlNmzZV3Xd0dDTatGnjcU1du3a1/zkqKgqdOnXC/v37AQD79+9H165dHWru1q0bCgsLceTIESQnJ2s+jry2SpUqoUqVKjh+/Lj9OF26dHE4jrwuwDaO6YknnsDq1atx9OhRFBUVoaioqFyrXKdOncod+4EHHsDy5cvxyCOP4Pjx4/j888/x7bffaq7dDI4fP47t27fbZ1QFI/laP6WlpYiIYE8/EbHbyyWLxaJL95OR9u/fb78AOOvQoQMyMzOxceNGfPPNN7jjjjvQp08ffPTRR4iLi3O777i4ONVg5SlpP0phTfwzBkp6PiIioty4qMuXL5fbZ4UKFcodQxqU7fx+Jc888wyee+45LFmyBFdeeSUqVaqESZMmlRvUrPQ7ct9992H69On44Ycf8MMPPyAlJQU9evRwe0wzSU5Oxpo1a4J6oHuDBg0QERGBoqIiHDt2DPXq1TO6JCIKAvxvUAjbvHkzfv31VwwePFh1m6pVq+LOO+/EihUrsHr1aqxduxanTp1CmzZtcOTIERw8eFD3un788Uf7n4uLi7F79240b94cANCyZUvs2LHDIZzs2LEDVapUsd8/qlatWsjNzbW/XlBQ4PGYjpYtWzrU4VwXAGzbtg0DBgzAPffcg7Zt26JRo0b43//+p2n/NWvWxMCBA7Fy5UqsXLkS999/v0f1mUH16tUxZMiQoB7HVKFCBTRo0AAAyg1oJ6LwxfATIqT/2R49ehTp6elYuHAhBgwYgH79+uG+++5TfM9zzz2HDz74ABkZGTh48CDWrFmDunXrolq1aujZsyeuvfZaDB48GJs2bbK3EH355ZeqNcyYMUP1WHIvvfQS1q9fj4yMDIwfPx6nT5/GiBEjAADjxo1DdnY2HnroIWRkZODjjz9GWloapkyZYu+y6N27N9555x1s27YNv/32G4YNG+YwEFuLsWPH4tChQ5gyZQoOHDiA9957z2FWHGC7b9WmTZuwY8cO7N+/H2PGjMGxY8c0H+OBBx7AW2+9hf3792PYsGEe1Uf64aBnInLGbq8Q8eWXXyIxMRFRUVGoXr062rZtixdeeAHDhg1THedQuXJlPPnkk/jf//6HyMhIdO7cGV988YV9+7Vr1+I///kP7r77bpw7dw6NGzfGE088oVpDbm4uDh8+7LbWJ554Ak8++ST27NmDK664Ah9//DESEhIAAPXr18cXX3yBqVOnom3btqhRowZGjhyJWbNm2d8/Y8YM/Pnnn+jXrx/i4+Mxf/58jy9sycnJWLt2LSZPnoyXX34ZV111FRYuXGgPYQAwe/ZsZGZmom/fvqhYsSJGjx6NgQMHIj8/X9Mx+vTpg8TERLRq1Soku1syMzORk5ODhg0bIikpyehyVEm3uWD4ISKJRWgZ/BAiCgoKEB8fj/z8fFStWtXhtYsXLyIzMxNWqxWxsbEGVUih5Pz586hXrx7eeOMNDBo0SHU7s/7uzZ49GwsWLMD48eOxdOlSo8tRNW/ePKSlpWHEiBF4/fXXjS6HiLzg6vrtDbb8EOmstLQUx44dwzPPPIP4+HjceuutRpfkF1WqVMEVV1wR1K0+ALu9iKg8hh8inR0+fBhWqxVJSUl48803ERUVmn/NHnnkETzyyCNGl+EWww8ROQvNf5WJDJSSkuLVLUrIP6Twk52djeLi4pANo0SkHWd7EVFIS0xMRExMDEpKSuy3byGi8MbwQ0ReGTp0KPr3729fnTtYRUREoGHDhgDY9UVENgw/ROSVb7/9Fp999hmKioqMLsUt+W0uiIjY+U1EXnn11Vdx4sQJ+zo6wYyDnolIjuGHiLwycOBAo0vQjAsdEpEcu72IKOSx5YeI5Bh+yC8sFgs2bNhgdBmKsrKyYLFYsHfvXqNLMa2zZ89i+/btOHDggNGlaMLwQ0RyDD8h4tixY5g4cSIaN26M2NhY1KlTB927d8err76K8+fPG10ehZiff/4Z3bt3R79+/YwuRRMp/OTm5uLChQsGV0NERuOYnxDw559/olu3bqhWrRoWLlyIK6+8EsXFxTh48CDeeOMN1KtXL2RvsWCUS5cuITo62ugyDCOEwBVXXGEPFcGuZs2aqFy5MgoLC/HXX3+hefPmRpdERAZiy48G586dw7lz5xxW7b106RLOnTtXbpqvtG1paan9ucuXL+PcuXO4ePGi2229MW7cOERFRWHXrl2444470KJFC1x55ZUYPHgwPv/8c/Tv39++bX5+PkaPHo3atWujatWq6N27N37++Wf763PmzEG7du3wzjvvICUlBfHx8bjrrrtw9uxZj+vKzc3FTTfdhLi4OFitVqxZs8bh9V9//RW9e/dGXFwcatasidGjR6OwsND+empqKiZNmuTwnoEDB2L48OH271NSUux3Y69SpQqSk5OxfPlyh/f89NNPaN++PWJjY9GpUyfs2bPH4fWSkhKMHDkSVqsVcXFxaNasGZ5//nmHbYYPH46BAwdi0aJFqFevHpo2bYp58+bhyiuvLPe5O3bsiMcee8yTU2U6PXr0wB9//IGvv/7a6FI0sVgs7PoiIjuGHw0qV66MypUr48SJE/bnFi9ejMqVK2PChAkO29auXRuVK1fG4cOH7c+99NJLqFy5MkaOHOmwbUpKCipXruzTInEnT57E119/jfHjx6NSpUqK21gsFgC2/63fcsstOHbsGL744gvs3r0bHTp0wHXXXYdTp07Ztz906BA2bNiAzz77DJ999hm+//57PPHEEx7XNnv2bAwePBg///wz7rnnHtx99932z3r+/HnceOONqF69Onbu3Ik1a9bgm2++KXc+tXjmmWfsoWbcuHF48MEHkZGRAcAWMPv164dmzZph9+7dmDNnDv7zn/84vL+0tBRJSUn48MMPsW/fPjz22GOYOXMmPvzwQ4ftvv32W+zfvx+bNm3CZ599hhEjRmDfvn3YuXOnfZtffvkFe/bscQhoFBwYfohIwvBjcn/88QeEEGjWrJnD8wkJCfbQNm3aNADAd999h19//RVr1qxBp06d0KRJEzz99NOoVq0aPvroI/t7S0tL8eabb6J169bo0aMH7r33Xnz77bce13b77bfjgQceQNOmTTF//nx06tQJL774IgDg3XffxYULF/D222+jdevW6N27N5YuXYp33nkHf//9t0fHufnmmzFu3Dg0btwY06ZNQ0JCArZs2WI/TklJCd544w20atUK/fr1w9SpUx3eX6FCBcydOxedO3eG1WrF0KFDMXz48HLhp1KlSnjttdfQqlUrtG7dGklJSejbty9Wrlxp32blypXo2bMnGjVq5PH5Iv9i+CEiCcOPBoWFhSgsLERCQoL9ualTp6KwsBBLly512Pb48eMoLCxEcnKy/bnx48ejsLAQr7/+usO2WVlZKCwsRIsWLXyuUWrdkfz000/Yu3cvWrVqZe+a2717NwoLC+3jH6RHZmYmDh06ZH9vSkoKqlSpYv8+MTERx48f97imrl27lvteavnZv38/2rZt69Ba1a1bN5SWlno8g6hNmzb2P1ssFtStW9der3ScihUrqtYF2Bbs69SpE2rVqoXKlStjxYoVDq13AHDllVeWG+czatQovP/++7h48SIuX76Md999FyNGjPCofjN68cUXceutt5brygxmXOWZiCQc8KyBUndSdHS04oBXpW0rVKiAChUqaNrWU40bN4bFYrF380ikloe4uDj7c6WlpUhMTLS3ishVq1bNoV45i8Xi87gk+b4AWxecc2Bz3iYiIqLc3dEvX75cbntX9Wq5u/qHH36IyZMn45lnnkHXrl1RpUoVLF68GP/9738dtlP6efXv3x8xMTFYv349YmJiUFRUhMGDB7s9ptnt2rULn376Kbp37250KZpxoUMikrDlx+Rq1qyJ66+/HkuXLsW5c+dcbtuhQwccO3YMUVFRaNy4scND3qqllx9//LHc99Ism5YtW2Lv3r0ONW/fvh0RERFo2rQpAKBWrVrIzc21v15SUoLffvvNoxpatmyJn3/+2WF6s3Nd27ZtwzXXXINx48ahffv2aNy4sUNLmCtRUVEYNmwYVq5ciZUrV+Kuu+5yaGUKVaNHj8by5cvRt29fo0vRjN1eRCRh+AkBL7/8MoqLi9GpUyesXr0a+/fvx4EDB7Bq1SpkZGQgMjISANCnTx907doVAwcOxFdffYWsrCzs2LEDs2bNwq5duzQfb+nSpbjuuuvcbrdmzRq88cYbOHjwINLS0vDTTz/ZBzQPHToUsbGxGDZsGH777Td89913eOihh3DvvfeiTp06AIDevXvj888/x+eff46MjAyMGzcOZ86c8ejc/Otf/0JERARGjhyJffv24YsvvsDTTz/tsE3jxo2xa9cufPXVVzh48CBmz57tMIjZnQceeACbN2/Gxo0bw6LLC7B1UY4aNQpt27Y1uhTNpPBz6tQpFBQUGFwNERmJ4ScEXHHFFdizZw/69OmDGTNmoG3btvbBxf/5z38wf/58ALbuoC+++ALXXnstRowYgaZNm+Kuu+5CVlaWPXBoceLECU0tI3PnzsUHH3yANm3a4K233sK7776Lli1bAgAqVqyIr776CqdOnULnzp0xZMgQXHfddQ5jqEaMGIFhw4bhvvvuQ8+ePWG1WtGrVy+Pzk3lypXx6aefYt++fWjfvj0effRRPPnkkw7bjB07FoMGDcKdd96Jq6++GidPnsS4ceM0H6NJkya45ppr0KxZM1x99dUe1UeBU6VKFdSsWRMAW3+Iwp1FaBkUESIKCgoQHx+P/Px8VK1a1eG1ixcvIjMzE1arFbGxsQZVSGYkhEDz5s0xZswYTJkyxeP3m/F3b+fOnahUqRIaN25sqsUeO3fujF27dmH9+vWmujErUbhzdf32Blt+iHxw/PhxPPvsszh69Cjuv/9+o8sJmJ49e6JVq1bIzs42uhSPcNwPEQGc7UXkkzp16iAhIQHLly9H9erVjS4nIIqLi1GvXj2cPn3adJ+Z4YeIAIYfIp+EUa+xXVRUFP744w+jy/AKww8RAez2IqIwwoUOiQhg+ClHr8X8iLTi71zgyBc6DMdWOyKyYbfXP6KjoxEREYGcnBzUqlUL0dHRqisQE+lBCIFLly4hLy8PERERppk19eOPP+Lxxx9H+/btMW/ePKPL8UjDhg0B2G54e+LECdSqVcvgiojICAw//4iIiIDVakVubi5ycnKMLofCSMWKFZGcnIyICHM0xP7555/47LPPcP78eaNL8VhsbCzq1auHnJwcZGZmMvwQhSmGH5no6GgkJyejuLgYJSUlRpdDYSAyMhJRUVGmamW8+uqrsWLFCtMGB6vVag8/V111ldHlEJEBGH6cWCwW1RuREpFtRfErrrjC6DK8ZrVasX37ds74Igpj5mhnJyLSCae7E5Hpwk9RURHatWsHi8WCvXv3Gl0OUdj5448/8Pvvv5v25qAMP0RkuvDzyCOPoF69ekaXQRS2HnnkEbRu3Rrvvvuu0aV4heGHiEw15mfjxo34+uuvsXbtWmzcuNHt9kVFRSgqKrJ/n5+fDwCm/R8rUTCIiopC9erVERcXZ8q/SwkJCQBsCx2ePn0akZGRBldERO5I/9botj6XMIljx46J+vXri507d4rMzEwBQOzZs8fle9LS0gQAPvjggw8++OAjBB6HDh3SJVNYhAj+ZU6FELj55pvRrVs3zJo1C1lZWbBardizZw/atWun+j7nlp8zZ86gYcOGOHz4MOLj4wNQeegqKChAgwYNkJ2djapVqxpdjqnxXOqD51E/PJf64bnUR35+PpKTk3H69GlUq1bN5/0Z2u01Z84czJ071+U2O3fuxI4dO1BQUIAZM2Z4tP+YmBjExMSUez4+Pp6/hDqpWrUqz6VOeC71wfOoH55L/fBc6kOvxWANDT8TJkzAXXfd5XKblJQULFiwAD/++GO5INOpUycMHToUb731lj/LJCIiohBiaPhJSEiwDz505YUXXsCCBQvs3+fk5KBv375YvXo1rr76an+WSERERCHGFLO9kpOTHb6vXLkyANtKs0lJSZr3ExMTg7S0NMWuMPIMz6V+eC71wfOoH55L/fBc6kPv82iKAc/OtA54JiIiInJmyvBDRERE5C3TrfBMRERE5AuGHyIiIgorDD9EREQUVhh+iIiIKKyEZfjJysrCyJEjYbVaERcXhyuuuAJpaWm4dOmS0aWZwssvvwyr1YrY2Fh07NgR27ZtM7ok01m0aBE6d+6MKlWqoHbt2hg4cCAOHDhgdFmmt2jRIlgsFkyaNMnoUkzp6NGjuOeee1CzZk1UrFgR7dq1w+7du40uy3SKi4sxa9Ys+zWmUaNGmDdvHkpLS40uLeht3boV/fv3R7169WCxWLBhwwaH14UQmDNnDurVq4e4uDikpqbi999/9/g4YRl+MjIyUFpaimXLluH333/Hc889h1dffRUzZ840urSgt3r1akyaNAmPPvoo9uzZgx49euCmm27C4cOHjS7NVL7//nuMHz8eP/74IzZt2oTi4mLccMMNOHfunNGlmdbOnTuxfPlytGnTxuhSTOn06dPo1q0bKlSogI0bN2Lfvn145plndLmPUrh58skn8eqrr2Lp0qXYv38/nnrqKSxevBgvvvii0aUFvXPnzqFt27ZYunSp4utPPfUUnn32WSxduhQ7d+5E3bp1cf311+Ps2bOeHUiX26OGgKeeekpYrVajywh6V111lRg7dqzDc82bNxfTp083qKLQcPz4cQFAfP/990aXYkpnz54VTZo0EZs2bRI9e/YUEydONLok05k2bZro3r270WWEhFtuuUWMGDHC4blBgwaJe+65x6CKzAmAWL9+vf370tJSUbduXfHEE0/Yn7t48aKIj48Xr776qkf7DsuWHyX5+fmoUaOG0WUEtUuXLmH37t244YYbHJ6/4YYbsGPHDoOqCg35+fkAwN9BL40fPx633HIL+vTpY3QppvXJJ5+gU6dOuP3221G7dm20b98eK1asMLosU+revTu+/fZbHDx4EADw888/4//+7/9w8803G1yZuWVmZuLYsWMO16CYmBj07NnT42uQKW5v4W+HDh3Ciy++iGeeecboUoLaiRMnUFJSgjp16jg8X6dOHRw7dsygqsxPCIEpU6age/fuaN26tdHlmM4HH3yA9PR07Ny50+hSTO3PP//EK6+8gilTpmDmzJn46aef8O9//xsxMTG47777jC7PVKZNm4b8/Hw0b94ckZGRKCkpweOPP467777b6NJMTbrOKF2D/vrrL4/2FVItP3PmzIHFYnH52LVrl8N7cnJycOONN+L222/HAw88YFDl5mKxWBy+F0KUe460mzBhAn755Re8//77RpdiOtnZ2Zg4cSJWrVqF2NhYo8sxtdLSUnTo0AELFy5E+/btMWbMGIwaNQqvvPKK0aWZzurVq7Fq1Sq89957SE9Px1tvvYWnn34ab731ltGlhQQ9rkEh1fIzYcIE3HXXXS63SUlJsf85JycHvXr1QteuXbF8+XI/V2d+CQkJiIyMLNfKc/z48XJJnLR56KGH8Mknn2Dr1q0e3aSXbHbv3o3jx4+jY8eO9udKSkqwdetWLF26FEVFRYiMjDSwQvNITExEy5YtHZ5r0aIF1q5da1BF5jV16lRMnz7dfj268sor8ddff2HRokUYNmyYwdWZV926dQHYWoASExPtz3tzDQqp8JOQkICEhARN2x49ehS9evVCx44dsXLlSkREhFQjmF9ER0ejY8eO2LRpE2677Tb785s2bcKAAQMMrMx8hBB46KGHsH79emzZsgVWq9Xokkzpuuuuw6+//urw3P3334/mzZtj2rRpDD4e6NatW7nlFg4ePIiGDRsaVJF5nT9/vtw1JTIyklPdfWS1WlG3bl1s2rQJ7du3B2Abi/r999/jySef9GhfIRV+tMrJyUFqaiqSk5Px9NNPIy8vz/6alCxJ2ZQpU3DvvfeiU6dO9hazw4cPY+zYsUaXZirjx4/He++9h48//hhVqlSxt6bFx8cjLi7O4OrMo0qVKuXGSVWqVAk1a9bk+CkPTZ48Gddccw0WLlyIO+64Az/99BOWL1/OVnEv9O/fH48//jiSk5PRqlUr7NmzB88++yxGjBhhdGlBr7CwEH/88Yf9+8zMTOzduxc1atRAcnIyJk2ahIULF6JJkyZo0qQJFi5ciIoVK+Jf//qXZwfSYzqa2axcuVIAUHyQey+99JJo2LChiI6OFh06dOD0bC+o/f6tXLnS6NJMj1Pdvffpp5+K1q1bi5iYGNG8eXOxfPlyo0sypYKCAjFx4kSRnJwsYmNjRaNGjcSjjz4qioqKjC4t6H333XeK/zYOGzZMCGGb7p6Wlibq1q0rYmJixLXXXit+/fVXj49jEUII33IaERERkXlwoAsRERGFFYYfIiIiCisMP0RERBRWGH6IiIgorDD8EBERUVhh+CEiIqKwwvBDREREYYXhh4iIiMIKww8RERGFFYYfIjKd999/H7GxsTh69Kj9uQceeABt2rRBfn6+gZURkRnw9hZEZDpCCLRr1w49evTA0qVLMXfuXLz22mv48ccfUb9+faPLI6IgF5Z3dScic7NYLHj88ccxZMgQ1KtXD88//zy2bdtmDz633XYbtmzZguuuuw4fffSRwdUSUbBhyw8RmVaHDh3w+++/4+uvv0bPnj3tz3/33XcoLCzEW2+9xfBDROVwzA8RmdJXX32FjIwMlJSUoE6dOg6v9erVC1WqVDGoMiIKdgw/RGQ66enpuP3227Fs2TL07dsXs2fPNrokIjIRjvkhIlPJysrCLbfcgunTp+Pee+9Fy5Yt0blzZ+zevRsdO3Y0ujwiMgG2/BCRaZw6dQo33XQTbr31VsycORMA0LFjR/Tv3x+PPvqowdURkVmw5YeITKNGjRrYv39/uec//vhjA6ohIrPibC8iCjl9+/ZFeno6zp07hxo1amD9+vXo3Lmz0WURUZBg+CEiIqKwwjE/REREFFYYfoiIiCisMPwQERFRWGH4ISIiorDC8ENERERhheGHiIiIwgrDDxEREYUVhh8iIiIKKww/REREFFYYfoiIiCisMPwQERFRWPl/p317F7xzWxAAAAAASUVORK5CYII=", "text/plain": [ "Figure(PyObject
)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "using Optim # Optimization library\n", "\n", "y_1 = zeros(length(y))# class 1 indicator vector\n", "y_1[findall(y)] .= 1\n", "X_ext = vcat(X, ones(1, length(y))) # Extend X with a row of ones to allow an offset in the discrimination boundary\n", "\n", "# Implement negative log-likelihood function\n", "function negative_log_likelihood(θ::Vector)\n", " # Return negative log-likelihood: -L(θ)\n", " p_1 = 1.0 ./ (1.0 .+ exp.(-X_ext' * θ)) # P(C1|X,θ)\n", " return -sum(log.( (y_1 .* p_1) + ((1 .- y_1).*(1 .- p_1))) ) # negative log-likelihood\n", "end\n", "\n", "# Use Optim.jl optimiser to minimize the negative log-likelihood function w.r.t. θ\n", "results = optimize(negative_log_likelihood, zeros(3), LBFGS())\n", "θ = results.minimizer\n", "\n", "# Plot the data set and ML discrimination boundary\n", "plotDataSet()\n", "p_1(x) = 1.0 ./ (1.0 .+ exp(-([x;1.]' * θ)))\n", "boundary(x1) = -1 ./ θ[2] * (θ[1]*x1 .+ θ[3])\n", "plot([-2.;10.], boundary([-2.; 10.]), \"k-\");\n", "# # Also fit the generative Gaussian model from lesson 7 and plot the resulting discrimination boundary for comparison\n", "generative_boundary = buildGenerativeDiscriminationBoundary(X, y)\n", "plot([-2.;10.], generative_boundary([-2;10]), \"k:\");\n", "legend([L\"y=0\";L\"y=1\";L\"y=?\";\"Discr. boundary\";\"Gen. boundary\"], loc=3);\n", "\n", "# Given $\\hat{\\theta}$, we can classify a new input $x_\\bullet = [3.75, 1.0]^T$:\n", "\n", "x_test = [3.75;1.0]\n", "println(\"P(C1|x•,θ) = $(p_1(x_test))\")" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "- The generative model gives a bad result because the feature distribution of one class is clearly non-Gaussian: the model does not fit the data well. \n", "\n", "- The discriminative approach does not suffer from this problem because it makes no assumptions about the feature distribition $p(x|y)$, it just estimates the conditional class distribution $p(y|x)$ directly." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Recap Classification\n", "\n", "\n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "
Generative Discriminative (ML)
1Like density estimation, model joint prob.\n", "$$p(\\mathcal{C}_k) p(x|\\mathcal{C}_k) = \\pi_k \\mathcal{N}(\\mu_k,\\Sigma)$$ Like (linear) regression, model conditional\n", "$$p(\\mathcal{C}_k|x,\\theta)$$
2Leads to softmax posterior class probability\n", "$$ p(\\mathcal{C}_k|x,\\theta ) = e^{\\theta_k^T x}/Z$$\n", "with structured $\\theta$ Choose also softmax posterior class probability\n", "$$ p(\\mathcal{C}_k|x,\\theta ) = e^{\\theta_k^T x}/Z$$\n", "but now with 'free' $\\theta$
3For Gaussian $p(x|\\mathcal{C}_k)$ and multinomial priors,\n", "$$\\hat \\theta_k = \\left[ {\\begin{array}{c}\n", " { - \\frac{1}{2} \\mu_k^T \\sigma^{-1} \\mu_k + \\log \\pi_k} \\\\\n", " {\\sigma^{-1} \\mu_k } \\\\\n", "\\end{array}} \\right]$$\n", "in one shot. Find $\\hat\\theta_k$ through gradient-based adaptation\n", "$$\\nabla_{\\theta_k}\\mathrm{L}(\\theta) = \\sum_n \\Big( y_{nk} - \\frac{e^{\\theta_k^T x_n}}{\\sum_{k^\\prime} e^{\\theta_{k^\\prime}^T x_n}} \\Big)\\, x_n$$
" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "##
OPTIONAL SLIDES
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Proof of gradient and Hessian for Laplace Approximation of Posterior\n", "\n", "- We will start with the posterior\n", "$$\\begin{align*}\n", "\\underbrace{p(w | D)}_{\\text{posterior}} \\propto \\underbrace{\\mathcal{N}(w \\,|\\, m_0, S_0)}_{\\text{prior}} \\cdot \\underbrace{\\prod_{n=1}^N \\sigma\\big( \\underbrace{(2y_n-1) w^T x_n}_{a_n}\\big)}_{\\text{likelihood}} \\tag{B-4.142}\n", "\\end{align*}$$\n", "from which it follows that\n", "$$\\begin{align*}\n", "\\log p(w | D) \\propto -\\frac{1}{2}\\log |S_0| -\\frac{1}{2} (w-m_0)^T S_0^{-1} (w-m_0) +\\sum_n \\log \\sigma\\left( a_n\\right) \n", "\\end{align*}$$\n", "and the gradient\n", "$$\\begin{align*}\n", "\\nabla_{w}\\log p(w | D) &\\propto \\underbrace{S_0^{-1} (m_0-w)}_{\\text{SRM-5b}} +\\sum_n \\underbrace{\\frac{1}{\\sigma(a_n)}}_{\\frac{\\partial \\log \\sigma(a_n)}{\\partial \\sigma(a_n)}} \\cdot \\underbrace{\\sigma(a_n) \\cdot (1-\\sigma(a_n))}_{\\frac{\\partial \\sigma(a_n)}{\\partial a_n}} \\cdot \\underbrace{(2y_n-1)x_n}_{\\frac{\\partial a_n}{\\partial w} \\text{ (see SRM-5a)}} \\\\\n", "&= S_0^{-1} (m_0-w) + \\sum_n (2y_n-1) (1-\\sigma(a_n)) x_n \\quad \\text{(gradient)}\n", " \\end{align*}$$\n", "where we used $\\sigma^\\prime(a) = \\sigma(a)\\cdot (1-\\sigma(a))$.\n", "\n", "- For the Hessian, we continue to differentiate the transpose of the gradient, leading to\n", "$$\\begin{align*}\n", "\\nabla\\nabla_{w}\\log p(w | D) &= \\nabla_{w} \\left(S_0^{-1} (m_0-w)\\right)^T - \\sum_n (2y_n-1) x_n \\nabla_{w}\\sigma(a_n)^T \\\\ &= -S_0^{-1} - \\sum_n (2y_n-1) x_n \\cdot \\underbrace{\\sigma(a_n)\\cdot (1-\\sigma(a_n))}_{\\frac{\\partial \\sigma(a_n)^T}{\\partial a_n^T}}\\cdot \\underbrace{(2y_n-1) x_n^T}_{\\frac{\\partial a_n^T}{\\partial w}} \\\\\n", "&= -S_0^{-1} - \\sum_n \\sigma(a_n)\\cdot (1-\\sigma(a_n))\\cdot x_n x_n^T \\quad \\text{(Hessian)}\n", "\\end{align*}$$\n", "since $(2y_n-1)^2=1$ for $y_n \\in \\{0,1\\}$.\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "### Proof of Derivative of Log-likelihood for Logistic Regression\n", "\n", "\n", "- The Log-likelihood is $\n", " \\mathrm{L}(\\theta) = \\log \\prod_n \\prod_k {\\underbrace{p(\\mathcal{C}_k|x_n,\\theta)}_{p_{nk}}}^{y_{nk}} = \\sum_{n,k} y_{nk} \\log p_{nk}$\n", "\n", " \n", "- Use the fact that the softmax $\\phi_k \\equiv e^{a_k} / {\\sum_j e^{a_j}}$ has analytical derivative:\n", "\n", "$$ \\begin{align*}\n", " \\frac{\\partial \\phi_k}{\\partial a_j} &= \\frac{(\\sum_j e^{a_j})e^{a_k}\\delta_{kj}-e^{a_j}e^{a_k}}{(\\sum_j e^{a_j})^2} = \\frac{e^{a_k}}{\\sum_j e^{a_j}}\\delta_{kj} - \\frac{e^{a_j}}{\\sum_j e^{a_j}} \\frac{e^{a_k}}{\\sum_j e^{a_j}}\\\\\n", " &= \\phi_k \\cdot(\\delta_{kj}-\\phi_j)\n", " \\end{align*}$$\n", "\n", "\n", "\n", " - Take the derivative of $\\mathrm{L}(\\theta)$ (or: how to spend a hour ...)\n", "$$\\begin{align*} \n", "\\nabla_{\\theta_j} \\mathrm{L}(\\theta) &= \\sum_{n,k} \\frac{\\partial \\mathrm{L}_{nk}}{\\partial p_{nk}} \\cdot\\frac{\\partial p_{nk}}{\\partial a_{nj}}\\cdot\\frac{\\partial a_{nj}}{\\partial \\theta_j} \\\\\n", " &= \\sum_{n,k} \\frac{y_{nk}}{p_{nk}} \\cdot p_{nk} (\\delta_{kj}-p_{nj}) \\cdot x_n \\\\\n", " &= \\sum_n \\Big( y_{nj} (1-p_{nj}) -\\sum_{k\\neq j} y_{nk} p_{nj} \\Big) \\cdot x_n \\\\\n", " &= \\sum_n \\left( y_{nj} - p_{nj} \\right)\\cdot x_n \\\\\n", " &= \\sum_n \\Big( \\underbrace{y_{nj}}_{\\text{target}} - \\underbrace{\\frac{e^{\\theta_j^T x_n}}{\\sum_{j^\\prime} e^{\\theta_{j^\\prime}^T x_n}}}_{\\text{prediction}} \\Big)\\cdot x_n \n", "\\end{align*}$$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "slideshow": { "slide_type": "skip" } }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "open(\"../../styles/aipstyle.html\") do f\n", " display(\"text/html\", read(f,String))\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "anaconda-cloud": {}, "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Julia 1.8.2", "language": "julia", "name": "julia-1.8" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.8.2" } }, "nbformat": 4, "nbformat_minor": 4 }