From 2c476f2ba2b827864c12644f659b0c1a1e90a955 Mon Sep 17 00:00:00 2001 From: Tom Birdsong Date: Wed, 17 Feb 2021 13:35:41 -0500 Subject: [PATCH 1/2] ENH: Added RegisterTwoPointSets example notebook --- .../RegisterTwoPointSets/Documentation.rst | 5 + .../RegisterTwoPointSets.ipynb | 814 ++++++++++++++++++ 2 files changed, 819 insertions(+) create mode 100644 src/Registration/Metricsv4/RegisterTwoPointSets/RegisterTwoPointSets.ipynb diff --git a/src/Registration/Metricsv4/RegisterTwoPointSets/Documentation.rst b/src/Registration/Metricsv4/RegisterTwoPointSets/Documentation.rst index a759d10f9..22755b10c 100644 --- a/src/Registration/Metricsv4/RegisterTwoPointSets/Documentation.rst +++ b/src/Registration/Metricsv4/RegisterTwoPointSets/Documentation.rst @@ -55,6 +55,11 @@ Python :language: python :lines: 1, 23- +Binder +------- +.. image:: https://mybinder.org/badge_logo.svg + :target: https://mybinder.org/v2/gh/InsightSoftwareConsortium/ITKExamples/master?filepath=src%2FRegistration%2FMetricsv4%2FRegisterTwoPointSets%2FRegisterTwoPointSets.ipynb + Classes demonstrated -------------------- diff --git a/src/Registration/Metricsv4/RegisterTwoPointSets/RegisterTwoPointSets.ipynb b/src/Registration/Metricsv4/RegisterTwoPointSets/RegisterTwoPointSets.ipynb new file mode 100644 index 000000000..7c5bb6287 --- /dev/null +++ b/src/Registration/Metricsv4/RegisterTwoPointSets/RegisterTwoPointSets.ipynb @@ -0,0 +1,814 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Register Two Point Sets\n", + "\n", + "Similar to image registration, an n-dimensional \"moving\" point set may be resampled to align with a \"fixed\" point set. An ITK point set metric may be employed with an ITK optimizer in order to register the two sets.\n", + "\n", + "In this example we create two `itk.PointSet` representations with an arbitrary offset and select parameters to align them with a `TranslationTransform`. We use the `JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4` class to quantify the difference between point sets and the `GradientDescentOptimizerv4` class to iteratively reduce this difference by modifying transform parameters. Our example also includes sample code visualizing the parameter surface with `matplotlib` and `itkwidgets` as well as a sample hyperparameter search to optimize gradient descent performance." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import sys\n", + "import itertools\n", + "from math import pi, sin, cos, sqrt\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "import itk\n", + "from itkwidgets import view" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Construct Two Point Sets" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# Generate two circles with a small offset\n", + "def make_circles(dimension:int=2,offset:list=None):\n", + "\n", + " PointSetType = itk.PointSet[itk.F, dimension]\n", + "\n", + " RADIUS = 100\n", + "\n", + " if not offset or len(offset) != dimension:\n", + " offset = [2.0] * dimension\n", + " \n", + " fixed_points = PointSetType.New()\n", + " moving_points = PointSetType.New()\n", + " fixed_points.Initialize()\n", + " moving_points.Initialize()\n", + "\n", + " count = 0\n", + " step = 0.1\n", + " for count in range(0, int(2 * pi / step) + 1):\n", + "\n", + " theta = count * step\n", + " \n", + " fixed_point = list()\n", + " fixed_point.append(RADIUS * cos(theta))\n", + " for dim in range(1,dimension):\n", + " fixed_point.append(RADIUS * sin(theta))\n", + " fixed_points.SetPoint(count, fixed_point)\n", + "\n", + " moving_point = [fixed_point[dim] + offset[dim]\n", + " for dim in range(0,dimension)]\n", + " moving_points.SetPoint(count, moving_point)\n", + "\n", + " return fixed_points, moving_points" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "POINT_SET_OFFSET = [15.0, 15.0]\n", + "fixed_set, moving_set = make_circles(offset=POINT_SET_OFFSET)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeA0lEQVR4nO3df5BeVXkH8O9DSHRxalbZGJJNYjJOhhZDivKKOqEdJMYADgR1WGmnhVE7caYyTpkOuKltTLGWNUybNlYdEak4rdLVgbhKMEYy6uCUkY2hIWAzpILDbvgh1gTHpIYkT/947yX3fXPvu+89995zz7nn+5nJbPa+e3dvLsvznvM85zlXVBVERBSWM+q+ACIiso/Bn4goQAz+REQBYvAnIgoQgz8RUYDOrPsC+jE0NKRLly6t+zKIiLyye/fuF1R1XtprXgT/pUuXYnJysu7LICLyioj8POs1pn2IiALE4E9EFCAGfyKiADH4ExEFiMGfiChADP5EZds7DmxZAWwabH/cO17teUQGvFjqSVSLvePAA7cAh6eAuYuA1RuBlSMzn/OtjwIvHW1/fvjp9udA73NNzyMyxJE/UZo4GB9+GoCeCsYzjcYfuOVUAI+9dLR9vIrz4mvljIFyYvCn5jMJjqbB+PBUvuNFzzN9k6LgMfhTs5kGR9NgPHdRvuNFz+OMgQwx+JNf8gYs0+BoGoxXbwRmD3Qemz3QPl7FeZwxkCEGf/KHScAyDY6mwXjlCHDlVmDuYgDS/njl1pmLtqbnccZAhrjah/zRK2BlBcm5i6I3i5TjvcTfL+9qn/hckxU6Juet3ti5SgiwM2PgqiTvMfiTP0wClmlwBMyDuE2mb1Kmb4omb8DkJAZ/qofJGnqTgFVkBO8LH2YM5BwGf7LPNHVgGrB8GMHbZnvGYPJmT5Vi8Cf7TFMHIYzibbI1Y2CdwEkM/mRfkdQBR/H1MnkDZp3ASQz+VI4803rT1AG5Ie8bMOsETiplnb+I3Ckiz4vIvsSx14rIThF5Ivr4mui4iMhWETkgIntF5M1lXAPVKO/6e9M19BXYtmcaq8Z2YdnofVg1tgvb9kxXel6QTHsRqFJlNXl9GcBlXcdGATygqssBPBB9DgCXA1ge/VkP4PMlXQPVJW/DkGlDUw8mwXjbnmlsuOdRTB86CgUwfegoNtzz6Iznmp5nep3eM3mzZyNZ5UpJ+6jqD0VkadfhdQAuif5+F4DvA/hYdPwrqqoAHhKRQRFZoKrPlHEtVAOTaX1G6mDbnmnctmM/Dh46ioWDA7hp7bm4+k3DPX98HIyPvnQCwKlgDKDnubft2P/yObGjL53AbTv2V3Ke6XV6L2+dgAViK6rM+c9PBPRnAcyP/j4MIJnwnYqOdQR/EVmP9swAS5YsqfAyqbCScvi2g/jBQ0dzHS96nul1AmZvik7JUydggdgKK3v7RKN8zXnO7araUtXWvHnzKroyypRn2l1SDr9XcOzFNBgvHBzIdbzoeabXWSTN5CUWiK2oMvg/JyILACD6+Hx0fBrA4sTXLYqOkSvyFnBTcvgPn/+3WLV9KFdu23YQv2ntuRiYPavj2MDsWbhp7bmVnGd6naZvit5igdiKKoP/BIDro79fD+CbiePXRat+3gbgMPP9jjHZ8XHlCHDjPmDTIWy7ZAeue/j1uUeqtoP41W8axq3vPR/DgwMQAMODA7j1vefPmE4xPc/0OovMGLwsLuedSbI4bKSUnL+IfA3t4u6QiEwB+ASAMQDjIvIhAD8HECfrtgO4AsABAEcAfKCMa6ASFZx2m+a2b1p7bkfOH+g/iMc/N29O/Oo3DRvlzk3OM73OhYMDmE4J9L3eFL0uLucpELM4bEza6Xi3tVotnZycrPsywrFlRWoBd+rkEN5/1hdnDFjLRu9LLfAIgCfH3t3zR3tf2KxAdyAH2m+KvWYbq8Z2pb5hDA8O4Eejl1Z2rdZl/K5i7uL2TDRwIrJbVVtpr7HDl06Xsn/LEZ2DzcdH+hpBmoxUY6Yj8SYzmTGYpoq8w+KwMQb/kPS7BUNi2n3y8BQOnjwbm4+PYOLkxQBmTuGYpm8oW943RdM3YO9mXtwqxBgf4xgKkxU8N+7DG/7v33Hxsa0vB/5YrxGkaUGUymNSXPZySalDW4X4hiP/UBg2zpiOIJm+qZdJqqhIE1ptTLqHuSU4AAb/cGTkQE8ensLEnmmmcBoo7xuwt3WCfruHuTKoA9M+ocjIgR48eXbPqT1TOOEw7bPwhkn/SoNx5B+KHit4jp7sPbVnCicMprM8b4rEXBnUgcG/CfrJY0afT31jAxbKL3FQO1fwOD+1p8qZ1Am8aibjyqAODP6+y5PHXDmC928fMl6DT82Xd5bnVZHY5PnDDcacv+9y5jFN95chSuNVkbiChwj5jCN/z+nhKUiO40X2wSHqVqSbuxZ5VgY1fEkog7/nnsMQzsEvMo6nYwGXypK3SOxFcTiQJaFM+3ju1mPX4IjO6Th2ROfg1mPX1HRFFJI8S4G96SAOZEkoR/6em3z1Goy+CNx85njHKp7dr15T96VRIPqdSXpTHA5kSSiDv6Pi6XHrxZ3YMOfrmI8XICm5x/a0+xgmjp3ae2dg9izcygIuOcab4nAgS0KZ9nFQPD2+8MWduHX2HTgHv4BkbMbGDlzyhTcdxIFsFseRv4Pi6fHNc8ZxlhzrfDFlMzYWcMkH3uwTlXezOE8x+DsongYvlBfSv6BhuUcKQ95lxrWuDJppSWgDloIy+DsoXjt9UIewKO0NoGG5RwpHv7NUp7eNaMhSUOb8HRR34W4+PnLaMs4m5h6JuvVaGVS7hiwF5cjfQaemx3Ow4UX0XO1D1EROrwxqyFJQBn/L+s1jnpoeXwrgVuvXSVQnp7eNaMhSUKZ9LErrcHzw3s/hyKd/F9g0CGxZkf1MXaKAOL0BYUOWgnLkb1F3HvOqMx7ELXIHzjoaLef0tHBEVDanNyBsyFJQUdW6r2FGrVZLJycn676MwpaN3ofk3X5wzkex6Iy01TyLgRv3WbsuIl95sVFcjURkt6q20l5j2sei7nwl1/ETmXN6o7i94+00rsPpXAZ/i7rzmAd1KP0LPSscEdXB2eWgcR/A4aeBjG1ZXMDgb1H3Pjx3zPkTHJ/1ys4v8rBwRFQHZ5eDetIHwIJvifrJP3Z2OL4b2PtG7wtHRHVwdjmoJ30ADP4lMW5H7/exckTUwdmN4jzpA2DapyQ9848eFH+IfOPsduae9AFw5F+SrDxj68WdwLf+1ftNoIhc1M9GcdaXg3rSB8DgX5Ks/OOGOV/PLv449stA1DS17Q7qQTqXaZ+SZLWjzwfX8hPVxbnloA6lgBn8S5KVf5SsIo9jxR+iJnJqOahj6/+Z9ilRav5x1sbOBz8AThZ/iJrIqeWgvdb/15Ai4sg/h217prFqbBeWjd6HVWO7+msjXzkCXLm1vV8PpP3xyq3O5wOJmsCp3UEdW/9f+chfRJ4C8GsAJwAcV9WWiLwWwH8AWArgKQAjqvqrqq+liEKFIw+KP0RN5NTuoI6t/7eV9nmHqiYrn6MAHlDVMREZjT7/mKVrMdKrcNTxi9SABzsTNUm/zw2u3Gq3UsB15fzXAbgk+vtdAL4Px4N/X4WjhjzYmSgUVnsAHFv/byP4K4DviogC+IKq3g5gvqo+E73+LID53SeJyHoA6wFgyZIlFi6zt74KR44VdIgoWy09AA6lgG0UfC9W1TcDuBzAR0TkD5MvavtpMqc9UUZVb1fVlqq25s2bZ+Eye+urcORYQYeIsjnXA2BZ5cFfVaejj88DuBfARQCeE5EFABB9fL7q6yiqr31EuKafyBvO9ADU1PhVadpHRF4F4AxV/XX093cBuAXABIDrAYxFH79Z5XWUZcbCkWMFHSLK5kQPQI11wqpH/vMBPCgi/wXgxwDuU9XvoB3014jIEwDeGX3uBKO1/DGu6SfyhhM9ADU++KXSkb+q/gzA76cc/yWA1VX+bBOlFIAcKugQUTYnegBqrBNye4eEvtfyE1Ej1N4DUGPjF7d3SDAqADm0Sx8RladQCrhfNT74hSP/hNwFIDZ1ETWStR6AGhu/OPJPyF0AqrFYQ0TVsdoDsHIEuHEfsOlQ+6OlgSNH/gm5C0Bs6iJqJGd6ACrE4N8lVwHIsV36iKgctfUAWNwYMui0T+GCTo3FGiKqTi09AJaf9BVs8I8LOtOHjkJxqqDDpi4i6ms7l7JZriEGm/YpbU0/m7qIGsl6D4DlGmKwI/8QCjpE5BHLG0MGG/yzCjd9FXTY2EUUnMqbvizXEIMN/sYFHctFGSKqXyk1wplYriFK+1kqbmu1Wjo5OVn69zV6hNuWFRnLOxe3GzSIqHFWje1KXfo5PDiAH41eWsMV9UdEdqtqK+21YAu+gGFBh41dRMFpYo0w2LSPMT6tiyg4hWqEJizUFYMK/qUUbNjYRRQcq01fluqKwQT/0go2bOwiCo7Vpi9LzV7B5PxLfVALG7uIgmOt6ctSXTGYkX8TCzZE1ECW6orBBP/SCjZs8CIiVNj0ZamuGEzwL6VgwwYvIkLFTV+W6opBNXkZNXUlscGLiOBP0xebvCKFCzZs8CIiNKOGGEzapxRs8CIiWGr6qri+2PjgX2pRhg1eRAQLTV8W6ouNDv6lF2XY4EVEsND0ZaHRq9E5/1Ibu2Js8CIiVNz0ZaG+2OiRfxOKMkQUIAv1xUYHf+s78RFRcCpp9rJQX2x08K+kKMMOXyKKVNbsZaG+2Oicf5yPK9TYlRRX4ONCTFyBB1gHIApQJXXFWMX1xUYHf6DkokyvCjyDP1FwKq0r7h1vx5bDU+1c/+qNpcaZRqd9SscOXyJKqKyuyHX+xZReiGGHLxElVNbsZWGdf2ODfyWFGHb4ElFCZc1eFrIMjc35V9bgBVSahyMiv1TS7DV3UcYOwuVlGRob/CsrxLDDl4iqtnpj58pCoDnr/EXkMhHZLyIHRGS07O/PBi8isqX0+mJT1/mLyCwAnwWwBsAUgIdFZEJVHy/rZ9y09lxsuOfRjtRP6bvuMf1DFLy4vhjHmri+CKBYOqjiLENdI/+LABxQ1Z+p6jEAdwNYV+YPqHTXPT7OkYgiveqLLqsr5z8MIFnNmALw1uQXiMh6AOsBYMmSJUY/pLJd99jsRUSR0uuLlrIKzi71VNXbVbWlqq158+bVfTmd2OxFRJFS64sWswp1Bf9pAIsTny+KjpWqkt32ADZ7EdHLSm30stDcFasr+D8MYLmILBOROQCuBTBR5g+obLc9gM1eRPSyUuuLFrMKteT8VfW4iNwAYAeAWQDuVNXHyvwZle+2B3C1DxEBKLG+aKG5K1Zbk5eqbgewvarvX/lTvNjsRURoZxlK2zbeQnNXzNmCb1Fs8iKiqpWeXrbQ3BVr7PYOlTd5xdjsRRSsyvYQsxBDGhv8S3+KVxo+2YsoaJWnlyvU2OAPVNjkFWOzF1HQFg4OYDol0OdKL9eUPWhszt8KNnsRBa3wGv8at4pp9MgfKLkS383isiwick/h9HKN2YNGB//KdtuLWVyWRURuKpRerjF70Oi0T+W77VlclkVEDVTjVjGNHvlbqcSz2YsoSKWklGvMHjR65F9bo9fecWDLCmDTYPsj9/knapTSmrtqzB40euRvrdEriWv/iRqv1OaumrIHjR75V/o0rywWt2Qlonr43NwVa/TIH7DQ6NWNa/+JGq9Qc5cjW8I0euRfCz7ohajxjJu7HHr+d1DBv7IneyXxQS9EjWecUnYoLdz4tE+s8oavGB/0QhQEo5SyQ2nhYIJ/pU/26sa1/0SUxqEtYYJJ+zShOk9E9SqcOnYoLRxM8K/9yV5s/CLyWimNXQ5tCRNM2qeWhq8YG7+IvFda6tiRtHAwI/9aGr5iDlX4ichM7tSx47P9YEb+QA0NXzGHKvxEZCZXY5cHs/1gRv61YuMXkfdyNXZ5MNsPauSfptInfcX40Bci7+V6apcHs/2ggz8bv4goj75Txw6t588SdPBn4xcR9WKcGfBgth908Hem8cuRXf6I6JRCmQEPZvtBB/9C27KWxYNVAUQhKpwZcHy2H/RqH+NtWcvkwaoAohD1lRlwfC1/L0GP/HNV76viwaoAohDNmBnwfNYedPAHamz8inmwKoAoRDNuCdNr1u5B8A867eMEh3b5I6JTZtwSxvNZe/Aj/zRWGr9iHqwKIApVz8yA57N2Bv8u1hq/kmZaFcCloETu8WAtfy8M/l2sNn71w/OiEpEPjGb7ns/aGfy7ONP4FfO8qETkuqzZ/vDT38Zb/uczvQO742v5e2Hw7+JE41eS50UlItelzfbXnPgBVvzkSwB+2z7QwBl3Zat9RGSTiEyLyCPRnysSr20QkQMisl9E1lZ1DSacaPxK4nbQRJVKm9XffOY4BuLAH2tY82XVSz23qOoF0Z/tACAi5wG4FsAbAVwG4HMiMqvXN7Gp1id+peFSUKJKpc3qF8oL6V/coBl3HWmfdQDuVtXfAnhSRA4AuAjAf9ZwLalqb/xK8ryoROS6tGauZzCEYaS8ATRoxl118L9BRK4DMAngL1X1VwCGATyU+Jqp6FgHEVkPYD0ALFmypOLL7J/VHoBYv0UlLgklyi1tm5eD592M4Uc/4e0yzn6IqpqfLPI9AOekvPRxtAP8CwAUwCcBLFDVD4rIvwB4SFX/LfoeXwJwv6p+I+vntFotnZycNL7OsnSvCgDa9YBa00Kx7iWhQPuX9cqtfAMgMtGAwZSI7FbVVtprhUb+qvrOPi/giwC+HX06DWBx4uVF0THnOdcDkMQloUSZumfs/3TeE41extmPKlf7LEh8+h4A+6K/TwC4VkReISLLACwH8OOqrqNMzvUAJHFJKFGqeMY+fegoFMCFL+7Eit1/HW3NoKeWcXq0HXMZqlzts1lEHhWRvQDeAeBGAFDVxwCMA3gcwHcAfERVT2R/G3dkrfWvrQcgiUtCiVJ1z9hvPnMcA3Ks84satoyzH5UFf1X9U1U9X1VXqupVqvpM4rVPqeobVPVcVb2/qmsom3M9AElcEkqUqntmHsIyzn5wS+ccnOsBSFo50i7uzl0MQNofWewlOm1mflCH0r8wsFlyodU+triy2qcftSwFzasBqxiI+tW9Su+qMx7Ep2ff0Zn6aejKuF6rfTjyL1F3YSneIGrbHocWM8VLQgMvdlEDZTxPt3vGvvvVa7Dvwr8LfpbMkX+JVo3tSt0UbnhwAD8avbSGK0qxZUXGAygWAzfuO/04kQcenvgCVvzkbzr342noaD4PjvwtcXopaIxLQqlhtu2ZxsLdmxu/EVvZGPxL5PRS0BiXhFLD3LZjPxak7cMDcFDTA4N/iZxeChrLuyQ0I49K5IqDh45yBY8BBv8SOb0UNJZnSSiLw+SBhYMD2Hx8BEd0Tsfxo3gF+1x6YMG3Bl4sBwVYHKb69bEsOV5lt+bED3DzmeNYKL/EMzgbBy+8GW+56sM1XbgbKtvYjfLLel4oAPfeAFgcpjp171Sb8SjFU1syz8EfHLrY7QGVQxj8LXN6Z9BucxdljPyZRyULMnaqPXL/RpzVNfp36gFMnmDO3zIvloPGTPYLYoGYypIxw3zlkWfdapz0FIO/ZV4sB43l3S+IBWKaSZ7BQcYM86Cejdt27K/m+gLC4G+ZF8tBk1aOtIu7mw61P/bqluz1QBmivIOD1RtPW8FzROdg8/ERN2fKnmHwt8yL5aCmWCCmXvIODlaOYPPsP8fUySGcVMHUySGMvvRnmDh5sZszZc+w4FuDPMUpb5aFAuYFYu4yGgaDwcEF716PNfe8HUePdT4329mZskc48neYF7uEJpkWiFkn8FPe4n7GIOBZDGX+Tjd6plwzjvwd5tWyUODUaD3PKJ4PnvdTn2vwO6ze2HkO2jn8v3/pGuzs0evCZZzVYPB3mFfLQmMrR/IFbdYJ/GTyph0df/aev8Lr9AUc1LOx+fgIJk5eDJx0eFDTUAz+Dls4OJD6fIBGFbtYJ3BD3vtp+qa9cgRv/+qrkLapjNODmgZizt9heZeFbtszjVVju7Bs9D6sGtvlbm0giXWC+pnczwJbg3vV69JgDP4Oy1Ps8q44HDN58HyRfgJ2IJ/O5H6mvGkfn/VKbPrN+2YcfHjX69JQTPs4rt9il3fF4SRbdQKTIqVvTNJhJvezq7h/ZOAcbPzN+/CNYxcB6L1h4amN2DxZwtxQDP4N4WVx2JRpnaDIyiLTGoPN80zf3EzvZ+JNe83YLkwf67y3vQYfXMFTP6Z9GiKoPKpJnQAoPmPIW2OwfZ5pOsz0fiYENfhoCAb/hjDJo3pZIAbM6gSAeZHSNKjaPq/AChyj+5kQ1OCjIZj2aYi8eVSvHiqTJm+dAEhtMqp0xmD7vCLPX0i5n3m2Frlp7bkdv08Ai7iuY/BvkDx5VK8LxKZMOpAB86Bq+zzTN7cUeQcHLOL6h8E/UMHmaG3OGGyfZ/rmlsJkcMAirl8Y/ANl0j3s1Q6jZTINqrbPi88tYdlqsIODgDD4Bypvjtb7GkFRpkHV9nkZ8r5xB7G1SOC42idQebfK7ZUGILeZdH+zC7f5OPIPWJ4cbZE0QLDpIkeY5u/jc/nfrZkY/KkvpmmA4NNFJTN5IzV942YBt9mY9qG+mKYBmC4qj+nmfWzAojQM/tQX08fpmY46ve0+ziHvv9H0jZT5e0rDtA/1zSQNYLqk1DRV5Et9weTfWCR9AzB/T50KBX8RuQbAJgC/B+AiVZ1MvLYBwIcAnADwUVXdER2/DMA/A5gF4A5VHStyDeQ2k7Z/0+7jut40TM41+TcWWX7J/D11K5r22QfgvQB+mDwoIucBuBbAGwFcBuBzIjJLRGYB+CyAywGcB+CPoq+lhjJJF5mOcE3TIkUehGN6rsm/kekbKlOhkb+q/hQARKT7pXUA7lbV3wJ4UkQOALgoeu2Aqv4sOu/u6GsfL3Id5La8o07TEW4VbxpVjOABs38j0zdUpqpy/sMAHkp8PhUdA4Cnu46/Ne0biMh6AOsBYMmSJRVcIrnKdIdI228aRc41/TcyfUNlmTHtIyLfE5F9KX/WVXlhqnq7qrZUtTVv3rwqfxQ5xnRlkWlapMhSSNNzTf+NRGWZceSvqu80+L7TABYnPl8UHUOP40QvMxnhmqZFiuxFX+RcjuKpTlWlfSYAfFVE/hHAQgDLAfwYgABYLiLL0A761wL444qugQJk802j6LlEdRJVNT9Z5D0APgNgHoBDAB5R1bXRax8H8EEAxwH8hareHx2/AsA/ob3U805V/dRMP6fVaunk5ORMX0ZERAkisltVW6mvFQn+tjD4ExHl1yv4c3sHIqIAMfgTEQWIwZ+IKEAM/kREAfKi4CsivwDw8wLfYgjACyVdTpPwvmTjvcnGe5PNtXvzelVN7ZL1IvgXJSKTWRXvkPG+ZOO9ycZ7k82ne8O0DxFRgBj8iYgCFErwv73uC3AU70s23ptsvDfZvLk3QeT8iYioUygjfyIiSmDwJyIKUKOCv4hcIyKPichJEWl1vbZBRA6IyH4RWZs4fll07ICIjNq/avtEZJOITIvII9GfKxKvpd6nkIT4O5FFRJ4SkUej35PJ6NhrRWSniDwRfXxN3ddpg4jcKSLPi8i+xLHUeyFtW6Pfob0i8ub6rjxdo4I/+ED5PLao6gXRn+1A9n2q8yJtC/x3Iss7ot+TeEA1CuABVV0O4IHo8xB8Ge3/L5Ky7sXlaD/HZDnaj6P9vKVr7Fujgr+q/lRV96e89PID5VX1SQDxA+UvQvRAeVU9BiB+oHyosu5TSPg7MbN1AO6K/n4XgKvruxR7VPWHAP6363DWvVgH4Cva9hCAQRFZYOVC+9So4N/DME5/cPxwj+MhuCGajt6ZmLaHfD9ivAedFMB3RWS3iKyPjs1X1Weivz8LYH49l+aErHvh/O9RVY9xrIyIfA/AOSkvfVxVv2n7elzV6z6hPQX9JNr/Y38SwD+g/dQ1om4Xq+q0iLwOwE4R+e/ki6qqIsL14vDvXngX/PlA+f70e59E5IsAvh192us+hYL3IEFVp6OPz4vIvWinxZ4TkQWq+kyUyni+1ousV9a9cP73KJS0zwSAa0XkFdHD4+MHyj+M6IHyIjIH7WLnRI3XaUVX7vE9aBfKgez7FJIgfyfSiMirROR34r8DeBfavysTAK6Pvux6ACHPuLPuxQSA66JVP28DcDiRHnKCdyP/XroeKH+fiDyiqmtV9TERGQfwONoPlP+Iqp6IzrkBwA6ceqD8YzVdvk2bReQCtNM+TwH4MAD0uk+hUNXjgf5OpJkP4F4RAdqx4quq+h0ReRjAuIh8CO2t1kdqvEZrRORrAC4BMCQiUwA+AWAM6fdiO4Ar0F40cQTAB6xf8Ay4vQMRUYBCSfsQEVECgz8RUYAY/ImIAsTgT0QUIAZ/IqIAMfgTEQWIwZ+IKED/D7rNmcFJZkmGAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Visualize point sets with matplotlib\n", + "\n", + "fig = plt.figure()\n", + "ax = plt.axes()\n", + "\n", + "n_points = fixed_set.GetNumberOfPoints()\n", + "ax.scatter([fixed_set.GetPoint(i)[0] for i in range(0,n_points)],\n", + " [fixed_set.GetPoint(i)[1] for i in range(0,n_points)])\n", + "ax.scatter([moving_set.GetPoint(i)[0] for i in range(0,n_points)],\n", + " [moving_set.GetPoint(i)[1] for i in range(0,n_points)])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Run Gradient Descent Optimization\n", + "\n", + "We will quantify the point set offset with `JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4` and minimize the metric value over 10 gradient descent iterations. " + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "ExhaustiveOptimizerType = itk.ExhaustiveOptimizerv4[itk.D]" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "dim = 2\n", + "\n", + "# Define translation parameters to update iteratively\n", + "TransformType = itk.TranslationTransform[itk.D, dim]\n", + "transform = TransformType.New()\n", + "transform.SetIdentity()" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "PointSetType = type(fixed_set)\n", + "\n", + "# Define a metric to reflect the difference between point sets\n", + "PointSetMetricType = itk.JensenHavrdaCharvatTsallisPointSetToPointSetMetricv4[PointSetType]\n", + "metric = PointSetMetricType.New(\n", + " FixedPointSet=fixed_set,\n", + " MovingPointSet=moving_set,\n", + " MovingTransform=transform,\n", + " PointSetSigma=5.0,\n", + " KernelSigma=10.0,\n", + " UseAnisotropicCovariances=False,\n", + " CovarianceKNeighborhood=5,\n", + " EvaluationKNeighborhood=10,\n", + " Alpha=1.1)\n", + "metric.Initialize()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "# Define an estimator to help determine step sizes along each transform parameter2\n", + "ShiftScalesType = itk.RegistrationParameterScalesFromPhysicalShift[PointSetMetricType]\n", + "shift_scale_estimator = ShiftScalesType.New(\n", + " Metric=metric,\n", + " VirtualDomainPointSet=metric.GetVirtualTransformedPointSet(),\n", + " TransformForward=True\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "max_iterations = 10\n", + "\n", + "# Define the gradient descent optimzer\n", + "OptimizerType = itk.GradientDescentOptimizerv4Template[itk.D]\n", + "optimizer = OptimizerType.New(\n", + " Metric=metric,\n", + " NumberOfIterations=max_iterations,\n", + " ScalesEstimator=shift_scale_estimator,\n", + " MaximumStepSizeInPhysicalUnits=8.0,\n", + " MinimumConvergenceValue=-1,\n", + " DoEstimateLearningRateAtEachIteration=False,\n", + " DoEstimateLearningRateOnce=True,\n", + " ReturnBestParametersAndValue=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "iteration_data = dict()\n", + "\n", + "# Track gradient descent iterations with observers\n", + "def print_iteration():\n", + " print(f'It: {optimizer.GetCurrentIteration()}'\n", + " f' metric value: {optimizer.GetCurrentMetricValue():.6f} '\n", + " #f' transform position: {list(optimizer.GetCurrentPosition())}'\n", + " f' learning rate: {optimizer.GetLearningRate()}')\n", + " \n", + "def log_iteration():\n", + " iteration_data[optimizer.GetCurrentIteration() + 1] = list(optimizer.GetCurrentPosition())\n", + " \n", + "optimizer.AddObserver(itk.AnyEvent(), print_iteration)\n", + "optimizer.AddObserver(itk.IterationEvent(), log_iteration)\n", + "\n", + "# Set first value to default transform position\n", + "iteration_data[0] = list(optimizer.GetCurrentPosition())" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "It: 0 metric value: 0.000000 learning rate: 1.0\n", + "It: 0 metric value: -0.043464 learning rate: 5594.753388446298\n", + "It: 1 metric value: -0.054787 learning rate: 5594.753388446298\n", + "It: 2 metric value: -0.062597 learning rate: 5594.753388446298\n", + "It: 3 metric value: -0.064588 learning rate: 5594.753388446298\n", + "It: 4 metric value: -0.064807 learning rate: 5594.753388446298\n", + "It: 5 metric value: -0.064815 learning rate: 5594.753388446298\n", + "It: 6 metric value: -0.064815 learning rate: 5594.753388446298\n", + "It: 7 metric value: -0.064815 learning rate: 5594.753388446298\n", + "It: 8 metric value: -0.064815 learning rate: 5594.753388446298\n", + "It: 9 metric value: -0.064815 learning rate: 5594.753388446298\n", + "It: 10 metric value: -0.064815 learning rate: 5594.753388446298\n", + "Number of iterations: 9\n", + "Moving-source final value: -0.06481531061643396\n", + "Moving-source final position: [15.000412861069881, 14.99997463945473]\n", + "Optimizer scales: [1.0000000000010232, 1.0000000000010232]\n", + "Optimizer learning rate: 5594.753388446298\n", + "Stop reason: GradientDescentOptimizerv4Template: Maximum number of iterations (10) exceeded.\n" + ] + } + ], + "source": [ + "# Run optimization and print out results\n", + "optimizer.StartOptimization()\n", + "\n", + "print(f'Number of iterations: {optimizer.GetCurrentIteration() - 1}')\n", + "print(f'Moving-source final value: {optimizer.GetCurrentMetricValue()}')\n", + "print(f'Moving-source final position: {list(optimizer.GetCurrentPosition())}')\n", + "print(f'Optimizer scales: {list(optimizer.GetScales())}')\n", + "print(f'Optimizer learning rate: {optimizer.GetLearningRate()}')\n", + "print(f'Stop reason: {optimizer.GetStopConditionDescription()}')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Resample Moving Point Set" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "moving_inverse = metric.GetMovingTransform().GetInverseTransform()\n", + "fixed_inverse = metric.GetFixedTransform().GetInverseTransform()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "transformed_fixed_set = PointSetType.New()\n", + "transformed_moving_set = PointSetType.New()\n", + "\n", + "for n in range(0,metric.GetNumberOfComponents()):\n", + " transformed_moving_point = moving_inverse.TransformPoint(moving_set.GetPoint(n))\n", + " transformed_moving_set.SetPoint(n,transformed_moving_point)\n", + " \n", + " transformed_fixed_point = fixed_inverse.TransformPoint(fixed_set.GetPoint(n))\n", + " transformed_fixed_set.SetPoint(n,transformed_fixed_point)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZxklEQVR4nO3dfYxcV3nH8e8vBkeIFxPqTWySuE5IgmpcGuxNilSI2iaFxGptQtXUoYJAUB2kIEFpTUNjgduCBKG8iEKBICLSFhIoEGEh0kIQhVZqwLvBGIcQsEOsjWt7NwQIKlZCyNM/5q5zdz2zu3Pv3PffR1rt7JmZvcd3x8+cec557lFEYGZm3XJS1R0wM7PyOfibmXWQg7+ZWQc5+JuZdZCDv5lZBz2p6g4sxcqVK2Pt2rVVd8PMrFEmJycfjIixfvc1IvivXbuWiYmJqrthZtYokg4Ous9pHzOzDnLwNzPrIAd/M7MOcvA3M+sgB38zsw4aSfCXdJOkaUn7Um3PkvQVST9Mvp+StEvSByTtl7RX0oZR9MGsSLt3fZQjO8/h8bet4MjOc9i966OlPt9s1EY18v8EcOm8tuuAr0bEucBXk58BLgPOTb62AR8eUR/MFpQ1AO/e9VHWT+5gFTOcJFjFDOsnd5T2fLMijCT4R8Q3gIfmNW8Bbk5u3wy8LNX+z9FzJ/BMSatH0Q+zQfIE4DPvejdP0aNz2p6iRznzrncv6dh5n+9PDVaEInP+p0XE4eT2EeC05PbpwFTqcQ8kbXNI2iZpQtLEzMxMgd20JskaCPME4FOj/+vv1HhwScfO83x/arCilDLhG70dY4baNSYiboyI8YgYHxvrW51sHZMnEOYJwNPq//qb1spFn5v3+f7UYEUpMvgfnU3nJN+nk/ZDwJmpx52RtJktKE8gzBOApzZs51gsn9N2LJYztWH7os/N+3x/arCiFBn8dwFXJbevAr6Qan9VsurnhcDPUukh64gsI9I8gTBPAL5g8zXs2/h2jjDG4yGOMMa+jW/ngs3XLPrcvM+v8lODtdtILuwm6Rbgd4GVkh4A3ga8E/iMpNcCB4Erkod/CdgE7Ad+AbxmFH2w5pgdkT5Fj0IyIl0xuYPdsGBAnNYYqzjxDWBaK1m1yDEv2HwNu+kFxFPjQaa1kqmN24cK4CSPXZV8DSPr86c2bGfF7LlKHIvlTG3cvujvODVmQP3al/apoXeuZpjWGFMbln6urBlGEvwj4soBd13c57EBXDuK41ozLTgiXSDA5AmEkD+AVyHPm1bWN8usb87WLI24pLPVU9bRYdYRad7Re1OV/akh65uzNYuDv2WSZ3SYN33TtNF7VbK+WeZJF1lzOPhbJnlGh3nTN7Z0Wd4s87w5e66gOXxhN8skz8qbvKtnrFhZV0Z5aWmzeORvmeQZHYLTN3WWNV3kuYJmcfA3YPiP607dtFuWN2fPFTSL0z6W6eO6Uzc2X97LYFi5PPK3zB/XnbqxtKyfBj1JXA0Hf/PHdRuJLHMFLiirjoO/5Z68NZs17KdBTxJXxzn/Fhr2oml5r1ppllXevRIsOwf/lvHkrTWJJ4mr47RPy3jy1prES4ar4+DfMp68tSbJOkns1UH5Ofi3jCdvrWmG+dTp1UGj45x/y3jy1trMu5ONTqEjf0nPBT6dajobeCvwTODP4fgQ9W8i4ktF9qXJhvmY29Vr3ls3OK05OoUG/4i4FzgfQNIyehu130Zv68b3RcQ/FHn8NsjyMdeTt9ZWTmuOTplpn4uBAxFxsMRjNp4/5po9wWnN0Skz+G8Fbkn9/HpJeyXdJOmU+Q+WtE3ShKSJmZn+hSBd4CIYsydkqUkZtuixK9TbT73gg0jLgf8FnhcRRyWdBjwIBPD3wOqIuHrQ88fHx2NiYqLwftbRkZ3n9P2Ye4QxVu3cX0GPzJpjTto0cSyWd6aIUdJkRIz3u6+skf9lwF0RcRQgIo5GxK8i4nHgY8CFJfWjcfwx1yw7p00HKyv4X0kq5SNpdeq+y4F9JfWjcXzpBbPsnDYdrPAiL0lPBf4ASEerGySdTy/tc/+8+zpjqUs4vXrHLBuvDhqs8OAfEf8H/Nq8tlcWfdy6c6WiWfF87aDBXOFbEecizYrntOlgvrZPRVypaFaOpaZNu3bBOI/8K+LrmJvVR5Z9MJrOwb8iXsJpVh9dTMM6+FfEuUiz+ujiklDn/CvkJZxm9dDFJaEe+RfA1xIxa5YupmE98h8xr983a54u7oNRyoXd8mrShd18ITYzq4uFLuzmkf+Ief2+WXu1qRbAOf8R8/p9s3ZqWy2Ag/+IdXHiyKwL2lYL4OA/Yl6/b9ZObasFcM6/AF6/b9Y+basF8Mh/SF7Db9ZNbUvplrGZy/3Az4FfAY9FxLikZwGfBtbS28zlioj4SdF9yctr+M26q221AIWv80+C/3jEE4kxSTcAD0XEOyVdB5wSEX896HfUZZ2/1/CbWZPUYQP3+bYANye3bwZeVlE/htK2CR8z664ygn8AX5Y0KWlb0nZaRBxObh8BTiuhH7l5Db+ZLaRJc4JlBP8XRcQG4DLgWkkXpe+MXt7phNyTpG2SJiRNzMz0H3GXrW0TPmY2Ok0rAis8+EfEoeT7NHAbcCFwVNJqgOT7dJ/n3RgR4xExPjbWf8RdNq/hN7NBmlYEVuhqH0lPBU6KiJ8nt18C/B2wC7gKeGfy/QtF9mOUvIbfzPpp2nW9il7qeRpwm6TZY30qIv5d0m7gM5JeCxwErii4H2ZmhWpaEVihwT8i7gN+q0/7j4GLizx2Hm26cp+ZlWNqw3ZWzNYBJY7FcqY2bu9e8G8iF3KZWRZNKwLzZi7zuJDLzNqijkVeteVCLjPrAgf/eVzIZWZFqFsBmIP/PC7kMrNRq2MBmIP/PC7kMrNRq2MBmFf79OFCLjMbpToWgHnkb2ZWsDrOJXY6+NdtAsbM2qmOc4mdTfu4mMvMylLHArDOFnm5mMvM2s5FXn24mMvMuqyzwb+OEzBm1i1Vzjt2NvjXcQLGzLqj6sKvzgZ/F3OZWZWqLvzq7GofcDGXmVWn6sKvwkb+ks6U9DVJ35N0t6Q3JO07JR2StCf52lRUH8zM6qrqecci0z6PAX8ZEeuAFwLXSlqX3Pe+iDg/+fpSgX2Yw0VdZlYXVc87Fpb2iYjDwOHk9s8l3QOcXtTxFuOiLjOrk6oLv0op8pK0FvgGsB54E/Bq4GFggt6ng5/0ec42YBvAmjVrNh48eDBXH1zUZWZdU2mRl6SnAZ8D3hgRDwMfBp4DnE/vk8F7+j0vIm6MiPGIGB8b658bG4aLuszMnlBo8Jf0ZHqB/5MR8XmAiDgaEb+KiMeBjwEXFtmHWVVPrpiZLabMeckiV/sI+DhwT0S8N9W+OvWwy4F9RfUhrerJFTOzhZRd9FXkyP93gFcCvz9vWecNkr4raS/we8BfFNiH41zUZWZ1VnbRV5Grff6bviUMlLa0cz4XdZlZXZVd9NXZyzuYmdVJ2fOSDv5mZjVQ9rxk66/ts3vXR5MiihmmNcbUhmp3zzEz66fsoq9W7+Q1p6o3cSyWe6LXzDqhszt5VX3JVDOzump18HdVr5lZf60O/q7qNbMmKqPSt9XB31W9ZtY0ZVX6tjr4u6rXzJqmrLnK1i/1dFWvmTVJWZW+rR75m5k1TVlzlQ7+ZmY1UtZcZavTPq7uNbOmKavSt7UVvq7uNbOu62SFr6t7zcwGa23wd3WvmdlglQV/SZdKulfSfknXjfr3u7rXzJqs6CrfSoK/pGXAh4DLgHXAlZLWjfIYru41s6Yqo8q3qpH/hcD+iLgvIh4FbgW2jPIAru41s6YqY86yqqWepwNTqZ8fAH47/QBJ24BtAGvWrMl0EFf3mlkTlVHlW9sJ34i4MSLGI2J8bKx//t7MrI3KmLOsKvgfAs5M/XxG0mZm1nllzFlWlfbZDZwr6Sx6QX8r8IqRH8QVvmbWQGVU+VZW4StpE/B+YBlwU0S8Y9BjXeFrZja8Wlb4RsSXIuK8iHjOQoE/K1f4mlmTFb3Ov7UXdivrmthmZqM2J3ORrPNfMbmD3TCyzEVtV/vk5QpfM2uqMjIXrQ3+rvA1s6Yq49pkrQ3+rvA1s6YqI3PR2pw/uMLXzJppasN2VvRZrTi1cfvI4lhrR/5mZk1VRuaitTt5mZl13ULr/Fud9gFX+ZqZ9dPq4F/GWlkzsyZqdc7fVb5m1iRFV/WmtXrk7ypfM2uKsjMVrR75u8rXzJqi7ExFq4O/q3zNrCnKqOpNa3Xwd5WvmTVF2ZmKVuf8wVW+ZtYMZVT1prV65G9m1hRlZyoKqfCV9G7gj4BHgQPAayLip5LWAvcA9yYPvTMiXrfY7xtlha+LvsysK6rYyesrwPqIeD7wA+AtqfsORMT5ydeigX+UZpdSrWKGk5KlVOsndxS6ltbMrI4KCf4R8eWIeCz58U7gjCKOMywXfZmZ9ZSR878auD3181mSvi3p65JePOhJkrZJmpA0MTPTfwnUsMpeSmVmtpgyq3rTMgd/SXdI2tfna0vqMdcDjwGfTJoOA2si4gXAm4BPSXpGv98fETdGxHhEjI+N9V8CNSwXfZlZnVSZis4c/CPikohY3+frCwCSXg38IfBnkcwqR8QjEfHj5PYkvcng83L/K5bIRV9mVidVpqILSftIuhR4M7A5In6Rah+TtCy5fTZwLnBfEX3ox0VfZlYnVaaiiyry+iBwMvAVSfDEks6LgL+T9EvgceB1EfFQQX3oy0VfZlYX0xpjFSe+AUxrZeGxqZDgHxHnDGj/HPC5Io5pZtY0ZVf1prnC18ysIlWmoju/h68rfs2srTq9h+9CvM2jmXVVp9M+rvg1szJVVdDVT6dH/t7m0czKUrdMQ6dH/q74NbOy1C3T0Ong74pfMytL3a4t1ung74pfMytL3TINnc75gyt+zawcVRZ09dPpkb+ZWVnqlmnofJHXIC7+MrOmc5HXkOq2JMvMbNSc9umjbkuyzKxZ6lTMNYhH/n24+MvMsmpK5sAj/z7qtiTLzJqjKZmDwoK/pJ2SDknak3xtSt33Fkn7Jd0r6aVF9SErF3+ZWVZ1K+YapOi0z/si4h/SDZLWAVuB5wHPBu6QdF5E/KrgvizZBZuvYTckq30eZFormdro1T5mtrgqd+caRhU5/y3ArRHxCPAjSfuBC4H/qaAvA7n4y8yyqFsx1yBF5/xfL2mvpJsknZK0nQ5MpR7zQNI2h6RtkiYkTczM9P8YVaUmzOabWfnqVsw1SK4iL0l30H9QfD1wJ/AgEMDfA6sj4mpJH6S3ofu/Jr/j48DtEfHZQceposhrIXNm8xPHYnkt/8Bm1l2FFXlFxCVL7MDHgC8mPx4CzkzdfUbS1hgLzuY7+JtZAxS52md16sfLgX3J7V3AVkknSzoLOBf4VlH9KEJTZvPNrBhtSPsWOeF7g6Tz6aV97geuAYiIuyV9Bvge8BhwbZ1W+ixFU2bzzWz0mlLEtZjCRv4R8cqI+M2IeH5EbI6Iw6n73hERz4mI50bE7UX1oSiuAzDrrqYUcS3GFb4ZNGU238xGry1pX1/bJyPXAZh1U1vSvh75m5kNoS1pX4/8C+LNYMzaqS2Xf/FOXgVwEZhZc7Vp4LZQkZfTPgVoy2oAs66ZHbitYoaTkmWc6yd3NHId/2Ic/AvQltUAZl3TpYGbg38BvBmMWTN1aeDm4F+AtqwGMOuaLg3cHPwL4CIws2bq0sDNq30q1qaVBWZt8MT/yWQZZ4P/Ty602sfBv0JeEmpmRfJSz5rq0soCs6q14TLMo+QK3wqdGjOgfu3tW1lgVqW2XIZ5lDzyr1CXVhaYVcmfsk/k4F+hLq0sMKtSl9bvL1UhwV/SpyXtSb7ul7QnaV8r6Vjqvo8Ucfym8JJQs3L4U/aJCsn5R8Sfzt6W9B7gZ6m7D0TE+UUct4mG2RfAy0LNspnasJ0VfVbWTW3c3qhr8I9SoRO+kgRcAfx+kcfpAk9Ymc01zGCoLZdhHqVC1/lLugh47+w6U0lrgbuBHwAPAzsi4r8GPHcbsA1gzZo1Gw8ePFhYP5vgyM5z+u4edIQxVu3cX0GPzKrjGpmlKWSdv6Q7JO3r87Ul9bArgVtSPx8G1kTEC4A3AZ+S9Ix+vz8iboyI8YgYHxvrn6/rEk9YmT3Bq3fyy5z2iYhLFrpf0pOAlwMbU895BHgkuT0p6QBwHtC+8t0Ra8u+oWaj4BqZ/Ipc6nkJ8P2IeGC2QdKYpGXJ7bOBc4H7CuxDa3hZqNkTvHonvyKD/1bmpnwALgL2Jks/Pwu8LiIeKrAPrZFlWajL2a2tPBjKzxd2aylPiFnTDLuUuU1X3yyKr+rZQV4dZE3iwUoxfFXPDvLqIGsSr94pn4N/S3lCzJrEg5XyOfi3VNYJMU8SWxU8WCmfg39LZV0dtH5yB6uY4aTkEhLrJ3f4DcCGNuwgwqt3yucJXzvOk8Q2Clknb716Z/QWmvD1Tl52nKsmbRQWnLxdIJgPc4Vby89pHzvOeVcbBU/eNoODvx2XJ+/qiWKb5UFEMzj423FZdxbzRHF7ZXlT9+RtM3jC13LzRHE75am69eRtPXjC1wqVZ6LYW1PWV9aJW/DkbRM47WO5Zc3xOl1UnizpG0/ctpuDv+WWNcfr67mUI+ubrCdu283B33LLOlGcd2TZtRVGWf+9Wd9kPXHbbs7520hkyfHm2ZpyzmRkMppdMbmD3bN9WcLzmzTXkOffm3VO5oLN17Ab5k7cbqz3ebKlyzXyl/Qnku6W9Lik8Xn3vUXSfkn3Snppqv3SpG2/pOvyHN+aLc/IMk/KKO9cQ55PHGWP3iFf+uaCzdewaud+Tvrbn7Jq534H/hbJm/bZR2+T9m+kGyWto7eN4/OAS4F/krQs2b/3Q8BlwDrgyuSx1kFZ00WQL2VU1RtHnufm+fc6fWP95Er7RMQ9ANIJnym3ALdGxCPAjyTtBy5M7tsfEfclz7s1eez38vTDmivrksA8KaM8S1PzLH/M89w8/16nb6yfoiZ8TwemUj8/kLQNaj+BpG2SJiRNzMz0H/VYd+UZzeZJg+QZgVc5enf6xuZbNPhLukPSvj5fW4rsWETcGBHjETE+Ntb/P6t1V56UUVVvHHlz71n/vWb9LJr2iYhLMvzeQ8CZqZ/PSNpYoN1sKFlTRnnSIFMbtrOizyUPpjZuX/T4eZ47229XzdqoFLXUcxfwKUnvBZ4NnAt8i16m9VxJZ9EL+luBVxTUB7OBqnjjcO7d6iTXhd0kXQ78IzAG/BTYExEvTe67HrgaeAx4Y0TcnrRvAt4PLANuioh3LHYcX9jNzGx4C13YzVf1NDNrqYWCvy/vYGbWQQ7+ZmYd5OBvZtZBDv5mZh3UiAlfSTPAwRy/YiVQxx0o3K/huF/Dcb+G08Z+/XpE9K0ubETwz0vSxKAZ7yq5X8Nxv4bjfg2na/1y2sfMrIMc/M3MOqgrwf/GqjswgPs1HPdrOO7XcDrVr07k/M3MbK6ujPzNzCzFwd/MrINaFfybsKG8pE9L2pN83S9pT9K+VtKx1H0fKbov8/q1U9Kh1PE3pe7re+5K6te7JX1f0l5Jt0l6ZtJe6flK+lDqa2eBfpwp6WuSvpe8/t+QtA/8m5bYt/slfTc5/kTS9ixJX5H0w+T7KSX36bmpc7JH0sOS3ljV+ZJ0k6RpSftSbX3PkXo+kLzm9krakPnAEdGaL+A3gOcC/wmMp9rXAd8BTgbOAg7Qu6T0suT22cDy5DHrSuzve4C3JrfXAvsqPHc7gb/q09733JXYr5cAT0puvwt4V03OV6WvnXl9WQ1sSG4/HfhB8nfr+zctuW/3Ayvntd0AXJfcvm72b1rh3/EI8OtVnS/gImBD+vU86BwBm4Db6e2N8kLgm1mP26qRf0TcExH39rnr+IbyEfEjYHZD+QtJNpSPiEeB2Q3lC6fervdXALeUcbwcBp27UkTElyPiseTHO+nt/lYHlb125ouIwxFxV3L758A9DNgbuya2ADcnt28GXlZdV7gYOBARea4gkEtEfAN4aF7zoHO0Bfjn6LkTeKak1VmO26rgv4DcG8oX4MXA0Yj4YartLEnflvR1SS8uqR9pr08+St6U+ihe5Tma72p6o55ZVZ6vOp2X4yStBV4AfDNp6vc3LVMAX5Y0KWlb0nZaRBxObh8BTqugX7O2MncAVvX5mjXoHI3sdde44K+KNpQfxhL7eCVzX3SHgTUR8QLgTfS2wXxGif36MPAc4PykL+8Z5bFz9Gv2MdfT2xXuk0lT4eeraSQ9DfgcvZ3zHqbCv2nKiyJiA3AZcK2ki9J3Ri+XUcl6c0nLgc3AvyVNdThfJyjqHBW1h29hogEbyi/WR0lPAl4ObEw95xHgkeT2pKQDwHnAyLYwW+q5k/Qx4IvJjwudu1L6JenVwB8CFyf/EUo5X4so/LwMQ9KT6QX+T0bE5wEi4mjq/vTftDQRcSj5Pi3pNnrpsqOSVkfE4SRlMV12vxKXAXfNnqc6nK+UQedoZK+7xo38M9oFbJV0snqbx89uKL+bZEP5ZBSwNXls0S4Bvh8RD8w2SBqTtCy5fXbSx/tK6Mvs8dN5w8uB2ZUHg85dWf26FHgzsDkifpFqr/R8Ud1r5wTJ/NHHgXsi4r2p9kF/07L69VRJT5+9TW/yfh+983RV8rCrgC+U2a+UOZ++qz5f8ww6R7uAVyWrfl4I/CyVHhpO2TPbBc+aX04vB/YIcBT4j9R919NbnXEvcFmqfRO91REHgOtL6ucngNfNa/tj4G5gD3AX8Ecln7t/Ab4L7E1eYKsXO3cl9Ws/vRznnuTrI3U4X1W9dgb040X00gJ7U+dp00J/05L6dTa9VVDfSf5W1yftvwZ8FfghcAfwrArO2VOBHwMrUm2VnC96b0CHgV8m8eu1g84RvVU+H0pec98ltapx2C9f3sHMrIO6kvYxM7MUB38zsw5y8Dcz6yAHfzOzDnLwNzPrIAd/M7MOcvA3M+ug/wfaLdn4SZ4TYwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Compare fixed point set with resampled moving point set to see alignment\n", + "\n", + "fig = plt.figure()\n", + "ax = plt.axes()\n", + "\n", + "n_points = fixed_set.GetNumberOfPoints()\n", + "ax.scatter([fixed_set.GetPoint(i)[0] for i in range(0,n_points)],\n", + " [fixed_set.GetPoint(i)[1] for i in range(0,n_points)])\n", + "ax.scatter([transformed_moving_set.GetPoint(i)[0] for i in range(0,n_points)],\n", + " [transformed_moving_set.GetPoint(i)[1] for i in range(0,n_points)])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Visualize Gradient Descent\n", + "\n", + "We can use the ITK `ExhaustiveOptimizerv4` class to view how the optimizer moved along the surface defined by the transform parameters and metric." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the new optimizer\n", + "\n", + "# Create a new transform and metric for analysis\n", + "transform = TransformType.New()\n", + "transform.SetIdentity()\n", + "\n", + "metric = PointSetMetricType.New(\n", + " FixedPointSet=fixed_set,\n", + " MovingPointSet=moving_set,\n", + " MovingTransform=transform,\n", + " PointSetSigma=5,\n", + " KernelSigma=10.0,\n", + " UseAnisotropicCovariances=False,\n", + " CovarianceKNeighborhood=5,\n", + " EvaluationKNeighborhood=10,\n", + " Alpha=1.1)\n", + "metric.Initialize()\n", + "\n", + "# Create a new observer to map out the parameter surface\n", + "optimizer.RemoveAllObservers()\n", + "optimizer = ExhaustiveOptimizerType.New(\n", + " Metric=metric)\n", + "\n", + "# Use observers to collect points on the surface\n", + "param_space = dict()\n", + "def log_exhaustive_iteration():\n", + " param_space[tuple(optimizer.GetCurrentPosition())] = optimizer.GetCurrentValue()\n", + "\n", + "optimizer.AddObserver(itk.IterationEvent(), log_exhaustive_iteration)\n", + "\n", + "# Collect a moderate number of steps along each transform parameter\n", + "step_count = 25\n", + "optimizer.SetNumberOfSteps([step_count,step_count])\n", + "\n", + "# Step a reasonable distance along each transform parameter\n", + "scales = optimizer.GetScales()\n", + "scales.SetSize(2)\n", + "\n", + "scale_size = 1.0\n", + "scales.SetElement(0, scale_size)\n", + "scales.SetElement(1, scale_size)\n", + "\n", + "optimizer.SetScales(scales)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "MinimumMetricValue: -0.0648\tMaximumMetricValue: -0.0153\n", + "MinimumMetricValuePosition: [15.0, 15.0]\tMaximumMetricValuePosition: [-25.0, -25.0]\n", + "StopConditionDescription: ExhaustiveOptimizerv4: Completed sampling of parametric space of size 2\t\n" + ] + } + ], + "source": [ + "optimizer.StartOptimization()\n", + "print(f'MinimumMetricValue: {optimizer.GetMinimumMetricValue():.4f}\\t'\n", + " f'MaximumMetricValue: {optimizer.GetMaximumMetricValue():.4f}\\n'\n", + " f'MinimumMetricValuePosition: {list(optimizer.GetMinimumMetricValuePosition())}\\t'\n", + " f'MaximumMetricValuePosition: {list(optimizer.GetMaximumMetricValuePosition())}\\n'\n", + " f'StopConditionDescription: {optimizer.GetStopConditionDescription()}\\t')" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "# Reformat gradient descent data to overlay on the plot\n", + "descent_x_vals = [iteration_data[i][0] for i in range(0,len(iteration_data))]\n", + "descent_y_vals = [iteration_data[i][1] for i in range(0,len(iteration_data))]" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAADck0lEQVR4nOz9eeAny1nXi7+fqu7P8t1nPzNz9pzDOUkgBDiExLDID1GvoqhXNnG9KGLCklz2VUBBFhUFVARBQQGBC1x2VCKLgoBJWJOQ5SQ525wz+8x3+yzdVc/vj+rqrq566jPz+Z7vIfnmTsHJzLy7P9VV1dXV1a/nqaeImXEn3Ul30p10J31gJvW+LsCddCfdSXfSnfTCpTuD/J10J91Jd9IHcLozyN9Jd9KddCd9AKc7g/yddCfdSXfSB3C6M8jfSXfSnXQnfQCn4n1dgDCdPHmS77///vd1Me6kO+lOupOOVHrTm950hZlPScferwb5+++/H2984xvf18W4k+6kO+lOOlKJiJ7IHbuDa+6kO+lOupM+gNOdQf5OupPupDvpAzjdGeTvpDvpTrqTPoDTnUH+TrqT7qQ76QM4fUAM8k9eexa/9q434vLOtZ7+7PZl/I/H34QLNy/19Ct71/E/3/1mPHn92Z5+Y7KD33jP7+Dxq0/19L35BL/5xO/hHZff29On9Qy//dQf4C3PvQthDKDK1HjzM2/F7z/7dli2rW7Z4veffTt+98LbUFvT6syMt116HG++8FbMTdXT33n1vXjzhbdgUk17137ixjN487Nvwe58v6c/s3MRb37urbgx3enpl/au4ncuvg1XJzd6+rXpTfzepT/Cpf2rPX17vos/uPJ2PLvXb7tJPcVbr74Dz+w+19PnZo63X38cT+1c6LVFbWs8fvM9eHLn6Z5u2eKJ7SfxxPaTSRs9s/sMnth5EiZqo4v7F/HkzpOobd279tXpVTy9+xTmZt7Tb85u4Jm9ZzAz/bbbrXbw7P4FTOpJVLd9XJw8h/2636YzM8PlySXsVXs9vbIVrk2vYK/a7enG1rgxu4bdqn8PLFvcnN/AbrXTawtmxk61Lep71S72BH1S72Ov2u21HTNjaibYr/s6AMzNDPv1Hgybnl7ZOSb1fqLXtsbUTJK2NmxE3bLBzExR2yrSragzW8zMDJWdJ3WbmxnmOd3MUt3KemXnmJkpOGoLp0+SNqrsHFMzEdqiwqTeT+pc2wr79R6qqG61rbFX72JuZv22szX2qp2kPxo22K22MYn63WGl9yvvmmXT/nyCT/33X4Q3vP03MSwGmNZz/N1X/RX8s7/0Rfj7P/r1+PHf+28YFgPM6jn+8od8Ar7n078WX/Hz347v+62fwLAYYF7P8Scf+kj8x7/+zfjO//FD+Oe/+h8w0CUqW+NDzz2KH/4b34r/9w/fgK/9r/8KhdKorcGDJ+7GD/21b8X/fvoP8IU/+60gAqy1OL12Av/pM74Fz9x8Dq/9qX8EYw2YGWvDVXzvX/3HsGzx93/yH2J/PgWIMNAFvvMvfjVOrm7hc37qa3FtcgOK3Dv3m/70F+KlZx7GP/jZr8WF7YtQSsNaiy/7mL+Hj3/wlfiCX/wGvOPqe1GqApWt8Tkf8en4lJf+WXzZf/+n+N2Lb0OhSlS2wqe++P/A53z4p+Prf/1f4deffjNKXaIyFT7xgVfjiz7qs/Adb/4B/Jf3/k8MdIHK1HjF2Zfhq175GvzgH/00fvLx/4qBcm3xkuMP4Wte+bn4pSf+B37o7T+FQhWobY1718/jKz/q8/D7l9+C//C2H4UCwcLixOgYvvgjXosLe8/he9/6H8HMYDDWylV83od+Nqb1BN/31u9H1Tw0A1Xis176tzEqhviBP/oBTOopiAgKhE97+NNwanwKP/iOH8ROtQMFBQbjk+77JDyweT9+/PEfx7XpVShSsGzxsec+Di898cH4uSd+Bhf3L0KTgmGDDz/1GD7i5Efgl599Ay7sPQ1FGpYNPmjzUXzU6T+B37r063hi973N+Rb3rN6LV575aLzl+u/h3TvvgoKGhcGp0Rm88vSr8Z7tx/Gu7T9y1wVjs9zEY6dehYuT5/DO7beBQGAwVvQKXn7iI7Fd3cQ7b74VAMAABmqADz7+YajsHO+8+TYwGGBAK40P2nwpCMC7tt8O2ww4RAoPrn8QBnqA9+y8E8YPOATcvXI/Vss1PLn7bsztHASAwbhrfB7r5RYu7D+BmZm2ZTo+PIVjwxN4bvIMpmbS6hvlFo4PTuHK7CImZg9o9NViDScGp3Gjuoq9erfJHxjpMY4PTmG33sZufbPJBxioIU4MT2O/3sV2daM9v6QSx0enMa0n2K6ut8+xJo3jw9OoucL12ZVWV6RwfHgKli2uz664NgJAIBwbnoIiwtXZpd4gfmx4CgUVuDK72LYdAGyWxzHSY1yePoeaq56+Uq7h0uRZzO3UNWjTFhvlMVyaPouJ2WvbaL3cxInBaVycXsBudbPVV4t13DU+j0uz53B9dqW7/8Ua7lm9H5enF3Fp8myrj/UK7l9/CNdn1/D03nsBIoAZQz3CQxuPYqCHz3d4bBO9P0WhfOyxx3gZF8q/85++Cv/5zb+Aad3N4FbKET7moQ/H/3rv7/Vmv+NyhFc/+HL89pN/gP1AHxYDvOK+D8HvX3h7Tx/oAh989mE8fu2pXj5aadx/7Bwu7l3FtO7e1ATCqbXjmNYTTOr+G3x9sAqlFfbm/ZnjqBxiY7yK65ObPX1YDHB24xSe3b3cm22MiiEePHEP3n39KdRBBx4XQ7zk9EN429XH24HTn//yM4/iD6++E7NgljvUA3zYXS/GW6719YEq8SGnH8G7br4H00AvqMAjx+/Hhd0LvfM1Kdyzdg435tcwD2YzBMKJ0THMeZrMclaLFZQKmNv+rHuoh1gbjDCNZjmlKrE13MRONCsuVYkz41O4PrvePvyurCXOrZ7FlellWNhAL3DP+t24Mr3cm6kVVODc6nlcjXRNGmfGZ3Gzut7TFRROjk5hr97p6QTC5mALMzvrDS4EwlivwLLplce19xCaVKJr0tCq7OXjrz3QQxjuzygJCqNihCpqU3/tuZ2lerGKmZ0BQdv58yV9pMaoeN5rawAYqjFqQS9pAMN1ohdUNn23ryu4l2Wsu6sjyccf4ajtfHtYmORcguoN8F5XpMU2KmiAmZ30rk0gDGiIidlL9KEaYc/s9crU6vVeUtahGmG/3kvu/0iP8ZKtDwURJXXLJSJ6EzM/Jh47qoP8rJpj80tfhVk9T44NikLsFGWhBRUodeFe4EI+eZ2Sz71RMQApSj7rRsUQSqkeinH5lCgLneiF0ijLood0ANdhykEBK9yzQVnKdS7lj7XhoBDbYlCWkPrWSBdQol5CK0quPSpKcQAb6RKlUskn8VAPMCpK1NEAVpDGsBgm5ysojIuhWOdRIc+CVsqRqI8ys6ahHooP2kCV7VdXmEoqoVSqF1SI5xdUoFBpX1XQ0EoLuhLPB4BCyfe5IFnXmfM1aVFXUGJbECg7GJH08IjaQVPaDlLbAGi/JiU97qOLdADdV1SU4r7bnW/Ea9c2fQkCrq0f2fpgrBSrYn5SWjTIH1kmP6mmySDrU+5GS4PjwpTtj+kA768bD/AAYNgmAzbg2L1UJrOwXrlCLVe3fFPIB3JXJcrNspB5eJAM2IBjt1KbWjDkichy5VyUlp/oZK6y1GAHEKUvx6ZE8oCUOTuX/0HSYU35DrNMcnr/mZzeKuXGo3xKJ4rPJx3ZQX5zvI67t86Ix7ZG60vp43IEJTygBWkUKp3ZMBgDXSa6ZVkHGFqY4ZWqgNRZC6XFQbhQGlqYESpSGGZmo6vFWNbLFVEfa3m2q0lDCQ8uszz7qy2jEHTLVtRz7y+V6aJEKjtTzM1Gc7PaUg1EXWfOl/oKAIBZLJNlK+rGGlFnyPn4o8ul5V480j0GXHvL+gs9mL/wSWoLZxOS66wy/UtjOT2XP8NidYlZ/K3SkR3kiQgPnLhbPHbXxklRP7N+IqtLs7ljo00Ym84uV4oRKpO+aTUp8Q3MzLBCPrWtxYeqtgbjIh14jLXYGKwlumWL4+ONRAeAE+MtUT8+2hT1YwNZ3xysy2hEj8QZeKG0qBNRw177ycKCWBgI2YgvWstWfLE5Y5f8Ast9/ub0sZZfkEOS9UINxDbSQvkBiAgH8Pw4Nwinv2FwdiApM74VBUmTEcfSxfMz+eiMnh0gDzTDv33s49l7ohOJbURE2brl2mJA8oRqoPLYb5nzCyqzL5KDpCM7yE+rGX793b8jHnvHZTmMwzuvPinqT954NpmRKCJcmV4XH9qpmWVnMCJ7VVqcyWulxU98BYVZZLD06cZ8W9Qv7l9JsiICnt57Ttb3L6bXJsKl6WVodI8QAdAAblQ3QdQfehQRKp5DERJdK6Ag96iTz4cIA6VRKg3d5OX1YVFCK9VrPwJhoEqAXLtSoBekUXOV6JoUZnaKQul2qFQglKrAvtmDJt0OoNTM+idmDwUVPb2gAlMzgUZ0PjQqnkNH57uvB3b8vXmsvK6VY++q+fogEApVoFDuYXYPtGsNRbp9yH1ZfT6loCuolvmHZVVQbvBSFJTVn1+6siKqMwowrNgWTCy2BYOhoh6j4Pq1Qv9ry3F9JeoFFb174/qRq5NvtbBfqOD/Qp3g+lDX8zrdD/RhHRQUoJD2i8abSmojgxoFlc3LxOXiDMqVeH7FddKmmnSTT9pfFBEm5vDcKY+sC+W0cn6x0mDLyOiZ89Gcr6NjhGj0Co4wbJoXOZ4e644r22RMdd4TRXI+U+MqJpZVNnRxc770E1HnJn/hgHsYkqvCnR61Ebn/iREGNfVI8BK5WbtS8lxPEfVmMUzOtVARAYFOwfl9vRnkSEHpzIwumiX5fiHhHHEGSD6fIjq3yz/8DTX5xPm7+9IMQVE7EaEdsG6l9+ostKqiNB9//XgW65+dtC1YbgssniXnUEhOT+4N/OtP6MCNFl7bT8mIVK9E7ks9rQM3/xe3ETN3/SIeJlm+n62enC+3qW3Giri/ALLd6qDpyM7kt1Y2cP+J8+Kx4ysycsjpa4NV2VtCFSIqABqPnChZthkmD5Gla6VFdCEya7gvAqk8ihRWCpmlr5cyitjI6CvFivheC2ccYVqWyRtrxfNzts+80VL29nBMPvf5nUEUGSaf81jJslR2149TjslzTl/I5JdNSzL5DELK20Zy+S+X8sbJF97Amn8R5doiw95znkw5Jp9pa4cc7zB5AMC5DTF8Mk6sbIl6bpA/Nl5PVsUBwMZoXWTyIz3orcb0SZESvWjYsjiYG2vFG23YYCi8LIy1WCtTHsxssTFMWT0AbA1lY/Nm5vz1war4WK2WK6JP8lANRfauiW7hJRRr8oBnWX4pWLYYiIMzZ10iRxmjco69j5R8fpalZgzpuYc8+wJbMETmWfdyRr+sa2WWsef0nFFRPv8gdZaTh32Smhu0ha8NyF8hwCK7xZJ2jkx/KZGbEOrsfT5IOrKD/GQ+xW8+8XvisXdl2Pu7rz0t6k9vXxQHtqsTmcnPbCXOPBUp0XhY6AI6mnkSEQrtWH2sEylUXKckhRjb1Z6gIwlL4NMz+5cExEK4MLmU+L0rAi5Pr0IJ7H2n2kkeKwWg5jkInOiOr/fRkiJq/ec9Sye4dvO+5zFjLagAKeoxU0XKMU02CUtVDZPXAquf2gkK6li9Z/sTuy/qMztFQZ4Jd7rhStSdfSFl/krpHj9uWbqvc8DqFXV8OtY1FY7hB0zb5aqhlAqu7a6gqQBJvJk0mNBj9f58JpmxM9nAdtDopAFikbGjyYeifPw9js+P6xXrYb9Q0f8hqEOfyXf3n0BQinp18Kyeyb+sYt1GbeS+Eg2l/c4zdne+6p1f81y0Z+TyUaQwNf2Fk88nHVkmPzdVy6GltByTl8/PMXlq/ieeMfgFFInuPdyjQc8ywzZl6l2/4dDLMPZbtUVyjD0DTesmnd/qSd4ZvfkfhT73d2SUE5bO4LZNezMrCsva6X4Ij3XK6GFbJBipZaYZ3hzrtJzetV3qltfyaaEB3WDZzwcNP+6f6uusRCy4qC0SxODbOsfeM/ksNXsWznf9RciHue0/8oy7f4U29IHA3qU6+EVSyfmhnmujpB/dQo/LnzsfSFY7P590ZGfym+N1PHTyHvHYqdVjS+kbwzXR93mgyqz7WynwWgsrsnq28stFKyW6bhYke91o0uJ1FSnR750AbAoul4BziZTSarkqvddQqgFk9ze5k9YZzFJbIzN5sTSLPuPdF4+UlvWHz+GXUsmf0zmWiiWZfJbVL7uwbeHR5fDIMl5jC/PJ3jc5Lb9g6PDSC87kc+stMvkwGOM7TN6lE6tbor6+gDdLaW2wIuKX1XIs+reXVIjWbwKJ51u24mBurBUfBZNZMMSwGAtL9pkZqxKrb+ogpdWMoXYlw61XipH4GA6olP3hoUQ7hy+vpMnGScFDxx3BIPPw5AbtQcbvOTeYF1mf8cz5qhDvpyI5nMaihUT5wVOYqUNm9Yt4c475Zu0HvJxBNvdyOayXTi6vRS+dZQdzvWQbFVm2n+l3Cwb/wzO8H2Fcsz+f4Lef/APx2Huuy+z9vTcuiPqFncsiFrk2uymeX6FGQxc6ly04bxlpYCvLEkowRGoiF+vEpC8MA2npP2MnCncLuDf/lel1EFFvACUiPLt/WazDc5MrvfL7OlydXWu7lz+mAOzWO1AEMLprKCIYqluiGOqlpmbwQWuuVQCGuoSn8d4YrUAYeE8my+1Lg0g5X3siaKhWV+Swh4VB0USN9DqBMOcZNPUXY7WsHroXbkFBYWpnbiBm27nskcKc+7rPx8KI+QDua4vR1cH5yfuHltswFp69+6iaYYTFsvHBdub6TldQTR+jXhs53uzsQZ3eDGsEaBbKSgzF3T3w1JmJobgrUzs4kvt7HJTLnxPrRASwq3NXh8b+xNyrs4+L4/8d1jnM29fBtRna/u4dAqhpS9XqXT6hx1KiE6C4i7PUnkc200ZWbFNDTb+I+p3hCpqKRvf3v7EpNWsNwrbwTH6cWdS3bDqyg3xtTYjrkiQz+VvNDKI8IMcW8Q9DMjdg91gmfu/MYEqv7Wy0LBpSkS2rwNf9EYm9Y3EbiVcQdWp8tNMM2/PjC7WMPUxNW4Cgw6tQ56PdY+xNnrH/fMjkiyWZfDoT83aRDD++zXwW+727ISK0gYTn90vDzQKzmL03htuIvR+szl5PGXiuTItnw7Iuza2z57fThbR/AfIMWub+Qh24u0aurFkmH7cFZ9YMLM3q5bZupgNJGQ+ajiyu2Rit4dEz94vHzqxlwhesHRf1Y+N1kcmP9CCDCmQfagMrs/oMk1cgMUCZJi0zfNKiH74Cif7wBODYUA53sJXR18s1ceAfZpi8tGAHAIzlLJOXzrdsM3A5/9mf48Q5VzbZ5RLZ2N15Jp9Zcp7Byu4LR2byop7JKBthcSHPXhaPHJaffE5//0u5QT/vrrocxsnFUsqxegZjpA9nFg8c4UEecCxdSsNMqNlRmYstMRCZ/LDIxCIhOS4LWI4gycgweZaZvF8Jl+iwGAgvEYaLfy/rmTpnBrBBRi91KQ4lBWlxJZPHE2IS5NxCoiyrR36wlV60QP5hyy4+WzYuC3VL6cMkzx39MVnLsffcDFjmyvkr5zj0sudnB39RXXxkmfMpc2QR25fbdNHdWc7wmlsbsDi+jzD5o8Nl8kd2kN+b7eNNT71VPPZkhr0/deM5Ub+4d1UczG/OdzOhgI28ACjniVMUYpxxRUrUQ3bbU5mxV6f+swzO2g8uTa8KGIVwZXqt18U8P78xv9H+G8GxSb3fukt25SdY1AC415EUgIFWTm9QDsHZIEbFoPWhD/NpY9dE/tDepz71k6bGOyX1h57bWc/j2j3ehJmdCZ/GhJmZRue7r6OK573IKQ4MqMZnvFvf0OIZz8F7ZVItk4/r4OOZJLpycenjmDaKtNgWnW94cL5vJeoPfB1j9wM3BXrnxx+3kcOT/R7jzw/923193Ey4P4AqUNMfVNJGbmCT9O5Yqqve16UvRbyYyJcvHujbnim0kfuLXGdLafRQAsGSCYymvkWd/7zURrbpR93kwN1nIoWZveMnv3xseORj2gAytxY580GuwZ7jxRfIxKjx8Fu6Ni9i8nKZyI3OcZFE9u5myL4zpvnkmDyANN6NcH5r0KN+rCAOftOL++LLEuvodMrqy/HjPKtdcH5y63zbpbweQFK3W/Hm7BqATJ0PK4bM4emuhLHuVLmN8r71fY94r+fuT68tOOzzcT4yw2dJR2fbS+rcMnaZyafn+zqn61UOcy+n5z2TJ6J7iOiXieitRPQWIvqCRj9ORP+NiN7Z/Hns+Re3S+ujVbz07IvEY2fX5XAHd63LIYhPrGyJnXSlGGfjt+fC30q6FBoBaN7+IpPvz2hCfSSEIFYg0e+dABwfbonXPj6S9c3BhtgWQz0SdT+7jFPN6aIdwO2qI+Vj2Igde5H/fI7J52LUZF0rl2TyOeyTS7lNSXJBqHILYQ7C5PP4JY81pHRYTH55XPPCp/wL6rCY/HIYx4XmkF2fD5IOA9fUAL6QmV8C4JUAXktELwHwZQDewMwPA3hD8+9DTdLCIwJQ6HxQMamL5ZisogU794gjkqz7CHhxssxiecRZP/qz7EQXsY8L9yulnEFZkxLrvGwccJKrDMEnqTmf5C8Rlh9CHzlQSqKHCfL3Oe9Lnos5s4jJC+dnZrrItcUh6f7oUroQlmNRLvkh5IUezJd/ecmOAx4mLpHXknXOvyyWe6EeND3v3Jj5WWZ+c/P3HQBvA3AewCcD+P7mtO8H8Jee77XCtDeb4PefeUevmfzfn74ps/cL2xcRkgtCE69l/5o4sO3V+zIbB4sGVq3TfTmBBUxeyUzeBn6zPZ1tskm4L8+N2bagA5dn1xMdAK5Or/e6t//7jflN8fxptKGxTy5uTcreB4WGY/XdB7sCYajLxHNBQaHUBZTqexUREZRWySDpBzWG7cXZcW6WhJqrHkv3fNr5zwc6uZg2cztvWKqPOaJalhrHkPFx412c80BvoEE84PovjpC9+79L8Xp6MW0S+4RuWH0Qu6axDeTaqJuR+vvgvhIXs/our/Dvqqd3DLzP0gmtP3wygHbllNsozSeM7dMrK1HjCRTnHzL8rs7+2mFOXfnTtnB/4fS66OwT/TZythqpbkymF8rY19eSbftmUGMQEWa2v6H980mHyuSJ6H4AHwbgtwCcYeZnm0PPARD36iOizwbw2QBw77333va1wlgU4hwpM5HoOHGqy+fn8smXTWTvC1m6wCxzfvIpBs7mc7t6UtzsJ3j+uq5NZfaqe52+O5b3V899Pnd6azeIzm8fvEWs9jbyT3Shv+TOD8vRlupW7H1J/Xba6MA6pDqEz5p8n3OMPX/+7eQTPucUqKE96zbqxu5/suwdaVvkmTy316X+D+S6BXrvKzPUM1+Nh5EO7buAiNYA/DiA1zFzb1rJOV7hjn03Mz/GzI+dOiWzdCmtDVfwsvMfJB67e1Pe+/XcxmlRP7V6XOzsuTjzWuls+FuRQ5uM5zPLBmQFeYavSIlhDQgk+r0TgJOjY9KVcSKjb2W2/xvrFbGNFMlhUY1lUa9yTN5a2dgku5hn49gD+Rg1w+z2bHIoh9xy9Nx1c5/lufAOh8XkF6fcCzuHCpacKGRf/bnSLHf+C50WveiWjV2zrJ5j9QAwzIS5Pkg6lEGeiEq4Af4HmfknGvkiEZ1tjp8FcOkwrhWmnEEzp4u+7QCsEAMeAKwQbsCdn+7yBOQt4uES6+gI5HdfhvlndfniLvclF9YssBPI+dx+zPhMMdtEmd/kWCrl8sosquIcb87kkyvP4qVH+dm9fHZ+Jn27uj8qptx9W/p98cIOzvm6HU5b5Gb8B0nLvmyzNcv2u8Nt6+eNa8i13vcCeBsz//Pg0E8D+FsAvqn586ee77XCtDvbx1uefZd47NmdTLyWnSviA3d1clO8E/tmJt5QJtmPvVBKnIUVOrO4gXxsi/5AaTODv2XGzMzT8oBxs9pN8wdwfXZD7DLXM371O3UXN751awQw5ylyj4kVBvpS63brvjANirKfeZMK7dg3R+2nyc+pqFcerXyclb5OINSoxZLWqJrZWfgLuH1qm7+HRwxSf2g0Zy3CH4TOYE/oYoTHL8pCFSClmhG3u7L38477nue3se7YNYkDN2WcB5wupUVDjNSHFw3C8iI5IkrawrddqrtfSZOMsJ3jMkp1I5Jebtzi0b7a7yP9lAk5QCzuJyGFHwfy/YvgmPxhedgcBpN/NYC/AeAPiOh3G+0r4Ab3HyWizwLwBIBPPYRrtelgb7vlfnOrN7bIRaWO5xlqLC8qD2Wen6ynxiL2Lmaff0CXZfULDBRL/Ub0pAnYe3gqQ9Rv77qLWPoy+eT125nRU09/PjP6zHqLpfORy3lL/QCz7du/hlw3P+inHF2+NmfzcUeT6wasPk4s5O+vcfv9Wi7nrfSDpMPwrvmfzEzM/DJmfnnz388z81Vm/gRmfpiZ/xQzXzuMAvu0Ohzjw+5+VDx279ZZUc+x+jNrJ8RG3Ryui4zd7dwje8VIN9kYI3LZHMZRkANNKSIxdDCBcHwks/RTIzlez6lxJo5Pxq9+tVgT26jbOaefrJUfEGl7RK/nPIpyLyOZaVLWHz7HOXNMvlyaycvlXNpPfkkEtji9rwb4BeU5xIU+zzctNkovy+RzrrXLxZkH8ms6DpKObFgDAJhUqTshAMzrFGkAQGUq8YbOTSX2ynk9lwcqznzGZ3aeygWoWvRAyfYDQmVr8fxKGDwJtGBQTfMBgMrk9Er+cOVaZObevTE9kEEdkD1G3Es286kuxrHPGzqlWP9AfvX0soOtn1vKKVPnpWfwy87wll9AtUQ2Sydq/0c8ctv6Ig8w8ctuQdvlF9wtac/K5ZQ1gCwPzA6SjmxYg53pHt528T3isYt78kfDpYx+Y7oj3qCZrTK+4W6gj5NSWvSfL4pCxDUgP4j187KZgGYMKw7yDMZuJs78jSr1nweAG/MdsSvt1bviA1TzHEqlxk5XrTTQWqmK5nO3nwa6FNu0ULr5OurTUK36PL69rmrW7UQMlEhlWWeNKtF83aSU21w897nuPb3jsjomnzJzt8mI8HUEmeF7n/g4n5ZnJ6xe5t8uL1n3v0yl5WbsecbuBluprEA6UPrBOdUh180XSKybPGznWD0IonFUsjUtSty4UKdXsOJ99n7y709M/n2StFq8BjPHy3M/EmfgC1gnQRyzm00l+gM9ESWbWgNdFLqU1ZJ4AWefW26f2kXlzz24UmoHXoFdLmhW4dp+w4+ojZqH+bYZO2f0g/LpQ2jTHA/O/YbbtnuB0ckh1O0g+R8or9vWM4x9aYa/JKtvLp37Yl+qjbJ9+P2Myb+v0spgjI+496XisQeOnRf1HKs/u35KbNSt0YaIHEpdQoq7nfvsN0aOWmnYQnZBJDFevQJhTdjCkEBZf/gzYzlez6mxHHP/eIbhr5Uyky/UQNRNjslzLc7kKyvrdsFLTdrLlZH3hx9lmHxuxpTzt88x+VxYg1x8eCvs/uX0w2Pyyxtdlw1BvJz+/payyIwWxZNf0h/+AH7yd5h8k7anqdsgAdiZ7Yvn7y7QpS45mU3EDlCZShz8LVuZT2cGf+lzFnAdTGLpRISpENYAACZCCGIAmFTy8uip4IoJADOThuN158t6zbKdI4c6iGW/d5UZFlwb5XSJyVPW3mAy9omcnmP42TUDuUB0Ge+KZQfCgwyc2bIusd6iObD0teV0WIO/zNiXf+nk1rAseqnmzl/u5Zxbt5NzPz1oOrK4Znu6i3defjLRGcCV/evib65Oboizwp35ntj35plZJyB7iWilRVavMwHT3OegAmI/+dxmImDxugzGXi0P5tu17D+/U+2I+n69715gUbUt11Aq1YlkRlnqQmbvunCxPxI+La8lUET+ItF1/b6d0fkgcSZMaPajFW6ngbwKl2FFv3TnJy/4qyuZgetM3TQV8vlU5Jk8JEy4iMlLfvVNm4r8W3YSWJRy/uo5f3ggfckcXJfKIw+sywydeXvFIgN7ut+tU618bXKrwqW6ze086/W1bDqyg3ypC7jn//Y75DLnAos/j+WHykVxjI2vjsmnerzpQZhPtrxLGn1yaZF9InttSfdc+Ta5onNMc8GZ4nyJSFhM0tD6hL3K1wgXId1W+W+hy9fw80iB+wptkStT3uZzC96clPAWdc72MfnbKZcWPT0vuD3gNvVbsvfo97d6lOR65f3n87/Jt98dJp9J43KEV9z7IeKxF524R9TvP3ZO1M9vnBab9MR4S2zsoR5kbr38WljM5OVuVgp7uRIIG4M1QQdOZxj72RU5Xs+ZFTlO0MmRnM96Ke8Jm2PyNvP05Nh7ZWtx9pTbCpHBKDLx3nP+8CMls/dlmXwuXn1uO8I8Msm4dB7iJs75wSKjH5qh9GikRUz+sGLX6Ex/yceTp2zfO0g6soM8AFzalV0ir+/LboPX92VEcWO6I97o7cmOyN6n9Vzc6s9aFiMfAvLb2qEFyb1NoTKpu58iwr7A5AmEHSGsAQHYm6eulQCwm9Mr2W6xV++JbTQ38lqC3OYgvryJRhn3NmRcAKFgWFgzAPepK6XKyi6UVeb8WsgfcMvRpZQbtJf98DrMWdyiGEWZA0vmv5z+PksLvtSWju+Uu/+Zl/NBFr0dbOGbnI7sIH9zsoP3Xn1GPHZtclPUb0zlwX9vPhEH2xpW9HtXGcOoVgqV5D9PSvS8sWTFsLe5m29hRSMhkzOMJomAXbMPb6Ly/yki7NbyYD4zsgGXIa9IVSoXu6YQ9UI1MW1ineTXY35FojyQuJdFrv3kBWO5hzDrQpfRc2XNbdCS27C526s1Lqu88fOyiEXCX4vy8WWU8eHi6y/KL9aW0W91BVFd5islO8bKoCtnB+gyk79G5ftM2YnHQdKRZfLDYrDwpuWCMi2TFr3hFZGAJEhk8kopaKUwN329aB7yuFi5FQCLfHcXPWvygCSfm/P0dftaij/IX1e4D7di8olht32oYvbq6yyz11zK+fBIqWXvgj1A1hez8dsvU0d3pfyXH1SXObK49Pln6PkP9Mvot57p3qad4xYM/6A2iuV+IT2fhzf/PrIz+VE5xCvvf5l47KGT8uYjDx6XWX3Of/7kihxnflQMl7Lq16YWZ+C1lZm8ZcZQ2MuVQNgcpnu5AsDpFdkf/tyKHK8nx+pz/vMbmTjzuf1Rc1/9yzJ5w0bs8AzO7sGaZ+8rGV1m+IMMS12ayWcGpFzc+IMx+UNi6S/Y4qVbp0Gxih8u/i7GxSWMixHGxSX8cPF3MSjStSGHlQ43ds1h7fFK2T52kHRkZ/LMjKeuXxSPXdq5KuqXd3P6dXHWeXOyDUUEEw0+k2qKQheJP3ZtDcqigI1YLpEcUlgrGeNoUqI/vCLC7jzFLATKoqgb0+XQ1c5st3lZ9dNutdu6LYZpWk8x0DoZmAwbDGiQoBCCPMNXpLJMXhoM3SdtytgJhCoTpmBmZDfTHMOvuRbHztrWICXbIaRFbNx+dqRlPby0HFfOGrQX7C+cczaQPU3yqz/jVH7SX4b9hI/Hf8Rfw+fiuwH4Qf1+fC6+BwqET7utnJZPzIu+2G8f+zk9x+Rl204uH242AD2s/nFkZ/I3p7vZvVxvzmTf8O2ZbGyc1FN5Bs4mw+SVuOBGK5UJIKYy3jVGfNBMZiGRhYUROhKDMcsMVPuZvSL3M+x9zrM2tED4n+94cVJK9ksvVTrwA85PXtJlD3A/s7x93ATID9tCVr/sTkyZ3R4WbSwuGt4F4z2Qj2Z4mC+FxehSIs6L9OeX7Cd8PIov+Qp81eZ3Yh/9Wfs+VvGV5//DoVxnGVyy/KYqwLLm5lz/usPkmzQuh3KUQxwik8/MwCw44xlDKJROjLKKHJM3CZMvMqs/M+9eRrb3yUNkPt2KrUrsUp75QZzhu5kfJe1HaDYxzq4lkGb+Qt0yuzwtijeUS7nBc3E88fy1b7dMuf6VS4fpcZFLJPwtd8ZhJfP6zweswTNfKiPBCxcKIBNcbvl0e/aP97V/UG7CcKC8Di2nP+Y0LAb4E/d/mHjskVP3i/rDJ2RWf18m1s2Z1RNil14pRgvQQjpbzDN5GQdYthgKe7kCwLFM3Pi7Mn7vd6/K9obzq3eJ+qlMrJutwZaoD/VIHHyY5Rdtjr3n9JqN2OEtbNaXOOf3Ps4yefn8XPyQZePMi4MzL8IB7zsmL62E7Yz9gn4ITJ5+9/egf+hHcC/SFewAcPc9z3/A7cqTri6VwggcJHZNjrFLMZYWnU8gFIfI5I/sIM/MePyq3CmeviGz+gs35W1mn9u+LHbKK/s3xAFmv5qgEB5oY2tR97PUOBWqEPej1UpjKoQpUKRwc5b6+isQrk5vJDoAXJnIdogr+5k1BtPrbgBi7v13s9oW22haTzIbedftku3wP7A8ACxCHVKoiLw/PGFaT8TrTjKIapph9Tm/+jqrZ9gr27RNwcuF+r1lWtbX+1YLtDj4b4G+RP6DT/or0N/27Z0wn6P8s38Rg498Nejx9+Ab8BVYQR+prmAP3/D03xKvsUzioMxhHbgXJDDUeYErcyYGEuT7n1tvsYjt3/GTB3BjsoNnt6+Ix3ILenYr+SGfmbns5ZL1k1eopIU4SmEuMXmQmE/NtciWTWZTEgsr+9uDMc8MPNMM25taOdBZZau2i6WPs2AkzjD5QmvR/1yrNMwwkPdMcsdkJp9zhhS/LHoPeXwsN3M+PCYv2jNy0SwP0WPl/SmZT/iTKL/kK6G/7dtBb/hlDB58FPqX3gB77iywt4fPwA/ju9Tn4O7N6yBY3L15Hd+lPgefaf/T+6C0LzyuucPkb5FWBqPsQ7U0k8/5jDeeBnFOtgn2n3BlQGbyWT/5Qtw0JLd4ptklI1OFTJ0zv1jW3znH5C0DhcTkmaEoNb4S3NaG8QtAkcoweelj2uEgsZwHeDhfeCZ/OOmPh8m/cC8S8/rPB6oa5Rd/eavVf+dvgp59DvXrPg/m9Z+PTwHwKQD0t30L1Bt+GdXP/iTm+LcvWJneHxPjDpMH4Jj8xzz44eKxF595QNQfOXm/qD9w7G5Rv2vtpPhYrZUr4hJ2BSXil4VMXkjGWozL1Hebwdm9XM+uyn7v96zJ8XpyrP6uVdmvfmt4TNRHepx5ucgz/NrKcXzyTL4WOzzDYpBh8jn2Ptayv3XOf34gxcBhYCAxec4z+dww/8fhJw/kw0VIjB2Uvkz8yy6+P90mGbfH6tVv/CaK7/+PnTH9sY9A/T3/BvOf/Qn3AgiSef3nY/6zP5Gp0+2kiLE3V037GMkTJPIOEGndFLTUctm48UVmLp1j8gp3YtcAcB3src+9Wzz23qsXRP3J68+K+jM3L4qd8vLeNXFWvT3bg1YaHP1fZWuU0mYfi5i8sc5XN/hPK439eYqWNClcm95MeDMBuLR3NckHDFzYvZScz8y4uCejrsv7sn5zdkNso0m9LzL52poFG16nrD5n5HJMvk7amkCYmZnI3qcZ9j4xe0k+DMbUyHhvLuTPYPG67v47ZBbnn4sb7lN8vlfT/26ly3nLv5IZuy9rXB7/4um1NW6T1e/toXz9F2PwcZ8IXLoMJoL56D8BetObob7tX4ptwEE+uWOLWjQuE7Nty9q/b7arW/h/zAEzj9qo0eM7kItp5HTpDtQZ3WQnAAdJR3aQvzHZweXMnq37mdjqOX1u5L1cLVhk6Tl/eKVUho3LTL7KeNfUNuM/z/IeQ5ZZtBEwgDlX4rAwzywYyi3ckGZ4wMGYvBhYLJP/IjiSS3lsspyeRWOZSbMSYrcDi3eSyn0FidfNlPIgKVfjw4Y16g2/jOHLPwrFd/wb2Fd9FOjGTVTf8g2Y/8p/RfUt39Ay+qOUloVmy24mcthM/sgO8quDMYrcYpLMU5jlXIuYvGgA5Sw3F2fySqEUyjrQhTzDzw4Ki5n8UnqW+cnn57Y2tJk2ys3OCfKgR1BiW1DWDCmnAzH5pWMaZfRcSOHsy2W5tOi1tmzK1/iQhvmbN1F+zudi+Gf+AlAUmP3yfwE2NlB9yze0aMa8/vNRfcs3QL/hVw7nmh9AKWeUP1Beh5bTH3MaFCVe/YDsJ/+S0w+Jes5//kXH5Jg2Z9dOiQ/oerkqz/ChxOiUtTGiXhk5HG9ta5HJW2acyOzlem5NZvL3ZZj8PRkmfzbjP388y+RXMrNR2YumtjbD5OWvqRyTt+CsH/s44/e+kmPyKsfkD8FPnuUhmBuEIJ0fopHkIBa8ZARNNsgzIOxsxeCFe7xKuYghKn7uFzD+0FdAf98PoPqi12P25v8F+zGvxuxnfhzmdZ/fW89nXv/5mP3Mj/e5h8+/wW/gflvk9K5U6YvVr5yWv5sk7y0Xu0ZqI2m3MAZDUyG0EaNA6vPOaFi9UGc6ZCZ/ZL1rmBl/cOGd4rF3XXlC1N999WlRf/LGs6L56LndKw1e6D+M27MdDIpB4kddcY1BUWIexYInRVDCitdSlaJBtlAFduf7wrZwGlcnN5LzCZRl7M/sZdYM7D0nzjyf27skDiXXptegdbpV2V61h1FRJGimMhXGOo1d44eGeIArqGjZdqzXQlAzTQpTO03KSSBMzEQY2Aj79Z7wcAL79Z6IYKZ2IgawmtupGJKgshVKVSbXsGzAmSBr7v+jOjA1X0jpIOadiqTBRNIdP44HYga3XDluj5xvOIsDvW1CcDAYuHoNw9d/KYof/lHYl74E8x/7QdiPfKzJ2JfPuvLEg5vPOacjGujbclKkZ/7OfWtH/5z0OfC8Ps4H6Nh7rDt/+LTn1ajFF3DNtWu/yFPMhVPJ27SWTUd2Jn99so3rEznIVn6T6txmEvJNACHD5LW4UEaRSgZ4ACBmGJMO5nNbyUyeaxFS5AMasbhgCIDI6gGINgXAMXkOur7/T9qvEwAKJcdjzzF5atwkEzsByQY12UAYGuak/+S5s+z/3x3JGfdiPTfbzTF2P/OL/28xk1/W2LgopbVeFjktzp6h/5+fxMrLXgH9Yz+B+dd8OSa//auwH/kR4rUPN73Q+d/ONRdf/1a2oKR/Ib+RzUHSkR3k14crC5h8/iGUUj62ep7J57j5QAh/q5QWt/Mb6FLk1m7xjFBO0aXLH1uOyefb4rCYvM2wdxlreD/59PwMk888z388TP5gD3Oi/3GNSULKO1wuZ/Sl5y5i+Kl/HaPP+Nuwd5/H5Ld+FdVXfxkwODzc8P/FdFizeOAID/KlLvHRGT/5D75LZvIvPv0iUX/RcTmmTY7JbwzWxJmzJtnrxjH5VK9MJQ4wtTVYLVOubNnixFhm4+czfvL3r8txeXL+82dXZFZ/fHRc1Ed6VZ6NCgueAPdlJDH5uakW+MkLWy3C5pl8ITP2VZ3ujwsAYyWz+qEaSVzEMXlBlx5Md5psmLbCZ7936vPnhOe3Mz+hnfKL4VLE5v4p+71D+GJjZoz/4qei+Bff2ctk+Nf/L6w8+FLoX/hvmP+Tr8P0198AftkHN7/P5I/0pXc7L8fky6qZOCdt4amP4A+fu5a0BsDpUts1fvKCrlGIdStQJm905obhx8Vv2u4OkwdgrcWbnnqLeHP+6OJ7RP2dl2X93Vefci6C0bEL25egVTpY3ZxsYzAoUZn+wF2ZCqXA5BX58MT9fApdoDbpIqCBLsS48ZoUrkyuJecTES7sXRLr9uTOBVF/avfZlvuG6bm9i+K07erkKgqV4oXdahejIg0fPLcVVtQwcZf0hqs6ivNRUuGgSWKHKFALRllFGjMzSR5ZAmG/SkNKMxh7Zld8yL3/fKLX+yBFCTOdmYnb4zfSKzuHVmXSfIZrKJaMbyF26l9fYvLOt9sz7T5Z9jttSfYA/9u+btP8ATAbMf/q4z8Woy/9akwB1H/pL2Dlk/4q9DveCXP/vZj8zI8DjzzsC9nLP25WH8c+fUe6esWpx957A71zJ44HaAsGGEndbItTonzYQzxK2trZG9I2yrWp79Mpq2/CpkSVds4GNulHhg2MraEFT72DpCM7yN+Y7GB7KseHz7H3WWbmLHF3wA2eoq+31skA786XmTwzoxbOn5saI53OCitrMCp0sr2gYdusnktnBbkFNzlWX9tafNgs5Lg5RPKm44VWYoz7QmkxMJNS6QAPOCYvPAfNwyyFMg5Zffyb5bBJbgMN3wzpC8Z/AMcvWyVeJRe7RkPWOyYv6el1D+r2KLWGD+URp/nrXgOAMfqSrwK+5Kuc9smfhOmPfD8og03dNXJMallWtSwiWy4f8dNswVWXyVs2+YbXiPoRVDNhOJzh+VBwDRF9HxFdIqI/DLTjRPTfiOidzZ8yZzhg2hitYiBwbiDPs3SmM+aYrDeOxcmyzdoDpDLpDJMf6hLSFmOFkpn8ongWWfaeC5e6pJ98jh8bZtlIzJy/DxKTh8rsXZup2WHGrjksJp/zk8/6zwOH45d+eHA/1xb6HY+j/JlfgA9JYD70QzD9sR9wO7nfSbdIy95jhs5sbXmQdFh36D8A+LOR9mUA3sDMDwN4Q/PvQ0uFLvDqB18uHssz+QdF/UWZOPPn1k+LM+SNwZro916oQpzh5/zk56YSF25Vps4y+ZOjLbGs59fkmDP3ZZj8vWuynvWTzzD5FbUizvDbIGTx10jGT34uIBlm5zUkvRQWMvlMLJock1/Rq+Ica6TGSzP5RG4Ya5zc15cBmEVM4LKMOHRm5udSPrZM+hUkM3wJW6CuMfjWf4nVxz4G+o2/AwZQv/qV0L//hyi/7V9l85G+RXI2Be/3zkJbMKP5r6/HeXZ/X9RGuReuzOSz/vAZ9q5RiHUrKMPqqRS6kfeTfz/b45WZf42I7o/kTwbwJ5u/fz+AXwHwpYdxPcAx+d9+4g/FY2+7+Lio/9Gl94gcOheX3jH51E/+xnQbg6JMjKxzMxf95JVSKJRK3ChLLfvJl7rAznw/mT0XpHApw+Sf3r0odtT37jwDy/1hhgE8vfuseP4zu8+BiJPzL0+uoFCpd89uvdv4yffbqDJzjIsBavQHejcTTDdXKaiAZZvWDap9AYRlIqI2Pnxc1km9J0Kc3Xq7HeA6lbFT7zYe4H19v94Tmfy0YfKxPrdzaKV7X08MzjN5z9hjY6cFrE4HSm+UlQZW24SnjgdA6YXhLiH7z7cxUxjQv/cHGH/O61D8zu/DvPTFUG/9I0y/+R9h/rrXYPAv/jVGX/Y1mAKoXv/apL3zsdibrS1ZOp+EOrh8UvYe6OJAn7L3Vue+7gfWtO26GDVhG9XoYhSFqeZmohLjV5ZDmc95DskWUcM5akj7BR8kvZBM/gwz+4hgzwGQp5oHTDcmO9gTgngBEGO6A27mLBp9mMUvKsfkc37yAm/OMnmLygg32VQik6+twVDr5DGp2WIocNzOIJcmX/746pWVY9m3i1UiXQnhhAG4l1emjSR7gMr41TO5QG3xV4G3BMSe4h0V7+s5LxN/bvhnq/uBNv7dLZh8rPtQF/EAtyimjRsU4pe2FnLvXk7SEVm/VZJehQyazjD8xn+O0T/7DvDJ49j94e/F8D/8EOZ/8zMwf90/AMDtn8V//1VUr3+tmNeiMr3Q+sG5ufTv22vXbP/K6IuZ/OxIDPJtYmYmkndaIKLPBvDZAHDvvTI2kdLGaBWjciAOqlJMd6enKzObMoidxTP5eNZpwShUkbhFMhgDLc3kNQqtE4OwY/LpSJtj8rl4OQQCkRIHek25QVjeXNwbD+PELA/0hhkDIW68ZUYptDejWcUaGWVbJh9dWuW2+M4w+YOkgzD53NoAktAS24xth+V4REs70C97fj4Vv/HbWP2c10O/412Y/Y1Px+Sbvw58/BjMX/7kZCI0f91rMH/daw/FqvCBkeTYUosmHnJiFEck1PBFIjoLAM2f4t57zPzdzPwYMz926pS8T6mUCl3gVfe/XDz20jMyk3/JmQyTz/jJn19bxOTr5H4WkJm8yTD5WV2JA3dlaqwJTN6wxUkhdg2Ds37yWSa/LsfQP7cif3AdHx0XO+qKWhU/zVXD5GP8Yqx1oYMjxtrGkw9OZ3bRNWMm77n1QA2F8ZExzjD2Fb0u+CV7Jp/+YqTHImMtqZT9p0mKRYL26zH1e7cih7botgvsn8+yvtAvPcfqIx/wnV2svO4rsP4JfxGYzrDzMz+C/e/5dvDxY03byD7juTjz4Z+3Ux5fr9tm8g2VSXX52nHZwpQbhFWmzjrjJ19kmLyG3F+KTD96v2TymfTTAP4WgG9q/vypw8zcWINff8/vJoydAbzlucfFydDbLr5H1N915QnnJx/d6Ke3LzZ+8hGTn+xgUBQuZEDwk8pWGBQlZomfvIImQh35RA6KErVJB8KBKrA930v0gnTD5PvlJwBP78qxaN6z/bSoP7H9tNh2T+8+K2KcK5OrKHSKcXbqHQy0RsjMGcCsnmNclgC6B5cB5yrpz2t156Fjvc9ycBEFlQ1eNqknAKG3Gx+zKxOhv42ehcVOdbPJM9CZsVNvN0rI0i12qx23EjfQLVvsm30UqoDiPnuv7AxaFWksElujVH3ezGBYZjDZRAcH4ZuDulm2MDDJF4834opcmetGD86H7dYvMFD+1/+Otc/9EqinL2D6ms/C/td9OWhtrb04g2FQI/Z7d/YG0+z93W8LZptc15XV2QLi29kx9nBg7WBcjtXHLw3/goyvbZtQGs4e0B1w7r9p32r7olDniudJHZw+c1Kkz+206eMq0mdg2KQf1Vy3cZAOIx3KIE9EPwxnZD1JRE8D+Idwg/uPEtFnAXgCwKcexrV8ujnZxbRy+5TGtygXl2We8ZPPfUrRAtQhXYNIJQM84DpMZdJ8cky+YoMhJCZvZCaP/Mwl6z9vjbi1IZpHIU5KyTsMaSXbLQqdbuUHOOQj6X4RS1q3jsmHyQ/sboYb6IG50iR5uX+bZPbkmbxgK2gGsn4dVDPwxds2Fr3rdHq3wUxfJ1EPY+uHx1omzyzrmX4cep30rnP1Gta+9Osw+sEfQ/3IQ7jxhv8X5lWvEPElN/0i5cdev919cm9V1uVir+dQVa7OaWvf+hqynuYUvpqk36fPYqPHNpzGT/79apBn5s/IHPqEw8hfSpvjNawMRphP0kG11DI2KXUphhdQJIfFBVhk2patu4bI5AeYR+xdKXUoTD4XL4dAzf6oEpOXwwsUwkpeYHkmbxkoBeOrYYsBDRL2znCfqbGuSWbylGPyGRLcekqIv8jt/SvXObt+guX9bnMGfKentbBNfJ+kTJyvg1ie2zwvzH/wEz+D9f/7K0HXbmDvy1+H/S/9AmA4zO7FkN/7dXkf8MNZG/D+l+R+uojJy79g8J3t/wDnvfHK+z5UPPaSTIyaF59+QNQfPJ6LJ39aHDg3huviTL6ARiXM5Bcx+cQgxw2TH0RMnp3P+Ekhdg2DcS6zN+v9mRg19+X85HNMfnhM7KgresXtWBXNjlVjXO0zU8BYV4+Ypc5tJTLZumHyYe6eWw/1kn7yRd5PXkojNU5qzAAGaiCzV9KI+YQf4B3C6Net/T+Byct8OhPTRuDQXo8xh7rwHDY//e9h86//fZjz53D9138B+1/zJcBwuJCx5/h0vOy/0xezeomxy3XmjI48qxfaQvq7E4CFfvIiY1+OvesMey8z+Rw2k6fcZ/77Ij322GP8xje+8bbOrU2Nk1/xMWKMl/FgKA6qw0KKbw4MikKcXBRaQSmVRGAstXZMPrpGqQqUuki2ABwVAxSFTlDOsCihdYp+SlVAFan3iyYFVaZGHwJhMEw7EgCMBqV7UYUzNGaMB0NY2N5bngGslAMAKasfqgKFVlE2jFK5tgg9e7ipw7iJaROy+pIKDIsyow+amCfdJ70mjVExaGfobZ2JHOpqWiAs7agcwfvj+2TBzWYi/XwYjFGRnu++yoYN21e9fAoqUCjde0EzM7TSKFQBilqVoBq9n79qdKCPXAjk9CCWvbdV+HziuhVKN189XataWBQNQiIGxj/wI1j/8n8Mms2w+9VfhP3P+3ugogzKZNtgcHFZ/b8p4srUNH/Mm7nJQ0X5+Hl83EYhe0dQh5xuAz3sR7Zh7HEbdc993HbdM9zT2cK7E/fbolvMF7eRD6GtItuO3w8hbiMfu8b1o/A+M+5ff3ip2TwRvYmZH5OOHdnYNdvTPdF9EoA4wDu9Xspdjigd4AHvJy/zZmmPV5v1k68xEsId1GwwEpi8YZtdbSdalBHE5YmOO4MZJdfgTAApx+Q5uYy0WAxwceYlP3kK/OTDrCzZ3uy2q1t/k+W2PByyXYlPSxgj5dyuzt1sW9Lj2DytP3xU785/Pj4/x967a8aM3d/T8Bd+YIzz6uos7wGs3/MkNl/7pRj+6m9g9jGvxM3v/Cbwwy8Sy+pRlNTHSGDvfjBPcacrU7ruwf9uEXuXEUZOzxPz22Pp/Wsc/PyWvSN2G/Z6hsnH/QgKczM/NGRzZHHN1ngda0P5s3yYiWkzLORGU0r13sphklwcrbXiXq4WjKFOr0FKodTp+TkmX6p0IAfQzNbS5Ji8zOv9TC69htxGi/aWlMpqWN7vtmYWr20zej52TW5f3uXj+BxWbP3sQ581fmf2JQi+WpJri3pusMkkY7D6Hf8OJ1/xp1G++fdx89u/Edd+/odhHpKxJZCv87JtmkvvP9zg+aTDaaOc7r8iDysd2UFeKYWPvOeDxWOPnJI7cW6P1weP3SM+uGfXT4qz1PXhWufXHaSCtPu6iLIyxqZfFywzeWbndbM6WIl0Rm2NuMcrg3F2VV5jcO+6HB/+3jVZv+sWTD6u86rE5BnQkJm8ZZnJ+/aUmLwineiWDQZq2HcNZOeWONYrgusesFKk/vAMdqxe+MFY5fzkZSbv44zHZSUgCdngv1o8f+/VAda5Vwo6BA7NnOrFW96O05/wKdj4sn+E+ce+Cpff9EvY/6zPBJRqyrSc33sujkvufK9LPDyrJ19xnr1LDJ/FtvA/zzJ5pH+X2gLAgrjxC+LJi0w+F9OmzPSXO0wegGPyJ75cZvKjUmbvw7IUB+1hmYlmqZW4JL3QCoUukrgzjsnrJKzCqBhCF2nIg2ExEJl8oTS0wOQVEXSZdhgCYVC6eOzhTIyZMWrYe5xGg0GPiwPu+Vgpyoax9vVSaRTKMXn/G7+qdaA1OIh3w8wYqALDsujiijf5DFSJYZHqpSow1J7VdyhGQWFYDgF00S7dw68wLgaOmVOHaCwzxuUQ3uPIn2/ZYlyMGk4cMnaLkR4CwfmubhZDNQCImvy9+6JFqcvWf96X1cJ5ymiloUl3OtuOvUeMXYGcXz36OpCydz8AFlRCNZtMh3Xz56uqxsY/+y5sfuu/gd1Yw41v/RpMPuUvtl9o/sWiqYBkt/BtELeR1/Osvq97LWb46B3r103aT9fhUhbP767bt0P41K9bRg9i18QM39sc4jYyzdqDuI38lqCKgn7Btt0PQUX9ouaqsYv1+xEAvGj90aVm8/+fY/LSAA9433DJ9a1vl/Qpu2cnyWETiEiMm2PZgBf5yQvlLwUmb5lRQOalLduNXgC5QFGGregn7zt8/CtFOSZP3Qyzp3esvseVKcPk2c1eu80d/A9cxs4MFjD5poQW3LObhAOBGCMHqX+749DdZhD9tki3PowfyLBuQNr/whg1vZwaLBffz3BP2PS6vr7poFW+8fdw6nO/EoO3vgO7n/JJuPrNXwl12q2Ejuvs6iegJwfZk7rdmtXHPamzE0gprV//JdcvUJ6x51i9nFdez+efP19qUyBl7C2TF/oXIPvJz+380JDNkcU1W+N1bI5ll7hRkdsWTtYd684xeTkWiRQ3npFh8uRm/nEaZP3kl2fyudjtku0AyNstdIbhS7MsIB83vrZ2KSZPJDP5XDz83GbawCIGmoutn9tnQD5/eSYvs3demsln0v4Ex7/im3DuT30a1I2beO5HvguXv/efwZ6Uw0MvTJnLLs/ql03vP0ThVmlpm0+uD2fOZzBGenSwwonXP6JJKYUPv/tR8djDJ+VYNA+dvE/U79s6C4k337V2Uoy/sjFYF/ck1VQ4PeqwxhqHdiJ+PK+bePKRXpkaq0XfT94z+eNRPHnPdmMm78t2d8Z//tyKHDf+zPhU7/c+HRvITH5Fr4ptpJudoVKWzpnYNXk/eUUq0U3L5GM+nTJ5z2pX9Yr06dKcnzLTkZb95EsaimXVJDN5gEQmb5kzus3qiDj06Fd/E/f+ib+Eze/899j525+Kp3/75zH5Pz6+LQdxuq2dxN47PcPqWda9j/lS7P0Q9PC/3vk5Jt9M+KX8pbZw+nJ+8rlYNJ7Jx9cokGPyCsUdJu/2Uz32Za/G/nyaHBuWpbicf1AWokvkQiav0htdaC26Dg50gUKnYYhHxQBa6cS9cqALt5Vg5NqpSaEoNYy1PYxEIBRNHSK3dwwHpTjDHJalyOqHg0HDYLtzLQPjogBRx7n9+QNdQinv49zoQMPkVYs2vD5QBYaFdtcO81Flo/fZe0Ehqw+4NSkMC2foDBk7cffFFrP3lXLsmGmP1VusFKPe+Wj0oR6253v0YNmtXiaiNhIm4D6tC1U28Yg6lmrYomiYfMxkHasvGvLaoR5FDb+PWL3/Mgv1ez/187Dzca/A9uf+bQBAcXMX5//q52D8xj/A/MF7cek7vh6zj34F+m6Wnr2n9glvN+ixerbQyvlsq/ZOO0zmWyD06W7tKtTn0y1jb3MKXTzdiyFl7GjsPZHO3LZJyuS9HtStjYHTP99y4D/fawsT5N8FsDZs2r4Yt5Fn8hT0Fz9R8W3U3X/jJjXghL3Xtm7tHWF/AYCHN1+61Gz+A5LJ70z3UZvbY2I+uUFzSSYvvBRyIQSISFwJa9mC4w1b4fZyVbrpcMF1DCxKLtqyteVsHhAgZeO5z13PsRNW3yyQSovlZ54RD6buaHikUN3er30mTyKTB4X7znZHLBkwpzFeHMtOGXtohA31PpPP45NejZkBSmPatOfH7D0Y2Ht6xn+eqBss4yR9HSkoMPX17Y/7SJz7mn/h+t7p47jrC78RamcPux/zGJ75kX8Fvbrmay1eI2ufiMNOePtEVFY/NKe82bVGzvaT1jo1ngL9+5bcf0FHT5frlmXvmbbIt52wJ7Fw33yZbpe9d6w+drDQqOz80JDNkR3kj61s4PjKJp7buZIcWylH2K/SGf5KOcKkniV62XjKSIYurVRiZLVsxR2gLDOGxSCJUeOYvE7Oz+1RWyp5g+ec7UDBrZCUXjBDNcDMphubD/VQ1AtVip3az3PiR9YwMBDi4NSWMS7SGDWWGUMhds2iePJSWszkFSSDn4Ico0hTkSxgAfKG9xxxduemcXCc3UJabOXiz6fGzGBBVJOuvuYzoXf3cear/lnzmgeu/YPPxKVv/OIs2z1IOiw/+XyJWDzqZ9PvT0k2KC/Sl+t3uXyYbbM6+3DSkWXyRISXnXu4p/k364Mn7u7926cHjsv6PZt3iW/mM6vHUVuJya+JTN77yXN/+t3GronZ27yukkHMM/mVMsfkNxLdgnF6fEJsi7NRnPlWX5HZ++nxSVHfGm65j+1IH+txwuSZGw8kgb0bzsWuqUUOvYjJDzNMfqRGsh6x947hj0U26vOP9TLDWF2sE8EfnjPsHTJ755jJM2Prh34GJ//tf27P23/Fh+LSN35xmzcE3uxZ9PNh8l091QvK2H05k2tDZu+cOR/BNeLrin9nAEJbOFVuC51ZD5Fj757h3+75AB0qkz+yM/l5XeHXHn9TMugwcxNPPr3Zb734ePf5F+jvvvpUO7kI9Qu7l7vPyEC/OdtxTD4a9GpbQ2sFGz1YmhQUqwZT+Ot3uzMln+uksDffFz4HgSuT62Kdn9656GZDIWO3jPfcfFrU37vt9N7+osx4avfZzk+ePfkELk+uQlEz0wvadme+i1I3sD7QK1M1s1RAsf8sZWhmsPWcF61eNuWwzAEZdYbDytaJzpaxV7s1EqFu2OJmdbP9CqMGPRk2uFFdd3rDRp37ZacrUi0asmyxXW2DqPFcCvR9s+f85AOd2WJuZlDKra3wMcJbd1HrOLh3j7Ts4ohb5ePDd6jHu5kSCKMnnsXdX/hNWP+V38LsnrPQN3ew91Evx8pv/S6OfecP4MprP9OhNZ5DofPdt7CBrzdBQTczR9u+lAkExaqdURoYsPXnd3oYW12hayN/jLzfPru74N1aWy7O3f3xGDDRGzdcfw98W/jrxueHTF7Kx+kye1eB3mPvQZ09L5fayNi6Y/Xc5V+x+1LXCHRYzOwMDE76UeXjzEf9xbLFpN7HuJBX9C+bjuwgvzffh7GZzYIzTN53yDjJH5D5z6lFTF7ynzfoNkgOx+fKGtch4vPZohRujX/LJ0yjrUU/f6fKemu8jV8YkP3kQ+Np+ItC+Tgr/VKFTD4m2p7Jh9zf+8knLL1h07HeZ/KdHrJd6V4wOImpw8wNk5dtPKneGdX6cmc863+ed/7zPbuCjzPPjJ7fOyuwrXHqe34M5/7JdwOkcO2TPg7Hfu7XcOHrX4err/lrOPGvfxDnvubbYGBx7TV/3V0XVuz7i7hysmag+d9k20bmJtaRbJ+I2XuOseeYvE+5Z3cxk0/zWcjYl7VbLMHqgXRNgi9frn/F/UWRiyc/xuEM8kcW12yNN3B6XfYDXhvIjbMo1o3IukmhULfvJ+88MlKdADGfgS4gDdjFAiYvxftWIJQZvj/OLKhYyawZGKqB/MKLdkjyybC8lqC2y8Wu8d4KcVqeyVP2N7m1BLkNk3PnL4wzL7Rezk/e7+YUp9FbH8ejf+7v456v/nbsvPrD8ZZf/0GoyRxPf93n4epr/hoAx+gvfP0XYP1Xf1ssy61S3qc7w9hzrH7J/XHlCcr7Oh2O37sSJmyLdTkf5wkmh78+SDqygzwR4dFMfPh7tlxclhjZ3LN5V/Nv9P48u346mFF16eR4S2Ty6+WqyN4VKcyNwKGtzTB5YWs2dkx+VIx6zwOzm5luDdZ7ZXFsl5O9X/21PKuP2yKnnxgedx/CUd02B5tCnBUXcz1l8u4TtLbLMfk4Bg5zw+QFNmrYoIjiunvOPdQpk0dT1pCB+r8PM3rsD+//XmQYq0KK3/zfc0w+bDuazXH+m74XL/3Ez8LgiQt493d9LR7/wW9Fdf4M3vXD/xQXP+fTete9+prPxHv+879YzOQz7D17PpNYZ3/+rfTw2DIs3f23zJoB2c7hj0n3IPd3yupKrLOKWLr/e8zqWz1nw8msq3AutIcHWY6sn/ysmmPrS1+FaZ16iJSFhrR/6aAoxHlEod2bVhG1S/3ZMsqicDOV3r6kaHT0ZtvMzltGaYJpttbz+rAYiFsGloFffXi+Vhplxk++bOLBxI7yZen85GN/+EFGd371STYoC+2YPFG3AbW1GBQFVBPHxfuPWGYMdQmtXLuF+kAVKBv/edXk4/WhdoOhavEGo1Alxk3dQl03fvId++3C8A6brxHd+G473uz84WM/Zhe7potpQ40NwLBp/e2977j/dB+qQcvkfVtYa9z+rqRaHXAL3goqoJSgN+crdLGQOv95d82NN74VL/q//ylW3vkELv/VT8RTX/95qE8ca/3bPcbwXx3et9p/6rexa6hj8pbdqmNJb/3km+u7wdQE/vMKRM6QbNi2bebrEaKh0K/f19n3rTCOT4yqFHy/8LHb0baRO7+L6a68zQR9vEGR7q/hbCade7L3VfdlBTzSq+FXHqtmRbllhgmZfNRf/DXi/uI9xkLdxa6pO/Ye2GS6mDbKvVDAbVleevzlWMlsciOlD0g/+f1qInJxoGO98YDuOHT+89JyHyy3sTqijCR/YXd+F8s+fHe6/UvT82troJsXTO98tiDk/OT9P+La9Wfk8W9i3bdFnE3L5N2Uq6sbuoE3rIlqXnbMffauQibfyweoBd1txtAZ7YKSwjN5FsoTM/Z4IU2cWvYafSWFMXXiFJ/fG7Si/iLpXU0iZk6A2t3D/f/k+3D23/80ZudO4Q//0zdg+xNeBa000Ly0pLrl1wYsyZsFt9uujUxP83UzwkYbC9sIqW45jtLvkmxX8P0i3JiG2z+lZ9F/AUh1i913vb64LdKy3q4Nx5c1tdV0ZQ2fHsfk5bhcB0lHFtdsjTdwbkMOr7s5lN+AGyOZc42LocgiXfTAlKcxWIwPn2PyYG4e2n5ycWXS3lMonY7hjS7Fblek5Jg5AFYK2d92tZTtE+PMAgzdzLziZDnH5OXYNYZlJo8sk8/wzGCVoFRWKS0bW39RHB8pLcvkt37pt/DhH/d3cfbf/zQu/F+fjDf/yvfg+ie8on3RitcQ1VwiLMubsxw6F3M/c/4Ln3LXXb48h8XeJSeKRToFq4fDZNlitVhPf3DAdGQHeSLCQz5Gje/5nrFvOt9wP1D6P8+ue70/u71r/USz4LqbHQBuwZXE5Ff1GHXE3j2Tr0Qmz+4TNtIr45l8qo/0oDeT9kx+fbDW46meWx5rYtr0eCCAk4IOACcGm6nOwLFS0AGsF2s9P3nPVwdqKDN5KNlP3tpM7JpabqNmdhXrtTVJrBA/4x8ogaUDGDT+82kdRi0r7nRGKeXDjAID8XzKMHnnwtnp5bWbeMnnfSs+9G9+DczqGL/7U9+Gd//j18KsrTRowYht6rRUt2wBm7YRN0v5RZaODHvP6YzM+fLfObpue8zKurU59m5k9m4NrJVi9Mv5OO6YthHQHIrq4JLcFirH6jN6HNOof75Uh/xE5SDpyOKaWTXH/3z37yQDPBj4o4vv6f4ZDPRvv/ye9vdhg7/72jOdn3ww0F/eu5acDwB79T5Y8Ku31rYbFoW6Vm6uEKIfBrst9ZKO5TwcJvW0Z1fwHfXa5Gbr/uh/Yy3jws4lNwOM9PduPwMALVP05Xxi55mGm3brACxbPL1/ESC07pXc5HNldr3nJ+8eDIv9er+1ZfiyWmthlAFsFNOGGUUTSiH0e3d+8hrW+pguHRpSUKise9C16nQwYWImQDvr9Q+kxc1qG4Rmts/+U9zgRnUj1dngxvyGu0+kWiOcYYud6mbHodvzLfbMbst8vW7ZYsbT1lfdzwING2iw85OHwrmf+h945Gu/B8XOPh7/gk/Hk5//GcBoANUiLDfr91tV+lmgRwY+dozPn0M+zYBSutO5i5Xu6tzxY2vt0joQsPSmjQxqwLN39rzZldX7kvvzQ3zifNI9Szfts9HX+37yKtTbvPpt4fdZdeshdHu+d3dsGb5n7DF7Z8FPPmoLw9Omv+hW94w9biPfLxx711DB+ZVx8eT9+gnA2TNAhL16F2vl4czmj+wgP6mmC31qRZ0XM/nbTwQWrk2UxjPx5fFlDctWWyOGIA65Xe8tH/yWoyp27F0uMUcHmklbOquCXw3ZRwM5P3lvt2DuX5yCevSuHbRR6N/uByvLEZMn/43FqINAO6p9Kfv8/WDQ6fK9SHXfL0zjpRInx1LDKtBiPeLKBMLomcv44K/8tzj9y2/G9Q97GH/wza/B7MUPNYt/+udr6DZ2TW9gDPaQTVm9uy8LufLz0H0buZl1eK67vjs/ZfXx+WEbSXaqnO5fsOG/m4IJ97npL9l8Usa+bBsBi/ze4zbiTu/l3eUT5qVRoBJCjhw0HVlcszlex73H5HC5x8Ybor41Xhfp2+pgLPqfF1SIn00MFuO0m4z/PLPMNPP+8PK7t2h2HYqTIsrGyl/LsPf1zFqC1UzMjIIKmcuy/GlZWUYhxeK3VtQ50xVVs3owTuRj3QhHlmXvA8GeASC7tDzHpxtv/75kLe79/p/Hx37iF+D4b74Fb/2av4P/9ePfiN1H7nNeWBKThczk4xdsfO1cicQ6ZNt7uTUGh8vql/nN4eUvldVFH80w9kz/yuqZubRkXwPc5GG9waaHkY7sIE9EuGerGeRDJgPg1Nqx5p99FHJyte8D7v88Md7q+dZ6NLIxXHOfktH5Yz0UfcAVZCZvrZV9wI0BCSy1NgalKt1sOri2sQar5ajHFj3z3Wz85+Oybg02RH2zlPWNci3VGVgtVnqs038BFDTo8WOvE0j0h6+ZxTUDdRO6ID2/7uXbtgUbKE5Zp19slTBWTuPAt3XAQGapkP2h471cuzr3/eFX3/U0PuZTvw4f8g+/F9c//BH82n/9Nrz3s/4CoN3vYz/5rm51lslLdgu3FN6IZYrZeKenbcrN552k21w+LOdjM/mbyO+91a3tPWtt3axcZ5uxW1i2LavvXcN9T6f9AsjWOfd3yqwlyLN6+XyPUKV7dofJA5hWM/zGu383GeDBjHdceiLAF92f77ryRPv7UH/vjWeS/C1bXN2/3rH64PxJPROZfDghCHWtHGk0se790BlAW173pTCr501H8Dff4Y0b053WZ7x70BjP7V5yOQTsna3FkzvPopFb3VrGk7uyfmHvUt8pg11bXJneACk3K2jZOzP2zaT1kQ91CwtjHcoJ2XvL5CO9VKqtY6jrBqN0dW7QAXPD5L2nTbdAZqfaCfSOjd6YX2/avdMNW9yovO5Yqo91s90wed3j0KZj8uj7PbNxPs+6snj0e34Wj37nT8KMh3jTt34Onvo/P65hsrY9X0P32Lsvpwub0fh7+1g0LeKqYbnz3faxTnyjeU5s2fZj0ZACs8OMFvaW7D3kzf7l0rYRB37yXLX5+EHOlbVh9eji9bRIovmy9XX2MXZi3TN5htvuMbRP+PYI90jtsAc5+0SwliCM4xO2Ueg/79vCn+9hXNwWtmHyqsfkDSZ2v9dG/nzvJ9/GLuJO9+XvxTSCwk61jY3B4czmj+wgP6vnefb+PJh8btAOk18wIfxY9INlyDE5DNv2BoTZLWLvYjkREI3MDxLZTdoS3Xf48GsB6BZN2SizLlZ6pKOLM99j79SxzD6T77aUiHX/oLvydLMvX96QsYa3zK8EjVPM3pm9PcXAiOfHep69H/+9x/HYV3wvtt7+FJ788x+F3/nqvwF7+mQToz0938+0Jd1fO9R71+XwfNcaWR/w58mb/bOT85OP/dJbP3kYuQ4ZP/ZUb/pX1EbcdGCGnM8i9r5MGwGL2Xvu/LCNwnsZ5tWy+qQOhJrv+Mljc7zehhSO08nVLVE/sSLrG8NVkcmXqsz4yct7pxq2ov88sxXZpVY6CuvU6VIqlBavq0hhtZRZ+uZAXjOQW0uwXq6K77YBlRn+Ku8vW2fYu22iMcYp9zKjDJMPVzTGvygpw96VbLfIbc5QKpnVS+XX+zO8/Bt/CJ/wKV+HwY1d/M/vej1+819+LmYnN7N+8iaj585flPJe40sy+QwmyPHmHLc+mIOD9BsWdW9svl29O5YmqS3a6KNCKjL9q6SMbSdzvstfsr1YbJZb4m8Oko7sIA907D1GNlsjmU9vDmV9Y7DW8OB+dqvlCoxN+WFJRY8rd6wOqE2qG8OoBZ/e2hgg8ht25xsUgm+tsabxn+/XwS2eGCc60C2Gul19tdnXtKezixsvsVSNQmTylklso5rlOD514Pcc/saFGU7ZaG0NpDgrltOYIG1Zl9Uz7B3cZ69nfv0t+KRP+ho88u9/EY9/2sfjF3/hn+DCn/rwoC0yHDpi751uRHtG++kf94uA1cdtYSQ+7VGIoIf/3ZYucG5/bUk3NsfkO7tCr7/YGkZcbyHX2c+Wb7fOAPJ2i2XbQuin7u+L1h6kbeG3hjysdGRxzWQ+xW+/5w/cksvGp9sD7nddebJ9l4cN+K5rTwJA711vrcVTNy8G7L378+rkRnA+t/qs2azbdxGvE6ngZdHpWqsG8aD7sOYgPoztG2EY1AY6C3XDjJuzvZ6fvO9YF/euujlPpD+ze9HVIWL4Xg+TtRYX9i8lszDLFldnN1r/eR/KxzJjaqZORxe7hplRkkXtwo704skXumj85Pv+8KUuYNnhANX44rsBW6FmH2+l0eEejomZAbAdS0UQBx6emfpPcRv4w/f1m5XT4xjqN1sm77m1++TmegeKFMY3J/iob/lxPPSTv4Eb95/GL/7HL8aVj3rU5dOGU7YoUKCyLoYPsWslZ29omDy6+CumiR9jrIFBZ4BzfBqAUoCPMcNuq0j2TJ7RxixvNyQh91Wgle7xZq1q+K8wVzfbGnB9GyHgzYYc4vAcGp43N51Cox+vxT9gPlZ6qzcpXkvQMvlIb5+Vhl37tuhsFKZXZx+H3yDQ2ce06dg7fJ3bHeHcffZ1NrZukYumotdGPftEE4uotnU/dk3A8B2T7/aV8P3O6TaJJ69I4+b8OraGcpTdZdORHeQrU3dcPOFumdR8+bkxOHtWm1pmKl6jG6xbjTsm39MRxNPpnd9/i3fn2/aCnClnLHNGD8uW/EAyHi/6Dfy7NNQDH+3g/HDxVxzJww9YKXsX9OBVWkf5O1/8vt9730/+dnkzWmaeZfhhjRl44L+8Ga/6hh/F6MYufvfv/Rn8zmv+HPTKKtwWiTl/e0baGpwMgNzWu8+PO3/4VPetlOPNteA/7wbqTEz0BX7y/RxcmxvIsWvieC1+ApGLfS/ZLQA5do1/PrN1Xoq9L2Ly+TaSfpPGrvG6heQpn7QRkOx58HzSCz7IE9GfBfAvAWgA/46Zv+kw8t0Yr+GRM/fhLc8+nhw7s3YCF3evJvqpteO9Vaw+bY7WsTPfTYypQz1wy+2FIEulsKeqn5EmurUoigLxIg+l5D1Ew8h9YSqUxkCXmEb71CoQ1ger2JnvReUEjg03cH22neR1bCTrm4M17FT7SbnGzZ6wsQ3B7c2aLkipLWOopb1fLYZFkeiW5Zezj84XJ8/k04GCMFADzIXFJCM9cl8eURoXY0wbT534/Jf8u59HuTPFIz/3Rqw/dx17pzZRj4fYeuISrrzkHvzid78W1158DwDAxXtMN5qxcLP5RGdnt0h1gwI6P1kRUo5QL9rXVhq8cnpu79+CCnFTDakdbpWk3/hZdly77sUm5XH7+QPu6yR+eRERCpSiATTXv5bVNRXNiyjtL8cGhzOLB17gQZ6INIB/BeATATwN4H8T0U8z81sPI/+NhrH3+gG5xU05/fJeN3Pzf66UQ2zPdnp4BACGusDMzBOdQC0PDHXHENNYN9baxNcXcPxeB/uXtuFMrcWASdapwKRx4WoRDAFDVWI7Oh9EGKqG4aMrK4EwoLzO2EvqVqLABLNEV+j84fuhExhGpXoNi4LTtjAwqK1qH1H/m5oNrO0Gwi7kq21Zfa/OaF4MiU7wLD0+n23a1u5+EsqdCR77d7/UvnLXLt0EA3jmFQ/jF77n84CyCOpswU0IizAfwwak3Aug3xY1iFODoGHjthGMTGYGFrAd2vF5WViQRaIzLAw0gLRu4cb1z0uHrNv2a5SSZ8ThvVi37bm9fBqXyNjI7vl6UudmVuyZdqgbcBtkr2sjRhgPKLw2wKJuMnq4ziOum5g/u/rFbWEgO2ocNL3QM/lXAHgXM78bAIjoPwP4ZADPe5Dfn0/wxif/MOAU3Z/vvvp0d2Kgv+da5w/ftjMDz9y8lBi5mRnXpm6mG84A3MDTvX3DwUoRtTc6jAmjC+dnbSK97dD+09N3OljUbIOXkWPshi225/0B2LP6i3vuyyWOUXNh9xKaSrSDOnPjD4/uZec58XP7V5pydTozcGO+0zGStu0YUzMHiNuBDWi+aApGZRxyCNn7qChgGEkbFY0fd90weeKOyRv2/uMdqweAeT13sT+UArgz+O1UuwCaFYXscQjj5vxmq4cvi+0qo+MmHvm5NybfVARg/akrqDSgmtjpPmaKanzY3UPqXszuq89t2h7Gky+URmWdn7luZvSunhq1tYDnyugGO1Kq4badT7eFbQ0lWnUc2rIFq0b3vuENDul4c6RTX/d4Q1El6n6A6uKyOLzRvNKacjZtwda1JTXfGOTtE0GMGu7rgMedph24QwOtbdoCbVt491zT2mQ6P3l373oMP3h5FVS0z4j3YZfqHDL5UPfhgdv9B7gzogN9+4TPx/aYvGsjTRrXZldxYnQSh5Fe6EH+PICngn8/DeCjwhOI6LMBfDYA3HvvvbedsbEy91yY/F2OU2YRU8fkU5aQu7SPadNj7M1sO9GDf4efkX5w7f4uXUc+kNN95+1dnPrX6ucj63E+fswP69Kd6mey/bq5l5l8PiNl+74kkg6gt9dvTxeYb6uHTdG8TKXYNevP3RDzWH/uOtLYNS5lfcAX+HSn/Fhi753e9z0P8hHZ+zL712b0pr8sjNfSUzomH8ZK9wsu4vg73WUWxOVZeg1A/3yfk8TqAWT1A7Wd8OjE7N0/H1Kc+VxcroOk97kLJTN/NzM/xsyPnTolx4eX0vpoFS8++6B47NyG/Aa8K6OfGG+JfvIrxUiML0GgzN6vLMadsVa+YSrDCLXqdtnp6aQwKlJfXAVq3UP75XR1k9LJjH5suCHaA1b0uP3U7ZdJjnFvrLw0u7ZG/BS18UuoTbJvM0Fl1hMQhtl9bTPxenKx9YsV7Ny1JR7buetYouVi3FuWY9EYm9/7dbnkP7HStMjXW5oQlFSKA5QY94fzPuDZvXkz5fRH00vIMzM/sbidPG51bWkNABG5dRLCRYZKXlcxzKy3GGTWW+RiKTEYx0cnxGMHSS/0TP4ZAPcE/7670Q4l+W3bYshe+psW620Y1oBDE6HQRTvDDPm0/4SOuTJAPf95n6wRoi7CDfLWqkSvrQUJurEGpe1Yb8deWeTNTN0naMjYQeQ+AyO9Ddsq6D5cbMLemeS2sIBVqe6YfNpGNWyze1ZUZzYwNjWwGTaojDPu9nRY1DZloAzAGImNEqxwvmtvmZkaw3jLn/swvPJ7f7k3PDCAt//5j0jqXFnrwjFH+VRc93CWT3N2S//jl6dhAxIWe1m4uPHxwG3Z4ZqUyTPmpmonMKFemar5kovrULdfrqIe1a22deummdgbkLL3ylagJsxHPx+H6eI6u7APKplI1NaAkW5M48MRSG3kQ0bEbeG+fuJ+1Pj6S23BVaYfyWzfh2FI2HtTpqTtGtx1WOmFHuT/N4CHiegBuMH90wH8tcPIeG+2j9956q19uN78+cS1C+7v7ZjBgGU8ef259kUfGoKe276c6mBsz3faISccrGrrb1rEoUvZk6HQjQ8xt1+rzd9dp66jPWE9E46Rj2GDvfkkQDie1TOu7F8PrshtPs/tXXFNEbH3ZzP65ck15wsfGctuVjttG4UsfW7nbXmc7K5NBRquzL02GmkNY53doecnr4rW7tCLJ0/OlmGtcewdnYGrINU80NQxeQC7tfMycoy9iaXDjJsBe+98wJ1fPYCWc3t9p9rGZGWA//V3Pg4v/cXfw8ZzN7B9Zgtv/fMfhmp1gJpN4DNu4bdUdPu/UlMfF2xuZqveIGZszOS7fU016dY9L97jFdT4gKvAf54toErUXHcve3YGWesC3qNo+LFp0ECPvRO1g2BFna+3tzE5X3SZN3dMvrE3BLp7gal2suQGLwWYJkYNKfdStqGfPLUx8T375sYfXwe6dzFt/eTRsffWTz5oozYvcK/Ojou7DlwEun+JSHX2jF03dggLz9j9WgLV9ItmkhLoYRv19n6lLu6PpgJXppdxenwGh5Fe0EGemWsi+lwA/wXOhfL7mPkth5T3bZzU/s/ilH1p5hzTuuyTcsgYu/vfkAW3FDpm7xGrzvJ3udSSnjD25is4Ze/uwKL27R/r2ijH5Hs6s7suYn/47k8jXJvhXxrdVX1Kz+8GXul+5PcGTvXf+lsfCwD41df+6VYb6GaDd8F90+djWdJZRDFuALx1bBkfyNiiz97DWV+OE+f8rnP8OI7X46HJQt68hB0ifz7DCnVL7Q3dZGy5tQH5OtQL2iLbRkL+OZ5+u/F6iHLo8mDpBWfyzPzzzPxBzPwiZv6Gw8p3bbSKDzn3sHjs7i35DXh3sy1gnE6vHhc/j9YGK9k9VUUmb2072wyTyTB5ygymWmnRRqBJiXHjCYTj401BB06tyP62pzP6idExuS3KFVEvM3HmmVlsu9ra7nMmSM74JHXsDJPP3AMCZfepXSvkeD05fbWU9wQeqkzMfSWvb/DeU3HKMfncQO3zklLu8z4bW58GYt/LxfcZCjozi+e7e5+LdZOzsch2hRyTz8/MFjD5jFtiKdgniAjDZrvIOI2lPRcYGOvUtsPMGGX2aBjooVheZsapkTxWHSQd2RWvQOM15meiAXvvXFnQTTQJ8NOrmL37l2jMJx3blRm755mxbm06cFtrYUnSWWb7wn6Xrjrc48ptO1CGNxOBGy+kHnvP6QDYWpG9m1q2N9TWgFVapto0bYT0fLeNX6SzhbYW8VIWA+ti/ERM3rLF3JgkfwacLjD5eRDrP6zbrJYZ66yu4LdgDPW5mcOSpFfiy7myFaCKZCCuuAYRkt94Jp9waDZgSg3aFhZzm25CYcGYmbnI5P3inPA+O4Y/E9n73Mwbm1Ck2zkYKVeu7ByO2PTrNrMz524ZlbWyFYjSoGC1x1mRQ4Pn6LHx3SGbNDieQzN10kaL2mLGuTaay/3FzEXGXuV0W7V+8mGaCwuknk86soP87nQPv//0O9LYNWA8feO57sSAAzwV+MN7BshgPLd7tePNrQ7sBSs/O3bddSR3fucho5Ryg3ZU1kI33jLJfXOzXc8Lw8UtrdEnuEbNBvtVs/crNaVlt+TlWhNnp825YdSX9q+J+sX9ZkWw59nNscuT6y1DDxn7br0fMXmnV7YGW19/al+KutSojOnpzk9+4AymGSZf+Rg1ja5AKMktMtPOWNDFPFcKPqYNNfzYtc2kadu+vt3EmU902m70vp+85Z2ANwc6rKgPVLP2AZ69e73EPGDvDkmZ1n8+jF1jGyZf2xo1uoHbR6zUSmHO8yimDYMVo2K/G1ezTR8YBfW5skcDRkV+8k296iBGjefNlg2q5nq68Xtv/ec9kw/WAJhg4Io5tHvz1+58pQF2/Zqa/lf7AZ07VMVg1KZu69z5vBOMMV1s/cbvnECYNfYJQhdznwFUzK3e+c/bpC0Md6y+H9+n4/vhqmsfmx5A0hbt+c1uZt7mU3s/+eb+d/2ixKXJJdy1chaHkY7sIN/zVDjIS+82eDZn9H420aw6m7XwGRlJnaFz8TW6c3Ksfgmd8/RvqXzgvkz6qVkokug59u6NVensvNX9DCo4EOr+GAu6T3k9F3MmZuyL9UW+27fD3n0t3KK4NBa7K+vhsXcppX7v/vzI792ff9vsnbvzexgzZO9VoDqMk64ZaCYavu2WqNth6UvbIZKYRr4tZFvNYaX3uZ/8QdPqcAUvv/tR8dh9x+Q34H3H7hKJ3dn1kyLTdHHmpb1ZC5E3W7agDJPPDoxCZ1BEIttXpLqQDUEiUDZW/l1rmTUDq/KahNNj2T93o1wT22ioMnHmbeoyCKBBL+npfoYcJyYSr0tQYux+ArVhlNM6bGT01J4BAOuDdO0BkGevhSqzbFxm74uYvPRi7L7u4rxzy+BzsfWHaij2yZEap+8/lmPuM7PoG87M2djqi7e1yzH5NOWnJjLbB+RY+USEgdAWRISxHotttKJTWw0DWNWpbYeZs/0l51fPwKF51gBHeCYPAFVdp88CNTrQn60Roaqrhuhw7wdz0+0jGqZ5XbuZbqRXpoZbMZ+ydKUkv3fbIItohmoZltL8GY5pS0y+CvhxVzXCvJL1WVWnOgjzet7iql6dq4xez5tZf1+fVTWgJSZv3EblcVuwgTYEjv3erUVNNvGH79z6Ip0sZnUtPuwzoY0AwrSai209reYte+/rM/Ehn9bTlsn39Rm0Sk2LUzPHQJeJPjcVSKdMvmYDECUTCe9WG28c4zxhGEU0MWBmTDCFRsqhp/UssSsxM6aYOHzBFJ0/bVFImCb1FESpYXm/nkBRutH31MygFCWLpea2arbJjFh94zMeM3YfLyluC4dATLLQyLCFMVWG4ddpWzR19u3S0+2kZeyhPjGpDgATM3WuwULbWbbp+olm688FNuSl0pEd5Heme3jLs+9KDzBwYfuSoDMubF8RGi72Me/SpJpkUIl1hsUoKVKoTd3wwrY4KLVqfpcO2s4I2s/HwMBy0ZobfD6WLaZ1OlAZNrgxdbw59Ld3deuYfOgTfzmoc2hvuDq7GTF5d2y3moh+8oYM2HByflmSm7Wj3xi6KFBb50nTZ/KOz1bNohivayJYpXv+877+hTKwgj5pjIehzgC2axfTRpH/2EePySvqokhaZuxiN2Ds1N4Dbpi8n0F7xjpQDMMpky9ViZmZB0zes/qiYfVofMndvSyogLEGNeqAvRsACqQIEzNrBz2HTxiFYtTG8+bmhQBGgRIVOt17MXk/+XYf3KYOPhRArNeZ81t/eL/HKwImzymfJiL3lRfZJwiOqae6AiNdM+CdJqZm2u6w5BdCKShMzLS1T4SLkUK238b3ab4gw7qF/vNJW3DH2P1aAgvbumy6tqD2Ptctk6egLbq4OWl/GeC5ybM4vyrvfLdsOrKDvOTFcMt0SG9GP4hLH4xECsyR8TWHHFQORYQzhNtLvjw52wEWHFte71krsue3bRTrzVslDO0cutDFIZ99a/TP75K43+6BdRmTyL77y+2dKjH50KgeiN1gE+jhHDdl76610+u6ci/iyoelS7htcbyeW7dRM91Y0KZ5O8RC+0SmDkvpWGCrEc+P2XtkV2iTHBLkoOlIM/mPuO+l4rEHT5wX9QeOnRP18xunxcH22MqmyDqHWubQ3sUsTsYYcWBhwR3SpS4EQ5gUKawPU75HIJxalVn6+XXZ3/b8usz87lqRGf6x4abYRiM9lH2frfwCM0b2hq+DGVWYmEn0qwdIjAlCIKyVq+ILamuwJeQDbA2OifqmcD4zsCr41TPnfcwJmbbIvGzCyIi9a2Re+YuY/EDJ/vCjBbw5+eJkYGUJH/CF/vOZeWXOlpFbPQEs4vJy0lQkdcv5wxOR2BYAsFasp84RzFgXbD7MjNVCYPjMGBcykwfoDpP3aXe6K0wfCXvTqTgFnkwn8K5aYdpv9Hgg3ptMHJ+O9GldiYzdGAsqBA5tLbRwvj+WMnNnD0jPZ0yqqdgh92b7InvfnU1EPr03mwjsnbA/n4pMfn/uNtWI85pUM5BO61aZCoVKHypjDUhL6Mqism7+GiZLjLmpEx0ETOoZ4hvNAHbmzoUyvASDsDufpOydIesAduf7gbtqp+9XE/dCj8alvXq/WRDVT5N66ph8zF7NHNBl6ifffN7Hi71qdpOFWHesfirwZsZ+tQ8d2UYsW+zV+66cEZPfq/fg3Yu7OjP26n0wJA7t+l1ct0k9ceExEj69326bF6apmTaLAKXYNZL/fN0Yu/tDmGWLylYodKpP7KTdO7l3PtdQEZO3YOxLbcGMPbsHH/8+bIv9ag/Wpox9r9rv4Safdud7PazYtlE1QRhn5/mmIzvIb0928fbn3pseYMbFnXRXKAB4budad18C3nJtsi2inKmZtfOIkI0D3PqAh0mRQtUYcVMmT8n7yMC2m0WHqUYNDSGgGRu3oCca2Awb7Mz2WjVk79cmN9vzEl2wT9yYbXedNGij/XoaMPnmbAastgGT7xh4UejWAB7qw6KENU388CAjauKgV1HsGk0MKJ3o1LS3Z/LhS1o1vsihzgB2sdce7+kZVr/Lu40/fLeIxTID2vvDo1lFyY1doePLvqyGLQa6bNm794e2jf+8j2sS+sMXpGGsRWX7TJ6IoJRqFjj5fYOb9RqqRN3TG18ccvFxvC+5bRa7+d+1vLlh9XWGyStyXlHOL76Js+IXQaH7mmhZfbMuJNEbRh5+fTh/e9WLx+51789eoerpgGPiMzvr+88zO93M2jp0KIma+D7ueWyZfNQWXk+YvG+jMEZN2y/6a2dCG07VPLN+nUO73qKxTcX6QA1wYe8C7l2//dDri9KRHeTdbGZJyB5P1W6RZPLsBrdwEVSbPQEaqjF6hTo1g0H/xaAyIYUJCkvXLSrfLfVM9r2Z1+3k4w0USGfnvo1i3QcdCzOiZiAFx+w9z+p9euH04Loik0c7QPh/e8raQzHNAUYaN6dj71HM9eZt02Py3Om9WCfBPVvIlTP8OHt+v6BOBwNL2iGW0YGYmVPQRrdeM3A7bbFM7KJb6kKbLrLheCOsL61n8n0dYsjyg6Yjy+RXhmN81AMvE489dEp+Az504h5Rv2/rLpEJnlqV48yPiqGIiZnlcbE2BtamHdou8J8fCD7gigibo5QHEwhn1mWWfs/GXbK+Luvn1mSGf3K8Jb4XVvRIbDuViVIirQj2uryYJkUBLhFGQtx4AmFjkDJTADgxPC7eoBND2Z5xfJTG92E4JpvoLMc0cShDjjOf3dRE2CoS8ANser77apENdQMt+8OvFDJvXi3WBSbPWNNror4q+YwzYyTE92HmbCyd3KyDM0w+t2ag94UYpULJMWokf3giath7ap/YHGyKbbFRbib3h5mxMZBZvbS/gfcaOntIq12BIzyTB4Ar29fETn9jb1tw6yDc2G82ro5u0M39HZHVb09228+0ME3mU2itEO88U9U1yiLl0NYyrOBXD2SYPBFmgt87oLA9Sdk7CLg52RXzEXUANyc7YnluTLZFJr89dUgj/s3efAIqUiY/q+YoyrQtamOcPSNm72xhBP95JsKkmkML05H9aio+6Ddnu0iYPAPXZzvu8zjSb8x2YW3K3m/MdkVWvz3fg4vfHuu7KHT6cptU0y5yZZBmZgbWg8QfvrJ1w95TP3nLs2QCYNhi3+4ni8MMW+zO99rtEX2ybLE930EcQsiFXd4BwFA9nbE93wYT9/zYmf1Wi31+zGBsz3ehVD+eDsPxaUUqWew3aVwh4zrPjYtpE9shKmPA5JBXmGprYGyVtpG1mFa7KIuiV2djDeZciW2xY3cSJm/Z4uZ8By7kRqDDtZFHT2Fb3Jzt9FyAfdqeu36a2nYmqG2NUssL2ZZNR3aQv7m/g8cvPyUO8ld2rwu/YFzZuyFOGG7MdkV9Zip5HkGu88VJa9XqPZ/xwseu5p7OZFBoncw6a1tDa5XqbECWuoEnYOz71X5SUmbH2OPEAK4Hehg0aXu+G/jJd7pfPNNWrtGtsrB1yOTd4KoLjaquAr1j8u1agkD3e5Maa/txechCFwPMTN9/nuBiBYVx/T1uIVOJTB61Mx47Nk6tvtfEn491YL9l7z5GjfNKn7dEpWOvjIHSsMYE7N2x+oF2/vBAtzjIskWhijY4VsiVCypg2DQ+3f3YNQNVYlLPIt3tMhTqPlZLoTQqwz3e7Pzn+7zZ67rl3Z69d9wYTX2dLcSzeiWuGSAQYGQm7/oxAiZvQVAwPG8CrbnNx90CIgXbxGQP6wwAxIR9O20Ze8jqJ6bfFn4Pgrqq2xhFxvbdPrNMvmkL30Ytk4/ayEebjetsbLe/bNgv/PVjW81QDfH03tN4YOMBHEY6soP8oFi8hFxc/k0yTole4rdOHO8U1SWtdPO53Wk+kJGB7em5bf7UbTD55ENl+Vp0eSUMPNWjH/TLkTlfwQW46unNSk6/YrE9N2DyHFROobO9pKy+K4OINzLld2xcZux9nZY63z/ICXtn9AJThbkT3MYnPa7cDKRuJWs0CDVt1Ec9IZ8+JN6cYe+LYuIvx+TjsrLI3gH3PFjEe7+qwG4h2A8kvUmLYhRJKdsWGf326yzrDM5uGXiQdGSZ/Hgwwqseerl47IPO3C/qD5+6T9QfOH5eHFJPr8lx5kdC/BgAIvIBXBgEIzD5XJx5hpvxpvkTjo1TvkcAzm3IsWju25LXBty3Ket3Z/znT60cF9torZDjzGto8ZVjjczeHZOXBur0MxdwD8ZKMRLZ+9ZwU9RPjU+IZTo9TtuOGTgp7LPJDGwNt0RdZqyNf7bEj1nmyqbZeSxOuUGIwYmbpE85f/i1QmLsLo5PjjeLrL6UfcZzvLkUBi/nEaOzjF3WO6+YOOUmf9KaAQJhVbRPUKbOwLHhcbEtjg2OiV/Tm2J/YawJ+xUwu6+uc6vy83mQdGRn8gDwzLXnEr4OIly+eTXVAVzZvi6+wS9tXxPZ+/XdbSii3o5EALA/3YfWGhX3kc28YfLJrK/5d+plAocuEl1jf576tysibAuMHUS4vHND7KhXRB24unvDhQiOnoerezdFJn9t/yaY0jrsTvdAhUoGoFk9R1lo1JFeWxdiN3lw2S0ai5l8G5dHpQ/nzdk+FHF0qwnXJjddPj2WClyZbKf2EgYuT264+ELBy8Qy4+rkJmp0oY87/QZIYPI3pjsotYYK6sBwfvvDouzpADCp5xgUZS8fwMVrYUbCrf3y/DR2jcVete88MoL6Gba4OdtGoVWfQ7PFjSZ8hYryuT67ATShAbo6W9yc34Rh24uDY5ldPujHdWcwrk9vQCvV829nMHbmOygEf/hJPYFSutfWgI/RnzJ8vwVnYoewBlM7wzDa8L62BhM7c20XtdG12U0USvXsELV1bcSI2DtbXJ9ed5OPqI2uz2642ExB3WzTFrHOYFyf3QTA6NszgN1qD3M7z25Iv2w6soP8jb1tPHH1QnqA+77hYbq2f1O0vO/O98Xz51Zm8kyO18epDJh8LxGJs7PKAkWh23js7XVNhaEaJB45Fi7KpTRzmtTTtqShP7wz7sR6Y5xsPmtDfafaExHS3Ar1AmALBpt0mzVVdIHfwlQUwVqCwD4BRVCsAybvdM/kjenrBi7oW2VTJu9q3t/gguGMe4DHJA3mgQum5fWQyXMdMXk0fJ87L8dw/9qB1pgZH5s8YPW6wMzME73ssfoun4I0atQtS3e649mlKjAx09Yf3rR66YJ/BTrg9i2t67qNy+NZvW1GtE5vltc3MWq8r7fXiWpZB7m2sAJ7F5i8Ji3aJ5yf/ByVpHOFyvR174Uyqftt4QZl1QRIc4zdtGiQXBuBgtj6zn+nsmjz8fYGn+I6eyYf64mffJN3P0ZNs7cCrHPICOLm+LoN9RBP7T6FhzYfSp6fg6QjO8iPymHGtW4Rk0/WTTbnL8vku4U48XULpROfV0UKShFMFNSs0Fp0NJTCGPv886mrxW2g9J59oqfnruGBs3BAam9GsyGK2EZFs2Kx07VfS8Am0nVb2J6dA929ZMiYbKGeqWL4FdYahpFn8mH+fvBHfD53X0xxPp432+i6nsmHfSxk8qIfu6Q3KReX52DxesKyhvqt+XTYT8KyctCPenWgpo24j/n85uBxnVXA6uOB2l/Zxp1JKM/COjflvN02upUeX5uZ5S0GD5iOLJMfDYb46Ic+Qjz26F0PivojZ2Rr9Ysy/vNn106Ig97acEUcPPxMKU61qUX+XhsjDjaWGeNS/lQ7sbIp6ndvyiz9wWNyJLsHj8l1vjfD6u9alf3w18s1SINEgUJk7MbKawZqKzNWY6176aV0B2vFijhanxhtiWW9a+WU+AI8u3ImzZ+BM+N0zQBzxn+egbVSjmlTCDHd/dgqM/n8BhxSYjj/c9kfXu6rOd68NTgm8ubjA5lDbxSp/YOZxX1zmTm7V2w7aAvHZKO6vMbER6GU3IyHwpoBIsJ6mfrDEwhbgy2xzidGJ0T95PCE2BbHhlJ/4Ux/YWhVHCqTp6wHxfsgPfbYY/zGN77xts5lZjz4hf8/vPfqM9ERwubKeuMr3U+b43XcnO8hfqJXh2PsmzR2+LAYoKZ0hlKQhip0uyTdJ00Kg7JEFSGbQhUoi6Jdtu1TqUsowX9ek4bSKvFi0KRRFDqdHRNhXAwwi/InEDaGKy4kQZAYwOZwDXtNfJFeW5RjzJr9OdvzGRjpAayyyYOoSKEotPhVMyiLhMkrEMqiRLzPvSbt/OdJaAuihMkrKOdmGjN8JgzKErbZUi6s86gYwnCdsPeVYgzD80Qf6xEMqkQfqFJk8mASmby1wKgoerHy3dcBMFQFVFQ3Y4GCFApNia4VJb7ktkFzpe7vI+vsDAo6inbqQicUbciG7nzT7KXa34Td8eQCFqbH2G3j+siwPT92v62iIko49NzUKHXRiw/PYMzqCoXSvXwY3b65sR3CryWIA/nV1oVLLqP4/cZYVGxQ6v5m647ts9gW/q6EdaibXahixu48dgiG6559onOVjFg9d6FRdNCRGICxjH/0yn+IcWbzGykR0ZuY+THp2JHFNTf3d/DMtYtIAxcybk7TAR6A0wXksFdNRb1q9mBMUpn3k5/V8/R8QuJWCbggXkM1QG27AYkBGGKM9SA5v2YDZeWvhRlVIjLZrSZpeQDsVs4OEV9jr5q2A3/I6ueqaqMyhLouVBCjptMHRek2XWnq74tWFKW4lgAKICo6X2PqmPywYfKh7ph8mZwPAGSoZfLoLo1pw8Vjhj8NWX3A8IFpi2xUoBPqtuyhXijd2mpCxj4sSszaeCeuoBZAqRo+bSUmb1DV/Zj7zlBaYmoqh/mocxMtVYFp7V9UXR8pVYHaeHbe6bZ54Xjdh8dWLDN51bL6uo1dEzL5ytZQoR7wQBWdPzNzzIA2brxn7zM7x8ymes3WbSoS+ck7f/i6aQtqZ/eKFKb1rNcWjuErTOt5ry38BKurS6f7Zz9h74E/fLjBd8vkbb+NwhdG2BbWui+P+P6P9AhP7jyNR449jMNIRxbXjAejxELvUy7WfJbhZ/Qsqec8N49nHIA7t9DpsvMy2DovJKN+E4SknJBtCovSorUEoh7N4Lt/yHq4UKb/wghiYge6Nx7689uXhdL9WZ/Xqb+WgMOBI7ya0DAxbQ4H/T5Kop7uj4TXDfUw9XXu6f4a3kjY6sF1uwBk3flEStR7ni3o9L6XB4uTgIR1R3pcNz9opRMHv9DKlzVsi/R8SfdtYSLO7usR6mEESsOdr3zYL2xQNxdOI2jrwBjqyxC3RU7PtUXb39u2i/W4zmE/CnW5v1i2onvlQdORHeSH5QAf84j4dYIX3/UiUX9JhtU/fFKOdXN+87Q4FG4MV0XerEmLxqfaGHEnqdhA65OxBisCk2cwTq0eE39z75YciyZnb3johFznBzKs/mwmps3mYF18sEsqZa5srbhmoDayb7hj8jp531q2WC9XxffwqbEco+bcasreAeD82l0CS21YPVL9pLAPLrO8hywzMKBBOlhyV484GWsSJAf4GWwazdTCNj7gyU9kH3BGljefHJ6UOfRQ5tDOZxyJvi7G93Hx59O2cBOC+KuZm4FPbKPGYyZxiV3A5EeFHDd+Q4hFAxBOjKS2AE6NT4lt4fS0zidHJ0SbjxzTxk3+zq0eXuyaI83k73vdx+Gpq88lx9bGqy2O6Omj1QZf9Ou8MhxjYuaIN60Y6BKW0sFYKw1VqITJF6RRlKWLEx6kUpUoCt26XrX5qxKk0h1sfOzxuDyq+SKIXySOyQ8x5xQhrZYrmNpZwqfXByuYmGlvRs9grBQj56cdtdFAlWCVLkzRpKAlJk/kQhgIqxfLsmxD2vpUUAFSBKZ+PhoapJBMR4jJLQBSUewPJpS6gCXb/yqxwKAYwEBi8sMmBG3IrR0yMzCJn3yhCihwj6W6Gb1CqVTK5Jkw0n327tlrzOSd7ri7Vv24OcYSFAFl0ffJqi1AYJS6HwjNNF8QhVIRb7ZQyn0XqIg3F1SAI/bu1zZY2J7ebUTOEZN3ITyUUj37gePQNQpVRByaUdWOZYf+8Ay3kFBLfvLGuSUWuojayDZxX1L7RMU1StXXXawbt99D2Bbu2XaRY1XC3l3cqpjJc6OHdTbMMBZgNkl/mRvXd0PTi4WzvXzrR3/9UrP5D0gmf2N/GxdvyHHjc37vu/N9kb3vV9NUBFAreRUmKWp9vcOkckweDGPqZCCc8RwjYRXevK4wKgeoI/95A9NsFB7lzoyZSZk80PmGx0e8b3j8m0k9E/FVrYxg/wB0EezlGqSB1pjVgt2i0G5T8CgZVUMHTL7VYTFUJeq6YcIBe1dKJaze/b37NO/ZFRrDtAGHZoKW1dfou9NOMQcRJbo7Ey1L9Uy+1MDcemzQ4aKhLpxRPD5faVHXSqO2FnPbLMShLl7RUJeYVFWf1cNFLZ3VdWJvcHYC07LzDv0EXBlo14N4g7hj79RDCCCgQt2Gq/C/BzyTd18Znskbaxr/80AHuXg9ka6gevaJVicl2i0AtOw9rFvH5Od9HT7O/LxXZ/c1gMBu0Z3fwZYmrr1vo6YfVTC9tvCOFfNmn2J/P62PvBrdZ8PumKtLp48bJv+S44/gMNKRxTWrwzGUyjH5JVl9hk9nj3F+D8aBEDlOK53sVOPOLcQ6aC2HBDjITjHLM/lMW2TNE+nONj7FkQOBxrtGsFsU0SzOp3jmdSsmH5ckfCHG8W5uqZOsh6nH5Lm7vg1+45btq+R8ipi81zUUVNN2FggG7K4tQoYbr5hs2X70lSZuQYk4mG/fPiG1RW5DPj843koP6xDq/Tbq6iDaJ0i354d1cxuWU6InTF4oZ9xGPcYulD9ui57O/RpLOgXXCnXDFpuDFHcdNB3ZQX5QDPCxGSb/wefllWIvPSvrj2Ri3dy9eUbstFuj9cS9EQBKKhIkA7hPTom/S7tLAW55/5oQH8eyxem11OcWAO7L7F/7QSful/WTD4j6Q8dlVp+LaXNsuCkOHoPA8yVMfpf6OFUZJl9bI26gYNliYyDE/gBweiVloABw95rA5Bm4b/2cqN+9elZk9WcysW42hxtpNgwM1EjgzQCDEvTG3GzzZ21ybd+HYt3tJjRMywpufMDTMh0XGbuL4yPy5pGsnxydFL8sN0vZf35crIhMvqBCZPK8gMlbpL7yzrVRjuMjrRkgEDYHW6k9g6mpW9pGZ1bOZG04kn5qfFJk8lvDzaSMzMBQD3B+7Y6fPJgZ5z/31Xj2xuXk2MpwnPiGA8B4MHb4IgLUo3KImU1xR6EKsE4fRE0KuiwSN8qCNPSgTJDNQJfQWiX+84NiAKJ4JxzPv5F0MMfkVbJhOMHx74rTF8xKOcLMVonXzEo5wtzOE31cDJuYPP1Z4ECVgErxjrcThA8os/tqksqkSaEsZCavNMFED3pJ2sXMif3noVyIB+pvKkIMFC2TD+pgGcNygBom1XXZoLBgxhXoMWMtSIMEJk+kUCrd94dnBykGuuizemYYJgyVBqm+u6e1bjs4Femwzm2y1P3vMGvd+oKB6vt618wAFAqFnl2hshYKGlpx7+upMs4fHrC9OlfWbUnIZHt+6bV1fBqRbtiiNq59St3XK2NQRl9t3me8VAXCjzkGoza2XTFOgW6Mm3UXurM3MIDauIiepdL9trAWlTUY6L5eGYPaMnQU+35uDJxfb9oWnsmXgV5bC5tj8gZwawliJu9sGmGdHXIifMef/EasD9LFUrn0Acvkr+7eEI/lGPvE+8NH77WpkTh645oleMVQqTAXeDNphVmV5mWZwcYgQuyY1XOMigFsdGDGcwzVMBn8LVtoJXutKCuHcvB1i9/lUzNzfxHawj0D0Qyp8eeNk9YQ7RNFUYrxfQpNot2CVY2SitRuAYtBWSYhITyrt8LXAkhuCzGuEIA5NazeZph8Ex+nzVE51u9is3R6qQnzNq5Jpw9UNyEIB+2QyYd6oTTYWrDtD/KOvZeYNuw91EtdYNrEsne82aVCATPDia6JURsATZx+/zVmKODKAW+2jUG8is4nuBdP1ZwfAhhjgXnL6ht7AynXRtbjDad7P/lO75h83bZRpzM7hFTVdWpvII2pqZo6d+GuHauvenYL/xJ2TN70zufGPjGP26iZoLiXZdemprn3cX/xu6FV0X22DEcEIn1cjPDE9tP44JOP4jDS88I1RPQpRPQWIrJE9Fh07MuJ6F1E9HYi+jPPr5hpWhuuiMwXyPuwZxn+AiYvhrllFq/tH8I4lUon0fIAYKAHYpmKZtVhWhYlqM21M1XI8fJl7Ra5xAvykjALkcq0USEz+chLwiedaYul7Svo3/8wz5h1t/mEXz9hRtznuJ1s2z7JvWPUrvwMdQXV9q9QD5l87/yYybe5k6jHx3Jx+hfx5laPWD1nzve6N776v4ftpNu26I4oqHbLwFAPI1n2mDz6K2bDjT/8TIejOkjnL6qz6rURRL3XXyJdOhbqxhocG6Uo56Dp+TL5PwTwVwD8WigS0UsAfDqAlwL4swD+NVHGUnnAVBYlPuaRjxSPfcjdHyTr5+QVZI9mYtrcu3WX2BmOjTZErlyqInGrBJzngczkK3HgqU2NtUEaj9uyxZk1eT/S+46dF/VHMuw9pz904j5RvyfD5I+PNsUvi4EatDObMFl2Lm7RKOh2tHegupfq5jM+Tq1xSuDQuTg7966fFdn7/Rt3y/raeVG/KxPrZmu4JfBmYKjHid87Ny+E2MXUM3mTYfJEqZ+8YYuRlpn8xmBD5MQnRJYO3DWW1wzkOLTsGw4cGx4T22hFr4ptVKhStE9YtkEUSa9zuw9u0hbwTD66Nvl9bRElwrHhMVE/PT4t1u3c6lnRPnEuZ8MRYiY5u8hmpr8McX7t/cxPnoh+BcAXMfMbm39/OQAw8z9p/v1fAHwtM/+vRfksw+SttTj7mlfh0va15NhoMMRUQAKjwVBEM4OiRIU6uRFaaUClTF4phUJi8qqALguRyZMQhthvDJIweV2CSfJJd/FdLPd5s4sTU8JYkzL2coi5raJVo5zRgVExQM11kk+hmtgyQC8fpRRKXaRtRAqDiL0zeyZfZJl8bHwrqOjaIiiThgJIAarfFrDkbAQCkx8UZRuiuD3dWgy0pHOjp6xeKxeHJmSsLrogNfFXguIww7KLUROyd6cThoVj4OHnOrOL+RIzeWbH3BVxj0MzOzhQxEzeMkCEQlGfyXvOTf19ZOfWNrNmk3BoN8u2CYdmaBBFOrsl+5qQ2C1qy40nFXptMbeMknSvzu58NPF3+m1hLKFmiyJqC2sINQzKoC0Ybi2B8/eP24JRs2fynT43BtwAq7gtmBUs1z17Q2UtrCVYrjGI9NoCYItS91n9vHZ9NO0vhO/+U9+MrWG6WCqX3hdM/jyA3wz+/XSjHVq6sb+NG/s74jFpgG914Yt9Lsw4AWfsE/m3yjB5knmzgQWb1BNgVs8xLAbJ18KsnmNUDpOvAosaA2qiO0a/0Ua1rLKXV8vko2s0L5zER99Uot3Cxb9OqoYCWmbyGfaei+/jmXw84zWYY1gOki8nC+sCnZn4A9itY3D16uu+PeN6VCTrc7iXYGwzaTYvbRYhhUzeBa2rIn2gih57b3VdtPch1Au/b2jM5Ikw1GXri9+/hvOHT/JSGnOTnq/gvprmET/mhslXke5j3cQ6NZuzxLpj2EBl4uuGTL6vzzi1T+jAT77/wnPXmNr+CxLwTL7fFp7Jx23k48y7ftTX/YQjbiP/hRrX2W8ulLRdg6uqOmprALZZSxDqK8UY791+Gi8/9RIcRrolriGiXyKiPxT+++TDKAARfTYRvZGI3nj5cuopk0vro1WRcwPoRYHr64fH5HPXjnekARyTlzj0MMfktbxdXI5DA3n/9hwvz8aszzJ5mWhbsJgXQ47jozJMfpBh8oMomqBPuXUKlCknkL/Py7dR9/f4yybP5HWiA10bhbom3USDjJg8RWsGMuXsBk4SdWD5NQO3y5ul1G+Lbs1AfL4W7BME1dp2Qr1UhRj3qYj6RfgSCS/Y5Z8v67K2mqwNZ0EbSf2ltrULzXFI6ZaDPDP/KWb+YOG/n1rws2cAhEFQ7m40Kf/vZubHmPmxU6fkfUql5Ji87Cf/sgyTf9nd8gqyR8/KMW3u2TorMvnj4w1x9jpUhZvhRz+pjRH95+emFl8itamxPkx9wA1b3LUu8+b7j8sfSo+ektn7o6fl+D45v/p7N86KD/GJ8ZYYZ2Wk0y8RsBsIEyYP1xYJk2fH6ktdJucbNtgcCvFRPJOPC8vAfRsyY3/Rxj2inmP1d6/eJeonhsfTlzMDY70i8mbAhcaIeXPNdRK1lJlR2bqNehgmwxZDne53a5mxmWHyp8Yykz+7InNlp3Oi37Uic+sTYmx1YK1cTb6OmRmlyOS5CZGQMvm6CTWctp0RY+szGKtlyuQZhOMjmcnftSrbIe5ZOyfrq7J+NpPPidExkcmPixHOr8mxqA6SXigm/1IAPwTgFQDOAXgDgIeZhRVEQVqWyZ/+7Ffg6t6N/gEGhsOhiASG5cB9HkfjaqkLVByuUXRJKwVoJfLmoiiSLfEK0iiKAtMwhgy7pejQGjPvtujLUwwcL4z49ECVLSoKZ9aKFJR2n/JhInJ+2Ib7+5QaZqw02Cdm7D6uTMzqPZ+Oz9c+/knP1xtQqokVEz1sOSZfKIVSl8n+uAMqAd3EQwmuXVABkPt07rUF3C5JNvaTt4DSGkzcZ+mGMSiKBUze9L6qrLUoden2NdWh3ny5KAaFcVZsEz9dqd7Uyb24XBuBOrTgCJDzIyfqc2iw9wu3Pa4MJoCox4+ZufWTD9k7A21YjCKIQc9YwOSNX3jWZ+zzxq8eMAlvZnZ9Yqj7/vO2WcIfxtkx1qIyLva8Dnz9ne5ZfXc/WyYPglL9sMu2cRmmmMlbahbQdV/uDOc/X1mLMmqLee1WZmhCry1mxoBZgQQ/ebYEC9NrI98Wzke/b58wlgF2Pv0+GWtRWfflF9onDLtvne//09+CE+NjuN30gjF5IvrLAL4DwCkAP0dEv8vMf4aZ30JEPwrgrXCBPl57qwF+2XRjfxs70z3xG1GOHwMXS0XgzVUmGiQTxF2MlM4w+UKJ9oBaOSYfX3dWzzEoB4gnwjNTYViUCXs3sBipQeJvz8yoyTmwxV8efjYdv8s7vX+gtrXzM0/atdm6LNKzTF5p+T4o1doJwjSnGqXSid+7Z/KOi/fbYlCUYCsspedmf82Yvd+SyffvddXEb/Fxc3xqbJBAvM2jbrZ+jPIvdCnafcog/nzvfIUm1ni/XoSG79dpnxzoUmTyJSnMjcW88aFvEY9iGE75MWd4M0HmzS4+PrCfYfJ1ZSK9Y/JhCtl7X9eYc3q+/028qtpNMBr7RLCinNktVpsZ08Ty6TAPgZwhNKiD85N3z2wV98nm34k/fPPQmGglu49dU9f9fCwIlg2q6HauFkM8sXNhqUF+UXpeLpTM/JPMfDczD5n5DDP/meDYNzDzi5j5EWb+hedf1H7aGK+J/BuQ/bMBiDHdgTx7BUjmzYfE5EfFUOTBA13KTD5jayDk/dtz7Dq3xkBl9Bztdkw+/Y37BBf85KEw1GkbuRWhaVsMhXYD8m2hiBasGViOvS/i/lLyXzCpbhNW7HPK7j8g6KUqRJuJu5fBl1Tz5yIm/8Lz5q6/9Agcd0y+VzZG6w8fJhfrSFpXUYptnXv243sZMnkSdGBRf+n74nd6zoaT0VnuS5WtcdfK7aPrW6UjG7um0AU++hF5j9cPvSfD5M/L+iI/eYk335LJR8lITJ7dTD4xBrLzn9+QmLw1IpNnAA8cl/dyfSTD5F98Ksfk7xP1ezfPiS+ek+NjYgz9cTFKeDPQsGUT2S1YZvLMjLmp3EAfM3lrRBczy4yzq3Ls+/s3BLsFAy/avE9k7A9u5Zi8HAPn1OhYugKXgVU9Rh35eoMdcqolJm8NapHJN6tKI722BqMMk98abcm+2xn/9ly8nnvWzmd4s+xXf3J0PPmqZAbWBmsCk3frKmpO28I2/D3WK1uJTN7hkjJtC1isl2syGx+nMfEBwrnM2oB7M21xz7qsn8vZLUZbyfMBdt415zL7NxwkHdnYNcYanPysx3Bjfzs5Niiaz+PoFV0Wpftkj/RCa9Q9G71LihQojqfCjZ+8Tpl8qQqookhwxFAPQJowtbE+BJPF3JpemQZUwqp0M2dNGlqrhNWD0e6zGjN2734Yeq4YazEeDFGZlNWX2u0hG/JpZm6YfH9PVc/kC+VikMf6UJc9FMZws/tCFz07gXvIXTx556NPrV6qHJMnt8VaxORhnZtmzOrZMsrChU0glequziGrd4zYgiMm33jKEENFrN75sGuQ5p7uYroXbp9a6q6rmKC1hiUT1Jmh2N1nA9Px5mYG7M6zPR1w8YxUwJUdh2aQIsebA27tYtEoUHK+gWrqVkasXjdx5mM/eRfs353vW89xaEApYBDolXVx5mM/ecOMumHyRF3dDLs4Li5sge0zeVbNVyS3dgjLDMPK2VECneFCLNRsne9+HOsGfSbvWH3trk3c+sMzXFgHhgJzyuQtO/zi9pEN6mwcUnTx/n3b+bg8BFLd16e1Foo1/tOf/6cu0N5tpg/I2DU393ewP5f3L537WXbMoU0tMvl4s+k2KYiRFBWlC5sAABk/+ZoN2KTXnfEMZVm60Sw4NsfcrRiNylU3Gw/46Hxhsla1LDFM/ovDRpywql1bxO94Y53TrhQTxm2Y0NcU0s1QABeRU+LNDNlmMucaJemEvc9sjskDg0IlOgAQd/tu9uvWRHE0fd2vRjaCLjN51w62jpl8sygsao5CKRdrPGpSpdIFdYBbQGSNFZi8xUCnm6O7gGIFKmsSflyyRsWM+CrU8ObaxpzYiDo3lYp5s2f1qU5uIIvyKUi1awl6dfZMHrGu21gx/fK4zTzi+gIGBel2LUF7PjuDeWVtry28n3zF3MurjX4pMHnvPx/vo+D7l6n75fV+8ibuR4Bow1kth3hy+8JSg/yidGRxzebKOsaDkXgsx8tzeqG0QMZcR80xeYmxE9yK0eS6qsCgkJm8xOSGGSaf4+jxTvO38xspvj2QZ90AiTyYweJvbKaNFEi0Wyxm8rffFoqUWE5ggX0io0u8HIDIlAHXFiJvBmevMRB4c57Jl5AMDi5vwYaTY8G4Pd58O3quLdx1BBsO52w4cnsTSGyjQYbJl1RIvhgu1IGgqwUrK3Ltl7uXC/ukeERuv8rUuHv98MIaHNlBXiuNP/Hwh4vHXnaP7A//IRlW/8iZB8QOcO+xszKTX9lyM7DoRwNVOg+e+EvB1qhNGq9lVs3TjsRu1enGMA0zWluDuzakWOaMB0/ITP7Fp+U1AC/J+cmfvF/U79tM/aQB4ORoK/FjBlzI4rmpZCZf93VewORnpnKGWoFDb40kJm9xdlU2Wt23KceNf9HWvaL+wKbA5LGAyQ+PpyukGVgrVrJMvrJxW4RMPuXQOSa/UowFu4XF8eFWojsmL3PiPJOXfcDPZ5j8qfEJkclvDNaFujkmL60ZcCEPbp/Juy3+SoHVy0zewu0JLDH53N6/96/L6y3uy8Q6ulvwkwc7e5bUX1bLcTb+0kHSkWbyx//Oh2N7spscK4syid0OAGU5EAOIKa3ERU9KKZCiZKBXSkE3S9h7+asCqtApky8GIFIuvK9/pTMwLIZgYveZGrzqSyphNSfX1c0s1fkHdzqz87eOubW1LkyvyOqLgdtgIWL1A124ePWCnzwR9crpY9c4dm2T80vtYtRQoLtokxpVEDfHL4ZRmlDbPpN3nBa9ujGj4d8yk1fahXiI/eELVYCpb28wxrliGrZQwWab1jgPqtR/vvGTT5i89z1XTShias8nuHqE8XR8PlppGPSZvNsZSrn1Ck0d2DrdrQ3o7AfMDLKuPEzo6bAAGm8jb1dgdmiFlIYK2DsDmBsfLdOiaLiyZ/XkY9dEvBmNn3yoOz9514cKQXf17jzCamthuNnViWzLzGtmsKWmTp3uWH1T7uB8H4vd7WHbMXbLLv68ASfsvWpCgMf2DM/kEbVRZWwT0cb2vPVqYwDrxoqi6PcvR33S/gKPB6P+oknhh//CP8fZJYyvH7BMflrNxGPSAN/qwueuNMADAAjiTF4pJXrXZJm8Nc4wGeHjWe2ZfF+vqEKpUyZv2HUsG814ARfUyMWuiTm0yyNh9VZm9aZ5IUhNYkUmT2KEzZzdwoJR+5dgcJHa1igEJj+3VmTyFrdg8uCEsRt29yD2h2+ZfJ22ERElsewt5Zh8MyGImqMQJgSAexGKth2FjskHeRkYt3DLWvRvm3UL0qxJ1gA4X3JObEuKLWpON6n3F4y5MmV07yc/S8538dfj8zUpVIJvuCadZfLis9ZEmzTWJIxdk8bcmGbjj+ja1ibs3UUD7dsP/FcEGGncJPYxbSL27mPXRJXzz17M5Nu8InmlHOGpneeWGuQXpSOLazZX1rE6TMPxArKvOuB2YpKSY/Lp4O9Wown80Mq8GWCMimGaP+mlmPxAD/JMXnhJtSsthTTIsHe5/Hk/YwKJDJQhx9a3bEV/eEVKXkuQZfL5eyYlN5tezj6RY+8SC/bXkDluJo4Py2sGctfQpMX7MMgw+XLB/gNpdJb8vQTy9z9X/sUxhAQbTpbJZ9ZV0CImL7VFuTSTP6z+sqiN5PFF5v5zU8nuvgdMR3aQ10rjlQ+9XDz20rvluPEvzez9+shdD4iD6j3H7nJv8Whae2J1s2HyfX2omq3/IvZurHFMHn19Xs/TwZkZVT3HZuwn37DXs4KfvGU+gJ+8zOofOnmvqN+7dVaMyHlifCxhrACwUowwM/OEQ1trBSbPHZMPktPnyYIoz6E3M0z+rrVTMjPdlB+cF0nsPcfkGbhnTbBPsPMNd3Hg+3VeLVbTNQOLmHyzJ3CeyUe6qbGi+0zerfi1OD7cFJg84+yK3Eb3rMgx9+9dle0ZOV/y0yvyHrKb5boYi2aoBt2eAoHOlkVWP2/aLbVn1C4sSNQWxhpslNL+A8CpUbonMIFwTyZG0f3rcr/I7RV8fu3Mcv2lXMHp1cPxrAGOMJOvTY3jf/PDsDPb609WLFCUpejWVxSF+zyN/OSVcj63sU4Nk+/FOG/O10on+5cWqoDWDZPvsfcBSFESy35UDGBVui1dqUtYaph8UCYFAlSzDVmgExOKQjX7jAZNwcDA+4ZHrH5YDlDbfqx0ZnYRMNlCByzVWAvdeCCFEw/P5LXSvZek96svVBH5vXO7A1Rl+7rbJYuc3jUdStIOXwR1YPgZWOMf3nJo9z9KNczcs9RmwPPRPVNWX8KiH6PGGIfGEpbacGvb5N/5N3NrJ+gzVgsFBa2UyORd3+vrRbM/ac2mx9j9rkc2sCvYxt8eRI5P66DOJt8Wilx8H+25csPqdWPn6DP8ps5wuq+zMY7JM9jFqPEsvXYDl7cFeN3x6WZVsO6YvHExHJoydfYG0/iMK0Iv5pD3JQcAQ51dyVoGcbOtH3Gv7dg0W0L26sZtHWK9rt3+q9L9h6W2v/TaoglFEfcjH3ojaQvTbH2o0dMLKvAjf/nbcPfG7Qcp+4Bk8tv7u92sOXpPSQM80BhHBD/51h85fptn/MVVM0jFqRdPPsgr5yc/recoB+mnaGUqFys9Yu8W3ISg7esuZk3TFNE1Wt/w6IBnsXGsdNv4ycc+wH4YjyfzCkrcJUsh5M3BC4C4858PylSZyrm/cf+bas6128uVuXe+BTDQGsw2Ye+q4akxe/dfIiKrBxJ/eNP4yccs1TYx1OP9f0krkckrjSyTl3izbbzP4xj+pmHvxtiencAADatP66zg7nGsM9zOU/GzYhr/+fh8vyAwbiN/b+cRh/Z7osZ2C03K7YYVryUgZ4xvK9QkIsLMlyWyQzgm328LZpdXarfo7BNh3RyT9y+/vu57YnL/WW4L/4yluo9dE7ed2wAo7hqDQYlndi4tNcgvSkcW12ytbmBjLO9mPipTLr5IH2g5JogCiVzO7zKUHmDxGgXlYtcMMn7yGSavtVxOUll+mLNPZO0WGVbfeo5Eya9ijZNl2+581U+yn3yxgMnn1wxIHJrEWCGEBew9U+ecPSDHWBkyY80z+cxagsz9dGxaZvK5/Qfy6z8Oxz6xaC2BdHU/CKc6y9dYkskPdN5PPlvOjH1i2f6S03XW5ieXa15XeNGxexL9oOnIDvJKKTz2og8Rj73knMzeHz0r8+mHTt/n3sLRbPf81hk34430YysbTdz4vj7QJWbVPNGNMaKf/Lyq3KDQ0x2HTmLXMKM2Bnetn0h0yxb3Hz8n1u2DTt4n6o+cul/UX5Rh++c3zsByurvV8fGWyJvHxQiziL27ekp+8oyqrmFtf2cr5yc/xyDjJ7856seT9x4RZ9fk+D73bMgLTB7ckuuc86u/e/0MEk8mBk6PjomMda1YSXizY/IalZH85GUmP7cVFGeYfJEyeWMtjg02E52ZcdeKHHM/x6HvWZP18yunRf3U+Hjixw4GNgdrYt0GatBy9l5ZFzB511/SthioItGNNVgv19K2AOP0+LjI5M9n1kPk2Huu7c6tnknsB2Dg5FhYY8LOT/7kyjEcVjqyTL6qKxz7Gy/H3mw/OaaLwiGEiLHrQrfL9kO99ZNPmLxj4HEb+bjuvU2YuZmNqnR59lAHTD5k9boEa+r7yXPj4aCQ2A+omaVaWMTwXRVF4iffi8sSs/pmCX4cQ927aMa82XNPrbqVr84fWTWxXCLerLSLCRTsO8vc7AClXbjZJHaNdgG7whglfl9bw2EcF+8nT03Mlo6lovk3Az1Wb61FWRRg2IBP+7jxru3itnC2Btu0RcdMtfL5p0yeCEDQRouYfNEwecN9f/iCnBdVuL+sDfzkQ999ZoZi1TL5UGfL7dqGmEM7+46LOkoU6KQafpzTVetia4xxM9TWPtG1kX+h6yLSbcOfVXd/jLFQTQz9Pnt3dQZiJs/Q7ILfxLpiV06LflvAoLW7+bawll0dnJGntU84Jt+FI47bgphaW1C/bhzYLajpXwZs+ro/36NPimw4BRX4sb/6L3Df1u172HxAMvmdyV7WH75lxDGfNtFmik1q/eQTJk/yIimtRN9wIhL9nntMPshuZiq3VN3G59ddeNVo5tHs6ZB8LfgXUeIn3/SkOAa959CxzcHXN7VFOCBQRxkVRK3/eVjYkqjlzf2JCgf7y3Z6ZQ0KavhxkI9FjaEauEElZPLM7tOcrbCpi38g+2XtmHy6/kDUPZOPdJ9r3EZtf4muSzpcSxAcU10bhX2gbgZmy32u7GLoFwljd6y+FNm7JgVrbcrkicEWsMk+wp435/So3zfPU9xGzmiZxlDX1Nktwrop1e3rENqJiDiILdRvCxfXJm2jQumkLZzrpos/H94397wQ2HKvLULbUNIWTT/M9a+EyXu2b+JnikQ/+cGgxIXdy0sN8ovSkcU1W6sbOLa6KR5bGYxlfSjHuhlqme/l2ChbFrkyM2NUpromLfqruzg3EpPPx64R2R6pbFyeccYOMS7ltnB2gjRp0hneLLNLy3IbAZmY+wv95NO2yLLgjO2AQFlmKq1tAFxoBinFe62GV1nGB5wzTL5QJHLrnI0g758t+/MT8v7wuWuMMnqO1WtSmTg+8rVzTD4fu6aQmbwQZhhY3EbSGhNCnrHn+ktOL5XcXxTyTP6Djt8v5nWQdGQHeaUUXv7Ai/tic3MfPnu/+JsPuqvRo17w4Kl7E+8NADi/dXoBk0/Z+0AVmFdVPHWFtSbd15QZ87oSYtc4fX0g+8mfWUuZPLPFfcfOJjqQZ+wvOiEbdu4/fr7x0onaYvMMjMjkN0UmP9QDzOq+n7xj5pCZvKndDjqRPjNz95BHXzSVrbE5TBmrZYszEZP37PXujTOJDjTsXUie4cd1Prd6UmSsJ8bHZMaqx6IPuPdASv3kjRhDfW7mDkcIbTfWQ4FDWxwbbghMHrhrfDLRwchy6HONHtfNs/2UNx9PbTgMbJZrqIW2KKkU28JmmHxl5f4yD7y0em1qDdbLVaE/Mk6OjqdtAcK5yN7g2+juNbmN7m7sFnGdz66eyvSX40J/YayWY5xY2cJhpSPN5Lc+82XYn0+TY0prceGOUi4Wu8jkJU8NpZqPVO4xcxfTJtr7ld1slBT1lzuz21sWKg15MChKQKWIR6sCpCMmz80qwoZDt69nhvMxLhu7gofvbqR2rF6IoV60ceOpJVhed7Fx+v7z3heePHdudKX8zLm/N2eh3MzfBLHPbbNSWGvvJ+/z6eLJx/7zhdIufpCN4rWgY/I60K0FyqKxKwR+9cyN/zK7vVl9nZ2ffN8O4d73HYeW9dQHXDWYJWa1nT88I/T11uT84WMf8IKct1fFdceVA7/6nv98w6GJ0Pruo3H9Zev6qmfpvo08YweFdXN4SxO1MXD8puHWWhBiJt/w6ZbJO9tJyOol33BiOHsSAX2/9zAuj2rPL8j1uzDef8jqazYBY3e8nODaNNR9eGluykoUMHnIdgvvPx/ef9v49Hv/+baNGjuEZdvT2XIQuyaw1ZjmGLuY//75tMbtOfyjn/ov8aLjt+9h8wHJ5Hen+yIXByAO8D09Gs+lAb5NAq9XJDN5EHWx7IPkvgbS0+eeyUfJ2IbJc/+67I2MjSFJKmYM332dOdJbVp+cz+1DH12huU7fA0ZBBe3KgU4dYw0qweD2pdZn8i52TTzjqUyNIQ0Exs5QDZOvIzbqf55GQWzqcJt2CG/IjnUOF8AEyTP5mNWSomBdQqdrhXZBXc8fntxExHuFhOUZkNvX1gT3zW0qrRv2Hvd91RoS+4mbATmus5utxnXgWzJ502PjlGkj5Zk89+usqFtLEP7GNovCgP75lpsFWpzWuRDawjN5G7H38HhYN9clfFvEsekbG07GTz5eP9E+miZ6oOGeB7bcs8mUqsClvatLDfKL0pHFNVurGzi5cVw8tjaSY9rkYt2Mi6Hof65zTJ5Z9AFnZjGevCYlxm/P8c9hkfcNVxI/JJL99gGsZmLurw3ktlgp5Xg6hSqy8bUle4BhFuvHyDN56R7k/OTdfqdpWbTS4j0jUJaZrmTsEyulbNsZqDJrG5EYq2V5T+Ack1ckr8/I7XdbqkLk0G43L4nJE8pMXjk7xEqRs2fJaz3cV8rtx8dh5PzkZTae6y+LmbzExeW2Jsh7QwD5/jLOtNFAy/1FkxJtL3NT4dGTsrv3QdKRHeSJKI1R09zdB07d8/9v78vDbSuqO39rT2ef6c7DmwceDx5vYBBBRASJCoh0azROoRMSk4h2YmM6nbRDzNRfTH8Z/IyaQGubqF+byXQbE0B5xDhAIlESDcigoLQKjQwRZXjvDmfv1X9U1d61q1bdd+/lKXnHs74P7n2/s+/eVWvXqV37t1b9VuPf5udxcwZv/tm26U2ab2t+sHF8Vubk876fJ8+6GpKHs9auKTx8ebCsXjsbC1TG0mAJvazt4YOywGzPyZ/VK9+t7u4444vJzaIvdkxt1v/kxs9tk6qWq8dDj82KnPykrnfr8o15nHnaNer1tcSSpF1T5ck38ZqT9znZsczn5IuywGzH51gZjI29ps68udZmzdW7vnBxYxu6M76PWNXsLATd+E6cY3mwBk5eax35vliuql41fFEMkMctZ3gpymE8kzh5xpzEQzOwoa00bRq+YKi6uRLf3JHjEzO63q3b1n7SDXDySYCTl7WOBoEYzpE4eT+Gw5jOJ3xfAJgPxC1cX9T4rOgjo+Pj4lPtcXHPQCdpY7I9jqNlxywnv7S8hPEfPYCFpUWBY4/D3LvQX6PRrsRZNFgarl7/sWUxRYDRmbc48zROAKJmHVnWK9eIPE2bNEpBsebkLY49RqRzxou6TQyACXESBeWPS7sPjCNw9bG4N8BwuCbnm6GoniZOGld6MK7GfcWlU7Mebcml2qkYGV2WmjNN46Ta5m9zr0mk6p2aUnzm/EohtKn7wvoLVmnOkMU3V7n+Dg9t4hBs9U1/+eNAn+08fJMPbbh6QHOsFa75Y0JD770oykqjxs31TiipUgobPHSVV2/7Tum1RHo8RrpeaKl1LqJItbtqU2m0a1SbyObebY7d8MpapjjSnLzhlZnNHoBa38XmoRXTofAoiurzs+K/VZ68WuAUZVFx7FVcgeu4BQCnz2XlU3ccxXrPQMFFpTmjeHE1ZkqUKiZDCq/jCqj7Vpo8eYPbvlB59SU34xMqBmLiEE0fVZx8FOm4hQoo25y8jbfiFB9+1TuDBXwkG0pO/onFwzXn5szbQe69ZOmNraYDfMpMAkCRxck3PqZao8TCi0Ce/LLh5BmObrgqzuC3J/xArp5dni/Y/CLjzvGsJ35fs976O+tBSWTx/tbxRFGti9JYwHCl823nMQ+4ACVmBVr/wXIxQBxl3kqIoTZKcemvtE1b3Idhla8s1P+U+mYeUC5ufnO17In0FgZX4z6q9xI0Ls3AcsU3W9eNyqpfti8GXCCteGVXoybxjld4BObSq2sLzcm7OBse2slvNwXcXV65xv08eYPbn0VUa+43/iay9JTs4xvxDItj13nybtxigAKJXskXDscemwevl6+ufT0I+cLVrtHjyOPeDe5+cblqvxsyMZy8jSdRgoef+PaaJvmV7Jilaya6Y5ifkMu8hTRtQngna4vcZRIl4ZzeACcv8c0RRY0qMsZW0o+RJq4kljn5iEjMzycAfVceQZsnm6Ct12qLGeBZnMqcPJOcJ1+Wwb0EonZNJOvyhPYMqNxj36JAfVQCBbnUUHyi56axasuTTNRlCWsIsciBhzj5kN674qxlX0gv5HEUifeMIOsHAWpLvWS9AB7iodNI3vcArJAnL8UJWM7FT6KkuetbmytpULdHHl9xYI8JIcy9h8ZLyHetpCWOlziSx8tSsYx987Jc+nrsmJ3kiQgnbtLBCWcVu21a54w7PPSW6Q0yPjkvrlzn+lMo3Px2qIfFssTJR5Hahes8xFfi5COQjxfL6AqcfFEUvqaFphY29We9vjGArePzYp+3jsn45v6cZnRYw+rnfG/a5+QZmMj7Yh5zprX1PY4VLHPyxUBt9RZ46ETgWJfLAfqZlPdcYqY9AZczZTDmtUa337cZGZeOZ2CmO4USPsc6kcscax61RI0aAol4URai5vpSsazS91xfFAO040w4T4kxp66peSMytV/dvs3peIaLz7anA7h8/GQ+LnDmQD+R8+QTJOK4CNZ4DXLyAySIPXxQFugmHfH8k1pzv3n/CTPtSdlHAV/MdWV8Op/0xgszY7I1Jo6XdpKLNZ7Xa8csJ7+4vIjxVx5QgmCORYLmOljRLBX3LuHueQL584nm/F1+P01SMPkl1VoaXxo4GjVG66Zo4nEUIYp12TMLB1mcvNsHOkLfTJEghhVvYMXVVziDYsXhmlxmPeODDLdr5ckbztxc0+RJl1zq7Ac0OHk1+WfVq3lTZz5WchGF1q7RqYsx1Ro4kW6r4aGTOEaJAmRx72WpuFLmWk/FcKYm17+ZA17vAYgt3HD45vjIOt68TblcvYlNRDqrpeL2Na9szsV6IjBvL4ZLN3iqfWr6bNI4Y8RVOqaNExPiOEbBTV80YinV8TYnX/vI4G5+u6KAmnsGDH9cFPW9jRxe2ejM1ziqfHvocWrjinpBpSevcK3jA1b1bq3xojh52xf1uFDl/ArEsQpQVXsGqKk508yTb/oihJdlraFfc++aJmN4PlqJkwdbefJWbCdPcnz4R9+Bk+Z2eXNPyIaSkz+8uODleBs7Et/s4UGuWzhY46G6sFL+fGMFbP3ZoBggjTIPL8oSkVXMweqA3WjZQn1zuMDq4eX4UC9i1C+Nj5Qv1ErQQqkudF5a/C6BRC17Jq606m284BJUag1yCy95gFhvViqdPPmk4uRdPlg/s1zuveJGGXbnaq6+iZs9A66ejv257YzIeqDYVyYdCAT8/PNal8W+bllN+A0cA5Unz/6egRhxNUm4JnH16l6uxCuHfOfzygrnpuZMAI8IVZ2ExveX7Bq8zfFSlIW+un3PtFCc0zdmXbCdGYOB22cVn/B9gbAvVvSRPI58H8m4eaC6MZyEYnz78He9Nq7Xjlm6Zrw7hi0zsqj+ZHdMxCcCeK/VkTnQOBbz2wGsiZMP1YrNAvnwWSzXqUziuK7kY1lEhE4gH3683Rfxybbsi/G8K+b0tpNMzAEnyHnyzCzmdYfq44b0vrMQJx/HMi8eRUqa2GsnBbnU8cCr8Xgu4+00F6+dkJyjDwTy5DnEyctcehaF4xOSKV0W4Z4RBbWLgjGcTPZFN+2IsZFWIIZDoCAnL20MZAYSkZNfYbwIC7CEwpy8tOclovB4CVEpIbydtuXxEsXiuF8slnDyhhPFc63HjtlJnoiwfVbOAd9gArLmXuufG8Zn5OPHplWGiINP9SZQ2Fy6XoH0UpX33Fxl62rwQp48m+o1zvGDipNv4kofvOWt4ouixGTeryiU6vzMmO1Oyn0zWjerxOe703oF08RnO5MoufDw8VZf5lgp8eq2Gs5c0nEJ5T0vF0rfx8cLtNMchmOuzl+WGG/1qpWwwRlc6YG43Pt0e1LEp3SussuxTrXGGxyrebvpO1rppm1ZoH4pGOE8+YCefIiTNyt8Gy/KEj3NQ7ttGs/UOJL65uEAJvMx8fjJVt/3NTPGsl5DW9/gnSQXx0sc0NYvucRAwJeLgafdb3BpvAy4QDsWxgsY/VTqAyrdH3+8BMZFAJ/MxwLjpdvYV2Halset4ANmPXbMcvILSwsYf8V+xXM7VvHTHm7lwxtj1Dy0Y0rrhp1cdVUrtoC/MUjpkrN67bSEWbIk05IHzbaqbJnYwyNSHKvHycPh5K02NTh5G684dqfP0SqOl3D7XKx2mJbqpbPJvesVeCOP3aziIzUBeZy89oWtURPp7IflQmvaaL6ZAMRxgrJUedIEVFRPEscAG50Ydc/LkpHGMRh1XMF8sSOLSzdWlrWGijkPIHPyBq84eRvXDzzj64ZueKwyhJTev9lLoHxhuPcGD60VNgflAKb6VeWLKG7ouzTjDWXtC42bPtu1XBVXr+6z2+eq3qmj41Jz8s16p2bo2Do+lQYOs6Ohrjh2BoOJG3gSqVRiW6+n4uqtuEXlC450fVyjgaQ3j3GdJ2/rzDf1d5q+sPtg+8KONzRxDhyP6vthjxd1v0trvKj7305yfPg/vAMHNpyA1dpQcvILS4veJGssxLE38uEbeMj0bWbhbwLXrrXs689LltVxBkWJTHjtLJlRoY0/5Oav0knZP8w/j+ML+3hGY3KvzYpPcBN2VzoGN/noHl4U+mWkxgsuQVxv8LHxJFYPEbMyMk3Q1W69vHdUeLMT5jc3rmBzqY1zMFcBN8nZQU7ewVXxdf0m17xwpcviauVTaeIfNT7gAlkU6YeTnTMOxDCTlt9O5nDudlgr3/Wd7KOqzR5uT2a2j1Q8SDXf9pG+j9zEzRsJw98zEEex0OcCESXWyrwem2qFL+m6Q8SP7At5fPnjpW6DfT+ruIUzXuIowqMLj3ttXK89KbqGiH6HiO4koluI6CNENGF99iYiupuIvkxEFz7pljo20RvHjnlZwGe6PynjvQkRH2/3ZC3zJA3qWoicPPSq3f9A5N4yvfKXcMmSOEaaCDm9RKJeDwGY6oyL55rpToj4ZGdM5OS7aS5yrDFFslY6y3yzWlELnDwi8ZkV4lhDedhxFIkaNQQKcqahLeRTAbyXdcR9FWmUhDn5NfgolCefRvJ4CWnDh+4NEdAJxSfytcVw+lkvEMNpyeMlkvcxAHJsgVnm30M1ftMoEZMikigR95jEFIl1mVceL4GYX8B3vawrj5dYHi9Lg2WcuukkD1+vPVlO/noA+5n5ZABfAfAmACCivQBeCWAfgIsA/CFRILL2JGzT5Jz6pVqiqR/VZG5utv45ZXRfnOOnuhNqZeTg452+Xl00z9NOWqpajI2zKkkm1XKFTslqruIUDx1x/fcGXy6KWq/FwouiRD/tNFeEhvczAVYLZ6DWpXZ9YSYwB5/O+/A4eVbcqxe3YKCftJv6PrrNCWKRhw5x8stlEdQiifQrvv3GMCgLJV7m4EVZoqfz511Ovp8p7RL3zcMEXl0udSwXcAD9tNeoX2qu1U06jfx2g6eUNvThDc6AGLcYlIXso2JQbYVv+mJQ7SVo+IJLtBOfh2YG+mnHWu3afZN91E8DPko79VuW1bdO2m5klRk8o0z0BTGJ+y0KDvsixMmbOrjNcVGgFaUerjh5IW6xwngZk3CuA6/euMg64njp6O+O56M4DW7EWo8dNU6eiH4YwI8w86VE9CYAYObf0p9dB+DXmPmzK51jLZz84cUFTLx8v5d7DkDndPv9MjnNHr4iJ++/2iWxqptacdr67SxLMxRceNvps0Tx04pjrxsbRzGSRG0asvtApHLAfU6+5l4bevIMUOzkyTd8gVVz9dW/CXVbG/nzTVzlZzNA3MDTJEUBzcmj5kxbaYYSpcpDtjl5vWdguViucahdpFmsxKtqLl01NIms/HnUX+xGnjzMw4WrTJNY88clG9665qEN91xqTRvo7tpaNKabRE2+udbcbx6fRKq2LAMVJ8ussklKsOZljeuULwhN7RouGZHeM7BcLNd7CVhx8kmU1PVxyfKFHsNGq6aaVG2O3fJRg2PX96DWqBH4Zm22jo+db+/iimIpAYefNpr+xn9mXJj8+WIVnLyJHcSRrhVs6fswqziEqq3g9Bl1LMX4zqafXB/VnHzTR+bBbeoPM1AtEl1cnZ8cTl75opO18Zc/9o41rea/X5z8qwH8uf59M4CbrM/u1ZjUuNcAeA0AbNu2bdUXW7J3UzozdJiTl23NjzlCrc9uL87BoniYGbzqoPoPirJELPSh8eB1G1eNqkDbPF9Yv0inXckpjQWAHZ9w/0jGK00bm7dkw7E2+6k4eRLwAmwVWGah4Z6GDGB9kXzztGsC52E96dkr3sbfWasz+7rS8YZCkFaqqp/WsVpjidnNAS+qIuIq15wdT/j8bo072V36H6vl2EM+whpxgh66ju+qNsIfdipP3uHkUVQ7Wxs+gnqjZm5y7HqaVfcn1OeQ74L7LQI+Mg8D4fylM16k60YU4Ymlw2Ib12NHpGuI6G+J6EvCfy+yjnkLgAGAD621Acz8HmZ+OjM/fXZW1qKRbLw7huM37RA/mxuflvGA/vxUd0zkzPK0JXJmEWQ9bi5lTRsCAhyuXMFK5VT7xydRLHL+EZGoy0MgX5pY21xf9sVMd0LkWPtZR+RY0yhZE9+slCsln8oZUSFdloRisTZnRJGoA09EQS61Sj91bCZQgm0878v1ReNMHi8k6/sw5Px59UIm33/JRwnFYq56TJG4V4FA6Af45lAcIlSObiIPxHCyPBAzCe0lWBsnHwcy6NRbk+SjRB4vkDWNIoqC4yU0LkL4WN4P1i+WYnXLxQCnbd4rnms9dsRJnpmfx8z7hf8+CgBE9BMALgFwKdeP5vsA2FHRLRo7qjbVGzeNbPw0XGp1r/XPvpkIneN7mdb3sHFWAbbS5eRZ6bIoTZsmToDm5Jt4WegqNfaixqw06qVf9TdFUSCl2MPLsqwHpIUzA31TEMXGwfXmFtdHmYyr49nDe5pXdH3Ujloqo8jpc4JI9AWXOrvG4SGLsqheq5t8c6Gb4/PNqcCxMpdop63GarrmiXMPB2qRLZdjzfVeBfc8rbglcqytuCVyrBGTiJdl2cgZt31RaH7a/szWa2n4iItKZ973kR+3YLDWn2fZF66PWPYRM1TuOXxf5HHeqPFq8ASx6Asu0eDkK19wWeH2f4OyqDj55ngZiOOo4EKMW5RQ980fR6oP8nhpiz7KY3kcdZKW6Is8zsQ8+YTiYOH09diT4uSJ6CIAbwdwHjM/ZOH7APwJgDMBbALwCQC7mdnf82/ZWjj5QwuHMfHy/Wrzkd8w63fUE2vokVaRrBbGms9m9vAkSTCANflry9IMA1YBRJsbVzy05uQtiyNVMUpx8hURCRBVueHuteMkFikhT5en6rPsiwZXb52/4SMnTlAdX71nQnPypYenSYoBCt2mmqtvpRkGKNWDwepzEseIoxiLg+Wq3iUzV1W1lguFE7QSA6taroOy8Gq5Gn34eocja7c6fDOzru2qcaNFU+FRdZ9qTZuyWpUZzRSwStdMY78OblmqzBrlo2YOeBarTJCmJr7CGSzkyasd2O5eArCJTwyqGq9c1nGIkktVB5cALqHr36Li2JUv0OTkrb0ENsXg+qhxvL4HNs9NVKnVaB2fqNoX0NgzECWKlqH6rVfVGUhRshov9l4CU6/AxHAIpp1qd/mgWG7kwxsNIcPtm3HU0OuJ6rhFU4rB8QWj4tht3Myl7v0Hq3cLFyfNybt4L+vgLy97B07fsg+rte8lJ/9uAC0A1+sbdhMzv5aZbyOivwBwOxSN87NHmuDXagN7Jb2SreUZ5hxrbr7Ii0u6OWwFbKyPQ3nyIa6++btw7ZX6tMrPWMBWPpfNyQeOdfBqAWH1R636Sh+3JhO7Hm3JDMSsuVS5sS7/WXPjgeNXy8mqOUPMw1fnkc9vVrn2acw4Yq/Pwp6BkqvJw8ZLve3MvzZX/w+1tXDaav61ov6T8FGY8w/7QvriS+ep6ss2fMT6e9LECy4Rw4/haKet6ItQn4PHO77gAG63WfoKu3joeCLCwvKi2Jb12JNKoWTm45l5KzOfqv97rfXZbzLzLmY+kZk/9uSb2rSxTh97tsgqbRsmZW5/Q0B/fro/KXJmnawt84GQNaiDtTxLazVrn4fiYE6vqLMRyTm9EUWiXg+BMD8mxyc2jsm+mOtNixzreN6T9VRieS8BAr4wGRZ+W+VJOYli8YuRBLjdmCJxSzgRBfcGzPcDMZxA3GKyHYjhBHLDVc1WmZMXc8YDE2bkbNc3FtJliSgS65RGRJjIpfESjk/MdWUfTbfHgzEcMdc/MF5CvlBvHULfAjEcFbfwTenGy+OlL9R+JlAwDjHfk8fFbHel8SLML2lL/C4sF4M1reKPZMesdg0AdHOfhwaAtpkIHbyaIKtHK/TxWovbOz5rPq31ajWJYj9/Xv9XcfjWf8ysqtE4eFkW6o3Aw8taZ966BjMjU+WTPDxPMrN8tfrAaMcBX5gvv4vHujBFo2/qeM9HzEgpFvLnVXyiFPqMglGaPQNOn1nwheH7PS5VS+96OHOAq2fNT6OBA0CL0sa/zc+UEg9nVkE8iZOPEYu54aZsnMc3l0WAby4bfHOFF0WDk7d5a8lHhgOXfJSQjzOAFGnjjcPwyqaCmYtLPLfilRORh46YRE0b20f2f2VZVhx+00eF6KOiLMACV1/qHdVSnxMIfYDqw1p8kYnjhZFQKvsCcWNvSI2THJBfpx2z2jVPLBzC5Ev31Ty3TWMEOfm1cfWVdo1jSZxoiVhu8NZZmmLZBFJsPMnABBU/sPAoUm8Ebp48iJCmicfhA5oDtzl//TdVnnzIFy73Hgd8QRCPb5zfbk9kxQisPqRJioHgi1aSooD64lrJ4So3vNLWr7l6IkKapB4O6NiIw1sDRu9f56tryo3ZcoXFHzNzlXWhPlefqOOpPt6iTwwfrs5jzs+15rqTJ53GiY7TUINvzpJU77do6vukcarz5P16tyInD5WZslwOPF9EVp1dw6WrtjbrA5i+hfpcDZUVceW7WhvH56ENL272Ehg8jfR3KhC3cPPhE4qrWgx+fCLCklUr2Dyc4yhW3H6FNzn51fbZfOb61PzJysdbvq44+aYvei3FyT9j+8lYrQ2ldo3Hq4aeVWvFV2Vcn8A9D/u4nX3QwNnm5IVzHKkJa/nsSR+/jjaJvrD4W5t7B9fCBlb/G6cM+CW0UFELLP8zs4KTjpf6aX+B/fPYk4GeONjJk9b5/2D2+GbTlmYf6rOWbnxC5OSbbQrhLPh1xb6tGV+N78yq1W+rGLew+tTQmUft+4YvWPadbauNpawOX+N4EXHnoQKq6yMfBTtm6Zp+p4d922XN5c3Tss785ql5EZ8bmxF5v37eFXNrkygRdeZVVSKJk5cDOjGF3hRiOceYInQyPwc8okjU5SEQNk3MidfePCH7YuPYrMixTrXHRb45T+TccEDWVDFKh15bhS89ACSBfOg4isQ8/IiiOoXWPj9RkHvfNCb7aGNfjlvMdCZFffBuYC9BHJHMQzPE8eUGSY0RycVqjDqla8Ec8ADfTKhLIboW8sVcd0ocL2OtkB5UJvLQIDkf3qh7SiaNlziKxPVITOHxItUNIBDmAtx7KJ61sT8j7leY6UysabwUZYEzth4Qr7EeO2ZX8oBcrAIAEuGLAwCR2Xhg3t31T6NOp94vUQ2SqMKbo4YAzR/Dw0sBZ2a1zc9duXAJKv1BySWDpOOh08vq90LdB10OzaxEK1wFiaU+V2WwbRyE2OQEujjXkq2w2kAlai7dtqJEaXLznL5V110lHgk4SjR9VLWJEZXUeEVWr8YAscobNK/mFT2h772N23yVOc5+GzP64DbOBTe4+qqpBYvS1INyoKgCBy+4AEu44abhHl8CKPzr6mu7viiJdEWtJg7tC993BJS+7wCq3pZ8H9U0ho2XBVeceaMPgjKp6RuRj6tNhL58SanHru8LBgu+cHc0176o5RBcX9RlFZq+KPVKvekj/XYn+KgsS79uMtRDXloIrdeGh5O3LeSgkN8CHH4UUaDGayq+TmWp4qGl2q+Anycf6Tx5TxOf4OfJ61MG8+QDNV7DcYgAV2//27p2SAMnjmx9e0M4AmmaYrmKW9Tpl1mi9X1KS3+n1PVu46jWIiKFExHSNPXq4AKoar+6146SuKoD0JhEK06+doaa21bDvdfUQVQdj8YkYTh5F091DIcAUBQBrB4GWZzqFEpWODRXH6cAqY11kYVHes/AUrFcvQGYxADDyat8eKomKJMnX+eG6+CxU7PXowwcvvnIXP3qfGHy5D0faU6e0ayDm8aJzpMvG5x8TDHiiLCk6wwQmdxz9Qbpah2ptybFydu4rZXv4tVXxWgaMVf+U18hOU++rv1rzlM/KGw5amLzEKtjNUVZot/q4iOv/n2cvfM0rNaGkpP/nj2c7NOS82/xoNV+IhxbzZqrbM+KbVpFA6TP1nN8EG9+UDW1gVu/u88qtn6upa3CtQ0uc/LyyVbm3iXcmeispohcLSBSdy7O0PQVmrnb9hlt3KZLXA19Y+vjodfiI7fPHMBDx6/QVg71oRlYrc5bXfvoxBtWricgHK86J56/WCV+NGe3Y5aT77W7OOU4Wd9h2+wmEd86s0mcUjdOzoqvR+PtMTk3PJF1oLmEeLwroWssolqUybYYcm3OKIrQE2pwRkSiLg+BsHVKjk9sm9oo4psn50WOdaYb3ksg572TyDebbAbXmEvRF6HXVlWbU8oBlzVHCIRNAS5126Tsoy3j8+J4mQ/sJRhry3sJkigRfaTeFiQfFeLDyazEXYtI9nVEJOeAE2FGyIcnAJvH5fjElvFADKc/K/posj0e1IMSa78SBXxh3gRdY3m8RJF4byKQuGcgjiJMdvzxEhEFufet40dnvPQDe0+YGWdu2y9eYz12zK7kAaAcFGapWINEKg8baD5NiVAWpXc4SOVzs7AarM/TvG6leuc+fbkAOPKf4iWDzd7uxh+QWsG5K2BCrXVj4wwUg0JYJVCde96ExUpVIFL69l4fCMWgVA0VOFMIKxR1XaFvBeulreOjUsZR6hWYhBe+jwBz/1eJEwWOJxRLhbfqJxCWB4PGytDY8mBQ8dC2LS0PxBVq6PhBWahdm65P2VfpBDSvW/rcO0PVIPBxwvJg2cMJwOKyjwMUwIGFJRlfHCyJbV1cWhLjE4vLyyJuaDrJRypU4IxHrWQK8n1XCr5mIiwuLQm+hogDhGXBFwTCkuBT0wdpvCwtLYv339wbr8963B0tO2Y5+ccPP4GpH94n8taNHHDb1hjMiOJA9ovWTHEtS1L1ZbO5cVb8NAgeJ0+kNj1IdWrTNA3kyUfNPHndvEq7xutEoM8hPAQbTt5rjxwjSJNEiW/Z52Tlo0IHnGy84uQLPz6RxSFOPqk5eQuPEqs2L6zPzL+PFLeojq/jFhVsnYcM5w9U/HfNv9bHK05+oPGaYsjitJLRdXFAFVKx8YiU+umSrbmvx2cWp2KevNGTrzJU9CRkqKWVtGiMG+w5InS8iVsAgnYN1RrqNSff5KGrPHnUbzf2HgMVV6jxWKuQNuITeuI38azI8oXRuxm4NYF1nVrjX5PB5MUnLF+Yj3ycRRya7gvlybvjYizv4iM/9U48e9fpWK0NJSffeP35vj+nwg+LKnpuL5yFLADAiDoJ51rNw+jfzrM5bFYwtIJAMFkczguYfIrV3OejhYc+4wAsrNrsz2xrfOlt3Apu+n/RPN72hTsJSXizPTK+Vi2alfhm2Uer56EBeHEIQl10pIFbw8I7PtAH869VaxdZbQ31bTU423hg7Nn3h1lOJ12vHbOcfLfdwdN2q1zSDVNz+NTvfRjzWrNmx5xc+3X73GZxet48PS/m4k52J1bI9RW4NEA8T1H6GTeAupkhHlrKw4+Iarlky4gIGyb8/GYCsH1ajk/snN0i4lunNsq8Yn9a5Fh7rbas1wI5178Q3oAAVKs118yX3Ds/EdoBHfBpQaOGQNgS2BuwY0r20Y5JOYazKbCXYKIt64ansVyPFszi+YuyELO6pIcEoO6zNF6IKJgDPt+X8+G3Tcixmu2TMr5pbE7sw0x3Urxvnawt5/RbK3vbVDzLP7/acCaNF4j3RtUZ8MdLHEWiDjyBgnGI7YHxEvKdyp9f/XgBgDO3j/LkAagSgADw1kuvwDn7z8RbL70CP/fuX8bColZwswcHKe5N4S5PuCSef3FpQd0al0tbXlI3zR19RalkcgVOXrpuRBBzzCMQisJX2SSKsLS85HP1RFhcXBR56MUlc3xz2bywcNjHASwtLkKRo80v0MLiAsA+vrS4LPqiLAaICH7WQMmgyPcFlfp/ni8gxhsiirUvXM6UA74AFhcWhT4TFhcXRHxhQcKBxcVFHVYQfMcCJ78U4pvLKqXOtgIMEoL1FfcOFycsLfu8MoFwaHFBfHsQcRAOCzgAHNI+cq99eGFBfKN5YvGw6IsF/Z2V4hYQfDHgopl/b/liUJbedUERlnQMxO3zwuKiGLeQfAECngj46NDCoTX56JCej9y/Obwkj5fFwbKK1xyl1fwxy8k/duhxTL9kP6bHp/C1D/4D2q0chxYP47gfOxsPfOdh+Y+OpCfvHh6s8RpjYAKXDve+LOTPp3ECJlJFNKzjlS7LGjn5JBYlUU1ert+JVXDvdlqm7QsLb3D+Vh8anLztiyRp6gppPEtSDIRMmjiqOVaXF89MnrxjtYaQ0+WImpy82zePk4eMh44PxHyMTkzjWsyNdkpcPcPSUOdavdPWyjf8d2Zx8oZXBoAsrrVrbDzRefIuV29urYKbm8dg4UCTYlgJl/Lk7XOquEVUjV8bT6PaFw08TlDot7zGeSjSNYGdOIT2q50nb1IvDSdvH2/2ANh7AfyHYrPPpeUrmauvKVph3VXd/yYnX+NjeRd/9dPvxHN2n4HV2kqc/DFL15in3FsvvcIKFkV466VXHMWrhB6AFvPH9WEUSmWLovpV1Do+0sUovLOHJuajaI0rhLq5Ei6SkU1c6rMdcLMtIssXto+C8Ymwj9j+ZaV2up+tFQ9fGbD4XZtCsGlf5Yt6U1M1eUDGlY/M6bkx2dSXleuL+ry/3R7peJlDXgm3d4javmhcm0nGAW/ydX1h4xHV3yk36NnwRUPvJuSLI+HOJG6fX/gb8yBwj6/P5V2igTMrBdyjZcfsJN/J23jZeZfgJy98BfJMyenmWQs/edErcP4pZyNLmzoVWZriwPY9aGdNXi6JY+zZsgudVlMTJqIIW2c21XLGlk10x9Bu+fxeFqWBGq8kc7WB/Ok8aaHX9vPhO1mODRNz3oOh22pj98bt3rnaWY79W0/w6tG20gwHtu1ByxlIaZxg35bj0XF8FEcRds1tRdfxkdGrd3EA6KQ5WgHt+1C8Qa6bmnrtBJSctKTX0221sWN6s+fvTpZj78Zdnp5OnrZw8uYTvfuWJSkObNpdy1Zri6MYJ87u8HxERNgysQFdQVtoLO+hnfm+SKNE74ZuWkTyeEmiWB4vcYZ+LoyXNBfjB500x+7Zbf54SVvYv+F4f7wkKU7eeIKXZ57GCU6aP86LjUQUYefUZtEXs70JUX8pz1rBfQ+STEkU0OtpJamYD99OW2I9gV6rgx1TwnhJc5w0f5xXdzZPMpy86QS03PESpziw8XhxvOye3Y6uOF7mRR91szbO3Hb0OPljdpIHgPf8/G97GyiSOMHf/Lc/xv7tJ6KXd9Fu5ei1u9i7bTcO/taHcNZJT0O31Ua7laPf7mH73BZc/7YP4eIzzke7laOj8fnJGVz/tj/Bjz/3pcizFjqtNvrtHqb6Ezj4m/8Lv/jSy5GnLXRbHfTbXfTbXVz963+M3331Lys876jrZzn+9Bffhfe/4ffQznL08g56eQd52sK7Lv91fPTN70U3V+fottrI0xbe+orX47pf+QAmOmPot7voaPw1F7wK1//qBzE7NoV+u4dO1kY7y/HiMy/Awbd+EFumN+jjc3RabZy390wcfMv7ceKmnVVbuq0OTt1+Eg7+8vtx+s796LY66GRt9PMujt+wHQff/AGcv++Z6GQ5OlmOft7F5skNOPim9+NHnnER2hY+05/EdW/8Y7zuuT+KPG1V5xlv9/Dx//pH+I2XXKF9pPBO1sZHrvgDXPUTv448baHX6qDXUr5430+9DX/xs79fncP44r+/7Bdw7Rveg37erc6Rpy385wsuw8Gffx+muxMaz9FOW3jVmS/E9T//Pmwcn63wTpbjgr3PwsE3/E8cN7sF/VYX7TRXX6YdB3DwP70XBzbtRq/VQTvN0Wt1cNKG43Dd69+Ds3edhm7WRjvN0W91sX1qE65//Xtxyb7z0E5zdDQ+35/G9T/3Xlz2jBcpX2h8sjOGgz/3HvzS816NPGmhm7XRb3XRa3Vw9Wv/AG9/8X9BnmToZu3q+n9y2W/jgz/2NrS1j7pZG3nSwrte+kb81U+/U5/D4BnefMHP4NrLr8REu49+q4tOmiNPWvips16Ca197Feb6UxXeTlv4d/ufg2svvwpbJ+bRb3XQ0f0457in4ZrLr8SJs9t1W1R7D2w8AVe/5g9x2pa96GZtdLSPjpvegmsvvxI/tPvM6hz9Vhebxmdx9eVX4qWnPL/hi+nuBK55zZV43dkvR55kGu9gLO/hb3763fjVC/8jWpYvOmmOP7/sd/HuH3lL5bte1kGetHDVy38Ff/rjv1O1pZu10Uoy/MYLXo+P/sy70W91VZ8z5YvXP/tSXHv5lWq8WL542akX4GOXX4WNY7MV3klzPPeEs3Dt5Vdi5/SWqi3drI3Tt+3HNZdfif0bd6OX1ePlxLkd+Pjr/geetfM0NeZ037ZNbsB1r7sKl+x7DtppqxpHc70pXPe6q/CTz3hxw0cT7T7+5vJ3ixvD1mvHLCd/JGNm3HDrP+KOb9yNPdt24dwDZ1Xc2+e/8kV88au3Y8f8FjzvtGdXDr3lnjtw0x3/jE3T87jw9POqVdZX7v0aPv2lmzDdn8TFZ5yPXD+Vv/Hgfbj+Czegl3dwyTOeV636v/XIg/jYP30KaZzgkjOeiwldcPzbj30H19z8dyjKAheffj7mdEbMY4cex9U3fwKHFhdw4WnnYsuMitIfXlzANf/0d3jk8e/ihw6cjV0btwMAlpaX8PEvfgbfeuQhPGvP6di37QQAasPS9bfciK8/dB+evusATt91oPLFp267CV+5/x7s33oCzj7x9MoXN931Bdzy9Ttx/IbtOH/fMytffOGe2/D5r96KrdMbccEp51Srvjvuuxs33Hkz5sdn8IJTz0WmV033PPhNfOK2z2K808cLT31O9WZ037e/hetuvRF52sIlp52PMZ0d9PBj38Y1//JpEAgvPPU8TPcmAQDfPfQYrv7iJ7E4WMILTj4XG7WK5hOLh3DNLZ/Go4cfx/P3no3tM5sBqKD5tbd+Gg8//gjOPeEMnLhhJwCl+3LdbTfivu88iGfsPBmnbN0DQGVrfOLOm3DPw/filC0n4sydJ1e+uPHuf8Lt938VezYch3N3P716Y/r8/70VX/jmHdgxvRnP3XNW5Ytb7/sKPnvPv2Dj+Cwu2vusagV814Nfx6fu+jymu+O4eN+5VbGabz7yLRy88x/Qzdq4ZP956OndqA88+jCuvf0GNV72nYcJvQPzkUOP4urbPo2iLHHx3nMwp1U0H1t4Atfc/hk8sXQYF+55VpU1dHhpAdfefgMeOfwozt99JnbNqCyz5WIZH7/j7/HAYw/jmTtPxb4Nx6vxUhb42y/fhG88cj9O37oXT9u6txovn/7qzbjroa9j34bj8cwdp1Q++sev34Jb778Lu6a34jnHn1GNly/eeydu/uZt2Dwxh+ef8Mzqbe3OB76GG+/5AuZ6U7hozznVG9M9/3ofPnn35zCe9/CCk86pVvf/77sP4eCX/x550sLFe59dKYo+/MQj+NjtN4KI8IKTzqkyqL57+DFce8cNWBws4cI9z6p2qR5aOoxrb78Bjy48jh864awqg2pxsISP33EjHnr8EZxz3NOwZ74eL9d/+bO477sP4Onb9uPUzfV4+eTdn8M9/3ovDmw8AWduP1D54u/v+QJu/9ZXccLcDpy3qx4vN3/jNvzzvXdg++RGPO/Eerx86f678A/3/As2js3gopPq8XL3Q9/AJ+/6HKY643ihNV7WYitx8kM7yY9sZCMb2Q+KDWXgdWQjG9nIRnZkG03yIxvZyEY2xDaa5Ec2spGNbIhtNMmPbGQjG9kQ22iSH9nIRjayIbZ/U9k1RPQQgK8/1e1Yh80ACGgpDK2N+jz89oPWX+DY7fN2ZharnPybmuSPVSOim0PpS8Nqoz4Pv/2g9RcYzj6P6JqRjWxkIxtiG03yIxvZyEY2xDaa5I+OveepbsBTYKM+D7/9oPUXGMI+jzj5kY1sZCMbYhut5Ec2spGNbIhtNMmPbGQjG9kQ22iSX6cR0e8Q0Z1EdAsRfYSIJqzP3kREdxPRl4nowqewmUfViOhlRHQbEZVE9HTns6HsMwAQ0UW6X3cT0Ruf6vZ8L4yI/oiIHiSiL1nYFBFdT0R36Z+TT2Ubj7YR0VYi+iQR3a7H9RUaH6p+jyb59dv1APYz88kAvgLgTQBARHsBvBLAPgAXAfhDIjo6FXmfevsSgJcA+IwNDnOfdT/+AMALAOwF8Crd32Gz90PdO9veCOATzLwbwCf0v4fJBgB+gZn3AjgLwM/qeztU/R5N8us0Zj7IzKbS9k0AtujfXwTgz5h5kZnvAXA3gDOfijYebWPmO5j5y8JHQ9tnqH7czcxfY+YlAH8G1d+hMmb+DIBvO/CLAHxA//4BAC/+frbpe23MfD8z/7P+/TEAdwDYjCHr92iSPzr2agAf079vBvBN67N7NTbMNsx9Hua+Hcnmmfl+/fu3AMw/lY35XhoR7QBwGoB/xJD126+oPLLKiOhvAWwQPnoLM39UH/MWqNe+D30/2/a9stX0eWQ/eMbMTERDmW9NRD0A/xvAG5j5UbIKyg9Dv0eT/ArGzM9b6XMi+gkAlwB4LtcbDu4DsNU6bIvGjgk7Up8Ddkz3+Qg2zH07kj1ARBuZ+X4i2gjgwae6QUfbiCiFmuA/xMz/R8ND1e8RXbNOI6KLAPwSgH/PzIesj/4awCuJqEVEOwHsBvC5p6KN30cb5j5/HsBuItpJRBlUgPmvn+I2fb/srwFcpn+/DMBQvcmRWrK/D8AdzPx266Oh6vdox+s6jYjuBtAC8K8auomZX6s/ewsUTz+AegX8mHyWY8uI6IcBvAvALIDvAPgiM1+oPxvKPgMAEV0M4B0AYgB/xMy/+dS26OgbEf0pgOdASe0+AOBXAfwVgL8AsA1KAvzlzOwGZ49ZI6JzANwA4FYApYbfDMXLD02/R5P8yEY2spENsY3ompGNbGQjG2IbTfIjG9nIRjbENprkRzaykY1siG00yY9sZCMb2RDbaJIf2chGNrIhttEkP7KRjWxkQ2yjSX5kIxvZyIbY/j928blxKwa0oQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Plot the surface, extrema, and gradient descent data in a matplotlib scatter plot\n", + "\n", + "fig = plt.figure()\n", + "ax = plt.axes()\n", + "\n", + "ax.scatter([x for (x,y) in param_space.keys()],\n", + " [y for (x,y) in param_space.keys()],\n", + " c=list(param_space.values()),\n", + " cmap='Greens',zorder=1);\n", + "ax.plot(optimizer.GetMinimumMetricValuePosition()[0],\n", + " optimizer.GetMinimumMetricValuePosition()[1],\n", + " 'kv')\n", + "ax.plot(optimizer.GetMaximumMetricValuePosition()[0],\n", + " optimizer.GetMaximumMetricValuePosition()[1],\n", + " 'w^')\n", + "\n", + "for i in range(0,len(iteration_data)):\n", + " ax.plot(descent_x_vals[i:i+2],descent_y_vals[i:i+2],'rx-')\n", + "ax.plot(descent_x_vals[0],descent_y_vals[0],'ro')\n", + "ax.plot(descent_x_vals[len(iteration_data) - 1],\n", + " descent_y_vals[len(iteration_data) - 1],'bo')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can also view and export the surface as an image using `itkwidgets`." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "x_vals = list(set(x for (x,y) in param_space.keys()))\n", + "y_vals = list(set(y for (x,y) in param_space.keys()))\n", + "\n", + "x_vals.sort()\n", + "y_vals.sort(reverse=True)\n", + "array = np.array([[param_space[(x,y)] \n", + " for x in x_vals]\n", + " for y in y_vals])\n", + "\n", + "image_view = itk.GetImageViewFromArray(array)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "becd8d61582a4e24ad9dd5ef4a6a36b2", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "Viewer(geometries=[], gradient_opacity=0.22, point_sets=[], rendered_image=" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Visualize metric over gradient descent iterations as matplotlib subplots.\n", + "\n", + "f, axn = plt.subplots(len(point_set_sigmas),len(estimate_rates),sharex=True)\n", + "\n", + "for (i,j) in [(i,j) \n", + " for i in range(0,len(point_set_sigmas))\n", + " for j in range(0,len(estimate_rates))]:\n", + " axn[i,j].scatter(x=list(hyper_data[point_set_sigmas[i],estimate_rates[j]].keys())[1:],\n", + " y=list(hyper_data[point_set_sigmas[i],estimate_rates[j]].values())[1:])\n", + " axn[i,j].set_title(f'sigma={point_set_sigmas[i]},est={estimate_rates[j]}')\n", + " axn[i,j].set_ylim(-0.08,0)\n", + "\n", + "plt.subplots_adjust(top=5,bottom=1,right=5)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAW6ElEQVR4nO3df5BdZ33f8fenCNNiKKzAyKpEYpI6Qx0SINnaaUtTt5ZlEdqRp0kImUy6BDweZkgJ7XQStW6qjh2myi9C0+ZHVaOgMIQQCIlFIVHWKqZMEojXVIANhjXYWFL1Y+11jI35YeDbP+5ZfHS5u3vvnpUlfN6vmTv3nOc893uf556793PPuXd3U1VIkvrrb5ztAUiSzi6DQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSem5dgiDJjiSfSnJXkl0jtj8lyTua7R9OclFr279v2j+V5Kr1GI8kaXydgyDJk4DfAF4KXAL8eJJLhrq9Gnigqv4u8GvALza3vQR4BfDdwA7gN5t6kqTHyXocEVwK3FVVn62qrwC/D+wc6rMT2N8svwu4Ikma9t+vqi9X1d3AXU09SdLjZMM61NgCHGmtHwUuW65PVX01yYPAs5r2Dw3ddsuoO0lyLXAtwPnnn//9z3/+89dh6Gvz8WMPrul237PlGetaYz08keYCaxvL8DjWo8Z6OFf2jTXOzHPkbDzPbrvttvuq6oLh9vUIgsdFVe0F9gJMT0/X3NzcWRvLRbveu6bbze152brWWA9PpLnA2sYyPI71qLEezpV9Y40z8xw5G8+zJJ8b1b4ep4aOAc9trW9t2kb2SbIBeAZw/5i3lSSdQesRBLcCFyd5XpLzGHz4e2CozwFgpln+EeB/1+Cv3R0AXtF8q+h5wMXAX63DmCRJY+p8aqg55//TwEHgScC+qrojyfXAXFUdAN4MvDXJXcAig7Cg6fcHwCeArwKvraqvdR2TJGl86/IZQVW9D3jfUNt/ai1/CfjRZW77BuAN6zEOSdLk/M1iSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnOgVBko1JZpPMN9dTy/SbafrMJ5lptb8hyZEkD3cZhyRp7boeEewCDlXVxcChZv00STYCu4HLgEuB3a3AeE/TJkk6S7oGwU5gf7O8H7h6RJ+rgNmqWqyqB4BZYAdAVX2oqo53HIMkqYOuQbCp9UJ+Atg0os8W4Ehr/WjTNpEk1yaZSzK3sLAw+UglSSNtWK1DkpuBC0dsuq69UlWVpNZrYMOqai+wF2B6evqM3Y8k9c2qQVBV25bbluRkks1VdTzJZuDUiG7HgMtb61uBWyYc5znlnj0vOydqSNJ66Hpq6ACw9C2gGeCmEX0OAtuTTDUfEm9v2iRJ54CuQbAHuDLJPLCtWSfJdJIbAapqEbgBuLW5XN+0keSXkhwFnprkaJL/3HE8kqQJrXpqaCVVdT9wxYj2OeCa1vo+YN+Ifj8L/GyXMUiSuvE3iyWp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknquUxAk2ZhkNsl8cz21TL+Zps98kpmm7alJ3pvkziR3JNnTZSySpLXpekSwCzhUVRcDh5r10yTZCOwGLgMuBXa3AuNXqur5wIuBf5TkpR3HI0maUNcg2Ansb5b3A1eP6HMVMFtVi1X1ADAL7KiqR6rq/QBV9RXgI8DWjuORJE2oaxBsqqrjzfIJYNOIPluAI631o03bNyR5JvAvGBxVjJTk2iRzSeYWFhY6DVqS9JgNq3VIcjNw4YhN17VXqqqS1KQDSLIBeDvw61X12eX6VdVeYC/A9PT0xPcjSRpt1SCoqm3LbUtyMsnmqjqeZDNwakS3Y8DlrfWtwC2t9b3AfFW9aZwBS5LWV9dTQweAmWZ5BrhpRJ+DwPYkU82HxNubNpL8AvAM4PUdxyFJWqOuQbAHuDLJPLCtWSfJdJIbAapqEbgBuLW5XF9Vi0m2Mji9dAnwkSSHk1zTcTySpAmtempoJVV1P3DFiPY54JrW+j5g31Cfo0C63L8kqTt/s1iSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6rlOQZBkY5LZJPPN9dQy/WaaPvNJZlrtf5rko0nuSPLbSZ7UZTySpMl1PSLYBRyqqouBQ836aZJsBHYDlwGXArtbgfHyqnoh8ALgAuBHO45HkjShrkGwE9jfLO8Hrh7R5ypgtqoWq+oBYBbYAVBVn2/6bADOA6rjeCRJE+oaBJuq6nizfALYNKLPFuBIa/1o0wZAkoPAKeAh4F3L3VGSa5PMJZlbWFjoOGxJ0pJVgyDJzUluH3HZ2e5XVcUa3tFX1VXAZuApwD9bod/eqpququkLLrhg0ruRJC1jw2odqmrbctuSnEyyuaqOJ9nM4J39sGPA5a31rcAtQ/fxpSQ3MTjVNDvGuCVJ66TrqaEDwNK3gGaAm0b0OQhsTzLVfEi8HTiY5GlNeJBkA/Ay4M6O45EkTahrEOwBrkwyD2xr1kkyneRGgKpaBG4Abm0u1zdt5wMHknwMOMzgaOK3O45HkjShVU8NraSq7geuGNE+B1zTWt8H7BvqcxL4+13uX5LUnb9ZLEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HOdgiDJxiSzSeab66ll+s00feaTzIzYfiDJ7V3GIklam65HBLuAQ1V1MXCoWT9Nko3AbuAy4FJgdzswkvxL4OGO45AkrVHXINgJ7G+W9wNXj+hzFTBbVYtV9QAwC+wASPI04N8Cv9BxHJKkNeoaBJuq6nizfALYNKLPFuBIa/1o0wZwA/CrwCOr3VGSa5PMJZlbWFjoMGRJUtuG1TokuRm4cMSm69orVVVJatw7TvIi4Dur6t8kuWi1/lW1F9gLMD09Pfb9SJJWtmoQVNW25bYlOZlkc1UdT7IZODWi2zHg8tb6VuAW4B8A00nuacbxnCS3VNXlSJIeN11PDR0Alr4FNAPcNKLPQWB7kqnmQ+LtwMGq+q2q+jtVdRHwEuDThoAkPf66BsEe4Mok88C2Zp0k00luBKiqRQafBdzaXK5v2iRJ54BVTw2tpKruB64Y0T4HXNNa3wfsW6HOPcALuoxFkrQ2/maxJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUcwaBJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRznYIgycYks0nmm+upZfrNNH3mk8y02m9J8qkkh5vLc7qMR5I0uQ0db78LOFRVe5LsatZ/rt0hyUZgNzANFHBbkgNV9UDT5Seqaq7jOCSdIffsednZHsK6Wa+5rEedc+lx7XpqaCewv1neD1w9os9VwGxVLTYv/rPAjo73K0laJ12DYFNVHW+WTwCbRvTZAhxprR9t2pb8TnNa6OeTZLk7SnJtkrkkcwsLCx2HLUlasuqpoSQ3AxeO2HRde6WqKklNeP8/UVXHkjwd+EPgJ4HfHdWxqvYCewGmp6cnvR9J0jJWDYKq2rbctiQnk2yuquNJNgOnRnQ7BlzeWt8K3NLUPtZcP5Tk94BLWSYIJElnRtdTQweApW8BzQA3jehzENieZKr5VtF24GCSDUmeDZDkycA/B27vOB5J0oS6BsEe4Mok88C2Zp0k00luBKiqReAG4Nbmcn3T9hQGgfAx4DCDI4f/2XE8kqQJdfr6aFXdD1wxon0OuKa1vg/YN9TnC8D3d7l/SVJ3/maxJPWcQSBJPWcQSFLPGQSS1HMGgST1nEEgST1nEEhSzxkEktRzBoEk9ZxBIEk9ZxBIUs8ZBJLUc13/Z7EkreqJ9j9+n2g8IpCknjMIJKnnDAJJ6jmDQJJ6ziCQpJ4zCCSp5zoFQZKNSWaTzDfXU8v0m2n6zCeZabWfl2Rvkk8nuTPJD3cZjyRpcl2PCHYBh6rqYuBQs36aJBuB3cBlwKXA7lZgXAecqqrvAi4BPtBxPJKkCXUNgp3A/mZ5P3D1iD5XAbNVtVhVDwCzwI5m26uA/wJQVV+vqvs6jkeSNKGuQbCpqo43yyeATSP6bAGOtNaPAluSPLNZvyHJR5K8M8mo2wOQ5Nokc0nmFhYWOg5bkrRk1SBIcnOS20dcdrb7VVUBNcF9bwC2An9RVd8H/CXwK8t1rqq9VTVdVdMXXHDBBHcjSVrJqn9rqKq2Lbctyckkm6vqeJLNwKkR3Y4Bl7fWtwK3APcDjwDvbtrfCbx6vGFLT3z+fR49XrqeGjoALH0LaAa4aUSfg8D2JFPNh8TbgYPNEcR7eCwkrgA+0XE8kqQJdf3ro3uAP0jyauBzwMsBkkwDr6mqa6pqMckNwK3Nba6vqsVm+eeAtyZ5E7AA/FTH8egs8F2n9K2tUxBU1f0M3skPt88B17TW9wH7RvT7HPCDXcYgSerG3yyWpJ4zCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknqu628WS+cU/z6PNDmPCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ6rlMQJNmYZDbJfHM9tUy/mabPfJKZpu3pSQ63LvcleVOX8UiSJtf1iGAXcKiqLgYONeunSbIR2A1cBlwK7E4yVVUPVdWLli7A54B3dxyPJGlCXYNgJ7C/Wd4PXD2iz1XAbFUtVtUDwCywo90hyXcBzwE+2HE8kqQJdQ2CTVV1vFk+AWwa0WcLcKS1frRpa3sF8I6qquXuKMm1SeaSzC0sLHQZsySpZdX/R5DkZuDCEZuua69UVSVZ9oV8Fa8AfnKlDlW1F9gLMD09vdb7kSQNWTUIqmrbctuSnEyyuaqOJ9kMnBrR7RhweWt9K3BLq8YLgQ1Vddu4g5YkrZ+up4YOADPN8gxw04g+B4HtSaaabxVtb9qW/Djw9o7jkCStUdcg2ANcmWQe2Nask2Q6yY0AVbUI3ADc2lyub9qWvByDQJLOmk7/s7iq7geuGNE+B1zTWt8H7Fumxnd0GYMkqRt/s1iSeq7TEYG+9d2z52VnewiSzjKPCCSp5wwCSeo5g0CSes4gkKSeMwgkqecMAknqOYNAknrOIJCknjMIJKnnDAJJ6rms8E/BzllJFhj8j+P19mzgPmtY41tgLNawxlp8e1VdMNz4LRkEZ0qSuaqatoY1zvWxWMMa68lTQ5LUcwaBJPWcQXC6vdawxuNQxxrWONM1JuJnBJLUcx4RSFLPGQSS1HdV9YS8APuAU8Dty2wP8OvAXcDHgO9r2i8CvggcBu4HHgXuAF4JLDTth4HfnbDGydZtv9Zaft1wDeB7WtsXgXuBh4CHgU8DX2ltPwycN8ZY7mhqfAn4BPC2NcynXeMUcHSN83m0qfdQ0+eLfPN8VntcDzePw+3N4/Noey4T7OOvNMvHgUfGnM8fteawCNwNfJzBvr67uf3H2/NpPe9mgHngQeDzNM9P4J5mvzwEfBn4K+B7l3lsxq3xmdY4HgY+xWP7eqnG3cAnm9s9wmCf3svge+xLdT7YtLXn9A9b87m7ue0xBs+Rn1nDfMapMc58HgWOtPocbub0labGQ8047hmezzI15ta4b1arMc5c5oGZ5nH+cLP9Xk7/ub1o1Fwmfr082y/YZzAIfpDBD+1yQfBDwJ8weLH4AeDDrReJ24drMAiC/96lRut2D69Wo7X9LcA1PPYi9t3NE+mSScYCbG7VeDpwAnhbxxqfBi5Zw3xOAc8efqwmeVxbL3zPHrVvJtg/9yyNZdz9M/S8eAvwI8AvAbua2ieAXxwxno3AZ5vrH2LwAvWJ1jj+G7CrWd8F/NaIx2bSGr/YLN8CTI+o8fymxlRrn/4H4CNDdR4c8Tgt1fl7wD9plp/b1Dg24XzGrbHafO5tHoep1n28D/iLdo3h/b5SjTXum3FqrDaXqWa5PZdX8s2vQ6fNZS2XJ+ypoar6PwzerS1nJ4N3jlVVHwKemWTzuVgD+Ouq+kiz/AUGQbBlkjpVdXypRlU9BPw/4Gkda3xy0nGMaT1qnLGxMHinNrxPdwL7m+W/Bq4eUesqYLaqFqvqfcAHGLz4Lnlpq8Z+YPs61FhtHHcC7wV2tPbpFPC8oTpPXaHOJ6vqA8As8JKmxpMmnM+4NVacD/B1Bo/Jjtb2FwN3dqwx0b4Zs8aK46iqBxg8HjtG9FtXG870HZzDtjA4dFtytGm7D/jOJIeb9o+1+vxYkpc0y/91jTXGHcfxZfpvZfAC/qtJvg78eVW9doKx/Dnwy8C3A9/Wap9kPks1XszgHfKk8yngz5IU8K4Rtccdx0bgz5rrZ7X3TVX9zph1NgOfTHIv8D+qavireyvVaNtUVceTXMTgsT2vqb+0f0bVOsHpP4PPA/4kyZcZHH1cADx56LG5e8Iam/hmI+fUjP3FDI6AnlpVS/vsBIMX5fcn+Rrw5aq6bJk6L2hqfHmN81mtxmrzKeBK4IokT2/25zMYnHZamssmBmH9jfkA71yuRtM26VzGqTH2vhnRb9jwvplIn4NgJZ+pqhcBND8c/6tpf0dV/fRSpyQ/vIYaa5bkaQwOSY9W1fdOcNP2WJ7G4J3K7zH4+vBa5rNU4/VV9fkkE80D+Pmq2pvkOQwOi48v1W7qr/RYtcexpaqOJfnXwG7gdc0R2Dg+U1UvatV4DjCb5M5VbzlGbWBLe05jenjocSha823a/t0aaozjPOAPgdcDf7u9oaqqqfNPq2qlv4GzdI78dcAblxnLavMZp8ZKXgK8ikFwvXZ4f7bmQns+Q+M4rQbwlDXMZdwa62W1fbOiJ+ypoTEcY3AucsnWpu2crJHkyQx+UP+YwQeEE9dp1XgbcFvXGlX17jXOZxGgqk4xeEf/t9ZQg6paanuIwbnUSyets1SjGcsfraVG42TrtNMGBp+DDBuudSHw1db6wlKN5nrUD/akNcYZx7cBP8bp+/SRoTpfX6lO87x4JfD+Vo2J5jNmjRXn0+zPrQw+fF3anw/SnNpaY42nTDqXMWuMs2/W8poysT4HwQHgX2XgB4AHW4fC52KNNzM4b/rmtdTJ4G37m4FPVtUbz1YNBi+Sf7Opdz7wjxkcmk9SgyTntw65z2Pww3P7hPP5Ro1mLNsnrTHUb6ZZfiZw04g+B4HtSaaSTDH40Pnh1vY/bdWYYXB+uGuNccbxCgYfprb36d1DdR7hm32jDvBWBkcSP7PW+YxZY6X5bEnyXAb78YM8tj//L4MPxdda49EJ5zJujXH2zfam7cyqdfiGzrl4Ad7O4Lz0owzOs70aeA3wmmZ7gN/gsa9xLX1qfxGPfaOkXWMRONSxxtI4vrRajdY83gL8RwbnHD/G4CtoX2TwjZGxx8LgUHWpxmEG32p4T8cah5txTDKfdzH4lsNHGXxN8Jeb2pM+rt/R1Pgog3dMf9m0TzKfdo07gOua7SvOZ2iffgH4TeBZDJ4fdzN4Yd7Y3H4auLE1/1cx+CrqQwzOUy89L+4DLm5qzAM3Ay9sxtmlxtI4bgP+eESNoyP26RuBG4fq3MvgG1qjxrJU40irxskJ5zNujZXmcw+DNxX3Du3P1zZ1v1Gj6XvFiHF8U42mbZK5jFtjtX1zF/BTQz8/r+QMfGvIPzEhST3X51NDkiQMAknqPYNAknrOIJCknjMIJKnnDAJJ6jmDQJJ67v8DSf8OztZCTpkAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Compare final metric magnitudes\n", + "\n", + "fig = plt.figure()\n", + "ax = fig.gca()\n", + "\n", + "labels = [f'{round(sig,0)}{\"T\" if est0 else \"F\"}{\"T\" if est1 else \"F\"}'\n", + " for (sig,(est0,est1)) in itertools.product(point_set_sigmas,estimate_rates)]\n", + "\n", + "vals = [final_optimizers[trial_values].GetCurrentMetricValue()\n", + " for trial_values in itertools.product(point_set_sigmas,estimate_rates)]\n", + "\n", + "ax.bar(labels,vals)\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 66ad2c1792c2ef03fd3a951529ada5ddf4ea6c87 Mon Sep 17 00:00:00 2001 From: Tom Birdsong Date: Wed, 24 Feb 2021 12:58:42 -0500 Subject: [PATCH 2/2] BUG: Added .binder folder for requirements.txt dependencies --- .binder/requirements.txt | 3 ++ requirements.txt | 67 ---------------------------------------- 2 files changed, 3 insertions(+), 67 deletions(-) create mode 100644 .binder/requirements.txt delete mode 100644 requirements.txt diff --git a/.binder/requirements.txt b/.binder/requirements.txt new file mode 100644 index 000000000..1c245d63e --- /dev/null +++ b/.binder/requirements.txt @@ -0,0 +1,3 @@ +itk>=5.2rc2 +itkwidgets +matplotlib diff --git a/requirements.txt b/requirements.txt deleted file mode 100644 index f1964dbde..000000000 --- a/requirements.txt +++ /dev/null @@ -1,67 +0,0 @@ -async-generator==1.10 -atomicwrites==1.4.0 -backcall==0.2.0 -brotlipy==0.7.0 -certifi==2020.12.5 -colorama==0.4.4 -colorcet==2.0.2 -cycler==0.10.0 -decorator==4.4.2 -defusedxml==0.6.0 -entrypoints==0.3 -imageio==2.9.0 -iniconfig==1.1.1 -ipydatawidgets==4.2.0 -ipykernel==5.4.3 -ipympl==0.6.2 -ipython==7.19.0 -ipython-genutils==0.2.0 -ipywidgets==7.6.3 -itk==5.2rc2 -itk-boneenhancement==0.4.0 -itk-bonemorphometry==1.2.0 -itk-core==5.2rc2 -itk-filtering==5.2rc2 -itk-hasi==0.1.1 -itk-io==5.2rc2 -itk-meshtopolydata==0.6.3 -itk-numerics==5.2rc2 -itk-registration==5.2rc2 -itk-segmentation==5.2rc2 -itkwidgets==0.32.0 -jedi==0.18.0 -Jinja2==2.11.2 -json5==0.9.5 -jupyter-client==6.1.11 -jupyter-core==4.7.0 -jupyterlab==2.2.6 -jupyterlab-widgets==1.0.0 -kiwisolver==1.3.1 -MarkupSafe==1.1.1 -matplotlib==3.3.3 -mistune==0.8.4 -numpy==1.19.5 -param==1.10.1 -parso==0.8.1 -pickleshare==0.7.5 -Pillow==8.1.0 -pluggy==0.13.1 -prompt-toolkit==3.0.10 -py==1.10.0 -pyct==0.4.8 -Pygments==2.7.4 -pyparsing==2.4.7 -pytest==6.2.1 -python-dateutil==2.8.1 -pyzmq==20.0.0 -six==1.15.0 -testpath==0.4.4 -toml==0.10.2 -tornado==6.1 -traitlets==5.0.5 -traittypes==0.2.1 -wcwidth==0.2.5 -webencodings==0.5.1 -widgetsnbextension==3.5.1 -wincertstore==0.2 -zstandard==0.15.1