{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "In the previous post, I showed how to solve inverse problems for coefficients of elliptic PDE using [firedrake-adjoint](http://www.dolfin-adjoint.org/en/latest/).\n", "The exact parameter field that I used in that demonstration was smooth in space and, to guarantee a smooth solution, I showed how to add regularization to the objective functional.\n", "Many geophysical inverse problems aim to estimate fields that instead have sharp discontinuities or interfaces.\n", "For example, the porosity of soil and hard bedrock are very different and there is no continuous transition between the two.\n", "For these media, the regularization functional\n", "\n", "$$R(q) = \\frac{1}{2}\\int_\\Omega|\\nabla q|^2 dx$$\n", "\n", "that we used in that demonstration would yield an infinite value.\n", "The inferred field with this penalty would have a more diffuse interface than the real one.\n", "\n", "Rather than use the integrated square gradient, we can instead use the **total variation** functional:\n", "\n", "$$R(q) = \\int_\\Omega|\\nabla q|dx.$$\n", "\n", "We can get some insight into why the total variation is a good regularizer for these types of problems by using the very wonderful [coarea formula](https://en.wikipedia.org/wiki/Coarea_formula).\n", "The coarea formula states that, for reasonable functions $p$ and $q$, we can express certain integrals involving the gradient of $q$ in terms of its contours or level sets.\n", "Let $ds$ be the element of surface area, let $z$ be an arbitrary real value, and let $\\Gamma_z$ be the $z$-contour surface of the function $q$.\n", "Then\n", "\n", "$$\\int_\\Omega p|\\nabla q|dx = \\int_{-\\infty}^\\infty\\int_{\\Gamma_z}p\\, ds\\, dz.$$\n", "\n", "The right-hand side of the last equation can make sense even when $q$ is discontinuous, provided we're a little careful in the definition of the $z$-contour of $q$:\n", "\n", "$$\\Gamma_z = \\partial\\{x \\in \\Omega: q(x) \\le z\\}.$$\n", "\n", "For example, suppose that $\\Gamma$ is some nice closed surface inside $\\Omega$, and we take $q$ to be equal to $\\alpha$ in the interior of $\\Gamma$ and $0$ outside.\n", "Then the coarea formula tells us that\n", "\n", "$$\\int_\\Omega|\\nabla q|dx = a\\cdot|\\Gamma|.$$\n", "\n", "This partly explains why the total variation functional is an effective regularizer.\n", "While it doesn't forbid a jump discontinuity as such, it instead penalizes (1) the magnitude of the jump and (2) the area of the surface over which it occurs.\n", "Gabriel Peyré has a nice visualization of the coarea formula on [Twitter](https://twitter.com/gabrielpeyre/status/985768327246237697)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Calculating total variation\n", "\n", "A new difficulty that we'll encounter here is that the total variation functional doesn't have a well-defined functional derivative like the mean square gradient does.\n", "It is a convex functional, so the minimum is well-defined, but we might be at a loss for an algorithm to actually approximate it.\n", "\n", "We've already encountered the mathematical concepts that we'll need to remedy this issue in a previous post on the obstacle problem.\n", "The obstacle problem is the prototypical example of an optimization problem with inequality constraints.\n", "To solve it, we reformulated the obstacle problem as an unconstrained convex optimization problem where the objective could take the value $+\\infty$.\n", "We then smoothed away the infinite values by instead working with the Moreau envelope of that functional.\n", "\n", "Many of the same tricks work for total variation-type problems because the Moreau envelope of the $L^1$-norm has a simple analytical expression in terms of the *Huber function*:\n", "\n", "$$H_\\gamma(z) = \\begin{cases}\\frac{1}{2\\gamma}|z|^2 & |z| < \\gamma \\\\ |z| - \\frac{\\gamma}{2} & |z| \\ge \\gamma \\end{cases}$$\n", "\n", "The Huber function looks like the $L^1$ norm for large values of the argument, but like the square norm for small values." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "zs = np.linspace(-5., 5., 41)\n", "γ = 2.\n", "H_γ = [z**2 / (2 * γ) if abs(z) < γ else abs(z) - γ / 2 for z in zs]" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXhU9dn/8fedPRBIIAlhC4RAAglhkwgCIgiCbLJUbcXW7VF5aLVqta3KpoK41Mdq1bbudam1rQsYEBUQxAUEArIkIYQQliQsCYQkBJKQ5fv7I2N/FIMzJDM5s9yv68rVmcxhzmew183JyTmfrxhjUEop5fn8rA6glFLKOXSgK6WUl9CBrpRSXkIHulJKeQkd6Eop5SUCrNpxVFSUiYuLs2r3SinlkbZs2XLMGBPd2GuWDfS4uDjS09Ot2r1SSnkkETlwvtf0lItSSnkJHehKKeUldKArpZSX0IGulFJeQge6Ukp5CYcHuoj4i8h3IrK8kdeCReRfIpIrIhtFJM6ZIZVSStl3IUfodwO7zvParcAJY0wv4BngyeYGU0opdWEcGugi0hWYDLx6nk2mAW/aHr8PjBURaX68HyosreSRZZnU1NW74u2VUsqlnl2dw/b8Upe8t6NH6M8CvwfON0W7APkAxphaoAyIPHcjEZklIukikl5cXNyEuJBZWMbfvtnP61/va9KfV0opq3y95xjPrt7D59lFLnl/uwNdRKYARcaYLc3dmTHmZWNMqjEmNTq60TtX7RrftyNXJMXw7Oo9FJw43dxISinVIqpq6pj/UQbdI1vxq9E9XbIPR47QRwBTRWQ/8E9gjIj8/ZxtCoFYABEJAMKB407M+V8enpoMwCPLsly1C6WUcqqX1uWx79gpFk1LISTQ3yX7sDvQjTEPGmO6GmPigOuANcaYX5yzWRpwk+3xNbZtXLa2Xdd2rbj7igRWZR1lZeYRV+1GKaWcYv+xU/z5i1ym9O/EZYlNOzvhiCZfhy4iC0Vkqu3pa0CkiOQC9wIPOCPcj7n10h4kxoTxyLIsTp+pdfXulFKqSYwxzP8og2B/P+ZPSXbpvi5ooBtjvjDGTLE9XmCMSbM9rjLGXGuM6WWMGWKMyXNF2LMF+vuxeEY/Cksr+dPne1y9O6WUapLlOw7z1Z5j3Dc+kZi2IS7dl0ffKXpxXHt+mtqV177aR/aRcqvjKKXUfymvqmHh8iz6dQnnhmFxLt+fRw90gAcmJtEmJIB5SzKor3fZaXullLpgf1yZw7GKahbPSMHfzyW35vwXjx/o7VsH8eDEJNIPnOC9LflWx1FKKQB2FpTx1ob93HBJd/p3jWiRfXr8QAe4ZnBXLo5rx+OfZFNy6ozVcZRSPq6u3jB36U4iw4L57ZW9W2y/XjHQ/fyER6f3o6KqlsdXnK9uRimlWsY7Gw+wo6CMeZOTaBsS2GL79YqBDtC7YxtuHdmD97YUsGlfidVxlFI+qqi8iqc+3c2lvaKYOqBzi+7bawY6wN1jE+gSEcq8pTs5U6vlXUqplvfox7uorqtn0fQUXNRReF5eNdBbBQXwyNS+5Byt4NWvXX4pvFJK/Zev9hSTtv0QvxzVkx5RrVt8/1410AGuSI5hXHIMz32+h/wSLe9SSrWMqpo65i/NIC6yFb90UfmWPV430AEentoXQXg4LRMXVsoopdR/vLhuL/uPn2bRdNeVb9njlQO9S0QovxmXwOfZRazMOmp1HKWUl9t37BR/WbuXqwZ0ZmSC68q37PHKgQ5wy4ge9I5pwyNpmZyq1vIupZRrGGOYvzSD4AA/5k9OsjSL1w70hvKuFA6VVWl5l1LKZZbtOMzXucf43YTedHBx+ZY9XjvQAVLj2nPdxbG89vU+dh3W8i6llHOVVdawaHkW/buG8/Oh3a2O490DHeD+CX0IDw1k3lIt71JKOdfTK3dzvKKaxdP7tUj5lj1eP9DbtQ7iwYl92HLgBP9O1/IupZRz7Cgo5e1vD3DjsDj6dQ23Og7g2CLRISKySUS2i0imiDzSyDY3i0ixiGyzfd3mmrhNc83grgzp0Z4nPs3meEW11XGUUh6urt4wZ8lOosKCuXd8otVx/sORI/RqYIwxZgAwEJggIpc0st2/jDEDbV+vOjVlM4kIi6enNJR3fZJtdRyllId7e8N+MgrLWTAluUXLt+xxZJFoY4ypsD0NtH153MnohJg23H5ZPO9vKWBj3nGr4yilPFRReRVPr8xhZEIUU/p3sjrOf3HoHLqI+IvINqAIWGWM2djIZleLyA4ReV9EYs/zPrNEJF1E0ouLi5sRu2nuGpNA13ahzFuaoeVdSqkmWbg8q6F8a1rLl2/Z49BAN8bUGWMGAl2BISKScs4my4A4Y0x/YBXw5nne52VjTKoxJjU6uuXvpgoN8mfhtL7sKdLyLqXUhfsyp5jlOw5zx+hexFlQvmXPBV3lYowpBdYCE875/nFjzPe/bXwVGOyceM43pk8MV/bV8i6l1IWpqqlj/kcZxEe1ZvboeKvjNMqRq1yiRSTC9jgUGAdkn7PN2SeSpgJuvWzQQ1f1xU+Eh7S8SynloL98sZcDtvKt4ABryrfsceQIvROwVkR2AJtpOIe+XEQWishU2zZ32S5p3A7cBdzsmrjO0TkilHvHJbImu4jPMrW8Syn14/KKK3jxi71MG9iZEb2irI5zXmLVEWpqaqpJT0+3ZN8AtXX1THn+a8oqa1h97yhaBwdYlkUp5b6MMfzitY3sKCjj8/tG0aGNtX0tIrLFGJPa2Gtef6fo+QT4+7F4Rj8Ol1Xx7Oocq+MopdxU2vZDfJN7nN9P6GP5MLfHZwc6wODu7Zg5pBuvf7Nfy7uUUj/wffnWgK7hXD+km9Vx7PLpgQ5w/4TeRIQGMnfJTi3vUkr9l//7bDclp86weIZ7lG/Z4/MDPaJVEHMmJbH1YCn/0vIupZTNtvxS/r7xADcNjyOli3uUb9nj8wMd4CcXdWFoj/Y88Uk2x7S8SymfV1tXz9wlO+nQJph7x7lP+ZY9OtCxlXfNSOH0mVoeW+HWl9ArpVrA298eIPNQOQum9KWNG5Vv2aMD3aZXhzbMuiyeD7cWsmGvlncp5auOlDWUb41KjGZSv45Wx7kgOtDPcuflCcS2D2X+R1repZSvWrQ8i5q6ehZO6+t25Vv26EA/S2iQPwunppBbVMErX2l5l1K+5ovdRXy88zB3Xt6L7pHuV75ljw70c1zepwMTUzry3Od7OHhcy7uU8hVVNXUs+CiT+OjWzBrlnuVb9uhAb8SCq5IJ8BMeSsvQ8i6lfMSf1+ZysOQ0j7px+ZY9OtAb0Sk8lN+MS2Tt7mI+yzxidRyllIvlFlXw4rq9zBjUheE93bd8yx4d6Odx8/A4kjq15eG0LCqqa62Oo5RyEWMM85dmEBroz5xJSVbHaRYd6OcR4O/HYzNSOHqyimdWaXmXUt5q6bZCNuQd5/6JfYhuE2x1nGbRgf4jBnVrKO96Y/1+Mg+VWR1HKeVkZadrWPzxLgbGRjDzYvcv37JHB7od91/Zx1belaHlXUp5mT98lm0r30rBzwPKt+xxZAm6EBHZJCLbbasSPdLINsEi8i8RyRWRjSIS54qwVghvFci8KUlsyy/l3c0HrY6jlHKS7w6e4B+bDnLLiB707ewZ5Vv2OHKEXg2MMcYMAAYCE0TkknO2uRU4YYzpBTwDPOncmNaaPrALw+IjeVLLu5TyCg3lWxnEtAnhNx5UvmWP3YFuGlTYngbavs499zANeNP2+H1grHjaPbM/QkRYND2Fypo6Fn+s5V1Kebo3Nxwg63A5D12VTJgXLT/p0Dl0EfEXkW1AEQ2LRG88Z5MuQD6AMaYWKAMiG3mfWSKSLiLpxcXFzUvewnp1CGP2qJ4s+a6Q9XuPWR1HKdVEh8sq+ePK3YzuHc2EFM8q37LHoYFujKkzxgwEugJDRCSlKTszxrxsjEk1xqRGR0c35S0sdcflvejWvhXzlmZQXVtndRylVBMsWp5Fbb1h4dQUjyvfsueCrnIxxpQCa4EJ57xUCMQCiEgAEA54XQdtSKA/C6f1Ja/4FK98qeVdSnmatdlFrNh5hLvGJtAtspXVcZzOkatcokUkwvY4FBgHZJ+zWRpwk+3xNcAa46UlKKN7d2Byv048vyZXy7uU8iCVZ+pYkJZBrw5h3D7SM8u37HHkCL0TsFZEdgCbaTiHvlxEForIVNs2rwGRIpIL3As84Jq47mH+lIbyrvkfaXmXUp7iz2tzyS+pZNG0FIICvPMWHLu/3jXG7AAGNfL9BWc9rgKudW4099UxPIT7xvdm4fIsPsk4wqR+nayOpJT6EblFJ3npy7385KIuDOv5g+s1vIZ3/jPVAm4c1p2+ndvyyLJMTlbVWB1HKXUexhjmLc2gVVCAx5dv2aMDvYkC/P1YPKMfRSereWbVHqvjKKXOY8l3hXybV8L9E/oQFebZ5Vv26EBvhoGxEfx8aDfeWL+PjEIt71LK3ZSePsPij3cxqFsE110ca3Ucl9OB3ky/u7IP7VsHMXdpBnVa3qWUW3ny092UVtaweHo/ryjfskcHejOFhwYyb3Iy2/NLeXeTlncp5S62HDjBu5sOcsvwOJI7t7U6TovQge4E0wZ2ZnjPSJ78NJvik1repZTVauvqmbc0g07hIdzjReVb9uhAd4Lvy7uqa+pZ/HGW1XGU8nlvrN/PLi8s37JHB7qT9IwOY/aoeJZuO8T6XC3vUsoqh0or+eOqHMb06cCVfb2rfMseHehO9KvLe9E9Usu7lLLSwmVZ1BvDI1P7el35lj060J2oobwrhbxjp3hpnZZ3KdXS1mQf5dPMI/x6TAKx7b2vfMseHehONioxmsn9O/HC2lz2HztldRylfEblmToWfJTp1eVb9uhAd4EFU5IJ8vfT8i6lWtDza/ZQcKKSxdO9t3zLHt/81C4W0zaE+8Yn8tWeY3y887DVcZTyejlHT/Lyl3lcfVFXhsZ7b/mWPTrQXeSGS7qT0qUtC5dlaXmXUi70fflW6+AA5kzqY3UcS+lAd5EAfz8WT+9HcUU1T6/MsTqOUl7rg62FbNpXwoMT+xDp5eVb9jiyYlGsiKwVkSwRyRSRuxvZZrSIlInINtvXgsbey9cMiI3gF0O789aG/ews0PIupZztxKkzPLZiFxd1i+Cnqd5fvmWPI0fotcB9xphk4BLgDhFJbmS7r4wxA21fC52a0oP99sretG8dzNylO7W8Sykne/LTbMoqa1g8wzfKt+yxO9CNMYeNMVttj08Cu4Aurg7mLcJDA5k/JYkdBWW8s/GA1XGU8hpbDpTwz8353HppD5I6+Ub5lj0XdA5dROJoWI5uYyMvDxOR7SLyiYj0Pc+fnyUi6SKSXlxcfMFhPdXUAZ25tFcUT326m6KTVVbHUcrj1dTVM3dJQ/nW3WMTrI7jNhwe6CISBnwA3GOMKT/n5a1Ad2PMAOB5YGlj72GMedkYk2qMSY2Ojm5qZo8jIiyc1pfq2noeXb7L6jhKeby/fbOP7CMneXhqX1r7UPmWPQ4NdBEJpGGYv2OM+fDc140x5caYCtvjFUCgiEQ5NamHi48O45eje5K2/RBf7fGdn06UcrbC0kqeWbWHK5I6MD45xuo4bsWRq1wEeA3YZYz543m26WjbDhEZYnvf484M6g1+OboncZGtWPBRJlU1Wt6lVFM8kpYJwMM+WL5ljyNH6COAG4AxZ12WOElEZovIbNs21wAZIrIdeA64zug97z8QEujPoukp7Dt2ihfX7bU6jlIeZ1XWUVZmHeWusQl0bed75Vv22D35ZIz5GvjRfwaNMS8ALzgrlDcbmRDNVQM685e1e5k2sAs9olpbHUkpj3D6TC0Pp2WSGBPGbSN7WB3HLemdohaYPzmJ4AA/Fmh5l1IOe+7zXApLK3l0ej8C/XV0NUb/VizQoW0Iv72yN1/tOcayHVrepZQ9u4+c5NWv8rh2cFeG9GhvdRy3pQPdIr+4pDv9uoSzaHkW5VrepdR51dcb5i3dSVhIAA9OSrI6jlvTgW4Rfz/hsRn9OF5RzdOf7bY6jlJu6/2tBWzef4I5E5No3zrI6jhuTQe6hfp1DefGYXG89e0BdhSUWh1HKbdz4tQZHl+xi4vj2nHN4K5Wx3F7OtAtdu/4RKLCgpm7JEPLu5Q6xxOfZHOyqpZHp2v5liN0oFusbUggC6Yks7OwjLc37Lc6jlJuY/P+Ev6Vns+tI3vQu2Mbq+N4BB3obmBK/06MTIji6ZU5FJVreZdSNXX1zFuSQZeIUC3fugA60N2AiLBoWgrVdfUs+ljLu5R6/et97D7aUL7VKkjLtxylA91NxEW15o7RvVi2/RBf5mh5l/JdBSdO8+zqPYxLjmGclm9dEB3obmT26Hjio1qz4KMMLe9SPuvhtKyG/53a6LIK6kfoQHcjwQEN5V37j5/mr19oeZfyPSszj7B611HuuSKBLhGhVsfxODrQ3cyIXlFMG9iZv36xl7ziCqvjKNViTlU3lG/1jmnD/1yq5VtNoQPdDc2dnERwoB/ztbxL+ZDnPt/DobIqFs9I0fKtJtK/NTfUoU0Iv7+yN9/kHidt+yGr4yjlctlHynnt6338LDWW1Dgt32oqHehu6vqh3RnQNZxFy3dRVqnlXcp71dcb5i3JoE1IAA9M7GN1HI/myBJ0sSKyVkSyRCRTRO5uZBsRkedEJFdEdojIRa6J6zv8/YTFM/pRcqqa/9PyLuXF3tuST/qBE8yZlEQ7Ld9qFkeO0GuB+4wxycAlwB0iknzONhOBBNvXLOCvTk3po1K6NJR3/X3jAbbna3mX8j4lp87w+CfZDIlrr+VbTmB3oBtjDhtjttoenwR2AV3O2Wwa8JZp8C0QISKdnJ7WB903PpHosGDmLNlJbV291XGUcqrHV+yioqqWR2ek6ILPTnBB59BFJA4YBGw856UuQP5Zzwv44dBHRGaJSLqIpBcX692QjmgTEshDV/Ul81A5b397wOo4SjnNpn0lvLelgNtGxpMYo+VbzuDwQBeRMOAD4B5jTHlTdmaMedkYk2qMSY2Ojm7KW/ikSf06clliNE+vzOGolncpL3Cmtp55S3dq+ZaTOTTQRSSQhmH+jjHmw0Y2KQRiz3re1fY95QQN5V19OVNXz8LlWVbHUarZXvt6HzlHK1g4rS+hQf5Wx/EajlzlIsBrwC5jzB/Ps1kacKPtapdLgDJjjK5+7ETdI1vz68t78fGOw6zT8i7lwfJLTvOnz3O4sm8MY5O0fMuZHDlCHwHcAIwRkW22r0kiMltEZtu2WQHkAbnAK8CvXBPXt80aFU98tJZ3Kc9ljOHhtEz8RHjoKi3fcja7RcPGmK+BH/31s2m4P/0OZ4VSjQsO8OfRaSlc/+pG/rI2l3vH97Y6klIXZGXWUT7PLmLupCQ6a/mW0+mdoh5meK8oZgzqwl/X7WWvlncpD/J9+Vafjm24eUSc1XG8kg50DzRnUhKhgf7MX6rlXcpz/OnzPRzW8i2X0r9VDxTdJpjfT+jD+r3H+Wiblncp97frcEP51swhsQzuruVbrqID3UNdP6QbA2MjePTjLMpOa3mXcl/19Ya5S3YSHhrI/RO0fMuVdKB7KD8/4dHpKZScOsNTK7OtjqPUef07PZ+tB0uZMymJiFZavuVKOtA9WEqXcG4e3oN3Nh7ku4MnrI6j1A8cr6jm8U+yGdqjPVdf9IM2EOVkOtA93L3jE4lpE8LcJRla3qXczmMrsjl9ppZHp2v5VkvQge7hwoIDWHBVMlmHy3lrg5Z3Kffxbd5xPthawO0j40nQ8q0WoQPdC0xM6cjo3tE8vXI3R8q0vEtZr6F8K4Ou7UL59Rgt32opOtC9gIiwcGoKtfWGhcszrY6jFK98lUduUQWLpqVo+VYL0oHuJbpFtuLXY3qxYucR1u4usjqO8mH5Jad57vM9TOjbkcv7dLA6jk/Rge5Fbr8snp5a3qUsZIxhwUcZBPgJD009d6VK5Wo60L1IcIA/j07vR35JJS+sybU6jvJBn2UeYe3uYn4zLpFO4Vq+1dJ0oHuZYT0j+cmgLrz05V5yi05aHUf5kIrqWh5OyyKpU1tuHh5ndRyfpAPdC82Z3FDeNU/Lu1QLenZVDkdPNpRvBWj5liX0b90LRYUF88DEJL7NK2HJd7oSoHK9zENl/G39fmYO6cZF3dpZHcdnObIE3esiUiQiGed5fbSIlJ21mtEC58dUF+q6i2MZ1C2CxR/vovT0GavjKC/WUL6VQURoIPdfqeVbVnLkCP0NYIKdbb4yxgy0fS1sfizVXH5+wuLp/SitrOHJT3dbHUd5sXc3H2RbfilzJycR3irQ6jg+ze5AN8Z8CZS0QBblZMmd23LL8Dje3XSQLQe0vEs537GKap78JJtL4tszY5CWb1nNWefQh4nIdhH5RETOu/KriMwSkXQRSS8u1pXrW8I94xLp2DaEuUt2anmXcrrHPt5FZU2dlm+5CWcM9K1Ad2PMAOB5YOn5NjTGvGyMSTXGpEZHRzth18qesOAAHp6aTPaRk7yxfr/VcZQXWb/3GB9+V8j/XtaTXh20fMsdNHugG2PKjTEVtscrgEARiWp2MuU0V/btyJg+HfjjqhwOlVZaHUd5geraOuYtzaBb+1bcOaaX1XGUTbMHuoh0FNvPWiIyxPaex5v7vsp5RIRHpval3hgWLsuyOo7yAq98mUde8SkemdaXkEAt33IXjly2+C6wAegtIgUicquIzBaR2bZNrgEyRGQ78BxwndG7WdxObPtW/HpMAp9mHmFN9lGr4ygPdvD4aZ5fk8ukfh25vLeWb7mTAHsbGGNm2nn9BeAFpyVSLnP7yHiWfFfIgo8yGRYfpbWm6oIZY5hvK99aMOW81z8oi+idoj4kKMCPxdNTKDhRyfNr9lgdR3mgTzKOsC6nmPvG96ZjeIjVcdQ5dKD7mKHxkVx9UVde/jKPPUe1vEs57mRVDY8syyS5U1tuHNbd6jiqETrQfdCcSX1oHRzAXC3vUhfgmVV7KDpZreVbbkz/q/igyLBgHpzYh037Svhgq5Z3KfsyCst4Y/0+fj60G4O0fMtt6UD3UT9NjeWibhE8tmIXJ05peZc6v7p6w9ylGbRvHcTvtHzLrelA91F+fsLiGf0oq6zhyU+zrY6j3Ni7mw6yPb+UeZOTCQ/V8i13pgPdhyV1asutl/bgn5vz2XJA+9fUDxWfrObJT7MZ3jOSaQM7Wx1H2aED3cfdPTaBzuEhzF2SQY2Wd6lzLP44i+qaehZp+ZZH0IHu41oHB/DQ1L4N5V3f7Lc6jnIj63OPsXTbIWaPiqdndJjVcZQDdKArxifHcEVSB55ZnUOhlncp/n/5VvfIVvzqci3f8hQ60BUiwkNXNZR3PZKWaXUc5QZeWpdH3rFTLJyWouVbHkQHugIayrvuHpvIyqyjrM7S8i5ftv/YKV5Ym8vk/p0YlajrFngSHejqP269tAcJHcJ4KC2T02dqrY6jLGCMYUFaJkH+fiyYkmx1HHWBdKCr/wgK8OPR6SkUllby/Jpcq+MoC6zYeYQvc4q5b3wiMW21fMvT6EBX/2VofCTXDu7KK1/mkaPlXT7l+/KtlC5tueESLd/yRI4scPG6iBSJSMZ5XhcReU5EckVkh4hc5PyYqiU9OCmJsJAA5i3R8i5f8vTKHIorqlk8vZ+Wb3koR/6rvQFM+JHXJwIJtq9ZwF+bH0tZqX3roIbyrv0lvL+lwOo4qgVkFJbx1ob9/GJodwbERlgdRzWR3YFujPkS+LH7wqcBb5kG3wIRItLJWQGVNa4dHMvg7u20vMsH1NUb5i7ZSfvWwfz2yt5Wx1HN4Iyfq7oA+Wc9L7B97wdEZJaIpItIenFxsRN2rVylobwrhfKqWhZ9rAtLe7O3Nuxne0EZ8yYnafmWh2vRE2XGmJeNManGmNToaL2+1d316diWX47qyYdbC1m7u8jqOMoFDh4/zR8+3c3o3tFavuUFnDHQC4HYs553tX1PeYFfj+1Frw5hzPlwJyeraqyOo5zIGMP9H+zA3094bEY/Ld/yAs4Y6GnAjbarXS4Byowxh53wvsoNBAf489Q1/TlaXsXjn2hvujd5d1M+G/KOM2dSEp0jQq2Oo5zAkcsW3wU2AL1FpEBEbhWR2SIy27bJCiAPyAVeAX7lsrTKEoO6teO2kfH8Y+NB1uceszqOcoLC0koeW7GLEb0imTkk1v4fUB4hwN4GxpiZdl43wB1OS6Tc0r3jElmVdZT7P9zBp3dfRutgu//XUW7KGMOcD3dSV2944if99VSLF9G7B5RDQgL9efLq/uSXVPLUZ7utjqOa4YOthazLKeb+Cb2Jbd/K6jjKiXSgK4cN6dGem4Z1580N+9m8X5es80RF5VUsXJbJxXHtuHFYnNVxlJPpQFcX5PcT+tAlIpT7399BVU2d1XHUBTDGMHdpBtW19Tx5dX/8/PRUi7fRga4uSOvgAJ68uj95x07xzOocq+OoC7B8x2FWZR3lvvGJxOuScl5JB7q6YCN6RTFzSCyvfJnH9vxSq+MoBxyvqOahtEwGxEZw66XxVsdRLqIDXTXJg5OS6NAmhN+9v11Pvbi57xetOFlVw1PX9MdfT7V4LR3oqknahgTyxNX9yDlawcLl2vXizt7dlM/HOw5zzxWJJMa0sTqOciEd6KrJRvfuwOxRPfnHxoMs/U7bHtxRRmEZDy/L5LLEaH45qqfVcZSL6UBXzfLb8YkMiWvPnCU7yS3SFY7cSXlVDXf8YyvtWwXxzE8H6FUtPkAHumqWAH8/nps5iNBAf371zlZdXNpNGGO4//0dFJyo5IXrBxEZFmx1JNUCdKCrZusYHsKz1w1kT1EF85bqsnXu4I31+/kk4wj3T+hNalx7q+OoFqIDXTnFyIRo7hqTwIdbC3kvXZets9J3B0/w2IpdXJHUgdtH6iWKvkQHunKau8YmMKJXJPM/ymDX4XKr4/ik0tNnuPMf3xHTNoSnrx2oxVs+Rge6chp/P+HZnw0iPDSQX72zVRfEaGH19Yb7/r2dopNV/Pn6iwhvpcvJ+Rod6MqpotsE8/zMQRwsOc2DH+7U8+kt6OWv8vg8u4h5k5MZEBthdRxlAR3oyumGxkfy2/G9Wb7jMG+u3y7wL3sAAAspSURBVG91HJ/wbd5xnvpsN5P7d+LGYd2tjqMs4tBAF5EJIrJbRHJF5IFGXr9ZRIpFZJvt6zbnR1We5H8vi+eKpBgWLs9iZeYRq+N4tZyjJ5n1VjrdI1vxxE90bVBf5sgSdP7An4GJQDIwU0SSG9n0X8aYgbavV52cU3kYPz/huZkD6dc1gl+/+x1bDmh/uiscLqvkptc3ERzoz5u3DKFNiJ4392WOHKEPAXKNMXnGmDPAP4Fpro2lvEGroABevymVzhGh/M8b6XonqZOVna7hptc3cbKqljduuVhXH1IODfQuQP5Zzwts3zvX1SKyQ0TeF5FGV50VkVkiki4i6cXFxU2IqzxNZFgwb94yhEB/P256fTNHyqqsjuQVqmrquP3tdPYdO8XLNwymb+dwqyMpN+CsX4ouA+KMMf2BVcCbjW1kjHnZGJNqjEmNjo520q6Vu+sW2Yo3brmY0tNnuPlvmyir1MsZm6Ou3vCbf21j074Snv7pQIb3irI6knITjgz0QuDsI+6utu/9hzHmuDGm2vb0VWCwc+Ipb5HSJZwXbxhMblEFs95K1w71JjLG8MiyTD7JOMK8yUlMHdDZ6kjKjTgy0DcDCSLSQ0SCgOuAtLM3EJFOZz2dCuxyXkTlLUYmRPN/1w5g474S7vv3durr9Rr1C/WXL/by1oYDzLosntv0tn51jgB7GxhjakXkTuAzwB943RiTKSILgXRjTBpwl4hMBWqBEuBmF2ZWHmz6oC4UnazisRXZRLcJ5qGrkvUyOwe9l57PU5/tZvrAzjwwoY/VcZQbsjvQAYwxK4AV53xvwVmPHwQedG405a1uHxnP0fJqXvt6H4H+woMTk7Sr244PthTwwIc7ubRXFH+4RrvNVeMcGuhKOZOIMHdSEmdq63nlq32UnKrhiav7EeivNy435pUv81i8YhfDe0by4g2DCQrQvyfVOB3oyhJ+fsLCaX2JDAvi2dV7KD19hheuv4jQIH+ro7kNYwxPfJrNS+vymNSvI8/8bCDBAfr3o85P/6lXlhER7rkikUXT+rJmdxE3vLaRstN6SSNAbV09v39/By+ty+PnQ7vx/MyLdJgru3SgK8vdMCyO52cOYntBKT99aYPP33xUVVPH7L9v5b0tBdw9NoFHp6fgr+fMlQN0oCu3MKV/Z/528xAKTpzm6r+uJ6+4wupIliirrOGG1zbyefZRFk7ry2/GJepVQMphOtCV27g0IYp3Z11CZU0d17y4gR0FpVZHalFF5VX87KUNbMsv5bnrBnHjsDirIykPowNduZX+XSN4f/YwQgP9ufbFDbz97QGfWCTjy5xiJj//NQdLTvP6zRdzld4BqppAB7pyO/HRYSy5YzhD4yOZvzSD295M51hFtf0/6IGqaupYuCyLG1/fRLtWgXzwy+GMTNCeI9U0OtCVW+rQJoQ3br6Yh65K5qvcY0x49iu+2F1kdSyn2n3kJNP//A2vf7OPm4fHkXbnpSR1amt1LOXBdKArt+XnJ9wyogdpd44gKiyIm/+2mYfTMj2+2MsYwxvf7OOqF77mWMUZ/nbLxTw8tS8hgXpZomoevbFIub0+Hduy9I4R/OHT3bz+zT427D3On2YOpE9HzzuaLTpZxe/e28G6nGLG9unAk9f0Jyos2OpYykuIVb9wSk1NNenp6ZbsW3mudTnF/Pa97ZSdrmH2qHhuHRlPeKj7L7t2praef6fn88yqHCqqa5k3JZlfDO2mlySqCyYiW4wxqY2+pgNdeZrjFdU8lJbJ8h2HaRMSwG2XxnPLpXG0dcP1NGvq6vlgSwHPr8mlsLSSwd3b8cRP+pEQ08bqaMpD6UBXXinrUDnPrs5hZdZRwkMDmXVZPDcNjyMs2PozibV19Sz5rpDn1uwhv6SSgbER/GZcIpclROlRuWoWHejKq+0sKOPZ1Tl8nl1Eu1aB/O+ontw4rDutglp+sNfVG9K2F/Kn1XvYf/w0/bqEc++4REb3jtZBrpxCB7ryCdvyS3lmVQ7rcoqJaBXIuKQYrkiOYWRClEuHe01dPZv3l7A6q4iVWUcoOFFJUqe23DsukSuSOuggV07V7IEuIhOAP9GwYtGrxpgnznk9GHiLhrVEjwM/M8bs/7H31IGuXGXLgRO8uX4/a3cXcbKqlqAAPy7tFcUVSTGMTepATNuQZu+jvKqGdbuLWb3rKGuziyi37WdEz0h+dnEs45M76iIUyiV+bKDbPWwREX/gz8A4oADYLCJpxpissza7FThhjOklItcBTwI/a350pS7c4O7tGNy9XcOR874SVu06yupdR1mTXQRLYEDXcC7q3o6ObUOIaRtCh7bBxNgen33+vaqmjqLyao6erOJoeRVHy6spKq8i81A53+Ydp7be0L51EOP7duSKpIafBFq7wfl75bvsHqGLyDDgYWPMlbbnDwIYYx4/a5vPbNtsEJEA4AgQbX7kzfUIXbUkYwx7iipYldUw3HcfOcnpMz+8Qal1kD/tw4Ior6ylrPKH3exBAX7ERbbi8j4dGJcUw6Bu7bTaVrWoZh2hA12A/LOeFwBDz7eNbVHpMiASOHZOkFnALIBu3bo5FF4pZxAREmPakBjThjsu7wVARXWt7ci7quFI3HYUfvxUNeGhgQ1H723+/9F7TNtgwkMD9Zy4clst+vOhMeZl4GVoOEJvyX0rda6w4ADCosPoGR1mdRSlnMKRLpdCIPas511t32t0G9spl3AafjmqlFKqhTgy0DcDCSLSQ0SCgOuAtHO2SQNusj2+BljzY+fPlVJKOZ/dUy62c+J3Ap/RcNni68aYTBFZCKQbY9KA14C3RSQXKKFh6CullGpBDp1DN8asAFac870FZz2uAq51bjSllFIXQvvQlVLKS+hAV0opL6EDXSmlvIQOdKWU8hKWtS2KSDFwwJKdN08U59wB6wP0M3s/X/u84LmfubsxJrqxFywb6J5KRNLP16PgrfQzez9f+7zgnZ9ZT7kopZSX0IGulFJeQgf6hXvZ6gAW0M/s/Xzt84IXfmY9h66UUl5Cj9CVUspL6EBXSikvoQO9GUTkPhExIhJldRZXEpGnRCRbRHaIyBIRibA6k6uIyAQR2S0iuSLygNV5XE1EYkVkrYhkiUimiNxtdaaWIiL+IvKdiCy3Oouz6EBvIhGJBcYDB63O0gJWASnGmP5ADvCgxXlc4qwF0ScCycBMEUm2NpXL1QL3GWOSgUuAO3zgM3/vbmCX1SGcSQd60z0D/B7w+t8qG2NWGmNqbU+/pWHVKm80BMg1xuQZY84A/wSmWZzJpYwxh40xW22PT9Iw4LpYm8r1RKQrMBl41eoszqQDvQlEZBpQaIzZbnUWC/wP8InVIVyksQXRvX64fU9E4oBBwEZrk7SIZ2k4IKu3Oogztegi0Z5ERFYDHRt5aS4wh4bTLV7jxz6vMeYj2zZzafgR/Z2WzKZcT0TCgA+Ae4wx5VbncSURmQIUGWO2iMhoq/M4kw708zDGXNHY90WkH9AD2C4i0HD6YauIDDHGHGnBiE51vs/7PRG5GZgCjPXi9WIdWRDd64hIIA3D/B1jzIdW52kBI4CpIjIJCAHaisjfjTG/sDhXs+mNRc0kIvuBVGOMJ7a2OUREJgB/BEYZY4qtzuMqIhJAwy99x9IwyDcD1xtjMi0N5kLScFTyJlBijLnH6jwtzXaE/ltjzBSrsziDnkNXjngBaAOsEpFtIvKi1YFcwfaL3+8XRN8F/Nubh7nNCOAGYIztv+0225Gr8kB6hK6UUl5Cj9CVUspL6EBXSikvoQNdKaW8hA50pZTyEjrQlVLKS+hAV0opL6EDXSmlvMT/A+pesew6PFcEAAAAAElFTkSuQmCC\n", "text/plain": [ "