Skip to content

Instantly share code, notes, and snippets.

@tomekkorbak
Created February 9, 2021 17:37
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save tomekkorbak/f498fac77e990fea5d93c39ea8dee84d to your computer and use it in GitHub Desktop.
Save tomekkorbak/f498fac77e990fea5d93c39ea8dee84d to your computer and use it in GitHub Desktop.
EM for Gaussian mixtures using einsum
Display the source blob
Display the rendered blob
Raw
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "EM for Gaussian mixtures using einsum",
"provenance": [],
"authorship_tag": "ABX9TyNhN/X90U6wayYxNTNDdMvB",
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/tomekkorbak/f498fac77e990fea5d93c39ea8dee84d/em-for-gaussian-mixtures-using-einsum.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "code",
"metadata": {
"id": "dDFPJymsTfa8"
},
"source": [
"import numpy as np\n",
"from scipy.stats import multivariate_normal\n",
"import matplotlib.pyplot as plt\n",
"import seaborn as sns\n",
"sns.set_style(\"white\")"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 271
},
"id": "ROXVDNLaVvF8",
"outputId": "9feb8195-70e1-478a-e082-bd031d149f7b"
},
"source": [
"def create_dataset():\n",
" x1 = np.random.normal(size=(100, 2)) + np.array([-1, -1])\n",
" x2 = np.random.normal(size=(100, 2)) + np.array([1, -1])\n",
" x3 = np.random.normal(size=(100, 2)) + np.array([0, 1])\n",
" return np.vstack((x1, x2, x3))\n",
"X = create_dataset()\n",
"plt.scatter(X[:, 0], X[:, 1], cmap='gist_rainbow', s=50)\n",
"sns.despine()"
],
"execution_count": null,
"outputs": [
{
"output_type": "display_data",
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD+CAYAAAAuyi5kAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2de3QUVbb/v1XV3XlcSKK8gyBBBgyIMN5E5eUj+SnImxnn5+goMBfkqjNrZs36eUccHQmgA/jz6nXduVcEFMjouLy/NQPhkRnAwCgSAmEEgSTDK0AgDyDhETCP7q6u3x+hmu5OPbtPd1V39ucfJd1ddepU1T7n7L3Pd3OSJEkgCIIgEgbe6gYQBEEQbCHDThAEkWCQYScIgkgwyLATBEEkGGTYCYIgEgwy7ARBEAkGc8P++9//HsOGDcPx48dZH5ogCIIwAFPDXlFRgUOHDqF///4sD0sQBEGYgJlhd7vdWLJkCQoKCkz9zuv14vz58/B6vayaQhAE0aVhZtjff/99TJ8+HXfccYep3zU0NCA/Px8NDQ2smkIQBNGlYWLYDx48iKNHj+KZZ55hcTiCIAgiApgY9vLycpw6dQr5+fnIy8tDQ0MD5s2bh6+//prF4QmCIAgTcNEQAcvLy8PKlSsxdOhQ3e+eP38e+fn5KCkpMe3GIQiCIDpDeewEQRAJhiMaB925c2c0DksQBEEYICqGnSAI+9LS5sHuQ3Wob7yBfj27YcLoTKQmO61uFsEQMuwE0YWoqG7C4jVl8EkS2t0iklwCPtp0FIvmP4gRg3tY3TyCEeRjJ4guQkubB4vXlKG13Yt2twgAaHeLaG33+v9OJAZk2Amii7D7UB18KklwPknC7kO1MW4RES3IsBNEF6G+8YZ/ph5Ku1tEfeONGLeIiBZk2Amii9CvZzckuQTFz5JcAvr17BbjFhHRggw7QXQRJozOBM9xip/xHIcJo0mVNVEgw04QXYTUZCcWzX8QKUkO/8w9ySUgJcnh/zuRGNCdJIguxIjBPbB+0UTsPlQbkMfen4x6gkF3kyC6GClJDjz+wJ1WN4OIIuSKIQiCSDDIsBMEQSQYZNgJgiASDDLsBEEQCQYZdoIgiASDDDtBEESCQYadIAgiwSDDThAEkWAw26D00ksv4fz58+B5Hqmpqfjtb3+L7OxsVocnCIIgDMLMsK9YsQLdu3cHAHzxxRf4zW9+gw0bNrA6PEEQBGEQZq4Y2agDwI0bN8CpqMgRBEEQ0YWpVsxrr72GPXv2QJIkrFmzhuWhCYIgCIMwDZ6+9dZb+Nvf/oZf/epXePvtt1kemiAIgjBIVLJiZs6ciX379uHKlSvRODxBEAShARPD/t1336G+vt7/7507dyI9PR0ZGRksDk8QBEGYgImPvbW1Fb/85S/R2toKnueRnp6OlStXUgCVIAjCApgY9p49e+J//ud/WByKIAiCiBDaeUoQBJFgkGEnCIJIMMiwEwRBJBhUzJogIqSlzYPdh+pQ33gD/Xp2w4TRmUhNdlrdLFMkwjUQtyDDThARUFHdhMVryuCTJLS7RSS5BHy06SgWzX8QIwb3sLp5hkiEayCCIVcMQYRJS5sHi9eUobXdi3a3CABod4tobff6/253EuEaiM6QYSeIMNl9qA4+SVL8zCdJ2H2oNsYtMk8iXAPRGXLFJADkH7WG+sYb/lluKO1uEfWNN2LcIvMkwjUQnSHDbgEsDTH5R62jX89uSHIJioYxySWgX89uUT0/i+fI6msgogMZ9hjD0hAH+kdl5Bd08ZoyrF80ESlJdIujxYTRmfho01HFz3iOw4TR/aN2blbPkZXXQEQP8rHHENaBKvKPWktqshOL5j+IlCQHklwCgI5ZbkqSw//3aMDyObLqGojoQncthhgxxI8/cKfh45F/1HpGDO6B9YsmYveh2gCXSP+oGkTWz5EV10BEF7pzMYS1IbaLf7SrB29TkhymDGmkRGNAj/U1ENGFDHsMYW2I7eAftVvwtisMMnYZ0An7QoY9hrA2xLJ/NNSw8hwXE/+o3YK3SoPMmqIjmDwuCxwQFUMf64Gkpc0Dt0eE1+tT/JwCngQAcJKk4qyLEefPn0d+fj5KSkpwxx13WNmUmKBkfGRDHO4Mt7Xda4l/dFvZWawuOqI6c1wwc2TMlvctbR7MXbJdM3DIoq8Dica9NHO+QKJ9biK+oBl7jIlGoMoq/6idgrdaAUUZpdVEuDPuWK9WlM4n4xA4zJ0yHPm5AyngSQBgZNivXLmCX//616ipqYHL5cKdd96JJUuW4Pbbb2dx+IQjUQJVdvL1ag0yociZI/17dQs7PsA6M0UPrfMJAg+XUyCjTvhhksfOcRzmz5+Pbdu2YfPmzRgwYADeeecdFocmbMyE0ZngVeraxtrXKw8yRmh3i6hpaDacC97S5sG2srNYt6UC28rOoqXNw2S1onRcNey0OiLsD5MhPiMjAw888ID/36NHj8Znn33G4tAJSyJkb8QyeKvXX1qB6VCSXAKut3gMzbjVsn6eGDsootWK2WwiO62OCPvDfO3m8/nw2WefIS8vj/WhEwa7pQhGQiw2txjpL6VBRg2e49A9xak7A9byoxfvOQ2VxYrmaqWlzYOS8nP4ePNReMVbA4uef94Oqa1E/MDcsC9duhSpqal49tlnWR86IbBbiiALohkzMNNfoYOMJAFbS09DktBpNXH+4g3dGbCWX1sCMHVsFopLzxhercgDlMcrBhn1QAJXC6GrlIWzc7C88IAlqa1EfMH0aVixYgXOnj2LlStXgudJhkaJWAfd4h2z/RU6yDz12DDF1URWZpruDPjzHcc0Z/UcB8OrFa2sltDj1jQ0Y+2WChR9eQocB3hFyW/EX52Ti0tXW2nrP6EJsyfi3XffxdGjR7Fq1Sq4XC5Wh004KAhmDr3+2lZ2BpIE1RiF2moi0HUj+nxwe3zgeUDgeSycnYOUJIchv7bR1YqRdEwAcDp4FJeegSdkA5LchmXry+NyVUfEFibT6hMnTuDDDz/ExYsX8eMf/xgzZszAz372MxaHTji0sjcSLQhmJutDDb1sl+M1V7G66AjmLtmOiuomU8ceMbgHXp2TC0kCBJ6DzwfwPIflhQdQUd3ENOvHaDqmx+vrZNQDIdVOwghMhv3vfe97OHbsGItDJTxdJQgWC71wmXBjFC1tHixbXx5kSEOPpTerN4rW7B/o2GTEcRw4AG4Nwx4Y2I33rCoiepAjPMZ0Bf1rPb3wy9daDc/klfpLDbOzWS33iOjzYfehWt1ZvVG0Zv8Cz2H+jHvwxNhBmkYd6HhWJAmYu2Q7VhcdwZ92nQx7xUIkLqQVYxFK+i6SJCXELExLQ8bl4CGhwzia0VeR+2tb2Rkcr7mq+r0n84ZgzpQRhtq5bksF/rTrpOrnj9zXHy/+cJSqBk1KkkN1haA0oz5d16ypLaPVbzLJNwe3NoXvaLWH6FrQE2ARoUG3RMpt1/Inh85IjbpR5P5qd4uorr2mmC4YGqOQc8bLqy6Ag4Sc7L7Izx3gHyz79ewGl5OH26M8S95zuB7D7rzddBaT1r3UyqLRczsluwRMHpeFrXtOm2oP0fUgV4wNYF0yz2rMbO+XMeJGqahuQmFxpWoOeGCMoqK6CbMLtmHVxiM4eOwivjl2Cas2HsFzBdv8LosJozOhtV7leQ7llRdMZTHp3UsAePyBOzFnygg8/sCdQQOZktvJIXBwCBx++OgQFBZMAgdQVhWhCxl2G5BotUu1/Mlq6Bkl2WAquSAAIDlJ8BvFljYPClbvRbun83fdHhEFq/eitd2L1GQnxo7sp9kmjpNMZTFFei/lTVYLZo7Ek3lD8OIPR+GPSydj7tQRQSmYRttDdE3IsNsAVrntLNILWaAWIHY6eDgdyo+cnlHSMpgOgcOcycP9Lqvdh+rgFdWDkF7R5zewI4f00jSUOdl9TaU8sriXsttJaVZvJ+E1wr6Qj90GsBB4spuPXklDJufu3nhhxU7FPG09o6RlML2ihKZrrUHfVXPXyN+XDaxe+ml+7kBkZaYbFjozcy/DSVm0umoWER/QU2ADIs1tt6v+jNKuTD2jpGbszBjMfj27wSFwqsbdIXD+7xsxlGaEzozey0gG4lgIrxHxDaU72oRIyqxZVaIu3E0yaqX8KqqbULB6L7yiD15Ruhk45FHw/BhkZaYZSjtUU08MJNkloLBgUpAhZFleUO9eapXxo5RFggX09NiESGZhLPVnjBrrSGacSjP5ljYPFq0KDnh6RQleUcSiVXvxh8WTOs2sHQIHSQKeGDsIkiQFtUnLFTN3yvCIyuLpoXcvSQiOiDZk2G1EuPK3LHz0LW0efP7FcWz66hSAW4qCSsY6Gq6fv+49o5jFAgDtHhE7D9RgyrjBWL9oIj7fcQxFN9sp+iRs3XMaxXtuyvOqHCOQ9cVV6NezG5YXlkctJqF1L0kIjog2lBWTAESaKVFR3YQ5i7fhz7tO3pwld8wm1XLpWadnVlQ3Yd3WSs3vlFdeAABIkoTi0jOd2tnmFg0ZdbmNb67dZ9m+AUpZJKINGfYEIBL9Gb38cKCzsdabcR4+2Wi47fL5jUZ6jMrfatHuFiH6lNMhY7FvgFIWiWhDrpgEQc2vK0kStpWdVfUjGzGUoe4BPaXC0sN1aH1ylCF3jFFDnZPdB4Bx+VstBJ6D6FM+Z7tbROnhWtRdip5eD6UsEtGGnqAEIhz9GSOGMtQ9MGF0Jj7ccFj1+xwHwwFAI+d3OQXk5w4EoD+oGIHjoKkR8+2JRvz9H5eiuheAUhaJaEKumATFqP6MEV2XUPeA3lZ8t8dnOABo5Pz/MnV4kFCWmhvD5RSQnCQYuh4tyQO9GIMaZnf+au0wjTZKbbXLzmUicmh6kKAYTanTVRQM0GAJZOSQXth7tF5x1msmADhhdCbWFB3R/M664krk5Q5ESpJD140xuH86Sspr8JfSMzh38bqi757jOUwZF1yIWuABNRUCIymIdtv5q4VSW1dvPAIJHasZu7ef0Idm7AmK0ZQ6NUVBgb+pKLhokuKLPWF0JgSVguVmAoCpyU48puOy8XjEoIBmqFDWgpkjsX7RRIwY3APVtddQWFylatSBjus/ef4annpsKEZk9cDQARmQJPUZvFGBMpZZNuHOnvV+p9bWdo8It0dMCHVRguGMfcWKFdi2bRtqa2uxefNmDB06lNWhiTAwk9sejr+XZQCw9uJ1zc9FHzoZVrVNTkYN0bcnLuHbE5cMtS8SgbJwNhyFO/s38juzWUXiTcG0eNgwReUCb8HMsOfn52P27Nn4yU9+wuqQRASY1Z8JZ3MUqwCgBH2J3+raZmwrO6v5srJIhVSi3S2iZ3qK6uesd/6Gs/nL6O/MZhW5vT6ca9AeeO1APLnCYgEzV0xOTg769VMPqBGxJVa1VVkEAHNvpjJq8c2xi7q1PVmkQqqxbP1+1ZWAVgDYIXD+QcmIO0VrcPJ4RZSU15j+XWBufjhFUJpb3Ka+H2sSrVANC8jHnsBo+aLNEO1sifzcAUhy6hsbvZc1HKNllDa3iM93HFP8TCtTxytKqoOSUr/qyRN/tOmo4sBmdNUQThGUtFR7uzMSrVANC8iwJziRzqgrqpswd8l2rC46gj/tOqk5aw53AEhNdmLxgjFIdgngDdgctZc1HKNlhqKvTikOKEqro1BCByW1fpUAzcHJK0qKA5tRmQK1YLkaTgePAX3TVD+3A6S90xky7IQqZpa4ZgYANUSfBJUNoUGovayBRkutUlMkSBJUZ3+Bq6N/vruXqrH0SRJKymtU+7V4z2ndiIPSwGZGpiB0JTdv+j1IVnUl8UwlDqKx+iPtnc5QHjuhitFsj0jVHuUapUqVlZTQelkDA7rnGprR3OJBWqoLgsBh657TkALaZhbRJ2kGEuXVUd2lG/j7P5QzbtrdIsorL6j2qwRg6rgsbPzylKbsQejAZjZLKTRYbqZKVLhEK8AZaaGaRISZYX/zzTexfft2NDY24qc//SkyMjKwdetWVocnLMDoEjfSdD+9GqWh6L2sahk+Tz02DCXlNVhTdFTVaOoRGEgMt9oToD64dBTQBubPuAcfbVIuFqI2sEWSpRRtiYNoVvki7Z3OMLvi119/Ha+//jqrwyUc8ZhjazQXPlIfp16N0sBzRvKydrhoBDgcPMQwZ+1yIFFr9qk3g8wZ3gcVp5s0+3X8qEwUFlfBK3b26WsNbOFq+kf6Wz2iXVyEtHeC6ZpXHWPiNcfW6BI30kIfejVKeQ74/rDeyMnuAwlAeWUDzl+8EdbgGElKpMBzqLlwA5t3V6OwuDJI6jh09qk1g8zKTMMfiqtUrpXzG6REmoUaGfwjnfxEc2CKN+Lr6YhD7FpoWkbrZTK6xI3UxynrxXhF5Rff5RIwfcJdTCoeRaIOKfo6UhcPn7ykOggFzj61ZpBGinqfv3gDj98/ENdbPeie6sTAvmlxOwvVG/wlCZi7ZHvcTX7sChWzjjJWFZo2gtEC2kYKPUdSjFv+/Rur9sIdUgUpySngN3NzsbzwAJPiz1qFpAPPyXEIO9A6dGAGHn9gkO6MU6uodyR9GYod3IBa/S5n5CgVe6Hi3uFBhj3KrNtSgT/tOqn6+ZN5QzBnyogYtqgDrRct3JfJyACg9/uS8hocqOoog5c7vA/ycgbiq4O1TAdH2XB6RV9QJo7LwUMQeL9K5O5DtSg9XItvTzQaigGEtiscY8z6vrAeJCJBrS1PjB2ErXtO23LyE6/QMBhlWBSajgbRCGZF6uNMSXJg6vjBmDp+cNDfWW9AUUuJHNC3e9BgpJe6qIWSu83IzJnlfbGbG1AtwPn5jmO0wYgxZNijjF1zbPWM5Z5v66JaHs4M0RgcjQ5Cej55raAvcMsY9+/VzVAAncUgJg8gew7XweNVPhaLTJTQ8xlx9Sj1u10nP/EM7TyNMrES4zKLnq7K4ZOXwt5BKsNql6GVxZ+1zp2cJGDe9HswdGCG6u/b3SJqGpoN7+DtoaEiaeTzwB3AB49dVB10WM2E9XYcG3kGqLg3e2jGHgPsmGOrVzkpsDwcYH7pHmmKZ+gscOHsHCwvPBDz1D+9zKARg3vA6RBwtkE9BnC9xWPYvdLQ9J1me7SkcMzo0bOYCeu5ehbOzjWUyUQbjNhDPRYj7JZjK79MSpkoSphZumu98AWr92LOlOFovNqqumyvqG5Cweq98Io+eEUJDoGDQ+Dxm7n349LV1ogGx3AyRPQGZj13W/cUpyH3SkubB8WlpzXb0ni1VfUzM3r0ajNhM/2jFw94c+2+oOC01iTBjpOfeIZ6rQuTlZkGlep2nTCzdNd64dvcon+rvNIMrqXNg0Wr9qI9YLDxihK8oog31+7DJ4ufCPtlD11FuJw8PtxwGGNH9sPIIb06GTGjRi412YmFs3Pw5tr98PkkiD4paMZ5/uINQz7k3YfqNNvvEDjNWbaRzVdaM2Gzqyy9eIDas6U2SbDb5CeeIcPehdl9qE61LmgoTgdveOmuZ2C03Dwl5eeCjHogbo8Pn/ylCs/PHNnpMzUjLP+9pqEZf9l7JmgGKRfi/vJgLcoqGoKMmBkjV1HdhOWFB8BxHZuYeB7w+ST85qf3Y8TgHsjKTDMUQNeTVpAkaPqbtYKQHICBfbtj0phByL9ZGDy0/8xm0GidT+A5UyJmBFsoeNqFMbO93uP1Gap0BJgveBEoQ1t+M4ddjc1fV+Nyc1tQQO6bf1xQDOBt/bra//dNu6s11SMDg5mXr7WqBjt/+2Ep1mw84g8EysqUre1e/0Dh83X01/LCcrS2ew0H0PX6bebDd2muVrSCkBKAhsstKCyuQnXttU6fh1OsQjPoyXNwOZXNC2W6RB+asXdhzGyvdzl4lFddMLRU1gvMhhI4g+OgvYSQJGD+WzvA85zfnSIb1MDjAcDKDUcMt0HGJ0lYX1ylUZ7Oh6Ld1XAIHNYUHcH9I/oq7piUjyW7HIz4kLX6LTlJwFOPDdNse2AQUvT5VPtFaQYeTpqlVtBTDnYDnQdTynSJPjRj78KYqTjk9voML5/NVulJcgm43NyOdVsq8E8pLt3je7w+vxEKNV6R0u4Wcf7idd3BzitKaHOL+Oqgetm1UIOoV81Ka2ZfMH+MKfnd8aP6Q1ApR6U0Aw+3WIVa+cX77u5jyzTfrgL1rs2Ipa6H0oxLDbPL59AZao/0FKwvrkRbe+dztLtFfP1tLdwen+ryPVYkuQTc0bs7zjboG3c99IKdSrDIDklJcuC27kmmfNyRbKRTC3pSpot1UA/bCLNZCSwGgcCXr6ahGcWlZxR90eEsn41U6QmdebOegZuF5zjMmZyNvUfqIz5WaLDT6P3Syw4xcpwe6Smqu2KVBulo5ZJTpos1kAiYTTAr/hQtcadoi0YFCoVdbm73z9RDcQgcJAmQIMHnA1xOHpLUkd3hNlhCTw2eA7L6Z6Cmodnvq09yCeAATB6XBQ4dwcbiCEvp/fDRIZg7tUPgjVW/GjmOvA9AzfevJSYWqZAbYQ+YGfbTp09j4cKFuHr1KjIyMrBixQoMGjRI93dk2DswI+8bDWXGwFlgj4wUcACarrVG9eXWU76c+dBgDOib5jcyOdl98MLyEtWdlfIKQDZ2cyZn+wOhof0qG/Ip47I6JHolBNVETXJ1SPdOGZuFk+ev4dsT5oTAkpMEFC6a5Bf/YnG/jBxHkiRNWeJkl4CC58fErfwvYQxmb+uiRYvwzDPPYMaMGSgqKsIbb7yBwsJCVodPCLReDDNZCayVGVnMJs2+9C1tHlxubr+Z89358ySXgAF90zpdh3oWRi4ar3XelZqXO1Cxzqnc18WlZ7ByYT5eWF6iWBGpuPQMVr6ShwXLShTz650CD8HRsbpQc1+wul9GjiNJUP2OQ+Awd+rwsIx6vFYB66owMexNTU2orKzE2rVrAQBTp07F0qVLcfnyZdx+++0sThH36L0YZhTuWMrYspB2NfvSB35fyagD6j59swE5vTqnPknC+q2VmgbzwD8uYvGCMShYsxdeb4DMgYPHq7NzUdf4HcorL4DjJORk9+20AYjV/TJyHHmAUcIrSpqSBGrYTf6X0IfJ3aivr0efPn0gCB1pTYIgoHfv3qivryfDDmMvhpmsBJYyp5HOJs2+9HpCVS4nD4HnNQN2ZgNyegZRK71RNpiPP3AnChdNChpQemWkYNn6YJGrytNXkJWZHjSgad0vl5P3p3oq7ZgNXAEZue+SBGbPhtyG0sO1MZH/JdhBw2wMMGo8jWYlsNR4j3Q2aXZg0Pq+wHMYP6o/XvjBvUxngHoGUSu9MdAYBg4oSv5utQEtJ7s3PtygfM1ujw97DtcFrXRCYwPy359+fKjqvZLvuyRJTJ6N0FWYGiQPYE+YJA3369cPFy5cgHizGLEoirh48SL69evH4vBxj1HjqbbZI9SdEanGe6BG9uXmds2t3/JsUk1L2+zAoPV90Sfh9rQk5st6Pb3vOVOGm9YDN7oFv6K6CS+u2InQowdu2AqVLVi54YiinMHHmytVr/HVOblISXIw0f8PXFUZERUjeQD7weQN6tGjB7Kzs7FlyxbMmDEDW7ZsQXZ2NrlhbmLGdWLUzRDu5g8lf7ha7njgxiE1v7lZt5CR77POvtDL0b49Ldl0DreRAU3WkVFKO/T5JLgcfMSpm0CHO+dSgO880o1BLOR/CWthNjUqKCjAwoUL8d///d9IS0vDihUrWB067olWeTyzvmYtfzgAv+5KoP6K/F81N4PZa9P7fs/0FMxdsp159oWesTNrDI0MUJ9/cVxDRybyfHwZt6ez3EMkG4OMyv9yAJ4YOwif7zhG6Y82gzYoxQg7VIvXy5Ufd28mfD4fdh+qU92OrlQ13uy1qX3/Vz/+Pv7vp39X3Pkabn5+tNDLKV/5Sh7+5c0dqv3IEqV7EsmqR+s5cQgcRn2vFwb1S+uU9x/r55lQxx5vSRfADroZeu6D7qlObN9Xo2mMQv3mLW0enL94A4/fPxDXWz3onurEwL5pmtem1Be9MlKw9ON9qtK6dsu+0HPvlFdd1CxjxxKfT0JOgKRypDnnE0ZnYtVGZWVMnufxi/89Gi+s2KmY90/pj/aAej+GWK2boec+uHK9XTWtLfB7st9ca6au92IrZZjo6aXbKftCb0Arr2zQLJoRiFyUQu6/x+8fiKLd1Z2+N2PCYGzfXwOv6AvqKw7AC8tLsGj+g8jKTGOSc642JnEA9hypZ7pBjmAPyfZ2ITRleiVgz7e1usZI9psrZU4EFqswUlBZxkiwjmX2RWBWkFq2jxYV1U3+Ah5Fu6ux53Adtu+rQf9e3QwXzZBxOXk8fN8d/iyola/kYfv+GsXvbt9fg//41cOd/u72+vz9XnLgnOmCGaHsPlSnadkPVF5gtkGOiA40Y+9CqLkPOHQE80SdWJ7TcWvj0Lays8xmbUaCdayyLyJ1UxjdkGW02IjA80F5+9vKzsKrciO8og+ff3EcvIbOejhGN0gnKD0FZRX1mscA2G2CIqIDGfYuhpJ/u90jYv1W9RxpoMOof/TaY7gtLRkAW1kDvUpOgQNKJLDYGm90Q5ae1r3aDtuahmZVl5TH60NNQzNTo2t0I1LgMXKH90HVmcuKn1P6oz0gV0wXJLSST9PVVs2X2iFwWPqvY/1GHQi/4o4SWi4ieUCRZ9NG3Sih32u61oqVfz4Ct0qhbKNuCjMDmjyITh2X1aEtc3NTksB3iIa9Oie30yrheou2W8gh8Jr9nju8j+HNVmY2IgUeIz93IFVHsjl0BwjNGbND4DBv+j2dDBDL3Hy9DBN5QDHqRgn9ntPB4/f/zwee63A5KRFolLVSBbX6iueBy83taGnz+L8vSRKKS88ExS5EnwTRJ2HZ+vJOq4SUJG2//J390nC24briZ7LRVSpoohTUNrMRySFwcDoE/zHskOVFqEN3oQuhZrC0jLTTISA/d2Cnv7OuuKNnKIy6UZS+J7s2tFLK5VWGXKTCK95ScVxTdMSvYa7VVz4fsOdwHfYeqfcPNma1dFoVSgcG4vH6OvW7XITk/uF98NXBWkwYnWnI6BqJbciM+l4vvDI7N2QQsmd1JNKNJ8Oe0AQ+4KEVgUJnu+EYadazNi1DYdRAmpmFBsJzHHLu7o3nl5UEuWu8ovNyt4oAACAASURBVASvKOKNVXvxyeJJQQOa6PN1kmMIHWz0XDd7vq3D+FG3DE/3VG0DlJbqDOr3IycvYc/hevA8hy8P1qKsosF/X/WMrl5sQybJJWDsvZlxMRsn3fgOyMeeoASm5P1p10n8eddJtLlF1dREowJkoYT666P18hv1bZuZhQId/m7ZN7zncL2qD97tEVFS3pGGKPfV+FH9IWhkqOw+VKub9nj45CXMXbIdFdVNAICBfdPgdCi/lk4HjwF90wB09Pv4UZnYV3EBHq8vrJRTzfTXAOIlIMoyBTfeIcOegJgJigUGDWNlpMPBaLC2R3qKqrENheeBh++7wz+AHahq0Pz+vqO3ClynJDlwW/ck1V268mCjZzy9ohRkeCaMzoRDUH4tHQIfZGCNKkyqoaQEGUi8BUQj7Y9Ewv53izCNGXeEmdREK32XRoK1FdVNKCyuNKzPkuR0BOWQS6q7cjr49mQjvvnHBVy62uYvxq2XWhjouvF4RdUNYGZ1+VvaPH4ddyWMBoND3Wk9M1IgSdGvdxsNWKbgxjvxcccIU5hxRxhNTQz1XToEDh/86VvMePguPPW/hkbdwOsFayVJwuI1ZapqioGoKRMO6N0NB49dVP2dJAGLVpf5jXmgCmYoge4L2XguX1+Ob1SOH2h4sjLTMHtytmq5PfleaMk/BAaD9XzOkiRBkjquTy5VKP/biEagXYKVLCuLxTtk2BMQo0ExwJj/VCnTRJ55/nnXSRSXnkbB/DFRD05pBWu1dsIG0uf2FDwwoh/+svcMir48BdEnQeA5fLjhsM58/RZyvwYadbm/1YLOKUkOjL03ExWnmzQNj5IhDiy3p1daUEYOBr+wYqdmJlF17TXVDUpGAo92ClZGSx47HiEfewJiJChmxn9aUn5Oc3bY1i7GLDilFgcwukrpnurCpt3V8Hh9fpeN6JPg8frC1kd3OXmMuzdTN+isV8kpJ7uPbvBPz83mEG4Fg8urLmr6nEvKazRjMXqBR7sFK1lUj0oUus6VdiG0NGGmjMsCx8FUxaWPNx/VFQezWtXP6CqluvYa83O7PT7cnpaEOVNGALi16zXUNaEr9Vt5QTf4pzeABeabl1c2aPqcD1Q1GFrlqN3bSAuhRwPaONVB17raLoTWAy77RPUq38gzMiPys+1uETUNzYoGLRYYEd1yCBx8USh8oSdlHOia0LoveoZY/r7aAOZy8kjvluS/rz3SUzR9zpLEGVrlqAUe7RqstOvGqVhChj2BUXrAzfhEzWTXuBw8ikvPgOc5S3ytgbPhUL1ygefgdPJ4YHhffHkwvJQ3gVdXv1SSMpYJ9GmvfCUP5VUX/Qb6R/nBQWcjwb/xo9QHMLfH58+UkVdoULl9PMchd3gfVZ+/0rlDoWClfYnYx15UVIRp06Zh+PDh+OSTT1i0iYgSZn2iZrJr3F5f2BtlWCHPhl/4wb2Y+dBg5OUMwMyH7sJLT45C4aJJGDmkF1xO/UdeToNPcglwOniMH5UJTiNm8eqcDteH1kAoij7M/90X/g1jq4uOBG1MAvR98BNG91f0IwdeU2D/t7lFSACSkwRFn3NezoCINigZaS9hDRHP2LOzs/Hee+9h1apVLNrDHLukYrEmnOsy6xM14rdOcgl+94aS3Gysfa1ay/Bb7hrtIOk9g3vCcXP3Z87wPmho/E7VHeVy8qi9dAMXr7Rie9kZ1b4KDcwq6dwY1d8Jdedcbm5XzWnnOGDO5OFwOQVFn7OWtLCerARrvSCCHRH3/NChQwF01EK0G3ZKxWJJuNdl1ieqLQ7GY8rYQRjQNw019c2KpdzUjhtLQgfAhbNz8Oba/aqa5w6BQ9XZy36X0tHqRtVcdaDD/fHRpqMQBN6UlIFM6MBnNPgXOICt21KhGyT99XO5igN/pBuUKFhpTxK291kUVYj0/NFYKURyXWZ9onozMnkQ2VZ21nJfq1J/n65rVmz7r5/Lwdt/OKBo3Dtm5rdm51pGPfA3XtG8UQeUBz6zwT+9ldW3Jxoxd8l21YE/0mAjBSvth65lmzVrFurq6hQ/Ky0thSDo13W0AitTsaK5Uojkuoxu4Ag1kitfycOBf1xUnZGFszGE5cCn1N9rio7AJyFI1Es2fO/+8Ru8/tMHsLywHKIYfv66EfRcVSwGPr2MoI6BxxuTCQ1hD3Tv8IYNG2LRDuZYlYoV7ZVCJNdlxCeqPChxmjKwZn2tLAc+rf5WwydJaLzWipWv5GH+774wdT4zDBuYgccfHOTfAapk2FkEGcPRo5FJ1BhUPBDNvk/YoduqVKxorxQivS69/PZwByUjefNyseTC4sogTZdIBr5w9NflAbC8qiO4GA2cDh6PPzjIf6+jHWQ0q0cDJG4MKh6Idt9H/ERt2bIFb7/9Npqbm1FSUoJVq1bh448/xpAhQyJuXCRYpRsR7ZUCi+tS84lGOigZyZt3CJypGaUa8mChlYmihjwA1l26YciHHg4erw+52X38/zYSZIx0BqenRwN0CHvJ57IyBtWViUXfR3znpk6diqlTp0Z6GOZYlYoV7ZWC2esyYyxYD0pa4mHhnqOlzYPPvziOTV+d0j2e1nl6pndkf/D8LUVDJQSeg8sp4LH7B6K49LTh87kcPMqrLgQNUtLN+3WqthnVtdfQ7haRnzsAqclOZjO4CaMzsaboiOrnRV+dgtfrQ3OLx3ZyAF3FLRSL+F9CD8lWpGLFYqVg9LqMGgv5haqua1adUYczKJl1k+idQ65HakSaV4/lheVY+UoeBJ6HT8Wyy4U4XvjBvZAkCTv218ArGtts5fb6Ork9Fq3ai/aAYO43xy5h3dZKvDY3F8sLDzCZwaUmOzF5XBb+vOuk4udeUULR7mrNAa3dLeLwycaYGvau5BaKRfzPfsnnDGlp8+Crg7Wou3QDfXt0bMeO9vJSaWegwHNwOngsnJ3D7Px61Y70dplebm7DtrKzeOeTA3h20V+xuugIDh67qDoj5QDTg5LZMnVaA598PSyMOtAxMzrwj4t4/af3q34nsBCH0n11COoOeoHncLm5HS1tHrS0eVCwOtioy7g9IpZ+vA+iipXVqvwji42t21KBbWVn0dLmAQBD8sNaqxQAKD1cF7Mdw3ZTiQxtm1IfR4LRamCRkLAzditnACMG98DC2bl4c+0+8HyHLKwgcFheeCBmMxCt2bJX9GH+WzvAccbytIEO32x17TVTbdfLr5ZXB0ZcZOEWqVZDnhn179VNsWCGyyl0ak/oSqlHegrWF1eirb3z9Yk+CV9/W4u9R+rxxNhB8KoJzdz8rtqAqjaD03q+zejxq8FxiJk7xo4qkUD0bEgsVvUJOWO3egbQ0ubB8sJyeLw+/8zI7fHFdAaiNVv23NR1MRM4bPeY11zX0hIReGDcvZmY8dBgQ4Wzjcz+HQKnOYsOJMkloEd6ChavKVPsB54HBvdP7/T3wJXS1PGDUTB/jGrNUPmeF315StM3L/v61doZOoPTe75zh/cxpAGjhdvji1pKcCh2VImMpg2JhW58Qhp2q4vaWn1+QHu5Fy5m265VLFn0AWUVDdi+rwb9e3XrlB0Suvw1cj0OgVctBB2KbPjU7pMkwdC1yrP4cfdmqhpnjgN4jQLbDoGDoPJjpRmc3vN1oOqCZpFqI8RSnTEWrgmzRPsdlp+bBTNH6hZoCYeEdMVYPQOw+vyAMX1ys4TTdvkBLimvwUebggt2KAUI1Za/C2fnaM5Ck5MEFMwfAwC6AVaHwOH+4X2wv/KC7n0ykqmRkuTAbd2TVP3WXlHq0IJXbQ+PV+fcj+WF5YYynfSer5qGZkgS8Pj9A3G91YOUJAHb99Wo6uMoETqgRDNjxY4l7WLxDkdTiiFuDbvWg2a1TrTV55f75v7hffD1t3X+EnAs6JmRYrodsgKhIPCKmiryDGj8qEzV/N7lhQfw6pxcLFt/y/g5BA6SBMx8+C489dgwvwGcM2V4p0EkEAnAlwdrNd02HYUogLlLthvyserd86njsrDl69OdAqgup4CC58eYyuDSOpdTQRef5zjMmzYC64urOg0ccyZnK/49cECJdrwqVqnJZgYnrTcmHrTm49Kw6z1oVs8ArDx/aN9oeAAU0dpABNza4GK6HTrpdfWNN3SXv5euthoyfo1XWzWvQbz5mdZ3OABbS08HBUa1UhD17vlTjw3DU48Nw84DNSivvAAAyMnug/zcgf7jGJ3BaZ0rdFYut3l9cRU+XJiP8qoLnfouL3egap/GaiNTtFOTzQxOLW0eFO85rXqscDLEYk3cGXYjD5rVOtFWnV+pb4xM1l1OHpIEjL03E9e/c6tuSQc65FzDaoeGF8Dl4P07QbWWv9vKzkCSoOsGCCcrJDBDhwOQnXU7Dp9oVPyuUqaG0Xs+ZdxgTBk32HC7lFA7l88ngUNn7Xe5zaEbpmS0BhS9wbakvAZOh8DERRMt14TZwWn3oTrNGfuUcVm235Vr79YpYDQ1ymqdaCvObzYlsENTPQsD+nb3t21b2VnVLelGl6Bm2+G+uf1+vwRNg3y85irONhzRdQOEE18Y9b1eyMpMgyQBW/ecxrcnLumuMEKJ5T1XOtfZhmZs+oqtLr6erzlQi96um4rMplPqZWBFS1+IJXFn2M0ENazWiY71+fUeSIHnIPqC88ZDX0AWbiSzG5OcN7ffGzHI8nFf+2AP5s+4B3k5AzrNEJVmtFouJrkItNvjw1/2ntENMmoNcLG856HnioYuvt7qJ1CL3q5aM2YDoVbHyFgQd+mOdkyNsgtafeNy8nj4vjt0U6uM5thq7cgzm2rpubn9Xis9MhTRJ+GjTUc71Q2VCU0nmzf9HiSrHFMuAr1pd7WhzBG71vOMRg1SrWOqEauUXqOYtRmJUMs17gx7InR6tNDeEMTjhR/cqypBEIhejm1FdRPmLtmuWpjZrDEIfLnkc8+dMlz3d15R0tww0mkz0fPBm4mUikBrwXNAsqvzjlS7EI2NL2rHFDSi8lZtKlLDrM2IxQaiaMNJEsN92mFw/vx55Ofno6SkBHfccYeh3yhFuNVcC12NaPdNS5sHc5dsVzWkWvnoasYzJcmBlQvzUV55K2PD7RHx8Wb1lMVAklwCFswcacgF0truNVQEWgmOA5IC0hPtSuA1svLzhx7T7RGxbmulqrvC6P2IFeG8F9Hox1gRl4YdiO9OjzbR7JttZWexuuiIoRc6tB0901MUN+Eo5VKLos+UJO+TeUMwZ8oIU9eybksF/qSigqhF4ADWVTE6wNuJrmQz4vaqrA6M2plo9k2kwevQTI6c7D54YXmJqdJ2oejFVpqutWL91irUXryO/r27Y86UbPRITwlbLEv0+SwTprILVqcUh0NXshn2633C1kSaMaCUyRGpaqNWbGXr19VYueFW4Ynj565i19/P4YVZI/FozoCwZBdiKZBlZ6xOKSbUifgOLF68GHv37oXL5UJqaipee+01jBw5kkXbCMaw0PtgvavWbGpkIHozxKZrrUFGPZCVG45g7L2ZirNOPTcQzyPi7KtYVQuK9nm60iw4nojYx75r1y6MHz8eTqcTu3btwltvvYUvvjBe+T1cHzthDpZBVTPH0jMsWj57LXgeeOS+Af5CGEq8+8dvsOvv51SPkZczAL96+r5Ovtd2j4j1KoFBoCPv/tMlT4Q9M41V8F/rPFmZaV2iDF1XhWnw9MqVK5gwYQIOHz4MXk3DNAQy7NGnpc2DOYu3KSoehhvoMhKIMmLAtIJweugFTP/Pf3yJ4+euqn4+bGAG3vnlw53+rtemxc+PwX139zbdXr1jsww6ap3H5RTA8x26P5RVlpgwzWP/9NNP8cgjjxg26kRs+PyL46oytuFuJtErzdd0rRW//bBUt1CBUs6wkcfHiD+/f+/ump/37fFPin/XKm+4ZEH4Rh2InVa/1nncHhFt7aLtytAR7NCdGsyaNQt1dXWKn5WWlkIQOh78rVu3YvPmzfj000/ZtpCIiJY2D4q+PKX6uZnNJEb9tRXVTfjth6WquzhD9TlCg3BG8suN+PPnTMnWdMWUVdSjorpJcZYarcBguDrfZn3l4cQurCxDx5JYxS/sjO5TumHDBt2D7NixA++99x7WrVuHnj17MmkYwYbdh+p0RYuMaKwblT2VlfS0tuYrGbDAIFxLmwd7j9Sr/j5w96fWS9wjPQUvzBqpGkBtd/s0dU2iERgMJ6soHD30cFI57bZjNBysrHVsJyL2mezatQvLli3DRx99RD5yG1LfeEN3o49elMVM/Uejyo5aX1Fyhcj1TH/46BAUFkzCiME9dKUNAGDK+MH4l2nDVXXpY61rYnZ7e7i1N8PReIl3rSWrax3biYgN+6uvvgqPx4Nf/OIXmDFjBmbMmIErV66waBvBgH49u2nqegD6Gutm/MJGXQBb95zWfNFC9Wpe/OEo/HHpZMydOsI/U9d7iWWhsq8P1anq0sd6lmpWhyRcn7zaeZJdApKcyoJYagOLmthbNAnnvHaoNWwXIg6/l5WVsWgHAfO+QSPfnzA6Ex9uOKxaHs/l5HVnaWb8wkZdABKg68+NpADE5zuOobj0jH9JroZD4NAj3Xi5PxaY8d9HUntT7TzVtdcM7RgN160RqY873PPaodawXUi4LWLxGjgx+zAb/X5qshOv//R+LFqtPAALPK8bhNQy1i4nj8vN7Vi3peKmREBvfLRJ3wUQ6Yum9xJv/PKUoVqvXlHC+q2VyMpMj6kP1qj/nvVOX8DYwCJnNQXGSozorUfq446kFF8i6KizIqHyEo34XO2IWd+g2e/fd3cfLH5+DJwO3p9K6HLyhmVItfy1spa53N8vrtiJOZOzkZLk0C0WHcmLpqWx7RA4U1Vu2tyibX2w0ZKp1kpXrahuwvO/+0I3qykUFj7uSNwpJOl9i4Qx7PEcODH7MIfz8N93d298uuQJ/OzJ0Xgybwj+dda9qsU2QlHy1yppmcv9LRdOnjf9HlX/fqQvml5w0IwyJGBfH2ystcHDzWoC2Pi4I3GnJIKOOisS5krN1jW0E2Yf5nAf/kjS98zkmsuFk6eOH4yszPSoKADKL3HB6r3w3tR26cic4TF5XBa27jmdMKl+LHPq9VyVRrKa1FZbLHzckbpTSJisg4S52ngOnJh9mK3yJQYODOu2VBjqbytetFFDeuIvpWdM/cbuPlgWOfVq/u+Fs3Nw6Wob6htvoLquWXdAVFttsXguWYjMkTBZArli4rkWqlnfoB18iWb6W09+IBxkl0GbW/S7XbyihDa3iOWFB/DqnFzFVD+XiVS/RELLVblodZk/LnXk5CXN4zgdvOpqi8VzSe4UNiRML7GWk40lZosW2KHIgdX9red6u3S1NaJUv3Cwc0aWnotFNvZasQmng8dHrz2G29KSFT9n9VySOyVy4rY0nhLxXgvVbOkuq0t9se5vM4ZRr6ydlvJjNPrNqJKlVYY/nDKADoGDV5RM31ern0siwQw7QA9VrGHV32YHCTO1V6ONESneqtNNeHPtfvh8EkTfLWMZ6N+OprEPR/f+n+/ujazMNHqP4pCEu1MUOIktLPo7nE0pVruCAtFzC33ylyps2l0d9Hf5+hatLoPLycPt8TEXrApcIfRIT4EZ5Zgkl4Cx92bSuxSnJEzwlIhfwsl/tlOQTS8ja/PX1Yqfybg9Pv93We27CN2st764EpLUUWRDaS9CKHaPSxHaJNyMPdGwc0BOJtI2hpuqapcgm1aan0PgDEkbBBLpvgutFVBykoA5k4ej6Vor+vXshl4ZKVi2vtyyIDwRHejO2Zh40JZm0cZI8p/t4HrTcgsB+rLIoUS670JrBSTP2gMDy3YYHAm2kCvGpsSDRAKrNtohLz8StNxC0x+6S9PloUSk+y7MroCisc+AsBa6gzYlHiQSWLVRK/954ewcfHWw1tauKEDdLSRJ0s1dsOraK6FEOpiRyiFBht2mxFoiIRw/Ocs2KhlGJf+v3VxRgai5heRBS/T54Pb4wPMdcsmzn7gbf9x+nLl/204ZQ4Q1kGG3KbGcdYXrJ2fdxtC6p6G54fJ5fvthKZ4YOwgD+6TZdgYfiFaQd+KYLOb+bTvsTCasJeE2KCUKRja9sNJcCfc80WyjkQ01sppjwfNjbDmDtxrarNd1iTh4+sEHH2DatGmYOXMmZsyYgeLiYhbt6vLEKk87Eg3taLbRSO1UWfTrjVV7bRFMthsUFO26RHynn332Wbz44osAgAsXLuCJJ57AuHHjkJ6eHnHjujqxyNOO1E8erTYarZ0KAG6PiJLyGkwdPziicxJEohCxhejevbv//1taWsBxHHw+4xkAhDbRztNm4SePRhv1csNDOXCzsAdBEIzy2D/77DNMmjQJs2bNwtKlS3HbbbexOCwRA+yaQ67k5iEIwhi6wdNZs2ahrq5O8bPS0lIIwq2X7tixY3j55ZdRWFho2LhT8NR67Cx3LAcAdx88j0MnGlW/98IPRmLKOJqxEwRgwBWzYcMGwwcbNmwYevfujf3792PixIkRNYyIHXbRXFFCdvOMH5WJ2QXb0O5RcBk5BeTlDLSgdQRhTyJ+c0+ePIkhQ4YAAM6dO4eqqir/v4n4wQ6aK1qkJjuxeMEYFKzZC683oHi1g0fB/DG2GIQIwi5E/Db853/+J06ePAmHwwFBEPD666/jrrvuYtE2gghixOAeKFw0yZYrC4KwExG/Ee+//z6LdhCEIey+siAIO0BTHSJuiQeteoKwAjLsRFwSD1r1BGEVpMdOxB3xoFVPEFZChp2IOyLRtyGIrgAZdiLuiLVWPUHEG2TYibhD1rdRgioEEQQZdiIOsau+DUHYBcqKIeKOrlwhKJ5SPOOprYkGVVAi4pauViHIzmJtocRTWxMRMuwEEQfEqlQiC+KprYkK+dgJIkq0tHmwrews1m2pwLays2hp84R9rHhK8YyntiYqNGwSRACs/MKsd8bGU4pnPLU1USHDThA3YWWMA3fGysiGbvGasrBcESxKGMaKeGprokKuGIIAW5mCaLgi4inFM57amqiQYScIsDXG0XBFKNWATXIJSEly2C7FM57amqhQDxME2BrjaLki7FzCMJR4amsiQr1MEGBrjCeMzsRHm44qfhapKyKeCo3EU1sTDWaumH379iE7OxuffPIJq0MSRMxg6RcmVwRhNUyesBs3buCdd97BQw89xOJwBBFzWMsUkCuCsBImT9ny5csxb948/O1vf2NxOIKwBNbGmFwRhFVEbNi//PJLXL9+HZMmTSLDTsQ9ZIyJREDXsM+aNQt1dXWKn/31r3/Fv//7v2Pt2rXMG0YQBEGEh65h37Bhg+pnBw4cwKVLl/CjH/0IAHDlyhXs2rULV69exc9//nN2rSQIgiAME5ErJicnB3v37vX/e+HChbjnnnvw7LPPRtwwgiAIIjwsD9GLYkfecENDg8UtIQiCiD/69u0LhyPYlFuux37gwAH85Cc/sbIJBEEQcYtSLQvLDXtbWxuOHj2KXr16QRCUCxQTBEEQythyxk4QBEGwhdQdCYIgEgwy7ARBEAkGGXaCIIgEgww7QRBEgkGGnSAIIsEgw04QBJFgkGEnCIJIMMiwEwRBJBi2N+x2LLn3wQcfYNq0aZg5cyZmzJiB4uJiq5sUxOLFizFp0iRMnz4dP/7xj3HkyBGrmxREUVERpk2bhuHDh9vmvp4+fRpPPfUUJk6ciKeeegpnzpyxukl+VqxYgby8PAwbNgzHjx+3ujlBXLlyBc8//zwmTpyIadOm4ec//zkuX75sdbOCeOmllzB9+nTMnDkTzzzzDKqqqqxukiK///3v2d1jycZcv35devLJJ6UFCxZIf/jDH6xujp/m5mb//zc0NEjf//73patXr1rYomB27twpud1u///n5+db3KJgjh07Jp04cUL6t3/7N9vc1+eee07auHGjJEmStHHjRum5556zuEW3KC8vl+rq6qRHH31UOnbsmNXNCeLKlStSWVmZ/9/Lly+XXn31VQtb1JnA93XHjh3SzJkzLWyNMkePHpXmzZvH7B7besYul9y77bbbrG5KEN27d/f/f0tLCziOg8/ns7BFwTz66KNwOp0AgNGjR6OhocFW7Rs6dCiGDBkCnrfH49fU1ITKykpMnToVADB16lRUVlbaZuaZk5ODfv36Wd0MRTIyMvDAAw/4/z169GjVwjxWEfi+3rhxA5xK0XKrcLvdWLJkCQoKCpgd03LZXjXsXnLvs88+w/r169HQ0IDf/e53tht8ZD799FM88sgjtjGidqS+vh59+vTxi9AJgoDevXujvr4et99+u8Wtix98Ph8+++wz5OXlWd2UTrz22mvYs2cPJEnCmjVrrG5OEO+//z6mT5/eSaExEiwz7HYuuafVttLSUgiCgKeffhpPP/00jh07hpdffhljxoyJmXE30j4A2Lp1KzZv3oxPP/00Ju2SMdo+IrFYunQpUlNTbVlo56233gIAbNy4EW+//TZWr15tcYs6OHjwII4ePYqXX36Z6XEtM+x2Lrmn1bZQhg0bht69e2P//v2YOHFiFFt1CyPt27FjB9577z2sW7cOPXv2jEGrbmGm/+xAv379cOHCBYiiCEEQIIoiLl68aFv3hx1ZsWIFzp49i5UrV9p6dThz5ky88cYbuHLlii1W2eXl5Th16hTy8/MBdBQcmjdvHpYtW4bx48eHfVxbumLsXnLv5MmTGDJkCADg3LlzqKqq8v/bDuzatQvLli3D2rVrmS7vEpUePXogOzsbW7ZswYwZM7BlyxZkZ2eTG8Yg7777Lo4ePYpVq1bB5XJZ3ZwgvvvuOzQ3N/sH6Z07dyI9PR0ZGRkWt6yDBQsWYMGCBf5/5+XlYeXKlRg6dGhEx40LPXa7GfZf/vKXOHnyJBwOBwRBwPz58zF58mSrm+XnwQcfhNPpDDJM69ats8UMBQC2bNmCt99+G83NzXA6nUhJScHHH39s6eB46tQpLFy4EM3NzUhLS8OKFSswePBgy9oTamMdKQAAAJRJREFUyJtvvont27ejsbERt912GzIyMrB161armwUAOHHiBKZOnYpBgwYhOTkZAHDHHXfgv/7rvyxuWQeNjY146aWX0NraCp7nkZ6ejldeeQUjRoywummKdCnDThAEQRjHvs4wgiAIIizIsBMEQSQYZNgJgiASDDLsBEEQCQYZdoIgiASDDDtBEESCQYadIAgiwfj/jSxTJ8g2hAUAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"tags": []
}
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "zj2K2dYBSKwd"
},
"source": [
"class GaussianMixture:\n",
"\n",
" def __init__(self, num_components, input_size, max_iter=10):\n",
" self.num_components = num_components\n",
" self.input_size = input_size\n",
" self.max_iter = max_iter\n",
"\n",
" # Parameters\n",
" self.pi = np.ones((self.num_components))/self.num_components\n",
" self.mu = np.random.normal(size=(self.num_components, self.input_size))\n",
" self.sigma = np.tile(np.eye(self.input_size), (self.num_components, 1, 1))\n",
"\n",
" def fit(self, X):\n",
" prev_elbo = -np.inf\n",
" for i in range(self.max_iter):\n",
" q = self._expect(X)\n",
" self.pi, self.mu, self.cov = self._maximize(X, q)\n",
" elbo = self.compute_elbo(X, q)\n",
" if np.allclose(elbo, prev_elbo):\n",
" break\n",
" prev_elbo = elbo\n",
"\n",
" def _expect(self, X):\n",
" n = X.shape[0]\n",
" c = self.num_components\n",
" d = self.input_size\n",
" pi = self.pi.reshape((1, c))\n",
" X = X.reshape((n, 1, d))\n",
" mu = self.mu.reshape((1, c, d))\n",
" inv_sigma = np.linalg.inv(self.sigma)\n",
" det_sigma = np.linalg.det(self.sigma)\n",
" distance = np.einsum('ncd,cde,nce->nc', X - mu, inv_sigma, X - mu)\n",
" Z = np.sqrt(det_sigma) * (2 * np.pi) ** (d/2)\n",
" p_x_given_t = np.einsum('nc,c->nc', np.exp(-distance/2), 1/Z)\n",
" q_unnormalized = p_x_given_t * pi\n",
" return q_unnormalized/q_unnormalized.sum(axis=1, keepdims=True)\n",
"\n",
" def _maximize(self, X, q):\n",
" n = X.shape[0]\n",
" c = self.num_components\n",
" d = self.input_size\n",
" nk = q.sum(axis=0).reshape((c, 1))\n",
" mu = np.einsum('nd,nc->cd', X, q)/nk\n",
" error = X.reshape((n, 1, d)) - mu.reshape((1, c, d))\n",
" sigma = np.einsum('ncd,nce,nc->cde', error, error, q)/nk.reshape(c, 1, 1)\n",
" pi = (nk/n).reshape((c,))\n",
" return pi, mu, sigma\n",
"\n",
" def compute_elbo(self, X, q):\n",
" n = X.shape[0]\n",
" c = self.num_components\n",
" d = self.input_size\n",
" q += 10e-20\n",
" pi = self.pi.reshape((1, c))\n",
" X = X.reshape((n, 1, d))\n",
" mu = self.mu.reshape((1, c, d))\n",
" inv_sigma = np.linalg.inv(self.sigma)\n",
" logdet_sigma = np.linalg.slogdet(self.sigma)[1].reshape(1, c)\n",
" distance = np.einsum('ncd,cde,nce->nc', X - mu, inv_sigma, X - mu)\n",
" constants = -d/2*np.log(2*np.pi)\n",
" log_p_x_given_t = constants - logdet_sigma/2 - distance.reshape((n, c))/2\n",
" log_pi = np.log(pi).reshape((1, c))\n",
" log_q = np.log(q)\n",
" return ((log_pi + log_p_x_given_t - log_q) * q).sum()\n",
" \n",
" def predict(self, X):\n",
" return self._expect(X)\n",
"\n"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 271
},
"id": "-d4KNRRaWBcU",
"outputId": "61b9b16e-a790-4e1c-b3b0-0e8cca1e3e4b"
},
"source": [
"gmm = GaussianMixture(num_components=3, input_size=2)\n",
"gmm.fit(X)\n",
"y = gmm.predict(X)\n",
"plt.scatter(X[:, 0], X[:, 1], c=y, cmap='gist_rainbow', s=50)\n",
"sns.despine()"
],
"execution_count": null,
"outputs": [
{
"output_type": "display_data",
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD+CAYAAAAuyi5kAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd5gUVdaH31tVHSeSc0ayoARBwQQqiCiYXdPqijkrJpQ1oWvENS1hAUVRQT9BFkEBQaKKIDnnnIbJ07mq7vdHTWq6exhkAIF6n2cemKmu27dqen731LknCCmlxMbGxsbmlEE50ROwsbGxsalYbGG3sbGxOcWwhd3GxsbmFMMWdhsbG5tTDFvYbWxsbE4xbGG3sbGxOcWocGH/6KOPaN68ORs2bKjooW1sbGxsykGFCvvq1atZtmwZderUqchhbWxsbGyOgAoT9nA4zCuvvMJLL710ROfpus6uXbvQdb2ipmJjY2NzWlNhwv7+++9z1VVXUbdu3SM6b9++ffTo0YN9+/ZV1FRsbGxsTmsqRNiXLl3KqlWruPnmmytiOBsbGxubo6BChH3RokVs3ryZHj160L17d/bt28ddd93F/PnzK2J4GxsbG5sjQByLImDdu3dn2LBhNGvW7LCv3bVrFz169GDmzJlH7MaxsbGxsYnFjmO3sbGxOcXQjsWgs2bNOhbD2tjYVAChsMlvK33oOnRp6yXJo57oKdlUMMdE2G1sbP6afPq/TB55exei8HvdkAy6uybP3FEDIUSZ59qcPNjCbmNzmjB1fi4PvrETfzB6W+3V/+6jSprG3ddUPUEzs6lobB+7jc1pwgsf740RdQB/UPLisL3YzdROHWxht7E5TVi5KZDwWGauQXaecRxnY3MssYXdxuY0IclTxp+7BK/bloNTBfs3aWNzmvD3PpVxOmI3SFUFendLxe2y5eBUwf5N2ticJrx8fy0a1HLidZeIu8clqFpJ48Nn7OTAUwk7KsbG5jQhPUVj6VfN+WRSJmO+z8YwJNddks5911WlcpotBacS9m/TxuY0Ismj8tBN1Xnopuoneio2xxDbFWNjY2NzimELu42Njc0phi3sNjY2NqcYtrDb2NjYnGLYwm5jY2NzimELu42Njc0phi3sNjY2NqcYtrDb2NjYnGLYwm5jY2NzilFhmacPPPAAu3btQlEUvF4vgwYNomXLlhU1vI2NjY1NOakwYX/zzTdJSUkB4KeffmLgwIFMnDixooa3sbGxsSknFeaKKRJ1gIKCArt/oo2Njc0JokKLgD3//PMsWLAAKSUjR46syKFtbGxsbMpJhW6evvbaa8yePZvHH3+ct956qyKHtrGxsbEpJ8ckKqZfv34sXLiQ7OzsYzG8jY2NjU0ZVIiw+3w+9u7dW/z9rFmzSEtLIz09vSKGt7GxsbE5AirExx4IBHj00UcJBAIoikJaWhrDhg2zN1BtbGxsTgAVIuxVq1bl66+/roihbGxsbGyOEjvz1MbGxuYUwxZ2Gxsbm1MMW9htbGxsTjEqNEHJxuZ0IxiUjB7vY+QXfvxBSa+L3TxxTxL165w8f1rBiM7nv25i9PwNBMIGfc9uwIMXt6R6qudET83mT3LyfPpsbP5iBAKSrv0yWL9Zxx+wfrZlh49PxvmZN7EqbVs5TuwEy4E/pNPtze9Zvy8Xf1gHYN2+XD6atYaFA6+iaY3UEzxDmz+D7YqxsfmTfPRJAes2lYg6QCQCeQWS2x87OZLz3p2+krV7c4pFHSCkG+T4Q/zj07kncGY2R4NtsZ8CbNki+e47CIehRw/o1MnOHzgeDB/rJxCMf2z9Jp0du/W/vEtmxNx1BCNGzM9NCb9vzeBgfpCqKe4TMDObo+Gv/ak7RZFSkp0NXi+43X9ehKWUPDkAhg4F0wTDgFcHw3nnSiZPPrqxbQ5Pvk8mPObQBHn5iY//VcgLRhIe01SF3EDYFvaTENsVcxyRUjJ0mKRWbahVG9LS4ZprJXv2/DkB+PJLGD4cgkHLWjcM8Pth/gJ44skKnrxNDN06OUmUXC2BMxodO7sppyDCo/9ZR1q/WWi9ZtDu3l+Z/GvGEY/TqWG1hMc0RVC/cvLRTNPmBGEL+3Hk1cEwYADs328JcTgMkydDp3MgJ+fIxf31f1lCfijBIHz6qbW5Z3PsGPRYCp44T0Vej+Cp+5NwuY7NE5MvYND5kd8ZNmUXeX4Dw4QVWwu46fUVjJi664jGevmq9nidaszPvU6NZy5vi0OzJeJkxP6tHSfy8yVvvBErxLoO2dnwZ8rXb9+e+JiiWAuIzbHjrDYOvhtdiVrVFZKTBGkpAq9H8MQ9SQx6LOXwA/xJPp2+m10ZQcKR6IXbHzIZMHwDwXCszzwRXc+owZh/XEglr5NUt4M0jwOPQ+WxS1rz7OXtKnrqNscJ28d+nPjlF3A4IBCIPRYIwDffWNb8kVC7NmzcGP+YYUC1xE/ZNhXEpRe42fVHDZauiuAPSM5u4yA56djaS5/P3Is/ZMY9JgT8sjqX7mdXLvd413VsRN+zGvDblgMEIwbnNKpGmtdZUdO1OQHYFvtxQo192o1C+xNL7IAnrQ3YQ3G54PrrISnp2G+erlghuetugy7n6dx+h8Eff5x494+Ux3cOiiLo0NbJ+Z1dx1zUwVq0yzxuHvn1OzSF85vV5NLWdWxRPwWwhf040a2bFbkSD68Xbr/9yMfs3x+uvcY6Xyn8TSYnQ+vW8PFHf36u5WXMZwZduhmM+VyycBF88ZXkgu4GH35UfldARSGlZMSnYRq2LkBNK6By/XwGvhykoMDkf1MjvD80xPc/RtD1E7/wHA2GIbmoXToup8Daoo1GNyTntbL7IJzuCHm8zZtD2LVrFz169GDmzJnUrVv3RE7lmDPiv5LHH4/2s7tc0KgR/LEYvN4/Z2EvWSIZ/zWEgnD55XDppZYVeSw5eFBSr5FBME4ct9sN61er1K9//MItn3ohyH9GRqLurdMpMZF43NZGtdNp3eOfvkuiTavDPEKVg/0ZBqPG57F0dZhG9TXu+VsqTRseu2zT0T/u4pnRGwiGTXxBa/GUpgCs++x1KbxwS2Oeu6nRMZuDzcmB7WM/jtxzt6B2LcmgQbBqtWVd33knvPTinxd1gPbtBe3bV+BEy8H4r2XCUD/ThC++MnnumaMXz/KwZ6/Jh8MjhEKlfyqtzUUB+QXWT0JhyC+QXHylj11rUoqjVnJyTRQFUlPK/wA7+7cAV961H8OQBEPg0OCjMfl8+FJl7rqp4tPwx8zYzSND18b41oUiUYWgepqTl29vwl2X16nw97Y5+bCF/TjTp4+gT58TPYujJzOLuNY6WNbxvn3Hby4/zNDRVAgd/qUAhEKSCZMj1KgBjw4sYP1my/pt11rjozeS6dyhbKs7GDLpd/d+fP6Sh92IDhFd8vBLWfTo5qFh3Yqz3E1T8tzoDXE3TIWAHu0r8cPgDnbHMptiKsTHnp2dzd13303Pnj258soreeihh8jKyqqIoW3+onRobz1xxCMlGbp0Pn4iI2Ucb7OQRR6KGPILYOqMMH1uyWXVOoNIxKrxsniZTverc1iyPHE2JsD/ZvhJ5MA0TcmocflHfA1lsTszSK5fj3tMSvh1TY4t6jZRVIiwCyHo378/06ZNY/LkydSrV4933nmnIoa2+YvSq6egatXYaB9FgeQUuObqihca05SYcSI+LuuuYR66XytFvL1FwNpsnr0gEjf01B+AZ171FX+/P8Pglfdy6XnzAe54PJNf/wixe79BKBJ/8HAEtu2KL8JR05OSSXOzuOShNbS+aRl3vrqJ1VviZJsBbodaZqSL22nHQNhEUyGfiPT0dDp37lz8/VlnncWePXsqYuhTkp9/llxyqaRmTclZZ0vGjIkvWH9lVFUwd5bKmW0soUxNhaQkaN4c5s9WKzTrctkKgx5XBHCk+3Gm+7mkT4DlK0uUvH49hTtudcSGfiaagoS9GYnFd/YCy2L/fWmIZhfs5Y2P85gxL8TYCX4u/VsGvywK4XTEH9zrFpzdpuxwQSkld722hVtf3MSsxXms3RZk7I8H6fyPVUyeF1sVslq6k7YN4yc8OTXBrd1rl/l+NqcfFb7Um6bJV199Rffu3St66FOCocMkfa6EmTNh/wFYvhwefAhuvuX4x18fLXXrCpYu1vh1nsonIxXmzlJZvVylceOKE/Ulywy6XRJk1hzTKnRmwszZJl0vCbJsRYm4f/yui5cHOqlaxcoJ8HoEvbpreNxW5BGA2wVeD3w+3JNw4xdAVayng6v7HyS/wNocBWtT2B+QTJ0ZItmrFIeYRp2rwh3XxRdh05QMn7Cfulcs5dPvM/AFSnzmhgH+oMltL20iFI71pY94rDXJHhW11Hs6NUGNSi6eu6lx4ouxOS2p8M3TV199Fa/Xy6233lrRQ5/05ORInngidtPR54Pvv4e5c+HCC0/M3I6Gtm0FbdseGx/vE8+G8cXxUPh8MGBgmJ++t7r8KIpgwCMunnzYic8HHo/1VLF7j8nIz8Ks3WDSuoXCXbc5qV1LofMIjQULY612RYErLnMyb2GIAn/8hTYYgu6t3KzZEiIjy8DQJZom0DTB96NrUDk9fjRQ/8Fb+fqnLPyhwveNe8sk037LZntWAUOn7GD3wSB1qrp54tqGLP6wC6+P28qPiw8WW+oDrm9I5RQ7ocgmmgqNY3/zzTdZv349w4YNw+ks34ftdIpjHztWcv8DUFAQe0wIuOMOGD3K3gQrQtcl7sp+jASJXULAoGcVzmyt0fcKDUcC90g8lq6IcP5VOfj9FG6ESoQCqckKi2dUYsnqEPc8nUVeQfw/j64dncz5tjqzfwuybnOEOjVUel/sTTiHlZv8dLlzDYGQCcJM6CZyeXWEJ0gwYlK8SSAkXrfCNV1r8tnTbct9jTanLxXmihkyZAirVq3i448/Lreon274fImzT6WEvNzjO5+TgbKsDiklr7wR4a4HAjRoVcCGjeXPeD27rYPffqjEJRdqKA4d4dBB1UmvbLBkVYh2rRxEErjhnU7o0sGJogi6n+fhgdtS6XtZUpkLy4SfswjrZslFxbswLUJI9ReKOljqL0AK/EGTCQv28cdG+0Nic3gqRNg3btzI8OHDOXDgADfddBN9+/blwQcfrIihTykuvJCEYXLJydC79/Gdz7Fk+3bJQ49GaNw8RKu2Yd56Ryf/CBtPaJqga5eyP6KiMAFp/wHJZf38R7QJXaeWwppNIVS15Jwduw3+8XgusxeEOedsJ86ocHRLkTUFHrrjyKo3RiKy1KIuSoYrfS3uUAJL3hL3YNhk3Oy9xT8N6wZGIkvB5rSmQnzsZ5xxBuvXr6+IoU5pWrQQXNJD8tPM6CqPmgbp6XDTTSdubhXJqtUm3S6ywgkjEQDJy4MNRn9q8PsvTlJTy+8yGfKGkwt7BePUnZeWS6PoOwnZ2ZLZ8wy6X1i+j/XQMT6yckx0o3C8QvwBeHpwHuvnV+eG+w+yeHkYXejo0noicKUqTJ5VwIO3pZU7frzXeem8P35/4YZpYSimKAzAF1bmqq6ULdKmCb6gwfQVe3j6yyWs2pWDIgS92tVmyK0daVrz2JUKtjm5sANgjzNffw0332zVU0lNtSI2zj8ffvv16MoK/JW46x6dvLwiUbcIBGD7Dhg4KMKDj4Zp0TbIOd2CjPpUJxxObGV3bK8y+wc33c4tHYViifqhmmqYsGlL+S3Y8ZMCBEMSMAsTmoq+TCSSTVt1fvqqOvUamqCUuHmyc02efTuTp/6VWe736toumQ4tknA7iyZtWeFOTaFtkyTWfd0OVS3795/iValSWXD1e3NYsTMHU4JuSqYu2805g6ay46CvzPNtTh9sYT/OuN2Ckf8V7NsLc+fAtq0wa6YgPR1GjpRce53krv6S+fPlSRf+CLB/v2TFyvjzDoVg6HCTkZ8YbNgoWfyH5NEnI/ToFSpT3Dt1UJk3w0Mw08sTjyg4XbGiDlaYYqMGSuF7Sb6aEOSF1wsY+kmA7Jz46fjFGaqlvwBfwCQQlEycUcDeA3qhVV+CPyD5eGwu+w/GOuIDQZPPpxzk6fd38tH4/WTm6Agh+OH95tx5ZTW8LgWPS+B1q/TvW51fR7WmUW0P3c+qnCAM06rL06CGhzG/riNwSCMNU0JBUOe1SSsT3kOb0wu7VswJIi1N0K6wQc2OHZIu50JenrXBKgSMHw/XXguffiJPqnTxggLLtRRKULhFymhL3u+Hpcslo8fo3Hd32fVV8vIla9cbRCKF/otDSEkR9LhIZfU6ne79cggEJQU+AMmDT+fToJ5gyKspXH2F1Zz5ystcrFgfjn0jYc1z3eYIPy8qSBj26HDA9Hl+bru6pOjXio1+ut+7jrAuKfCbeFyCZz7YydjBjbn64sp8/ExDhjxen8xcnSppGq5SWaMfPNiCzg8vxBfQ0aPWIUHvzlV49c7GdHtlZ9y56Kbku8U7GX5XlzLvoc3pgW2x/wW46W9w4IAl6mCJis8H335rCfzJRIMG4EoYFBVfIP1+GD6y7IiWcFjS9VIfs+YWva4ktETToEplwY8TvUgJPa/PJSOzSNShyBTfvlNyy325/HuYdaBZE63MRKWJPwTiJiGVpvSiq+uSng+uJyvPoMBvKXMgJAmEJLe+sIVd+61FxOVUqF3NGSXqAM3qJrF06Lnc0bMO1dIcVE93cmuPmqwd1ZXJr3QgLclRZpSQSJhqa3O6YQv7CWbHDsnSpfG74vh88N6/yz9WICDZtk1SkCD2+nigaYIXB6lxOzsBlusjDnl5ZY874X8R9uw1CYc5xFiXCCGZNcVDm1YqM2ZHyC/j+oNBeP71AvLyTdJSFZLK2NdwOAQ3XpFCcoLX6Dr0PL/kQqf9los/GN/Hb5ow8ruMsi4RgIY1PYx4vDX7v7mYfV9fxGfPtKV5vSTrWLUkqqe6456nqYJrOtU77Pg2pwe2sJ9g9u0rSXmPx969iY8VEQpJHn7EoGo1kzPbmlSrbnLzLQY5OSdG4B96QOWlQSrJyZCaAslJULsWOF3xKy6qKlx0Qdkfxf/9oJeywInyiWsOWPCbtTJu3mYQSdglyTrBoQlmzg3TvasrYZu5JK/g9uuS6HdpEs0aOQs7FpXg9QgeuzONalVKsky37g4RTvDeoYhkzdaSUCgpJfNWZTH4i028880WtuyNXwAsavZC8PGd5+BxRme2KkKQ6nYwsG+bw45hc3pgC/sJpkmTxP5oIeDMcvytXnudyejRVuSJz2eNN2ECXHChiWEcf3EXQjDgCY0Du53M+NHBgrkOdm518tD98S15lwuefrLs7R5PfEMVAEWAs1B4G9VXcGhlVP9CIrHqpycnKbwxMBWvJ/r1bhc0a6xxfR8PDodgzld1eOCWVFIKe8jWr63x/qCqvPZklajzGtdx4Uzw3i6noHUjq/xBQUCn2xO/ccULi3l57EYGjdnAmffO45mR6w67Yd77rDp8P+BiOjSyNlo1VXB1x3osGtybulWSyjzX5vTBbo33F+CWWyUTJsTWkPF6Yfo06No1sbtg+XJJ164m/jglaFNS4PPPFK666tj4Xnfvliz41cTjhkt6KHg8Zb+PYUieei7CsP8auJygG1CpEoz9xMkF3SwrdObsCO98EMYfgBv6ObjnHw4cDsHM2Tr9bvYV70OUxu2CbatTqF5NwTAk9dpmsu9AvI+1JexuN2xfWo1qVS27ZtK0AC++k8vaTTqpyQr33JLEwEdSSPLG2j2mKRO2HTQMSb3ey9mXGVvP3eMSbJjYljrVndz+1nL+b94+QpFot02SW2X0k2dy3fm1yryPpeciBCfV5rrN8cGOivkL8N8RkJUFc+ZYVrqqWv7bD94vW9QBZs6SMaF4ReTnw9SpJlddVf4WdaYp2b7divioWzeRb1lyz/0648abOJyWg8M04T8fqtx6S+KPlKoKhrzl5MXnJctXmqQkC85qJxBCYJomnS/288fSErGbt8DguZeCrFuSRPcLVS7qpvHzXD1qEfN64dnHXVSvphAISD4c5UdzWH73YpNF0a1JSnC7FO77u7dY1AH69vTQt6enXPenrF6yqiqY/nEzut+3nmDYpMBv4nUpSOCr1xtTp7qTPF+Eb+OIOljJR2+O31JuYT/WfW1tTl5sYf8L4PUKfpgKq1dL5s+3LO0+fShXhqbTGdvsoghFAVcZLoxDGTrM4OVXTfLzrcic+vVh+FCVCw/xfz/znM7X/2cSDFFc0hbg3gcNGjRQOL9b2R6+tDRRbKEX8fizoShRL6LAB+df5mfzyhQmfunloxEh3v0oxIEMHUWVVK+hULmKJDfXpMe12azbqBMIFiZ0KiYUZXMW3koTg3M7qWRmG7z731y+mFhARIfeF3sY+FA6jesfXUu7Nk297JjSjgk/Z7Nmc4A6NZzcdFllKqVaf2p7skI4NEEwQZOmrfsO72u3sTkctivmJGfXLskZzcy4fnqvF36aodClS9kLxNx5Jn+7xYi7UevxwOyZKp06WWLt80lq1A3HSfG3uLSHYNrUIysC993kCNfeEigzlG/7miTq1VVZuDjCZdfkEgxZTzVg1VivWhUO5hilrHlpZYvGuXS3CypVNcjKNQgVhrGrqtUkY+7/1aZty2NXxC6nIELtv82Ka7EDtG2UwtKh3Y7Z+9ucHtibpyc5desKnnjc6l5UmqQk6HMFlGpsFZc//pBcfkV8UQdrQ/b5QSUitG27TPiEALBs+ZHZCZ9/FebWuwIJi6MVsXadiZSS2+/Lp8BXIupg1XbZudeI3mdIEFYJ1sbp/gxZLOpghZvm+yT3Dzx4RPM/UtKTHVzZpTpOLfZPL8mt8tQNdtMMm6PHFvZTgMGDFUaNFJx5plUlskkTeOstwZdfKoX+68St9174pxG392dp5i8oObdKZWHFkifgcAk9pdF1yZPPhg77/gCtWiqs32iwZ198Szd2YUjczNowwNTjH1yyKkRWTvnL//4ZRjzWhpb1k0j2WCukpoLHpXD7JXX420Xl86/b2JSF7WM/BRBCcMMNghtuiP754sWSp57WmTff2pTt3h3efVulTZsSUZs3//AWtqOU27lmTUG7toLfFxVWKYxCkpkFO3dK6tU7/P7A6rUmodI1YuJXCqB+PUHdOiqLlkTQEn1iY84trAuQoPZKQoteMXnqrX34/Cbnd/RyW790UpPLv/lcHtKSHCz+qCszlhxkxpKDJLlVbrigFq0T9DW1sTlSbGE/Rfn9d0n3S4woX/iMGXBeN4NfF6i0bm0pnuMwnwBNgxuuj1bHy3sJfl8E0WpqCaWiwlfjdZ4ecPhNSCGIrTJQZHqXOjh2lBWx0qq5FuWCKY0iFFTNLGmOIQGlaPDS87dcSUa8ErlamIii8/mkIKYJP8zN5+WPMpjzRSNaNikjiwwr4WjushwWrs4jJUnjmgurUaNyYl+9ogh6dqxGz47VyhzXxubPYLtiTlEee9yIu8Hp88Gzz5WI2o03iIRWsAAqV4aXX4y2WKPdHtHtgMJhq0l3eWjdsiT2XSBBGIVt4wpL6UoToZj0f8iPlJKkJMGj93niJjm5XYJvRqVxVhsNtwtcHlkSEVM8R2ueqmbGxtwrhtVFiZIuV76AJDvP4OoHd5SZOJSVF6HTXYvp++wK/jlyC8/8ZxNNb/iVj7/dVb4bcZwJGzrjN/7O0798w5Bl09nryyEzmM8ry77lzElP0e5/z/DGyknkhO0ywCcrtsV+ChIISBYtjn9MSpg23SoJLITgxX+qTJqsk5UVXXURLNEf8q5KzZrRItihvUJKikF+nN6tKcnQqWP57AVVFXw0xMXf7gxgStOyqwUULxSF/+zZa7JwsUGXThqvDLRU/f1hATRNYBiSalUVxgxNoUtHjUUrgmzeEcAXtEr7yiJxl6LY/eJ2CR65O41P/y+fg1kGiiIJESKe915K2Juhs2RNkA6t48e63/ryatZu8xWXEwiErJEGDt9MuzOS6dY2vVz3I/p9JQcLgrgdKinuiovS2Zizn+6T3qYgEqIgEsKlavzz9+9wug0iaoigYX0IXlsxkeHrf2LhFYOp7kmrsPe3OT7Ywn6aU6uWYNkfGq+/YTD+a4muW66WQc+rNGsW30/eu5dC5SoCf0BG1VpRFEhOgWv6lf9B8Jq+DkwKLcND306U/LN5qyXsiiIY/EISzz7mZeVanZRkQesWKkIIbnsok0k/BvAHihpmiJLqjaV86v6g5N+jsjAVHcUJwUL/TrRTqQRNgT0HdDq0jp3/zv1B5i3PjVsjJhAyGfLVjiMW9i9+Xc8LExZysCCAKaFr01p8dOsFNKtZ9jhhQ+fjRb8wdPEvZAX9nFm9Fi+cfwk9Gp1hXZeUXDnlA/b78wsLK0DIsK495JPg0Yuf4YNGhP2BXF5YOp4R591zRPO3OfFUmCvmzTffpHv37jRv3pwNGzZU1LA2fwKPR9ChfeLjl1wSnYZeo4bg/fc09u12cHC/g8/HaAlFHSxLe+5PTlq3Eni9VqGvpCRo3kwwb6azuG5LeShu2VnGKbohSUoiqu5NcrLg3E4O2rTUEEKwcUuEiVP9lqgXb47Gd58YhsQXMAj6BcGwHtt79JDXB8OSlo3jW81b9wZwJWhiLSWs3XZkCUcj567hwbFz2Z3jI6SbRAyTORt20/X1b9mVFecRqRDdNLj8y5G8PGc623KzyQuFWLBzG9d+PYbRS38H4Lf9W9jvzysW9Rgi0fsiEWnw1dZfjmj+J4oJgSl0yLiEtH1NaHrgHD4sGIkuE2zInAZUmLD36NGDL774gjp16lTUkDZHwXtD4hfcSkqCN/919FEe9eoJli1yMf9nJ6NGOJg9w8mqpU4aNz6yj5SmJfbxW0lGOiFd57b7cqnd4iAf/dcf19/909xgiSJHRbzE7M4CZuHCVqp70iGIUq/XDQMlwS2rW91NKJLY/96odvlTfyO6wcBvf8UfjhYkKcEf1hkyfVnCc/+3YQ1L9u7Gr0f70/x6hCdm/A9fOMy2vINl1J8XYMb+7oJG5C/fyWtw/hD+nvsQy/RVFEgfW40dDMx/jeuy//GXn/uxosKEvWPHjtSqZcfg/lU491zBT9NVzjvXcpGoClx8Ecyfq3LmmUdeYyQclhw4IAu7F5VwVjuFa9YyrsQAACAASURBVK9W6dBe+dPFqK66PF4EjQRVRxSKdCAAmdmSga8U8PYHsVawpgkUUcqZUtpXH7V5WnpT9XDzLex/qoV4fXj8DK7GtT2cdUZynH6lEpc3wk5fBk3u/IG+Ly1gwero5KesghCjZ6/n3SkrmbduH6t2Z6Efmm8gTHCEiCgBxi1ei2HGj+P/bPlifJH4CQaqUJi1bRON06qRIJ2B6PtS6u2BA8HcRCedcPYbGfyr4H38MjoZwk+AWeH5zAmfHE8cFY0dFXMK06WLYP48jYBPJeBXmfmTRrt2Rya+gYDk4UfDVK0ZoHGzANVqBhjwdJhQqOIsoY/e9ZKeFidKJY417Q/Aa0N8BALR73/FJR6MYtUSJZGY4hBRL2pYfRgkEiFMFFcIieSLHw6SmRP/0X7cK22oU9VZnHCkqhItyQ9uPxt257PrYIAf/9hPnxcXMGraVgDGzN1Ao0fG8eTYhbz4zWL6vjOd2z+eg1m6zLIjBB4/OCIIZ4QsPZd2/xpNRkHswuY/dOf7EIJ6hHOqN6JucqVSC+AhOA4dQ+J2qOwOZJU59olkcmg6qoj/OOWXAb4MfHucZ/TXwBb20wCHQ6AlrFGeGCkll/cJ8ckYHb/fKitc4IPhI3SuuS62OM3yFQaDXgry5DMBps3QE2a7HkqN6grzZyTTuZNAVS0BVpTE56qK4I/l0SJUu6bKY3enWLXV5SHXWtyUoyhGvtQCYChxXfECioVOCMsv/+LH8cMX61RzsfbLLox4pgX39qvNzb2q4fIaRIxoC9gfMnhixHLmrtnHo2N+JRgx8Id0IobEF9LZeqAAI1goUqpuCbqg2H0ikWzLzOX2z76PmcNVzVrh1eLnDoQNg271GiGE4H9XPEwtbzopDstF5FEdKAJwBQot9tKLoInEpEFSxcXam9Ik18ivMP93RIYxZYJsZCRBmaDZwSmOLew2CZkz12TpMjOmTnwgCPMWmCxabIXESCm554EA53f38daQMO9/FOHGW/106uojN/fw4v7Lwgjn9shmybIwUuooql5m7RgpibtBO/jZNIa+WYlmTRxoioIQVv9Vqwdr4YCqXsrPbBbGzuvEuGscIas6ZCnGfp+ZcE5Oh8L13avz0RPNcXkjBCPxyxJoqsILXywjFOd4xDBBKrg0DbRIXH94xDRZsGUXu3Pyo35+e7uOpLrdqCL6T9qrObj1zPbUSrEabjdOrcbGW15jVPc7GNihN/8691o+u+wOvF7FKppW6sulaVxRtwNVXEefEWtKkzcPjKD22q7UXtuVqqvP4eHdr1BgHF2sfHfX+QmPJYsk+rgvO6rxT1ZsYbdJyA8/xk9yAqtL0/QZlvB9NjbCuG8i+AMlvVsLCmDdepMHHyu7EEwoJLnqplx8vsLQyegk0bi4XNChXeyOqxCCW65NYs3cWgS318e/qQHjhlbn7UGVeeGxVLzJRmEnpsINVFcQ4QyjaAZCDYMSBjUMzgBCjbUC4/UzzffrFASirc9cXyThwmQYJtsz/Al93R6Hyqt9O+N2JX7CcmsqO7Kjm8Smutz8cudDXNywCS5VJdnhJNnp5JFzuvHR5VdHvdahalzduD0vn9OPB8/szo2Nz+OZNv1wqw7cmgNVKCQ5XLSr1IARXe5NOI8j4b7d/+RfB4aRZeQSkTp+GeTT7An02PL3o7Lem2tNudzdAw/Rm9ROHNRWatLPffnRTv2kpMLi2AcPHsz06dM5ePAgd955J+np6UyZMqWihrc5ATgclhsgnkgpSkkNmXf+Hb+MbygE303Syc+XpKTEF6rvfwyXxMLHe8khNWA8HhjxXmqczcp48xdceWlJaNCT96RzRveNhAzDcnWoRrFVXOTykEgrOiTORqK71FPCnGUHGfCfNazbboUgnnVGKu891IaOLdK5tH0NfvxjH75gPKtdcEYdL/s25cSdc8Qw6XlmfaZvqcPMDdvjviakGzSsHJs0VDc1nak39ycr4Ccr4Kdearpl/ZeDgW2u4fbGFzJhx0L8eohu1VvStVrzCunOtDW8i3E5U2LcIiEZZmN4G1Pz53BVao8/Pf4X6UN5Lu81/hv4HABd6vRxX8bQtLdwimNXgvmvTIVZ7C+88AJz585lzZo1LFiwwBb1Q5g50+SqqwzattW59VaDJUv++mFY11yt4U4Qradp0LewM9PuPfF9nEWv2x+3TZ3Fzt0G4XDscXHI/90uuPh8Bz/+XzpX9S67bksiDuYYBPUIQjMQpUQ9+n2LHdqHIKlf21rJ5i7L5OrnF7FqSz66IdENyeJ1ufQc8BvLNuZy4wV1SfU6UA/pcOR2KnRtXYXnr2uL1xUruIoQnFEzjZZ10nmyxzl4nbE+c6eqcOEZ9amVlpzwOit7vDStXLXcol5EXW8VHmnRm2fbXE236i0qrOXejPz5KAkikApMPxNypxd/H5ERMo3sI7LincLJu2kvc6DGGlZWncP+Gmv4utJIqiiVj3ruJyu2K+Y48PzzBv36mUyZIlm9GsaPl1xwgcHYsce2PGx50HXJDz8YDB0W4aefjKgNz/ZnK/TpHRsP7/XC325Uad7M+vg0rJ/4Y6QbULNGYoFofoZa4i+PkygksHzki36uRN/eCs+/lkmPq/cyamwegUDiBSUeR9ZJrrTP3cThDuNJDXD7S8u5960VBOK4ZfxBg0Gj1uF1a8x9+yI6N6+M26mQ6tVwOxT6nVuHbwZ24eLWtXm0V2vcDhWtcFLJbo3qqW7GP2pZrj2aN+S5S7vg1lRcmrWAJrscNK9ehU9vveKIrvtEkyhqpQgNFZ/p5+GM56m2rQ0Nt3eixra2vJD5JhFZdrRPadzCTQOtHilK4kXvdMHuoHSMWblScu658Wueezywe7dK2qGhfseJpctM+lwZJBCwGleoKlSqJPhxqotmhaJtGJKPh+q892+dffsldWoLnnpS4567tWKLbvw3Ee5+IBDjjnG74abrHYwclrifqGFImrTNYvc+I27ZXrcLLrnYycr1+WRml9xHr0dQv67GnMm1SEstn30ipaRFz/Vs2xUu3DRNXNZXaoV9UoWJ5i1AUWQpv7hECgla6Q5NElQDVTX57PkOXHFOHVwOlR0H/OzPCdK4ZhJVUqOfNNbuzmbsvE1k5Ae5oGVNrjunEW6nZWWv2LOPyavWkVkQwJQqlT1eujapS/dmDY7Ykt6St59J238nbOpcXPtMOlVtclwbYO+O7KfF+p6EZGycfbLi5at67zE4/12Wh1dHuWs8ws1lnov4puaI4zbXUwVb2I8xAwYYfPBBdE2VIpKT4f33Fe644/g/OBUUSBo1CZBziKtXCKhZAzZv8uBIkCp/KFJKnnk+xH+GhzFNa5HwuKF9e5UpE714vWWPs2GTzqV9c8nJMfEHSz6ODgdcc6ULgxDfT/OXlOQtxOWEe/6ewtsvVynXPAGmz8/nhke2WRZ3nNZ5Hreg1/mp+CMRfl9VQIh8wmbsZqgsSujRrMga4Q6AkAgg2auhKYKJL17IOc2rlntuYDUT7z9+Iv9btY6wrmNKidfpoEPdOky86xbch6uzXHqOUvL0758xav1MTCnRTQO35qRT1SZMuPQZPNrx8z8/vectRmSNj0ok8ggXHT1n8kzNf3DTgfsokLEbNR7hZkGdybRxNj9ucz0VsF0xx5iMDOKKOlglbg8V1uPF+K+NmGqOYG2UFvhgytTyu4mEELz1upuVfyQz+CUXLw1yMW2Kl1k/Hl7UAZo11diyojJfjk7ljReTeGVgEhPGprF5WRVGf5zKlBmxog4QCsOYcYnrp8Tjsm4pTB7RiM7tkihSdU2FJI+Cyyno2yOdz95sxP8+bMGaSe2QSvzQS1Gcgi9LRL3Qd1QQ0MnxRbjqxdnk+kqs1Iy8ACNmreDNyYuYsXJ73Dj/Yb/8zuRV6whEwhjoSKHjC4dYuH0nL0ydcUTXOm7LfD7Z8DNBI0LY1DGR+PUQCzM28tTCz45orKPlzVpP8U6tZ2jgqI2CoLKazuNV72Rqo5FMCcyMK+oAhjSY7p99XOd6KmBXdzzGXHihYOJEiS9OuK7DAZ07nxg3zLJlRtw5gRWqOH26QSQsadpU4eyzy1dbplFDhSce/XMbm6oquKKniyt6Rp+fm2eWGdPu8x35A+f5HZOZ+1VTpJSEwpJZC/Px+Q26nJVMvZolVmxWbhiHphBOEJMuAKnqJaJ+CIYpGTd7G/de0YzP5q3mibFzEAKCEYMkp4Oa6Un8+Mw11EwvaVj7/txf8OvB6KgcYRCSBp8uXMK/+lyGo6yms6V4Z8Uk/Hpsgk7QiPDl5nm83fn2CrPapZQYGGgivqQIIehf5Qb6V7kh5piGikDELU4mEKhUbAer0wHbYj/G3HSTICkptheo0wnNm0OXLidmXvXrC1wJNVjy+dgw9z8YoPslPjp0KmDnziPbqCxi+w6Dpcsi5OX9ufNTUwSVKyX+mDZrevhOTYkQQuB2KfS+II3re1WOEnWAOtXdZS4qaSkaqiYTVpzxhwxWbstlxY4MnvxiLsGIQSBsWE9FoQjbMnK56cPo6LF9eXkgzOLwy9KZpyEZJCcQjPNOJazJ2skdMz+k6ecPsOHgXoRB3HwARQgOBI6+Bky+4eOJ7a9TbUln0ha3p/nySxmTMeGIim9dl9wHr4i/DyOEoE/SJUc9z9MNW9iPMV6vYP58ldatrWiStDRrU/H882HaNPW4bmKV5tZbHAkr/UlpJQ7l54PfD2vXmlx6mS+qbO7h2LRZp9vFmZx1zkF6XZlNg2YZPDYgL6aIWCK2bo8waYqPBb8Fef7xNKtUwCF4PYKXnznyJhblxeNS6d+3Lh5X7J+J163w2YvtGDHgLJI98a1Ul0OhUY0kPpy2NG6mqW5KVu/OZN2eklosDs2qOimEUfxlfW/59bMTZYwB8/es5eLv/snELQvZH8gtrAogUHQRI+6mNKnqPrqM0pAZpvva2/gk41v8ZgCJZGd4H0/u+Bev7v6Y1cEN3LfrWbpu6sctOx7mF1/87i9dXB242HMeHhEdW+sVHv6efANNHA2Pap6nI7Yr5jjQpIlg2TKN1aslu3dLmjUTNGx4YgS9iFq1BCOGafS/J4wesdwdQoApBeIQ14JhwMGDkhk/6fTqeXgLOSvb5MJLs8jOlkgJwcIN0c++CLB3v0HXczV8fknXLk7O7+qKWtzyC0xuv3c/cxcEcTitRSbJq3BTPy9fTfQXJ0UZBrz2fCWuujwp3hQSkp0XYceeEDWrOqlR9fBuiFfvPYP9WSG+m3MATRUoAiK65NV7m3HpOVXxBdN5fPgfcc8VQnBLj0Zc+/4KzAQWrENV2Lw/hxa1rZhrhyYJF1ZwLEoOE8LElApOVWNXbi7NqsfWbpFScvfPQ/Hr0ZEnRS4OYYIs9Gi4FI1rGnYhyRGbpOA3gszI+I0cvYCOaS1pndIk4b2ZmD2dbaFdMdEufjPIO/tG8e+c4USIYGCwLLiGafmzGVDtPp6t/mDMffq6xgjezx3JB7kjyTAyqaPV4um0B+mfenPC97dJjC3sx5HWrUVxE+m/Ahs2mDg0Ez1S8rgvCiM7Dg099Plh2TKTXj0PP+4nYwIEAjLGjREIGkye6mPaT1YbPq9XcEYTjanfVSM9zbKKb7tnP3MXBAiFIVjoHi4oMPi/SQX8Oq0OW7frqCp06+y26puUE1/A4LHBG/nup4O4HAqhiEmXdqm8/GhDqlZ2UK+mO+7Tk0NT+GRQW17qH2DesizcTpVLO1chLdlaYZLcGl8PPJ/rX5uHYZoEwyYuh4IiBP99rDM1K3k4o2Y6K3ZmxHXr6IZJg6pWHZcNGQcI6qGoRbXo/wpW68AGlSrFvb612bvICsXfSBYIhCkwVJNkzU2ztNq8d+6dMa+bvH8Od68cjEBgShOJpENaS8a3f5NULXYB/TpzKj4zfskIXeroetiKGsJ62vDLIG9nDOPqtF40d0UvGJrQeDL9Pp5Mvy/ueDZHhi3spyl79pi8OyRM6JC9NVGYlFNUB10WeuvcbkHVquVblH6cHooTt29atcWxooHA2vhcsy7CPQ9l8fXnVdm8NcLcX4KE4pQV9/lMnn/lIP/3Wa24Arx9V4i1GwNUq+Kg/ZnemNfc+OhqFi7PIxSWhMJW0+x5Sw/Q/c4DuF0K1Ss7eefpFlzW1QpPPJgTYvjEzUyasxshBNf3qEv/vo1JT4m18tufUYnHr23IqB834w9B24ZpvHVXR9o2tsIwH7z0LKYs20rgkAYaihA0qp5Gm3rWe366aFFCf70QUL9SGk2qxg/t9OmhmAJgFlZYZrLm4e+tunNZ3bO4qFZrlENeu7ZgC3eteIWAGf2B+D1nNf9Y8RL/1/7tBDM7MiJSZ2z2BF6t+VSFjGcTH9vHfpoydapl+UZj1SAvcsUIAQITgYlhSK69pnwblXEThhQzrk8/HIYZM4NkHDRYviqEHqd3qIVg2iw/AwYdAGDPvjBrNwTIyIxw7V0b6XT5au5+citX3r6Bdj1WsXKtn0jEZMKPGVx2+3LmL84lVFS6QJgINVL8RBIMmezYG+Tvz65g9u+Z7Drg57y7ZvLh+I1s3FnAhh35vPPFerr2n0lGdrTw5RSEuOCZb/lg8nIyfbkE9DyWbtvJ1a99z/YDVqGuTk1qMqhfZ9wOFYdq3Ztkl4MaaV7GP1ySRborJwejjE3HWzucnfBY68r1MKLK10qEI4jiKUB1+Qk7spl/YBF1ktJjRB3gg63jCMdJ4w/LCLMz/2BnYF/MsesqX06Skij5TEKcQmo6Ovv1jITXYVMx2Bb7aUo4UqrfaCFClPh2S35m/duvr0qlSuWz2O+43c3seaFDwikTC5bLKdix00BVi/qaxq8GJk0Y81U2v/yRzeZtIRyawB+2uiwZJsXNP3z+ML1vWU/TZoINW/0UBA2rRHuRu0nR4y4ygZDJoA820qiJICsvROly6qGwyf7sIH0G/EyH1sn06FCHK8+rzwP/+Zkd+/KiPFfhiEGOIXli5Hy+HdgbgEd6teeqDk0Y9+t6MvIDdG5Si74dmuBylKyunerXZ/qG9QTiJBgkOZ10adAw4T30ai4ea9eHfy+fjF8PIxxBhFZynYY0WZ2znZ7TnmbxVcOo6o4uIrYkby2GjB/W6VIcrPdtp56nZtTPr63ck3f3jmJzaAfhUqn/buHE0MJE4tzjJMVL16ROCa/DpmKwLfbTlO7d1UPk0xLFRJEyv/xS/qJMvXu56HaeM6rGjEjoZIBQWFKrhkL1qkoZPTkBTCJmkDXrg4RCkgKfVdvGiBNJ6Q+arFgdwhdTT6bsDkqrNxbww697445pGLBhm5/xs7bw4JD5tLpzHNMWby+8vlItVwHTNJm3ag++YIngNayWxrNXncO7t1zIDV2aRYk6wC3tO6AdGhcLaIpC3bR0zmvYMOG8AZ5tfzUDzu5HksOBosUuXhJJ0IgweuMPMefWdlVPOK4uDWo4YwtquRQns1p+zu1V++FR3Cgo1HZU5+36z3KGtx6OQ+xGgcAr3Fyf1qfM67A5emxh/4tyrCs9tGiu0rOXhidxGZcocnLKPx9FEXw7Lp1330jlzDYatWsptD/LGbdSpKZBWiq07byZ3v12AqW7+FDq/xKhhaPt/jIEWtdBLzIZj+BWCkXEPMnEIK3yutlxNitLa6kiiPGrl0Vlr5dJ/7iLqklJJDtdJDmdeJ1OmlWrxqR/3HXY0FghBE+d3ZePL/pH3IgXgKARZuaeJTE/v7/BdXjV2HMEgjqu6rRJaRp3vDQthQ8a/pOD7X8ns8MiNp01k7uqX8/Uxp/T0dsOj3CTqiSTJDyc4WzET43H4U3ovrGpKGxXzF+IUEjy1hsRhv4nQlYm1K0reOpZjXvudRyTePexn3sYNCjEiP9aNV4O3UgtTfPmR5b9p6qCO273cMft1h+xlJKBL+YyfFSBJbo6JHkFwZBBbl6gOGLECs4zKLE5SoRdiSm1K0is2iWLg0CULJRCgFSQxPr8hYBLzq3M+gM62/cmiBcvTB6Sceq1l54VQHqyi8rJCeoeJ6BD3XpseHYgszdvYl9+Pi2qVad93bpH9Puv7E5JWCYXIM0RG+FyadUu3FK7N1/snkrAtPq8ehU3TsXBF2cPLtei4ipV+7y6VoWZjcexIbSFDaGt1HXUpJ271QnL2zjdsIuA/UUwTUnvnkF++80kWCqixOuF2/6u8f6HsWmiC3/TGTY0zM4dkg4dVe5/wEnDRkf+EBYKSTIyJEOHBRk2LBRTpdHrhc8/S6J376NPP1+/IcL/TfSTly/ZvCXA7Ll56KVcuxIwEHHc7BLNHUISvUGIasZ1ySsKmFp0ezuJRCqAkKjOsOUyKeqYp1r10p/qX4f8YJj/TNhCMBzrxlGcAYRqIBUDXOEE1SGtReK9ey/ktu4tAMj2+fni10X8tHoDyS4Xfzu3A73ObIkax/VSFlJKftuzllnbl+FUNK5o2plWVRtEvebX/Su5esbz6DICCAwUZOGGaZLmZnjXJ+hTL37K88KcVYzZNZmMcDYXVu7ArXV6k+44+tZ4NscXW9j/Isz8yeCGa4Nx67e43bBslYeGDUtEYPArQf49JEwgYCWxOBzW11dfe7n0sj/3IJaZafDww36m/qAXd0+KRODVVz089OCRWZ7locVZm9m3P9ZVIbHa39Vv6GTT1khxUlL7szWWr8knECj1kRWFVRZLCazXo3DhucnMWZIRp266RNWCJCfrhAwVl8uJEKBoIQxZgKqaKIpASggLDRMFf9Cao3AEUTS9cBQJnmBCYe9xdj2+fs7aON2WkckVQ4bhD4cJFlYz8zqddGxUj7H3/b3ctV/8kRA3TXqdVQe34o9Y4Y0OVePKpl14/5L7UYTCmA1T+efiEQSNUJQzy0Tg0pLoVqMN4y4eFDcy5nCY0iRHz8WrenErf64mkM3xocJ87Fu3buXGG2+kZ8+e3HjjjWzbtq2ihj4tmPCtnrAolxDww5QSs3b5MoP33rXa0RUty5GIlf5/69/8xZme5SU31+SO23NoecZB5swKoBGmbRsY8q6HbVvTyyXqu3bpfDI6j1Ejc9m6tXzNERI9lQtAUwVfjazD7rVNWTSrIdtXNmXyFw1p08KDxy1KvVbg0lR6dEulXSsPvS5OY+zHTfhqaFNeG9AIj1sh2WtVYbTcOWFULUQwZCD1MG4lSIP6fsJGDqFIBH/QoMCv4wvoOKXOyIEdcCeFUNwFxaJe9L5EtBhdVwRUS/Mw/OHuxT976PNvyPEHikUdwB8Os2jrDr78JX6afTxenv85yw9sxh+xfGaGNAnqYb7ftJDPV/1EZjCXFxYNI1BK1IvvJ/BEm2v58qIXjljUpZQM3f0pzRd1pfXii2i08Bz6r3+Sg5Gsw55rc2KoMGF/8cUXufnmm5k2bRo333wz//znPytq6FOCrCzJ4JfDnNnCT6sz/Dz7VIi9e0usSbOMKrlSRvcdHT0qXJzkE49pP5Z/w840JVf0ymbK5BChEBTkSyJhWLkszKej/KSmln2+lJKBAw/SqeNOXnghk0GDsuh63i4euP9AwtoyBQUG739wgIg/hJMwKjqH+srTUlWaNnGQ5FVoWN+J16PgdCp8/2VTnn+iFo0bOKlSWaXnxan8MK4ZE0afwdzvWjF+eFN6dEtFCME/bqjFptmdubZ3Ol5vAKc7D4czELWghIwI23YXoMeZq2GYbNqey8C/n4nTESdaxXTy6JVtqZaq4FRNHKrg8k4NmfnGNVQq9K3vzcll1a69cUsKBMIRRs/7tewbXGqeX6+bQ8iIXTQDeoj/LJnMlB0LEiQpgSoU/BEfmnLklRJf3P4Or+/8gCw9h7AME5ZhJmdO59IVN+AzEteusTlxVIiwZ2ZmsmbNGvr0scKY+vTpw5o1a8jKsld0gIwMSeezAwx5O8LmTZJtWyXDPtbp2C7Als2WuPftp5JcRkevnr1K/iD37InfuAOsVnQHD5bfYv95VpgtW4yYhSIUgjWrI/z6S9nW92ef5fPZmHxCIUkgIAkGra9Jk3z8+73YYvMFBQaX9drEu+8dIC/PuggFEwcRRKH/3OMWvPtGjbgbbW6XwsP9q7P051ZsWXwm40c2pkM7b8zrikjyqjgcOoYMF2fTRl1nRMeIUxcdIBg2WbI2k4euac2r/TtSKcWJQ1NQFUGrhqlc0FHj85kL0Cgg3ROgSXWTZ68/mzpVSn6R2b4ADi2xmGb54qfkH0p2IL/M4/t82eRF/ITN+Iu6Lg2yQmWPcSgRM8KKgjWM3DsW/yGlA3R0Dkay+CZj8hGNaXN8qBBh37t3LzVq1EAt9BWqqkr16tXZu3dvRQx/0vPKi2EOHJAES1VcDYchNwcef8R6rL60p0qr1kpMSKDXC9ddr9Gkacmv6rzz1IRhigJod1b5rbLZP4fxFcQXNr8f5swu49EAGPJuDn5/HGs0IPn449yYZhJDhx1k565IlLuoKAZcQ+esti7GfVaH3j0rrm9l4/pJuONUaARwOlQ0NYGVqwjqVLfKE/Tv04J1Y29g0Yh+rPviBurUCrBs825CEYP8QBh/KMKug7nc8OqXZOWXWLENqlZC143iDN7STyYCSd2qTh7/ZiQDvv2EORtWYRbGWmb58/hk8Y+8/vMXfLd6Pl6HK2EhMYBayZXpWLUFTjV+dnCS5qFbjbaHu1WA5Ut/e+tHtFxwHr2X/o2QGf8z4DcDTMr8sVxj2hxf7HDH48D4r/S43YpME2bPMgkEJB6P4McZbgY9H+aT0TrhEKSkwmOPOxjwdPQf6+13OHjzX7H1WBwOaN5CoUOHw6/XWZkGhmFtzGqaFX54KJoGnjjlckuza1dit4/PZ1KQL0kt1dP1y3HZxRmih+J2Cb78pA61a//5GuvxuKF3Xd4asSHuMYdwoDgiBMOxj0AOTeHWPiXx26qqULd6Mht2HWTJpj2E9dhzIobBV7OW82Dfc5FSMmLmz2imD4FZFI1PBIe1ZH9iLAAAIABJREFUmZmcw7bcHNZlWh+OH1b/Qds6Dbmx87k8++MIQBLUIyQ5rAJlUo2AEpuZ69GcPNS+L52rt6Zpal3W5WyLstxVoZDq8HJlg27lul/PbXyVr/dNImAGMYUs0/xL1FjD5sRSIRZ7rVq12L9/P0ahf8AwDA4cOECtWrUqYviTnuBhnraLLHmvV/Duey4ysrzsO+hlz34vzzznRFWj/5CrVFH4YUYStWsLkpMhJcVqjN2+g8J3k2MLYJVm4a8hzj9nP62b7KNts31884UPNUFMtpQwbmwWZ7fawf39D7BubazlVqmMJhiaJvAmRc8lEJMJGv16n8/6DGVnRzhwIFQhiVpVK7kY9a8OeNxqseXudaskJ2mMe78Lo1/thtet4Sr0o6uKwONSefrONrRolBYz3sqt+1CV+Pc4GNb5dc0OAD6cNoORs37GKCrDCyhIXIRxu/1omiXcRfjDIZbu3MKAiSMJ6uHiY75IkIJwABEUCFkS119k/ber1ohbWndHCME3l7xO1xrtcKlOUh1JuFUnbSs3ZcrlQ3Cphw9XPRDKYNzeiQRM60NZRg4YSYqXG6v1PeyYNsefClluq1SpQsuWLfn+++/p27cv33//PS1btqRy5dg05NORszsoLFoYX9Bq1hSkH9IrQlVFmf52gLPPVtmwJZlfFhjs3ydp2UqhVeuyXTBL/whzbZ/MqHDBnTusGi1uFwRDlrAJTDRFR5EG/8/eeUdHVa1t/HfalPSEhBB6B6mCIKIUFcGGih3sHa9evWK/dhHLtXMt2MDCByICSlVEAemISJfeQyrpyZTT9vfHTBozQ0ACojfPWlkL5pyz9z4nmfe8+y3Ps2t7wNB+O9Vkzswyxn/VgH5nV8aB7rgzjv+OLgqpxHE6YdiwGFS1ugE8s3c0c74rDtvdqShQUqJzySWb2LK1DFmWqFdP4+mnWnHhhaEc5EeDAWfWZ82MAUyde4A9Bzy0axHD5QMbERMd+AosG38xX8zYwdqt+TRpEM3Nl7amQ6vwIh7x0eEpfiFgvJPjo/AbBmN++BF/uK0QAkwbM0zG3G8aYCkQbtMiQLNshCqCJL6B/IRt+yrWk+iM5avzRnGgLJfdJRk0jEqhZVzDGp9POZYX/Yoma/gtPXg/ErIlsBWqbRSckpPW7uYMSuxPnpFHgpqAItVJ2J0sqLV91HPPPcfjjz/O+++/T1xcHP/5z39qa+i/PEaOcnDFpb6Q0Ik7Cka+9Me7SmVZok/fI/8VPv90UfUa8CAsC9IaQHySwu4dHmw94ClWdZYtU+C14J47c1i/tSly0GMdMSKRZUt9rFnjr9AfjY6WaN1a47nnQylmH3qwPj/NLwlZh9stcdUV0Qy9dg2GYVMuOJeRYTPiwS0IIbjoosh8JkeCxHgHd1zTIuyxBsluHr2t8xGN07dz84h9nS6nxvUDuvHpgoURjDqAFDDeEXC4TYokJCRs5Cqxeq8R2jKc6k5kf+k+9pTsIt7hop7ryJwsh6RRacEDc6hBPh5LlZBkiQQ1nmtSLiHX2keXtd2QkHDLLu5J+wf3NLj7D9XI16F2UdegdIIw41uD++7R8XoC9duKAqNednDbnbUbTz4cUuMOhI2lQyCePvr9WJ546GBIMrSi0UWSiI6RmDYzjW6nVWZ5bVvw889evplWimXBJZdEM3BQVEgIqRzLlpVy/4h0cnPNinPO7iOz+OcsLLs6UYBNgD+4YUMnS5f0Omla0hes3ck//jsdw7Qwg4xhUU6NK/t24rmbBtD7iSfJ8/gjFusLyYLYyLQFIjr0mECAy6hGh+uQNW7vNpjHz7yp4rP5B5bwyPKRWLYFEhiWwZAWF/J8z0dQ5cM7AmVmGZ2W9cVne1A5hD8eF++c8joDU87m3I3nkeHPxKAylOSW3QxLHsqoZs8fdo46HH/UZT5OEC4dojH4UpWNG2xMEzp3kdG0E2ukFCV8krQcY/4bvsKlwtAKgSJLlAaraDIO6Lz9n0y+nxWofjl7QBwP/TuNlq0P39B05pkxrFrRjh07/Ph8grJSH7ffth67SgNp+ZxBVhfy8gyysnTS0k6OjsczOzTlg/sG8+3yLWxLzyM1KZZbBnWnX+cW7MjKwoxUjwogBJJqBjzzQ/4EXKqGoXkJ92sq56iRMdEkP7JkI0sap9ZvXnHOpvyt/GvJU3it6qLXM/bMxa24eKrHiMPeV7QazfDGN/H+vndD50cnR8/g27zp5BoHqxl1AK/tZULuRP7V8J+kaMcWOqvDsaFuz3QCIcsSXboqdD9NOeFGHeDCwW4iUZO0baewNUxy9FDoOnTp6iB9v86F/bYwdVI+xUUWpSU2c2YUcsmArWz5vebabEmSaNPGRefObr74PB1/SOt/FZsnBEKIP+WZHQrbthkz5QfOvnskT7w7nmW/rCRRKeXpYX3p36UlkhQkKBYCBSu0uyyY/JSxkKyAyLRTVXFrDhyqyrCe/WgQV9kpG7wICcHgjj2JVi3cShmKbCJLNuDnifkv8+GvEwAYs+kz/FZoaMZr+Zi081tKjQjtzVWw27M9LM2yjc3ove8xM38WHjv8bkOVVJYUL6txjjocX9QZ9v8hPPNCHLFxUjXlJFkOVOMc2F1z84o7SuLGW2KJT1B4deQBSoqtao1Stg1lpTbPPrb/qNa1a5fnsHFlgNato0g+AvHpmpCTW8rbHyzmmlvHc8s9k5jx3SYM4zDe9SF4+8vv+HzOIrx+HY9PRzdMft+Vzs3Pvc/BwsAzbJmaSrTLhQyoWAG5waCBl7BRJANZklANhabuZJ6+6FqeuXgoix58mXpxNmW+fOIoIZ5i4ikihlKiKGNv9g4cSnlSsxI+088Hq8eTXpzJ2oMbg4WVoVBlld3F+0I+z/BmMGX/FL5N/5Z8PZ8VhSsDYZ8wMGwDI0z3azkkJNS6JOqfjrpQzP8QmjVXWbi8Pq+8UMzsmT5sGwYMdJKWajF5vI0CWBxS/hCEosAtd8TzzMhAEu6HOUURectX/1KG12PjPkKx6ZYto9ixPbJxd7lkRo1qe0RjHQ7bdx3k9n9ORtdNDNOu+OybWRv58K0rcTgO/3UoLvPy1bxl+I3qgRIB+HWTL+cu5b5rL0CWZZ65+ioe+WI8PsNApTrDJKqNLEk4VY3Xr7uNni3bVBz/eu132Ja/mnCHErx+b/FeJHf4h2QLmznb5xPviCPLG156zrRN4p2VHBGWsHh6w9PMPDATRQoIr+i2jl81IcI71MLi0nqDWetdS1kYr90UBv3j+oW/uA4nDHUe+/8YmjZTef+TJPZmN2R/bkM+m1iP3Cwd3R/Y7kvVhC4CP8kpMqs3N+P5F+tVJDvD9OZUQ2Tt0lDcfmdTnBE6Q+PjVL75pjvduwUMkmnYrPsth9W/ZOHzhk8YCCHYsD6Tryat5bs5WygrDXi5T436njKPXmHUAXw+k207cpk2c2ON69y0az+aGt7466bJ4rVbKv4/6NRTGXPXXXRo0hgJm2i1lCTHQZK0XJKkEno2aciUBx6rZtQBir1FQPVXa0XewQx4zOFg2CbF/lJuaHsV7giCGc1jm9A0plHFZ+9vf4/ZB2ah2368lgeP5cEUJrIhkI3wv79W7pYMqz+Udu52OKXq+Q637ObhRg8Rp9ZAMFSH4446j70OtGjlQNMCDJEqNgIbi/JGHnj7vVTS0qr/qZx+RjTLFocqCAE0b+UiNu7It+OnnRbPgw+34I3XdiOEQNcF0dEKSUkaX03pTmpqwIDMnb2LV19YGaApkAL6qHf8oys33d6pYqyiQi/33z+dfXsLME2Bqsq8/PJ8zhnYnJ27D4ad3+c3+Xr6eoZeeSoQeDEUFXtQFJnYmMqafZdDO0zDlMB5iMd/Zvt2TG79ILe9/SR7c4rQy18owmR/xnaWbfyFNg2GVJzvNXzYthm2lFICsCBKc+MxQnMYUZqbbmmd6NesF3P2/cS6vE14zMB5DtmBS3Hw5pkjK843bZPPdn+Gz/aFjCUhoRkyftWuqOpRUHApTl5r/zKqpPJ1+694N+NdvsidQLFVTAtnCx5uNILBSReHjFeHE486w/4/hOJCk2ULSrEtwel9Y0hODZRaXndrAp99WIBhlCsOURE+iI1ROHtgaLfUY880ZOil20Pq0V1uiadeaBRyfk24/Y6mXHhRfWbPyqGwwKBb93jOObdyh/DL8gxGPbMMv6/6VuGT99cSG6dx+dXtAPj349+xc0ceZtCIGoaJpBr8MG8TQlYjlh+WlgUSjit+3cw7H88kKzsfIQTNmtTnsX9dQ/s2TejcuinqIYRequTHpZQhSzbp+wr499v/5f7rh5GWEqgKmb9uBekHs9DN6p62z/Az9oepXHHWQGLdAUWjvNKCwz4jCYnU6BT2Fx3ArCI8rUgK9dwJ9G92Bqqs8Ok5bzN3/wIm7ZhOqVFKv4ZncmObK0l2V/YV5Ov5Eb1/CHScKpbAUgLzdontwOiOb9PcHRD1cMsuHmn8MI80fviwa/4zYAubHDsLp+QiUf7fbJKsC8X8jbH+Vw+fvZPL15/m8c6LmfRptYkn/rGPp+7dzzmn/M5Ljx1ACEGLVg5efCsVp0vCEdxdu6Mk4uJkvpjWJKR7FKBr92g+/7o1bdu7cDglnE6JJs0cvD+uBWcP+GNb8YYNXdx5V1MeeawV5w1MrlYHP2b0mhCjDuDzWXz07lpsW5CeXsSmTVkVRh0ASSDJwQRmBEiSxKmdGrJ81Waefmk8GRnZ2JaBsE327M3g3odGs3rddlRFYeTwq3E5NCRAk3xEKSUokh1UY7JZsW49dz33AgcLAkb6+9+W4NXDaw4qsswvWzdU/H9d+sbDfiET3LF8MeQtuqd1wqk4iHVE45Q1WsUnM6BJJ77e8CV5noOossrFzQYyfsC7fHPBZ4zoclc1ow4Qp8VFTJCWQ7PBZQhchmBg4jkVRv1kxtdlE+ma2YqemR3olNGcC7P7sUnfUPOFfzPUeex/Q5QWW9w5ZBeb13uxTEACPSjz5q/S+j953EHqp6nc8UAq19yQQJ+zo5k8oYiM/QadTnVxxbVxxMRGDqmc3juGH5aeQm6OgW1B/QbqcWsg2rYlMgV0SbFOYYGPvXvyUTUFv98iUFJoI8kG5eXiigiGmA5Zo9OhcNsNPXn2P+Mwg12cVc+wheCpF8Yx66sX6d+9A588NZwPp81j/aZQLnVbCDw+H5O+n8s/hw2tYGsMB6/u440po7GsMgadNgAhBE5ZwWuZoTsLIWiT3Jh6UYl8NuQtDhRn8VvmGt5f9jolnn188/sWHIqDD1a+wyN9n2DwKYfncHEpLgamDuKHrLkYorrnLgiqUgXhVtz0SgwvpXcy4cvSz3m8aAReUZnUXW38wuDcc/gpdSUt1VZ/4upOLOo89r8hHr9rHxvXePF6AvFq3R/euHg9go/fyKmg1m3YWOOBx5J59d00broj8bBGvSpS6mukph0fwe1yuNyRfRDbBrdbJaV+DJYZaOCJppRoSlEks8JGqsJCEXaV2nJBlEvjtZGDSa0fRW5uwMs+9C4kwO/X+WV1IDnaqVUTRgwbhNsZvmvYtCx+/nU1AOed2huXI3JTlddXzFtT3+X7VfPo2bwbMjYyVdYY/IlWVC7vemHFdQ1jUxm/+gOKfAV4jIAh0y0d3dJ5fckr7M7feZinGcAznZ6hUVRjopRD+OwlEWSRBIfkoHV0G3ok9qxxvD8TlrB4vviJaka9HD7h5e3iV/6EVf15qDPsfzMczDb4eW4JegRq3EPhKbUpKjiyOm7db7FsXg4/Tssgc9+JVc4ZPKQ1WhgVI1mWOL13Gu4ojTZtkqkXrxIninGio2ESY/mJMnzIdiBcogkLp22g2SYxiuDjt6/gzF7NURQZIQLnSMJGKn8BBCGASV/PxbZt9KAqyeGecHmDz8BuZ5Ecl4imHPpiEqgExD/8hp+P5nxKckwSg7sMIlpzoGJV/EQpEmnxKQxo37fi6vVZa8kpzQ4bTjEtg683TKrxmSY4EpjZdyYjO7/AwNSBDGowiK5JXdA0lRg1BofsYED9gXzW44uThsohEnaZO/CL0EQwBEo0f/LNPcEr+nNRF4r5m2Hfbh2nUzpiw44E0TE1v98Xzc5i1D3rALBMC8uwadc1ntcm9yIu8dgbh2rCHf/oypKf95OT7amItTscCtExGo89EwgT/DR7Hf6cDJRyrtmgLdKEjWb6KFWdWHKgXtvtkOnSuRHt2weopWNjokiMcaIX5lUQbAnALzkwZA0Jm127NnPDrXfj9XnRHBq4wt+3pqoMOON0AFwOB2P/NYrR079gzq+LKC8h1SQ/qlQZAvHpXjLzsnjs/PtIi0tl/PIv0Q0/tiRz3iln89B5/8ChBuZbtWcpT3/3MLrwhxXTtoTF7oLdgXsQgvXZq1i45zssYdKn6UB6NOxTQdTlUBy0i25JcWxnVEmlX6PzcWkxZPuzSHM1JF4LpS0OB0tY6LYPl3x42ujjBU3SsIkc9tKkE8fJdDKgzrD/zZDaUEPXwxj1MLwkqgaDLo3HEaGGvBzbNhTx3F1r8XvNQCt80PBt+S2PIe2/44XPTuesC48v935snIPxX1/C9GnbmTN9J6Zpc+7AZlx9XXsSEl1kpRfw/iuzK8sRq9xr+T+jTB0jNhZsk3ZNwFG2mxceeZdBl/alYbNUZE8RSjWFI3AJHWyB7PAhy148wUpD3a+j2TZ+TanmM8uyTExUFNecP6jis/joGJ657h7Wb19KXkkw3HPI78K2BZqqsSNjC6vWzCDOk48syThUJ73T2hDnjgVgT94Onp31ID7Lh9CqtjFVQpEUWia1xLQNnvzpbtZm/4Lf9CIQ/LDzW1oktOXN88ejyAr/XnYXa3KXY9omkiTx7vpRJLtT6Va/F5e3vJHO9Q4fgik2inhr1zPMyf4aUxgkaPW4o+mDDG105wk18M2UFqTKaeyxdoUcc+DkyqhhJ2wtJwPqDPvfDI2aOujYzc26XzwRdVEB3NESyfU1nny9ZkbN/xu9E8NvVRj1qt9XyxQ8f8cqPl5wDs3axtbCHURGVLTGsBs7MOzGDiHHvvvmV+wgy2IkTl0FwbVXdeTHb78nb79JbjC3sP7XLahOBUsJbXgKGHcDFG+IMVYNk0RFwZVSj8z8fDRVoX+PHgy/5ioSw6iAn3faOUxbMgPTCp2nQVIquunhiU/vwVelTt1v+Pj4+7cRQnBhzyv4ctU4DEtHFoFQuECEvCVUReXqzkP5cuPHrMlagb8KIZjX9LAj/3c+XP0quOC33GWVx4NvqGzPAb7fO5WFGbMZ0uJG7u/yXNjnqdt+bl5zPunePRgiEJ7K03MYvWskWf50RrQaGfa64wFJkngrcQzX5Q3BKyqfn4ZGPbke98Q8cMLWcjKgLsb+N8RbnzcjJU0jKhhikWUJp0uiz4BYzr4glj7nxfLUa42Ztao9ScmHf7cbusX6RVmotg8NPZDYO/Qcw+brD2pO1h1PZOzLr17mGBaChbPmofuNEC1Ww28iIjJfChDhvyrC56dzagMWjPuYeR9/yFPD7yQlMZGCwjx2791OaWlxxbnDzrmahOh41CrxdkmScGpOHrrqfr78eSy6GVoa6Td8fPHTB1iWyabMddjBXIDTEBWEY+XL1GSNx/o9RfPElkz9/fNqRr0cuq0zZ8fXzNg1MexxkEBI+Cwv3+7+P37LDU/qNS9nOlm+AxVGvRw+28OXBz4mXw9PbXC8cJarPzNSfuJc5yCipWgSpSRuib6Ln1JXUk9JPqFr+bNR57H/DdGgsYN5G9rz3bQili8oIS5B4fLrkzila6CLctuaQlZ8l8Xkt3PpfVED2nQNrxR0MMPDAxf8hC/fi6MiTWdjCRmzSszStmDbukKy9xaTsauIemnRNG1/YhtDWrZrwK/LdqAbRtiwE0BUlIynLHyCTUJCNiQsJXJnaURIUkXYobAonzEfvcy27RtRVQ3TNDitex/uuOVB4qPj+WjEu0z4aRI/rlmIYRp0bdWZWwfdQOtGrXhl8i/YIvzLybB0MgsOEO9KILMoHQh4ZS5dYMlgS+CQFEZe8ArdGgXCJ0X+yCWipm3UWMcO4Le8TN35Kd1Tzgw59n3ONLx2eLZIVdJYWbCIC1OvrHGO2kRXR3cmpcw4oXOejKgz7H9TOJwylw1L5LJhiRWfmYbNyBtX8duCXPRgAvKrt3bQc2B9nv68B4pa3St98bZl5GV6EXZlRyoESKlsYWFXYfHLSy/i/r6T0ZwKlmGT2jyOxz8bRIPmR5Z8O1acP6Q7Uz5fGtGoOxwqzVq42bs7shepCIElQkMbAFIE8U+Xy0m/s3oDYJoGz794P3n5udi2hRHsNl29ZgmFRXk89dibJMTEc+9lw7n3suEhY6lq5CS0bds4VAeXdR3Knvkv4gvSBUiAaoMq/ERJFm/P+idC2JzSqAdJrmQOenPCjudU3HhFzRS+AkGuLyvsMbkGFseTu47m7426UMz/ECa+to3V83PxeSxsO1D/7fdarJqXzaS3dlQ7N2tvKTvWF2Bb4YU3FCoD+LIssIoLMPwWnmIdv9dk/9YC/n3xdPyeyG3rtYmk5Fieev1a3G4nTkf1Coj4hCjufXwwiYnaYXTnAvzpzirsZpIUCGO1aN2IaMUOc60gyunkjNN7ALBq9RKKS4qwD9EyNQyDXbu28H8T32LsuOeZPeczioryQlYwoOtFaBEEp1MT0qifkMaA9hfRtUkPXFoVDhssXOgIYWFYfkzbYGP6SsyiYpxhCMGciourTrmF85pcikMON5+oULFWJY2OSd3DrunC+lfglqPDHjNtgzOSzgl7rA7HH8fssU+fPp1PPvmEnTt38sQTT3DDDTfUxrrqcAQoyPbh81jUb+IO8bYPhRCCaWN24feGEVD22kx9byfXP1JJjZub4UVzKOhhBDAAJASyDKomo4lSZLt6bFjYAp/HYMm3OxlwXfs/cHdHj269WvF/cx9i+YIt5B0sISk5hvZdmtCgYQKyLJOXtY8Nv22L6JU7MMAEd0oipmXRqXNbhl4/mO+/+hx9TyGq7KBUdWMFSwVdto5alENm+h6K8nNZuXI+fn94kRHd8LPg5+lIwoemOZg951OG3/UiXbv0qTjnirOuZ+GGuRSUHKzgPJckCYfq5P7LngRAkRVeuGQ0S3cuYNaGKZT4CsnJXcuh718hbFSvSYvktuwq2xlIE2AjSwpnNjmXm0+9D7/lY2fRZvaV7sZrlnvvQVbPcsMuq1zd6vaw9zQg5RI+3f9fdpdtC5ReBuGSo7i5yT9J0P43eVpOBhyz5um2bduQZZmPPvqILl26HLVhP96apxtWeZnwXgHpuw06nubihn8m0qTl8a+7Pp7YvbGYt/6xll0bilEUCc0lc9PT7blkePOIJWa6z2Jwg9lhPfByzC24pOIFkZfp5ebuszAidK3G13dz+sVNiY8ymPfZWkw9/Hn9r2rDv9479yjvsPZg+HRK8ouIToiltNTLiOufRjcMLKlKAhOBCx8KNrIs039Ae/Zt3U5sQgKupGg2rP6lmrduE9BiNRQJt0NHVQWaw0GR0PFEjE4IVNmPKldmaB0OF6++Mp2YmMpwVam3mClLxjN/3RwMU6dLix4MO/sOmqeGb4ffmrmG56bdglcPz7SZltCcZ676jKX7f8QSNr0a9ad5QuvKexE2K7N+ZmH6HFbkLKBQP4giqciyjCzJvNTrY3rU7xt2bACPWcqYPa8wLXM8ZVYJjVzNGN7sUS5pMPSkb2r6O+OYPfa2bQNenhxJc+1PxIcvH2TMi3noPoFtw7oVXr76qJDRkxvR/6JQxsK/ArL3enjw3CV4SgIGwgB8HotPnvwdQ7e48v7WYa/TnDJRsSqlheFDI3H1HNW8/nppbk7tl8ran7MxDjHaTrfCM5/1ptMZKcyftJUFqhzWsMuKRFy9w+uf1haKcwtZ+91ySgtKaNKxBW3O6MTsdyfy66xFSLKEbdl0GdCLex6/mQ9GvYdAx0YO7DwCktmBNaOz9uclCNvmYEYGpivMfUFACUn1AzaWCZZpoMiAOzKDpCIdUnYjBCt/mcuAc6+p+CjGHcctA+/lloH3HtF9K7IaQiUsCxuHbSELm7KCPWzcMZ8hXa5HDRPmkSWZ3mnn0DstEDbZWbSFbYUbiHcmcXr9fqjy4Rt7otQYHmo9iodaj0IIUWfMTxL8bZOn2zf5GTMqD18VWlnDAMMQjBh6gKVZbY5Y4edoIIRg3sRsxr+yl6y9PpIaOBg6ogmXDW+ELB/7H/3kN3eEZTn0eyz+b9Q2LhneAocz1G2UJInL727J5NE7QsIxTrfMlfe0rPZZyUEPN45oTX6WhwO7yjD8FqomI2zB7c91pdMZAVraXhc256PHFoddq6rJnDu03R+91SPGiikLmPHq/4EAUzdwRrkwTRNbM7GMyhfZuh9XkL0rnfMuP4+FMxdihSh7C5yUIYLEXbYiIiZjJdlCkqobfdUGh2Gja8oh1wgUSQ+x97rhJy8v8w/fN0CLlA44VCe+oJapYlu47MB9SYBt+fly4ShWbp3Jv6+ZhCIf/ivfKr49reL/WOiszqifPKjRsF9++eVkZGSEPbZs2TIU5eTUN5wytrCCX/xQSBIsmFnKRdfWvtLL+4/tZPqHB/B5Al/6zN0+xjy+k3WLC3luYqcarq4Zv3yfjRXhvgD2bCyh7Wnhyxevf7Qt29cVsfbn3Mq2fJdCj3NTGPpgQMmnONfDp//4ji2L96M6FIRu0ad/axr3akFskpM+gxuTkFLphUfHO7n79X588MhiTL9VUR/ujFK5+I5ONOtQL2QdxdlFfP/aHNbNWotlWDTv0ZwLHxtMs+7Nj/p5ZGzZy8zXJmD6Kw243xMoaRS6qCbxZhkmB/dnMWj41ZQVl/Lrol8Dhl8IBBIuuQxNPsJkr1xdoFoTFpIQyH4Fny1Qk5KQHBoxUdEU5O/BMkPHdTrdNGx4bIwxhcwrAAAgAElEQVSDiqxwz3kv8tZ3I/AbXlxhhDp008vurPWs3DKTMztcfkzz1eGvgRoN+zfffHMi1lHryD5gEqbBDwDTgPzcIxcwPlJk7vHyzZgDIQlHn8dm2Zw8Nq8q5pSex/YyOVz7v20LHK7Ix1VNZtTkXmxZXcDyOVlIMpx5URptuwVeBKZu8fLAieSnF2ObAtMfeEa7F+1ENXwM++aqsOOefXVbmrZPYvqYdezZmEf9JrFcMrwzXfo1Dq7LpjizCM2pYgt464LX8RSWYQcbinYu28GH177HrZ/eSZs+R6dtunjCXEw9YmdRICBe5ZH4PT62rVjP8KfuRn3seTasWIuEjSbrGHIgyFIO2ZYiCkOXf6xZBu6qghW2gdtSOL3H2Vxxy/0Yhp+HH72EsjBLVBSVnj0GHPG9RsIZrQfxzOWfMnb+s2Rnh+ce9xseFqyfWGHYhRAs+n0is1aPJr80g6SYhlzc/T76d7yhzvM+Adhq/s7HntFsMTfRUmnDnVH301nrVmvj/21DMd3PcvPz7FK8ntAvpqLAKadGplL9o1g8/WBE6TS/x+anydnHbNgH3tiECS9tC1utEpuo0axDzW397U9LpP1piSGf/zZzOyW5HuxD9EoNn8mOFQfYvyGHJp3rhx2zZedkRrwfaqR+m/wLc1+cgb/Yi23baAnReIo9FUa9cg6DqY9P5rHFTx6xYTH9Opnb9lWETsJDompzkSRJqMFyyN2rVhMt/EhSoPTTUKq/FCUhIVkgQvS9BZIlIas2btsI8ZAVYZG/83cANM3Jww++xxtv3Ydp6hiGjqo6UFWNEQ+Mxul0Uxvo2Ph0burzMO/N/CdevSTsOR5/5eefLXiYpVsm4zcDLJ1ZhTv5v0VPsCNrNXec93atrKkO4THFO4GHi4ej48fC4jdjJTN8k3k29jVujbqnVuY45iDzrFmz6NevH99//z2jR4+mX79+7Nixo+YLjzOG3BSP6gg1EKoKjVs66H5W7XyhqsIyBREaBxECjAilg0eDS+9uSUpjN1oVz12SAgnNER+cekze1sYfd+MvCx+KsG3BtqXpRzXe6kkrmPH4V5TmFGP4DCzdojivJMSol6M4u4iC9MjdkuXQPT6+e+o9Rp9+I/lbdnN4At1D6wBt9s5fQHFWDqIKmY4MKLYIqVVXLAWH7Awa70ApoIqJGx13mPAKBN4B6Vs3UJwXaA5q0qQNr786i1tvepLBg67n/L6XcN0FN4CnqMIRMAw/y34cx+inz+PVR3sz+eMHyMnYXuOzqIpWDbphWnr4gwL8pXnsyfiNxevGs3jzlxVGvRx+08PybVNIz9sSfow6HDPy7TweKr4TLx6sYC+IhYUXL8+VPMwBa3+tzHPMHvvgwYMZPHhwbaylVhGXoDB+QVPuviSdkiILkLAsQesOTj6Y2fi4bDd7DEhk3EgJM0wM3B2jcObgY+eriIpVeWdJP756Yzvzxu/H77Xo0DuJm55qHxJbN3wmS8etZdnYNfhLdVr3bcrAh3vToH3oOmzbxhmtIclSRadpVciKhCPqyP9cbMvm+xemY3iPokFJkrBq4HvJXL+NKcNfxFtYElAckkx8WlSYBGewHruq6yIEDqHjOXCAL+98lMYd2pG+aVPFYadl4UfBKr9GkmjQtCm3PP04v8ydzrI532JbZrUO3Eh/RYrm4GDGPuLq1ce2TJZMfYdfZn+GEaxz36DIyE4HMQn1GfrIR0wZN4KczG0YeiA/sOHXWWxe+wM33T+OFu2OTL0oNiqJfp2uYeGGL7HsKrGf4MujJH8vr316EaZbxlDCU/7qppdFv0/kur4njsArsETB8pxZfL37bXJ8+2kWcwrDWj5C56Q+NV/8F8IM32SkCP60jeBr73geiHnimOf524ZiANp0cvLWF0msX+7BGe+iW58Y2nU5fuV3bU6NpVv/RNYsLMDvrTRQmlOiSRs3PQfWTsNGdLzGbSM7cNvIUJbDcph+k/8OmkDm77kY3sCXfPXk31k/fRvDv7ma1n2aUpZbypL3lvHrF79SlluGM9ZJjKRQIgJhKo1AmEFHw7YEXS8MX0oZDnm7czF9oUZdFoE/4HAlga4YJ/WahSZby7H8nUms+OBrdNOquF4VNlGmH48aDK1JEggbGYFTLsUnHFgoKNi4hQ+X8CMEeAoK6T10MAc2b8Fh+FCxsSQZSQiEJIGq0u28c7jq8YcAGDD0JjYu/5nSwsKKRKhARiJ8rsbwe8ndsZGm7Toz68Mn2PrLPEy9kqdGWBaWx0uhfx+fPX8tHqUUw6hy3LYwdC9Txj3Ew68sCeuIeMry2b1jObIs06JNH1yuWG4a8ALL1k7EW2WnIguBy7CQIUB1YBkQqeZBwJIN47n2rGdqrKCpTYzZ/Ahz0sfhswLVPdnevazLW8Rd7V/ismb/OGHriITaKuXMtw/iI0ITG35y7PD0DUeLv61h/+2nfF64fmMg/CGB4RdsH5LCY592qJF//Fgw6utOjHl8J7PHZWDbAssUtGghc8ejKUHP6cQkppZ/vp6szZVGHUBYAt1jMP7mb2jcysmelfurhY78JX7cioyGH6coQ67Cj9L87FOJSzlEQu0wUBxqCIMiBMIdthL6DDS3xiXPDonYD5G9aSerPpqKYYTqgbptA6du4ldUbFmgSiYOdCQsouzwpF+GX8ebmUUzRcfn8SAFKbEEEkUxCaS0acWlD1TWkkfHxfOv0WP5eeqX/LZgLpZp0qBZEzJ3rccMEasWSKbOLxNGs3LSO3hUX0hVjBScUQiBt7QI02mGNbaeskKyD2yhQeNTKkcXgnmzXmLZwg9RFAdIYFkGAy56jDP734nq9xGDwJYCs8iH/BpUE4yIPXoC2/SxfudcurW5ONJJtYqdxeuZvf8T/HZ1g+e3PXyw5TH6N7iKBGfKCVlLVei2nw9yX+LLgvcpsvJJVRtxR/JjXJd0zx828p217kRLMZSK0DxItBRDd63XsS4b+Jtyxezf5uHJIesoyjXwlFh4ii0Mv83SGbm8cdfm4zq3wylz/5utuf5mB41EOk3FXowtO/j0zkU82e1bSvPCG5raxvJP16J7Quu0ZUzMnGz2LN8bNuQiLBvFNpEBGVHxk7loI6s/C0/fGg6JTZKIbxhadikBbskiOdmJooCiKSQ1rcfQt2+g22WnRRxv/ZffYemRwzoyArdtEG37cKIH6WyJGH5XNZVtM6eil5QgE6C/lQlwtte3/dz2xis4XNV3d9Fx8Vx069089cU3PDtxJne99B49zrsMzelCCtIMIASSEMSYHgyfB5+/JGypY+BZSEh28JoIEShZktF91WPhKxaPZfmijzFNP35/CX5fCabhY/53r7J5/XfERNULineHGnUA2Qr8hDwbIZAEmIaPbelLwy/oOGBexgQMO3xuQEZhSfb0E7aWcgghGL7vYj7Ne4MiK5D3yTYP8Gb244zMPLLmsXA413EBSVI9lEPe4hISUVI0l7jCV54dLf6Whn3yG3vDdkLqXpufp+RQkBMhwVRLWDV1D0u/2A5+HYLiD75Sk+ydxXxyR/hmntqGv7TyHmUs4iiiHnnESsUBPc8aYInqf3iGx2DhqO9qqECphCRJXP76MDS3VrFJUYRJfXJIMzJILt5HKzK54Kq2PLLgcbpc1PWw45VkHkTYNsphGDBkVUV1BzslhcBt+knSS0j2FROnl6FWIeeSbG+FElTI2oVg988Lj+geB9/5CLe/8AEtT+mCE5soy0eCWVpNiSlyxWTggKyoSHL42IhtmaQ2rmwYEkKwcO5bGHrodt7Qvfz03asM7H0PDjVccUBwTyKB2wuqQTXBbMUGhyVQZY0oV/heiOOBUqMQO0JIyxIGXit8lc/xxIqy+Wzw/oJfVH/OXuHh28LPOaDv+UPjKpLCt0k/007tSBTRxBBHtBRDC6U1MxMX45Rqp1rvb2nYNy4tjFjDrjlldq0Pz6tRW5jz+gb8YQqXLd1m47wMSg4ef6+97TnNkVUZGYsECivi5QI5Yst7TfCX+CnJKq75xCBa9G7NXdMfoN2AjriiVRpIOUFvWmCU+bB8Otu/ns+ih8dUuy5v3VbWvDKW1S9+RPaK9QghSOvWHsWpBcWmQ6tXVKdG16sHcsqFfdFcThKtYmKEF03YKAhctkmSXooTG9XlpG2/Hhje8LFOw+slb0elcEjhrh1k//YLvoLwFTuN23SkZctWRBuluET18kelhnYJIYOiOlDcoa37msNN7/Nuw+mqZFD0+0rweosijncwZye9OlxOk8S2yJKCLKvIFbHySnIvCXD5Iao0INjhNAWaHdy5SAq92ld6jraw+X3XD/yw/DWWrh1LmbfmyqWjQbd65+BWwlN8qLKDjgmhXPDHG3OLp+CJwDUPsLBk9h8eu5HShAX11jIjaRFvxn/ElMQfWVZvCy3UI89h1YS/ZYw9PtkBeMIesy1BXL3jK2ybty/yi0N1yhQcKCM2+fglcYUQnHVTJ36dtAnNLA7I2R3lGIoUapFs08IRdXQeRWyMysC7erGvQzQbPsrA9FaPR1tenZ3Tl3D6v68nqn4Ci+8ZxYH5K7H8Aa3RbZ9PJ6lzW3qPfpxfP56GhYEiRJCIC5AkFE1l4DN30fmKQB19UryTLZ9PoKoTWH7/CaYXzfSS88NcJIcW1plWnE5iGjSgYMdWFj/xAGXZmciqiqX7aTbgQs749wsozurPISGtKarTjXkIu6MkJBy2guWUqoVkBAJJklGdDq59aAyO6GimjHuQovxMZFnBti3OGnQH515SXdJNc7iRDvPbdMoKY54/DQEkCz9mlAvNEUVK69PYlrEc3Qh8L2RJQZJkJIfApnJ351CjGNjjHuonBigmCksy+O/E8ykpy0E3PKiqi2k/PcrQC96jZ8ehEddxNOibejljtz6FbvuwqshYaZKDFjEdOSXh9FqZx7RN8vUsotU4otXD95NYtoVsVb6YbYVAtZQU+N0diUhJTeisdavVpqSq+EsadiEEW5blseDzvXiKDU67MI0zr2mM0x3Yzg65tzHbfivGVxYaNkhI0Wh96vElAKvfOo7i3PBeuanbJDU5fvP/PnElS56egTevlGRLxScpVNWIkIWNhYwklQciDk3oCiRspBAJPIE7zokr4cjq/z1ZeSy67SXyN+xEdqiUeEwsEd4gyZpK1srN6NmZAaNexfibHh8H125h60dfc+VnI5l+94tYfh0hBMIWJLdtxpCPniaqXiVDYt7q34go+GoaUObHIQk8WnjCLglodnpPvr/9KozSQBigfLR987/H8vno98p/q13Trt9gFo17NeyUThz0v/0pls8aS0HWPhRVIy4ljXa9z6fHgOuITQw0fT3wwk8czN6N31dKSoPWOMI0LymKRsdTB7Nx7QzsQ7alqqwiWX7MYKmjBmglPiTZwJmRy22XfciCXz+iqDSbFg17MPCMeyny5jJ31WiyCnaQEt+cQT3vo1Pzykazj6ZeRX7RXmwReAJGsPZ90vf/pHH9LqSlRK7MOlI4FCf/7f0zo9beyLai1WiyE8P20SN5EI92GXvM1ShCCL7c9zoT9r2CYevYwqR74rk83O5D6rtCGWVN22BX8RpceuW3Q1gBlSq/MxCC6xtzwTGt6XjjL2fYbVvw31tWsfKbDHSPhRDw25wsJj6zif8sP4ekhm76XVmfBV9ls+qHvArjrjklNIfMM5M6H/eW6cGPduH96xeGJC9Vp0y3S5oSk1T7Xa8AG8YtYcEjUzE9AQ9Mw8SQnAipMuKmEDDsiECstUpPJopDwYUnmKSUqDD4wbBHbNSReSm2aTF38COUpecgLBvLp4OkQoQ4MoAa5WTtx1OqGfWK8fx+dkyYSef7r2P40s/Yt3w9nrxCkts2o36HliHnG2U1KQNJyAhi/AalzuDuTZJQHA4kSWLA88/x+2fvY/pCQzWW30/60gWUZh4gJq1RxeeumDguffI9Zrx4L0IITL8P1elGkiSGPPMBTbr0ouu5NSfGklNb1HjOxVeMYv+eXyktOYihBwyt5ojCJYHpC3UohG2Rl7WVhnEteeD6b6sdS6M97ZuGp+U9kLOBnPztFUa9KixLZ+Hq9xh6/rukZ6/B482nYf0uxEaH70yuCcmuRrx9xnyyvXs56MsgLaoFSc4Gf2isQ/Hhzn/zzYH38NmVu/hV+fMYvroXX5y+iVitej5hfNYb7PUGiizKLYVEoFTXbSgMTLmaZs7aC5scD/zlDPviiftZ+U0G/rLKPzZfqYXu9fL2jasY+VM/ZFni2cmdWT7zIDM+TKc4z6D7uUkM+WdjUhodfxrZboObcuGDnZjz+gaEEFiGjcOt0qhjIrd9cNZxmdM2LRY9Nb3CqJdDwcKsYqQlwCFMTElGyApalINGpzXmnEfOJaFhNBP6jMQwDHQ0DBEweiomDgySWjU9orWkz12J72ARwqr0+lUR0EoN5yFbPp3GfbuwJPfQ2K3AiR8HOvhh3umDSe7Tk26vP42rfuQtbMOzzqB4735sI3yipXw34rBtEnx+dLeLhE6daNT7DBq0asGGkY9S7C2NuN2WVY2839dXM+wAzbqdxR2fLmTz/OkUHNhDUpOWnHLOZbhiapdsLjqmHvc9tpB1q6ex/rdvkGWVU3tezaofR5OTURh+zYpGYd4eUtKOnLkxJ39HRPk7W1jsObCSF8a0xeMtQJZlTNPPqadczTUXvIeq/jHnJdXdjFR3sz90bTiUGAVMPfAO+iFlrzYWZWYxszPHMbTpg9WOTcp9B58IfalLgGJLPJv2fq2t73jhL2fYZ7y1rZpRL4dtwdZleRRk+Uhs4EKWJc66LIWzLjvx9a8AVzzbnf63tmXVtD34PSbt+zWg7Vmpx223kPd7JrYR+lwcwsRErRZtkQC3S+GKz2+g3cUdq51fr31Dctbtw2kbOKmMCWtRTk67d+ARrSVn5SbMsupfDAUbGRu7qnEP7gSihBdhmrhTk/FkVGp0uvGiUtnpaesGuYtW8vPgmzln7gQcieErNzreNJTtU6aHMeyBcs+qvwJZQAwyHfr3Ja5BfVb9+19Yfh+oMigRagsk0KLDh9PcsQl0v+zm8NfVIhzOaHqeeSM9z7yx4rPdG78nN3MzIkzVk2UZJNRrflRzJMU3jSiuLSGRl7cF267+jNdtmQoIrhs89qjmOl7YULQUTXKgE7qT8dseFuVODTHsRWaobGE5FBR8dhlu5ch7Ov4M/OWqYgoyIleUqE6ZwqwTUyd+JKjXNIYLHujEZU+cSrs+DY5rCEjWlLB6nhLgFj7cUQqqO0AZ0LBHE4ZNvT3EqANc/MXdRNWPQ4sOeFySIqG6NTrd2pfmg46MdtiZEIusVfcZJMBlm2jCRBY2krBRMYkVZUS5ZAo37aLjPdeiuIPzYlcz6uWfaaYHkbGPBV37sWLwdRSuXhcyf3RaKhd8PoaE1i1QXE60mGhkTUVVJdRDxC5kYeEuy2XP2/9h3WMjUEs8qIaNaoXTOA2uQ5JJPa12GklqE73OuQcljKcsyQrJDdoelbcO0LRBdxJiG4ZN1ioREriG6WXt5ikUl9ZOB+WxQpXCJ8jLocmhO/g0R+QdgyY5iVVDCfRONvzlPPYmneIozA6vNG/qNvWbn9xv0uOFpPYNcMZHYZSF1uhrTpWe9/flzKdr5vSJb5bMbeteYuu0X9m3cDPuejF0uK439btUhmF8mQfZ9c5XZM9ZgqypNBp2Ps3uuBwtLlCW1/zKs9nw1iQ4pC8nEAay0ISnmtdsWxZabBTtbrqU/PXb2DNjAbKv7JBr7WA3aflFgqK1G1g19A56fjWWhO5dqp1fr0M7hsz4kuK9+9FLSolp1IAfb7yFsgMHsMvFN4QgRniQEdhBiuJAYw+4dRtTkYOVEMFZJQnF6eSs519D0U4+ecVGzU9jwGXP8dO3zyAQWKaOwxmDOzqRa+6ccNTjSZLEXVdOYfSEgfiNMnSjDEV2IEkyca54ikvD6zSoqpMD2euIi6mdGPmxoGtCf0RIIUAALjmaC9NCd1e3pT3Bf/bdh++QckenFMWw1PtRpZPfbB6z5umx4mg1TzfMz+GlS5fh91QPOzhcMmcNbcJ943ocr6We9Ng9dxMzh32MWYV4S9YUolNjuWHlE7iTwivKHw3Kduxn6bl3Y3p8iGAnqOxy4kqrR5+FH6MlBmiD1785iY1vfxWSDJUxUA8ppYxqkkr/sU9StHw1SpQbUlJYOfwJJLsybKLhr970UwUJPbvR65svaly7XlTMb6+9zt45gUYrtypwGqUII7QzVAC2ZCJU8GgObElCE3DK0Jvo+PAzNc71Z6K4MINNq6fiLc2jUfOetOl0PrLyx42RYfpYs2UaezJ+IS46lZ6druOrOcPZsXdh2POdWgzDh86meaOTY1czJ/NT3t52P/4qyVOH5KJ5dAfeO20JDrn6LkcIwTvp/2ZSzjvIyBVkE+ckXs5zLT6rM+xHgj8iZj3nvZ188ch6ZFWivJmwY/9kHp3au6Lk8X8VGct3seS5GWSu3I3q0mg/tCdnPnUx7uTaKbFcMfhf5C1eGxKmkBwaze+6nA4v/bPis6wl6/j9/WmU7M4krlVDvPsP4NmbjuUJGHtJVZA1hWat4vDt2IUwLSRVQRgmhixTWtFMFUiiRgxkyTIDd/6KrB1Zf4Jtmlg+H+sevp/cRQvCniMQaJgoVYrhBSC53XR44U0aXDTkiOb6u2Ldlml8OftOdCO0AikuugHP/HMXsnTyRHpX5c9j3O7n2Fm6nmg1jksa3sl1TR/FdZhYeY5+gMVFs7CExRlxA2nqanMCV3xs+EsadoDSAp1VMzLxlZl06JdMs07xNV9Uh2OCUVTKvJaXIiJUm2hJcQzaMyvi9ZZusHvSXHZ+PhOzpIz6fbrhOriPgvmLEf7qISQB+BQVvxJI/NZo2HetRlbDe1LCsvDu2IGkqrhatqzIdax7/EEyZn4bIZYu0DBQqmzjy6lnXM1bctacE8elcjLCFjafTr2G7XsXVBh3WdZQFY07r55OqwgllHU4MTj59xQREJPo4Jyba68sqg41w/L6kBQZEYGLK1wNelUoDo3WNw2m9U2BWL9ZVMyyjmeHGHUIxuMtEwkbp9AD8W4pfLlkUu8eEY169oQJ7B01EuEPrE1NSqLVG2+S0P9smlx1LdnzvseKQC0gHxKbLZ/Zu3cXwraRIjBR/i9AlmRuveIr1m2dxtLfPqLMe5BWTfpy9un/Ijnx2HRc63Ds+Msa9jqceDjrJ6HFxeD3hecKSTjt6KoufPszkB1agD7gUAiBA5NYM3DMFDLFmjvgXVcx7kqUm3bPPhp2/MyxY9n7/LPVPjNycth66y10mDyFxNNOp8EFF5P1/RwsbzD+Kstg26iESt5VzOl0HbNR92UfoGDVQiRZIenMgTgSIvPQn6yQZYVup1xNt1Ou/rOXUodDcMyG/fnnn2f58uU4HA6ioqJ48skn6dy5c22srQ61CH96NhkfTqFo8VocqUmk3XkFCQNOP6oSTEmWafPU7Wx+/B0sT/WyUtntpN3Td1b7TFgWxfMW4dmwGbVeIolDLkRNqqw9d6SmYIej4hWCKLzVOG402yZe91KmOjAVBRSFen170+6JEcR2aBdmCMG+F18Iex/CNNn99JN0+W4unUe9Sv2zB7Bn/Kf4c7KJ69AJz8bV6Ol7Ij6HlIF/nKdc2Dbb/jOCrFkTkYJhJts0aXbLQ7S48/E/PG7YuYSgYNsqDiz7FmGZNDj9IlI6968Tq/4fwDHH2BcsWECfPn3QNI0FCxbw4osv8uOPPx7x9X80xl6HI0fxig1sGHwfQjcrK1miXKRcez5t3v/3UX/Rd7//NVtHjQ2KFQmUKBdd3nmE1IsqZcz09Ey2XXgdZl4BdpkH2eVECEGz918h6arKssv11w6nYNFKMCvj9k7hD6lhh8oYt60qpNx4LS1fe/Ew97yC3685TAu/LHPGnn1hD+0f/wm73nwR2x+mJ0LT6LN4A1rcH8vp7P30DfaMew37EI512RVF+6feIfX82vF+bctkxSvDyFnzE1aQ4ldxRpHQqhtnPDaB3Qu+YPeC8QjLoGHPwbS79AGi6jWqYdQ6/FVQq8nTgoIC+vbty/r16yMq4RyKOsN+fCFsm+VpA7GKQhkn5WgXHSb9h8SBR6apWRWWX6d4/Q5kTSGuS5tqoQkhBL+fcTH61u1g2QG/O/jykNwuTln8La62gTisnpvHmotuwMjNwyrzgATRdlnEMIgAbFki8eLzafPpmAhnQd6smWy/5zCSaocx7Jbfx5obh1C2Yyt2Fe4VyeGg68eTSOzZO/K4h4GwLJYMaoFZHL7tP6pZW3pN+fUPjX0otk59g81fvojlr/4CkTQnRDmxTT9WUKpPVh0oTjcD/7OE2IZta2X+Ovy5qNXsz4QJEzj77LOP2KjX4fhj9yNvhjXqAHaZjwOjJ/6hcRWng8SeHYg/tV1IvLlg4jTszVtQLRMVGw0LRVggBMIwyPlgfMW5jpR69Fw2g7Zvj6TB9VfQ8JZrkJTDl6zKLhexZx6+Rjq6Y8eg/mkYv0UINBl8m3+PcG8uuo3/llYPPU1021NwNWpC2lXXc/r0hX/YqAMYJYVYYYjFyuFN3xXxmOUto2jdMkq3reNIfLEd098JMeoApuXD8BRVGHUA29QxPMWsGnNPjeP+FZCRu56l68bw6+//h9cfmbv+74waY+yXX345GRnhO8yWLVuGEvwSzp49m5kzZzJhwtF3uNXh+MD2+sj+ZBqh1LyVKF665oiFer2r11M0dTboBjEXnEP0uWeFGHXjQCaZ/3oyIIZRBXLA18YyLbybtlY/pmnUH3IB9YdcgBCC9TNnYuYejLgOyeUiZeiVh12rq0VLojt0oGzTpuoJVyGQEDh8ZewecjGtf16G1iAt5HrF6aLx9bfR+PrbDjvP0UCNigmKnUQ4HhfKfSOEYO9HL5A+cTSSoiJsCzUmnsDWNjMAACAASURBVHbPjSWxx9kR5/IXhe/OFpF8LiHI27oSvawQR/SJU0+qTfj1UsZOH8KezOVAgHN+8o/DufLcd+nV6dY/eXUnFjW61t988w0rV64M+1Nu1OfNm8dbb73F2LFjSU5OPu6LrsORoWTF+gCHTEQIJEOn7JcNhx1H2DbptzzAnvOvI/+/48h//3PSr7uX3f2uwC6t3qCS98HnCDuUjCygKRowsK7WkalpJUkibcS9yFGhXOQCUFJT6DBrMkpsLJ7Fi0m/bAg7W7dhz+m9KPzkE0SVWH378RNwNW6MKiykID+NJkyiLD8yIPw+8sd9fNh7r03IDicp512OpIbSEchOF42uvCPk833jXiF94mhsnwerrBjbW4aem8Gmh66kdHvk35srKfRlVRMkWcH01UR5fPLiy7m3sTtjKYbpxTC9+I1SDNPL1Pn3sSdjxZ+9vBOKY46ZLFiwgJdffpmxY8fWxchPNkgSsiwFtT3D+4mKCt5NOw47TMHHEyiZ9SPC44Wg5qldWoZ/01YyH3q+2rlli1dAhAYmAMm2cDc8/Ms/5c6bSbnleiSnEzkmGjkqCknTSPn/9u40Pooqa+Dw/1b1nj2BQBAJKouAIigOoOiwzQBKWBRl1XEEdUDEDWfgVWdUVBZFQcQVBUFARQEFlGUEZRREcGBkM6xCJIBEA9l6rar3Q4eQ0N0hkCbdCff5wC/pTt86Scjpqlv3nvPXIbTatgF7k0aceP99sgcOwvnNN+i5uXj37iXnqafJHjiopC+ruVYtLnvuORw2M7Gam1jNja24UTeA4fGQv7riN/rDofGjE7FdlI5qP1XeQbXHEHd5KxrcVbbKoO52kTVncsCN1pPPHZw5MeRxmtz6KKo1cFelKGc3qMkehz0p8A1B87jwBYnhfCrK/5U9Wz/lQOa/0Xxn7lGcX/Qr2/ctwacF7qXw+px8uTH0z6omqvRyx7Fjx2I2mxk1alTJY7NmzSIpKforoEUj3e3B8PpQYytfzCyuXUt/A+ji5hqnp3YTPlSzHXNa+Yn2t6nv+JP6aQy3h7wFS0mb8gyK3V8lT62VXO5YZsNH3pTpJNzSE2uz4DfqhBDUf+Zx6jxwH/lr14GiEN+xA6biMr16QQHHxozFOG1jkeF04tywgcJVqzClJHN80kRcGzdgLSzAB/hUEbDBSYmNKzfecDPHJ/GHees4tvpTfv1yMcJkpm73/qR06BZwb8GZtSd0IjZ0Tmz+JuRxLrvpPo7v+x9ZX3+AoWlg6CgmC8mNWvHbgc0lK2VOUq0Orhz4VJmptRP7tvDjWw+Tu3MdAHENr+TKYZOpdeUfy/0eNZ8Hd14OltgkTJaKddw6Sdd8fPnJSHZsnINqshTPIgq6DXiHxlf1Dfm6Y7m7Mam2oIkdDLKPBVYBrckqndi/+y66LnEMXcdz7ASqw4oprvpUenTvPsAvD44nf/UGACwNL+KiFx4lIaNTyNd49hwg952P8B44hO3qFiT+tR+mlFNvqIrNSoPnHuTAP15CuDwne9QX/2ugoiPMJhL+VH6zYF+IapoAwtA4MXk6prq1sd/cjZRhQyj8ZgNGYeAZnijedGR4DI6/OZM6U8aXe1xz7Vok39or4PGir76GEDdYjcJCcidPRtu5tSTxK/jbxKmGgdtUdveqkbkT9+YfsLa+ptxYwkmxWKnT/XbqdL+93K9THXEYWohtvoDqCF3/RygKbUa9QdNbHiH7+6WgadRp043ES1qSvWkZG18b7p92Kf5ZXDHgX1z256Elr/9953d8+8Sf0Eu9AeTt28J3T/Wk7b8+o3bLwP+Xms/Dlrlj2Lt6Boah+5e3XjeANndNwWyv2Bvo2qVj2LlpLprPheY7dYP3i7l3Epv4JWnpwfufJsTWQ9NDn9knxF5YSzmrba2YYA6+u5LMJ+bgPV6AoRsk39CClq/fT0yjemGK9vzwHDxM5tX9/KtXSv06hMNGg7efJql/j4DX/P7WfI6OHo/h08DrRdhtCFWlwbJ3cFx3ddmvXbyafQ88j+9YbvEjOqrNjGIy0WTJa8S1b1VufHta/QnPrtNWbBgGFtz+NwerFaEoGLpO/GMjKdybRd5ny0+VCij+nsz4Sq4b7DdeR/3PF5zFT+mUvI8/5tdHR2MUBF/tI0wqapAzNwPwqAJN9bcGVAwDq66hxMZSf90PmILcRI20TYPaULR3e8DjitVO+j1PcvEdD5/TuIauc/zAVnSfh8SGLVHN/gqHusfNjumj2P/lTHQRvNxtXPoVdH51S8DjX03qxdFtq8tcDSgmKwkXN6f7c9+fcbeu113I60/WwecNtnJIcGmLm+kz7NOQr395Xjuyjm4KaDRiMccwuPt7tGx8S7nHr0lqzLrE/dM+Y/uoN3EfyUV3eTE8Pn5b8yPftH0E1+HgW+CjxdGJM9AKigKW5hlFLg49Mqlkzvgk9+6fOTr6eQynC4pLzhpOF3pBIVl97wsoQ5vcpzNtsv5Ny43zqTviNlIyOnLR2HtpuWPpGZM6QK3HhiNOu5l5skCWAHC7/WfHbjf5L71GUq8/U+eJhxEKCMPfOcmC99RkkNmEpcXZlR8ozd6uXcn3HcBqRajBV/gIwKQbCMPArOtYdc1//eL1kleFN1HPxuX/moHqiEWYTlWuVGwO7OlNqHfbfec8rlAUki65ipTG15YkdYAfX7iL7NXz0EPUMAco+CUTb0HZtfi5B/4XkNQBdJ+b/MO7Ofy/FWeM6XjOHhQlVIVOgyMHy1/j/5ebPyDWnorF7L9/IYSCxRTD1U0HcmWj0NM4NVGNqBWjub1kPjGnpBxsCd1AK3Sx/+XFNJsUvmVr4Xbi0zXgC1xJAqAXFOHedQDb5adWkuS+/YH/TD0Iw+ujYPla4jK6BDznuLIJDV86+23rCYP64tqyndx35vvXovt8mHyBO0MBjCIneZNeIe3b5RTMnod3z/6SG64nCZOJxPvOffmZuX59Ynv3pmDJkrLz7KqKYrej6G7wBO+kpRoGdu20m7tuN65v1p5zPOdTbNNWXDNvE1nvTyH3u1WoNgd1e/+Vur3uQrWd3fz1SVphPrrbiSmpdpllrs6jB/h13RJ0rwvK7SNi+GvqlHLkx3+jn/5zLeZzFXDov8uo1zrwyrM0myMZTQs9nWJzlH/fLjmhIY8P3c3mn+bz04FVOKxJ/KHFX2hY79z3HlRXNSKx523ZF7TqH4Du8XF40fqoTuzlXaIahoE4rfemd39WyJUnhk/De+hoeOMTgrovPEnyiL+Qv2QV2pGjON9423/FEIR3zz6EEFy0aC6/dLsF7fgJjMIihN0GhkGdN6dgaXRppWKq88pU8HnJX7SoZMGPomkkXN+egi9XltMOLUSD6ihepmtLS6fxYy9XehxX1h72Tbif/M1rQVEwxadw8fBxpGb4uwjl7ljvvzLwulB00BWCbn+Ib9gSs6Nsg26hmv3TccHON4RAqUDHqbiki6mV1oKjv/w34OrVZHbQqsP9ZxzDao6h3ZXDaHdl4NLRC0mNSOyKWcXQQ/8pl7+WO/ISb+/GsekfBJ1eMKUkYmnUoMxjtqtbULB8LYYrcB5ZqArW5o3OS5yWSxqQMmooel4+B197K+TXqXVSATCnX0zDbespXLEa94/bUWunEHdLBmpy5VdMaYcPw7IlxHjdJTeEFcC7cgWWyy7Fc/BA2bP5cjZgCUcM8XcNDfl8TeA5ls3Wu9qjFZyA4jlob042+yeNpGjvVgyfl/xDuzCKz7pVrTixn7a3TbU6aPm3qQHj12/Tiy3zxwY9tmqxk96+/BvFJ/UYMocPpl6P1+tE8/pPHMyWGOo2+ANXtr+wk/XZqBFz7PFXXYJqD35GoNgs1P9L4LRENEl97G5MSXEBKz2E3cbF058M2BWadPftwVeFKAqmurVx3HDt+QwXJT4Ox01/giAdi4TDQfyoU3O/wmQi9uY/kzL2YRKH3RmWpA5QOP1VDI+nZONTyX9kpxN27Sb+7mGImFhETCxYrVivaUOtadMRDgeUmq8WjhgcPW7G3qlrWOKKVtnzXkZ3FZYk9ZMMt5Mj86ZwdMFrFH67Ct1dvJLIALMXFJ2Si5yEhi25fvxqki8PnNqITW1I4673oVrLtl9UrQ7qXdWdlEYVa5OXnNqUu8bupG3X/6Nug7akN+nKnwfM4Na/LUdVK9YhS6pBq2KOLF7Pfwe/gO48NUcnLCZsF6Vw4+ZpmOOje+mj59BRDj/+Csc/WYXh8eJo04J64x8i9sbgPVwL135PVt+/+ee8PV6E2YyamkL6yvewpJ//pV3a77kc6dwL36HDGAX+ZXPCbsferTO133/znOqVa7t34/70Mwy3G0unjpjatw9Z6uBY+7b4dgSv9SLi40l6bw7m667Dd/AASkIipjp1APDu28uJN17FtXEDptqpxP31Hhzdb6p0KVvXmpXkvzIJ397dqA3SiRvxCLYevcuMqx07iueHDQibHWv7GxFWazkjhtfmvk1wHwpdi+YkjxmKHJQ9S7fFkNy6C1f/6+MzThvu/88cti8aT2HOQewJdbj85odp3G0EihLdV801TY1J7AC/rd1G5hOzOb5pN6rdSv07O9P4yYFYkqt2E0pV0d0eCpatwZd9FGuLxjg6tqvSWtuG10vRkhUUfbEKxeEgpn9frO2vPesYDMOg8LG/43p3Jmiav4Svw4Gp5ZUkfPYpIiawCXdOj254160LOp6IjSX5syVYrqmaxub5UyaSP3U8RtGptfvCEYPjzntIfHoShqaRO/YBCj+ajbCcTOYGSS++SUyvqmlSseW2FrgOZJ75CwGfCm67glInDWtSKul9HqBel8FnLM4mRY8aldil6sn1wYcUjHwAik7b1GS1Yu1/O3FvBJbndX68gBMPPIBRFFjbRKlbl9SdmVXSuk7L/oUj7ZpDsNrtNjupK9dTMH8mBe9OB99p91CsNixt2qL9vA8lKYXYu0cQ0+8ORIg2f2fL8Hgo+uEbDLeT3O3ryf5gCkawOINQbA6uWrAda52LwxKLVLVqxM1TqXpzvjg5MKkDuN24P1pA7OQXA87abX1voWjePDzr1516rcmEsFhIfGdmlfUjdS5dFPpJn5fCD+dQMOOVgCWfALhdeL79GgRo2Vkcf+IhnMsWUmvW4kqfHeetWMiRx4eVrC4xPG4cMRaKTGaM099gghKYk1IrFYMUOTXi5qlUvWkHgze8AEBV0Y8GLt8UqkryRwtIeHkK5muvRW3UCPvgIdT65lusHToEGej8MAoLwBti7bXPh3fLpuBJPdhYziLc3/0H58ollYrJuXk9h/9+J3r+CfSCPPSCPAyPG4vTR3KDlqhxiQiLFVuDJgiLLeD1wmontc8wFEvV3QOQwksm9iim//Y7RROmcLxzL/L6DsGzdEXALtRI836yiPzW13IiPpm8ho1xTXwRw3PmanylqfXLudmr+VBSg585CpMJx4CB1Pr3alJ/2EziK9MwXXZ+lnqGYml/AyLERiERE4vu84YuwA4BzxlFhRTOnVGpmHJeexYjSEMPw+VEz/yJa5Zl0fabQq5asJ20gaMQFhvCakeYrShWOwltu9LggfLr+EjRTU7FRCnfT7s40aEHhssFxRuBPF9+jbnLH4n/+L2ouJHlmjAJ94QXSqZCjOxs3M9PwPflamKWL63wdIj9kUcoeOjh4HPsfW9BxIYudhVplrbXY2rSDO+OH6H0G5rJjFo3DdPFDfBs/DZ4r5MQCV/Pz6tUTK6tG0M/KQTeX/ZjbdQcIQQN7n+eOv1GkLv2Mwyvh4S2XXFcdkWlji9Fnjxjj1L5g+7BOH6iJKkDUFiE98uvcc89t+JZ4aTn5OB+bkJgMnY60TZuwrdiZYXHsg4ehHVAf7DbT63Pj4nBdEULYqa8FMaow08IQa0FX2D7003+OjVx8f43pBs6Ueuzr3Bk9EOxFS+1LZ3IT358erK32bB17FapmJS4cjogeb2o8WX3Eljr1KfubSNIG/SQTOo1hDxjj0La3v1ou/YE79dZWITz1bex3TkgrMfUjx7F9+HH6MdyUK9pjannTeWuzvAtX+nfoOQOUv+6sBDv/A8x9+heoWMLIYh7dRr2EcNxL1wELhfmLp0xd+yIvn49rpnvYhw7hvrHjpjvugsRZbX+lfgEUt79CC3nGNovB1HrXYSaWhcAW9ebMDdqguenHRil5uKFEBhCQOliW0IgLDZi77y3UvEkDhpOzpQnA6djFAVr81aYUqOviqUUXjKxRyE95zeE2RyyFotRTj/Qc+GZORv3A8XlX10uvHFxiMQEHF+tQklvEPxFPm/ALsYyMQZL+Gdgat4cU/Pm/tcbBu6HHsI3Z7Z/N6lhoH29Fs+E8VjfnYmpQwdEQsJZH+N8UmvVRq1Vu8xjwmSi1uI1nHjq7xQteB80DRETQ+w9o/D+mk3RgjkIiwW8XkyXNiZl+hzUlNohjlAxSYNHULByEa6dWzCc/uWgwmpD2B3UmzirUmNL1YNcxx6F9Nzj/F6/BbiCJHYhsPTtSfyCWWE5lvbjVoo6dILTOyQpCkqzy4n5X/D5Wv3AQfKvaBX8jD02Bvv0aVgG9j/nuHwrluMaNAgKT65TL1WnXgiwWFB798b66muI+Pjgg0QZw+fDKCxAxMWX3H/Q807g3ZuJmpSCqeFl4TuW10veFx9xYsE76M5CYjtlkDjwb5iSo7fYmRQ+8ow9CilJiVgH9cM9/xP/2WppNhuOMQ+F7Vjeaa+BO8gqFl1H//lntM1bUFsH1mxX0htg7n8b3o8Xlp1nN5tR0tIw39KncnFNnx6Y1A1/uS8MA1wutIULce3eje3b9VW2br0yhMmESCg7/63EJ2BtHbwrUKWOZTaT0GswCb0Gh31sKfpV+q/h9ddfJyMjgz59+tC7d28+//zzcMR1wYt9dRKWXt3BaoW4WIiPQ8TFETv7dUzXnLk5RkVpO3/yb+MPRlXR94auL2J/8zWsD42CuDj/jU+rFVOf3sSsXV3pOij6oUNlHziZ1Cl1v9HrRd+2DX316kodS5JqmkqfsQ8ZMoThw4cDcPToUXr06MH1119PQpTNf1Y3wmolft4MtIO/4PtuEyIuFnPnG8JeOEq5vCn6xh+CJ3dNR7nkksDHT8aoqtieehLr42MwcnIQCQn+6olhoF7TBl9mpr9uTHFSL72A5OTHhseDd9FC1K41uzqjJJ2NSp+xx8WdKrBVVFSEEAI9yjbRVGdqg/pYb++DpUfX81IN0DJyOFiClDxWFJQGF6NcfearA1E8/RKupA5gefhhqEBzBgBj756wHVeSaoKwTEzOnz+f7t2707dvX8aNG0dSlC1Hk0JTW12F9aVJYLP5p30AYmMRdeti/+yTKq0WWZrSrBm2ObP9sQSp+36SAJT6slCVJJV2xlUxffv2JTs7O+hz69atQy21AzIzM5PRo0cze/bsCid3uSomOujZh/F98BH6r8dQ21yNqXdGuQm1qhguF95JE/FNGB98ushsxjJ3Pqbevas+OEmKUmFf7jh06FBuv/12unWr2O45mdilMzF0HdfVrTB27Sqb3BUFkd4Q27btYSt1K0k1QaWnYvbsOTW/mZWVxc6dO2nUqGoLMUk1m1AUbKu+RLnhBv+UUUIC2O0o7dtjW/OVTOqSdJpK/0VMmzaNPXv2YDKZUFWVJ554gssuC99GC0kCELVrY1uxCv3AAYz9+xHp6eWu2JGkC1mlE/vUqYEdyyXpfFHS0yE9PdJhSFJUk9ewUrVl7N2LNuNtjL17EVddhXr3UESaLHAlSTKxS9WSNvs9tJH3+zcweb0YX3yOPmkipoWLULrIzUrShS36C2xI0mmMrCy0++/319HxFvfvdLmgqAhfv1sxgvVPlaQLiEzsUrWjvTer3JLB+qKFVReMJEUhmdil6ufnn4OXCwb/mfuRI1UajiRFG5nYpeqndWsIVZfGbkc0a1a18UhSlJGJXap21CF3QLBNSUJAQgKiW8Va8klSTSVXxUjVjkhIwLR8Jb6Mm8Hj8a+MUVWoVQvz8pWIUvWLahLD50NftgRtwzpEcgpq/0EoF4doXRhhhtOJa8Ec3B/OBq8HS6/bsN1xD0pCOY22pbCRrfGkasvwejFWLMfI+gXRtAmiY6dq0UnpXBiHDuHu2gHjt9+gIN9fallRMP3rWcwPPhrp8MrQC/I50eM6tIP7TnXXstkRCYkkrtqImnZRZAO8AMgzdqnaEmYzomdGpMOoEu6Bt2D8knWqCJrH387QN+6fKG3bo7a7LoLRleV8cRza/t1lb3C7nBheDwWPDSfh/c8iF9wFomae3khShOkLP8FzdUs88XY8Devje2EihidIb9mKjLUrE2P71uBli51OfFMnVzLa8HLNfzf4qiVNw7t6hdxnUAVkYpekYobPh754Edrdf0Ub/jf0NWs4l5lK38Tn8d19J2zb6l9+mX0I/dmn8WXchHEO3cWMn/eH7iZlGBi7d531mOeTUVAQ+klFYBSW87wUFnIqRpIA4/hxtI43+tfIFxSAEBjz5iL+2BFl4aIKlwY2cnLQnxvnT+ilOZ0YGzdgrFyB6N7jrGITl14G3hBn+0IgmkbX8k612RVo//sh6HMiNh6RUquKI7rwyDN2SQL0hx6EXbv8SR3AMKCwEOOrNeivVLyCqf7F58GXYgIUFKDNe/+sY1MaNUa0bBV8XLsd84OPnPWY51PMmHFgD7LPwO7A8eiTNfYGdzSRP2Hpgmc4nRgLPiq5IVlGURHGtFcqPpjP639TCCXUjtkzsM77BJHeEGKLm8fbbGC1YXpmAsof2p3TmOeLpWsPYsZPg5g4RFw8Ii4ebHbsI/+ObdjISId3QZBTMZKUmwvlnUUeO1bhoZROXdCC3eQEiI1F6d33LIPzE2lpWLf8hL5yOdrG71CSklFv7Y+oV++cxjvf7IPvxnbrILwbvgGvF1Pb61Hi4iMd1gVDJnZJql3bv8EplLPo1CQaNkS5rT/6JwtOreEG/7rzumkot/Y75zCFqqL2uBm1x83nPEZVEjYblj/KEsqRELapmA0bNtCsWTPef//s5xAlKZKE2YwYPgLs9sAnHQ6U/3v8rMZT33oH5cFHIC7OP6bViujdF/N/1iOs1jBFLUmhheWMvaCggBdffJEbb7wxHMNJUpVTnhmHvn8/xtIlxQ8ooGmIUQ8iBgw8q7GEqmJ6ehzGE//0T+MkJiJCFS2TpPMgLIl9woQJDB06lK+++iocw0lSlRNmM+oHH2JkZmL8exUU72qtzBy2MJshSufApZqt0on966+/Jj8/n+7du8vELlV7omlTRNOmkQ5DkirljIm9b9++ZGdnB31u+fLlTJ48mZkzZ4Y9MEmSJOncnDGxL1q0KORzmzZt4tixY9x2220A5ObmsmbNGo4fP87IkXK9qiRJUiRUaiqmTZs2rF+/vuTzMWPGcMUVVzBkyJBKByZJkiSdm4ivYz+5meOI7FMpSZJ01urWrYvptHITEW+0sWnTJgYPHhzJECRJkqqtYE2KIp7YXS4X27Zto3bt2qg1tKWZJEnS+RKVZ+ySJElSeMnqjpIkSTWMTOySJEk1jEzskiRJNYxM7JIkSTWMTOySJEk1jEzskiRJNYxM7JIkSTWMTOySJEk1TNQn9mhsuff666+TkZFBnz596N27N59//nmkQyrj6aefpnv37vTq1YsBAwawdevWSIdUxqeffkpGRgbNmzePmt/r/v376d+/P926daN///78/PPPkQ6pxMSJE+ncuTNNmzZl165dkQ6njNzcXO655x66detGRkYGI0eO5Pfff490WGWMGDGCXr160adPHwYNGsTOnTsjHVJQr776avh+x0YUy8/PN/r162fce++9xpw5cyIdTom8vLySj48cOWK0bt3aOH78eAQjKmv16tWGx+Mp+bhLly4RjqiszMxMY/fu3cZjjz0WNb/XO+64w1i8eLFhGIaxePFi44477ohwRKds3LjRyM7ONjp16mRkZmZGOpwycnNzje+++67k8wkTJhhjx46NYESBSv+9rlq1yujTp08Eowlu27ZtxtChQ8P2O47qM/aTLfeSkpIiHUoZcXFxJR8XFRUhhEDX9QhGVFanTp0wm80AtGrViiNHjkRVfE2aNKFRo0YoSnT89/vtt9/YsWMHPXv2BKBnz57s2LEjas4827RpQ1paWqTDCCoxMZG2bduWfN6qVauQjXkipfTfa0FBAUKICEYTyOPx8Mwzz/DUU0+FbcyIl+0NJdpb7s2fP5/33nuPI0eO8Pzzz0fdm89Jc+fOpWPHjlGTRKPR4cOHqVOnTkkROlVVSU1N5fDhwyQnJ0c4uupD13Xmz59P586dIx1KgMcff5xvv/0WwzCYMWNGpMMpY+rUqfTq1SugQmNlRCyxR3PLvfJiW7duHaqqMnDgQAYOHEhmZiajR4+mffv2VZbcKxIfwLJly1iyZAlz586tkrhOqmh8Us0ybtw4HA5HVDbaee655wBYvHgxkyZN4u23345wRH6bN29m27ZtjB49OqzjRiyxR3PLvfJiO13Tpk1JTU3l+++/p1u3bucxqlMqEt+qVat4+eWXmTVrFrVq1aqCqE45m59fNEhLS+Po0aNomoaqqmiaxq+//hq10x/RaOLEiRw4cIA33ngjqq8O+/Tpwz//+U9yc3Oj4ip748aN7N27ly5dugD+hkNDhw5l/PjxdOjQ4ZzHjcqpmGhvubdnzx4aNWoEQFZWFjt37iz5PBqsWbOG8ePHM3PmzLBe3tVUKSkpNGvWjKVLl9K7d2+WLl1Ks2bN5DRMBb300kts27aNt956C4vFEulwyigsLCQvL6/kTXr16tUkJCSQmJgY4cj87r33Xu69996Szzt37swbb7xBkyZNKjVutajHHm2J/cEHH2TPnj2YTCZUVWXYsGHcdNNNkQ6rRLt27TCbzWUS06xZs6LiDAVg6dKlTJo0iby8PMxmM3a7nXfffTeib457miYAqQAAAKdJREFU9+5lzJgx5OXlER8fz8SJE7n00ksjFk9pzz77LCtXriQnJ4ekpCQSExNZtmxZpMMCYPfu3fTs2ZOGDRtis9kAqF+/PtOnT49wZH45OTmMGDECp9OJoigkJCTwj3/8gxYtWkQ6tKAuqMQuSZIkVVz0ToZJkiRJ50QmdkmSpBpGJnZJkqQaRiZ2SZKkGkYmdkmSpBpGJnZJkqQaRiZ2SZKkGub/Ad0GWkGIoW/JAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"tags": []
}
}
]
},
{
"cell_type": "code",
"metadata": {
"id": "SMnU-9e-WZgv"
},
"source": [
"def test_e_step():\n",
" num_components = 3\n",
" X = create_dataset()\n",
" n = X.shape[0]\n",
" gmm = GaussianMixture(num_components=num_components, input_size=2)\n",
" q = gmm._expect(X)\n",
"\n",
" # Reference implementation\n",
" true_q = np.zeros((n, num_components))\n",
" for i in range(num_components):\n",
" mu = gmm.mu[i]\n",
" cov = gmm.sigma[i]\n",
" pi = gmm.pi[i]\n",
" true_q[:, i] = multivariate_normal.pdf(X, mean=mu, cov=cov)\n",
" true_q[:, i] = pi * true_q[:, i]\n",
" true_q = true_q/np.sum(true_q, axis=1).reshape(n, 1)\n",
" assert np.allclose(q, true_q)\n",
" \n",
"\n",
"def test_m_step():\n",
" num_components = 3\n",
" input_size = 2\n",
" X = create_dataset()\n",
" n = X.shape[0]\n",
" gmm = GaussianMixture(num_components=num_components, input_size=input_size)\n",
" q = gmm._expect(X)\n",
" pi_hat, mu_hat, sigma_hat = gmm._maximize(X, q)\n",
"\n",
" # Reference implementation\n",
" mu = np.zeros((num_components, input_size))\n",
" sigma = np.zeros((num_components, input_size, input_size))\n",
" pi = q.sum(axis=0)/n\n",
" for i in range(num_components): \n",
" mu[i] = np.sum(X * q[:,i].reshape((n, 1)), axis=0) / (pi[i]*n)\n",
" error = X - mu[i, :].reshape((1, input_size))\n",
" sigma[i] = np.dot((error * q[:,i].reshape((n, 1))).T , error) / (pi[i]*n)\n",
" assert np.allclose(pi, pi_hat)\n",
" assert np.allclose(sigma, sigma_hat)\n",
" assert np.allclose(mu, mu_hat)\n",
"\n",
"\n",
"def test_elbo():\n",
" num_components = 3\n",
" input_size = 2\n",
" X = create_dataset()\n",
" n = X.shape[0]\n",
" gmm = GaussianMixture(num_components=num_components, input_size=input_size)\n",
" q = gmm._expect(X)\n",
" elbo_hat = gmm.compute_elbo(X, q)\n",
"\n",
" # Reference implementation\n",
" elbo = 0\n",
" for i in range(num_components):\n",
" pi = gmm.pi[i]\n",
" mu = gmm.mu[i]\n",
" sigma = gmm.sigma[i]\n",
" inv_sigma = np.linalg.inv(sigma)\n",
" logdet_sigma = np.linalg.slogdet(sigma)[1]\n",
" constants = -input_size/2*np.log(2*np.pi)\n",
" for j in range(n):\n",
" distance = (X[j] - mu).T @ inv_sigma @ (X[j] - mu)\n",
" log_p_x_given_t = constants - logdet_sigma/2 -distance/2 \n",
" elbo += (np.log(pi) + log_p_x_given_t - np.log(q[j,i])) * q[j,i]\n",
" assert np.allclose(elbo, elbo_hat)\n",
"\n",
"test_e_step()\n",
"test_m_step()\n",
"test_elbo()"
],
"execution_count": null,
"outputs": []
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment