{ "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": {}, "outputs": [], "source": [ "using Pkg;Pkg.activate(\"probprog/workspace\");Pkg.instantiate()\n", "using Random; Random.seed!(1234);\n", "IJulia.clear_output();" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAG2CAYAAAB/OYyEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3wU9b3/8c+ySIAWYhGBIOFiDxaEKgKeaoUS4BS8FFGOttqitnjjAVo4PNqjKJiYRKhFqbeCYD3WFrX+qsHq45QWWhCwHB8iEY8XsLZeQIiHYDVRxCDJ/P6YTnZ2MrNz2ZmdmZ3X8/HIY8nu7Ox3N9F55/O9pRRFUQQAACABOoXdAAAAgEIh+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMSITfA5evSoLFq0SIYMGSLdunWTE088Uaqrq6WtrS3spgEAgJjoHHYDnLr99tvl/vvvl4cfflhGjBghL774ovzgBz+Q0tJSmTdvXtjNAwAAMRCb4PM///M/Mn36dDnvvPNERGTw4MHy2GOPyYsvvhhyywAAQFzEJviMGzdO7r//fvnrX/8qJ510krz88svy3HPPyV133WX5nJaWFmlpaWn/vq2tTf7xj3/IcccdJ6lUqhDNBgAAeVIURT7++GPp37+/dOqU5ygdJSba2tqUG2+8UUmlUkrnzp2VVCqlLFmyJOdzKisrFRHhiy+++OKLL76K4Gvv3r1554mUoiiKxMBvfvMb+fGPfyzLli2TESNGyM6dO2X+/PmyfPlyueKKK0yfY6z4NDU1ycCBA2Xv3r3Ss2fPQjUdAADkobm5WcrLy+Wjjz6S0tLSvM4Vm+BTXl4uN954o8ydO7f9vtraWlmzZo3s3r3b0Tmam5ultLRUmpqaCD4AAMSEn9fv2Exn//TTTzv066XTaaazAwAAx2IzuHnatGly2223ycCBA2XEiBHy0ksvyfLly2XWrFlhNw0AAMREbLq6Pv74Y1m8eLGsXbtWDhw4IP3795dLL71UbrnlFunSpYujc9DVBQBA/Ph5/Y5N8PGDkw9OURQ5evSotLa2Frh1xeeYY46RdDoddjMAADHnZ/CJTVdXIRw5ckQaGhrk008/DbspRSGVSsmAAQPki1/8YthNAQBARAg+7dra2uTtt9+WdDot/fv3ly5durDIYR4URZHGxkZ57733ZOjQoVR+AACRQPD5pyNHjkhbW5uUl5dL9+7dw25OUTj++OPlnXfekc8//5zgAwCIhNhMZy+UvJfCRjsqZgCAqOEqDwAAEoPgAwAAEoPgk6eqKpGaGmfH1tSoxwMAgHAQfPKUTovccot9+KmpUY9jjC8AAOEh+ORp8WKR6urc4UcLPdXV6vFxtWLFChkyZIh07dpVxowZI1u3bg27SQAAuMJ0dh9oYeaWW7K/Fyme0PP444/L/PnzZcWKFXLWWWfJqlWr5JxzzpHXX39dBg4cGHbzAABwhIqPT8wqP4UKPQMGDJAVK1Zk3bdt2zbp3r27vPvuu768xvLly+XKK6+Uq666SoYPHy533XWXlJeXy8qVK305PwAAhUDw8ZE+/JSUFK7Sc8YZZ8j27dvbv1cURebPny/z58+XQYMGZR27ZMkS+eIXv5jzy9iFdeTIEdmxY4dMmTIl6/4pU6bItm3bgntjAAD4jK4uny1eLFJbK3LkiEiXLoXp3jrjjDPkl7/8Zfv3v/71r2XPnj2ycOHCDsfOnj1bvv3tb+c83wknnJD1/cGDB6W1tVX69u2bdX/fvn3l/fff995wAAAKjODjs5qaTOg5ckT9vhAVnxtuuEE++eQT6dSpk9x0001SW1srPXr06HBsr169pFevXp5ex7gSs6IorM4MAIgVurp8pB/T09JiP9vLL2PHjpV0Oi319fXyk5/8RI477jiZNWuW6bFeurp69+4t6XS6Q3XnwIEDHapAAABEGRUfn5gNZM4128tPXbt2lVNPPVXq6upk9erV8swzz1juOealq6tLly4yZswY2bBhg1x44YXt92/YsEGmT5+e/xsAAKBACD4+yDV7q1Dh54wzzpB77rlHvvWtb8nkyZMtj/Pa1bVgwQK57LLLZOzYsXLmmWfK6tWrZc+ePTJ79ux8mg0AQEERfPLkZMp6IcLPqFGjpHPnzrJs2TL/Ty4i3/nOd+SDDz6Q6upqaWhokJEjR8rvf//7DrPGAACIMoJPnlpbnU1Z1x5vbQ2mHY888ojMmTNHvvKVrwTzAiIyZ84cmTNnTmDnBwAgaASfPLnZdNTvSk9bW5s0NjbKgw8+KG+88YasXbvW3xcAAKDIEHxibMuWLTJp0iQZNmyY1NXVSWlpadhNAgAg0gg+MVZRUSFtbW1hNwMAgNhgHR8AAJAYBB8AAJAYBB8AAJAYBB8AAJAYBB8AAJAYBB8AAJAYBB8AAJAYBB9EQlWVuu+ZEzU17lbMBgBAQ/AJQkODemVuaAi7JbGRTqubuNqFH21T2HS6MO0CABQXVm4OQkODyK23ipx/vkhZWditcWT/fvW2f39/j3XKyQ72WuhxsiksAABmqPig3f79mVBjtGXLFpk2bZr069dfTjghJX/4w1O+v/7ixWqoMav8EHoAAH4g+EBE1OpN//7W4efQoUNy4omnyoIF94mISK9ewbTDLPwQegAkBkMlAkfwKQIDBgyQFStWZN23bds26d69u7z77ruOz5Mr/Jx66jly2WW1MnPmDD+anJM+/JSUEHpgwIXBHT6veNGGSvDzCgzBJ18NDSL19R2/RMzvD+CX+YwzzpDt27e3f68oisyfP1/mz58vgwYNyjp2yZIl8sUvftHy66STvih7927NCj/av7Vg5Cer2VyLF4t06SJy5Ih6u3gxs7nwT14uDEm++Id5IfXyuSf5Z4WCIPjka9UqkTFjsr+uvlp97OqrOz62apXvTTAGn1//+teyZ88eWbhwYYdjZ8+eLTt37sz5dfbZY9srPzt2BBd6RKxnc9XUZELPkSMikyczmysS7C5KUb1o8Ve0vSB+dl5Dar4/qyB/D8P6HY/qf1sxxKyufF17rTp7S6++Xg09DzwgMnp09mMBzPI644wz5IYbbpBPPvlEOnXqJDfddJPU1tZKjx49Ohzbq1cv6eVggE63bup/X4oikkoFE3pEzGdzGcf0TJ4ssnGjyKRJdHeFzm7GYlRnNDY2Zt+io6j+7LwI8r2E9TkV088nZASffJWVWf8Sjh7dMfgEYOzYsZJOp6W+vl7+9Kc/yXHHHSezZs0yPXbJkiWyZMmSnOdbt26dfPnL49tDj6Jkqj5B0IefZ59VQ44WempqMqFn40b1e8IPXDt4UL19802R//kf9Q8W7b/bhga1Equ/D+HQfhZJ0dDQsYKjHyqht2tXYdqUAASfItC1a1c59dRTpa6uTlavXi3PPPOMdOpk3os5e/Zs+fa3v53zfKnUCVndW/rxPkGGHy30aJUdY+VH+147Hgng5sIgkvsPERGRjz7q+Fez/i9pkXiHIL8/r0JraBBZvTrsVhTOqlXq754ZbcgEfEfwKRJnnHGG3HPPPfKtb31LJk+ebHmcXVeX2UDm/v1FDh36RDZt+pscf7x639tvvy07d+6UXr16ycCBA/Nuv7GyU1Kiju3Rz+ZyssghiozbC0NlpRpajBf/3bvV23feUW/1fz3ru7+C7iIJOlR5+bz8mjHgJXRpz9Pofy5xCGr50oZKNDaK1NWJzJghsnev+VCJXbtEZs4Mr61FhOBTJEaNGiWdO3eWZcuWeT5Hrtlb+/a9KDNnTmz/fsGCBSIicsUVV8gvf/lLz6+paW3NhBwt9GizufS071tb835J2LG7kDU2ZrqQRER691b/p60/Ts/LRcvLGLpcF/8HH1Rv9ReQa65x1yavCjFGw+nnpf3sevfuGE68/uy8hC6R/INakFUu/bm1cDJ0qD/n1h9fX69Wuq69Vtr/uuzWLfvYw4fdvy7duOaUBGlqalJERGlqaurw2OHDh5XXX39dOXz4cP4vtH+/olRWqrcFMnHiRGXevHmen79vn6Js367e+nGconj7TKurFUVEUbp0UW+rqx0/FX6rrFR/CH59VVbm1x7tv6s//EE9344d5se99JKirFmT/XX55epzzjtPva2pyTz2yCOZ8+3Ykfvc+Qjy3G5f1+3P1snPbv/+zGeofT3wgPr8Bx7o+Nj+/R2fs2ZN5jWtnmPk9r1MmOD8/835fk5OrwX6n5H2b6+vq3/NsH7nApDr+u0WFZ8glJUVZMGZtrY2aWxslAcffFDeeOMNWbt2bV7nczJlPagxPiKM6Ykcu+pBebl5xcfLjEYnf5lqFZM1a3K3+6mnrKsI//3f6q3+F2r8ePV248bMa2tdLtpf0sX2l3MQs1HdTvQwq9RoVQ2RjhUPq/O7eS9ad1FDg7P3pD+39txFi0Rqa519TvlU+dasERk+3P49GV/XOGYNHRB8YmzLli0yadIkGTZsmNTV1Ulpaannc7kJNEGEH7NtKRjTEzIvMxa1MrzbGY1+dgNpF6tVq5wNlN26Vb398Y8z92ldYddco55v165guqmMgapQASsCs1Fzdo2JdBzPYjUeKcj3YnbuYcOcnzufZRSGDzc/f6F+PkWM4BNjFRUV0tbWFnYz8pZrLy7CTxHTX+T9Pt/o0epFUjv373+v/vJceqnIY4+p902cKLJpU+5zrl7tbZaR03EnWqA688xM8EnKWi25KjUiHSseQXweQQdNrSKqr4za/W6Ul6uBu7FR/d6PAd2NjZn/Hor998oBgg9C5WQDUsJPTLhdg8VrSb53b/Wvf7tuBe3xVavU54iInHxy5vgbbhC5445MF8aaNWpXix8zatwO9K2rE5k61fn53SgrM/+8wmZ3QbeqePjJKmjmCifaDEH9DDQ34cTp74YWuJ3OvNParLVr165M1+Fzz6ndc0OHZsJksc2Qc4Hgg1DpZ3PlwmyuCDLrpvFrDZZcF569e9ULlf4Yu3FB996rfn/ssdmP6y+s+uqCVXfCqlXqRchuvNLHH4v84Q+ZGTpa+42hSgtUMwLc/LdAYw4tX9tL6LrmGn8uysZF/4xVN31IMIZlq3BSW6ve6oOw3TIKu3dnXvPMM9WQ3bt35vfD6/gq/X8Dxq5dffus2pzQzQ8JPgiVm//uqPREgP5CZvXXst1FK9dfptrFIdf4HKtpzmaVI21sxUcfZd+KqH8FH398dju0ypDVBXP1apFvfKNjF4yxWrB8ucj3vud8jIY+IIVh506RuXP9Cxwar6Er3y6ZsjKRCROsq3TG36GZM7ODQK5uuGXL1BXAZ8zI/NzsllGorc2ED41Z8Mg1fscsRNqNkzKaMUPk5psz50sogg8A5/QXMqvNEu0uWsb/WesvTmahRhtgbPdXsVl76urUW+OIeZGOF6OZMzNr+uTq1nI66NZIC2FaqGpsVMOXiMjater9WoVg7Vr1cX0gCrJr4rXXRJ5/XuS664rjglhWpo7lMqsa6n+H9N2ckyZlP9/qc5g0SeRHP+p4/wUXZNb40axfL/KrX4lcfrnIlCnZj40Y4f49GX/P7Ga0deumvj9tr59rrmFgtBB8AF9UVak7xzupStXUqF12iaoy60vyWql/9241eCxapN6v/XvYsOxuAONFaPTojkGnocF88PBXv6re3nuvyNe/roaJJUtEtmxRX+vCC7NXzRVRqzpWU4lFvA+61UKYWagyVgOM34sko2vCz/FIZpU4bYq8dqtVGg8f7ti96jZo5lpG4Ve/Ur/0KitFRo1yfn4zZm3UArZ+OYCWFvV2796OCyAmcKxPrILPvn375IYbbpB169bJ4cOH5aSTTpIHH3xQxowZE3bTELKwg0c67WzwtX4wd6wZu6pEcq+W+8gjavePGf1FXvu33UXe7eDh9evV4HP88SLf/KYafPr2zfz1qw0sbmhQX3vSJOuLgddBtzNmZIcqreKjD3z6MDhuXMeKT7ELcjyS/nfGGD792M7DrPqydm3m53nhhdmPBfXzNAvY2n9XQW9bEhOxCT4ffvihnHXWWTJx4kRZt26d9OnTR/7+97/LscbBikiksIOHk5lnTmawRZbxr2H9GBwnF5EFC0R27Mi+z6wkr4UCu4uC2UWmqkrkmWfMj3/mmY6P/b//p4YhPe3Cq68giVjvIaXf+kE/UNV4nGbBguxQdfzx6kXpwgvVMFVfn/09/HPttWpXlP73LN9FG43HG5+j/d4MG5b75+lnpUsfsEXcL7qYALEJPrfffruUl5fLQw891H7f4MGDcz6npaVFWrQSn4g0NzcH1bzIcrOretA7sAcpCsEjVxtiHXpE3A+i1MblaOwWmdM4raaYnW/uXDXc6LuicoWhrVtFjNVi7a9fP3bNtvrrOuwLzc6d6pgevfXrs2/1RozIv0smCsrKMr8Xxt8zJ4sCBrksgJ+VLi2Aa39AVFaq3csiLH74T7EJPk8//bRMnTpVLr74Ytm8ebOccMIJMmfOHLk6x/+Eli5dKre6+Z91kXISaPQblPr1mkePOjvWr66nKAQPszbEPvSIONsWwDhQ1O8LhN2FR/8/fO1/7mZhyEn3g/H96tfxsdrsM9+pyYUyf77I5s3mj5mNRZkwQeTZZwNvVkHkE168hBNtlqB2Wwj696i12az6mGCxCT5vvfWWrFy5UhYsWCA33XSTvPDCC/LDH/5QSkpK5PLLLzd9zsKFC9t3ERdRKz7l5eWBtrO1tVW2bt0qDQ0NUlZWJuPHj5d0Oh3oa+aiBZlc4SfXruz5OHBA5PPPcx/jd9dTFIKHvg21tepO87EOPSLutgVw0lVldn67C5KXC49ZGNIGf44bl3vqsL4tZWVqFWv1aud/NUf1r+u77jKv+Pg1+yjKCr2mkfb7V8jlCsJctykufNg0tSCOOeYY5cwzz8y67/rrr1fOOOMMx+cIenf2J598UhkwYIAiIu1fAwYMUJ588knP5/SL1a7qTndbX7p0qXLyyScr3bp1U4YOHao88sgjtq/57ruHlXXrXlfuvNP8M9V2Yvd7B/bKSkWZNMnZLu/V1flvGm5Fe+0uXYI5f+iMOz+73Qna752jzc7n9D4/2+vm/Mbdu53u5u03bVf0NWsK+7phC3r38rB+nkZFsEu7n7uzdwo1dblQVlYmJ+uXmxeR4cOHy549e0JqUba6ujq56KKL5L333su6f9++fXLRRRdJnTbSPiRaNUer7oi4q/Rs3bpVfvazn8mrr74qM2fOlMsvv1zeeuutnM/p00ektFTknnvUqoteEFWYqir1vOm0utF2Oq1WW7p06fgaNTXqVk233KIe57eamsxrHznS8f1D/B8zEdWtGXLR/jrX2mz8HsEK+ncmKj/POP63ESQfglhBXHrppcq4ceOy7ps/f36HKlAuQVV8jh492qHSo/9KpVJKeXm5cvToUdfnduKEE05Qfv7zn2fd95e//EXp1q2b8s4772Tdr1V4XnzRWaXHzAcffKCIiLJ169acx2mf6Z13Hs6quuRT6amszF290c6rVXxSqY6vpR0XRLXJ2A6z74uG8a/IqPx1q2fWpihVfKIiqRUfxIafFZ/YBJ8XXnhB6dy5s3Lbbbcpb775pvLII48o3bt3V9a4+A81qOCzadMmy9Cj/9q0aZPrczvx7//+78r3v//99u/b2tqU008/XVm4cGGHY2+77TalW7cvtH994Qsdv7Zs2WL5Wm1tbcr3v/99ZeTIkUpLS0vOduk/U+3ib9f1ZMcuROhDjfalhaDqauvQ41eXl1X7ijL8RDHoOOG13U6fF8fPheCDiEtk8FEURXnmmWeUkSNHKiUlJcqwYcOU1atXu3p+UMHn0UcfdRR8Hn30UdfndmLZsmXKiBEj2r9/+OGHlb59+yrNzc0djn311Q+Uuro3lbVr31Tq6t5UnnvuTeXNN7O/Pv30U8vXmjVrlnLSSScp7733nm27jJ+pX2NecoUIfbAZMqRj+LEKPX6EEqehrKjCD4rDSy8pyoQJ6i0QQYkNPvkq1orP1q1blU6dOikff/yxcujQIeWEE05QHnjggQ7HGQcyOx3YrHn55ZcVEVF2797t6PggKj4asxChDz1a0DEGHu2+XOfxqz35HAcAyPAz+MRmOnuUjR8/XgYMGCD79u0TRVE6PJ5KpWTAgAEyfvz4QF5/7Nixkk6npb6+Xv70pz/JcccdJ7Nmzco6RhvI/MQTS+S++5a036/FgVRK/RIRWbdunWlb3377bRER+cpXvuKqfStWZA9k1gY2i3gf2Gyctq7/t/F10ml1nSAR9T1u3JgZbGw1wNrt2kKtrc4GamuPa+0BABRY/jksPoKczv7kk08qqVRKSaVSHQY2p1KpwKe0jx07Vpk3b57SrVs35U9/+lPWY/rKzgcffNCha+u557K7vay6uj788ENl+/btjtt0+PBhZdu215VBgw4HNuZFO0+nTubn0ypM+q8hQxSlosK+u4yqDABEA11dHoWxjk95eXlB1vG57rrrlFQqpUybNi3rfqfdWU6Oq6urU77yla84blOh1vGxGjtkHOica5yP320CAPiHrq6ImjFjhkyfPj2UlZtHjRolnTt3lmXLlnV4zMk6PU5WbG5qapI33njDUXv271dXbi4tFZkzx/wYJ/tr2TFbL8fYnaatCm3s9nr22SLbTwsAYM+HIBYbQVd8wjRx4kRl3rx5YTej3b59asXHyWfqdSq51Xo5+sqOvnKj3Z9OdxzoTKUHAKKLig9ERKStrU0aGxvlwQcflDfeeEPWrl0bdpPa9e8v8tlnIv8cD52T10qPsTqzeLFaxdm4Uf1e/1hNjXr/pEnZtxs3inTu7HxwshNVVWplycm5/NqgFQDgTGy2rEBHW7ZskbKyMlmzZo3U1dVJaWlp2E0qCKsuKS3cDB5sffyf/6zeauFHRA0eZttaeJVOq69nt02FftYZAKAwqPjEWEVFhbS1tYXdjILKNQ5HX7XRjtMqQMbKkEimMiSSPT4oX07GLjGeCADCQfBBrOTqktJ3F+m7vayOf/bZ7IHP+a4tpJcr/BB6ACA8BB/EipuxMOPHi1RUmC9OaBY8Nm4MPvwQegAgXAQfFC2zkGQVPBYv9m9VaT19+KmtVbvUCD0AEB4GNxsoJltOwJsofJZVVZlBxnbVFm2GVXW1s8HJTi1enFlnyM9B1AAA9wg+/3TMMceIiMinn34ackuKx5EjR0RECrKAoxX9DKtc44P0M6wWL1aP82s/LbNFFgEA4aCr65/S6bQce+yxcuDAARER6d69u6S0XTvhmrbGUPfu3aVz5/B+zfRdTXahx2zmV76M5/a7Kw0A4A7BR6dfv34iIu3hB/np1KmTDBw4MPQAGdYMq1yBivADAOEg+OikUikpKyuTPn36yOeffx52c2KvS5cu0qlTNHpTCz3DKte5CT8AEB6Cj4l0Oh3quBQEo1AzrJwEKsIPAIQjpURh6k2BNDc3S2lpqTQ1NUnPnj3Dbg5CUlKSGWzc0uL/+dmrCwD85ef1Oxr9EIgd/TRxOzU14VzYzdpoNcMqVxvdvlcR5xWcxYsJPQBQSAQfeBKHjTiNbdR3QbW0ZNbrmTw5dxvj8F4BAA4pCdLU1KSIiNLU1BR2U4pCdbWiiKi3Xh4vBK0NkyaZt0W7f9IkZ+eJ8nsFgGLl5/Wb4IO8WF3woxQErMKNXSgyisN71VRWOm9PdbV6PABEFcHHI4JPMIwX/igFAatw47XNfrzXQoQSr+8HAKKI4OMRwccdNxdoLVh06RKtC6n+PWgXeas2Og0Zdudx+vygQwndcwCKBcHHI4KPO24v0Ol0JhAUitNwVlmphjMnbXQSgLTQ4/W9FiqUxKl7DgCsEHw8Ivi45/QCHVbFx+kFXGufXRudnC/fio/da/kdSqLcFQkAThB8PCL4uKdVSnJdoPWDh7UKjJMLq19jXZyGM+2rutr8OW5Cj18holChxK+wBgBhIPh4RPBxz25wsD70uL2A+znWxa56YhZ4rP7t9TXyDT9Bh5J8u+cAICwEH48IPt5YdWdZhR7j8/IJP25ChVX1JFdlx+qYINtpVumyCiVeZ3VZtY+KD4A4Ivh4RPDxzjiAWbu1WwPHr0qK024x7TmdOuUONMYg4Eeo8XqcVSgJaoAzY3wAxA3BxyOCT360kKOFniFDnF/o7aoWdhdnt6FCCz+5jtdXWczaqIUtJ6FLf5ybKfF2XYh+hx67+wEgigg+HhF8vDNekPUVH02+g5XtumPcDGK269Jx0vXjd6XHyGpFaafbaNjxs3sOAMJE8PGI4OON8QKphQUt/HitypgdZzcA1+kgZifHOml3UOEh6IpP0KENAAqJ4OMRwcc9u/EobmZz5Zoab3buigp3bcon5HgJP/mGniDH+LBXF4BiQvDxqBiDT5AXOKfjbpyGn1xdOPrn5AoyxuPNBicbP5NcVRT9Z6Jvg9U5zN6j2881jFldABBnBB+PijH4BNmlYbbPlVXVQ1u80Or17MKHMfQY/21GH3oqKnK3VR+69MHGqi2DB1s/pq/MaPdZVaeccDLeCACSjODjUTEGH0UJfhCr1/M76cKxCzp257aaVm98ntbNVlFhP4amujpznPEY40wwu8qUHaddcQCQZAQfj4o1+CiK/+NQ3D7f6jizLhz9WB9j6DF2g+lDiDEcaFUWp11uuapOud6TsY36SlMQn23Q4YfxPwDihuDjUTEHH0UJpnqQz0XSruKjdSc5GSht7J6yqgy53VrDyXsyhhz9fV6nnAddpfPy2nbjmoznIBABKBSCj0fFHnwUJTrjRXKFMH3FR+uq0ipCZkGlstK8SmT1msaFFY2vpe/+cvL5mO3s7jREmb1GrkpPrnFSfrLrVrQLPXTHASgkgo9HSQg+ipL/ZpR+LURodcHXAoMWUIyBxPh4dbV96NFoxx17bHZQMQYsY3vsdn43fuXqTrMLL3ZhyDioPFdlJZ/KS67wQ+gBECUEH4+SEHz8qPjk0xViFyaMVRntVgssWveX/n5FMR9srL/o6ytC+kHMxtBjbJuTKfa5go/xPXsZs+O1i9KPEGL12mGMPQIAKwQfj4o9+Pg5xsdLV4hdANAYKzrGrihjV5LZ9HHjaxgDjOm/WscAACAASURBVDG0aPcbp7Xnaq9xBpixjfrjvW7W6vZ+p4+7YRaW/fxdAoB8EXw8Kubg4/UC6vacuSoCdisz65+vVV2077VAkUpln9+s284q5BgHMOsDi9X9ubqsjMcaA1uuEKidw8naRn6FpXyYfc5+VA8BwA8EH4+KNfgEWR3wqyskVygyBpSuXbO/z1XxcTI+SN9mfaXHKtQZ1wbKdX63n4XbMFGIykuuNuU7XgwA/EDw8agYg08hxoP40RWiDyxm95tVfMyqKvrX0J6TTmefyzg+yNjtlatSY1V50tPCj9Uxdp+F2zARZOUl18+Rig+AqCD4eFSMwcduBpZfM4T86Aqx607SQoc2I0v/mHFPLuNz9NUkfXVJ3zaz8Tn692w21sjqPRkDknExRbsA6DZMBFF5satKBV1pAgCnCD4eFWPwseNHRUh7rFMn510hViHKbjyOcbaXPogYt4kwG+OjZ7a9hHauwYOzjzWOGzKuFm32HrRzDRlivfdXrs/Y7c/Gz8qLk5+320oWAASF4ONREoOPouQ3Bsis68OuK8Tu9axmdeUanKwPGU4uysa26R8zhhyrKe1m78N4n9VYI6efcVBhKRe3P2/CD4CwEXw8SmrwURRvF17jY8bVk3N1W+VaFLCyMhNgtK4irUJkPJ9WndG/5pAh6uMTJphflPXdXvrz6HdQN47Tsdp6wi4IKErHsUZ2n6PTx72GJTtW3aNWQcfpApUAEBSCj6IoS5YsUUREmTdvnuPnJDn4KIq76kGuv/atqi5mM6pytUMLHWaBwdilpf1bH5TM2phrnR7jsdq5Uilnn1uuEKIfa2T3OeZ6DaeVnSDCB5uXAoiqxAefF154QRk8eLByyimnEHxccjJeJNdF1W4tG+OMKjNme3Xpj9e/vll3l1lFp7o6e8FB4+KDxmOdVnz0bdVXjMw+J7P37iVMmFXa8qnQAEDcJTr4fPzxx8rQoUOVDRs2KBMmTAgt+MT5r2O7GUJ2F1rjGBxjIDFbr8esG8tszI32WEWFeaXFrjtL321m1g5jJUn/+nbv2ex+Y7vsql1OeN0lnS4oAMUq0cHn8ssvV+bPn68oimIbfD777DOlqamp/Wvv3r2+fXBeuzDC5qTiY8UsvOjHtxiDgNmqxXaLGGpT2XNtSqo/l/YaxgqP3esYZ3W56fZz+70fwuj6AoCoSGzweeyxx5SRI0cqhw8fVhTFPvhUVlYqItLhy6+urrhdjPy+QJttCeHkeKvjjAOezdqmha9cY5CM08ydvr7ZOY2vZwxU+ucEXX2xOmfUfs8AwG+JDD579uxR+vTpo+zcubP9vjArPpq4XIz8bqddxcfqeLsdzHOFHmN7zdquVbKMY3e0mWF24Uzf1Wb1esb1gyoqrLvVzGZE5dPtWYjqEgBETSKDz9q1axURUdLpdPuXiCipVEpJp9PK0aNHbc8R1ODmqF+M/K5MGUOMFjZyjZXJNcbHeD79zC279ppVWYxr/pit5WP33vQDmc3Cj76r0K6SZHaefOTTXQkAcZTI4NPc3Ky88sorWV9jx45VZs6cqbzyyiuOzhHkrK6oXoz8Hotk1+Vjtzu78TxWocltWDO+vnHgtX71Zi8BUP+Y2YrQTnel9+v3gs1DASRJIoOPmTBndZmJysVIX12xm31mrJgYFwU0HmsMF24GMBtpn5dVN5nToKIfsKw/3mxVZf1zrNpuRTtGf04nYc3v0BPVkA0AQSH4/FOUgk+ULkb5VHlyPVc/RdxurE6ulZv1x1kt/Of0vVhtLGpVATJ2TTl5T8bXMn7l+vyCCj1R7VYFgCAQfDxK0hgfr906Vuvf6J9jty2F20qN124iY3gybmJq1Q1m1hVmHLBsfC2zhRutus+CCsGFqigBQNQQfDwKIvhE+WLkpW25uqvcVEf0XWhWr51rwLOT8+nbYxwrZBZGtLV+jGN1zNYgMusyM4YeY3VJ/3p+d3vmMz4JAOKO4OOR38EnDhcjL9Uoq/DjdgyPXVvcPm53vNkK0GbvW18VshtcrS2oaLY/mdk59StO+1Xxyafr0i9xXqkcQPwRfDzy84OLwsXIKS8XYmPIcTJl3cn5/Pq8zI4zdtNZvW99SDF7r8YgZTWex+y5xgHWfvz8oxA64vT7DqD4EHw8SvJeXV66XoxjaMxmXzn9HLRQ4uRzsPu83FSO9O/bzfgl43vWjwMystuVvljCQBwqnACKE8HHo6Tuzp5P14vdthR+VwKcTL/XV3SsApJZd5bdzvHac7QFFPWrSNtVfJyeO+6hwOp9FMv7AxBNBB+Pkhh8vIzxMT7XblsKPysB+sqLGbO9s3INota6s+yqMhr9WB6t68o4LshqzJPTKfxxDwf5/E4BgBcEH4+SFnzy+evc7RgfPysBTjcTddMWJzPSrHaIN17gjdWzuHV7+sHvAdwAkAvBx6MkBZ98qjBWwcJt4MinEmAVfryEHifPt6oIWb0nt+OlilFUVioHUPwIPh4lJfjkM+7GabDQP262Ho9VJcBNxcMYfvINPVbvwbipaa7z6oNPkisdVHwAFJKf1+/OgqLT2ipSXS2yeHHu47THW1uznztpksjGjSIVFdbPueUW9bjWVpF0Wv1+40aRTZtEamtFjhwR6dIluw01Nepx1dXO3sef/ywyebJ63s6ds9tm9v7M3ndVldo+/X369715c+Y9DBqkfhnPq3/PIpn2b9yYuc/usy4m+p/j4sWZ70WS9TkAiCkfglhsJHk6u1Nm1Y1c3WHGKohxTJAf3V76KeVuz+HkPZhtOGp1rDZY2sn5ixGzugCEga4uj/z64HKtB2PkdcG/sBgDnV1w0LaBqK523zXlhLFrKZ12fw59G83elzG4GUOqn7PW4ozPAUBYCD4e+fXBOb2g281Oigsnf+XbDRr28hlYjcVxey59sDFWsvT/NmtrPuOligmfA4AwEXw8CmLLCrup0XEPPZpcM7asHstnmwu/gpRx3R+7AGTccsJJl6b+Ney6NOPU7alXrF27AOKB4ONRUJuU+lnliDKzmTxWf+Ebpzq7qQT4GSrNKjpW78HJ2Ca71/DjOABANoKPR0FMZ/erOyYIQfyVbgw0Zq9hNdXZaUXE725E/cBkY/gxhh79c9xULRj/AgDBIfh4FNQ6Pn4MwA2C35UIq0CT61z67512G2nhxDiDyshqLaFcbde6soxfVu/Fj/DjZ+ihywlAEhF8PApyAcN8plwHya9KhL4SY/U845gZ7aLrtIqjMZuFlatdWhiwe46+2qPfgNSsauQ1rOQKfn6gWw1AEhF8PAoq+Jh1c0XpouO2EmE1pd1qC4uKityDnXOdw2lbvbxPsxWlrb7009jzDQ1OKmP5oFsNQNIQfDwKIvgEsXZNENxUInKFGLPHnA4Mtvts/Lhgm7VPX+nR/9vsMb9CQ9D7WBWiWw0AooLg45HfwcfrRpphsapE5BqgbFa50bqitMe07528X6uuKT8v2PpzWYUe7TNwOt7Hy+sHvY9V0N1qABAVBB+P/Pzg7GYVRTX8mM3KsmqjcayN/j0bQ0I+a/UEERC0cxvHXpm1WX/c4MH+vG6hwkihQhYAhIng41EQKzc7PS4Ks2vMLpJOu5+0ncv1O5j7cdENsktIO7c+/Fh1SeqPy3dsT6G7n4LuVgOAsBF8PPJzry6nFzGtaydsTsbqWIUfLexoM6H0j+dz0S1ExUffpWW2gapxnJLXkBrWgGMqPgCSgODjUZDT2aPMSSXCKvxo32uhRx9wvF50zbrXcrXRbRCxGoBdXZ0d1HId5/X1/DjO6+syxgdAsSL4eJTE4OOmEmEMP1rXj1bxybVVhfH7XFUxs5lwfg14tmqHPtjoq0BWA7fdvGYYFcCwutUAIAwEH4+SFny8VCKMXT/G7q1clRGz81gdYxwobTx/PoOlzZ6jb7fZ4Gyn53EriFAUVrcaAISF4ONR0oKP14uucSaUfhC3XZdQrvBj/N5qZpyXTV6dhJ7q6uxzO31OPvzuBgurWw0AwuTn9buzoGhVVTk/dvFi9XbyZJHWVpFUSo03Q4aIbNwoUlOjPn7LLSLV1Zl/65+r/7d2XHV15rjWVvX7xYvV823cKDJpUub8xvvHj3fefv259WpqMm1ZvFg9TkQ9t76tVu9DO96rXK9h1cZcrN6nXlWVSDqtHmfX/poa9Rg3vysAEGs+BLHYSFrFxy2tGmKcuq7fNNQ4JsZqJpS+gmRX+QlqZlKUqiOFHJMTpfcNAH6gq8sjgo81fReQftaV1Vgc4/OcXmStwk0Qa9FEbdmBQs7CYhwQgGJC8PGI4GPO6dYbTqegW7EKN0FVfKKokO+VmV8AigXBxyOCT0f6cGP3uL4i5PbiaXXBL2QVJCoKudJyEj9fAMWH4OMRwacjJ91BxrE8XkOP0ynrxXxxDqO6laSKGoDiRPDxiODjndeLp1WIcbrJazFdpMOsvrCfF4A4I/h4RPDJj9uLp90YE7uB0cUUfsIcb0PFB0DcsY4PCm7iRJEjR0S6dFFvtXV3zGhr8Tz7rPmaM8b1fKzWuPFrLZ2w5Vqnx8k6P36+dq7PGwCSgOADWzU1aogREVm0SL21unhqF9aKCuuF9vSL5dmFm7hfnJ0sThhU+DF77aCDFgBEHcEHOekvniLmKzJrF083KxDrFfPF18lKyyL+V7fCrDIBQJQRfGDJ6uJpFn7093MhzfCybUi+wqwyAUDUEXxgyurimWsvLkJPNIRVZQKAOEgpiqKE3YhCaW5ultLSUmlqapKePXuG3ZxI0za6zDWA+ZZbMoOdKypENm0qZAsBAEnh5/Wb4APPSkoyM71aWsJuDQCgWPl5/e7kU5uQMDU1Hae3B62qyvnr1NS4G18DAEgGgg9c04//aWnJjPMJOvyk085eR2tfOh1cWwhhABBPDG6GK2GuDePkddxOqbcby2Q8d2tr5jlO3q9xOQA3vLYNAJBD3ms/xwhbVuTHbnuFQm0x4ef2D06fY3Zc0J9HPm0DgGLCXl0eEXy8i9pF2M8NP/MJMEHvwRWVsAkAYSL4eETw8a6y0t2O7JWVQbYm8zp+bb6ZT4AJetf1MDc4BYAoCD34fPrpp8p7773X4f5XX3017wYFieBTfNzuGJ9LPgEm6B3Qgw5XABBloQaf3/72t8qAAQOUU045RfnqV7+qPP/88+2PnXbaaXk3KEgEn+ISRNjI55x+hjC/2wYAcebn9dv1dPba2lqpr6+Xl19+Wf7rv/5LZs2aJY8++qg2UNqnIdfmli5dKqeffrr06NFD+vTpIxdccIG88cYbgb4moimoKfWLF2fWJurSxfkMtUKsa+S1bQCADNfB5/PPP5fjjz9eRETGjh0rW7ZskVWrVkl1dbWkUinfG6i3efNmmTt3rjz//POyYcMGOXr0qEyZMkUOHToU6OsiWqym1PsRfrwEmEKtaxTGopEAUHTclogqKiqUl19+Oeu+lpYW5ZJLLlHS6XTeJSg3Dhw4oIiIsnnzZkfH09UVf0HOcvIyjqZQA48Z4wMgyUIZ49Pc3KwoiqLs3btXaWhoMD3mueeey7tBbrz55puKiCivvPKK6eOfffaZ0tTU1P61d+9egk+MBTml3kuAKdRUc2Z1AUi6UILPqaeeahl4wtDW1qZMmzZNGTdunOUxlZWVioh0+CL4xFNQU+q9BJhCrWvEOj4AEFLwufLKK5WBAwcqu3btyrq/vr5eOeecc/JuiFtz5sxRBg0apOzdu9fyGCo+sOM1wBRiXaOoLRoJAGHxM/g43qvrF7/4hdx6660ybtw4eeqpp6RPnz6yaNEiefLJJ+X888/3d+CRjeuvv16efvpp2bJliwwYMMDyuJKSEikpKSlgyxA3ra3O9vXSHm9tVW/d7InldfaV17YBAKylFMXdHPSlS5dKdXW1tLa2ytSpU+XWW2+V0aNHB9W+LIqiyPXXXy9r166VZ599VoYOHerq+c3NzVJaWipNTU3Ss2fPgFoJAAD85Of12/F09oaGBvnhD38oNTU1cvLJJ8sxxxwjl1xyScFCj4jI3LlzZc2aNfLoo49Kjx495P3335f3339fDh8+XLA2IFqqqpxP666pYfdyAEg6x8HnxBNPlK1bt8pvf/tb2bFjh9TV1cmcOXPk9ttvD7J9WVauXClNTU1SUVEhZWVl7V+PP/54wdqAaEmnna2Zo621k04Xpl0AgGhyPMbnoYcekksuuaT9+6lTp8qmTZvkW9/6lrz77ruyYsWKQBqo57JXDgmgjW+55Zbs7/XMFjwEACST4+CjDz2a0aNHy7Zt2+Tcc8/1tVGAG7nCD6EHAKDnOPhYGTx4sPzlL3/xoy2AZ2bhh9ADADDKO/iIiHzpS1/y4zRAXvThp7ZW3c+K0AMA0HM9nT3OmM6eDCUlmc08W1rCbg0AIF+hTGcH4oAdzAEAuRB8UDT0Y3paWtRbJ1PdAURYQ4O6AFdDQ9gtQZEg+KAomA1kXryY8IOI4OLtXUODyK238tnBN74MbgbClGv2lpN1foDAaRfv888XKSsLuzVAohF8EGtOpqwTfhC6xsbsWwChIfgg1tjBHLFw8GD2LYDQEHwQa242HaXSA0RcQ0PHsTz19dm3emVldB3CNYIPAPjJ7OK9e3fm1ngB5+KdsWqVOhbKzNVXd7yvstLdXz+AsIAhAPirqsr64m2Gi3eGVcXn6qtFHnhAZPTo7McIjYnh5/Wbig8A+OmCC0SGDs2+b/16kV/9SuTyy0WmTMl+bMSIwrUt6nIFmdGjOwYfwAOCDwD46amnrCs+v/qV+qVXWSkyalTw7QIgIgQfAPDXtdeq6/XorV2r7py7aJHIhRdmP0ZXDVBQBB8A8JOxu6ahQeT119V/DxtGdw0QMrasAIAgNTSI1NWF3Yr4KitTuwOpjMEnVHwARFZra6ts3bpVGhoapKysTMaPHy/pdDrsZqlhZtUqtVvLzQW5d+/g2lSsysqY9QZfUfEBEEl1dXUyePBgmThxonz3u9+ViRMnyuDBg6UuCtUTrxtnHn98MO0B4BjBB0Dk1NXVyUUXXSTvvfde1v379u2Tiy66KBrhx61rrqG7BogAuroAREpra6vMmzdPzNZWVRRFUqmUzJ8/X6ZPnx6Nbi89/QJ8jY3q2B5tTZ/TT++4QB8L8AEFR/ABEClbt27tUOnRUxRF9u7dK1u3bpWKiorCNcwJtlwAIo/gAyBSGhyOm3F6XN7cbJx55pkif/iDOpZn1y6RmTPVtXtqa623XABQUAQfAJFS5jAMOD0ub/lWcYYNU2+LccsFr7PbgBARfABEyvjx42XAgAGyb98+03E+qVRKBgwYIOPHjy9Mg8xWYrbbODNI+YYNP8OKNrvt/PMJPogNgg+ASEmn03L33XfLRRddJKlUKiv8pFIpERG56667Cjew2cnGmfruMO3fWjfY7t3q7a5dzs5pJ9+wQVhBwjGdHUDkzJgxQ5544gk54YQTsu4fMGCAPPHEEzJjxoyQWmZh1SqRMWOyv7RusNpa9XbmzMxjq1aF19Z8NTSoXXmFGmMF+IyKD4BImjFjhkyfPj2aKzcb5eoOmzdP5O67Re69V+TrX1cfi3OlRV8xAmKI4AMgstLpdOGmrOcz9iVX19Xpp6u3X/96vAc3a114Wpfdrl0ihw+r/zbObhNhjSJEFsEHAETcjX0pxMaZDQ0id9whMmVK9lYXVlPptXYZd4a3WlDRbVgxzm6bOTPzb9YoQowQfADArUJsnNnQILJ8ufplxknY8HNBRa07T1ufaM0ateIT1uw2wCOCDwD4RV9h0bqEtFldXruD1qwRGT48872bqfT6sUdeFlQ0W7xR6946fFikWzf139ot3VuIAYIPAPjFrMKizery2h00fLj52CAnCyKaBRE3Cyo6rRhp3V50byEGCD4AksfNNhQizisZUVvsMF9276dbt0y31/DhmffDis6IMIIPgOQJajNRJ4sdmnETxPQLIVrJdT43CyravR+NsSrFIomIMIIPgOSJWmXGbRDL53z6BRU1dsFOX8FpbFTva2wUOeWU4Ge3AT4j+ABIHq+VmaC4CWKNjSLr1+cOG7nOt2yZyJtvisyYkZkmbxdc9BWcgwfV+w4eLMzsNsBnBB8A8JOX8S1ug9jUqd7PN2mSyI9+5KxdXmndabt3x3vRRhQlgg8A+Mk4vqUQix0WUu/e2bdm44k2bszcarPINEx5R8gIPgAQJLPuoDjPetK6x7TbXOOJHnxQ/dJjyjtCRvABAJHCVmaiPOvJuM3FwYOZrqu1azPHrV2rzgw76SR1OruIWgU6/niRFSvUwHPllSJz5mSfP2rvF4lD8AEAkegN1PU7iDk9n5MZYcZ/a7RqzuDB6veDB/s/xifO1TJEAsEHALwKaiFE7Vh9EMv3gu802OlnhOkrPrW16nYXIpl/Dxsm8tFHIh9+qN5/7LEijzwi8vzz6vfPP69+rzdihMioUe7br4lytQyxQPABAK+CWgjRTKEu+GbhrL5eDTsXXqh+r/179GiRigqRzZvNz/Xf/61+6U2YIPLss363GnCM4AMAXjlZf0ermoioY2Dq6/2rCkXBXXeJvPZa9n2PPaYGnvPOE7n00uzHRowoXNsAEwQfAPDKyfo7VVWFqwoFQetiszJqVMeuqw8/VIPP2WeLfO97wbYPcIngAyB+4jTANWrbY7jV0CCyerW753zpS9m3+bx2UGOokFgEHwDxE6cBrl62x4jiBf+aazKvUahp/4UcQ4XEIPgAQNQqSGFe8K0+C/33hQoXca+WIZJiF3xWrFghy5Ytk4aGBhkxYoTcddddMn78+LCbBSDO/Kwg+bH+TpgXfD8+ixEj1Nlb+Q5kjtpmsigKsQo+jz/+uMyfP19WrFghZ511lqxatUrOOeccef3112XgwIFhNw8A/FkI0ekFX1+dyYd2ngsuUFdjFsnceuliGzWKKeuIrFgFn+XLl8uVV14pV111lYiI3HXXXfLHP/5RVq5cKUuXLg25dQACEcXxLvnya1VmvypV2nn0A5lnzsw+hjE1KBKxCT5HjhyRHTt2yI033ph1/5QpU2Tbtm2mz2lpaZGWlpb275ubmwNtI4AAFOMA16htj6GZMUPkG99QQ8+aNSLDhwfbxWY3tipqY69QFGITfA4ePCitra3St2/frPv79u0r77//vulzli5dKrda/Q8TQDz4Pd6lGCtIXmmfhdatdfCgyOHD6r+1227dMrd+fxZ2FSuzxwu5mSyKUmyCjyaVSmV9ryhKh/s0CxculAULFrR/39zcLOXl5YG2D4DP/B7gGtcKUhAXfONnoe/eMn4WM2dG47OIarUMsRGb4NO7d29Jp9MdqjsHDhzoUAXSlJSUSElJSSGaByAu4jxF+vzzsytWXipV+uefeabapZVrE9IPPxS5/nr1uEmT/H9PQIHFJvh06dJFxowZIxs2bJALtY3yRGTDhg0yffr0EFsGIFZjMeI6RdqvSlWu89TWdvy39nkMHx79ny3gQGyCj4jIggUL5LLLLpOxY8fKmWeeKatXr5Y9e/bI7Nmzw24akGxxWkk5rvyqVNmdp1u3zOBmEfXf+lWbvbIbW6XfzFVE3dB1797s4/SKeewVAhWr4POd73xHPvjgA6murpaGhgYZOXKk/P73v5dBgwaF3TQAhZTEAa5+VarszqMZPjzzbz8qeW4rVnaPR2G8EWIpVsFHRGTOnDkyZ86csJsBIEyFHuAap668fDU2qre7d4ts3uzfee0qTeXl5hWfOIy9QqzELvgAgO/sKkj5dOXFNTT94x/ud2XPxUvFSuviivLYK8QOwQeAO8W4Dk6QFaQ4jX/S1vMREfnoo8y/4/pzBUwQfAC4E9d1cOxErTLjpD35jnXS79E1YUL2Oj6LF2f+HeefK2BA8AHgThzWwfESYqJWmXHSnnwrVfrXeOyxzCrOM2eK1NRkwk9Ufq6ADwg+ANyJwzo4+YSYxsaOXTtx78qzog1k3rZN5EtfUv+9e7d6++abmePef18NRL17ixx/vH/v2a5ilcTZewgcwQcA9OrqrAf12nX5OFmrpq5O3QzUzwDhVV2denv99R0f+9WvMv/Wd3uJ+NfNZVexYnsKBIDgAwB6M2aoXWR6TrvynI5/0oJV2ONkZsxQ23LvvdkVn9pakcsvz4SfmhqRIUOyKz5ATBF8ACRPrsrM3r3qxV1P29zYrivPzarI+i0gCjFTzuw1tJWRu3bNXrBQROT00zPB59xzo9GFCfiA4AMgf2GOxfASGtzOTLvmGmdtcbMqsv57t+1ZsEDkzjudtcnra2zapN76sV0FECEEHwD5C3Mshpfp9W5npjU2+ruYn5HT9mgzrqZMCf41Jk5UxwAZZ8ZFbdo/4BLBB0C8eZle73ZmmlnlyI5WidIWBdy1S+TwYevzWbXJrD3Grjina/64ec/amB+jqE37B1wi+ACIt0JMr/fSlWesROkXB/R7QUA/w4j2XkeMYCo5ihLBBwDseOnK0ypR+inscdh0U/9eR40KtSlAEAg+ABAEfSVq6lT11s2mm9rigvr9s6wGbeuPAZATwQcARKK3SvD69eqtvotMY9ZV5oWTXemLbUNaJB7BB0DxcRpijIOCo7RK8I9+JPK972Xfpw3aXrRIZNiwzP3aooNuw4jdey7WDWmRaClFUZSwG1Eozc3NUlpaKk1NTdKzZ8+wmwMgbPX1ImPGiOzYUZgF+py8Xq4ZWtrz3ch30LRZxcduxhzgMz+v31R8AKBQnFSinMzQWrNG3T4i6EHTcdiQFnCJ4AMAheJXd5q2vcTq1WplSFvXhzAC2OoUdgMAAAAKhYoPgGQohhlK+q4y43sB4AjBB0AyRHGGkpcwprUprOATtWn/gEsEHwDJ4GVPr6B52SX+2mvVf+sDUnm5+lhjo/p9kNWqqE37B1xiOjuA5Cr0dHYjp9PFV61ytzs86+mgyDCdHQCKgdPp4lVVaqUnbvt+ARFE8AGAqNMCUn09h+w2VAAADodJREFUU9iBPDGdHQAAJAbBB0ByMUMJSBy6ugAkVxRnKBHGgEARfAAgSrQw1tDQcS2fsKawA0WE6ewAEEVVVdZr/JhhCjuKGNPZAaDYRXHBRaAIEHwAIIqcrvEDwBVmdQEAgMQg+AAAgMQg+AAAgMQg+ABAXLDGD5A3BjcDQFxEccFFIGao+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMQg+AAAgMSIRfB555135Morr5QhQ4ZIt27d5Mtf/rJUVlbKkSNHwm4aAACIkc5hN8CJ3bt3S1tbm6xatUr+5V/+RV599VW5+uqr5dChQ3LHHXeE3TwAABATKUVRlLAb4cWyZctk5cqV8tZbbzl+TnNzs5SWlkpTU5P07NkzwNYBAAC/+Hn9jkXFx0xTU5P06tUr5zEtLS3S0tLS/n1zc3PQzQIAABEWizE+Rn//+9/l3nvvldmzZ+c8bunSpVJaWtr+VV5eXqAWAgCAKAo1+FRVVUkqlcr59eKLL2Y9Z//+/XL22WfLxRdfLFdddVXO8y9cuFCamprav/bu3Rvk2wEAABEX6hifgwcPysGDB3MeM3jwYOnatauIqKFn4sSJ8rWvfU1++ctfSqdO7nIbY3wAAIifohnj07t3b+ndu7ejY/ft2ycTJ06UMWPGyEMPPeQ69AAAAMRicPP+/fuloqJCBg4cKHfccYc0Nja2P9avX78QWwYAAOIkFsFn/fr18re//U3+9re/yYABA7Iei+lsfAAAEIJY9Bd9//vfF0VRTL8AAACcikXwAQAA8APBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJAbBBwAAJEbsgk9LS4uMGjVKUqmU7Ny5M+zmAACAGIld8PnP//xP6d+/f9jNAAAAMdQ57Aa4sW7dOlm/fr08+eSTsm7dOtvjW1papKWlpf37pqYmERFpbm4OrI0AAMBf2nVbUZS8zxWb4PN///d/cvXVV8tTTz0l3bt3d/ScpUuXyq233trh/vLycr+bBwAAAvbBBx9IaWlpXudIKX7Ep4ApiiLnnnuunHXWWbJo0SJ55513ZMiQIfLSSy/JqFGjLJ9nrPh89NFHMmjQINmzZ0/eH1zSNTc3S3l5uezdu1d69uwZdnNii8/RP3yW/uGz9Aefo3+amppk4MCB8uGHH8qxxx6b17lCrfhUVVWZVmT0tm/fLtu2bZPm5mZZuHChq/OXlJRISUlJh/tLS0v5JfRJz549+Sx9wOfoHz5L//BZ+oPP0T+dOuU/NDnU4HPdddfJJZdckvOYwYMHS21trTz//PMdQszYsWPle9/7njz88MNBNhMAABSJUINP7969pXfv3rbH3XPPPVJbW9v+/f79+2Xq1Kny+OOPy9e+9rUgmwgAAIpIuqqqqirsRtgpLS2VPn36tH+l02m5++67ZeHChXLSSSe5Olc6nZaKigrp3Dk247oji8/SH3yO/uGz9A+fpT/4HP3j12cZi8HNRk4HNwMAAOjFMvgAAAB4EbuVmwEAALwi+AAAgMQg+AAAgMQg+AAAgMRIZPB555135Morr5QhQ4ZIt27d5Mtf/rJUVlbKkSNHwm5aLKxYsUKGDBkiXbt2lTFjxsjWrVvDblLsLF26VE4//XTp0aOH9OnTRy644AJ54403wm5W7C1dulRSqZTMnz8/7KbE0r59+2TmzJly3HHHSffu3WXUqFGyY8eOsJsVO0ePHpVFixa1X2NOPPFEqa6ulra2trCbFnlbtmyRadOmSf/+/SWVSslTTz2V9biiKFJVVSX9+/eXbt26SUVFhbz22muuXiORwWf37t3S1tYmq1atktdee01+9rOfyf333y833XRT2E2LvMcff1zmz58vN998s7z00ksyfvx4Oeecc2TPnj1hNy1WNm/eLHPnzpXnn39eNmzYIEePHpUpU6bIoUOHwm5abG3fvl1Wr14tp5xySthNiaUPP/xQzjrrLDnmmGNk3bp18vrrr8udd96Z975ISXT77bfL/fffL/fdd5/s2rVLfvrTn8qyZcvk3nvvDbtpkXfo0CE59dRT5b777jN9/Kc//aksX75c7rvvPtm+fbv069dPvvnNb8rHH3/s/EUUKIqiKD/96U+VIUOGhN2MyPvXf/1XZfbs2Vn3DRs2TLnxxhtDalFxOHDggCIiyubNm8NuSix9/PHHytChQ5UNGzYoEyZMUObNmxd2k2LnhhtuUMaNGxd2M4rCeeedp8yaNSvrvhkzZigzZ84MqUXxJCLK2rVr279va2tT+vXrp/zkJz9pv++zzz5TSktLlfvvv9/xeRNZ8THT1NQkvXr1CrsZkXbkyBHZsWOHTJkyJev+KVOmyLZt20JqVXFoamoSEeF30KO5c+fKeeedJ//2b/8WdlNi6+mnn5axY8fKxRdfLH369JHTTjtNHnjggbCbFUvjxo2TP//5z/LXv/5VRERefvllee655+Tcc88NuWXx9vbbb8v777+fdQ0qKSmRCRMmuLoGsYa2iPz973+Xe++9V+68886wmxJpBw8elNbWVunbt2/W/X379pX3338/pFbFn6IosmDBAhk3bpyMHDky7ObEzm9+8xupr6+X7du3h92UWHvrrbdk5cqVsmDBArnpppvkhRdekB/+8IdSUlIil19+edjNi5UbbrhBmpqaZNiwYZJOp6W1tVVuu+02ufTSS8NuWqxp1xmza9C7777r+DxFVfGpqqqSVCqV8+vFF1/Mes7+/fvl7LPPlosvvliuuuqqkFoeL6lUKut7RVE63AfnrrvuOvnf//1feeyxx8JuSuzs3btX5s2bJ2vWrJGuXbuG3ZxYa2trk9GjR8uSJUvktNNOk2uvvVauvvpqWblyZdhNi53HH39c1qxZI48++qjU19fLww8/LHfccYc8/PDDYTetKOR7DSqqis91110nl1xySc5jBg8e3P7v/fv3y8SJE+XMM8+U1atXB9y6+Ovdu7ek0+kO1Z0DBw50SOBw5vrrr5enn35atmzZIgMGDAi7ObGzY8cOOXDggIwZM6b9vtbWVtmyZYvcd9990tLSIul0OsQWxkdZWZmcfPLJWfcNHz5cnnzyyZBaFF8//vGP5cYbb2y/Hn31q1+Vd999V5YuXSpXXHFFyK2Lr379+omIWvkpKytrv9/tNaiogk/v3r2ld+/ejo7dt2+fTJw4UcaMGSMPPfSQdOpUVMWvQHTp0kXGjBkjGzZskAsvvLD9/g0bNsj06dNDbFn8KIoi119/vaxdu1aeffZZGTJkSNhNiqXJkyfLK6+8knXfD37wAxk2bJjccMMNhB4XzjrrrA5LKvz1r3+VQYMGhdSi+Pr00087XFPS6TTT2fM0ZMgQ6devn2zYsEFOO+00EVHHnm7evFluv/12x+cpquDj1P79+6WiokIGDhwod9xxhzQ2NrY/piVKmFuwYIFcdtllMnbs2PZK2Z49e2T27NlhNy1W5s6dK48++qj87ne/kx49erRX0UpLS6Vbt24hty4+evTo0WFc1Be+8AU57rjjGC/l0n/8x3/I17/+dVmyZIl8+9vflhdeeEFWr15NNdyDadOmyW233SYDBw6UESNGyEsvvSTLly+XWbNmhd20yPvkk0/kb3/7W/v3b7/9tuzcuVN69eolAwcOlPnz58uSJUtk6NChMnToUFmyZIl0795dvvvd7zp/Ed/mncXIQw89pIiI6Rfs/fznP1cGDRqkdOnSRRk9ejRTsD2w+v176KGHwm5a7DGd3btnnnlGGTlypFJSUqIMGzZMWb16ddhNiqXm5mZl3rx5ysCBA5WuXbsqJ554onLzzTcrLS0tYTct8jZt2mT6/8YrrrhCURR1SntlZaXSr18/paSkRPnGN76hvPLKK65eI6UoipJ/RgMAAIg+BrYAAIDEIPgAAIDEIPgAAIDEIPgAAIDEIPgAAIDEIPgAAIDEIPgAAIDEIPgAAIDEIPgAAIDEIPgAiJ3HHntMunbtKvv27Wu/76qrrpJTTjlFmpqaQmwZgKhjywoAsaMoiowaNUrGjx8v9913n9x6663yi1/8Qp5//nk54YQTwm4egAhL5O7sAOItlUrJbbfdJhdddJH0799f7r77btm6dWt76Lnwwgvl2WeflcmTJ8sTTzwRcmsBRAkVHwCxNXr0aHnttddk/fr1MmHChPb7N23aJJ988ok8/PDDBB8AWRjjAyCW/vjHP8ru3bultbVV+vbtm/XYxIkTpUePHiG1DECUEXwAxE59fb1cfPHFsmrVKpk6daosXrw47CYBiAnG+ACIlXfeeUfOO+88ufHGG+Wyyy6Tk08+WU4//XTZsWOHjBkzJuzmAYg4Kj4AYuMf//iHnHPOOXL++efLTTfdJCIiY8aMkWnTpsnNN98ccusAxAEVHwCx0atXL9m1a1eH+3/3u9+F0BoAccSsLgBFZ+rUqVJfXy+HDh2SXr16ydq1a+X0008Pu1kAIoDgAwAAEoMxPgAAIDEIPgAAIDEIPgAAIDEIPgAAIDEIPgAAIDEIPgAAIDEIPgAAIDEIPgAAIDEIPgAAIDEIPgAAIDEIPgAAIDH+P4rQ9p08IAnrAAAAAElFTkSuQmCC", "text/plain": [ "Figure(PyObject
)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# 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, we create a prior and likelihood for the class labels $C_k$, whereas in the case of discriminative classification, we create a a prior and likelihood function for the parameters $w$. Both are proper Bayesian approaches. \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" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### 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)." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### 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*}$$" ] }, { "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)$, evidence $p(D)$ and the predictive distribution $p(y_\\bullet = 1 \\mid x_\\bullet, D)$. Unfortunately, these expressions are not analytically computable. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Many methods have been developed to approximate the integrals for the predictive distribution and evidence. Here, we present the **Laplace approximation**, which is one of the simplest methods with broad applicability to Bayesian calculations." ] }, { "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": { "slideshow": { "slide_type": "fragment" } }, "source": [ "- Note that $\\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", "\n" ] }, { "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)$. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "### Working out the Laplace Approximation \n", "\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 \\log f(z) \\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." ] }, { "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 \\tag{B-4.151}\\\\\n", "\\mu_a &= w^T_{N} x_\\bullet \\tag{B-4.149}\\\\\n", "\\sigma_a^2 &= x^T_\\bullet S_N x_\\bullet \\tag{B-4.150}\n", "\\end{align*}$$\n", " 1. 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", "$$" ] }, { "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 [**precision-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": 5, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAG2CAYAAAB/OYyEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzde3wU1d0/8M9mgYQockciARIqCCIFEfmBSglQQesV5KnaUqug4MPlgaI+XgokJilo8YK1RYHiBbFaL8FbhXoBJUitCmKxAkLlEiA8RIWES8j1/P4YJ5mdzH1ndnZ2P+/XK6+Q3dmZsxt0PpzzPeeEhBACREREREkgxe8GEBEREcUKgw8RERElDQYfIiIiShoMPkRERJQ0GHyIiIgoaTD4EBERUdJg8CEiIqKkweBDRERESYPBh4iIiJIGgw8REREljcAEn9raWsyZMwfZ2dlo2bIlevTogfz8fNTX1/vdNCIiIgqIZn43wKoHH3wQTz75JJ599ln07dsXn332GW655Ra0bt0aM2fO9Lt5REREFACBCT7/+Mc/cM011+CKK64AAGRlZeGFF17AZ5995nPLiIiIKCgCE3wuueQSPPnkk/j666/Rq1cvfPHFF9iwYQMWLVqk+5qqqipUVVU1/FxfX4/vv/8e7du3RygUikWziYiIKEpCCBw7dgxnnXUWUlKirNIRAVFfXy/uueceEQqFRLNmzUQoFBLz5883fE1ubq4AwC9+8Ytf/OIXvxLgq6SkJOo8ERJCCATAiy++iLvuugsLFy5E3759sWXLFsyaNQuPPPIIfv3rX2u+Rt3jU15ejm7duqGkpARnnHFGrJpOFGg1NTXo0KEDAGD37t1o166dp9ebP38+HnzwQdx888147LHHPL0WEQVDRUUFunbtiqNHj6J169ZRnSswwadr16645557MG3atIbHCgsLsXLlSmzfvt3SOSoqKtC6dWuUl5cz+BBZJIRAOByGEAIHDx5ERkaGp9fbuXMnNm3ahO7du2Po0KGeXouIgsHN+3dganxOnjzZZFwvHA5zOjuRx0KhEFJTU3Hq1ClUV1d7fr2ePXuiZ8+enl+HiJJTYILPVVddhd/97nfo1q0b+vbti88//xyPPPIIJk6c6HfTiBKeHHyUQ8dEREEUmODz+OOPY+7cuZg6dSoOHz6Ms846C1OmTMG8efP8bhpRwmvRogUAeB58qqur8dZbbyE7Oxv9+/ePfvYGEZFKYGp83MAaHyJnunbtiv379+PTTz/FoEGDPLvOzp070atXL6Snp+P48eNcdoKIACRpjQ8R+Sc1NRUAPK/xqaysxP/7f/8PaWlpDD1E5AkGHyIyJQcfr4e6fvzjH+Pjjz/29BpElNw4gE5EpmJV40NE5DUGHyIyFaseHyIirzH4EJGpWNX4XHvttRg+fDg2bdrk6XWIKHmxxoeITMWqx+cf//gHDh8+zGnsROQZBh8iMhWrGp/XX38du3fvRq9evTy9DhElLwYfIjIVqx6fIUOGYMiQIZ5eg4iSG/uTichUrGp8iIi8xh4fIjIVix6fLVu2YPfu3ejfvz969Ojh2XWIKLmxx4eITMWixmfFihUYN24c/vSnP3l2DSIiBh8iMhWLHp/MzEwMHjwY/fr18+waREQc6iIiU7Go8Zk9ezZmz57t2fmJiAD2+BCRBVy5mYgSBYMPEZniXl1ElCgYfIjIlNc9Plu2bEH37t0xfvx4T85PRCRjjQ8RmfK6xuebb77Bvn37kJGR4cn5iYhkDD5EZMrrHp+f/vSn2LhxI+rq6jw5PxGRjMGHiEx5XeNzxhlnYOjQoZ6cm4hIiTU+RGSKs7qIKFEw+BCRKa9rfJ599lm89tprOHbsmCfnJyKScaiLiEx52eMjhMDUqVNx8uRJ7Ny5E61atXL9GkREMgYfIjLlZY1PZWUlrrjiCuzZswddu3Z1/fxEREoMPkRkyssen/T0dLz00kuun5eISAtrfIjIVCz26iIiigUGHyIy5XWNDxFRrDD4EJEpL2t8pk2bhqysLDz99NOun5uISI3Bh4hMednjs2vXLuzduxehUMj1cxMRqbG4mYhMeVnjs2LFCvznP//Bj370I9fPTUSkxuBDcSEvDwiHgblzzY8tKADq6qTXUGwoe3yEEK72znTu3BmdO3d27XxEREY41EVxIRwG5s2TQo2RggLpuHA4Nu0iiVzjI4RAbW2tz60hInKOPT4EwP8eF/m68+ZF/qy+7rx5QH6+tXaSe+QeH0Dq9WnevLkr5921axfeeust9OvXD6NGjXLlnERERtjjQwDio8dl7lwp1Gi1g6HHX8rg42adz0cffYTf/OY3WLBggWvnJCIywh4fAhA/PS5a7WDo8V+zZs2QkpKC+vp6V2d2nXXWWRg/fjzOP/98185JFGilpcCSJcCUKUBGht+tSUwiiZSXlwsAory83O+mxK38fCEA6buVx71uR4sWsb0u6UtLSxMAxO7du/1tyMGDQuTmSt/JHD+vYNm0Sfqf3qZNfrckrrh5/+ZQF0XQGm7yssclL097eG3uXKBFC6C6Wvou9/xwJpd/vFzLx5bSUuD++6Xvdl6Tl2fvNYnCyefl5rXtfu7J/LuimGDwoSaU4Sc11dthJr3aooKCxtBTXQ2MGsXZXH5rWMtn0SL9m5LNm1ZdXZ1LrTPh580/KLwIHE5DarS/Ky/Dk1/BjIHQNQw+pEmrx8Wr6xj1MFVVASNHAmvXSt9Z4+Ofhh6fJ580Dj4Wb1rV1dVIT09HdnY2jhw54mZTmyori/xOTSVSOPTyvfj1OSXS78dnLG4mTeoel4KC2BQ0f/CBFHLkHqaCgsbQs3att+0gYw37dbl0vpKSElRXV+PQoUNo06aNS2fV8e230vedO4F//COycJTFpPFD/l0ki9LSpkFm8+bI77Jt22LTpiTA4ENNqGt65J8Bb8OPHHrknh0/2kH6Gnp8XDpfdnY29u/fj9LSUv2VoO3cGAApuBiFl6NHpX81X311ZPCRHwOCHYLc/rxirbQUWLrU71bEzpIl0t89LbfdFtu2JBEGH4qgVchsZaq7G9dV9uykpko9TbFuB+lrqPFx6XwpKSno0qULunTpon+Q3RtDbq4UWtQ3/+3bpe979kjflf96Vg5/KUOQ24EgFj1LTj4vt2YMOAld8utkyt9LEIJatKZMkf6ulZUBRUXAuHFASYn0u1q2DBg4sPHYbduACRP8a2sCYfChBkazt7wOHXV1jdeVQ49WbZH8c6xqYpOa6kaW+sNWFVWAdFMqK2scQgKADh2k/2nLz6s5uWnJNwalzZu1bwzyNYxu/suXS9+VN5DJk+21ySkvQ5XM6ucl/+46dGgaTpz+7pyELiD6oOZlL5fy3HI46dnTnXMrj9+8WerpmjIF6NhReq5ly8hjKyvtX5fDuNpcmF4fGFzHR5/VdXq8Xs+H6/fEkdxc6Zfww9dPAAFA/FXxmK2v3NyI069YsUIsWrRI7Ny501p75PVo1qwxXufk88+FWLky8uumm6TXXHGF9L2goPG5559vPJ+Xa6j4tT6L1nVVv1u7vztNBw82foby17Jl0uuXLWv63MGDTV+zcmXjNfVeo2b3vQwfbn1No2g/J6trKCl/R/KfnV5Xec0EWhPIzfs3e3wIQGSPixEve1xY0xNnVL0HqVOnAv/8p9Tjs2wZ0LWrdo+PUW+MwuLFi/Hxxx8jMzMTZ592mvm/TOUek5Urjdv92mv6vQh/+5v0XfkXatgw6fvatY3Xlodc5H9JJ9q/nJ30pJkx6u0YOLDpObV6auReDaBpj4fe+e28F3m4qLTU2ntSnlt+7Zw5QGGhtc8pml6+lSuBPn3M35P6uuqaNWqCwYcA2Bvm96rGx4/aIjKgutGk/tAFXw1o38iAxm54vecVrrzySmRmZuLcc891dxhIvlktWWKtULa4WPp+112Nj8lDYZMnS+fbts2bYSp1oIpVwLIbUrxgNDQGNK1n0atH8vK9aJ27d2/r545mGYU+fbTPH6vfTwJj8CHf+VlbRNa5PavrtxMnSje/Nm0i/6XvlDI0DBwo3SSnTJGee/tt6S/PjTcCL7wgPTZiBLBunfE5ly51NsvIat2JHKiGDm0MPl7XAcULo54aoGmPhxefh9dBU+4RVfaMmv3d6NpVCtxlZdLPbhR0l5U1/veQ6H+vLGDwIV9Z2Q6D4Sc+mK7jY3cNFqdd8h06SP/6NxtWkJ9fskR6DQCce27j8XffDTz0UOMQxsqVUgBzY0aN3ULfoiJgzBjr57cjI0P78/Kb2Q1dr8fDTXpB0yicyDMElTPQ7IQTq3835MBtdead3Ga5Xdu2Nf6DYsMGaXiuZ8/GMJloM+RsYPAhX8VDbRFZ06THR2uYxmLvSFVVFcK1tfr/AzK68ZSUSDcq5TFmdUGPPy79rF4oUXljVfYu6A0nLFki3YSMbhilpcCxY8CaNY0zdOT2q0OVHKjGjdM/X7QyMvzb5M5p6Jo82Z2bsnrRP3WvmzIkqMOyXjgpLJS+K4Ow2TIK27c3XnPoUClkd+jQ+PfDaX2V8r8B9dCusn16bU7SzQ8ZfMhXftcWkXUN6/jk5BgPy5jdtEpL8efHHsOshx7CRABLlP8ylW8ORvU5etOctXqO5NqKo0cjvwPSv4I7doy8+ck9Q3o3zKVLgZ/8pOkQjLq34JFHgF/+0nqNhjIg+WHLFmDaNPcCh8xp6Ip2SCYjAxg+XL+XTv13aMKEyCBgNAy3cKG0Avi4cY2/N7NlFAoLG8OHTCt4GNXvaIVIszoptXHjgN/+tvF8SYrBh4gsaejxueiixuCjxeymtWQJ9j74IGoBnAZE3py0Qo1cYGz2r2Kt9hQVSd/VFfNA05vRhAmNa/oYDWtZLbpVk0OYHKrKyqTwBQCrVkmPyz0Eq1ZJzysDkZdDE//+N/Dxx8D06YlxQ8zIkGq5tHoNlX+HlMOcI0dGvl7vcxg5ErjzzqaPX3tt4xo/snfeAVasAG66CRg9OvK5vn3tvyf13zOzGW0tW0rvT97rZ/JkFkaDwYfIFXl50s7xVnqlCgqkIbug9TI31PhUOShvVnbJDx2KB1aswMzPPkPKH/4gTQ8GpBAyZ440a0Y5DKC+CQ0c2DTolJZqFw/36yd9f/xx4KKLpDAxfz6wfr10rbFjI1fNBaReHb2pxIDzols5hGmFKnVvgPpnIDmGJtysR9LqiZOnyMvf5Z7Gysqmw6t2g6bRMgorVkhfSrm5wIAB1s+vRauNcsBWLgcg/zdbUtJ0AcQkrPUJVPA5cOAA7r77bqxevRqVlZXo1asXli9fjgsuuMDvppHP/A4e4bC14mtlMXfQNPT4VFU1LaIEjFfLff55afjnBykAGjaqUN7k5T+b3eTtFg+/844UfDp2BC69VAo+Z57Z+K9fubC4tFS69siR+jcDp0W348ZFhiq5x0cZ+LZvb/z5kkua9vgkOi/rkZR/Z9Th043tPLR6X1atavx9jh0b+ZxXv0+tgC3/d+X1tiUBEZjgc+TIEVx88cUYMWIEVq9ejU6dOuE///mP97s6UyD4HTyszDyzMoMtbpWWIvW77wAA1QcPSv+jlGtwrNxEZs8GNm2KfEyrS14OBWY3Ba2bTF4e8Oab2se/+WbT5156SQpDSvKNV9mDBOjvIaXc+kFZqKo+TjZ7dmSo6thRuimNHSuFqc2bI38m90yZIg1FKf+eRbtoo/p49Wvkvze9exv/Pt3s6VIGbMD+ootJIDDB58EHH0TXrl3x9NNPNzyWlZVl+JqqqqqIbvmKigqvmhe3/O4JiZV4CB5GbQh06AGAJUuQ+uSTAIAq+V+URuS6HJnipnDs2DHMmzcP2c2aYQaAkHpmlZUbvtZNZto0Kdwoh6KMwlBxMaDuLZb/9evGrtl6/7r2+0azZYtU06P0zjuR35X69o1+SCYeZGQ0/r1Q/z2zsiigl8sCuNnTJQdw+R8QubnSTDKAix/KXNhCIyb69OkjZs2aJcaPHy86duwoBgwYIJYuXWr4mtzcXIEf9hdSfiXTXl1+7cGVm2v9XPn51rYCsnourffh9R5jRteK5bU9c/CgePSOOwQAccOYMdr7MMn7LK1cabg30eeffy4AiA5t2tjbI8ts3yOtc8j7eq1c2XidOXOkx+bM0d8LyuoeUmvWSM+tWWN9fyqjNsdib6Xhw+3vbZUo1H+HvP685b9/a9Z4c34tWv+dJMCeXW7u1RWY4JOamipSU1PFvffeKzZv3iyefPJJkZaWJp599lnd15w6dUqUl5c3fJWUlCRd8BHC/MbrxY3Zz01P4yF4JOJmq4sXLxYAxNixY6UHHN60v/nmG3HnnXeK//31rxuPt7qZoxGt6xuFITs3o4MHhZg82drNw85Nxo/gY7SJ6003NX3u88+9a4vfvP684yVwxEs7opCUwad58+Zi6NChEY/NmDFDDBkyxPI5knl3dj96QvwIXEJI98+RI60FDzd7m9Tka7do4c35Y2358uUCgLjiiiukB6K9abv9P2Orwcfpda2+zs751YHPjQDohLK3Lpl4HQj8+n2qMfhECEyNT0ZGhrSZoUKfPn3w6quv+tSiYNGqP0m0mhe5nikcljbaDoeB6mqgRYum1ygokI754ANvZlgVFDReu7q6cRmNIIuY1eUGt2sm4nVrBiPq2g4/V1lORl7/nYmX32cQ/9vwUIrfDbDq4osvxo4dOyIe+/rrr9G9e3efWhQ8c+dKN/l584DU1NgU2yqvWVAgPRZN6MnLazyPmnJm18iRUqF2KNQYPGTy9eXQ4/b7V76/qqqm7z+oTNfxsfg/12+//RZ1dXWNNwU3g4/6fPwfPhlx++9gvEqW92mVCz1QMfHJJ5+IZs2aid/97ndi586d4vnnnxfp6elipY2u2WQe6lLyYwjGrZoXq8Nnyi952Cs/P/J55TncGvKKh+Jqr7z++usCgBg8eLD0gMNu/L59+4pmzZqJDz74wP1GWuF0+MHq6+JleMOOZB3qosBIyhofIYR48803xXnnnSdSU1NF7969TWd1qTH4+Ft061bgMgoRymCTnd00/OiFHjc+C79qmmJlzZo1AoAYMGCA43PU19eLdu3aCQBix44dLraOovL559LsrUQuZKZAS9rgE61kDz5+znZyO3BptV0ZeuSgow488mNG53GrPdEcF4/WrVsnAIg+ffpEdZ7a2lpRUlIiampqXGoZESW6pCxupuho1dVYWfTPi2vLP0dzTXXblX9WXycclup9AKnmZ+3axnobvVoju4s51tVZqxeSn5fbEyRR7dWlEA6HkZmZ6UaTiIjscyGIBUay9vj4OQTjdc2LfJ6UFO3zyT1Myq/sbCFycsyHy4LYK+Olzz77TAAQXbp08bspRJRkONTlUDIGn3haSNCra+rVDqkLnY3qfNxuUyLaunWrACA6duzo+BxFRUXiN7/5jXjnnXdcbBkRJToOdZFlfg3BWJmy7sZQm956OcrhNHmdHvWw1wcfJNh+Wh5zYx2fNWvWYOnSpTj99NNx6aWXutU0IiLLGHwSnJ21s9y82ccicOnVDn3wgVTHA0S2QX5cDj9r1wKjRgHvv8/QY4UbNT5XXHEFTj/9dOTk5LjUKiIie0JCCOF3I2KloqICrVu3Rnl5Oc444wy/m0NR0Asqo0Zphx75+JEjpefl70BjEHIr9MgrSFs5l90iaj8dOnQIGT8sgFZfX49QKORzi4goWbh5/w7Mys1EMr3QI29DkZWlf/z770vf5fADSMFDa1sLp+QVpM1WalbOOgsCeagLAGpra31sCRGRcww+FChGQ1Jyr83u3Y3bRIwapT2NPz8fqK9vfK16W4toaG3TYed9xCtl8HEy3HXq1Cns27dP2q6CiMgnrPGhQDEaklIOF82d21jTo3e8coPSefPcXc8o1hu0xoJc4wNIwef000+39frPPvsMw4YNQ8+ePfH111+73TwiIksYfChQ7NTCDBsG5ORoL06oFTzWrvU+/AQ19ABAs2bNkJKSgvr6ekc9PocPH0azZs24eCER+YrFzZRUzIKHF8FEPqc85T6IoUeWnp6OyspK7N69G1nqYioL6urqUFFRgbZt27rfOCJKWCxuJrIoL6+xzsZK6JGH0qwUJ1s1d25j6HGziNoP0a7lEw6HGXqIyFcc6qKEJs+wAozrg7RCkZuLOWotshhEbu3XRUTkF/b4UEJTzrDSW1tHbwNXN9bWUZ67qsr93qRYi6bHZ8qUKbjjjjtw6NAht5tFRGQZe3wo4fk1w0ovUOm1JQjk4FNdXW3rdTU1Nfjzn/+M+vp63HHHHV40jYjIEgYfSgqxnmFldO4ghx+nPT51dXVYtGgR9u7di86dO3vRNCIiSxh8KGkoA0dhoXczrGK1QasfnNb4pKWlYcaMGV40iYjIFgYfSipz5zaGHq9mWMVig1a/uLFDOxGRn1jcTI4op4mbKSjwZxNOrTZqzbCSH9dro933ClgPVG4VUceK0xqfvXv3oqSkhNtVEJHvGHzIkSBsxKluo94MK3k/L702BuG9xorTHp/f/va36NatGx5++GEvmkVEZBmHusgRKzUqfm/PoGyj1r5dyv28Ro6Mrh7H7/caK05rfKqrq9GsWTNkZ2d70SwiIssYfMixIGzEaRRuCgoaH1+71nhhwSC8V6W8PP11i9TkFautDLk57fF56aWXUFtbiyTaIYeI4hSHuigqygUCrW4NEUta4Ubdxvfft7awoFvvNRb1UV4Nzzmt8QGkTU6bN29u+3VERK4SSaS8vFwAEOXl5X43JRByc4XIz7d27MiRQgBCtGghfbf6Oq8p30N+vnEb8/Ol482Yncfq681eZ/U4p693cv5bbrlFABALFixw1igiIgfcvH8z+JAuuzfocLgxEMSK1XCWmyuFMytttBKA5NDj9L16EUrsnMfp+adMmSIAiLy8PMuv+eijj8R1110nHn30UXsXIyL6gZv3bw51kS6toR01eahk5EipTkQ9TdxrVod0iouloS6zNloZ+tGbEm+H0Wfr5lCh20ORTmp8tmzZgldffRXr1q2zdzEiIi+4EMQCgz0+9sk9JUa9BvLzI0c29sBY6U2wM5Rm1AtjtfdE/srP136NlXarj3F7OMqtnh6960Q7FHnXXXcJAOKOO+6w/Jovv/xSPPbYY6KoqMjZRYko6XGoyyEGH/vU4UZ9g1aGHrs3cDdrXcyGdLQCj96fnV4j2vDjdX1UtMNzQggxZ84cAUBMnz7dvYYREZlg8HGIwccZdciRb6B6oUf9umjCj51Qodd7YtSzo3eMl+3U6unSCyVWC67NuBWu8vPzBQBx2223Rd8oIiKLGHwcYvBxTl3ALH/XCz3q10Xbk2J1WEx+TUqKcaBRBwE3Qo3T4/RCiVcFztGc98EHHxQAxE033WT5NR999JEoKSkRdXV19i9IRCQYfBxj8ImOHHLk0JOdbf1Gb9ZrYXZzthsq5PBjdLyyl0WrjXLYshK6lMfZmRJvNoQYb7O6Fi1aJACIG264wdLxFRUVAgD/uyOiqDD4OMTg45z6hqzs8ZFFW6xsNhxjp4jZbEjHytCP2z09asqhQiuP2+Xm8JzsiSeeEADE2LFjLR3/zTffiOzsbHHmmWdavwgRkQqDj0MMPs6ob5ByWJDDj9NeGa3jzApwrRYxWznWSru9CA/K13nV4+NVaFu+fLkAIH72s5/Zag+HuYgoGgw+DjH42GdWj2JnNpfR1Hitc+fk2GtTNCHHSfiJNvR4WePj1lIBaitXrhQAxKhRo+w3iojIIQYfhxIx+Hh1g5OPt1J3YzX8GA3hKF9jFGTUx2sVJ6s/E6NeFOVnomyD3jm03qPdz9WPWV1uefnllwUAMWzYML+bQkRJhMHHoUQMPl7WoWjtc6XX6yEvXqh3PbPwoQ496j9rUYaenBzjtipDlzLY6LUlK0v/OWXPjPyYXu+UFVbqjeLFG2+8IQCIwYMHWzr+jjvuEOPHjxcbN270uGVElMgYfBxKxOAjhHd1KNGe38oQjlnQMTu33rR69evkYbacHPMamvz8xuPUx6hngpn1TJmxOhQXL9asWSMAiP79+1s6vl+/fgKAWL16tcctI6JExuDjUKIGHyHcr0Ox+3q947SGcJS1PurQox4GU4YQdTiQe1msDrkZ9ToZvSd1G5U9TV58tl6Hn2iGR9etWycAiN69e1t6/erVq8WiRYtEaWmp7XYSEckYfBxK5OAjhDe9B9HcJM16fOThJCuF0urhKb2eIbtba1h5T+qQo3zM6ZRzr3vpnFzbrK5JCCE2btwoAIi2bXvEVe0RESU2Bh+HEj34CBE/9SJGIUzZ4yMPVck9QlpBJTdXu5dI75rqhRXV11IOf1n5fJTXVfb2WAlRWtcw6ukxqpNyk9mwot61N23aJKQFCbvE7XAcESUeBh+HkiH4CBH9ZpRuLUSod8OXA4McUNSBRP18fr556JHJx7VpExlU1AFL3R6znd/VX0bDaWbhxSwMqYvKjXpWzJ43YhR+9D7nadO2CgAiPb2D6fn37NkjiouLxaFDh5w1kIjoBww+DiVD8HGjxyeaoRCzMKHulZG/y4FFHv5SPi6EdrGx8qav7BFSFjGrQ4+6bVam2BsFH/V7dlKz43SI0o0eIb1r67+PrwUA0apVK9Nzy/t6/eIXv3DeQCIiweDjWKIHHzdrfJwMhZgFAJm6R0c9FKUeStKaPq6+hjrAqEOL/Lh6WrtRe9UzwNRtVB7vdLNWu49bfd4OrbCs93dp9uy9AoBITU01Pe+jjz4qsrKyRC6LgYgoSgw+DiVy8HF6A7V7TqMeAbOVmZWvl3td5J/lQBEKRZ5fa9hOL+SoC5iVgUXvcaMhK/Wx6sBmFALlc1hZ28itsBQNrc9ZKxAdOnRIyJuO1tfXu9cAIiIDDD4OJWrw8bJ3wP5QiPZ1jEKROqCkpUX+bNTjY6U+SNlmZU+PXqhTrw1kdH67n4XdoUg3e/HMrqHVJnUg+v777xuCT1VVlXuNICIywODjUCIGn1jUg9gZCrHS06P1uFaPj1avivIa8mvC4chzqeuDtGp19Hpq9HqelOTwo3eM2Wdht/jcbliyw+j3qHXdEydONASfiooK9xpCRGSAwcehRAw+ZjOw3JohZKjebBAAACAASURBVHUoxIjZcJIcOuQZWcrn1HtyqV+j7E1S9i4p26ZVn6N8z1q1RnrvSR2Q1IspmgVAuyEm2pl6Rm3R65XS+l3l5dU2BJ9vv/1W99ylpaXiggsuEDfccAOHxIgoagw+DiVi8DHjRo+Q/FxKivlQiPI1RrO69Opx1LO9lEFEvU2EVo2Pktb2EvK5srIij1XXDalXi9Z6D/K5srP19/4y+ozt/m7c7PGx8vvW68lKSQkLAOLAgQO659+wYYMAILp37x59Y4ko6TH4OJSMwUeI6GqAtIY+jIZCrFxPb1aXUXGyMmRYGV5St035nDrk6E1p13of6sf0ao2sfsZehSUjdn/fTd9/SwFAfPPNN7rX+O6778Rrr70mXnrpJecNJSL6AYOPQ8kafIRwduNVP6dePdlo2MpoUcDc3MYAIw8VyT1E6vPJvTPKa2ZnS88PH64dSpTDXsrzKHdQV9fp6G09YRYEhGhaa2T2OVp93mlYMqM3PKoXdNS/x7S0NgKAmDFjm7MGEBHZxOAjhJg/f74AIGbOnGn5NckcfISw13tg1Nuh1+uiNaPKqB1y6NAKDOohLfnPyqCk1UajdXrUx8rnCoWsfW5GIURZa2T2ORpdw2rPjhs9P2pWV+w+88wzBQBx++1fuHdxIiIDSR98PvnkE5GVlSV+/OMfM/jYZKVexOimaraWjXpGlRatvbqUxyuvrzXcpdWjk58fueCgevFB9bFWe3yUbVX2GGl9Tlrv3cn2H1o9bdH00LitW7duAoD45JNPdI959913xYYNG8SxY8e8bQwRJYWkDj7Hjh0TPXv2FO+++64YPny4b8En2v2s/GQ2Q8jsRquuwVEHEq31erSGsbRqbuTncnK0e1rMhrOUw2Za7VD3JCmvb/aetR5Xt8ust8sKK1uDKJ+T/2550Quk5eyzzxYARHFxse4x2dnZAoBYv369t40hoqSQ1MHnpptuErNmzRJCCNPgc+rUKVFeXt7wVVJS4toH53QIw29Wenz0aIUXZX2LOghorVpstoihPJXdaFNS5bnka6h7eMyuo57VZWfYz+7PbvBj6EtP3759BQDx3nvvaT5fX18vLr30UtG9e3exf/9+7xtERAkvaYPPCy+8IM477zxRWVkphDAPPrm5uUJec0T55dZQVzzdjKxw+wattSWEleP1jlMXPGu1TQ5fRjVI6mnmVq+vdU719dSBSvkar3tf9M4Z679nAwcOFADE22+/HZsLElHSS8rgs2/fPtGpUyexZcuWhsf87PGRxcvNyIzb7TTr8dE73mwHc6PQo26vVtvlnix17Y48M8wsnCmH2vSup14/KCdHf1hNPcQZ7bBnLHqXzAwdOlQAEKtWrYrdRYkoqSVl8Fm1apUAIMLhcMMXABEKhUQ4HBa1tbWm5/CquDkebkZG3O6ZUocYOWwY1coY1fioz6ecuWXWXq1eFvWaP1pr+Zi9N2Uhs1b4UQ4VmvUkaZ0nGtEMV7ph+PDhAoB48cUXY3thIkpaSRl8KioqxNatWyO+Bg0aJCZMmCC2bt1q6Rxezury+2akx+1aJLMhH7Pd2dXn0QtNdsOa+vrqwmvl6s1OAqDyOa0Voa3uSu/W3wsvtrCwavTo0QKAWLFihebzv//978WgQYPE0qVLY9wyIkpUSRl8tPg5q0uLnzcjJWXvitnsM3WPiXpRQPWx6nBhp4BZTf689IbJrAYVZcGy8nitVZWVr9Frux75GOU5rYQ1t0OP3yH7qquuEgDEsmXLNJ+fMGGCACAeeOCB2DaMiBIWg88P4in4+H0z0mqLk14eo9cqp4ib1eoYrdysPE5v4T+r70VvY1G9HiD10JSV96S+lvrL6PPzKvT4Oaw6fvx4AUD88Y9/1Hz+66+/FqtWrRLbt2+PXaOIKKEx+DiUTDU+Tod19Na/Ub7GbFsKuz01ToeJ1OFJvYmp3jCY1lCYumBZfS2thRv1hs+8CsGx6lEy84tf/EIAEA8//HBsLkhESY/BxyEvgk+83Iy0OGmb0XCVnd4R5RCa3rWNCp6tnE/ZHnWtkFYYkdf6UdfqaK1BpDVkpg496t4l5fXcHvaMpj7JbbfccosAIBYsWOD9xYiIBIOPY24Hn3i6Gelx0hulF37s1vCYtcXu82bHa60ArfW+lb1CZsXV8oKKWvuTaZ1TueK0Wz0+0QxdukUZUm+//XYBQOTl5TU5rqysTIwf/7y49Vb97SyIiOxi8HHIzQ8uHm5GVjm5EatDjpUp61bO59bnpXWcephO730rQ4rWe1UHKb16Hq3Xqgus3fj9x8P2KMr3MXPmTAFA3HvvvU2O+9Wv1ggAolOnvu43goiSFoOPQ8m8V5eToRd1DY3W7Curn4McSqx8Dmafl52eI+X7tlO/pH7PyjogNbNd6eMh/LpBfh8XX3yXACBmz56t8fxakZU1XEyaNMmnVhJRImLwcShZd2ePZujFbFsKt3tyrEy/V/bo6AUkreEss53j5dfICygqV5E26/Gxeu7ECD9zBAAxffp01ePBf39EFJ8YfBxKxuDjpMZH/VqzbSncrN1R9rxo0do7y6iIWh7OMuuVkSlreeShK3VdkF7Nk9Up/EEPByNHFggA4oILbhNCJM77IqL4xeDjULIFH70bkpMCZ7Man2iupWZ1M1E7bbEyI01vh3h10FH3ngVt2DNav//974W04e9NcbFuFRElPgYfh5Ip+ETTC6MXLOwGjmh6AvTCj5PQY+X1ej1Ceu/Jbr1UIlm0aJEAIFJSrm/yOZx77rli8ODBYs+ePf41kIgSDoOPQ8kSfKKpu7EaLJTPa63Ho9cTYKfHQx1+og09eu9Bvamp0XmVwSdZezqeeOKJH3p8xkb8no8dO/bD4xBHjx71u5lElEDcvH83AyWcujogPx+YO9f4OPn5urrI144cCaxdC+Tk6L9m3jzpuLo6IByWfl67Fli3DigsBKqrgRYtIttQUCAdl59v7X28/z4wapR03mbNItum9f603ndentQ+5WPK9/3hh43voXt36Ut9XuV7Bhrbv3Zt42Nmn3Ui+fvfUwEAPXtW4euvG3+vdXVp+Pzzz7Fv3z60bt3a51YSEelwIYgFRjJPZ7dKq3fDaDhM3QuirglyY9hLOaXc7jmsvAetDUf1jpWLpa2cPxFJ7/d5AUCMGjVK9XjyfA5EFFsc6nLIrQ/OaD0YNacL/vlFHejMgoO8DUR+vv2hKSvUQ0vhsP1zKNuo9b7UwU0dUt2ctRZk8vu8/vqXBQBxySWXaD6f6J8DEcUeg49Dbn1wVm/oZrOTgsLKjC2zomEnn4FeLY7dc6mnouv1amm1NZp6qUSifH9vvPGGACAGDx4cccy7774rxo//iwB2J+znQET+YPBxyIstK8ymRgc99MiMZmzpPRfNNhduBSn1uj9mAUi95YSVIU3lNcyGNIM07Kmk/Bz+/ve/CwCif//+Ecdce+21AoC44oo/BvI9ElH8YvBxyKtNSt3s5YhnWjO29Ho61Ftk2OkRcTNUavXo6L0HK7VNZtdw47h498EHHwgAonfv3hGPz5s3T/zkJz8Ra9eu9allRJSoGHwc8mI6u1vDMV7wogBbHWi0rqEVkKxew4thRGVhsjr8qEOP8jV2ei2SqQ5o48aNAoDIzs72uylElCQYfBzyah0fNwpwveB2T4ReoDE6l/Jnq8NGcjhRz6BS01tLyKjt8lCW+kvvvbgRftwMPfEwm3DTpk0CgOjSpYv7Jyci0sDg45CXCxhGM+XaS271RCh7YvRep66ZkW+6VntxZFqzsIzaJYcBs9coe3uUG5Bq9Ro5DStGwc8N8TCs9uWXXwoAokOHDg2P1dfXu38hIqIfMPg45FXw0RrmCkL40Xtcb0q73hYWOTnGxc5G57DaVifvU2tFab0v5TT2aEODlZ6xaPg9rLZz504BQLRq1arhsVdeeUV06tRJTJo0yZuLElFSY/BxyIvg48XaNV6w0xNhFGK0nrNaGGz22bhxw9Zqn7KnR/lnrefcCg3qWii3xWJYTc/evXsFANFC8eYWLlwoAIgbb7zRuwsTUdJi8HHI7eDjdCNNv+j1RBgVKGv13MhDUfJz8s9W3q/e0JSbN2zlufRCj/wZWK33cXJ9r3cu93pYTc+hQ4cEIO3JJQ9xVVRUiM2bN4t///vf3l6ciJISg49Dbn5wZrOK4jX8aM3K0mujutZG+Z7VISGatXq8CAjyudW1V1ptVh6XleXOdWMVRmIVspSOHDnSEHyqqqq8vyARJT0GH4e8WLnZ6nHxsKCb1k3S6vCTvHO5cgdzN266Xg4JyedWhh+9IUnlcdHW9sR6+MnrYTW1kydPNgSfioqK2FyUiJIag49Dbu7VZfUmJg/t+M1KrY5e+JHDjjwTSvl8NDfdWPT4KIe0tDZQVdcpOQ2pfhUc+9HjU1tb2xB8ysrKRH19vXjggQfECy+8IE6dOuV9A4go6TD4OOTldPZ4ZqUnQi/8yD/LoUcZcJzedLWG14zaaDeI6BVg5+dHBjWj45xez43jnF43VjU+QggRDocFAHHgwAFRWloqAIhQKMShLyLyBIOPQ8kYfOz0RKjDjzz0I/f4GG1Vof7ZqFdMayacWwXPeu1QBhtlL5Be4bada/rRA+jXsJosPT1dABDffPON2L9/v7jpppvEdddd5+1FiShpMfg4lGzBx0lPhHroRz28ZdQzonUevWPUhdLq80dTLK31GmW7tYqzrZ7HLi9CkV/Dakpt27YVAMS2bdu8uwgR0Q8YfBxKtuDj9KarngmlLOI2GxIyCj/qn/VmxjnZ5NVK6MnPjzy31ddEw+1hML+G1dTOPPNMAUBs2bLFmwsQESm4ef9uBkpYeXnWj507V/o+ahRQVweEQlK8yc4G1q4FCgqk5+fNA/LzG/+sfK3yz/Jx+fmNx9XVST/PnSudb+1aYOTIxvOrHx82zHr7ledWKihobMvcudJxgHRuZVv13od8vFNG19BroxG996mUlweEw9JxZu0vKJCOsfN3BQBSU1MBANXV1aitrUWzZvxfCREFhAtBLDCSrcfHLrk3RD11XblpqLomRm8mlLIHyaznx6uZSfHSO2J0DS+uHYv33bNnTwFAFBcXi5/+9KeiU6dO4q233nLWYCIiExzqcojBR59yCEg560qvFkf9Oqs3Wb1w48VaNPG27EAsZ2F5XQd03nnnCQDivffeEz169BAAxIcffui8wUREBhh8HGLw0WZ16w2rU9D16IUbr3p84lEs36uXvUwDBw4UAMTf/vY38f3334tNmzaJ48ePR9dgIiIdDD4OMfg0pQw3Zs8re4Ts3jz1bvix7AWJF7Fcadmrz3fo0KECgFi1alX0jSQiMsHiZnKNWbGssjhXLjgOh60X4wJNi3flnz/4QCpkVp7HSjFwkBUUANXVQIsW0ne5qNsrys+zsFC6ptXfmxG5uLmqqirKFhIRxZgLQSww2OPjnNMhGr0eBqubvCZSz4+fvVtu9zKNGTNGABBz584V8+fPF+vWrXPnxEREGjjU5RCDT3Ts3jytLmLo50J8sRLLWV1613azruiqq64SAMSVV14pAIjbbrst+pMSEengUBfF3IgR1odo5LV4PvhAe1hFvZ6P3rCWW2vp+M1onR6vh/b0hhmjvZY81JWRkYFf/epXGD58uAutJSLyHoMPmSookEIMAMyZI33Xu3nKN9acHP1aEuVieWbhJug1PlYWJ/Qq/Ghd261rycGnd+/emD17dpQtJSKKHQYfMqS8eQLaKzLLN087KxArBT3cGLGy0jLgfu+W171MLG4moqBi8CFdejdPrfCjfDyRg4xdTrYNiVYseplatGgBADh+/LjDVhIR+YPBhzTp3TyN9uJi6IkPsehlknt85s+fjxdffBG7du1CKBSyfyIiohhj8CFNRjdPZfj54R/+yMlh6IkXsehlkoMPAKSlpTH0EFFgpPjdAIpPeXnGN8W5cxtneLVoAaxbF7OmURyQg8+kSZPw+uuv+9waIiLrGHzIEa0ViL2Wl2f9OgUF9no+yB65xicUCuHss8/2uTVERNYx+JBtyvqfqqrGOh+vw4+8VYbZdeT2hcPetSXZQ5jc41NdXe1zS4iI7GGND9ni5dowZqxcx+6U+rw8KSBZ3XOsrq7xNVber3o5ADucti0W5OCzZcsW7N+/H5mZmbG5MBFRtFxYSTowuGVFdMy2V4jVFhNubv9g9TVax3n9eUTTNq89+eSTAoAAIDZt2hS7CxNRUuJeXQ4x+DgXbzdhNzf8jCbAeL0HV7yETbXly5cLAOKss84SR44cie3FiSjpMPg4xODjXG6uvR3Zc3O9bE3jddzafDOaAOP1rut+bnCq5/nnnxcAxKhRo2J/cSJKOr4Hn5MnT4r9+/c3efzLL7+MukFeYvBJPHZ3jDcSTYDxYgd0t9rmhVdeeUUAEJdccok/DSCipOJr8Hn55ZdFZmam+PGPfyz69esnPv7444bnzj///Kgb5CUGn8TiRdiI5pxuhjC32+a2F198UQAQF154oX+NIKKk4eb92/Z09sLCQmzevBlffPEFnnrqKUycOBF/+ctf5ELpqIutjSxYsAAXXnghWrVqhU6dOuHaa6/Fjh07PL0mxSevptSrF2a0OkMtFusaOW2bF1566SUAwN69e/1rBBGRA7aDT01NDTp27AgAGDRoENavX48lS5YgPz/f82XrP/zwQ0ybNg0ff/wx3n33XdTW1mL06NE4ceKEp9el+KI3pd6N8OMkwMRqXSM/Fo3Uc/ToUQDgVhVEFDx2u4hycnLEF198EfFYVVWVuOGGG0Q4HI66C8qOw4cPCwDiww8/tHQ8h7qCz8tZTk7qaGJVeBxvNT4bN24UAET37t39aQARJRVfanwqKiqEEEKUlJSI0tJSzWM2bNgQdYPs2LlzpwAgtm7dqvn8qVOnRHl5ecNXSUkJg0+AeTml3kmAidVU83ic1bV58+aG6exERF7zJfj0799fN/D4ob6+Xlx11VWGs0pyc3MbFllTfjH4BJNXU+qdBJhYrWsUr+v4fPnllwKAaN++fWwvTERJyZfgM2nSJNGtWzexbdu2iMc3b94sLr/88qgbYtfUqVNF9+7dRUlJie4x7PEhM04DTCzWNYq3RSNlO3fuFDfeeKMAIFq1ahWbixJRUnMz+Fjeq+vPf/4z7r//flxyySV47bXX0KlTJ8yZMwevvvoqrr76aveKjiyYMWMG3njjDaxfv95wj6DU1NSGPYWItNTVWdvXS36+rk76bmdPLKezr5y2zWv/+te/8MILLwAAqqqqYnNRIiKXhISwNwd9wYIFyM/PR11dHcaMGYP7778fAwcO9Kp9EYQQmDFjBlatWoUPPvgAPXv2tPX6iooKtG7dGuXl5TjjjDM8aiVRYtu6dSueeuopLFq0CABQX1/P2V1E5Ck379+Wp7OXlpbif/7nf1BQUIBzzz0XzZs3xw033BCz0AMA06ZNw8qVK/GXv/wFrVq1wqFDh3Do0CFUVlbGrA0UX/LyrE/rLiiI3e7liaxfv37Izc1t+Lm6utrH1hAR2WM5+PTo0QPFxcV4+eWXsWnTJhQVFWHq1Kl48MEHvWxfhCeeeALl5eXIyclBRkZGw9df//rXmLWB4ks4bG3NHHmtnXA4Nu1KdMohZAYfIgoSyzU+Tz/9NG644YaGn8eMGYN169bhyiuvxN69e7F48WJPGqhkc1SOkoBc3zJvXuTPSloLHpJze/bsQadOnRp+rqqqQqtWrXxsERGRdZaDjzL0yAYOHIiNGzfiZz/7mauNIrLDKPww9LhLCIG+ffvi5MmTCIfDqKurY4EzEQWK5eCjJysrCx999JEbbSFyTCv8MPS478iRIwiFQgiFQmjRogUqKysZfIgoUKIOPgDQtm1bN05DFBVl+CkslPazYuhxV7t27XDs2DF8//336NmzJyorK1njQ0SBYnuTUqJ4Fk87mCeqUCiE9u3bNxQ4s8eHiIKEwYcSSjztYJ7oGHyIKIgYfChhKGt6qqqk71amupN1f/zjHzF16lRs2LABLVq0AMDgQx4rLZUW4Cot9bsllCAYfCghaBUyz53L8OO2N954A0888QR27drV0OPDGh8LePN2rrQUuP9+fnbkGleKm4n8ZDR7y8o6P2TdlClTMHjwYAwePJhDXXbIN++rrwYyMvxuDVFSY/ChQLMyZZ3hxz3XXXcdrrvuOgCs8bGlrCzyOxH5hsGHAi1edzBPBqzxseHbbyO/E5FvGHwo0OxsOsqenugcOXIEZWVl6N69O1JTU1njQ+4rLW1ay7N5c+R3pYwMDh2SbQw+RGTJm2++iV//+tcYNWoU3nvvPQ516dG6eW/f3vhdfQPnzbvRkiVSLZSW225r+lhurr1//RCBwYeILDp27BjS09ORnZ0NgDU+uoxu3oWF0pcSb96NpkyRCsCVNm+WQs+yZcDAgZHPMTCSAww+RGTJtGnTMHXq1IahLdb46Lj2WqBnz8jH3nkHWLECuOkmYPToyOf69o1d2+KdUe/XwIFNgw+RAww+RGRZKBRq6OlhjY+O117T7/FZsUL6UsrNBQYM8L5dRASAwYeIHOJQlw6t4ZpVq6QhrjlzgLFjI5/jcA1RTDH4EJGpmpoaXHPNNcjKysLDDz+Mli1bMvjoUQ/XlJYCX30l/bl3bw7XEPmMW1YQkal9+/Zh9erVePrpp5GWlgaANT6WlZYCRUV+tyK4MjKk4UD2jJFL2ONDRKbatm2LZcuWoaKiAqFQCEBsanzq6upQXFyM0tJSZGRkYNiwYQiHw55dz7LSUmn21pQp9m7IHTp416ZElZHBWW/kKgYfIjLVrl073HrrrRGPeT3UVVRUhJkzZ2L//v0Nj2VmZuKxxx7DuHHjPLmmZU733urY0bs2EZElHOoiIke8DD5FRUUYP358ROgBgAMHDmD8+PEoCuLQ0eTJHK4higPs8SEiU1u3bkVaWhq6d+/eUNvjVY1PXV0dZs6cCSFEk+eEEAiFQpg1axauueaa+Bj2UlKu2lxWJtX2yGv6XHhh01WduWozUcwx+BCRqcmTJ+Pjjz/Gyy+/jPHjxwPwrsanuLi4SU+PkhACJSUlKC4uRk5OjqvXjhq3XCCKeww+RGQqLS0tYrsKwLuhrlL1PldRHhc1OxtnDh0KrFkj1fJs2wZMmCCt3VNYyC0XiOIEgw8RmVq3bl2ToSevgk+GxTBg9bioRduL07u39D0Rt1xwOruNyEcMPkRkiTyNXeZVjc+wYcOQmZmJAwcOaNb5hEIhZGZmYtiwYa5eV1e8bZwZbdhwM6w4nd1G5CMGHyJyxKsan3A4jMceewzjx49HKBSKCD9y+Fq0aFHsCputbJypHA6T/ywPg23fLn3fts3aOc1EGzYYVijJcTo7ERn661//iiuuuALLli2LeNzL6ezjxo3DK6+8gi5dukQ8npmZiVdeecX/dXzUliwBLrgg8kseBisslL5PmND43JIl/rU1WqWl0lBerGqsiFzGHh8iMvTZZ5/h7bffRq9evSIe93oBw3HjxuGaa66Jz5Wb1YyGw2bOBB57DHj8ceCii6TngtzTouwxIgogBh8iMjRhwgT06tULffv2jXg8Fnt1hcPh2E1Zj6b2xWjo6sILpe8XXRTs4mZ5CE8estu2DaislP6snt0GcI0iilsMPkRkqH///ujfv3+Tx2OxV1dM2al9icXGmaWlwEMPAaNHR251oTeVXm6Xemd4vQUV7YYV9ey2CRMa/8w1iihAGHyIyBGvh7riWiw2ziwtBR55RPrSYiVsuLmgojycJ69PtHKl1OPj1+w2IocYfIhIV2VlJTZs2IDs7GycffbZEc8ldfDRo+xhkYeE5FldToeDVq4E+vRp/NnOVHpl7ZGTBRW1Fm+Uh7cqK4GWLaU/y985vEUBwOBDRLp27NiB0aNHo0OHDigrK4t4LhY1PoGj1cMiz+pyOhzUp492bZCVBRG1goidBRWt9hjJw14c3qIAYPAhIl2VlZXo06cPOnXq1OQ5ucenpqamYfPQwLCzDQVgvScj3hY7jJbZ+2nZsnHYq0+fxvfDFZ0pjjH4EJGuoUOH4quvvtJ8Tg4+gFTgrPw57nm1maiVxQ612AliyoUQ9Ridz86CimbvR6buleIiiRTHGHyIyBFl0KmqqgpW8Im3nhm7QSya8ykXVJSZBTtlD4485FlWBvz4x97PbiNyGYMPETki1/gAAazzcdoz4xU7QaysDHjnHeOwYXS+hQuBnTuBceMap8mbBRdlD86330qPffttbGa3EbmMwYeIdI0bNw7V1dVYuHAh+ihnFgFISUlBs2bNUFtbmzhr+bjBSX2L3SA2Zozz840cCdx5p7V2OSUPp23fHuxFGykhMfgQkSYhBN555x2cOHECDz/8sOYxqampqK2tDV6Pj5fU9S2xWOwwljp0iPyuVU+0dm3jd3kWmYxT3slnDD5EpEkIgaKiIuzevRtZWVmax6SmpuLEiRMMPka0hoOCPOtJHh6TvxvVEy1fLn0pcco7+YzBh4g0paSkYPTo0YbHJNRaPrHsmYnnWU/qbS6+/bZx6GrVqsbjVq2SZob16iVNZwekXqCOHYHFi6XAM2kSMHVq5Pnj7f1S0mHwISLHEmq/rngr1HU7iFk9n5UZYeo/y+TeHLmHMCvL/RqfIPeWUVxg8CEiTV988QXKyspw3nnnoXPnzprHJP22FV4thCgfqwxi0d7wrQY75YwwZY9PYaG03QXQ+OfevYGjR4EjR6TH27QBnn8e+Phj6eePP5Z+VurbFxgwwH77ZfHcW0aBwOBDRJqeeOIJLFmyBHPmzEFBQYHmMUkffLxaCFFLrG74WuFs82Yp7IwdK/0s/3ngQCAnB/jwQ+1z/e1v0pfS8OHABx+43Woiyxh8iEhTx44d0bt3b5xzzjm6xyRUjY8TVtbfkXtNAKkGZvNm93qF4sGiRcC//x352AsvSIHniiuAqZUt+wAAIABJREFUG2+MfK5v39i1jUgDgw8RaSooKNDt6ZElVI2PE1bW38nLi12vkBfkITY9AwY0Hbo6ckQKPpddBvzyl962j8gmBh8icsy3oa4gFbjG2/YYdpWWAkuX2ntN27aR36O5tlc1VJS0GHyIyDFfg09QClydbI8Rjzf8yZMbrxGraf+xrKGipMHgQ0RN/POf/8SkSZNw0UUXYanBv/YTpsYn3nqQ/Lzh630Wyp9jFS6C3ltGcSlwwWfx4sVYuHAhSktL0bdvXyxatAjDhg3zu1lECWXnzp3497//jU6dOhkelzA1Pm72ILmx/o6fN3w3Pou+faXZW9EWMsfbZrKUEAIVfP76179i1qxZWLx4MS6++GIsWbIEl19+Ob766it069bNtevU1dWhpqbGtfMRGWnevDnC4bDfzYgwZswYrFmzpiHY6En66exa3FgI0eoNX9k7Ew35PNdeK63GDDR+dzLENmAAp6xT3ApU8HnkkUcwadIk3HrrrQCARYsW4e9//zueeOIJLFiwIOrzCyFw6NAhHD16NOpzEdnRpk0bdO7cGaFQyO+mAJCmso8x2wEcMQo+8VjvEi23VmV2q6dKPo+ykHnChMhjWFNDCSIwwae6uhqbNm3CPffcE/H46NGjsXHjRs3XVFVVRfwPuaKiwvAacujp1KkT0tPT4+YmRIlLCIGTJ0/i8OHDAICMeL9hq8SkxicRC1zjbXsM2bhxwE9+IoWelSuBPn28HWIzq62Kt9orSgiBCT7ffvst6urqcOaZZ0Y8fuaZZ+LQoUOar1mwYAHu1/sfpkpdXV1D6Gnfvn3U7SWyqmXLlgCAw4cPo1OnTnEx7PXCCy+gU6dOuOiiixrapyUmNT5u17skYg+SU/JnIQ9rffstUFkp/Vn+Lv/+W7Z0/7Mw67HSej6Wm8lSQgpM8JGpe2GEELo9M/feey9mz57d8HNFRQW6du2qeaxc05Oenu5SS4msk//e1dTU+B58ampqMGHCBNTX1+PAgQOWgo+nPT5uF7gGtQfJixu++rNQDm+pP4sJE+Ljs4jX3jIKjMAEnw4dOiAcDjfp3Tl8+HCTXiBZamqqaXGmGoe3yA/x9Pfu2LFjGD16NA4ePKi7OakskMXNQZ4iffXVkT1WTnqqlK8fOlQa0jLahPTIEWDGDOm4kSPdf09EMRaY4NOiRQtccMEFePfddzFW3igPwLvvvotrrrnGx5YRJZZ27dph9erVlo6Va3yqjx6V/hUehFqMoE6Rdqunyug8hYVN/yx/Hn36xP/vlsiCwAQfAJg9ezZ+9atfYdCgQRg6dCiWLl2Kffv24fbbb/e7aURJqaHHp7w8OCspB5VbPVVm52nZsrG4GZD+rFy12Smz2irlZq6AtKFrSUnkcUqJXHtFngpU8Ln++uvx3XffIT8/H6WlpTjvvPPw9ttvo3v37r61KS8PCIeBuXPNjy0oAOrqODxNiaMh+MR6AcNkLHB1q6fK7DyyPn0a/+xGT57dHiuz5+Oh3ogCKVDBBwCmTp2KqVOn+t2MBuEwMG+e9Gej8FNQIB2Xnx+bdhE5dfvtt+Ojjz5CXl4errvuOsNjfQ0+sbzpJdO06rIy6fv27cCHH7p3XrOepq5dtXt8glB7RYESuOATb+SwYxR+lKHHSs9QvOJ2Iclh69at+PLLL1FXV2d6bEONT22t183yllkPUjQLBQY1NH3/vf1d2Y046bGSh7jiufaKAofBxwVG4SdRQk+stgsh/61YsQK7du3CgAEDtA9Q1GqkHjwIAKj67jvpuaDWYnjZgxSkneTl9XwAQLmCfVB/r0RaRBIpLy8XAER5eXmT5yorK8VXX30lKisrHZ8/P18IQPqu9bNXunTpIv70pz9FPPbRRx+Jli1bij179rhyjcGDB4vbb7894rHevXuLe+65x5XzJzs3/v7FTG6u9BcbEK8CAoC4+IefNb9yc/1usTUHD0ptPXiw6XObNknvZdMm++d1+lqj9tg5xso1Pv9ciOHD9X+HXv1ezT6baD53SihG92+7UnzOXQll7lypZ2fePCA1NXY9PUOGDMGnn37a8LMQArNmzcKsWbOaFH7Pnz8fp59+uuFXcXFxxGvk7UJGjx4d8bjRdiGUwKZMATZtAjZtQuqiRQCAKvnv2bJlDc81fEW7gaYTpaVSD456FpHZa+T9quKBlfbIPVVOe17ka9TXAy+8IP2+5NlcBQWNx8XL75XIBRzqctncudLyF9XVQIsWsRneGjJkCJ555pmGn5977jns27cP9957b5Njb7/9dvz85z83PF+XLl0ifnayXQgF09dff43i4mL069cPgwcP1j5IMcTR4ochruofan3iphYjmuGlsrKmQzuJuqWFXMi8cSPQtq305+3bpe87dzYed+iQNAzWoQPQsaN779mstioZZ++R5xh8XFZQ0Bh6qquln2PR43P33Xfj+PHjSElJwX333YfCwkK0atWqybHt2rVDu3btHF3HznYhFEzvv/8+pk6diquvvhqvv/666fG+zeryUlGRflGv2bRqK2vVFBVJm4G6GSCcKiqSvs+Y0fS5FSsa/6z+n5hbU8nNaqu4PQV5gMHHRepCZvlnwNvwM2jQIITDYWzevBnvvfce2rdvj4kTJ2oeO3/+fMyfP9/wfKtXr46YreVkuxAKps6dO+Oyyy7D0KFDLR3fEHx+2OsuIYwb13QYx+pCgVbXqpGDld9r0YwbJ7Xl8ccje3wKC4GbbmoMPwUFQHZ2ZI8PUUAx+LhEa/aWlanubkhLS0P//v1RVFSEpUuX4s0330RKinb5lpOhLm4XkjzGjh0b8Ts2E9geH6OemZIS6eauJG9ubDaUZ2dVZOUWELHYMV7rGvLKyGlpkQsWAsCFFzYGn5/9LD6GMIlcwODjAqMp67EKP0OGDMEf/vAHXHnllRg1apTucU6HurhdCGlpWMenrs6/WgwnocHuKsKTJ1tri51VkZU/223P7NnAww9ba5PTa6xbJ313Y7sKojjC4BMlK+v0xCL8DBgwAM2aNcPChQvdPznic7sQcp/duq2IoS6/hmycbN5pd9+rsjJ3F/NTs9qebdukHiPVDEtPrjFihFQDpF54MagLMhL9gMEnSnV11qasy89bWAzXkeeffx5Tp07FOeec480FEH/bhZC7jh07hi5duiArKwuffPIJ0tLSTF/TEHyqqrxunj4nm3faXUVYq+fIjNwTJS8KuG0bUFmpfz69Nmm1Rz0UZyWM2H3Pcs2PWpAWZCTSwOATJTv/yHW7p6e+vh5lZWVYvnw5duzYgVWrVrl7AUoqu3fvxrFjx3DgwAFLoQdoDD41NTWor6/XrS3zlFubd5pdw+5QnronasKExj+7vemmm2FEfq99+3IqOSUkBp8AW79+PUaOHInevXujqKgIrVu39rtJFGB9+vTBtm3b8J28/YQFco0PIIUfOQglHCfTquWeKOUU9iBsuql8r3rblhAFGINPgOXk5KC+vt7vZlCCaN68OXr37m3rNcqgU1VVlbjBxwllT9SYMdJ3O5tuyosLKvfP0ivaVh5DRIYYfIjIMWWPj691Pm6It1WC33lH+q4cIpNpDZU5YWVXeq+n2RPFGIMPEQGQCuSrqqowevRoZGZmWnpNSkoKmjdvjpqamvgKPlZDjLooOJ5WCb7zTuCXv4x8TC7anjMHUPbOyYsO2g0jZu/ZyYy5ePoMiTQw+BARAODBBx/E1q1b8fbbb1sOPoDU61NTU4PqeFrE0GqIiccZSsowpjccVlio/bjbYcTJjDmiOMfgQ0QAgEsvvRQZGRno1auXrdelpqbixIkT8dXjE6+s9ERZCWMrV0rbR3hdNB2LGXNEMcbgQ0QAgIftrgT8g7hYyyco3BpOk7eXWLpU6pWR1/VhGCEy5cOiG0SUSBh8iChI2ONDRKitrUU4HLa1XYWsYb+ueKrx0ZIIM5SUQ2Xq90JEljD4EBEeffRRFBQUYOrUqXjggQdsvTYwPT7xOEPJSRiT2+RX8Im3af9ENjH4eIGb+FHAyNtVhMNh268NTPCJxxlKTnaJnzJF+rMyIHXtKj1XVib97GVvVbxN+yeyicHHC/E4RZbIwMMPP4wZM2agVatWtl8bmOATjzOUrIaxJUukQmb5S0kZkOTnuJ4OkS4GH7Jk/fr1WLhwITZt2oTS0lKsWrUK1157rd/NIpe0bNkSfeSZQjYFpsYnHlkNY3l5UkgK2r5fRHGIwYcsOXHiBPr3749bbrkF1113nd/NoTgSmB6fIJMD0ubNnMJOFCVOZ08AmZmZWLx4ccRjGzduRHp6Ovbu3evKNS6//HIUFhZi3LhxrpyP4sd3332HefPm4bnnnnP0egYfIgoS9vhEKw6myA4ZMgSffvppw89CCMyaNQuzZs1C9+7dI46dP38+5s+fb3i+1atXY9iwYa62keLXjh07UFBQgO7du+NXv/qV7dcHOvhwhhJR0mHwiVYcTJEdMmQInnnmmYafn3vuOezbtw/33ntvk2Nvv/12/PznPzc8X5cuXVxtH8W3Nm3aYPLkyTjjjDMcvT7QNT7xOEOJYYzIUww+0YqDKbJDhgzB3XffjePHjyMlJQX33XcfCgsLNWfotGvXDu3atXO9DRRc5557LpYsWeL49YHu8YlHchgrLW3aa+zXFHaiBMLgE604mCI7aNAghMNhbN68Ge+99x7at2+PiRMnah7LoS5yG4OPR6z2JnMKO5EtDD4JIC0tDf3790dRURGWLl2KN998Eykp2nXrHOoitRMnTiA9Pd3RdhUAg49n4qA3mSgRMfgkiCFDhuAPf/gDrrzySowaNUr3OKdDXcePH8euXbsaft69eze2bNmCdu3aoVu3bo7aTPGhX79++O6777Bu3ToMdNBDGegan3gWB73JRImIwSdBDBgwAM2aNcPChQs9Of9nn32GESNGNPw8e/ZsAMCvf/3riMJqCpa6ujqUlJSgtrYWnTp1cnQO9vgQUZAw+HjBh1kZzz//PKZOnYpzzjnHk/Pn5ORACOHJuck/4XAY5eXl2LNnD8466yxH52DwIaIgYfDxQoymyNbX16OsrAzLly/Hjh07sGrVKs+vSYknPT0d5557ruPXM/gQUZAw+ATY+vXrMXLkSPTu3RtFRUVo3bq1302iJMQanxjiGj9EUWPwCbCcnBzU19f73QwKsFWrVuGLL77AmDFjMHToUEfnYI9PDMXjgotEAcPgQ5TEioqKsHLlSqSlpTH4EFFSYPAhSmJjxoxBy5YtHYcegMGHiIKFwYcoiU2YMAETJkyI6hys8SGiINFe3peIyCL2+BBRkDD4ECWpqqoqHD16NOrzMPgQUZAw+BAlqQ8//BBt27aNqr4HYPAhomBh8CFKUgcPHgQAdOjQIarzsMaHiIKExc1ESermm2/Gf/3Xf+H48eNRnYc9PkQUJAw+REnstNNOw2mnnRbVORh8iChIGHxcVldXh+LiYpSWliIjIwPDhg1DOBz2u1lEnmHwIaIgYY2Pi4qKipCVlYURI0bgF7/4BUaMGIGsrCwUFRX53bSoPfDAA+jbty/S09PRq1cv/OUvf/G7SRSladOmIS8vD0eOHInqPKzxIaIgYfBxSVFREcaPH4/9+/dHPH7gwAGMHz8+8OGnuLgYjz76KL788ktMmDABN910E7755hu/m0UOHT9+HIsXL8b999+PlJTo/jcg9/jU1NRw7zgiinsMPi6oq6vDzJkzIYRo8pz82KxZs1BXV+fJ9TMzM7F48eKIxzZu3Ij09HTs3bvXlWv87W9/w+jRo9GjRw9Mnz4ddXV1DbOCKHjq6+vxu9/9DtOnT0fr1q2jOpccfAD2+hBR/AtE8NmzZw8mTZqE7OxstGzZEj/60Y+Qm5sbN/+TLS4ubtLToySEQElJCYqLiz25/pAhQ/Dpp59GXG/WrFmYNWsWunfvHnHs/Pnzcfrppxt+GbVTCIE77rgD5513HgYPHuzJ+yHvnXHGGbjvvvvw+OOPR30uZfBhnQ8RxbtAFDdv374d9fX1WLJkCc4++2x8+eWXuO2223DixAk89NBDfjcPpaWlrh5n15AhQ/DMM880/Pzcc89h3759uPfee5sce/vtt+PnP/+54fm6dOmi+9ytt96KjRs3Yu3atQ21HZTcmjdv3vDnePnHCBGRnkAEn8suuwyXXXZZw889evTAjh078MQTT8RF8MnIyHD1OLuGDBmCu+++G8ePH0dKSgruu+8+FBYWolWrVk2ObdeuHdq1a+foOv/617/w1FNPYfv27YbhiOLfwYMHkZ6ejjZt2kR9rpSUFDRv3hw1NTXs8SGiuBeIoS4t5eXlpjfwqqoqVFRURHx5YdiwYcjMzEQoFNJ8PhQKoWvXrhg2bJgn1x80aBDC4TA2b96MBx54AO3bt8fEiRM1j41mqGv37t0AgHPOOceT90GxM23aNLRt2xZPPvmkK+fjlHYiCopA9Pio/ec//8Hjjz+Ohx9+2PC4BQsW4P777/e8PeFwGI899hjGjx+PUCgUUeQsh6FFixZ5tp5PWloa+vfvj6KiIixduhRvvvmm7kydaIa6hg8fHlFLRMFVXl4OAOjatasr50tNTcXx48cZfIgo/gkf5ebmCgCGX59++mnEaw4cOCDOPvtsMWnSJNPznzp1SpSXlzd8lZSUCACivLy8ybGVlZXiq6++EpWVlY7fz6uvvioyMzMj2t+1a1fx6quvOj6nVdOnTxehUEhcddVVnl2jqKhInHPOOZ6dP5m58ffPrmPHjolTp065cq6MjAwBQHz++eeunI+ISKm8vFz3/m2Xrz0+06dPxw033GB4TFZWVsOfDx48iBEjRmDo0KFYunSp6flTU1MjZpx4bdy4cbjmmmt8Wbl5wIABaNasGRYuXOjZNcrLy7Fjxw7Pzk+xdfrpp7t2Lg51EVFQ+Bp8OnToYHln6AMHDmDEiBG44IIL8PTTT0e96JpXwuEwcnJyYn7d559/HlOnTvW0/ubmm2/GzTff7Nn5KbgYfIgoKAJR43Pw4EHk5OSgW7dueOihh1BWVtbwXOfOnX1smb/q6+tRVlaG5cuXY8eOHVi1apXfTaIAWL9+PZ566ikMHz4ct9xyiyvnZPAhoqAIRPB55513sGvXLuzatQuZmZkRzwmN1ZKTxfr16zFy5Ej07t0bRUVFUa/AS8nhk08+wbPPPovq6mrXgg/36yKioAhE8OEQi7acnBzujUS2DR8+HIWFhTj33HNdOyd7fIgoKAIRfIjIPRdeeCEuvPBCV8/J4ENEQRGfFcJEFCgMPkQUFAw+RElECIEtW7Y0LGDoFtb4EFFQcKiLKIn83//9H84//3ykpKSgsrLStY1m2eNDREHB4EOURMrKytChQwe0bNnStdADMPgQUXAw+BAlkX79+qGsrAynTp1y9bwMPkQUFKzxIUpCaWlprp6PNT5EFBQMPkQUNfb4EFFQMPgkkVAohNdee83XNmRmZuKPf/yjr23QU1tbi1AohLfeesvvpnjmzjvvxC233IJ//etfrp6XwYeIgoLBJ+BuvvlmhEIhhEIhNG/eHGeeeSYuvfRSPPXUU01WdS4tLcXll1/uU0spHqxatQrPPPMMjh496up5GXyIKChY3JwALrvsMjz99NOoq6vD//3f/2HNmjWYOXMmXnnlFbzxxhto1kz6NXu1oWtdXR1CoRBSUpijq6urXZ0t5bb58+dj586d6NOnj6vnZY0PEQUF71QGhBA4ceJEzL/sbryampqKzp07o0uXLhg4cCDuu+8+vP7661i9ejWeeeaZhuOUQ13V1dWYPn06MjIykJaWhqysLCxYsKDh2KNHj2Ly5Mk488wzkZaWhvPOO69hCOiZZ55BmzZt8NZbb+Hcc89Famoq9u7da7m95eXluP7663HaaaehS5cuWLx4ccTze/bswdVXX43TTjsNrVu3xg033ICysrKG5ydMmIDx48dHvGb69On46U9/2vDzJZdcgt/85je444470LZtW2RkZKCgoCDiNTt27MCwYcOQlpaGvn37Yu3atU3aeuedd6Jnz55o2bIlevTogdzcXNTW1jY8P2fOHAwaNAjLli1DdnY2TjvtNDz11FPo1KkTampqIs51zTXXYOLEiZY/Jy9cf/31mDNnDjp27OjqednjQ0RBwR4fAydPnsTpp58e8+seP34cp512WlTnGDlyJPr374+ioiLceuutTZ7/wx/+gDfeeAMvvfQSunXrhpKSEpSUlAAA6uvrcfnll+PYsWNYuXIlfvSjH+Grr75COBxueP3JkyexYMEC/PnPf0b79u3RqVMny237/e9/j9/+9rcoKCjA22+/jRkzZqBPnz4YMWIE6uvrcfXVV6Ndu3YoLi5GdXU1/vu//xs33ngj3nvvPVufwVNPPYW77roLn3zyCTZs2ICJEyfikksuabjOtddeiy5duuCf//wnjhw5glmzZjU5R+vWrbFixQpkZGTgiy++wG233YbWrVtj9uzZDcds374dr732GoqKipCSkoKzzz4bs2bNwltvvYWxY8cCAA4fPoy3335bM1wlAgYfIgoKBp8E1rt3b90i1n379qFnz5645JJLEAqF0L1794bn3nvvPXzyySfYtm0bevXqBQDo0aNHxOtramqwePFi9O/f33a7fvKTn+B///d/AQC9evXChg0b8Oijj2LEiBH4+9//jm3btmHPnj3o0qULAODZZ59F//798fnnn+P888+3fJ2BAwdizpw5AICePXvi8ccfx/vvv48RI/5/e/ceV1O+/w/8tbtHFxRdxqRNKkm68nANIbdm3EWIZDjhuDzGnZEZcs0MMprG0TExCA+XQWhIDKcpKRqVjhRG+sY0ooaaav3+cNo/W6HYWa16PR+P/XjYa3/W+rxbttbLZ33WWr1x6tQpZGZm4ty5czAzMwMArFq1Cl5eXkrbWL58ueLPlpaWSE1NRWRkpFLwKS0tRUREBJo1a6ZY5u3tjfDwcEXwiYiIgFwuR48ePWqyq1Tq3r17yM/PR+vWraGvr6/SbTP4EJFUMPi8QaNGjVBYWChKv6ogCAJkMlmVn02aNAn9+vWDjY0NBgwYgCFDhqB///4AgOTkZLRs2VIReqqipaUFBweHd6qrS5culd6HhoYCANLS0mBpaakIPQDg4OAAPT09pKWl1Sj4vFqfmZkZ8vLylPqpCD1V1QUAkZGR2Lx5M27duoWioiKUlpYqBRwAkMvllZZNnToVXbt2RW5uLkxNTREeHo7JkydXu/ba8K9//QsrV66Ev78/vv/+e5Vum3N8iEgqGHzeQCaTvfcpJzGlpaVBLpdX+ZmzszOysrIQFRWFn3/+GaNHj0bfvn1x8OBB6OrqvnXburq6rw1V76JiW28KaxXL1dTUKs2DenU+DQBoampWWr/iSrfqzKO6dOkSxo0bh1WrVqFv374wNDTE7t27sW3bNqV2VX1H3NzcYGdnh4iICLi7uyM9PR2+vr5v7bO2GRkZvfY78T444kNEUsHgU0+dO3cOKSkpmDt37mvbGBgYYMyYMRgzZgxGjhyJAQMGID8/Hw4ODvj999+RkZHxxlGfdxUXF1fpva2tLQDAzs4OWVlZyMnJgbm5OQDg+vXrKCwsVFyJ1Lx5c9y6dUtpG8nJyTU6fWNnZ4fs7GzFiExVdf3yyy9o06YNFi1apFhWk0nc/v7+2L59O27fvg1PT0/FzyOWwMBABAYGVrrNgSow+BCRVPCqrnqguLgYubm5uH//Pq5evYqgoCB8+umnGDJkCCZOnFjlOl9//TX27duH9PR0ZGRk4MCBAzA1NUWTJk3g7u6Onj17YsSIEYiOjlaMDJ06deqNdRw8eBD29vZvrTc2NhbBwcHIyMjAli1bcPjwYcyePRsA4OnpiXbt2sHHxwdJSUmIi4vDpEmT4OHhAUdHRwAvJm7HxcVhz549+O9//4tly5YhPT29RvvM09MTrVu3xsSJE3H9+nXExsYqzecBACsrK2RlZSEyMhKZmZn45ptvcOzYsWr3MX78eGRlZWHnzp2iX831stq47QCDDxFJBYNPPXDq1CmYmZnB0tISAwYMQExMDLZs2YKjR48qXYn1Mj09Paxbtw6urq5wc3NDdnY2Tp48qTgoHjp0CG5ubhg7dizs7OywYMEClJWVvbGOx48fVyuALFiwAP/5z3/g5OSEtWvX4ptvvoGHhweAFwflY8eOQU9PD927d4enpyesra2xd+9exfqDBw/GkiVLMG/ePHTq1AnFxcXw8fGp7u4CAKirq+PIkSMoKiqCm5sbpk2bhqCgIKU2w4cPx6xZsxAQEAAnJyfEx8dj6dKl1e6jadOmGDZsGAwMDCpNmq5vOMeHiKRCJtT0pjES9uTJExgaGqKgoAAGBgZKnz1//hxZWVmQy+Uqf4AjNVy9e/eGk5MTNm3a9MZ2tf39y87Ohp+fH+zs7GrlkSEXL15Ez549YW1tjZs3b6p8+0TUsL3p+F1TnONDVAvy8/MRFRWFixcvqvwKqndx69YtxMTEIDc3t1a2z1NdRCQVDD5EtcDBwQGFhYXYuHEjrKysxC4H9vb2iIiIeO2pz/fF4ENEUsHgQ1QLfv/9d7FLUGJqaorx48fX2vY5x4eIpIKTm4novXHEh4ikgsGHqAG4cOECUlJSai2YMPgQkVQw+BA1AKNGjYKDgwNu3LhRK9uvCD6lpaW1coNEIiJV4Rwfonru77//xscff4zS0tJaeVwF8P/n+AAv5vnwlhBEVFcx+BDVc5qamrhy5Uqt9lEx4gO8ON3F4ENEdRVPdRHRe3t5xIfzfIioLmPwIZWTyWQ4cuSI2GVUKTs7GzKZDMnJyWKXUq/IZDJF+GHwIaK6jMGnHsjNzcXs2bNhZWUFHR0dmJiYoHv37ggNDcVff/0ldnkksnXr1sHDwwP79++v1X54Lx8ikgLO8ZG427dvo1u3bmjSpAmCgoLQoUMHlJaWIiMjAzt37oS5uTlAo7sMAAATy0lEQVQ++eQTscusV0pKSpRO7dR18fHxOHfuHIYOHVqr/Whra6OwsJAjPkRUp3HEpxqKiopQVFSEl5/nWlJSgqKiokq/5CvavnxJ799//42ioiI8f/78rW1rKiAgABoaGrhy5QpGjx6Ndu3aoUOHDhgxYgROnDih9FTwgoICfPbZZ2jRogUMDAzQp08fXLt2TfF5YGAgHB0dERERAUtLSxgaGsLb2xtPnz6tcV0PHjzAwIEDoaurC7lcjgMHDih9npKSgj59+kBXVxdGRkb47LPPUFhYqPi8V69emDNnjtI6Q4cOxaRJkxTvLS0tERQUBD8/P+jr68PCwgJhYWFK68THx8PJyQk6OjpwdXVFUlKS0udlZWWYMmUK5HI5dHV1YWNjg82bNyu1mTRpEoYOHYo1a9bA3Nwc1tbW+PLLL9GhQ4dKP7eLiwu++OKLGu2r2rZs2TL88MMP6N+/f632w3v5EJEUMPhUg56eHvT09PDo0SPFsg0bNkBPTw8zZ85UatuiRQvo6enh7t27imXbtm2Dnp4epkyZotTW0tISenp6SEtLe6e6/vjjD5w5cwYzZsxA48aNq2wjk8kAAIIgYPDgwcjNzcXJkyeRmJgIZ2dneHh4ID8/X9E+MzMTR44cwfHjx3H8+HHExsZi7dq1Na5t+fLlGDFiBK5du4bx48dj7Nixip/zr7/+woABA9C0aVMkJCTgwIED+Pnnnyvty+oIDg5WBJqAgAD84x//QHp6OoAXwXLIkCGwsbFBYmIiAgMD8fnnnyutX15ejpYtWyIyMhKpqan44osvsGTJEkRGRiq1O3v2LNLS0hAdHY3jx4/Dz88PqampSEhIULS5fv06kpKSlMJZXeDk5IQJEybAxsamVvth8CEiKWDwkbBbt25BEIRKBzRjY2NFWFu4cCEAICYmBikpKThw4ABcXV3Rtm1bbNy4EU2aNMHBgwcV65aXl+Pf//437O3t0aNHD0yYMAFnz56tcW2jRo2Cv78/rK2t8dVXX8HV1RVbt24FAOzZswfPnj3DDz/8AHt7e/Tp0wchISGIiIjA//3f/9Won0GDBiEgIABWVlZYuHAhjI2Ncf78eUU/ZWVl2LlzJ9q3b48hQ4Zg/vz5Sutrampi5cqVcHNzg1wuh4+PDyZNmlQp+DRu3Bg7duxA+/btYW9vj5YtW8LT0xPh4eGKNuHh4XB3d0fr1q1rvL/qA87xISIp4Byfaqg4BdOoUSPFsvnz52POnDnQ0FDehXl5eQAAXV1dxbIZM2Zg6tSplZ6MnZ2dXantu6gY1akQHx+P8vJy+Pj4KP73nZiYiMLCQhgZGSm1ffbsGTIzMxXvLS0toa+vr3hvZmam+JlqokuXLpXeV1xJlZaWho4dOyqNUnXr1g3l5eW4efMmTExMqt2Pg4OD4s8ymQympqaKeiv6efnv7dW6ACA0NBQ7duzAnTt38OzZM5SUlMDR0VGpTYcOHSrN65k6dSr8/PywadMmqKurY8+ePQgODq527R9CXl4efv31V1hZWaFdu3a12hdHfIhIChh8qqGq00haWlpVTnCtqq2mpiY0NTWr1bYmrKysIJPJFKd2KlSMOLwcqMrLy2FmZqYYDXlZkyZNlGp9mUwmU9kjCF4+7fZqWHu1jZqamtKcKuDFXKlXvaneV9evSmRkJObOnYvg4GB06dIF+vr62LBhA3799VeldlX9XXl5eUFbWxuHDx+GtrY2iouLMWLEiLf2+SFdunQJw4cPR6dOnSr9TKrG4ENEUsBTXRJmZGSEfv36ISQkBEVFRW9s6+zsjNzcXGhoaMDKykrpZWxsrPLa4uLiKr23tbUFANjZ2SE5OVmp5kuXLkFNTQ3W1tYAgObNm+PBgweKz8vKyvDbb7/VqAY7Oztcu3YNz549e21dFy9eRNeuXREQEAAnJydYWVkpjYC9iYaGBnx9fREeHo7w8HB4e3srjS7VBRoaGnB2dkbHjh1rvS8GHyKSAgYfifv2229RWloKV1dX7N+/H2lpabh58yZ2796N9PR0xem1vn37okuXLhg6dChOnz6N7OxsXL58GcuWLavR4wxCQkLg4eHx1nYHDhzAzp07kZGRgRUrViA+Pl4xednHxwc6Ojrw9fXFb7/9hpiYGMyaNQsTJkxQnObq06cPTpw4gRMnTiA9PR0BAQF4/PhxjfbNuHHjoKamhilTpiA1NRUnT57Exo0bldpYWVnhypUrOH36NDIyMrB8+XKlCctv4+/vj3PnziEqKgp+fn41qu9D8PLyQmJiYqWr3WoD5/gQkRQw+EhcmzZtkJSUhL59+2Lx4sXo2LGjYiLx559/jq+++grAi1NAJ0+eRM+ePeHn5wdra2t4e3sjOzu7RnNqHj16VK0RkZUrV2Lfvn1wcHDArl27sGfPHtjZ2QF4MVfq9OnTyM/Ph5ubG0aOHAkPDw+EhIQo1vfz84Ovry8mTpwId3d3yOVy9O7du0b7Rk9PDz/99BNSU1Ph5OSEpUuXYt26dUptpk+fjuHDh2PMmDHo3Lkz/vjjDwQEBFS7j7Zt26Jr166wsbFB586da1RffcMRHyKSAplQnYkQ9cSTJ09gaGiIgoICGBgYKH32/PlzZGVlQS6X8wGLVG2CIMDW1hbTpk3DvHnz3nk79eH7N2zYMBw5cgShoaGYNm2a2OUQUT3ypuN3TXHEh+gd5eXlYdOmTbh//z4mT54sdjmVCIIAe3t7eHh44OHDh7XeH0d8iEgKeFUX0TsyMTGBsbExwsLC0LRpU7HLqeThw4e4ceMGUlNT3/t/SNXBOT5EJAUMPkTvqK6fJTY0NMSFCxfw4MEDxWhMbeKIDxFJAYMPUT2lra2NHj16fND+AAYfIqrbOMfnFaq6WR9RTdSH7x2DDxFJAUd8/kdLSwtqamrIyclB8+bNoaWl9dq7CxOpiiAIKCkpwcOHD6Gmplbl3cDfVUxMDAoLC+Hm5gZTU1OVbfd1OMeHiKSAwed/1NTUIJfL8eDBA+Tk5IhdDjUwjRo1goWFBdTUVDcIu2bNGkRHR2Pnzp0f5KozjvgQkRQw+LxES0sLFhYWKC0tRVlZmdjlUAOhrq4ODQ0NlY8w2tra4tGjR7CxsVHpdl+HwYeIpIDB5xUymey1DxUlkpItW7Z80P4YfIhICji5mYhUgnN8iEgKJBd8iouL4ejoCJlMhuTkZLHLIaL/4YgPEUmB5ILPggULYG5uLnYZRHXa3r17YWVlhQULFnywPhl8iEgKJDXHJyoqCmfOnMGhQ4cQFRX11vbFxcVKv4QLCgoAvHjYGVF9lpKSgszMTOTm5n6w73vFBQFFRUX8N0ZEKlXxO0Uld8wXJCI3N1f46KOPhISEBCErK0sAICQlJb1xnRUrVggA+OKLL7744ouvevDKzMx87zwhE4Q6/sAhAIIgYNCgQejWrRuWLVuG7OxsyOVyJCUlwdHR8bXrvTri8/jxY7Rq1Qp3796FoaHhhyi93nry5Ak+/vhj3Lt374M8ALO+4n5UHe5L1eG+VA3uR9UpKCiAhYUF/vzzTzRp0uS9tiXqqa7AwECsXLnyjW0SEhJw+fJlPHnyBIsXL67R9rW1tat8OKOhoSG/hCpiYGDAfakC3I+qw32pOtyXqsH9qDqquMmrqMFn5syZ8Pb2fmMbS0tLrFq1CnFxcZVCjKurK3x8fLBr167aLJOIiIjqCVGDj7GxMYyNjd/absuWLVi1apXifU5ODjw9PbF//3507ty5NkskIiKiekQ9MDAwUOwi3sbQ0BAtWrRQvNTV1bF582YsXrwY1tbWNdqWuro6evXqBQ0NSV3QVidxX6oG96PqcF+qDvelanA/qo6q9qUkJje/qrqTm4mIiIheJsngQ0RERPQuJHfnZiIiIqJ3xeBDREREDQaDDxERETUYDD5ERETUYDTI4JOdnY0pU6ZALpdDV1cXbdq0wYoVK1BSUiJ2aZLw7bffQi6XQ0dHBy4uLrh48aLYJUnOmjVr4ObmBn19fbRo0QJDhw7FzZs3xS5L8tasWQOZTIY5c+aIXYok3b9/H+PHj4eRkREaNWoER0dHJCYmil2W5JSWlmLZsmWKY0zr1q3x5Zdfory8XOzS6rwLFy7Ay8sL5ubmkMlkOHLkiNLngiAgMDAQ5ubm0NXVRa9evXDjxo0a9dEgg096ejrKy8vx3Xff4caNG/j6668RGhqKJUuWiF1anbd//37MmTMHS5cuRVJSEnr06IGBAwfi7t27YpcmKbGxsZgxYwbi4uIQHR2N0tJS9O/fH0VFRWKXJlkJCQkICwuDg4OD2KVI0p9//olu3bpBU1MTUVFRSE1NRXBw8Hs/F6khWrduHUJDQxESEoK0tDSsX78eGzZswNatW8Uurc4rKipCx44dERISUuXn69evx6ZNmxASEoKEhASYmpqiX79+ePr0afU7ee/HnNYT69evF+Ryudhl1HmdOnUSpk+frrTM1tZWWLRokUgV1Q95eXkCACE2NlbsUiTp6dOnQtu2bYXo6GjB3d1dmD17ttglSc7ChQuF7t27i11GvTB48GDBz89Padnw4cOF8ePHi1SRNAEQDh8+rHhfXl4umJqaCmvXrlUse/78uWBoaCiEhoZWe7sNcsSnKgUFBWjWrJnYZdRpJSUlSExMRP/+/ZWW9+/fH5cvXxapqvqhoKAAAPgdfEczZszA4MGD0bdvX7FLkaxjx47B1dUVo0aNQosWLeDk5ITvv/9e7LIkqXv37jh79iwyMjIAANeuXcMvv/yCQYMGiVyZtGVlZSE3N1fpGKStrQ13d/caHYN4D20AmZmZ2Lp1K4KDg8UupU579OgRysrKYGJiorTcxMQEubm5IlUlfYIgYN68eejevTvs7e3FLkdy9u3bh6tXryIhIUHsUiTt9u3b2L59O+bNm4clS5YgPj4e//znP6GtrY2JEyeKXZ6kLFy4EAUFBbC1tYW6ujrKysqwevVqjB07VuzSJK3iOFPVMejOnTvV3k69GvEJDAyETCZ74+vKlStK6+Tk5GDAgAEYNWoU/P39RapcWmQymdJ7QRAqLaPqmzlzJq5fv469e/eKXYrk3Lt3D7Nnz8bu3buho6MjdjmSVl5eDmdnZwQFBcHJyQnTpk3D1KlTsX37drFLk5z9+/dj9+7d+PHHH3H16lXs2rULGzduxK5du8QurV5432NQvRrxmTlzJry9vd/YxtLSUvHnnJwc9O7dG126dEFYWFgtVyd9xsbGUFdXrzS6k5eXVymBU/XMmjULx44dw4ULF9CyZUuxy5GcxMRE5OXlwcXFRbGsrKwMFy5cQEhICIqLi6Guri5ihdJhZmYGOzs7pWXt2rXDoUOHRKpIuubPn49FixYpjkcdOnTAnTt3sGbNGvj6+opcnXSZmpoCeDHyY2Zmplhe02NQvQo+xsbGMDY2rlbb+/fvo3fv3nBxcUF4eDjU1OrV4Fet0NLSgouLC6KjozFs2DDF8ujoaHz66aciViY9giBg1qxZOHz4MM6fPw+5XC52SZLk4eGBlJQUpWWTJ0+Gra0tFi5cyNBTA926dat0S4WMjAy0atVKpIqk66+//qp0TFFXV+fl7O9JLpfD1NQU0dHRcHJyAvBi7mlsbCzWrVtX7e3Uq+BTXTk5OejVqxcsLCywceNGPHz4UPFZRaKkqs2bNw8TJkyAq6urYqTs7t27mD59utilScqMGTPw448/4ujRo9DX11eMohkaGkJXV1fk6qRDX1+/0ryoxo0bw8jIiPOlamju3Lno2rUrgoKCMHr0aMTHxyMsLIyj4e/Ay8sLq1evhoWFBdq3b4+kpCRs2rQJfn5+YpdW5xUWFuLWrVuK91lZWUhOTkazZs1gYWGBOXPmICgoCG3btkXbtm0RFBSERo0aYdy4cdXvRGXXnUlIeHi4AKDKF73dtm3bhFatWglaWlqCs7MzL8F+B6/7/oWHh4tdmuTxcvZ399NPPwn29vaCtra2YGtrK4SFhYldkiQ9efJEmD17tmBhYSHo6OgIrVu3FpYuXSoUFxeLXVqdFxMTU+XvRl9fX0EQXlzSvmLFCsHU1FTQ1tYWevbsKaSkpNSoD5kgCML7ZzQiIiKiuo8TW4iIiKjBYPAhIiKiBoPBh4iIiBoMBh8iIiJqMBh8iIiIqMFg8CEiIqIGg8GHiIiIGgwGHyIiImowGHyIiIiowWDwISLJ2bt3L3R0dHD//n3FMn9/fzg4OKCgoEDEyoioruMjK4hIcgRBgKOjI3r06IGQkBCsXLkSO3bsQFxcHD766COxyyOiOqxBPp2diKRNJpNh9erVGDlyJMzNzbF582ZcvHhREXqGDRuG8+fPw8PDAwcPHhS5WiKqSzjiQ0SS5ezsjBs3buDMmTNwd3dXLI+JiUFhYSF27drF4ENESjjHh4gk6fTp00hPT0dZWRlMTEyUPuvduzf09fVFqoyI6jIGHyKSnKtXr2LUqFH47rvv4OnpieXLl4tdEhFJBOf4EJGkZGdnY/DgwVi0aBEmTJgAOzs7uLm5ITExES4uLmKXR0R1HEd8iEgy8vPzMXDgQHzyySdYsmQJAMDFxQVeXl5YunSpyNURkRRwxIeIJKNZs2ZIS0urtPzo0aMiVENEUsSruoio3vH09MTVq1dRVFSEZs2a4fDhw3BzcxO7LCKqAxh8iIiIqMHgHB8iIiJqMBh8iIiIqMFg8CEiIqIGg8GHiIiIGgwGHyIiImowGHyIiIiowWDwISIiogaDwYeIiIgaDAYfIiIiajAYfIiIiKjBYPAhIiKiBuP/AW2h2Pz9Y/EIAAAAAElFTkSuQmCC", "text/plain": [ "Figure(PyObject
)" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "P(C1|x•,θ) = 0.9999999132138306\n" ] } ], "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": null, "metadata": {}, "outputs": [], "source": [] }, { "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": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "anaconda-cloud": {}, "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Julia 1.6.3", "language": "julia", "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.3" } }, "nbformat": 4, "nbformat_minor": 4 }