{ "cells": [ { "cell_type": "markdown", "metadata": { "cell_id": "9F60A887F43F4A678C6B6447F3327B83" }, "source": [ "# Implicit functions in pytorch\n", "\n", "*Thomas Viehmann, tv@lernapparat.de*\n", "\n", "Sometimes, we do not know the mapping of functions we wish to apply, but only an equation that describes the mapping. In mathematical terms, we wish to apply $f : \\mathbb{R}^n \\rightarrow \\mathbb{R}^m$ of which we only know that $F(x,f(x))=0$ for some function $F : \\mathbb{R}^n \\times \\mathbb{R}^m \\rightarrow \\mathbb{R}$.\n", "\n", "This is the realm of the [implicit function theorem](https://en.wikipedia.org/wiki/Implicit_function_theorem). \n", "Under reasonable conditions (smoothness, the right sort of nondegenerate derivatives), you have the following:\n", "Given $x,y$ such that $F(x,y) = 0$ there is a neighborhood $U\\ni x$ and a function $f$ such that $f(x)=y$ and $F(x,f(x))=0$. And if $F$ is nice enough, we can also compute the derivative of $f$ at $x$, namely $\\frac{df}{dx}(x) = - (\\frac{dF}{dy}(x,y))^{-1} \\frac{dF}{dx}(x,y)$.\n", "\n", "There is an example below.\n", "\n", "Nice. As computing the entire matrix Jacobian with respect to $y$ is not that straightforward in pytorch, though, we will stick with just the scalar case.\n", "\n", "Let us not be lazy about it and get our hands dirty. First import everything." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "cell_id": "5802BF8F7F724585BEFEA2528B60BBD2", "collapsed": true }, "outputs": [], "source": [ "import torch\n", "import numpy\n", "from matplotlib import pyplot\n", "from mpl_toolkits.mplot3d import Axes3D\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "F3ABBD4BB63E42C18685C2BDDDA7421E" }, "source": [ "Now we can turn to the implicit function. In evaluating $f(x)=y$, we need to actually search for the solution $y$ to $F(x,y)=0$. We do this by searching for a minimum of $(F(x,y))^2$, but we also need to provide a starting point $y_0$ near which to look. (And indeed, in the circle example below, we will have two solutions and it would typically be a problem when these are close.)\n", "\n", "We use the pytorch LBFGS optimizer for a limited number of steps. Note that LBFGS also has other stopping criteria, so this really is a bound, in fact, achieving the other critera can be considered success, hitting the maximal number of iterations might be considered failure to find the minimum. The scipy optimization documentation has more details. In order for LBFGS to work, you need to provide a function that re-evaluates $F^2$ and the gradient as the `closure` argument to the `step` call.\n", "\n", "~~As we need to call `backward` on `F` for the gradient computation and I ran into locking problems when doing so in `Implicit`'s backward method, we compute the required $\\frac{dF}{dx}$ and $\\frac{dF}{dy}$ in the forward. In fact we precompute $-\\frac{df}{dx}$ and save it to the context `ctx`. In the `backward`, we wrap the saved result in a variable and multiply by the `output_grad` to do our step in the backpropagation.~~ PyTorch is ever improving - for the 0.4 update, I moved the derivative calculation in the backward.\n", "\n", "We compute the required $\\frac{dF}{dx}$ and $\\frac{dF}{dy}$ by calling backward in the backward. I don't know if all `detach` are needed or if some of the tensors are pre-detached, but I'm going to play it safe here." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "cell_id": "3DB49B74E00149CE847EF4F7F5E625AA", "collapsed": true }, "outputs": [], "source": [ "class Implicit(torch.autograd.Function):\n", " @staticmethod\n", " def forward(ctx, x, y0, F, max_iter=100):\n", " with torch.enable_grad():\n", " y = y0.clone().detach().requires_grad_()\n", " xv = x.detach()\n", " opt = torch.optim.LBFGS([y], max_iter=max_iter)\n", " def reevaluate():\n", " opt.zero_grad()\n", " z = F(xv,y)**2\n", " z.backward()\n", " return z\n", " opt.step(reevaluate)\n", " ctx._the_function = F\n", " ctx.save_for_backward(x, y)\n", " return y\n", " @staticmethod\n", " def backward(ctx, output_grad):\n", " x, y = ctx.saved_tensors\n", " F = ctx._the_function\n", " with torch.enable_grad():\n", " xv = x.detach().requires_grad_()\n", " y = y.detach().requires_grad_()\n", " z = F(xv,y)\n", " z.backward()\n", " return -xv.grad/y.grad*output_grad, None, None, None\n" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "331983CBCD354721B4F41718553388C3" }, "source": [ "So now that we have a cool new `autograd` function, let us apply it to an example (and I must admit, I just took it from the Wikipedia page, the application I have in mind needs more context).\n", "\n", "The unit circle in the plane can be described by the equation $x^2+y^2 = 1$ or, equivalently $F(x,y) := x^2+y^2-1 = 0$. So let us define $F$." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "cell_id": "1A05EB10B00F41B6876792D5E32F3608", "collapsed": true }, "outputs": [], "source": [ "def circle(x,y):\n", " return x**2+y**2-1" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "1F30165C0DD34AF28F2DFC0732845D1F" }, "source": [ "Everybody loves pictures, let us plot $F$ (the blue grid). In black is the circle $F(x,y)=0$." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "cell_id": "9CF19E5055C841DDA0E7A7A7955C72A5" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADuCAYAAAAOR30qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzsvXd4HNX1Pv5uk6wu2VaXZVvd6lpZtiyBwfTQQu8x/ZtACMkvED6AAWMSEyCQkEASktiBmE4g1IADBmyDZGFZvctqVrF6L6ut8/tjudezs9N2tbIse97n0aPdmTt3ZnZ33jlzznvOUTEMAwUKFChQsPBQL/QBKFCgQIECOxRCVqBAgYITBAohK1CgQMEJAoWQFShQoOAEgULIChQoUHCCQCFkBQoUKDhBoBCyAgUKFJwgUAhZgQIFCk4QKISsQIECBScItC6OV9L6FChQoMB1qOQMUixkBQoUKDhBoBCyAgUKFJwgUAhZgQIFCk4QKISsQIECBScIFEJWoECBghMECiEr8CisVitsNhuUOtsKFLgOV2VvChTwwmazYXZ2FiaTCV5eXlCpVFCpVNBoNFCr1fS9SiVL/aNAwSkJhZAVzAkMw8BsNsNgMMBms0Gr1UKj0VAL2WKxOIxXq9UOfwAUklag4HuoXHy0VJ5DFVBYrVYYDAaYzWZoNBpqAet0Ot7x5LfG/s2pVCoHglasaAUnKWT9qBVCVuAyGIaB0WjE7OwsAFAiBSBKyHzzsP8TKFa0gpMQsn7AistCgUuwWCwwGAywWq0OREzgyg2eTeLs7RmGcXB1qFQqGAwGBAQEKFa0gpMaispCgSwwDAODwYDJyUnYbDbqovA0CNlyLeTKykqYzWZqmc/OzsJsNiuqDgUnFRQLWYEo2EE7hmHmjYjFQPZHyJkvYMgmcrY/W4GCxQSFkBUIwmazwWAwwGQyUaKTA4Zh5pUMpVwdhKgV2Z2CxQaFkBU4gWEYmEwmtLe3IzAwEP7+/m4TWfLWPQCApm3neOS4hI6DS9KK7E7BYoTiQ1bgAKvViunpaczMzGBmZgYWi8UjhEWI2V24egx8vmiVSkWtaJPJhNnZWRiNRphMJlgsFsUXrWDBoRCyAgB2i3J2dhaTk5OwWCz0Ud9ms7k9pyetY3KMcwGXpAnJW61Wh4Ch0WikAUPiClGg4HhAIWQFsFgsmJqagsFgcAiKuUvIQtbsXKzk46XoYFvRHR0dOHr0qGJFKzhuUAj5FAbDMJiZmaFSNq1W60B87hIym7CIdUz+z4WUjwcRskl6dnYWVqsVgD3AqcjuFMw3FEI+BUGCdhMTEzCZTNQ9wcVcCJmPeOfiuliIwBsJIgq5Oogv2mg00j+z2awQtAK3oRDyKQabzYaZmRlMT08DgKiu2FVCNhqNqK6uRnFxMQDg4xtjMTAwAJPJ5DDOXSv5eJOckKpDKHmFuDoUK1qBu1Bkb6cIiFVsMBgAiBMxgVqtdpKNCc199OhR9PT0ICEhAeHh4cAXXyI4OBhjY2Po6OiA2WxGQEAAvrg9CefubEby1j0uWczEt3s84YqeWpHdKfAEFEI+BcBXlU0ONBoNjEaj6JiZmRk0NTXB398fubm50Gg0SHn8SwBAWFgYwsLCANgt86mpKYyNjdFtk7fuwcc3xiI4OBjBwcHw8vIS3M9CkZa7+3WlTodS7U4BgULIJzG4VdlcTXsWc1nYbDZ0dXVhYGAASUlJCAoKchrDtoLVajUCAwMRGBgIoJmO4bOiCUH7+/s7+LYXwkL2JISsaKvVSoOHADA9PY2AgABotVrFij7FoPiQT1JYLBZMTk46SdlcgRAhT05Oory8HDabDbm5ubxkLKSq4GqTT/9zNZKSkrBu3Tps2LABq1atAgB0dHSgpKQEBw8eRHNzM0wmk5Mver5xPFLA+SzkxsZG2n1FSV45taBYyCcZGIbByMgIpqenERQUBK3W/a+YS8hWqxUdHR0YGxtDSkoK/P39nbbJ2L6fvm7ado6o2oKsJ5Y024qOjY0FAJhMJoyNjaGvrw+1tbWwWq2iVrQnMd+EzAeVSkUliOxgoc1mc7CilTodJycUQj5JwK7KNjY2hsnJSYSEhMxpTjYhj46OoqWlBREREdDr9aIXf+PjZzu8FwvgCZE2gZeXF8LCwtDV1YW0tDR4eXlRX3RHRwempqag1WopQUv5ol3BQhAyYHcHsW8ycgKGSrW7kwMKIZ8EIFI2s9kMtVoNnU43p5RnAqKyaGpqgsFgQHp6Onx8fCS3YxMZ2woWg5Tqgl2CU8iKluuLlouFImSGYUSPV6l2d/JCIeRFDCEpm0ajcXi8dRfj4+MYGRlBYmIikpKSJC9mtruCDbYVzEe6XNeFEIR8p8SK5lN0zLcVfaJAkd2dHFAIeZFCTMo2V0I2mUxobm6G1WpFYGAgIiMjZW9bs2Wj2/sFhC1lV9UhnrCiF8pC9gQU2d3ihELIiwxypGzuEjLDMOjr60NXVxdWr16NoKAg1NfXuzRHxvb9aNh6luB6PsJlW89iro25qAvcsaIXMyHzQUp219/fDx8fHwQHBytW9AJBIeRFBBK0s1qtojI2dwjZYDCgubkZ3t7e0Ov10Gq1VGYlB0LuCjbECJdN0nyk7WlCkGNFT0xMwGw2IzQ0dN4VHQsBLkGPj49Dp9MpVvQC4uT5dZ3EYNefYBjGqSobF67UoGAYBl1dXaitrUVsbCxSUlKoVM7VWhbEXbFm21cOy4W0yHwQqwo33/pbYkUTXfSyZcsQExMDwFkXzVejY7HDYrFAp9Px1ulQqt0dHyiEfAKDBO0mJydpXzs5FppcIp2enkZFRQVMJhP0er2TTI5oYqXAto4JKXMJVagMp5hOmb1uIWpZAEBAQABiY2ORmZmJgoICZGdn0+zCyspKFBUVobq6Gp2dnZiYmJizusVmsy2Y9WmxWHh160q1u+MHxWVxgoIrZZPbYBSQfry32Ww4cuQIhoeHkZycjICAALfmYYMdzKt/bBNSn/haMEDHlcKJjXHnWOYT863o4GqQjyeIK0wKiuxu/qAQ8gkGd6qyuYKJiQk0NTUhLCwMer1+Xi5+i8WCxsfPRsrjX4r6jOWU4WST+olY7c3TuuiFJGQhC1kOFNmdZ6AQ8gkEq9VKG4t6moitViva2towNTWFtLQ0+Pr6ytrOFe0xSa0eHh52GFPyi1yXiYZN1slb9+Ddq8Jlb+spuKuykLKiJycnodPpeK3oxWAhy4Eiu3MPCiGfACBSttbWVvj6+mLZsmUe+VESQhkZGUFLSwuio6ORkJAwLz/48fFxNDc3IyIiAvn5+dBoNCjPmYH+mQPIf74ML5/vC41Gg6CgIISEhCA4OJhuK+baIOuvercfpfet9Phxi8FTsjdXrGhfX1/YbLYFI+b5LqbE/s+W3ZWUlCAvLw/AqW1FK4S8wGBL2Yj14CkSMBqNaG9vh9lsRlZWFry9vT1wxMfAto5Pe7EKB3+ZBx8fH2rl6J85QNcXFBTAbDZTEurs7AQAvHNFKK75zyCSt+5Bw9azeEmIuDfynjvosQ7WcjCfOmQhK7qvrw8zMzM4cOCAoBV9soBN0FarFWq1+pS3ohVCXiDYbDZaYlGlUkGr1VLtrydgsVhQVVWFVatWISwsbE4/YDFi2nnuEkRHR+OCV1p561xw06JDQ0MRGhpKXRJxcXH44vYQnLuzGQcOHICXlxdu+mQcgP1mpdPpAAD/uSYSV7zT6/Y5uIvjdeETK5rEENLT0+etRseJBrayRMjVcapUu1MI+TiDXZWNFJEhPyStVksz8NyF0WhEc3MzLBYLMjIyvi8I7z6IhI7tW7RYLGhtbQUAZGZmYsmSJQBakbF9P6/sTahWBXltV3k047bPDah+eC3wyTf46IYVKC8vp+nbpAegq62f5oKFyNRjuyrm4ot2Z78LRWjsGy8fTqVqdwohH0cQKRtxS3ADKHOxkBmGQW9vL7q7uxEfH++yVE4IXEIeHh5Ga2srbv/CfuOwk7Fd9paxfT8ytu9H/WObAMCJfKW0x8lb9yDzyW/sY5KTAdg/s4mJCTQ2NtKxyVv34NPNqxESEoKgoCCPBaK4WAgtrVhgbT4r3c1FYTFXSBEyF+7I7tiGz4kMhZCPA7hSNqEfh7uETPra+fn50bTnwcFBj1R8I4RsNpvR0tJC/dH44junQkKElMWIzBXtMdl/cHDw95b+FJXTXbirHe9fa0RTUxMAOAQLlyxZ4tGg6PGEq8E8T1nRi4mQ+SBlRXt5ec3bjduTUAh5nkGkbCRoIXaBazQalwjZZrOhu7sb/f39Tn3tPFWCU61WY2hoCD09PVi5ciXCwsKoFSuEtF/v5V0uV3vMB/YFR+a5/O2jaNp2DiwWC8bHxzE2NoajR49idnYWfn5+lIACAwMXTT3kuaor3LWiPSl5cxWeIGQu2L+XhXTHuAqFkOcJ3Kpsch6ZtFqtbBKdnJxEc3MzQkJCkJub63QRe4KQSdq21WpFdna2LP8ksZKlIOUL5q6/7K0eh/Vcv/SyZcuwbNkyAPbPfnp6GmNjY+ju7sbExAQ0Go2DlShHcbJQhOxpYpRjRZPzHBgYOO6KDpPJNO/7WyyEvPhDtCcgzGYzJicnMTs761KDUTkuC6vVitbWVjQ3NyM5ORlxcXG8FtVcCXlgYACVlZXw8fFBfHw87wXDR7xyq77xgevO4FrTxffm8M7DHadSqeDv74+YmBikp6ejoKAAOTk5WLp0KSYmJlBVVeVQg2JyclLQzXK8L2TyJDWfIFY0qdFRWFiIVatWwcfHZ95qdIhBykJ296lqMUIhZA+CBO2mpqbAMIzL2XZSJDo2Noby8nJ4eXlBr9fzNhmVO5cQTCYTampqMDQ0hJycHJqoQEAIl/iP+Qi47tEzAUhfSFINUOXMITYXG15eXggNDUViYqJDh2uGYdDW1obi4mKUlpbi8OHDGBoagtlslrVfT2MhE0KCgoIEO4AfOHBg3irdmUwmQUKeKxkvtiJHCiF7AOyqbEajkUZ2XYVQRTPS166jowPp6elYsWKFrBoLrhAyKU5fWVmJyMhIpKamQqfT8VaOI2QsRsp8hMqVwrGXyan6JnRxuUrewDErceXKlcjKykJhYSEyMjIQEBCAoaEhlJWVYWpqCnV1dejp6aGlT+cbC0XIVqvVIajHZ0XPV6U7s9ns9AQmVXzKFSwmCZxCyHOE1WrF1NQUZmZmAECyVrGrGBoaQnl5OQIDA5GVlSWrySg5DrmEbDQaUV1djbGxMej1eixfvpyuYxOylDuCrCfExUeUfDpkKeWFHLhDylwsWbIEERERSElJQX5+Pvz8/BAVFUVbWhUVFaG8vBxtbW0YGRnxSNCUi4UKrpH6KWLg1ov2lBXNdVl4iogJFgsZA0pQz22QoN3IyAiOHDmCtLQ0j37xJpMJhw8fBsMwbqU9y6mJzNYuJyQkYOnSpZLzCEnduK4MAqmu064oLwpfqHTqLCK0L09czCqVCiEhIbRONMMwMBgMGBsbQ19fH5qamuijPltyNxcslIVssVhk3+wJPKWLZhOy3HZerkAh5JMcbCmbl5eXx+pPAPYfT29vL+1rFxoa6tY8Uj5kg8GApqYm+Pr6Uu2y2DyEcNnZeARylRXukKTURcmd09OkzIZKpYKvry98fX0RFRUFAFRyNzo6iu7ubhiNRiq5CwkJQUBAgEsEu5CE7Akdsju6aLJvPst4rt/hYuuLqLgsXACxkIi/TK1WQ6fTeaz+BGmNMzIygpycHLfJGBAmZIZh0N3djZqaGqxatQpJSUmiFyLbQpbTUZphGPT396Orq4sGN+W6E8Ss6E9uWilrDjnzeRJarRbLli1DQkIC1q5di4KCAsTHx0Oj0aCzsxMlJSX47rvvZD/Gz4fsTQ7mKzFEji96amoKKY9/CQAovW8dbDbbKaWsYEMhZBkg9ScmJiYcgnYkd36uwR5CktXV1fDz88Pq1avnLJTnI+SZmRlUVFRgdnYWubm5DiUwhaBWq3H+yy0Oy8Ss4cwnv6H+1ZaWFhQXFwsG9tiQIm2VSoX9d2fQMXIt5+N9YatUKgQEBGDFihXIyMigkjs2ARUXF6OmpgZdXV1OkrvjIXvjw/H0XbN90Xl5efhZ0bHzJb5oAPj4xliPKDoWUwGmxXOkCwSbzQaDwYCpqSkAnu/gQfraEZL08fHxiMXNJmSGYdDZ2Ym6ujrEx8cjISFB9sVHfsxiygryeue5dv/p1e8NYNWqVcjOzkZhYSEA4Osfp+L9a+2P+SUlJWhoaEBfX59s5QX5zKUeZ8XmWyhwg2H5+fmIjY2lmvLi4mIcOnQILS0tMBgM86r5FcJCpU4Tyxiwf1/EigYwr70LT1QohCwAErSTK2Vz1Uq22Wzo6OhAQ0MDEhISKEm6oo4QAyHk6elplJeXw2KxIDc31yG9Wg7O3dnstIztujAajQCAf19pbwnFdWsQMoyKikJqaioA4ObPphAaGorJyUkAwD/P86EX26H71ztsxyZdMfUGG2JKjhMBarUaQUFBDjet9PR0+Pv7Y3Z2FvX19Thw4ADq6+tx9OhRzMzMzLvkbiEImXwnb/5wKe9yTyk6FpMPWQnq8cBqtcJgMMBsNlOSFAPJsJPrZpiYmEBzczOWL1/u1NfOUzWRVSoVZmZm0NDQINrIVA4+vD6GdznbSk5JSXFYxw6q8VV9K3yhki4rKCjA5OQkRkdH0dbWhn+e54PbPrcXYhoeHkZQUJDDRcV1gUgFfuYz0OcpEMldd3c3MjIyoNVqaX2OxsZGzMzMONXn8KSL4Xi6LNjfX+l965D33EGnddzvyV1Fx2KDQsgscOtPuJryLEXIpK/d5OQk1qxZAz8/P8G55oLJyUk0NTWBYZg5NTIlhPvDN7tRsyWOLjcajXj9kmDc+PEYAGEpnBypG7nwyMW2cuX3bZo+t68fGBhAc3MzjEYj/P39ea1kQrSuNk09EUGCxRqNBkuXLqVSRIZhMDMzQwsoNTQ00Ep45G+ukrvjYUlyXUqk/6I7Mjc5io7IyEgnY+FEhuKy+B4WiwVTU1Mu158AAJ1OJ5lqOzIygrKyMvj4+CAnJ4eXjIG51aCw2Wxoa2tDc3MzUlJS4O3tPeeAxjc/zQQAWlazr68PVVVViI6OFt2u8fGzRdfLdSVc9lYPNmzYgBUrViAwMBBnvlQHACgqKkJVVRW+uD3JYR4xv/KJ6L7gQkhloVKp4Ofnh+joaKSlpaGgoAC5ublYvnw5JicnUVNTQz+TI0eOnJB+Vj7/PnE1CCUPuQKuoqOgoODYDX6R4JQnZGJ5TE5O0ovBVUtBzKo1m81oaGhAV1cXMjMzERMTIzq/uxbyxMQEysrKoNFokJOTM+fHNWIdq9Vq/PtKuwVSW1uL0dFR6PV6WlmNPdYViJEiH4FqNBqHcyooKEBcnN1qJ8cHAC0tLRgaGnL6DN2tkXG84UqpSJ1Oh+XLlyMxMRF5eXn0M1GpVA5+1sOHD2NwcHDB6nOIqWLIDbZp2znz8uSyWArTE5yyhEzqT0xMTMBkMrldfwIQtpAHBgZQUVGBkJAQVqsjcbhKyERe1tLSgtTUVKxcudJjMp+aLRuhUqlo+6QbPhrFmjVroNVqJYsMsaPn7kjduAR6/sstDu4KIi8j1hDBRa92YHBwEIcOHaKBMQC0OQB3vyciMbtLINzPpLCwEFlZWQgMDMTIyAjKy8tRXFyM2tpadHd3U534fNYLZt9chW6K++5Kd3jvSSwmMgZOUUImVdkI0cxVysYlZKPRiJqaGgwODiI7OxsRERGy53eFkMfHx1FeXg5vb29BN4g70XlCriaTCa2trTAajSh/YIPDOkC6yBD7IpRbZIgPZMxZf28QPW4y7rK3epCfn49169bh8rePAgDq6+tRVFSEyspKHDlyBAd/mcd7PCcbvL29ER4ejuTkZKxfvx75+fmIiYmBxWJBS0sLioqKUFZWBpPJhOHhYY8lOQHy3EgAHGIv85FduZhwShEyCdp1dnaio6NjTlYxG4REGYbB0aNHUVVVhaioKKSlpblceFuOD9lqtaK5uRltbW2i1d/k1LMQQ2VlJcLDw+Hv7w+dTida3Y0d2OMjZkBekSEhcvz8tkRZx8x1c5Blubm5NItOrVZj3e9LAQCvXRxExy/UI/3xBAkEsiV3CQkJ8PLyQn9/P0pLSx0kd6QZrysQ+h759OX/vjJs3orTL7a0aeAUImSiyZ2ZmYFarcbs7KzHviydTgeDwUDTQLk+VlcgZSGPjo6irKwMfn5+yM7OFi0I406AkE2mt38xi+XLl/OSulAaNdmeG9RzRzvMhpDsje+9WOF6kkVHxmVmZtIMwPTf7MOBAwdoyU2bzbbo6um6CpVKBa1WCz8/P6SmpmLDhg3Iy8tDeHg4DAYDfbKoqKhAR0cHxsbGRG/yXKtYyOIly81mM9Y++91J/ZTiCk56QmYYBrOzs5icnKQlBr28vDxWYJthGAwNDWFoaAhxcXGStSGkIETIpCZyZ2cnMjMzER0dLXlDcYWQSQ0KwJ5VRwhXyMqWspL5IPdxlO/iJAkqQmTLV2RIzvzkkZ6Mv2X3NK54p5fWFSkuLkZlZSU6OjowPj5+wikXPAFuUgipzxEfH4/c3FwUFhYiMTERWq0W3d3dtD5HU1MT+vv7YTQaJesXc5+OyPsbPhp1WO5pLKa0aeAkJ2SLxYLJyUkYDAYHKZuXl5dHHk+npqZQXl4OhmEQGBjochYcH/gIcHh4mNZElhscJHPJIWSTyYTa2lqcs8PewZldDzn7qSKnOcRcFwTsoJ4rELsoicTN091E+HDhrnb4+fmhoKAACQkJ0Gq1tFgQyRDztHJhoaxxqSw9vpZYer2etsRiN7394vYkTE5OSsYMFPDjpCRktpSNYRinovFzTb4get+mpiYkJycjPj7eY8EQdtcQs9mM+vp69PT0ICsrC5GRkS63hJKy6Pr7+2mXEMDRyiWvSdYcX0CPD0Lr5BQZkjtWLikL+afFHqXZ1nLK419SMiLFgki1Mq5yoaenZ05pzgzDLJrCQjqdDqe9WIWLXztCl5Xetw4A0NbWBsDun5cK4L58vq/DOk9iMXUKITipCFmulE2oVZIcjI+Po6ysDFqtlva1m2vwjA+Dg4OoqKjAsmXLkJGR4XKBekDcZUF65w0PDyMnJweb/lbPO44vC48NviJDZJ9CROqO8uLLO1OcvjO5bgmhfQqNB+x9AXddGMC7jmSIsZUL0dHRMJlMaGxspD7X9vZ2SZ8rGwtV6c2dOhZ8cjaSlJGVlQUAuOmTcYdtyI2LoGHrWbj1fzNzOXRRLEb//0lDyK5K2VwlZYvFQpUNaWlpiI2NpfN7ulPI7Ows+vv7kZ2djfDwcLfnFyJk0lGa9M7TP3OArhNyQ4j5jNnrin6WjfLycrx9ud3tkbx1Dw4fPiwr8CZWFIjv+5JrRQu9FxvPjdCLKUDUajVCQkKwevVq6PV6FBYWIikpCTqdzsnnKlYIZzEUp+fzFYsFWdnBvfz8fFz57z4AwFuXLUNRUREA4LOb4+atJdZis5AXfS0LYhUT4b+n608Adh9ua2srYmJikJiYOC9fMsMwGBgYwJEjR+Dl5YU1a9bMudgLl5BJbziVSoWcnByHc2eTK7crCLvYDxfcFk4tLS1IS0tDUFAQmrJVSN66hz7W/usH/jhw4ABteVT10GnI+u23smpeCGG+alqwCZnP7SG2LUlz9vPzQ0yMvTATKYQzOjqK9vZ2WK1WBAYG0s4ivr6+C9pPT6p9k9wgqtjns2bbV/R1Tk6OvaTtnhIsWbKEtsQCQGtzhISEuF2fYzFax8AiJ2RuVTZXiFKn04m2HweO9bWz2WySfe1IAM2dC8poNKK5uZmmPdfX18tqOikFNiEPDAygo6PDqS0U1/Lla8fEJmO+Fk4E//3RKidNNJsg8/PzYbVaacujnp4eOo4Eg/z9/Xm/x7P+3oDdt8QLniu3d59UTQupLD2GYZweuV0lZjb4CuFMTExgdHQUzc3NmJ6epsHm0dFRBAUFHTdrmdtxmg0hdQR3vZQriM9qzv1dCQB7aVZuS6yxsTH09PTMqSXWYvQhL1pCnp2dxZVXXonXX3/drUw78uPna85JJGCdnZ1YtWoVvYjEQFo5uUKipFhPV1cX4uPjqXaZWO/u+I3Z0Gg0MBqNqK2t5bWKCfgIlku8fO4Jm82Gzs5O+v6iVztQsyVW8pjYVczwsf1CPXdnM965IhRTU1NYsmQJgoODceGudgDHiOCCV1rRtG01AGEXhVyLWsqi41aVE9qXq8RMwK7URvY3ODiI9vZ2Ws1No9E4WIvzlUDB97sVC4BKfSZyPgsyP2nNRUAkd+RaYBgG09PTGB0dRWdnJyYnJ+nnQhrLztfnshBYtD5kb29v9PT0uF08hFjIXKKZnZ1FdXU1RkdHkZOTI4uMAfsPyRUJFNnP+Pi4UyLJXCq+sTE1NYXu7m6Eh4cjLS3NiYzZjUvZ4LovAMdHwIzt+2mnE7Kc+znyPTJKJYVc859BFBYWIjU1lT5Cv3KBH0pKSvDGpSFOc8gN6olhLhI6brBwLskNRI7p7+9Pq7kRadn4+LhD6yd2HQpPgO1Ddvc8pJ5I2DdCobRpPhDJHbclVkhIiGRLrMVmHQOLmJDn+mHzaZFJX7sVK1ZgzZo1LvW1k1OCE7ATVU9PD91PSkqK0+PiXGV5ZrMZtbW1mJycRFhYGG+zVL4nAzbYy1+5wI8qBaofPh0AkP98GZKSknDpG11O27AJXm7pS/Z6Hx8f+gi7YcMG5ObmIiAgAK9eFOgwh5DMTErFwQZ3DHmf9dtvecez5+d7Pxdi5rq8dDodQkNDnVo/kToUxcXFKCsrQ1tb25yCYtf8ZxC5vysRLQREzo39Wuy7FXpyYI/5360JbvWOFGqJZbPZHFpiebIux/HConVZAICfnx+mpqbc6obh5eVF++SV/CIX+c+X4Qf/akPlg4Vu+W7lkKjBYEBjYyP8/Pyg1+sF/XZzIWTy2Lt69WpotVoMDQ05jeEjYzHfMJH1EVnXB9dF47K3eujnznVtCNWyYD/m87kXhHy7Op2Odsho2raOrifHQ7pphITYrWibzea0L/ZxcCFEqG9dJi/9nY+Y3En7vbbWAAAgAElEQVSGkFJZkNZPJAGJdEEfGxtDb28vmpqaoFKpHB7nhdxeclLYhSDlzhCrY8GG2WyGr6+v5P6kwP5cSP3jmZmZOTcKXggsakIODw/HwMCAW4RMXBYdHR0YGhpC0c+yUfhCJbKfKpLV7p5vPiELmVjFvb29SExMlOz27A4hm81mNDc3g2EYZGdnw8vLCxMTE4JWEx+JskmZTaqbP53EJ0uPYmBgAElJSd8TQo9oph47os6GlO+VfUGz1/E9Edl9yudQPyMJEpWUlECn0+GTm1bi4teOuGSxsvd/3QfDaMo5tk5sHq5kz51MNVdlbyqVCr6+vvD19XUIihE1R2dnJ0wmEwICAhASEkIzMbnHUlxcjFv/NyNoDXMhZhGLkTF7u6Zt56Curm7e/L9LlixZdGnTwCJ2WQBAREQEent73drWZDKhv7+ftjkKDAyUlRIsBCESnZmZoV2l9Xq9JBkDrvuQSRJJaGgo0tPT6Y+cbx45fmP2MpJ9dfFrR6DX66l1JjdTj4+A3KlpQXTIfBc38TOS5QUFBcjMzERAQAA+uC7aYb6+vj7aokuuPM4da5eMZf+x5+N7EvCE7E2r1dKi9Td9Mo7bPjfg6vcGKBn/87xj8jYpN4eYqkTMwhbyF3NhNpvn1YpdjD7kRW0hR0ZGYmBgwKVtrFYr2tvbMT4+Dl9fX6xevZp3nNgjPB90Oh1NSgHsVnFXVxf6+/uRnJyMwMBA2XNptVpKGmLgs4rZ4BIy21XBpzfmaop7e3vR1dWFXRcGYPOnkw5kwc3QE0uXdictls/lwFU+iMnXvL29ERERgYiICAA9qHv0TKT9ei/O+GstXrs4iCZn7LkjGVNTU/Dz83O4gA/dvx7t7e24+r0Beh5C5+fqebG35d++1a2bltz94vM92PuTNKr9JdlyfX19Dr355GiP+W5WYlI39jKFkJ2xqC3k8PBwWqFMDkjpyiVLliAnJ4c3Si1W11cMbJfF9PQ0ysvLYbFYkJub6xIZA/JcFkJWMRtsQub6jaVuNq9fEkwVIOQC5VrXNVs2ij5VyMnkkpPlRaxgdt0E7hjua+58Wq2WrieWI2B3E7ADZGQ8CWKyLT6+4+Xbr6skSfbxxe1J2HNHMp1Dzh/fMXKtcr7PhWRpbthgbzxQ9LNsTE1NoaamxuH4D92/3mlbKe2xELhW9nwT8mLEoiZkuRay2WxGY2MjLV0ZExMDtVotKBtyh5SJ7K2jowMNDQ1ISkpCXFycW34sjUYj2qOvrq6OplaLyfLkuD6Ezu/Gj8eoAkSj0eDbe7IcxvNplPnmFAukCRE2H8me9mKV4HxyFRzc/QLAypUrkZ2djYKCAqxZswYA8J9rIlFdXY2hoSGa7szdj5SaQg5Jc5exfch87g7ue779SFnWfPsE7BX+Lnq1gybDlN63Dl/cnkQLBb1+STBaWlpol2ju/vi+Azk+6fm0YhUL+TgjIiJCkpCJJRkcHOxUulKsnoWrpGw0GjEyMgIA0Ov1bgUaCbRaLS+RDg0NSVrFbJCbjtDxs61bdpIMadfEnsdqtUpa1ez91D+2yWGdkFUlRwfsCsm4KjsjY0mADADS0tKQmZmJ5cuXIyQkBKf/uRqAPfj1n2si8fWPU2XPz4aUP91qtfK2qXJXESFnOxJ85Y7lFgpKT0+3lyP9UwUd88F10aLBSz5ydsd95SoWa9o0cBL4kIVcFuy6DXz+VcC1ehZCflKbzYaOjg6Mjo7C29sbq1atcvk8hI6LwGw20xRuoXMRAnk0l/IbkwJDfFI4OdXsuP5nNuQqK+QG2aSITWou9nEIjWUYBhqNhj6BNG07BzabjaZ9v3FpCC2u7q7umAu2m4QLT5CvmKuIq4DgjuPqsusePRPj4+MAevD6JcEO85Xet86BFPmCe/PZWBVYnGnTwCK3kMPCwpx0tgzDoLe3F5WVlYiIiBDtayfVOUTsURwAJiYmaCnOnJwcjxWGYbssiFW8bNkyWVYx3/EKNSO1WCxoaDhmkQmlShNC5iu1yQZ7+9HRUd5Oz0KQYzl9eH0MANeqvElpj4VcAHz92EhVt7i4OErGZb/Kx547kvHvKx1dR+yGoXJ1uXw+crnbso9dCtztuccuBPZnRlKcAbt7i42Ojg7aoOCTm1bSa5R9fCaT6aRKefYUFjUh63Q6WK1Wejc2GAyoqqrCxMQE9Hq9Q+cLoe2lsuv4SNlqtaKlpYVWNouNjfWo5pH4fuvr69Hb20vLcLoCKTfF6OgoysvLaTIF3zZk7CWvd+LsfzTSZUKBPPb7whcqaT+2qqoqh5oXUiQjN9DHXcYHVwNsxGImhCzlSiFpvZmZmQ7rCv5UgbRf76Xb77srHUajkXefbJT8Ilfw2KSWscF2EfEFAMn2Jb/IdUib5puDbCfne2nadg5VpgBAQEAACl+oBACH5qlTU1OKwoIHi5qQVSoV9W92dXWhtrYWq1atQnJysqz6rnJ763FJOfupIlz+9lHk5OQ4ZRp5olD90NAQDAaDW1Yx4GwZs4mS+KY3/qUGmZmZtFedEMkKBe+448n/necuwd6fpAEA7fQcFxdH97Pz3GOSKm7xdrFA34Gf651kb1JkK5bAwAXXnbLhj+W47K0eh2VywBeMA4Az/lqLzCe/cSBHvnoU7tzYhciXfUxC4Gqf+dxLfE8Q3PVCUrcz/lpLX69btw4REREwGAw4fPgwBgcH56Vf4UJ1XvEEFrUPGQCWLl2Kiy66CH/4wx+g1+tdchvIrT8hhMwnv3EgKVLxzd1HMeIrtlqtWLJkictWMeBMxq9c4Idbdk8jY/t+fHtPFpqbm7H7lnhc8Eor8p476DCWL+1ZDGQ8eTR/84dLkZycTC0f4o8kAc6mbfasugq9ATlPF2P9Hw7h5fN9qQ52aGgIwcHBvL5drsUjlnTAZ9nJsajZVt7w8DANYAn5VeWCva3NZqOBNJJZyMa635eK3nBckZyJKTDI9tzi9Fz/utD+5YA7ll3pj5Q9iIiIwOjoKI4cOUKtZlLhLjg4+JSTxS1aQjYajXjyySdRU1OD3/72t0hNdT3yza5nIYaRkRHsPHcJbv/COVmDHSQj0jd3CJkUwV+5ciXCwsJw6NAhXj+mEPjcBzVbNkKj0aD0vnXIe+4gWltbkZ6eDh8fH9RsiRbtDsJnMQsFNnOeLgZgj8QT1D+2CalPfO10URI1A7ngCRnvvzsDg4ODOHz4MADQehmAnRxGR0cFo+dCgT4CKYtaKKjHN0ZILSCXqNiyNrIfg8F+g/r7WTr8v6/Moq4BAi65unPDaNp2Drq7u78vc3qs1CnfHHw3OSH3kRwLnVwn/v7+1O0D2K/rsbExDA8Po62tDVarlTY0CA4Oho+Pj6xrQnFZHGfs3bsXAQEBuOGGGxz8oK5AykK2WCxobGxEV1cXMjMzBWVfRF1ALGRXYDab0dDQgKNHjzq0bHIlfZqb7sw+TpvNhqoqu4b35s+maFlLOVI4Nj6+MdZhudVqxeHDh2kFNi4YhpEMrrGXhYeHY82aNdiwYQPy8vIcKtQlb92DtrY2TE1N0epuXIJ05TFdzmM92z0i5koRm0MOMbLldv7+/qh66DTUbNmIb35q90n/6wf+ePl8X7x/bRTvflxVV3CPj8QG+I5fDuEKuS6kIJQU4u3tjfDwcKSkpGD9+vVYv349oqKiMDs7i8bGRhQXF6OiogIdHR2i/QoVQj7OOP/883H//fcjKirKpWw9NsR8yMPDwygvL0dQUJCDfllMi3vVu/0uuUCGh4dRUVGBkJAQJ1+x3AJDYp2gSd3iuLg43hrHUn5jsvzrH6fCZrM5LC8vL4ePjw9+9N8J3mPhwpVHXVKPAThGCjd8NAqr1UqbiALA57clYmJiAo2Pn033IcdHLIdcSSIKG65I6/i2EdonWUYSQ7y8vKj2OT8/H+vXr/8+BRx449IQfPvttw6BUnL+YuckRNCfbnYuHSCmLWa/lnNzEroxSnXrIdBoNFTZotfrUVBQgMTERGi1WnR1daGkpAQHDx5Ec3MzBgcH6fWnEPL3uO222xAWFubw+MrG3r17ERQUhOzsbGRnZ+OJJ56Y0/7cqWdBQCq+sWE2m2kkOCsrC5GRkU5frhgpb/pbvaQflsjNenp6kJWVhYiICKd9iGXrAeLEOj09jVcu8AMA3L2foUTPF5STKqjE1iHbbDZ8dIP90TI1NRU/+FcbHXM8Uqhv2T0NvV5PXUcqlQrt7e0oKirCP8/zwduXL4fNOAObcQZTU1OYnp7GzMwMZmdnwVjMqH3kDMljYoPbzYK9LVeV4QqEnha41d7IOOJ7BeyB0sLCQodAKan/29raCsAxKUdKbUF+Y3zHJOZ35/u+pPzPbLjr2iOFpGJiYmjB+uzsbAQHB2NkZATl5eVobGyUnugEhcd9yLfccgvuuecebN68WXDM6aefjk8++cQj+4uIiEBJiXNwRA64GXGkljDx487lLitUAJ74imNjY0U7Sgtl64kRJyHpl8/3RXJyMmpyA6hMT+62fFb2t/dkwWg0oqKighJD/vNlTufHnqfu0TMFg2tCUfvkrXtQ/fDp6O3txWx3A/5wUTRefPFF9Pb2Yl1nJz6p6ETQmw+DMRuwOlCNTa9NYWrK/sf9rJY/73yeAc/Zn4qCgoIQGBiIpYGBqB+2Qu3ti2X//T1u3piC5cuXIzo6GoaOQXR0MIiOjqYuBSE3h9xAm1SQjBCy1DiVSuUQKAVA6yIDwMGDB6FSqXDzZ47xET6iZDcY4O5HTOom5irijuGDJ+tYcPsVWq3WRWshe5yQN27ciI6ODk9PK4iIiAj09fW5tS1JnZaT1ceFXEUCe8w7V4TCbDZLNkwFnF0WQq4JsnxmZgb/+oE/bv5sCrf+bwY1a4+lbhPpEfu4hTL22HOSZae9WIW/bFQhKyvr+zKlqyXJnYDvMTfx4f/CMtoD81AXtm0rQlNTE4JaWlB7uANBTx9zgdzwuvDn0+Re1VWYTCYMDg5icHDQad1fah1rOF/ytv3/0qVLER0djYGZJbjlvLX4xz/akZqaitTUVHpjEnIJCN2A+MCWa7mqPfbx8UH2U3ZXzi27px3WvXpRICwWC2pr7b+D8gc2OASMy36VL+pT55O6cc+TYK4+5LnClUD4iYgFUVkcOHAAWVlZiIqKwrPPPou0tDS354qKinLbZcEwDCwWCyoqKhAXF8fb6kgMrsrErvnP9yTw8Xd0eyHYFQY9DsuEyDRj+3788zwfu1W8JdDhmEiRdvb+pCxiMsZoNOK1i4Nw0yfjuHs/g5rTAx3GkNd855H2670A7G2WvvvuO5SVlaGhoQHqhgYcae8AYA+a/fZjwY/ANWi0UGl08PPSgGEYu3rBZAXAgGFs0KlVsFqtbmldR0ZGaJ2Sv/611HG3/stw1gY9HnxwD9asWYPs7GwwNitUag2v5I7PenRF58ydR4xIybp169Y5yO1Il2s/P7tbi8jvuOoJ7jxCkKusYMMTtZ+FsFjTpgFA5WIhDlmDOzo6cPHFF9M7MhsTExNQq9Xw9/fHp59+ip///OdU6uQObDYb9Ho9vvnmG5e2MxqNaGpqoll9c2klI0TKFf9XgJyni/H6JcFO6aVysPuWeERHRzuQJHt/pfetQ2NjIw2scdfz3TD4rGs+/++eO5Jx5MgRJCQkwMfHB+t+X+owJ9++zGYzampqcPmv34DxaCOMPQ2wTjhboqJQqaHxC4bGfxnOz0tBZGQkoqKiEBwcjO37BqHy8sG7P92Ia1+pgdrbB6WPXgB/f39kbN8v6QdtfPxszM7OYnx8HAW/2Q2bcRpPbgpBX18fhoeHYTab8cp3PbBMDME6OQTL5BCsk8OAzTXljMrLB+eccRoKCwtRUFCAvLw8ar26okSQAh/5ypHD8R0DWw9e9DO7T5bcVNnb8FnLUrI8PhQXF6OgoEDyHF0F4TN2EbETBLLuEMedkLlYtWoVDh06JJnmLHhADIOcnBzs379f1l2RYRj09fWhq6sL8fHx6Ovrw6pVq6i14C6kLGU2kfERpthrMj95zTAMMp/8BjvPXYLExESEhITwWs0En25ebe/aK2ARs/dX/sAGWmio/IEN0Ol09ObFvqmQeVpaWnD2/X9BpvUwDh486FC/QhgqaIMjoFsWg7sv24jk5GQkJydjxYoVOPOvNWj+9flOBDM1NYXc35WIyrmECJlLGnyWILsN1L+vDMPIyAjUajXu2DML6/QYrJNDeP7CKBw+fBh1dXVoaGhAVW09YJVW1eh0OqhD43HP9Rfh/PPPR0FBAbRardMx8LVSYp8D97zlkq/Q5wXYiXjDhg00TXz/3RkYHR3FD9/spmP2/iQNZ75U57StEDlLkTHDMDhw4MC8EbJKpZJ0CS4AZBHycZe99fX10bvYwYMHYbPZaJESd0Dap3PrBPBhdnYW1dXVtPD6smXLeJUW7kBuaUqudSmUmkxKPHLdA0ajETU1NQCA27+Y5a1Fwd7myztTeLuG8B3b3p+koaKiAl/emQLgWIt2kp4OAIzVgr9sVGHLli3Izc2FXq/H6Fc7sG/fPl4yVmm94b0iHffffz927dqFyFtfQOx972Gqrx1hV23Fu9pN2Lx5M9avX49Nf6uHSm1/jOUG+/iy2thqBz5wSVjM98luA5WZmYmEhARERkZi310Z0PovhXdkEv6vwh87pnLwzDPP4JtvvkHsL99FdXU13nzzTTzyyCO47LLLoPFf6jS32WyG8WgjnnvuOZx33nlYuXIl7rzzTswcLsHMzAwdRyxUoXPgnrsUhMiYuy27ZsfGv9Q4kHHxvTkOx8g9NleUFQTz6a5Y7D5kjxPy9ddfjw0bNqCpqQkxMTHYuXMnXnrpJbz00ksAgHfffRfp6enIysrCvffei7feemvOH6BU5xDSZLS6uhorVqyghdcBe4R2LunT7kDMd0teazQavHd1BF1Gxq999jtER0fzugz43p/9j0ZKpkJ98wi6u7uRlZVFo9XscWVlZRj65Dl0v3AjLrroIrzwwgu8riZNYBiuuuoq/O53v0NxcTEmRocQccNT+LfmTDxatxReYavRvP1CANKFgriE+/bl4k9Rch6dxSxG9nuGscsFIyMj6XZENnf6n6uRsX0/VGoNzGYz1q9fj/vuuw9vvfUWou/+F+rr67Fjxw7ceuutSE5Odjqm4eFhvPrqqxj8z28QHR2Nq6++Grt27YLNyE98UuQrpr3mk7pJzcvWEC9btuz7TD47SI9FNj67OU6WsoJgviu9LWZCnheXxfHGHXfcgRtuuAF5eXlO6wwGAxobG+Hn54f4+HinO3NPTw8YhkFMTIzb+7dYLGhtbcXs7CztuCAFOe6Kty9fjmvft5cufPOHS6FWq+l7rpUtldwht+A+e9ya//sPpmr3IKTrW0HljI+PDzZt2oQLLrgATzcEQhu4HDVbNkKlUjlE0cUuWDHrTWo9n4XGJno5bgy+8STpgmh9XQl0fXF7Eq3XoFKpkPCr9/DcaVrs2bMHn3zyiWBjXpWXD27ffCPuuOMOZGdn856DlLtCrgyNbEPcJELjpWRwXOy6MADBwcE01ZmPeMfHx9HV1SWYqzAX2Gw26HQ6WcXFjjNOTJfFfIBP+kaajNbW1mL16tVISkrifUyaq4VMxOiBgYHIzMx06OorBjmuC6vVSl0X1384gtTUVF5iZZMu1/IuvjfHYQx7foKDv8xz2OdHH32EK6+8Er1/uw3j37zmRMbR0dHwz7kQYVdtRUdHB9566y081xMHbeByui92WVQ23AkAAfYAqdT2cpI1+MhdaDyxtMQsST63ybk7m5H7uxIUFRUheeseaHyDcMYZZ+CPf/wjWltbEfGj5xCYfzVtGUXAmAzYsWMH8vPzcfrppyP68gdgMznXT+EjaO45iFnD7G0IGfONF7sRCsny9Ho9li5dirGxMVRUVKC4uBi1tbXo6emhae9ys/TcxWK2kE8aQma7LGZmZlBRUQGj0Qi9Xo/g4GDBbd2t+GaxWNDU1ISuri6HjD6NRkPJQwpi/lzAni7c29uLsl/lO60XSjzhjlOr1XjrsmW86wCg+uHT4ePjA4vFgqnarxAan46bbroJX3zxhYNETO3thzvvvBN79+61ZzLuews+8XlY9/tS3htC1m+/hclkgsVioS4TdoozG3Iy9lQqFT0PsUQMKb8yGSM0nswp5ouUuimQ1+xiVOv/cAhrtn2FQ4cOwTsqGfXv/gGlpaWoqanB9u3boV3q+IRWWlqK4c/+hPGdt+HRRx/FyMiIZNadO/I5bi0Svs9WrrICsF9PoaGhSEpKwvr165Gfn4/o6GiYTCaa9t7S0oLp6WmPltw8WXBSuCzee+89HDp0CA8++CC6urowMDCA5ORkWd2eZ2Zm0NbW5tLj0+joKA4fPoyYmBin1OqamhoqFXNFo+wKhFQYfO9L71uH5uZmJ9ld8b05KPhTBRibFQ/FHcWzzz5LU2/ZWLIyCy8++jNsqw9B3ePnOqzjWt02m42SL1Fq1D16JpVP1T5yBlQqFX3f+PjZNKAk5E8m6wwGAxoaGqDX611ycXDn4tvGVQmXmFZXzL0g9rjPMAyM3XWYrPwMM01FgNVRaqfy9kNg3mVoef95mqU3F6kb+UzZcjy5n4WYm0gKDMPg8OHDmJ2136wmJyfh5eWFkJAQhISEICgoaE7uBpvNBm9v7xOxHvKp47KIjIxEV1cXduzYAavVitzcXFlkDPDXsxACsYpJ9+qoqCgnK4qdYSdkvXLrSLDHvXvVsRrIpJg731xiPmL2+7znDsJqtToEY6oeOg3+/v54JKEXR//xE9x1110OZKzy9kPghmsQ9eMdGKj5Btdccw3UOm9RTbPVaoXFYqGqF1JPgZBv/WOb6GdFgmMpj39JL2ibzQaGYdD4+Nm8BMZuSCsl5SJj+JQVcgJZcpUbQnMILeOzzOsePRMA8P61UXj5pkx8/MJj2PO/3QjZdBvi4+PpeMY4jfFvX8eaNWvw/PPPU0WLnM+Cu45sQ26ehFjlkKtUkosUSFOJ8PBwZGZmorCwEBkZGfD398fAwABKS0tx4MABNDQ0oK+vjxK3K1BcFgsIs9mM9957D/v27UNmZibi4uJcujsK1YzggrQ8CggIcOpezQbXBcJHyly3Q8b2/bTc4lXv9qP64dN5t5ciYfY69vrNn07SYvQAkPbge7jssstw5513wjJ2LMAUHByMhx9+GEcONyBk42bogo+pPIhvnKsGIch+qghardbBOuF2niYZVGq12mFd7SNnUEK3Wq2wWq1o2HqW07lxwUfcUhaumLtECHwuAlfnEDoeYg3GxsZi+fLl2LBhA27/YhaB667AK6+8gkceeQQrVx4rcjQ0NIQHH3wQGRkZMLSXC87N51PmgkgJxdwTQtvOBVyVxZIlSxAREeFQfjUsLAxTU1OoqalBUVERqqur0dXVxdtlhWAxd5smWPSEfO211yIoKAjh4eHIz8+X3oADtuXFB6vViubmZlGrmA1SpN5VnP7nahz8ZR49JsDRBymWVSe0jm35V/xfASofLMRdIbU4uvNufP3113RdcHAwgk67ETU1NXjwwQdp6Ucu+PZd9qt86jPndiZOfeJr3tdcpP9mH9RqNT1vdpovYCcHYkGzayEDrhGImBKBD3t/ksbrn5Xr1uBCbBkpLMSuZZGXl4ctW7agpKQEyy78BZXgAXaJ4sA7j+Gee+5xaLLgSnBPzrFLKSvcIWupOhakgWpCQgLy8vKwYcMGrFy5kvayLCoqQnl5Odrb2zE6Ourgh17MadPASUDIb775Jh599FFYLBa375BCpDw6OoqysjL4+/uLWsVs8BWpl0oaISDpyUIKDD7LWMyHzE6mSH/gHZx33nl48MEHwZjtSTQqlQp33303amtrEVx4PU57sUqU+Lk1kQH7xaPT6ajFS4iX/K9/bJPTOu568p5tQZN15JE+5+liBAQE0M+2/rFNlKDl+GgJXFFi8F3Y3O3n4sJgL+OW3iRQqVTIf74M/hnnoKmpCX/6058cEql27NiBiPg0rLz1DwDsPnspvzr7uNn1kN2RubkDVwsLqdVqBAUFYdWqVcjOzkZhYSFSUlLg5eWFnp4elJSU4LvvvkNzc7NgEstiwaInZG9vb3rhuEvIXKuWWMVHjhyRZRWzIaTakEvKBBnb92PXhQH0tVzLmGu91mzZiJ+FNqL31ftQWnqsMI5ueSzCb/wdnnrqKSd/O1+tZLVaDYZhnKLiOp3OyUXBJVu+dexlQmQO2DWrxF1y6RtdSP/NPgD2m4PFYqF/fDWAhd7LsXj/+6NVgt+5kNUsRtJSy0j2mphPPGP7fjzXEwef6/8In6QNdJ1lvB/9bz6MH4x/Qn2wAH+na/b8X96Z4rBPV8nYXVfGXCu9kS4r0dHRSE9PR0FBAXJychAcHHwiBvNcwrwdvVSheoZhcO+99yIhIQGZmZkoL3f2h7mCwMBAjI/LS8rggq1FJlaxn58fsrKyXC5SItbpQ4qUues3fzopuF6IhHeeuwSf35YIwO4jfPPNN+1WsdFeklGn0+Ghhx5Cb2M5vKNTHLYXy+S7Zfe0vXTo924JEpjjuiK4fmO+dXzuC+662kfOQGtrK1paWpCVleVE7uTxnpAmezvATi7EggbkuTHYZMqVvUlZwlIkLTQPGc+1kLnbsrfT+AVjpOorvPzyyyxJJ4OXXnoJb7zxBm7+1F5sampqCtXV1Q7blj9wjMitVivOf7nF6Vikzpc71lXMR1doUhN5LkXCTgTMGyHfcsst2L17t+D6zz77DIcPH8bhw4fx97//HXfdddec9sfVIrsCnU4Hg8HgYBVHR0e75YsSspAZhkFHRwdePl/4ByOWvMEmTi4JA6A9127/YhaRkZGofvh0jBW9iR//+Mf0BqHxX4plN/wODz30ELy8vGR3EKl66DQAoIFBQohirgjuawI5hA0A5eXlUKlU0Ov1tA8ge0V3jLIAACAASURBVF4yN9tnSJQcxM2R8viXSHn8SwCgpTcJQYv5WQF+HbIrfmJusgr7NZdok7fuwel/rqZlUrnkzkf2KpUK119/PcrKyrBkZRZd/o9//ANDH/0OjMWMi17toJmjB3+Zh923xKOpqQkA8NZly3DeP4+lvns6cLcQWOx1LIB5JOSNGzfS7hJ8+PDDD7F582a7jyw/H2NjY4IppXIwF0ImLgp3rWI2+IJ6MzMzKC8vp6VC5bgvMrbvp7UshAoQ1WzZiAM/1wOwu26IOiNj+378+te/xvi3xyq8r1mzBjUHv4V3RAJvESK+mwHJuONa/CUlJVSWRCwuritCyBoWI2z2+82fTjooZsR80uztCEGzx5Ex5FzI+TRsPUvQir7k9U6c9mKVW35iwDnBgk9Wxybar/7fGqf5xNwXBGf9vQFhV2/DtddeS5fNNH2LgXcfd0hQCgoKwurVq3H9h/a6ztd9MEzX/fM8H5SWltL9id00xc5ZDmw226InzfnEgjlcenp6aOtvAIiJiUFPT4/IFuJwh5AJEY+PjyMyMtJtq5gNdlCPpG/X1dUhMTHRgWDkkPKV/7ang4sF+Uj/sAt3tUOlUqFmy0ZMN36LZ599ls6zZGU2ps5/HDExMbxBQL5jIVaxzWajwUFyod72uQERERGYmZlBfX099fEKuS/kBPrIZ/bOFaH0xsK2hNnz8c3NRyLcdem/2eeQmMKV2hFCanz8bHx0wwoU/SzbaU4CKZ+wq5lzZ/29gXe5HIJWabQ4sPJGBOReQpfNHqnCL37xC4c5hJQot31uoJb0+9dG4bvvvhM8zrla0vPVKYRA8SG7Cb4A3FzI0FVCHhsbo77iuLg4t/fLhUajgdVqhcFgQGVlJU3f5ktUcSXQxyXif/3AXiry1v/NOCyvq6vDzBd/otude+65+N+uP0Lt7cfr6uDuAzgWMAOOZdxxA3CFL1TSTsDr1h1LOnn9kmAcOHAAtbW16O7upgkpUoG+zCftDQbS09Oh1WpFrWzu9nLG8CWqcKV2gN0HbbFYYDQanX6jUnpkAinSEiJWUvZUDriWe/MT52Hg23cQdNqNdNmuXbswVfuVLE0yweVvH+UtA0pw9OhRt5I1COabkBc7FoyQY2Ji0NV1rMFid3c3oqKi3J5Pbvdpq9WKw4cPo729HRkZGYiOjoa3t7dHaiITmEwm1NTUYPXq1UhISBCt/SqHlNmEu/sWe/ZWfHy8k/vCNjuF039wOZX+aEMisXPnTofgILGIxeRyVqtVNH2VS5RE+QAAN348hvz8fKxcuRI2mw0tLS0OGYdsYT/DMA6dnfmsXKnHZ77j4b7mm4tY38TFQc6BYRg0NjbCYrHQFOXaR85A3aNnov6xTbQeB8CvrJByYRDwWb7coJ7Q9kKWeMrjXyKo4Dpcf/31dP3I53+BebjLyaftTuLHgZ/rYTAYUFtbKztZg4v5JuTF7g5ZMEK+9NJLsWvXLjAMg5KSEgQFBTkI311FRESEJCETq9jHxwfZ2dk0WORugSEujEYjqqur7Y/5ubmiRY3YkFvcHgAueMWe4nzai1VO66/WHIRlzO7mUHn5IPTyR2iSx7f3ZDnNxd6vzWajRYz0zxygagohK5XPZcC2QgMCAhAbG4usrCyHBJd1vy9FSUkJysvLUVRU5GBt8ZGoXJ8z+3jE3BjcY2WPP/jLPBw6dIhW7iOfAduKJk8OfAQNuCd/A+wyu01/q+dVZbjiCml+4lx8G3EFLVbEmGcRWvoSkh773OkYXr0oELWPnCFbv7106VLEx8dj7dq12LBhA1atWkUNHHayxtjYmGDRIIWQxTFvxYWuv/567N27F0NDQwgPD8e2bdso6f3kJz8BwzC45557sHv3bvj6+uLll1/G2rVrXTz8Y+jr68ONN96IDz74wGmd1WpFW1sbpqamkJKS4hS1J73g9Hq9W/tmGAb9/f3o7OxEQkIC2tvbkZ2d7XJXBKliREKtdADANjuF7r/eCsZkr3Gw7OL70P7GVjrnvrvSsXTpUkGJW9mv8qHRaGixGYDfJSC1jLucLwEEAHac443Q0FDanYLdd7D2kTMcgnly5HJC+3dlOWAnZdI5hL1v7ti6R88EwzAOTwfEDQI4duFwBUIFe5K3OhcNEvMvmwbaMfLG/VSHHHb1NvjE5TqMYffRE9Mdy7GgGYbBzMwMRkdHMTo6iomJCeh0OoSEhGDp0qW0aFB3dzesVqtDOrinYLPZsGTJkhOVlBeup95CwGKxYO3atU7NTsfHx9Hc3CwatGMYBocOHeItcC8Fk8mEpqYmaLVaJCYmQqvVoqqqCikpKW719ZJbIY5b8S324nswtn8XAEC7NAZRd/wFtY+c6TAnl4SrHjoNVquV+on5IETAUvpjPuuVFGe6+r0Bh3nqH9sEg8GA0dFRnLPDLst69aJAh+at7JubnEAe37HxESzRVr97VTiuetc5BsF3XmI3CKKDJkTdsPUs6hYRIloA+OC6aISFhWHZsmWCZCunIht72Q8Nn+NPf7LHE3zi1iLs6sfpeqGqbXxwN5A3OzuLsbExjI6OYmxsjH4OgYGBiI+P92jXEKKWWeyEvLhDkixoNBoHCRN5lCKlNWNiYgS/KKl6FkIYHBxEZWUlIiMjsWbNGup3nYsLRG6gjx3kS9+2B5q6/9J1Lz75CFQqNSXe//5oFd2GvY+s337rFJUmj/NimXVSadHssQTj4+M4dOgQQkJCqE6YvY2Pjw+ioqLodoSM37s6AqWlpTh06BBaWlowPHxMriUWyBNyS3CPibglUlNTnbbjbiuW+k3+q9VqSsbsokls6SBJWiEgqdNSmXpciPmlm7adg3dn00F4YLa9DJ9uXs2baDIfZAw4Fw1au3Yt7X8pVLx+LljsdSyAk4iQVSoVzZIbHx9HeXm5k69YantXAhN1dXXo7+9HTk6OU8dssWw9KVitVodeemIgpGzqb6VEpQlYjidbIpzSnv93awJ9bzabqX6YuCjkELAQ5CSIbPhjuew0dDYhXvnvPuTn5yMzMxOBgYEYGrK3sCIqEyHCZM/FJViSpEOSJMRcG1J1OMT82uxsQqLuIGCXHiUp6ezCQgTcBBK2VSuUdEL+60KiqJuCYRjs2rXLabwYPJ0sQuqerF692qF4vdFoRENDA4qKilBVVYXOzk5MTk66TNCLnYyBk4iQASA0NBRbtmxBQ0ODpFXMhdwqbcPDw6ioqMDy5cuRnp7OG6Bw10KemJhAeXk5fH19HUpwiiFj+348vO6Ya+Sy886ASqN1sKB/8K82WK1Wqi3WP3OA100hRcpC/mChceUPbMCrFwXioxvsevO1z37nsJ1U8gh7PUmNvfztowCA7Oxs7LsrHR9cF+2wzd6fpDkpZrjzpf16Ly7c1Y5bdk87fQbc8+G+FptXaHv2b5BdMAmwW9Gzs7OYmJigmZBCQUG+7D82+JQTvmuOPXG98P43dN1cpXvugh3UU6vVCAkJQVxcHHJzc1FYWEglqK2trSgqKkJZWRna2tqcqrrx4WQg5JPGh/zdd9/h0ksvxVVXXYXt27e7HFCrq6vDqlWr4Ofnx7veYrGgpaUFJpMJycnJov7h7u5uqNVq2TI+m82GI0eOYGRkBGvWrKH5+AzDUH2uGEb2/B2TZR8BAB555BE88MADggWJCLhWozsBPLFlX/84FW1tbUhJSUFISIjDOqFtxI5PaB3fmH/9wB82mw1BQUG44p1eus3IyAiampqQkJCAM/5aK7ovMZ84WS/15DAX1D5yBnV9ECUHSQOXqm7HXj/bXYf+1/8PAOAVkYiJjhoHlQVxDfFtPx84ePAgcnJyZCktGIahsQUSKNRoNLS7SEhICHUT2mw2aoGfoDh1fMglJSV46KGH8MMf/hBnnnmmy2QM2IuTCGmRSXH6oKAgZGRkSAbrXHFZTE9Po6KiAgCQk5PjUBxFpVI5aHj5ULNlI8zDx/Tcf6uzida8IIXwXZGQCbkixFKYN/2tHmvXrqVkzB7PBdetwB3HXc8Fl0Rv/mwKt/5vBmFhYQDsnZD37duH2tpaREdHU32xkI/ZFRcG33Fz/7O35TsXdjcXsp7dXYXtg65/bJNDrIRASGesDQyjYywTgw4kfrzJGLAbNnJbNAlVdQsKCsLw8DAOHTrk0F3kZChQf1JYyORR5o9//CO8vLywefNml+fo6OiAr68vvYiBY3K56elppKSkyK5xMTw8jLGxMYcWPFwwDIOenh709vaK9v87ePAg1q1bJ6q+iC/9Pb76yl7QPfSKR+GbuN5BhVH5YCEYhkHO08UAgFcu8KOP6/vvzkBwcDC1xOdiFQP2uggrVqzAuTub6TK5gT85pMi3XkyeBtgJmcivxsbGcNGrHQDsadrE0iLuFPb8rh6PXBUKWcdVdoh9zmIBTLVa7WA9s90RBoPB4aYY+8BHUKn47bDjUWCouLgYBQXymgDLgcViwdjYGMbHx5GcnOyWMXaccOpYyCR4EhkZib6+Prfm4PbWYwcGXS04JGUhz87OoqqqCgaDQTCtmgsx9cXBkWMWO0kMydi+3yFw5+3tTS/uW3ZPU3/yxr/UoLy8HK9cYHfVuGoVc0GyH4VUDnICZWLrxKx3NtjktvnTScTHxyMkJASrV68GYLc8Y2NjYbFY0NzcTJ9ESOlSV9wRUv5mPqL98s4UquwQm5fvHMk2dY+e6ZCsAji7MVJ/+Rp9rQkMXVAyng8LVqvVYvny5fR7XeyYd0LevXs3kpOTkZCQgKeeespp/SuvvILQ0FBkZ2cjOzsbO3bscHtfcrL1hEBqIttsNrS2tqK1tRVpaWkuBQYJxIJ6fX19qK6uxsqVK5GYmCh5R1er1bTnnxApa1m97yxjvZSI9c8cEHzEZnf5uGX3NHJychzmfOUCP3xwXTT23ZXusD0bbLL48Hp7Zhi7dx97f0Kv+YKFYuvIazEis9lsdLvKBwsd5iLL1Wo1AgMDsXLlSmRnZ9NsQqHAEXt7oc+Ce8x86wj4slK58/MRMdv6Jinf7OQUooMG7ORsGmin773C+Gu2HK/Sm0JdUTyFkyGoN6+EbLVa8dOf/hSfffYZ6uvr8eabb6K+vt5p3LXXXovKykpUVlbijjvucHt/cyXk6elplJeXQ6fTOflzXQEfIZvNZtTW1mJkZAR6vd7hMVIM3CasfKSsCz52cRt7GpDzdLFDAXkpWRoAp1TpW3ZP00d8tvU8MDBArWuCb+/JQmJioiz5GxtS1qHYei7I3NPT05Sg6h49k3bAFnM5sF+T1HS2EuLty5dj57lL8Mal/N+ZEHGS19wxbMtXzMLnm0dov2Qsm5wBwNR/rJu4V5izFXk86yBzm5sqcMa8EvLBgweRkJCAuLg4eHl54brrrsOHH344b/uLiopyqyayzWbDwMAARkdHsWbNGsTGxs7pbst1WRCpXHh4OFJTU2UHNfjmApxJ2Ts2HSqt/Ydu6m+FsbvO4cKUIki+IBpgd2ckJiZi3bp1VIZ35kt1lLz/fpY9os2uq+GKm4ILsUCjkBuFOze7iL7QdygUzGPPx9722veHkJ+fjxs+GqXL/nmeD167OIgm3XCPW8jK5VvGF+wT2p7PfST02dpMBkzXHVvmFZHosH7PHcm09ChfoNDTUOpYSGNeCVluzeP33nsPmZmZuOqqqxwqwLmKgIAAl5scEpWDWq2Gr6+voOzNFZAfhtVqRVNTE3p6epCdnY3Q0FCX59JoNLz+aHbFNo1PIPzSjl3ME4eO3fSkXAbc90K+XO5N5H+3JiAlJYVqjFOf+BqNjY3o7++nbgKhfXLJRoywuWTFR0gWiwXvXBGKty8/lqAjRGRsEnNFucFO7Kh/bBPy8/ORmprqEFt49aJAfHxjLC3k5Cqk/Ojcz4vPTcI+r8mKT2Ez2JUUmsAw+MQdq9Xy+W2JiIyMhM1mc+pPOF8EPV+EfDKoKwjmlZDl1Dy+5JJL0NHRgerqapxzzjm4+eab52W/fGM6OztRX1+PpKQkxMXFObgG5gqr1YqysjIEBAQgIyPD7Uc1rsuCwGazOWTcBaz9IV1naC6BabADgLPLQug1IE7KZBlxX5z/cguWL1+OhIQEOu6Kd3oxOTnpUKAo9YmvHYKlcmRjQo/yfEQ1OTmJQ4cOYenSpbR/oztyNvYyuWoJ/TMHcPY/Gun7H/13Av7+/g5PDB9cF02lhkIBP6l9Cu2fD2Qum2kWEwf/Q5cHbbgGKo2dDA/8XE9jIyQgzq0NPR8EPZ8ui5MhbRqYZ0KWU/N42bJlVNd75513oqyszO39Ed3i9LRwBhYAGAwGVFRUwGw2Izc3FwEBAR77MklQ0Gg0Ij093aWO1Xzgc1mQ+ggqlYr6SL2Wx2LJalLNi8HgB0/BZrQ/LbhLylxi+OC6aKxdu1bUWr30Dfv3XbNlI7UUs58qQuoTX1O9KBuu+JW5SH3ia9TX1yMjIwNRUVFOVqyQ1S1XnSHXvcLe16a/HYuRVD5YiODgYFoCle9cxVwS3ONhH7uQuwmwGxujX+2AbcbeBUQTGAr/jLPx2sVBqHywEEFBQbznREhtvgh6Pl0WJwMZA/OsQ7ZYLEhKSsKXX36J6Oho5OXl4Y033kBaWhod09vbSyPO77//Pp5++mmUlJS4shsHXH755di6dSuvBphhGBw9ehRHjx5FUlKS0w+ztLQUa9eudfvLnZqaQmNjI0JDQzE8PIyMjIw5/wC7urqg1Wrp4yUhYq1W6+RGSLj3X+jb9f+BMdtLLvqmnI7llz7g0ASUzxITs8r23ZWOlpYWqkIQ0srKIRI2/nmeD4KCgmibKr4bg5DlXvXQaaivr6f94VzZr9h4sXOQmgewZycSQnZ1W7INX/DPVUyUfYzRPX+j75decC/2PnMH4uLi5kRc7MYCBITE2X9CaG1thZ+fH/7/9s48rKkz/fvfJOw7QQUFkX2nWoSqXRToojO26sz46+LVulR5q291tPPrgm11tJ1eVqd9p61ardZR26kC1lqdqljAta2Fuo2ssgiILLJkISGBJCfn/YN5jgmE5ASSEOV8rsurJTlbwuE+z3M/9/39BgSw02ox57r4fL69Lxiy+uKt3hhy4sQJrFmzBhRF4eWXX8Y777yD9evXIzk5GXPmzMHatWtx7NgxODg4QCgUYseOHYiJYW9l05fly5dj3rx5eOSRR/Re7+7uRkVFBdzc3BAeHm6w3Ozq1auIj483+xdLvPNaW1sRExMDDw8PlJSUICwsbMi25E1NTaAoCoGBgdBoNODz+XB0dBywfCj4ufXo+OFj5mff9Ax4pcw1uO1AlK1PA0VReo0lTk5OA07xdffTfc1YikCr1UIqlUIsFjNuy4Qrb05jcrN9A9Ivf34Q5eXlCAsLg7+//4DXYOj6Bvqsxj6LqWOQz9K3usHYuQYTZNmirLuG1pz1AN1bvucel4qS43v1Gp4sBZsATV4DgPLyckZi1JLcA23TgL0EZFuzceNGhIaG4g9/+AOA3pulpaUFDQ0NiIyMNFpuZkrPwhBKpRIVFRXw8vJCaGjoXcWuGzcQEBAw4PSQLXfu3IFcLkdwcDAEAgEEAoHJWs41a9Zg9+7dzM8+qYvh9dCfDI5eBhop753phnHjxmHm3mrmPYD9qNjYAlXfbfu+//VsL6jVanh5eTGddLpiSJffmKqn4EdSL8ZGl+aOWI1dt+5nUyqVjAksm89ragT886pJeGTrNZPXZ4ju2+Vo+3YDtD29KTunsZGoLiyweAAcCFMBuqSkBCEhIawaoczhfgrI90Wnni7+/v5M6RvxtpNKpaxqf43pWfSFpmk0Nzcz3nnh4eF6gXIoEpwErVYLJycntLa24ubNm+jo6GClIrd582Y99xXJ2X0QF+wCre2/ONg3j0w8+5acUmDm3mpWC3/GcsrGyt8MVVsAvYtjDz30EMaNG4fu7m4mGJNuOt0AOFAu19A5+8KmysLY9nHvndG7loEw9v30RTcY6zZ5mEJeegZ3stYywXj06NEoPvuDzYIxYDwHLRaLIZVKmUVqS1dy3Otu04T7boR87NgxnD9/HhkZGbh9+zbCw8NZ35SG9CwMoVKpUFFRAScnJ0RERBisK7516xacnJwGnS+jKIrpbOLz+ZDJZIzqlVqthre3N4RCIXx9fQ2mWCQSCZ5//nk9BxW3qIfh9/T/gu9oXByJWCiZGg0bGwWbm1c2NYI9tyIBTk5OEIvF+P1XtXrvkREzm+qEgQKtoQVOcxnK/oPNF9O0FtIL30B6MZt5zdfXF//+97/7dV8OF83Nzbh16xYSExPh6uraLwjrBnLyszmQgYsd61gAIzVlkZeXh8zMTMycOROvv/66WdOYpqYmaLVaBAUFDbhNW1sbamtrERYW1k+Y3txjGUKr1TJlbo6OjgZvMoqi0NnZCZFIBLFYDI1GAx8fH2Z6TwJ0T08PMjIycPjwYWZfR79g+P1+NZzHRZt1XX0xli9mm1c2NmLUzctefevhfgp7AwWvvKVRjLCRuUHOnHTDUDD1MGMLpZBCdGo7FJW/MK+FhITgiy++wAMPPGDR6qHBoNVqUVVVhe7ubsTHxxscuBhLc7AN0FqtFs7OzvY+Sh55Afn48eN444034O3tjR9//NH0Dn1oa2uDXC43KFSi0WhQVVUFjUaD6Ohokwt/xo41EGQqZ2rhztB+ZIFMLBaDoihmBO3t7Y3169dj27Ztd3fg8eGVMg/ejy4A39G4aJKx4DHYvPJAxyc/983Lsj2WLnuedIGnZ6/Km67y3ECYCowD5aRNzRwGO2o3tj9N0+gqzof4zD+h7ZYxrz/++OP44osvmBSBTCaDi4sL86D29PS0WdAi6UKhUIiQkBDWD4bBBGguINspW7duxbx58/CHP/wBZ86YP+KQSqVoaWlBdLT+6FEsFqOqqgrBwcHw9/dndXNJJBK0trYiKirK5La6o2KBQDDkxQmKohhzSeK0kJ+fj08++USvk9HBdyz8frcaLuMTWB/bVFA2tg3B1Ot9j2PoPVP7l65LhUwmw9RP7ta1Z/9hFBOc3Nzc+tkqDXR83eOSipOBrsmcYxnaxtQx1R0N6Di1HT0NJXqvr1ixAh9++GG/UaiuwLtMJmOcoH19feHl5WWVaX5nZyfKysoQERFhdBbJBjaVHHZubkoYeQEZ6P3FTZo0CRcuXDD7F6RQKBhTVOCuHrJcLkdsbKxZEpxyuRy3bt1CXFyc0e1IbbG5o2JzIAG6uLgY69atw5UrV/Ted42aBp+HX4CTv2E1sMFgbl6ZoihUVlYydckDHXMoJW4AkLs4HGKxGAqFgqmtNgdjKRc2jR+DhVJ2orPoCDqLjgDau4vFgYGB+OyzzzBr1ixWx9F1giaLbCTd5ePjM+QArZsvHmrJpyH6BujKykp8/PHH+Prrr7kRsjnk5uZi9erVoCgKy5YtQ2Zmpt77PT09WLhwIS5fvgw/Pz9kZ2cjJCRkUOeaNGkSzp49a/bNpdFocP36dSQlJUEmk6GiogIBAQGDkuDs6elBRUUFJk40rGvQd1TMppxtqKhUKpSUlODUqVP49NNPIZPJ9N53jZwKn0cWWCwws80rsz3OQKmAweSKdbc/OFeoF6CL/pICd3d3kyNocxlMgNbI2tH52/eQX8sFrdZ5gPD4WLXyVaxbt25I+isqlYoJ0BKJBHw+Xy9AsxXC0mq1qK6uhlKpHDBfbGny8vKwfv167NmzR6+qyE6xn4BMURSioqKQl5eHoKAgpKSk4ODBg3qjx88//xzXr1/Hzp07kZWVhSNHjiA7O9vIUQcmPT0dX375pdliPjRNo6ioCP7+/ujo6EBMTMygb3aKonD16lWDN4pux91AC3eWpr29HVVVVYiIiMDo0aPR1NSEN998E0eOHOm3rWvkVHhPew7OYyMNHMl0QDQn8OguwpF9dRlqTtrY9egei6ZpowG4dF0qeDweq8oQS4yW1eJmdBYehrwkH6D0yycTEhKwa9euAR/2Q0GtVusFaADw9vZmArShtRPyoPf19TUrXzxYtFottm3bhpMnTyI7O9vinX9Wwn4C8sWLF7FhwwacOnUKALBp0yYAwNq1a5ltZs6ciQ0bNmDatGnQaDQICAhAW1vboH65CxYswKpVq5jUA1sUCgWKiooQHByMkJCQIY9Yf/vtN6SkpOi9RhbuBAIBHBwcrD4qJivdCoUCcXFx/aoVrl+/js2bN+P777/vt6+TfzjcE5+Ae1wqBK6eRs/DJuAYGhkf+tMY/M/hVoPbEcwN9paqYih8LRlT/nEJWfP88Pz3HazPZe75yDFoSg3lzcvoKjkNRdWvTLcdISwsDK+++ioyMjJsVuKl0Wj0FoyJeSwJ0CqVCqWlpQgPDx+UmqG5dHd3489//jOcnZ3x+eefm/S3tCPsJyB/++23yM3NZdxAvv76axQWFuqt/CckJCA3N5cpEwsPD0dhYeGgFgX+8pe/YPr06Xj88cdZba+rcUFRFKZMmWKRp7xuQGZTzmZp5HI5ysrKEBAQgPHjxxv9TMXFxUhf9L9Q3Pi5/5sCB7hFToNH4hP45v88jGX5PazOP1AeWSqVYtqnV1htO9SqB2P7mDpHX37734cYveU9T7ro5aANjaBNUbY+DbEbT0PVUoWuktPoKj/PyGXq4jwuBuv+74tYsmQJfHx8WB/fGuhW9Ny5cwdKpRKjR4/G6NG9/oTmrLOYS0tLC1566SU899xzWLlypb3njPtiPwH50KFDOHXqlF5ALioqwtatW5lt4uPjcerUKb2AXFRUNKhOoy1btsDb2xsLFiwwuS3J9bq6uiI8PBzFxcWIi4uziFAJESsiwdhWo2KapnH79m00NTUhLi6OcVlmQ0lJCdJfzoTixs+gNf27FgWeo+EWNQ2u4SlwGZ8AnoP121XZNHWY2s7U8Y6/FIKOjg4kJCTAxcVlyItxfQO07jlpmoa6owHKql8hLzkNjei2wWO4hDyI7etWIjIyckgSrpZG553dRgAAIABJREFUN18cGxsLhULBjKB7enrg5eXF5KFdXV0tMri5dOkSVq5ciY8++ghPPfWUBT6FzWH1JVg/8w52Mpxkm6CgIGaaJBQKB3W+sWPHoq6uzuR2ra2tqKurQ0REBHMuYr9kiZufz+eju7ubUWaz1cJdWVkZXFxckJycbPZIPCEhAa1FPyDm7X+jq+IC5NfzoGq+m+OlZG2QXT4G2eVj4Dm6wCVkElzDkvHDXxcw6mvE5fmZb24x+13LfERPJ3kgDI1YjdX16v7MtvPPULDVaDRISkrS61A09zoJOX8c3S8fHf32v9Fddw3dtVegrL0KStZmcF+B52i4x6fCIz4d2YsT4OXlhfDwcLsp6VKr1SguLu51Kn/gAUYC1sfHB6GhodBqtZDJZJBIJKisrIRSqYSHh4deuaE5n4WmaeTk5GD79u04fPgwIiMNr2vcL9hkhMxGhnP79u0oLi5mFvW+++475OTkDOZ0yMvLw9GjR5lcdV/UajUqK3uDTFRUlF7db1VVFUaNGsXa824gtFotSkpK4OzsjDFjxsDX19fqaYqOjg5UVlYyC3eWIO69M1C11UFenI+uktMGp9R98X74BXhMmomjLycYLWMDjC+4sW2wMLRv389gLsYCOpvUh1bVDdWdanTfKkb3zcvoaa7slxMm8Jxc4Rb1CDwS0uEcnIBf1yQzdby2yMuyRSaTmZ0vpmkacrmcGUErFAq4ubkxAdrDw2PAAE1RFDZu3IgbN27gX//615CFuoYZ+0lZAKZlOLu7u/HSSy/h6tWrEAqFyMrKQljY4MqviouL8d577+Gf//xnv/dEIhGqq6sxYcIERr5RF7Z6FsYgAkA0TaOzsxMSiQQSiQQCgQC+vr5MB52lRstkCimXyxEfH2+VhQ6y6NR9qxjKmt+grLkEjaTZ6D4CDz84j4uG07jo3v/6R4DvZLkc42AX0UrXpTIj2MtvTGUlEGQMWqOCqrUWqpYq9DRXQ9VSCXXH7QEDMADwnN3hOmEiXKOmwS1yGvhOLihdl4rGxkY0NTUhISHBKnW8g6WlpQV1dXVITEwcUpkdTdN6KQ65XA5XV1cmxUG6CaVSKZYuXYoHHngAH3zwgb3rVLDBvgKyLWlvb8ef/vQn/PDDD8xrFEWhpqYGSqUSMTExAwatwWpQAKbL2VQqFcRiMUQiEaRSKRwdHSEUCiEUCgfd1trV1YXS0lL4+/sP2ZyVDSTw0TQNjaixNzjfvITu+v+Y2BMAjw8H33FwFAbCURgIh//+11EYCL6bj8lrH0zJm6UoW58GiUSCyWtzoJE0QyNuxtMTaBQXF+M/xSWAASU9fXhwGhsJ19AkuIQmwXlcNHj8u/fHniddmFRZVFQUfHx87GLRiqZpVFdX97p5JyRYvL6Ypmmmm1AikSA7Oxvnzp2DSCTCiy++iLfffvteqqQwxsgNyFqtFg8++CB++qnXHbmzsxM3btzAuHHjTFoqtbW1QSaTmT06H0w5W3d3NxOgOzs74eLiwgRoY1M5oPdGbmxsxO3btxEfH2/Wwp0l6Bv41KJGSH/JQlfpGTiPT4CqpVq/kcEEnp6e6PEIgIPPWAjcfbH00VB8fYOGwN0HfHdfCNx8IHD3RvnGuws6bCow2IyiaUoNStGJwwujMeeTAmgVUmiVnaC6JNBIWhDn3oXLJTdYpWsYeHw4+gXBaWwUXEMehEvIJAjc+k+5y9anQaFQoKSkBKNHj4azszMkEgk6OzuZNmehUAgvLy+bB2iSL/b29h6y2whbTp8+jY0bN+Lpp59GU1MTLl26hG+++WZIphV2wsgNyKR9+ty5c6ivr4dYLEZsbCyrKeBAehYDoVvOZshWyRyUSiWj4CaTyZhcm1AohLu7O/MHoVKpUF5eDkdHR0RHRw/rdG7AIKeloG6/hZ6mG+hpugFV0w2oOxow1FuI7+IBvrM7eA7O4Dk4gufgdPefwBE8R2fGzJPWakBTGoDSgKbUoCkNaK0GD45zx5W6dmi7u0AppKBV5jmVG8LBdxycAiLhPDYSTgERcPIPB9/JdcDtyYOira0N1dXViIuL65cjJQ9ssVisF6B9fX0tmvIyxGDyxUNBq9Vi586dOHr0KHJychhbt/uIkRuQAeDRRx+Fu7s7tmzZYlaTh0KhQE1NDRITE01ua80mD5JrIwG6q6sL7u7ucHZ2RkdHB8LDww3mwIcLNqkCbY8CB+ePw7zNR6EWNUItboRG1Ai1qBG0SmmDqxwaPAcnOHgHwMF3LBx8x8LRd9x/A3EEBC4erI5BAjFN07h58yakUikSEhJYVfX09PQwAZroUOgGaEs9mO/cuYPa2toh54vZ0tPTg9deew00TeOLL76wai3zMDIyAzJpq9y4cSO2bt2KZ555xqz9dfUsjJ3D1joUFEWhoqICEokELi4uUKlU8PT0ZFIcw3kTNzc3o76+HrGxsf0aPthA0zSoLjE0okZoOltBdUmg7ZKAUkhAdUlAdYlBKSTQKjqNLpQNGh4ffFdPCFy9wXfzgsDV67//9YbAewwcfcfCwWccBJ5C8Hjm/577Vn2QVuOhlrT19PTotTkPNUDr5ovj4+NtYol0584dLFq0CPPmzcOaNWvsIm9uJUZmQM7Ly8OpU6fQ1taGRYsW9WtdNgVN0/jtt9/w0EMPGXzfFupsfSELd2PGjMGECRMYyUGZTAaRSASRSASVSgUvLy/GRcQWCyHkIUFRFOLi4vTSNdZYXKO1FLRKGbQqJWiNqvcfpQKtUYPW9GBNkgv+3y8dTEMLT+AInsABEDiAJ3AAj6///yf/kg4/Pz88uu2ayUA7mAVDQzZNUqnUaiVtZNGYjKAFAgFrJbfhyBdfu3YNy5cvx+bNm/G73/3O6ucbZkZmQCa88847iI+PN3uEDABFRUUGAzKxVbLVqJi0dDc0NCAuLs6oOaRWq9VzEVGr1QZdRCwFyTGOHz/e5EKpJYOzpd07+h7bEscdyC/v9u3baGxstJo0ZV8MKbmR+0E3QMvlcsYl3Rru1H2haRqHDx/Gp59+er8s2LHBvgOySCTCc889h7q6OoSEhCAnJ8dgM4ZAIGDyucHBwTh27Bir42/duhUajQbLli0z+9pIy7OuSaOt1dnUajXKysoGvXBHNAdIgNZqtfDx8WFG0INdfCTVHY2NjYiPj4eHB7vcKcFaZWnDjTHDVDKToGkasbGxw7YIq1ar9UbQ5H6Wy+VITEy0SeMFRVH429/+huLiYnzzzTdDbsC6h7DvgPzmm29CKBQiMzMTH374IcRiMTZv3txvOw8PD8jlcrOPf+jQIVy9elVPUY4t165dY/QsdM1GbaFDAfQ+rG7cuIGwsDCLLdxpNBo9FxEATAUHW2FytVqN8vJyODg4DLm6434IzMaCMEGpVKK4uBhjx44dlK62taBpGpWVlZBIJPD29kZnZ29Jn+4I2tI1xzKZDBkZGYiKijLobnKfY98BOTo6GmfPnsXYsWPR3NyM1NRU3Lhxo992gw3IP/30E/bu3YtPPvnE7H1LS0sRHBwMZ2dn8Hi8IZezsUWr1TIr7/Hx8VZdqCO6tyKRiJnOktGzoaYEqVSK8vJyhIaGWqW6414I0GwCsC5Eg9pQSdtwolarmUVF3XyxIS1k3Rz0UBb5amtrsWjRIqxatQoLFy60mweTDbHvgOzj48P80oHeJzMZueni4OCASZMmwcHBAZmZmZg3bx6r49fU1GDNmjU4cOCA2ddWWVnJTO9ttXCnUChQWlqKUaNG2UTkuy8DdRGS31N7ezsSEhLg6jpwba0lsYcAfWCOLyOOQx5WbNTLBlPSZivMyRfrzqokEglomh5UgD537hzeeust7Nq1C1OnTrXEx7gXGf6A/MQTT6Clpb+4zAcffIBFixaxCshNTU0YN24cbt68ifT0dBQUFCA8PNzkuRUKBVJTU1FQUMD6ekk5m0QiQX19PXg8ntnTenOhaZrxIYuNjbWbkVRPTw9aW1tRW1sLrVbLVHCQNu/hHOFYOlgbG/kScRySi1cqlYybtVAo7PeAIqNPT09Pu1JpA+7WFyckJJid+wcMB2giVu/r69svQNM0jd27d+PQoUPIyclBYGCgpT7KvcjwB2RjsE1Z6LJ48WI8/fTTmD9/vukLpWk8+OCDOH/+PKs/ChKMyUIHn8/XWwQh4kAkKFmilZXkZAUCAaKjo+0qp0aU4yIjIzFq1ChGb0AkEkEmk8HV1ZX5LnS7CK0NMUNVq9X9Su1sgW65oVgsRnd3N7y8vJiAVFNTY7NqBbbQNI2amhrIZDIkJCRYrL7YkLu5j48P6urqEBcXh7///e/o7u7G7t27rTKzOnToEDZs2IDy8nIUFRUN6Ktnys/TRth3QH7jjTfg5+fHLOqJRCJs2bJFbxuxWAw3Nzc4Ozujvb0d06ZNw9GjR006ORMmTpyI8+fPmwycbBfuSKcUmdaboz3RF7FYjIqKCoSGhtqVJ5huHjshIcFgPbOuYpdIJGK6CMm03lzNW7YQzQd7WiAj+r91dXUQiUSMNjAZQQ9315ktR+wkQG/ZsgXHjx+HWq3GvHnz8Pjjj7NONZpDeXk5+Hw+XnnlFXz00UcGAzIbP08bYT8C9YbIzMzEs88+iz179iA4OBiHDh0C0OsMsHPnTnz55ZcoLy/HK6+8Aj6fD61Wi8zMTLO+SLJ6PJDtjbm2Ss7OzggICGACKNGeqKurg1wuZ4ISmcoauvlJwJNIJJg0aZLNcrJsUCqVKC0thVAoRFJS0oB/vDweD+7u7nB3d0dQUBBomkZXVxdEIhGqqqqYab1u3nWokOm2PaV1gLvuLAKBANOnTwePx2NG0GVlZejp6dGb1tsyQJN8sbUWYvsiEAjQ2NiI8+fP47PPPkN6ejp+/fVXVFRUWOV8sbGxJrcpKipCREQEIxb2/PPPmzWoszXDFpD9/PwM5neTk5MZq6eHH34YxcXFgz6Hv78/7ty5YzAgEx2KoXTcubq6IjAwEIGBgXpBiTgl6HbOubi4QKlUoqSkBH5+fpg8ebJdjPAIra2tqKmpQUxMjNm1oTweDx4eHvDw8EBwcLDetL6iogLd3d1MUBIKhWZ1ERKTVqVSicmTJ9uknZctA5W0eXt7w9vbm3HQ6OzshFgsRmlpKdRqNZPiMPe7MIfW1lbcvHlz0Plic6FpGkePHsXf//53vRFoeno60tPTrX7+gWhsbMT48eOZn4OCglBYWDhs12MK+0laWgESkHWV2/o2eVgqB9k3KJGprEgkQmlpKRQKBSiKwoQJE+xmug30PpiqqqrQ3d2N5ORkiwQ8Ho8HLy8veHl5ISQkRC8olZSUQK1Ww9vbm3lYDVSFQB5gY8aMQVRUlN18ZwD7kjY+nw8fHx89iyPSUan7XZAR9FADtG6+2FYPMIqi8OGHH+LSpUsoKCgYtPWaIYwVBsydO9fk/oZSsvZ0H/Xlvg7IAQEBuHPnDvOzLXUo+Hw+vL294ebmBrlcDicnJwQEBEAqleLatWugadrqFRymIBoZAQEBiI6OttqNaigokcWgW7du6XURknIqIksZGxs77E7LupCSNolEgsmTJ5td0qb7XQC99yTpqLx9+zY0Gg2rh5UhdPPFkyZNskngkcvleOWVVzBhwgQcP37c4ous+fn5Q9qfjZ+nPTEiAnJfzWJb6FAAgEQiQUVFBSZMmMDouxJBGY1GA7FYjPb2dlRXVzP2Tn5+fjYRI29qasKtW7eGRdyeNKEIhUKEh4czi0EikQg3b95Ed3c3BAIBIiIibH5txiABz8PDw2iO3Rx09SWAuy3vYrEYDQ0N0Gg0rDRJbJ0vBoD6+nosXLgQK1aswJIlS+xy5JmSkoKqqirU1tYiMDAQWVlZg+pNsBX3dUAeN24cSktLzVq4swRarRa1tbUQi8WYOHGiwUUtBwcHjB49mgnQKpUKIpEITU1NqKiogJOTE4RCIfz8/Myu4DCGRqNhFlmSk5PtotROIBDAz88P7u7ukEqlCAwMhJeXF8RiMerr65mgRbwIh2M20dnZibKyMquXtOmWVgL6miS6swnd2l9b54sB4MKFC3j99dexY8cOPProozY5Z1+OHDmCVatWoa2tDbNnz8akSZNw6tQpNDU1YdmyZThx4gQcHBywbds2zJw5k/Hz1DVXtjfuW7U3oFeKc9WqVXjvvfcwffp0o2pploJUKvj6+iI0NHTQI11SwSESiSCXy+Hm5sb8oQ62rIwEleDgYLubtpGcrKFFRV1ZSYlEYnNrI2KVZSuVNmPo1v6KRCIolUrw+XzG2cPaOWOapvHPf/4TBw4cQE5Ojt6CGYdR7LsO2RZIpVJkZ2ejvLwcP/30E1xcXPDYY48hLS0NKSkpFm9pJc68MTExFs176lZwkD9Cc8TpaZpGQ0MDWlpaEB8fbxMXCLaQRajOzk6znDNIY0ZnZyczm/D19YWXl5fFZhMUReHGjRvQarXDqtJmCLVajdLSUuZBPdT2ZjaoVCq89dZb6OzsxJ49e4b94XSPwQVkXWiaRmtrKwoKClBQUIDffvsN48aNQ2pqKtLS0hAfHz/okRZJA9A0jZiYGJuMUshKvUgk0qtaIPobBPKH6+LigsjISLsKKj09PSgpKYGPj8+QRNEH6iL09fUddLrHXlXagN7F2OLiYoP5YkOqfroLpoNNUbW3t2Px4sV4/PHHsXbt2vvZ2cNacAHZGGRklp+fj4KCApSXlyM2NhZpaWlITU1lnDlMQVTQdBfubA2pWiCjRjJKcnJyQnNzM8LDw+2qlRe4KzEaFRUFPz8/ix2X2MqThxXpIiQpDjbpHpI+sbcKD+BufTHbxViyeKyr4GauxGZpaSkyMjKwYcMGq3TcjRC4gGwOFEXhP//5DxOgW1pakJKSgtTUVEyfPh1+fn56f8g0TaO2thYdHR2Ij4+3q+mbWq1GRUUFxGIxU2tNRs/Wdis2BfnexGKx1SVGyflIukcsFkOhUDDqbX3FgXRL2hITE+1KpU1XQS4xMXHQszBDsqu6I2jdGRRN0/jhhx+wadMm/Otf/0JCQoKlPs5IhAvIQ6GnpwcXL15EXl4ezp49C5VKhUcffRSpqakYO3YsDh48iCVLliAsLMyupm8kDUC80fh8fj9pTZJztbVym66553B9b6SLkHwfRBzIy8sLra2tjPGoPf1ONRoNSkpK4O7ujoiICIv+vvoKaPH5fNTX18PBwQGlpaUoLCxEVlYWRo0aZbFzEqztGmRncAHZkkilUpw7dw47duxAYWEhJk+ejIcffhipqalISkqyi5ZeMtU2lQbo7u5mpvQymcwiFRymIGJKRD3OXtBqtWhpaUFVVRWcnJzA4/EG3ZhhDUi+OCQkxCYiVCqVCidPnsS2bdtQU1OD6OhopKam4tVXX7V42svarkF2BheQLc3+/ftx8uRJ7NixA93d3cjPz0d+fj6uXLmCCRMmYMaMGUhLS0NMTIxNR1harRbV1dWQy+VmC6IT5TYSoBUKhVkVHGyOX19fj7a2NiQmJg67+llfSElbQkIC3N3d9TrnxGIxKIrqV/drK9ra2lBTU2PT5p2GhgYsXLgQS5cuRUZGBjo6OnD+/HmkpaVZ3P/O2q5BdgYXkC2NWq2Gg4NDvxEkEcDJy8vD6dOnUVlZicTERKSmpiI9Pd2kK/NQIE4jo0ePZr0QaQxdYSCRSASVSjXkVl4y1banNADbkra+mr9kwXSoVQvGIHl2ksu21UPgl19+wWuvvYZt27ZhxowZVj+ftV2D7AwuIA8XFEXhypUryMvLQ0FBATo6OvDQQw8hLS0N06dPh4+Pj0UCNJGktHTdsy66I0aRSAStVsuMFo25V5PqE3sTaweGVtLWt2pB11XGEl2EGo0GpaWlcHV1tdlDjKZpfPXVV9i3bx+ys7MREhJisWMPp2uQncEFZHtBqVTi559/Rn5+Ps6dOwetVss0qEydOtVsvWAyutNoNIiNjbXpNJrUuZJVeh6Pp2eOyuPx0NDQgDt37tjUg48tli5pI4ti5PsgFS2+vr5mV7R0dXWhpKQEEyZMsJlpgVqtxttvv407d+5g//79Nm0asrZrkJ3BBWR7hKZpiMVinD17Fvn5+bh48SJ8fX2RmpqK1NRUZmo2EHK5HKWlpYwO83A3LOhWcEgkEqhUKri6uiIyMtJiMwFLoFtuZ82SNl1XGdJFSEbQnp6eAwbo4cgXi0QiLF68GNOnT8e7775r85SSLVyD7AguIN8LEMcJUv987do1hIeHMwE6MjKScUypr6/HnTt3hkWhzRREJyMwMBACgWDYvfd00VVps3VJG6loIW3efW2/ADD5Yls6VJeVlSEjIwPvvvsu/vjHPw7L76WjowPPPvssbt26xbgGCYVCPdegX375Rc81aM2aNVi6dKnNr9UCcAH5XkSr1aK8vJwJ0MS26Pbt25gxYwbWrl1rFwptBPJAaW5uRkJCgl6DjKEKjoGaMqyFrVTa2KDbRSgWiyGTyaDRaODu7o7o6GibPbBOnjyJ999/H/v378fEiROtfj4OAFxAvj+4ePEiFi9ejMmTJ6OpqQmdnZ2YOnUqUlNT8dhjj1lUTMdcNBoNysrK4ODggOjoaJMLWn0rOIjfHAnQlh4d9i1psydIvtjf35+ZUZAHFklxDOTLOFi0Wi3+8Y9/4MyZM8jOzmakXzlsAheQ7we2b9+OWbNmMavKXV1d+Omnn5CXl4cLFy5AIBAwC4RTpkyxmkdbX2QyGUpLS4ek4UHsjDo6OvRqfsmi2GBnAmTRk6IoxMXF2ZWgEgDGlKBv6ommacjlcuaBRboILeFgrVAo8Oqrr0IoFOLTTz8d9oaXEQgXkO93aJpGR0cHTp8+jYKCAhQWFmLMmDFMg0piYqLFgxFN02hqarLKyFPXOUQsFg+qpMyeVdpomkZdXR1EIhGrhUVdX0axWAyVSsUY55ozo2hsbMTChQvx0ksvYcWKFXb1nYwguIA80iB/8CT/XFxcjJiYGCZAD0XiEugNmOXl5eDxeIiJibH6yLNvSZmjoyMTjAylajo6OlBZWWmXKm2kvpjIoA5mYbFvFyGxdyIzCkPlj4WFhVi9ejU++eSTYXV/5uAC8ohHq9WipKSE6SBsaGhAUlISU8ExevRo1gGalNsFBQUhMDDQyldumO7ubr2SMl3d49bWVrtUaQN60wXFxcUIDg62qERr3y5C0rQjlUoRHh6OEydO4Msvv0R2djZCQ0Mtdl4AyM3NxerVq0FRFJYtW4bMzEy993t6erBw4UJcvnwZfn5+Fm84uQfhAjKHPiqVCr/++isKCgpw5swZKBQKTJs2DWlpaXjkkUcGLKVrbm5GfX29XZXbkYqFtrY21NXVgaZpxoPQVhUcbCD54ri4OKtbiJGmnZ07dyInJwcSiQRLlizBrFmz8Nhjj1nsQUVRFKKiopCXl4egoCCkpKTg4MGDerXBn3/+Oa5fv46dO3ciKysLR44cQXZ2tkXOf4/CKiDbj7jAMHPo0CHGNeTSpUsDbpebm4vo6GhERETgww8/tOEVDh0nJydMnz4dGzduxPnz53HmzBnMmjULFy5cwNNPP42nnnoKH3zwAX7++WeoVCrI5XIcPXoU7e3tSE5OtptgDAA8Hg8URaG5uRmxsbGYMWMGQkNDGfeWixcvoqysDC0tLVCpVDa/PtKIUl9fj6SkJJv4ORJH9cuXL+OFF15ARUUFHn30URw7dgxKpdJi5ykqKkJERATCwsLg5OSE559/HkePHtXb5ujRo1i0aBEAYP78+SgoKICZg78Rif0UtA4zCQkJ+O677/DKK68MuA1FUXj11Vf1RgZz5sy5F7uGAACenp6YPXs2Zs+erWdxlZ2djVWrVkEul+OJJ57A8uXL7UoYCOjVN2hoaEBiYiKzsOjp6QlPT09MmDCBqeAQiUS4ffu2XgWHpb3m+kLKAZ2dnfHggw/a7Lu7ceMGli5diszMTPzP//wPeDwe5s6di7lz51r0PI2NjXrmpkFBQSgsLBxwGwcHB3h7e6Ojo8OupFftES4g/5fY2FiT2+iODAAwI4N7NSDrwuPx4O/vjwULFsDLywtXrlzB5s2b0dzcjI8//njQFleWRqvVoqKiAhRFITk5ecCFReKEQfz6dCs4amtrAUDPRcVSC5QkXzx+/HibOnv/+OOP+Otf/4q9e/ciKSnJqucyNNLtey+w2YajP1xANgM2I4P7gaioKJw5cwbe3t4AgOXLl+tZXL322msmLa6sASlpCwgIwPjx4806n0AggJ+fHyPcTyo4WltbUVlZyVRwEOfqwYxqbZkvJmi1WmzduhW5ubn48ccf+5meWoOgoCA0NDQwP9++fbvfw4dsExQUBI1GA6lUCqFQaPVru9cZUQHZmBQgm2ndSHnqR0VF9XtNIBAgKSkJSUlJePPNN/Usrnbs2KFncfXwww9bvDPO0iVtjo6OGDNmDNNO3dPTA5FIhMbGRpSXl8PFxQV+fn6snKt164uTkpJsVuWhVCrx5z//Ga6urvjxxx9t1hSUkpKCqqoq1NbWIjAwEFlZWThw4IDeNnPmzMH+/fsxbdo0fPvtt0hPT78v/1YszYgKyPn5+UPan83IYKTg7OzMlM8Bdy2u8vLy8N5778HT0xMzZswYssWVrkpbUlKS1YKOs7Mzxo4dy5SlEQ2O2tpaxrla1+aKQFEUSktLbZ4vbm5uxsKFC/H8889j5cqVNg12Dg4O2LZtG2bOnAmKovDyyy8jPj4e69evR3JyMubMmYOlS5fipZdeQkREBIRCIbKysmx2ffcyXNlbH1JTU/HRRx8hOTm533sajQZRUVEoKChAYGAgUlJScODAAcTHxw/DldovNE2jpaVlyBZX9uI4QlqaxWIxOjo6mJZmDw8PNDc3Izg42KYP5kuXLmHlypX4+OOP8eSTT9rsvBxDgqtDNocjR45g1apVaGtrg4+PDyZNmoRTp06hqakJy5Ytw4kTJwAAJ06cwJo1a5iRwTslojboAAAIXklEQVTvvDPMV27/DGRxNWPGDKSnpxvUdSZaGaGhoTbJi5qDVqtFQ0MD6urq4OLiomfrZE3fPZqmkZ2djc8//xwHDx5EZGSkVc7DYRW4gGzvjDAbdAZDFldTpkxhFghzcnIQHByM6dOn251KGzFt7ejoYLoCKYrSs7kCwGhw+Pj4WKSCQ6PRYOPGjaiqqsLXX3/NLLhy3DNwAdneGWE26ANCLK5yc3Nx4MABCIVCPPHEE3jyyScHZXFlLUi+2MnJCVFRUQOmUNRqtZ7NlUAg0NPgMDf1IpVKsXTpUkyaNAnvv/++3anXcbCCC8j2zgizQTeKVqvFU089hWeeeQYvvvgizp07NyiLK2tBSu6CgoLMzhf3tXVydnbWcw0xtiBXVVWFJUuW4PXXX8cLL7zAVSrcu3AB2d4ZYTboJmlvb+/XyUXTNBobG5n889WrVw1aXFkTUnIXFxdnkVQBcQ0RiUSQy+V6FRy6ovQFBQV49913sWfPHoOLzEPBlDjQvn378MYbbzBCUitXrsSyZcsseg0jDC4g2wOcDbplMWRxNXHiRKSmpiI9PR3+/v4WG0WSfHF7ezsSExOtUnJH0zS6urqYAF1VVYVDhw7Bz88PFRUV+P777y2qEAewEwfat28fLl26hG3btln03CMYVjfliKpDHg6M1T77+/ujubmZSVkM5PlGpshhYWFITU1lRokjET6fj/j4eMTHx2P16tVQq9W4dOkS8vPz8fLLL1vM4oqiKJSVlcHR0RFJSUlWG4XzeDx4eHjAw8MDwcHBCA8Px8mTJ1FdXQ1PT088/fTTWL58OTIyMix2zvtZAuBex74UY0YYpJsJAPbv32+wW1AsFqOnpwdA75T+559/5v5wdHB0dMS0adOwbt06nD17FhcuXMDcuXNRVFSEefPm4cknn2TU7cj3aAqlUonLly9DKBSyrpm2BC0tLZg/fz6mTJmCCxcuID8/H7/88gtmz55t0fMYkgBobGzst93hw4fxwAMPYP78+XoNURzWgxshDyOZmZl49tlnsWfPHsYGHYCeDXp5ebmeDXpmZiYXkI3g7u6OmTNnYubMmXoWV0eOHEFmZqZJiytL54vZcvXqVaxYsQJbtmzBrFmzmNednZ0t3nTCRgLgmWeewQsvvABnZ2fs3LkTixYtwunTpy16HRz94XLIIwTO4cG4xdWMGTOQk5ODBx54AE899ZTNdCFomsbhw4fx6aef4sCBA4iOjrb6OS9evIgNGzbg1KlTAIBNmzYBANauXWtwe4qiIBQKIZVKrX5t9zFcDpmjFzY6znv27IGvry+qq6uRlZWFt956675zeODxeAgNDUVGRgYyMjIYi6vjx49j9uzZEAqFaGhogEKhQGpqKsaMGWPVMjOKovD++++jtLQUBQUFNvMBZCMORNY2AODYsWOs5Gk5hg6XQx4BcA4PhuHz+UhMTMS5c+ewceNGXLp0CRkZGairq8OiRYswY8YMvPXWW8jNzYVMJrPouTs7O7FgwQIAvd+9LU1ZdcWBYmNj8eyzzzLiQKQL9LPPPkN8fDwmTpyIzz77DPv27bPZ9Y1kuJTFCODbb79Fbm4uvvzySwDA119/jcLCQr2SpoSEBOTm5iIoKAgAEB4ejsLCwhHh8CCXy+Hh4dHvdZlMhvPnzyM/Px8//fQTnJ2dGQW7lJSUQcts1tbWYuHChVizZg1efPFFrtljZMClLDh64RwejGMoGAPGLa5ef/11jB07FqmpqUhLS2P8GE1x9uxZZGZmYvfu3ZgyZYqlPwrHPQ4XkEcAnMPD0NG1uFqwYAFomkZNTQ3y8/Px0Ucf6VlcpaWlITg4WO+BRtM0du3ahcOHD+PkyZNMBxwHhy5cymIEwEbHefv27SguLmZs27/77jvk5OQM41XfW+haXBUUFOhZXE2bNg2bNm1CT08Pdu3aZTdiSRw2hd10k6Zpc/5x3KMcP36cjoyMpMPCwui//e1vNE3T9Lp16+ijR4/SNE3TSqWSnj9/Ph0eHk6npKTQNTU1Qz7nyZMn6aioKDo8PJzetGlTv/f37t1Ljxo1ip44cSI9ceJEevfu3UM+p73Q3d1Nnzlzhn7nnXfokJAQeuXKlTRFUcN9WRzDB6sYy42QOawCp5dwF5qmLZqPf/nll/HDDz9gzJgxKCkpMXi+1atX48SJE3Bzc8O+ffus7kTNYRJWNwBX9sZhFdiU2o0ULL04unjxYuTm5g74/smTJ1FVVYWqqirs2rULK1assOj5OawHF5A5rAKnl2A9pk+fbnTB9ejRo1i4cCF4PB6mTp0KiUSC5uZmG14hx2DhAjKHVTCUCjOkl1BXV4fr16/jiSeeYBpTOIYG24chh/3BBWQOq8Cm1M7Pz4/RjMjIyMDly5dteo33K2wehhz2CReQOayCrl6CSqVCVlYW5syZo7eN7jSa00uwHGwehhz2CReQOawCp5cwfMyZMwdfffUVaJrGr7/+Cm9vb4u7jnBYB67sjeO+434vC3vhhRdw9uxZtLe3w9/fHxs3boRarQYALF++HDRNY+XKlcjNzYWbmxv27t1rcU8+DrOxiqceB4fdw+PxpgOQA/iKpukEA+//HsAqAL8HMAXApzRNc8ISHMMOl7LguO+gafo8AJGRTeaiN1jTNE3/CsCHx+Nxc3qOYYcLyBwjkUAAukXPt//7GgfHsMIFZI6RiKF8Hpe74xh2uIDMMRK5DWC8zs9BAJqG6Vo4OBi4gMwxEjkGYCGvl6kApDRNc73FHMMOJ1DPcd/B4/EOAkgFMIrH490G8FcAjgBA0/ROACfQW2FRDUABYMnwXCkHhz5c2RsHBweHncClLDg4ODjsBC4gc3BwcNgJXEDm4ODgsBP+P0IEy1e3MTQcAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x = numpy.linspace(-1.2,1.2,30)\n", "y = numpy.linspace(-1.2,1.2,30)\n", "t = numpy.linspace(0, 2*numpy.pi)\n", "xx = numpy.repeat(x, len(y)).reshape(x.shape[0],y.shape[0])\n", "yy = numpy.tile(y, (len(x),)).reshape(x.shape[0],y.shape[0])\n", "\n", "\n", "fig = pyplot.figure()\n", "ax = fig.gca(projection='3d')\n", "ax.plot_wireframe(xx,yy,circle(xx,yy),cmap=pyplot.cm.coolwarm_r,\n", " linewidth=0.01, antialiased=False)\n", "ax.plot(numpy.cos(t),numpy.sin(t), 0, linewidth=3, c=\"black\")\n", "ax.view_init(elev=40., azim=30)\n" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "678E9C7603434BC68617DD76EE238D5E" }, "source": [ "Now we can pick a point $x$, say $1/2$ and seek the matching $y=f(x)$ on the circle, starting from $\\frac{1}{2}$. We know that actually $y = f(x) = \\sqrt{1-x^2} = \\frac{\\sqrt{3}}{2}$. Of course, we would run into trouble for $x$ close to $1$ or $-1$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "cell_id": "1D879F5FE3D14C618F9AA2ED54FC2C1E" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.8660251031723489 0.8660254037844386\n" ] } ], "source": [ "x = torch.tensor([0.5], dtype=torch.double, requires_grad=True)\n", "y0 = torch.tensor([0.5], dtype=torch.double)\n", "y= Implicit.apply(x, y0, circle)\n", "print (y.item(), (1-0.5**2)**0.5)" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "B75F8DDBE73D4A0BB15F41914462984C" }, "source": [ "Works! Let us compute the derivative. We can do that by hand, using the implicit function theorem, we have $\\frac{d}{dx} F(x,y) = 2x$ and $\\frac{d}{dy} F(x,y) = 2y$, so $\\frac{df}{dx}(x) = - x/y$. That is about the technical complication I can handle.\n", "\n", "If we didn't like using the implicit function theorem, we would have to do this by $\\frac{d}{dx} f(x) = \\frac{d}{dx} f(x) = \\frac{1}{2} \\frac{-2x}{\\sqrt{1-x^2}}$ and plugging back in $x$ and $y$ we see that $\\frac{d}{dx} f(x) = -\\frac{x}{y}$.\n", "\n", "But of course, we can also let the autograd do its thing:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "cell_id": "E9A5A59E194646BE9DAE7B9C44721367" }, "outputs": [ { "data": { "text/plain": [ "(torch.Size([1]), torch.Size([1]))" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y.backward()\n", "x.grad.shape, y.shape" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "4AACABC131B4409880261154E7AA8A17" }, "source": [ "Awesome. In fact, we can also use `pytorch`'s automated checker (and indeed this is why I used `DoubleTensor`s, as the gradient checker can be too strict to use single precision floats). It seems that gradcheck does not like the 1d $x$." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "cell_id": "74596AC8D005455F851D5B8DE48DC51E" }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "torch.autograd.gradcheck(lambda x: Implicit.apply(x.unsqueeze(0),y0,circle), x)" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "54903BD23B58404A889C1A2C31C311B1" }, "source": [ "So this toy example works great. But is the implicit function really useful?\n", "One thing is that the limitation to scalar $x$ and $y$ is quite severe.\n", "\n", "Quite likely, this is only the case for very specific use cases. In general, more common alternatives are\n", "- Manually solving for the implicit function.\n", "- Adding a free parameter $y$ and introduce a penalty $\\lambda F(x,y)^2$ into the loss function.\n", "\n", "However, when there is just this last transformation you want to add and a free $y$ isn't a good option for you, this might be nice to have.\n", "\n", "I hope you enjoyed this little feature, your feedback is very welcome. I read and appreciate every email you send to .\n", "\n", "Thomas" ] }, { "cell_type": "markdown", "metadata": { "cell_id": "0F9A83A4955547DF85BDC676EBA4CBD0" }, "source": [ "P.S.: And here is a **thing to try** based on a recent [PyTorch issue](https://github.com/pytorch/pytorch/issues/7698) (thank you for the idea!):\n", "\n", "HIPS autograd has a [fixed point with parameter function](https://github.com/HIPS/autograd/blob/master/autograd/misc/fixed_points.py).\n", "Write a PyTorch equivalent. Note:\n", "- Careful with the naming: $x$ in the HIPS autograd fixpoint function is $y$ in the implicit function above and $a$ is $x$.\n", "- You replace the iteration in the forward by the simpler fixed point iteration, but you can write the fixed point equation $f(x,y)=y$ as an implicit function problem to compute the derivative.\n", "- In particular you get to store the result and I don't think you need to iterate again." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.6.5rc1" } }, "nbformat": 4, "nbformat_minor": 2 }