From 2854807cfdeaf2767ead164c261bed184765a611 Mon Sep 17 00:00:00 2001 From: CosmoMatt Date: Wed, 9 Mar 2022 20:02:12 +0000 Subject: [PATCH 1/2] add fb algorithm, pnp-prior functionality, abstract classes, and keras gradients --- notebooks/custom_operators.ipynb | 9 +- notebooks/fb-proximal.ipynb | 298 +++++++++++++++++++++++++++++ optimusprimal/__init__.py | 2 + optimusprimal/ai_operators.py | 144 ++++++++++++++ optimusprimal/forward_backward.py | 103 ++++++++++ optimusprimal/grad_operators.py | 124 +++++++++++- optimusprimal/linear_operators.py | 57 +++++- optimusprimal/prox_operators.py | 82 ++++++-- requirements/requirements-core.txt | 1 + 9 files changed, 793 insertions(+), 27 deletions(-) create mode 100644 notebooks/fb-proximal.ipynb create mode 100644 optimusprimal/ai_operators.py create mode 100644 optimusprimal/forward_backward.py diff --git a/notebooks/custom_operators.ipynb b/notebooks/custom_operators.ipynb index 7ba4128..2ea344c 100644 --- a/notebooks/custom_operators.ipynb +++ b/notebooks/custom_operators.ipynb @@ -96,7 +96,7 @@ }, "outputs": [], "source": [ - "class custom_phi:\n", + "class custom_phi(linear_operators.LinearOperator):\n", " \"\"\"A custom linear operator e.g. a custom measurement operator\"\"\"\n", "\n", " def __init__(self, dim, masking):\n", @@ -173,8 +173,7 @@ }, "outputs": [], "source": [ - "g = grad_operators.l2_norm(sigma, y, phi)\n", - "g.beta = 1.0 / sigma ** 2" + "g = grad_operators.l2_norm(sigma, y, phi)" ] }, { @@ -222,7 +221,6 @@ "outputs": [], "source": [ "h = prox_operators.l1_norm(np.max(np.abs(psi.dir_op(phi.adj_op(y)))) * reg_param, psi)\n", - "h.beta = 1.0\n", "f = prox_operators.real_prox()" ] }, @@ -242,8 +240,7 @@ }, "outputs": [], "source": [ - "# Note that phi_adj_op(y) is a dirty first estimate. In practice one may wish \n", - "# to begin the optimisation from a better first guess!\n", + "# Note that phi_adj_op(y) is a dirty first estimate. In practice one may wish to begin the optimisation from a better first guess!\n", "best_estimate, diagnostics = primal_dual.FBPD(phi.adj_op(y), options, g, f, h)" ] }, diff --git a/notebooks/fb-proximal.ipynb b/notebooks/fb-proximal.ipynb new file mode 100644 index 0000000..af8ed8a --- /dev/null +++ b/notebooks/fb-proximal.ipynb @@ -0,0 +1,298 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "collapsed": false + }, + "source": [ + "
\n", + "\n", + "# [`Optimus-Primal`](https://github.com/astro-informatics/Optimus-Primal) - __Basic 1D FB__ Interactive Tutorial\n", + "---\n", + "\n", + "In this interactive tutorial we demonstrate basic usage of `optimusprimal` for a 1-dimensional noisy fitting problem." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "How to run a basic 1D unconstrained proximal primal-dual solver. \n", + "We consider the canonical problem $y = x + n$ where $n \\sim \\mathcal{N}$. \n", + "This inverse problem can be solved via the unconstrained optimisation \n", + "\n", + "$$\n", + "\\min_x [ ||(x-y)/\\sigma||^2_2 + \\lambda ||\\Psi^{\\dagger} x||_1 ]\n", + "$$\n", + "\n", + "where $x \\in \\mathbb{R}$ is an a priori ground truth 1D signal, $y \\in \\mathbb{R}$ \n", + "are simulated noisy observations, and $\\lambda$ is the regularisation parameter which acts as \n", + "a Lagrangian multiplier, balancing between data-fidelity and prior information. Before we begin, we \n", + "need to import `optimusprimal` and some example specific packages\n" + ] + }, + { + "cell_type": "code", + "execution_count": 106, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from scipy.stats import norm as normal_dist \n", + "\n", + "import optimusprimal.forward_backward as forward_backward\n", + "import optimusprimal.grad_operators as grad_operators\n", + "import optimusprimal.linear_operators as linear_operators\n", + "import optimusprimal.prox_operators as prox_operators" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we need to define some heuristics for the solver, these include:\n", + "\n", + " - tol: convergence criteria for the iterations\n", + " - iter: maximum number of iterations\n", + " - update_iter: iterations between logging iteration diagnostics\n", + " - record_iters: whether to record the full diagnostic information\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 107, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "options = {\"tol\": 1e-5, \"iter\": 5000, \"update_iter\": 10, \"record_iters\": False}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we simulate a standard de-noising setting by contaminating a known\n", + "signal $x$ with some Gaussianly distributed noise. Note that for simplicity the\n", + "measurement operator here is taken to be the identity operator.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 108, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "size = 2048 # Dimension of the 1D vector\n", + "ISNR = 20.0 # Input signal to noise ratio\n", + "sigma = 10 ** (-ISNR / 20.0) # Noise standard deviation\n", + "reg_param = 1e-1 # Regularisation parameter \n", + "\n", + "x = normal_dist(0, 0.5).pdf(np.linspace(-2, 2, size)) # Ground truth signal x\n", + "y = x + np.random.normal(0, sigma, size) # Simulated observations y" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the unconstrained problem with Gaussian noise the data-fidelity constraint\n", + "is given by the gradient of the $\\ell_2$-norm. Here we set up a gradient operator\n", + "corresponding to a gradient of the $\\ell_2$-norm.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 109, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "g = grad_operators.l2_norm(sigma, y, linear_operators.identity())\n", + "g.beta = 1. / sigma**2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We regularise this inverse problem by adopting a wavelet sparsity $\\ell_1$-norm prior.\n", + "To do this we first define what wavelets we wish to use, in this case a\n", + "combination of Daubechies family wavelets, and which levels to consider.\n", + "Any combination of wavelet families available by the [`PyWavelet`](https://tinyurl.com/5n7wzpmb) package may be\n", + "selected.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 110, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "wav = [\"db1\", \"db4\", \"db6\"] # Wavelet dictionaries to combine\n", + "levels = 6 # Wavelet levels to consider [1-6]\n", + "shape = (size,) # Shape of nd-wavelets\n", + "psi = linear_operators.dictionary(wav, levels, shape) # Wavelet linear operator" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we construct the $\\ell_1$-norm proximal operator which we pass the wavelets\n", + "($\\Psi$) as a dictionary in which to compute the $\\ell_1$-norm.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 111, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "h = prox_operators.l1_norm(np.max(np.abs(psi.dir_op(y))) * reg_param, psi)\n", + "h.beta = 1.\n", + "f = prox_operators.real_prox()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally we run the optimisation...\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 112, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2022-03-08 16:47:37,249 - Optimus Primal - INFO - Running Base Forward Backward\n", + "2022-03-08 16:47:37,262 - Optimus Primal - INFO - [Forward Backward] 0 out of 5000 iterations, tol = 0.305357\n", + "2022-03-08 16:47:37,328 - Optimus Primal - INFO - [Forward Backward] 10 out of 5000 iterations, tol = 0.109128\n", + "2022-03-08 16:47:37,381 - Optimus Primal - INFO - [Forward Backward] 20 out of 5000 iterations, tol = 0.073595\n", + "2022-03-08 16:47:37,448 - Optimus Primal - INFO - [Forward Backward] 30 out of 5000 iterations, tol = 0.049439\n", + "2022-03-08 16:47:37,519 - Optimus Primal - INFO - [Forward Backward] 40 out of 5000 iterations, tol = 0.033118\n", + "2022-03-08 16:47:37,703 - Optimus Primal - INFO - [Forward Backward] 50 out of 5000 iterations, tol = 0.022142\n", + "2022-03-08 16:47:37,824 - Optimus Primal - INFO - [Forward Backward] 60 out of 5000 iterations, tol = 0.014785\n", + "2022-03-08 16:47:37,939 - Optimus Primal - INFO - [Forward Backward] 70 out of 5000 iterations, tol = 0.009865\n", + "2022-03-08 16:47:38,039 - Optimus Primal - INFO - [Forward Backward] 80 out of 5000 iterations, tol = 0.006579\n", + "2022-03-08 16:47:38,135 - Optimus Primal - INFO - [Forward Backward] 90 out of 5000 iterations, tol = 0.004386\n", + "2022-03-08 16:47:38,201 - Optimus Primal - INFO - [Forward Backward] 100 out of 5000 iterations, tol = 0.002924\n", + "2022-03-08 16:47:38,285 - Optimus Primal - INFO - [Forward Backward] 110 out of 5000 iterations, tol = 0.001949\n", + "2022-03-08 16:47:38,404 - Optimus Primal - INFO - [Forward Backward] 120 out of 5000 iterations, tol = 0.001299\n", + "2022-03-08 16:47:38,492 - Optimus Primal - INFO - [Forward Backward] 130 out of 5000 iterations, tol = 0.000866\n", + "2022-03-08 16:47:38,570 - Optimus Primal - INFO - [Forward Backward] 140 out of 5000 iterations, tol = 0.000577\n", + "2022-03-08 16:47:38,707 - Optimus Primal - INFO - [Forward Backward] 150 out of 5000 iterations, tol = 0.000385\n", + "2022-03-08 16:47:38,816 - Optimus Primal - INFO - [Forward Backward] 160 out of 5000 iterations, tol = 0.000256\n", + "2022-03-08 16:47:38,938 - Optimus Primal - INFO - [Forward Backward] 170 out of 5000 iterations, tol = 0.000171\n", + "2022-03-08 16:47:39,103 - Optimus Primal - INFO - [Forward Backward] 180 out of 5000 iterations, tol = 0.000114\n", + "2022-03-08 16:47:39,280 - Optimus Primal - INFO - [Forward Backward] 190 out of 5000 iterations, tol = 0.000076\n", + "2022-03-08 16:47:39,475 - Optimus Primal - INFO - [Forward Backward] 200 out of 5000 iterations, tol = 0.000051\n", + "2022-03-08 16:47:39,625 - Optimus Primal - INFO - [Forward Backward] 210 out of 5000 iterations, tol = 0.000034\n", + "2022-03-08 16:47:39,873 - Optimus Primal - INFO - [Forward Backward] 220 out of 5000 iterations, tol = 0.000023\n", + "2022-03-08 16:47:40,004 - Optimus Primal - INFO - [Forward Backward] 230 out of 5000 iterations, tol = 0.000015\n", + "2022-03-08 16:47:40,170 - Optimus Primal - INFO - [Forward Backward] 240 out of 5000 iterations, tol = 0.000010\n", + "2022-03-08 16:47:40,207 - Optimus Primal - INFO - [Forward Backward] converged in 241 iterations\n" + ] + } + ], + "source": [ + "alpha = 2 / (g.beta + 2)\n", + "best_estimate, diagnostics = forward_backward.FB(x_init=y, options=options, g=g, f=f, h=h, alpha=alpha)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "...and plot the results!\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 113, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEKCAYAAAAiizNaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAACDMklEQVR4nO2deXxU1fXAvzfLQBbIRgg7CUlYgiJC2HfBBaXYVnFprbgV7aLUWqu01Fq1xV/ditYqtC5Yd7Qi4g6yI0tQQAwEsgFhzQpkIZPl/v54817evLyZTMKQhOR+Px+YzJv73rvz5t577jn3nHOFlBKFQqFQKJqLgJaugEKhUCjaF0rwKBQKhaJZUYJHoVAoFM2KEjwKhUKhaFaU4FEoFApFs6IEj0KhUCialQYFjxDiFiGENP0rE0LkCiE+EEJcJ4QQTbmxEGKoEOJhIUR0U873ct04IcSzQoh9QogKIUSBEGK7EGKhEKKDqdwa1/d5zeYad7g+izcde9XyHMqFEN8IIW47i7peIIRY5KqfUwhh69suhOgkhHjSVedTrvtPbsR9XhFC7HGdWyqE2CmEuFsIEejlnLFCiFrXvYJ8uEe8q+wtpmMPW55ZtRDigBDiJSFET1/rb7nPZJtrHhRC/EsIEdWUa7Y2XM/yYSFEv9Z2f1fff7UF6jRUCPG+67euFEIcFUKsFkLcYymnt4t6/VII8boQItdyLNfSnkqEEF8KIcY3sZ6dhRAPCSE2CSEKXdfbJIT4oU3Z8a5xZberHefWv6LH+/g0JtiM39Z/3Rq4z2TrtW3GQqcQIksI8ZQQItKX+jdG45kFjAGuBP4EVAJvAV8KIUIacR2docCfAb8JHiFEZ2ALMAN4Bq2udwKfAD8A7Or5UyFEio+3yEd7BmOAG4ES4CUhxDVNrPJwVx0PAmleysUAtwHVwJdNuE8I8Bzab/hjYCWwEHjarrAQIhhYBBxvwr3sGI/2zKYAfwOuAj4WQpyNxn2P65qXAf8F5gD1JhHnKfFofaNFBE8D9/8R8GhzVkYIMQLYDHQBfg9cDtwPZLjqY8efhRAOH2/xOVpbGgfMBRKBT8wTz0bQB/glsBa4Cbge2Ad8IIT4laXsVGAC8D2wp5H38XVM+Ji6MUv/NxYoBLZJKY818r465rHwUuBfaGPtf306W0rp9R9wCyCBJJvPrgFqgecauk5jrtvUf2g/hAQusvlMAML0fg3wLdrg+r6l7B2u68Sbjr0K5FnKhQPFwGdNrG+A6e/HtJ/Dtpy53tNcdZt8ls/qLeC0h8/+AOwG/uq6V5AP14t3lb3FdOxhu/NNz3dQE+o92XXuNMvxf7uOd/NXe2qpf56+o6e2ATha6v7N9DxeA44BHWw+C7C8l2iCRAJ3Wz57Hci1HMsFXrccG+c6/8Em1DUMCLU5vgo46KnudnVr6Hc3/d2oMQFN2EngV41oC5NNx+qNha7jjwE1QFhD1z2rNR4p5fvAh8DPhRCh+nEhxF9cZqhTQjN1fSWEGG36/BbgFdfb/SaVLd71+a+FEF8LIYpcqupmIcRVPlRJ157qSXHpwnK4DG0G/mMhxHAfv7b5mqVos5k+jT3XdX6tj+XORXqJQrTZkhtCiERgPtqsrcruRCFEqMu0Vegy3S0HejXi3qdcr8GNrLM3vnG9Gr+FECJICDFPCLHXZZ454jIHdDSfKIQIE0I87jIXVAohjrnMOnGmMiOFECtd37dMCLFKCDHScp1XhRB5QoiLhRDrhWaO3S+EuMtSrpsQYomrPrrZaIUQoqvLpLHaVfRLU9+Y7Do312Uyuk0IsRdwAlfZmURc5XVTS7zl+M9dfbRCCFEshFgrNPOqL/d/1XItvz0bD0QDxVLKSusHHvrQNmAZ8EfzuNQI6rUlX5FSlkkpy20+SgN6WMr61P893OdsxoTZaO3mLfNBIUSsEOJN17hdIrRliMhGXPcUmhXNowlfxx/OBZ8AHYBU07GeaKauq9E0mxPAOiHEha7PP0aTjlBnwhsDHHUdiwf+4/rserQfbYUQ4ooG6rLV9fq2EOJyIUSYD/V/Ec3U9VcfyrohtDWS3kCW5fga0Qh7bXMgNIKEEJEu0+Bs7E1tLwJLpZTrvFxuEZrW8jSa6S4DeNNL+UDXvUNcAv4PaOaF3U35Lh6IR5tt5ZqOvY4mRN9EM+8tAG4H3tALuMwxXwJ3o83kZgC/BoqAKFeZIWimkyi09nwz0BlYK4S4yFKPzq77vY7W/rcBLwghppjK/Betvd+PZqa4B8gDQtEGPd0ko5sTx1A3GIJmsvwt8BfgCmBXQw/HjBDiSWCx65rXoZmE1qENtL7c33wtfz8bO7YCA4UQL7qEXINrjmi/e6zrOzSWeNertV9Lq9BtBBOBvU08128IbVlkFrBCSllk+fh/aO3/D2jjbjWaid7TtYJc/8KEEBPR+s1nUspTns4x8EHVugUvJjE0e6sErvfweSAQhDY4LfT1umZ11HX+F8CHPtT3ITRpLl0PLg3N5BNpKbcG2OD6+3ZX+Qmu9x5Nba66BKHNXp5D05pG2ajVmb6qzSY1VfpQrkmmNleDkq5/tcDfbMrchDbgdnW9fxiLqQwYgDbAP2g59wU8m9qs//YAiY2pv+mak13XuMz1O3QCfog223rSVE43J9xsOf+nruNDXe918+xML/d8D209L9J0rLPrWf3P0kYkMMV0rAOadrnYdKwUuMeH71jP1IUmWMuxmBSxMYlY+lm8632S6/d7+izu/+q5ejYe6hMCfGBqP+Vo48HPsTe1Peb6+7+uekS43nsytb3haksOIAVNkO4Doixlq4GXmtBm57jq9VMvZRplarOc6/OYgLY2Xa+9o02AJHCD5fin1mubfkvrv6+BLr7U2R8aj+7VZqh+QohpQvM40c05VUB/tEGr4QsKMdxlejhuOv9SX86XUj6CNnO7A63hxaAtlO42m04svIrW0P7WwOV7uupSBRxGmxneJqXcYqnDVCllUkN1bWbWAyPQGunjwO+EEIaWJzTvwqeBP0gpT3i5zii0ycC7luNvezlntOveo9Bm2GXAF15+D1/4HO13OIU2KK1D0yB0rkCbgLxnmpnpExjQZqCgCbBjUsrlXu41EW2GWKIfkNqsbjkwyVK2XEq52lSukvrm2G3A/UKIuUKIC4VotGfoZtn0ReFpaL/f4iaeb8Xfz6YeUsoKKeWPgMFov/GnaBaWxcCnXp7fn9HWYe/38LnOT9DaUiWaJn4B8AMpZbGlHkFSytsbuJYbLhPls8BrUso3vJduFmajWaA+sRwfgzYhed9y3FO/PoHWp0e4zp2N5vzxqfDB2cwfgqe36/UogBBiGNqXKkXTJPRBZyfQ0e4CZoQQvdE0hmg088dY1/mf+XI+gJTymJTyJSnlrVLKBDQVsCceGqCUsgZNUxovhJju5dL6wx6F1lhzgJeFEAN9qVdLIqU8KaVMk1KuklL+AU3IPijq3JofQ/sN33WZ4yKpe94RJrNld9er1ePNmwfcdte9t0opl6KZvRLQzEVN5VfUCdJ3XNf8k+nzrmgz2DLqJgtVaL8haBMS/fVwA/eKps4MbOYYLnOciWKbcpW4t93r0Qbm36OZyQ4LzQXX1/5oVxdf0b933llcw4y/n41HpJTpUsonpZTXoFkcXkebONiu/0ops4GXgLlCiFgvl/4UrS2NBX6DpmH9T1jWAhuL0LzxlgNfoU2EWxQhRHe0/vKmlNK6vtsdbR3Nuq7rqV9Xufp0mpRys5TyNbQxMRVNy/aKL7bShrgKOANsd72/Bk1L+bH5SwgtxqLEh+tdAUQA10kpjc7RxEVCAKSUzwshHkVToz3xLvAg2gD8oocyVVJK3e15qxDiG7SB4yk8NP5WTBraxCMBbeBNAYagmT6sFKA5kfyQukEmDsg2lfFZe5FSHhdCFLju11T26b+FEOIr1/3nCSFekVIeQvseZ9BMbnYccb0WoM1wvVEE2MU7dMN+MPWKS6P8FfArIcQAtNniX9BcVF/w5RI2x864Xq0uxDGW9wWu155o5u+zxa/PxleklGeEEE+gmYdTgBUeij6K9nz/4OVyRaZ+/bUQ4iSa89PdwBNNqZ9rPftzYAdwjc2A3hLchLb0scTms6NAlBAi2FLXxlglvne9Ntivz0rjcS1SzwRelHWeHKFoKpvZ9HYJ9dVp3UPFqpbpAsYstPqjuTg2VJ84u1mjS9JH4GWmKDXj5XxgGJrwbBApZQbwPHCla3ZzPjEJ7TfShcdv0Batzf/0BjoN7dmAFidVi2YyM3ODrzd2/R5d0Abas8b1292LtmbwoOuwriFHmGZm5n+64PkC6CaE+IGXW6xF+407mb5DJ7TYsDVnWfcMlwZaTJ0A9NQ3vHHA9WoVotYJ0Uq032+Ol2s15v7n7NmYrtfdw0e6pcFbvz6C1kd/ge+el0vQnCnub8qEVwiRjOawkg3MkFJWNPYa54ibgV1Syh02n32NJpSsY5/P/Zo6gdNgv26MxjNUCNEFbUbVB22xehbaA55nKvcZ2iD2qhDiFbS1nT9R35yR7nr9lRBiCZqg2YXWMaqB14QQT6GpgH9B8zxrSFD+DJgjhHgDzROm3HX/+9Ds/c97O1lK+bEQYiOaw4SvPI7WiR9C62wIIVYBfRta53E16itdbwe6jl3rep9rmoXhMgGGAbpn4CTX71EmpfzUVC4TOCClnOp6fxVwK/AR2jPsBEx31XmRPgDbNUZR55q7VlfNpZQZQog3gUdcQn4bmrnjSuv5JkYJIWrQfr++aCbPGkyapahz471VSvmql2vZIqXcIYR4H7hdCPFXKeUaIcRbaGs8T6O1h1o0j6UrgQeklPvQzDU/B94SQixAE6yd0NrAP6SUe9FmzTOAVUKI/0MT2A+gTZIeaUw9hRARaG38DTQvpyo0D68o6taf9qH1gduEEEVogiBDSnnay/c/KoRYi6b1FaCZFG/CEgQqpcwSQjwD/NYlIJaj/RYjgb1SyncaeX+/PRsvLBZacPj7aJ6QgWimsd+jeZ590MD5eh+dRJ2A9oiUUgohHkLTon6BZtFACFENLPG2ziOE6Io2JjrQ1phSLEtQ37rWtnCZ//R1sD5AqKn/p0sp013lJqEtP9zmMmnp9/JpTHCVHYY2KbnPw3f+UgixAVjkusZ+NJOwJ2uAQ9SFyAShjbN/BE6jrZl7pyHvA+q8YvR/FWg/3gdogkfYnHM32vpHBdrANA1t9rPGUu7PaAJJ15B0z5vr0DrlGTT17QbXl/Hq9QEMQnPj/hbN1FKFNht6DxhmKbsGl1eb5biuCRj1kXWeHPWCplyf/c1V/mLTtRv0UKEu6NLu36uWsrkeyuXalFtjej8QzU3yENoAchzYgObdFdBA/R7GPgA0FM0kVIS2lrecuqC7W2zO1//Vopm4PgJGWq55lavMFQ3UaTKePa4GudrSQtf7ALRI9J2utnTS9fffcXk6ucqFo5lUDqBNUPQ209VUZhSawChFWzdaZfMdbNsIpraPppUtQmvXpWjOEduAn1jOuRNtxlyNyasIm4BH0zm9XM+2BG2N5W/YeGi6yt6FNtGrdP2Oa4AxPt7f2jb98my8/OaXo2khGWgDWyWawPknEGcpK3F5tdmMNZ76i6fnucnVFkJM1361gbpOxnOfto4p3so+bFPuFpu6NzgmuMouRBsP47zUPRZXYLmrDb2GNikyfn/Tb2m+XzVa33kLGOjt+ej/hOtCCkWLIoT4G5rZ9kKpGqVC0aZR2akVrYVJaLFFSugoFG0cpfEoFAqFollRGo9CoVAomhUleBQKhULRrPgjgPSc0qVLFxkfH9/S1VAoFIrzhu3btxdIKb1la2hRWr3giY+PJy3N2x5pCoVCoTAjhGgwXqklUaY2hUKhUDQrSvAoFAqFollRgkehUCgUzUqrX+NRKBSKplBVVUVeXh5nzpxpuPB5SseOHenVqxfBwf7cRf7cowSPQqFok+Tl5dGpUyfi4+Np/F57rR8pJYWFheTl5ZGQkNDS1WkUytSmUCjaJGfOnCEmJqZNCh0AIQQxMTHnpUanBI9CoWiztFWho3O+fj8leBSKc0RRmZNFa7MoKnO2dFUUilaFEjwKxTliadohFny6l6Vph1q6KopWwMMPP8yTTz7p8fNly5aRnp7u8fO2hHIuUCjOEbNSe7u9KhTeWLZsGTNmzCAlJaWlq3LOURqPQnGOiA5zcOekRKLDHC1dFUUL8de//pX+/fszfvx4MjIyAPj3v//NiBEjuOiii7jmmmsoLy9n06ZNLF++nPvvv5+hQ4eSlZVlW66toASPQqFQuPDnutz27dt5++232bFjB5988gnbtm0D4Mc//jHbtm1j586dDBo0iJdeeomxY8cyc+ZMnnjiCXbs2EFiYqJtubaCMrUpFAqFC31dDuDOSYlnda3169fzox/9iNDQUABmzpwJwO7du5k/fz4lJSWUlpZy+eWX257va7nzESV4FAqFwkVzrMvdcsstLFu2jIsuuohXX32VNWvWnFW58xFlalMoFAoX/lyXmzhxIsuWLaOiooLTp0/z0UcfAXD69Gm6d+9OVVUVb7zxhlG+U6dOnD592njvqVxbQAkehUKhOAcMGzaM66+/nosuuojp06czYsQIAB599FFGjRrFuHHjGDhwoFH+hhtu4IknnuDiiy8mKyvLY7m2gJBStnQdvJKamirVRnAKhaKx7Nmzh0GDBrV0Nc45dt9TCLFdSpnaQlVqEKXxKBQKhaJZ8ZvgEUK8LIQ4IYTY7eFzIYR4VgiRKYTYJYQY5q97KxTnGpX+RqHwH/7UeF4FrvDy+XQg2fVvDvCCH++tUJxTVPobhcJ/+M2dWkq5TggR76XI1cBrUltU2iyEiBRCdJdSHvVXHRSKc4VKf6NQ+I/mXOPpCZini3muY/UQQswRQqQJIdLy8/ObpXIKhTesbrbK9KZQNJ1W6VwgpVwspUyVUqbGxsa2dHUU7ZCGBIsyvSkUTac5Bc9hwGyn6OU6plC0OhoSLLNSezNv+kBlelN4pKSkhH/9618tXY1WSXMKnuXAzS7vttHASbW+o2itNCRYVOZpRUN4EjzV1dUtUJvWhT/dqd8CvgYGCCHyhBC3CyHuEkLc5SryCZANZAL/Bn7pr3srFC2BWudReOPBBx8kKyuLoUOHMmLECCZMmMDMmTNJSUkhNzeXCy64wCj75JNP8vDDDwOQlZXFFVdcwfDhw5kwYQJ79+5toW9w7vCnV9uNDXwugV/5634KxbnElyzF/sxkrDjHPBxxjq570uNHjz/+OLt372bHjh2sWbOGq666it27d5OQkEBubq7H8+bMmcOLL75IcnIyW7Zs4Ze//CVfffXVOah8y6GyUysUaNrL0rRDzErtTXSYwyf3aeVirWgMI0eOJCEhwWuZ0tJSNm3axKxZs4xjlZWV57pqzY4SPAoF9bUXfQ3HilVAnY2mY72W4hziRTNpLsLCwoy/g4KCqK2tNd6fOXMGgNraWiIjI9mxY0dzV69ZaZXu1ApFc+Orl5o3b7fGrvno11qyKUetFbVBrNscmImLi+PEiRMUFhZSWVnJihUrAOjcuTMJCQksXboUACklO3fubLY6NxdK41EooJ724kkbsZrXzOUau+ajX6PcWaPWitogMTExjBs3jgsuuICQkBDi4uKMz4KDg3nooYcYOXIkPXv2dNv24I033uAXv/gFjz32GFVVVdxwww1cdNFFLfEVzhlK8CgU1Bc0dkLEThiZyzV2zUcXdkVlTkIdgWqtqA3y5ptvevzsnnvu4Z577ql3PCEhgc8+++xcVqvFUYJHoaD+Go9ZiOgCp7DUyeL12ZQ7q5k9NoGlaYeYlqLNYqelxBnvG7tuc7ZrRQrF+YYSPAoF9bUVszBYtDaLBZ/uZXxSjKu0qCeo9DKbswtZnZFvHFcoFPVRgkehwLvWoQujaSlxrEw/zrSUOJbvOMzcqcnGZyPio0mMDWP2mHhG94vxajZT3myK9o7yalMocPdIs3qn6UIpMTacOyclsjL9OAtXZQKSpWmHKCpz8txX+8nKL2PJ17mGo4FKMKpQ2KM0HoUC9zWecmcNC1ftp9xZw+yx8SzZlAMIZo+Nd31ezdypyQDGOfNnpADpzJ+RYlyr3FlNqCPITbMpKnMa5ytnAkV7RQkeRbvCFzdpTdCArtFo2g2EOgIBWLgqk3GJMQzuEWEIkOgwB6/cOhKAqFTtunZu0vr1pgxQ230o2i9K8CjaFZ5ibcxrPLPHJhiaCmgaDgjjve5AsDGrkLlTk+oJMm9u0iPio4kOC2Z1Rj5L0w4pB4Q2zrPPPssLL7zAsWPHeOCBB3jwwQdZtmwZ/fv3JyUlpaWr12IowaNoV/gSa2MWHEvTDjF7bIKbqWxIrwh6R4WSXVBGRVUtC1fVN6t50qye+2o/RWVVJMaGuQWhms15yuGg7fCvf/2LlStX0qtXL+PYsmXLmDFjhhI8CkV7oTExM9bg0KVph1zrP5qpbENmAcP7RjJv+kA3s9q0lDjmvJZGVn5ZPYFkXgsyB6GazXlKC2ob3HXXXWRnZzN9+nRuu+02srKy+MlPfsLy5ctZu3Ytjz32GO+//z6Jie3v91aCR9Hm8MVd2a6MtwzVSzblsHBVJnMmJDBv+kCmpcQxut9xNw1HN6vd9+4OsvLLSIwNA4QR3/PUdUNJjA031oJ0ZqX2rmfOU/iXC5dceE6u+93s7zx+9uKLL/LZZ5+xevVqIxfb2LFjmTlzJjNmzODaa689J3U6H1CCR9HmaOpeOnYZqnVNp6JKyyQc4tJedNOYjtk8l9w1nKoayV+uHkxUqINdeSWszsjnvnd3MH9GCivTj9dbE7r30gHn5FkoFK0RJXgUbY6m7qVjTZOzZFMO2w+UsCGzgLlTkwxN5753dxjZCXTTmK4tlTtrWLw+h3nTB5IYGw7AU9cNNZ2TrjIbtADeNBNF86MEj6LNYdY+Fq3NsjW5mdd6zCY261oOwJQBsYaDwaK1WazOyCc+JpTLBndjWkoci9ZmGbE/uoAyC7ToMAdPXTfUyOWmm+ispr2mmggV5xfetktoL6jMBYo2i12GALsMBUs25RrllmzKYcGne6lwVjNnQgLjk7pw9yXJLNmUwzNf7mNaShxTBsSSW1hOTJiDlenHXeY5ybzpA5k9NsEtx5un7AfmDNh6/XzJaKCyHpz/3HDDDTzxxBNcfPHFZGVltXR1WgSl8SjaLHbmNPM6DmiZB8xaypJNuYC2lhPqCGRDZg7BgcIwj+3KK2H+jBRDa4H6cT5mDWjdvnxS46OZObRHvbUda/2mpcSxbl8+hS6haKfRqO22zy9yc3MBuOWWW7jlllsAGDduHOnp6S1XqVaAXwWPEOIKYCEQCPxHSvm45fM+wBIg0lXmQSnlJ/6sg6J9YmeC0p0DlmzKocJZS4gjkEn9Y5kyIJYR8dGs3XeCuVOT3WJnZg7twa68EmYO7UFUqINyZzUVzlqSu3Yi/egpVmfkM7rf8XrrMwtX7QckoY4gQ5hNGRBrBJrqDga6e7WecNRc35Xpx9mYVcjGrEJibNy+lZlN0Vbwm+ARQgQCzwOXAnnANiHEcimlWbTPB96VUr4ghEgBPgHi/VUHRfvFkyebOUYGMATAwaJysvLL6qWuWZl+3E24hDqCWLhqL/OmD+TZGy82Bn6zEKjzbhNurtHzZ6QwpNdhKpyaR9yQXpEAHrdPaMiturE7nCoUrRV/ajwjgUwpZTaAEOJt4GrALHgk0Nn1dwRwxI/3V7RjpqXEsTm70NiYTUcfzHWNZ+bQHuieZYmxYfVS11hNWfrrlXFF7H73Faad2kHEzmM4T53g+qpq5Ncx/DquHxMHJpPcN4LOocmG1qNrQKCZ3XRzXqgj0M3JQKcht2plZms8UkqEEA0XPE+RUrZ0FZqE8FfFhRDXAldIKe9wvf8ZMEpK+WtTme7AF0AUEAZMk1Jut7nWHGAOQJ8+fYYfOHDAL3VUtF30jdjmTR/YoDaQlV/KYyvSufuSZLblFtULIjXS14zuRXTWMtiyCI7u8K0i0f0ou+hWniocywe7iygur2Lu1CSf0uko/EtOTg6dOnUiJiamTQofKSWFhYWcPn2ahIQEt8+EENullKktVLUGaW7nghuBV6WUTwkhxgD/FUJcIKWsNReSUi4GFgOkpqaenyJd0azYaQNmIWJe3NfNaUN6Rbg0kjp009ykgJ3M3vEelGleRxWBndgUPJoNNRfwdWkcPxw/lO+PlpKelUOSOMLogHSuC91OWFE2Yav/xC9FFMWVN7A2dAozh/ZkZfpx4x5LNuWycNV+1u8v4NkbL7YVPko4nT29evUiLy+P/Pz8lq7KOaNjx45ueeDOF/wpeA4DZhtAL9cxM7cDVwBIKb8WQnQEugAn/FgPRTvD0yBtXt/R13YAkzdajVs6m+gwB7OGRDHquzcZWrACyuBkh+68HXojz5R0pjboEAGhWcT0zGBvSBZTJ06jX99JVDhrKXEEUjm6F7VZX3Jm5d+ILd3LM44XONl9D+9sCWPBhpOUO6tdpjRtLrUhs8Aw81m/g7f1HCWUfCM4OLieJqBoHfhT8GwDkoUQCWgC5wbgJ5YyB4GpwKtCiEFAR6DtTkcUzYKnQVp3T07sGk7HoECG9Ip004h0D7bVGfnc/eY3XNb1JDflPMDQkhwqZTAfRN3MwxVhBIasITiywDjvFLDq0G7WHf6KO3q/xeyxiYYAWHTyAhYU/IlFF+xlcu4/iDi8hp+e2MVK8SsqnP1YtDaLSf27kpZbTGLXcMqdNYZmtnBVJuXOamaPTfC6WZz+fc0CU0cJJcX5gN8Ej5SyWgjxa+BzNFfpl6WU3wshHgHSpJTLgfuAfwsh7kWb9t0iz9fVMUWL4S2ZpxndPdkRFMDqjHzGJcYYn2nuz3D3JcnszCshMOcrfpT3LIGigurYwbzU+xesrP2I4JN7AIgL7UZc4HC2ZIQwa+hAVpx4jKraKhbufJwDpTNJCL+ASwd3dwmM/myvSmT+7m68H/cKfU6m8XbHv7GyWDJnfQpTBsQa9Xrt6wOuDebqPON0TW3e9IEeY3l0rzjrnj7K801xPuDXNR5XTM4nlmMPmf5OB8b5856K9oevg6s5MFP3ZNM3b9M9zO57dwdjKtaz0PFPgqilNPFKvplwG//d9BClVaWEBcbwm2G/4dqBMzhVUcvSGE3ghezKYOm+pTgi0/ikMI3aY51ZmjuSnAPJ/H7SpUAA+UTxv8HP8hveJmDTs1yW+RgL4u5gwCXzGd0vpp5nm3XTOE/ea+YUPNYyyvNNcT7gN6+2c0VqaqpMS0tr6WooWhFWjccXjzZPm60d3/AasSvnEkAtL1bPYE2/H5Ie+CgSSWLoGHZ8ewVzpwx2XUU7FzThNyz5DOuPfs6n2Z9zrLwuMiAutBtTe19O7emLuf6ikXy08whR37/GzcXPEyAkafFzSL3liXP4hBTtHeXVplD4EV3oTEuJc0vsCQ3P8vWMAYbQyv2EuJX3AJKiEb9l6Z6JHDyxho7dJRN7TeSRUU/zXtc8I/2Ndo1AQAsCncdAfjvpt8we8Cue3/QVInwnaw+v5Hj5Md7MWAIs4ePjPTi2/xpqK8dyoFMwf6paSGruYsrX9OTMiF+q9RhFu0QJHsV5hXlh3Rz5r2eV9hQrY3Zh3pBZQFH6ah4s+AMCCZP/wFJxLVkFexmQ7OAIEFwbw++W7nTtGAppuUUM7hnhJtzMmalf/qqWedNn8cW1v+eb49/wSc4nfHngS0oqjxAcsZPKEz155fQIOkbezQNnFhK65s9sPVrDgp0Dje+gULQXlOBRtGo8ORJY10esnl7114E0k3JK9070rTnAr478CSEqOXnBbN7mGqalxFHurObTQ5pZd+P+k+QfzCfzxFb6xoSxMauQif1jDc3kzkmJhonPvGZUUl7N9owofp36IAOjB/Lo5kcRwUVc3CccR6CDF3JG0a3bL5ld8i8mZvyVJ4c/xyVqPUbRzlCCR9GqsdsVVNcOEieFG+Wsnl5mAbVobRYzh/YEoEPVKe4ve5RgUc5nNSN4Kvdq9qdlUF5VQ1puMbkni+kQC32jI8g/CIeKKzhUXMGUAbHGtey86aLDHHxzoJg7XttGUVkVAP0TtDxwwZ13czz4DwSfGYpwjKBo8K2k5Z4m9fB/mZHxIB0vGQ1hdRqPcolWtHWU4FG0avRca3q8i6eB2OrppQsoXSspd1bz3aEibs55gODAg1THDeF/jj+zf/8pxiXGsP1ACZuPryek53oApiYPYFinBHbmneSi3pHcNSnRqxAsKnNy25KtlJRXExUazKzU3nQK6cO9w+/lo6yPyCzJhMB1hPfbQGHHA3zV+xaKD37PpXwDb14PP/8KOmppDJVLtKKtowSPolUTHeaol3TTkyZgFgQ65iwFQ7NfZHLQTmpDogm68Q0eqIqmbNluyquq+L7sA0J7rQQhqTp5MfL0cGLCg9mSU2Q4FOjXGhEfza2vbGX+jBRje+ulaYcoKa8G4NJBdQ4Mt11wG7ddcBtfZe/kL2sWURy4gQ+z36Nz4CZODfsNF+//C10K95P58u1E3/w60eEdGnSWUBqR4nxHCR5Fq8ea4gYa1gTMg/OdkxI5/d3HdApahhQBBFz7MkT2YeXaLDbmHCKk5zt06JoBCO4Y/Cs6lE7juhF9AOoFat45KZFbX9nqcmxId9vSWt/SAKhXz/2Hwjmw7yoCHKl06beUUzV5rDj5N9ac+SFf8DxJJ77g0/eepP9Vc3lo2W4Su4azZFOu4b5tdgVXGpHifEcJHkWrx7yhmzmNjLeZv+7FVu6s4ZaLwuiwXEuSvjXhlyR3GwdlTnJL0+k26AXKavPpIMJ5dNwCpidOdruOXaCm5ummZbe+790drM7IN5waAF5ck8n4pC6MiI821oTq9tpJZnTidH6/7veUiB2c7LOMOXlX83r1O0zNfYbH/pfAxpwINmYVuu4mScstNt6bg0xVkKjifEUJHsV5gTmNDLhvLw3uM/+iMidpuUXaG1lL8Vs/p19VEd93uIgb0kfzQJ+D7Dr9MesKX0GIGi6IuYAnJz9Jz/CexvnWtSIzibHhvHLrSBatzaq3rw/A4vU52q2lZGNWoZEcVN9rZ9HaLA7tnUXCwE4UiPXs6JnGc/kTuLt8PfMrnuRQwtP07taFqFBNmOpCZ3xSFzfHBruAWIXifEAJHsV5gXmWr5uazG7MZpamHWJjViFTBsQyJ3QNYcUbKZFhfDXgEe6NimZ39T9ZX7QKIeCapBv44+jfExwY7HZ+Y1Py6FsuAIbJraKqxiU0hO15lwyazG9WzuNg1TqW9i5lTkF/OhTu45m+H/FK6BwA18Z1EquAMWfeDnUEKpOb4rxCCR7FeYFZ87BuEW3NzlzurGbOhH70qD5E6Jo/A7DtgocYNaoLf958P3mleYQFhTEh6lfcM/R6N6GjX9/86kudolIdhpY0e2wCS9MOMXNoD2JMbtd258VW3sQBsZ5iZyHv9X2I6wvmELnzP2xxdmNzbYoWuNojghCXg4O5jt62yVYoWjNK8CjOO7x5uumawCX9Y7gs94+IgDNw0U+49Ed3cc2Hs8grzaNf5/6EltzK0nWCpLBD9bQFq5u0Nw8y/XOz2Q/qOxd4utafZlzItR93oppTjBo3nh1Fd5Cau5hFnV7insjnWZtVaJjaduWVGNsgNLRNtkLRmglo6QooFHYUlTlZtDaLojKn7eezUnu71nsECz7da6yv6Mf/r+82UgP2URsWB1csAODAKW1fwo6Fc/g6Q7gFher3sd5XN7vp17dSZ5aThtlPr4PZCWLR2iyWbMplwad7ue/dHcb1E2PD6RquuWRX1lSysfstHA0dQETlUe6uftW4T1RosNs60tk+P4WiJVEaj6JV0tA6i66VFJU53Ty9osMcXJ8MYS/9FYCAGU9TVBvK0rVZVNdqg/DgbrFckhzhtl6k38f6viGzmzV7gY7dHjlzpyYxZUAsqzPyWbIpx0haWlaprQG9vOk73t8cwApxOx875pFauJwRYigZHS6guLyK8Uld6mVPaOrzUyhaEiV4FK0SuwHfaqoyb3VQXO76bHgvTr37SyJrKtgXM40ufS5zuTyfoNOgSgAC6UBhqZMlm3KZ1D+WcYkxFJY5KSpz1ruvnVebmYY+t/suumlOzy1XcCqO4MjDrCp4AcQcouOH8E3QrYw+9B8WBP+HK88sAIIZ3jeSlenHjUwM1mBa8/PQnBKUy7WidaIEj6JVYjeg22knumeXvoV18rEVXFKymRIZxqe972XHuztYve8wfQd+QhEQEhDBvzfkGtdcsesIWfllbMwqJMZ1T39qCHbrOmZNbVpKHEu/uZulR+/HGXyIhAGf869r/0FMh4up+ddqkoqzuCvwI7b0uR0Qrk3t7INplaeb4nxBCR5Fq8RuwLZqDtNS4li3L5/BPSO4fkRvJvbez6Q0LVB0S//fURvWlTXbviVuwOsUiVxCgkL408g/k9mtDxXOWtKPnmJDZgHjEmNIjY/22YzVGDyZvMyCdc74IZSufYAV+X+kQGxg1eEPmdbrh2zo+wAzi+cw1/Ehr3W/gb+sKkbPsl3hrDWCac17FJU7q6lw1jaY206haEmU4FG0SpZsymHhqkwj+BLqa0Er04+zMasQR1AAUaEObq18Hc4Usbl2ELk9ZzJ9cAc+yH+V4qo8eob35LExT7P++0BAcNfkOi1Bj8NZvuOIbUBqQ3jzfPO0RmQ+Z2naIV5fV8OQAT8jJ+AlFmxZwJvrK9mRGU1k9KVMLP+Smwr+gfOKhRSWV7F4XTYA86YPdNuBVTe/geZhp7QeRWvFb15tQogrhBAZQohMIcSDHspcJ4RIF0J8L4R401/3VrRFhOW1PrNSexuL9V+t/hLSXkaKQA6O+gtjBsLctXdQXJVHTHBfnpv0Etv3d2DhqkwWrtrP0rRDhiDT103Mnmm+YPVWs/M4MwtLO2+5JZtyKHdWMz6pC7sykhnQ8UqqZTX7xb8QQSXcU3QtFcGRBB/cwJ3R3xASrHXZUQlRhlZj9fBr7PdQKJobv2g8QohA4HngUiAP2CaEWC6lTDeVSQbmAeOklMVCiK7+uLeibTJ7bLybtxq4awmgDd7zZ6QwOuEoV+/7BchaKoffxX5HEM+vuYOiykJCa/uRu/smVsSWMXtsvG3QpSfPNDusmoq3DApmPHnL6fE/c6cmMSG5C+v2R1Ht3EdQWCZdEt/kB7F/ZX3J3VyW+Sh88Sdm3/41oY4gU9yQdN1BcygwPzOrFqayWitaC/4ytY0EMqWU2QBCiLeBq4F0U5mfA89LKYsBpJQn/HRvRRvEzrnAvH318L5RdWaxiC1wZBuEx/FU4EjeOjQPEVhBr45D2LPjGpAdAGkEXVoHYF8803TMAqQxAsuTt5zZHTw6zMGI+Gh+s/R2nJ2eoYyD7DizmFd2X8Hm2AuIO72b6LSF3HnpX4zzzIGrablFpMZHU1zu5LEV6W5bg1vrrkxwipbEX4KnJ2C2M+QBoyxl+gMIITYCgcDDUsrP7C4mhJgDzAHo06ePn6qoOF/RBUWFswaADZkFlFVW0TsqhML849Ru+TMBwIcXX897hx9HBJ6hT8cRPDnpCT6LLkDPc6bjaQD2RSOwChtfB3BPZa3Ht+UWcbBAMOei+Sw78SD7y9YzdkQsVQMeh/dmIL9+HjHsZuiopcwpLqtiVEI0oCUT3ZhVaHj4TRkQ61GzUyhakuZ0LggCkoHJQC9gnRDiQillibWglHIxsBggNTVVWj9XtC/MJq25U5PYfqCEDZkFAPTc8S9kUAGP9BnE0rxlAETUjOb7b3/Aum4ltmllPA3AvmgEjRE2DWEXl1TurGbu1GRmj43n1Lp7WXF8AbtK/8dfsyOZUj2J64LWwmcPsrTX44brNMCcCQkEBQjDw290v+P1BKg/665QnA3+ci44DJh7cS/XMTN5wHIpZZWUMgfYhyaIFAqPuA/GCdx76QCevfFiRsZHkSzyuCloJV+GhrE0sAwhgzhz7Gp6VN3C3KkDKCx18syX+9zSxjTkgdaci/LWdDx6HE6oI5DoMAe/n/hjpsVq7uEbil9m+4hLqHV0gv1fMDP0O8YlxnDzmL7MnZpMiCPIiEVKjA03tuY2f2+VQkfRWvCXxrMNSBZCJKAJnBuAn1jKLANuBF4RQnRBM71l++n+ijaGe/LNTKYMiDU+iw5z8OLPUsl77iECz0g+CRgI5BNZM57hXWcyf0YKK9OPs3DVXtcZ0tB8vGk1za0RWDUvu3Wgf1z1c5Z87+DJtCf5/OSrXD/6Foate47Qr+azrfgxJva/0DZ1kBW1vqNoTfhF45FSVgO/Bj4H9gDvSim/F0I8IoSY6Sr2OVAohEgHVgP3SykL7a+oaG94Ss4J0nCZNrsrRx9dx5Az2ygXoayVgwGY0r83r9w6kqhQB+XOamPtw+yS3dxajTd0QadrJtb3OrMHz+by+MupkTWsDAuHLgOIqDjEksHf1nNW8LY21Vq+t0LhtzgeKeUnUsr+UspEKeVfXcceklIud/0tpZS/lVKmSCkvlFK+7a97K85/rGanaSlxTBkQy8yhPXnquqHMmz6wLpP06XL4fD4AFaN/S5eYMACCAxxGXM3CVZmM7hfNvOkD3RwLGhqgWyvD44YD8NmeoxwZ/ScAxuS9QjSngTrBnZVfamtSO1+/t6JtojIXKJoNXyL8deFS7qxxuQOn89R1Q5mV2tuV7DOfQUfeZ2L+HojsyxIxiuPyfgTw8c5Cjh6sH1fTFmJXnFWa1na45DS/39GV1xMvgayvYO3/UTTpMePZrNuX77bdtkLRGlH78SiaDW9729hlETCb2JamHWJ1Rj7Tk0MZf2gxABnjfsFHxQ8jgsoJklEcOzKIUQlR6Mk0l6Ydst0D53zAanr89qCm2QhRQ0r3TnDZYyACIO0lvli33nCfHtwzAoDtB0rOq++raF8ojUfRbPgSR2IuU1zuBNKZlhJHVKimrVx38mUCDuWzoecQ7tu3hPLqMi7qMowLguby4t6jBAcGsnDVfiOWxbwHztK0+ruNtlaszgBj+8WxqgD6dw/hrslJEOaAi38G3yxhxP5/MHfqU8weG09xuZNVe46zIbPA2PPnfNf2FG0PJXgUzYYvXmPmMks25WpbHXQ9REy4tsFb+H8W80lYKPMcp6itrqXq5BDG9p3HDSP6EdUxwkj4OS0lzohlgTpzmxm7OJrmNst5uqfV9Ni1WzAA3SMD68pP+SPOnUtJLFpHUtk3RIf1Z2naIbLyy1xegEJ5silaJUrwKFoxWuxw+tGTbMgs5KIt/+Z4aDB/6hKDpJbrk2+ii/PHXDeiL9FhDiN/mj6IJ04KN66kuxybtz1oaPfR5qChbRP0zNM/naBZxTcdXcf6M5m8dkjy5OQnqO15M+MOvsj0I/+kqPSHbgGoGpJyZw1Z+aWGQF6ZXj+4VKFoTpTgUbRaZo9NMLaHfvP99zlWupWHukQjBdx98d3MGTLHKFtU5jQW2MFecDS0rbU/Uso0VmvydWvtcmcNZ45fSce4TwjseIRTNfDAugc5se/npEV8QPiJ79i54kUW7kh2ZarWtpbQszzopsfN2YVuz0glDlW0BErwKFothtlNSi4JfIlfdYlGCsFvhv2G2y+83a2s7nxgzU8GdcJA373T07bW/gggbazW5OvW2lqA6ByqQ1J5cc8jABRX5TFsQAgFfR8kfN29TDj0An+67H1+5NL89JQ6vaNC6B0dytypycwc2qOeCVKZ4xTNjfJqUzQ7ntK3eDp+9NtX+X1gATVCcOvAn/Kjfj+rV04PkHzquqH1Zu764Loy/bhPsSxnk17mXAVq6gIosPxit+M7a//C/blB7KztR1DZMW4P+NgwO45PigHgUHEFr319AKBeOh0VWKpoCZTgUTQ7ntyqrceLypw8+/kuHtn6JKcCAxksevKzwXNtz/e04Ro0fnD15vbdEOcyULOozElFVQ0XhFxLp8A4BkReQEBQGfuDnmbZhVdphTb+A04dJTrMwbM3DmPu1CRTBof6+XZVYKmiJVCmNkWzY/XY0tcXZqVqqf71nTWXph1i7/ZH2dAzgJBaSMu6iX+vzeH7I6eYM7GfrSDRhcbm7EJD+2msCa21bR9gzVs3b/pN3Dnpz1TXVvOP7f9gSfoS3i9fzoXJY/nx/k2w+jG4+nmiwxzMHpsAwOh+MW4ZHBSKlkRpPIpmRxcyDy3bbWz/rB8PdQSxcNV+7nt3B5clBHM0drf2Wf5wJiT2J/3oaTZmFbJqz3Hba5u3w26KxqLXozVpAea8dWbNLSggiN+N+B1zh81FIvlzdR5LO3eGb9/gZM43/O3jdH78r41uGa+tqKzVipZAaTyKc4Y3j6mlaYfYmKXliNWj7HWB9NXe46zOyOfOshc5HqrNja5IncsPLkzmna0HOVBYRlZ+GUs25XLvpf3r3eep64baxu2cr3jb6bSozElN0RR61BzlSOC7PBodSUhNNQnv/JbFJb8DBFGhwYZjhfm8Oi3KtZNrE50LlGecorEojUdxzvC2VjIrtTdzpyYxPqkLGzILWJp2iKIyJy+uyWTvsdP0E0dILVhGtSux9M8n9Gdl+nEWr8+hT3So6yramoWeFmfJplyg9WksvuBN8/D0fXQX8gWf7mVS3DV0q/4xUsAjXWJIqvyWWZ20neeLy6tYmX7cLZGofl6Fs5opA2LrCabGcDZrYor2idJ4FOcM80zdOiuODnNw76UD3I4vTTvE4vWa2e0fHd8miBqqRCAgOVVRawRHzhzawwiC1JCW1/OPprg1m13Irx/Zh/0rruJI9edUBJVxWgTw1/B36JX6GrUBwcbz1de/VmfkMy4xhvSjp9mQWcCQXoebnF6nta2JKVo/SvAozhlWTzNPEfr6oDgtJY7C0kpkzgamFKRxWnakWmjCZNk3R10L6wNJjA0nKtVhCCw90PR8HvgaM3jbxSXpQihmYAeclOGM6oujMJO5UZsoGjybpWmHGBEfzZQBscweEw9AclwnFq/Lrpdex5oBoiHUltqKxqIEj6JZ8Dawmmf7f5g+EP59FwC7+90CfEKACOD6EfEEiLodNq0awvk+8PkyeHtal8nKL2XdvnzmTOzHhopOHCotonLcPbD8Xiq+fIy3C4fx97XHDKcLgNUZ+QzpFeHmrKDvYKqCShXnGiV4FM2CeWC1mt3MbtSlaW8SfnQndOrBxbMegKWfECSC6g3MbdG809AivS4QzPsNFZU5mfNaGln5ZTiCAgiJ6QCAs+9YjkQMo8fJbxiZ9zJzp95DhbOWIb0i3bIXmO+jP9+2+GwVrQvlXKCw5WzcbK3netrWWnco0HOKLVq1G+fnD2sXmfoQxytrAaiVAe1iR82GFun1QNjZYxOM765no06MDWP+jBSiOkQB8GH2csJm/h8Aw4+9S1zNcRavzybUEWhsMWFGuVUrmhOl8ShsORtzi6esz+XOaiPpJ9StJeg5xf7UeSXRznwOOJLolPQj/vzFMwBUlsf4vLfM+eza2xRNw7p/UcmRyQQGpvHf9P+SMiGFGUNuQOx6m2uK/8Pp6X/2aEozHwOUqU1xTvGr4BFCXAEsBAKB/0gpH/dQ7hrgPWCElDLNn3VQ+IezMbd4yvZcWOZk4aq9lDtrmD023lggL3dWE+os5NYdywB4sPR6hm/ZwY7S9wBIDr4eX/eWOZ/XJxpa57F+N6uQve/dHaRlRDEw+QYOB73BY18/RurUF4hLX0aHjA+5cfiddA6rn5Xb/LfdMYXC3/jN1CaECASeB6YDKcCNQogUm3KdgLnAFn/dW+F/zsaUZT1Xfx8SrDW3zdkF3PPWt0biznsvHcCdte8SUFXG9o6jSR51JWsK/kONdJIcOp6nZ86iwlnN+KQuDcabTEuJO+u4lNaKNeec1TQ3f0YKUwbEMqHbTKpOD6Ksuoy/7F7M9p4/BaBixQMgZb3fR4+fWr+/gOJyZ5s0YypaF/5c4xkJZEops6WUTuBt4Gqbco8C/wec8eO9FecBs8cmMGVALFtyitmQWVC3hcGJPfDNEmoI4Pcnr2Fl7joOnNlCkOjIv678sxE4uiGzwC0Q0m49YmW6lvVgZbp9Sp3zGatAsAqixNhwXrl1JDeM7MOQDrcRLMLYcHgDSyLjKQuOIe7Ud/D9/+pdV4+f2pBZwGMr0hush1oPUpwt/hQ8PQHzqmie65iBEGIY0FtK+bEf76toQczR8NbByDpARYc5mD8jhXGJMcyZ2I/5M1JYuu0gVR8/CLKW0gtvplf/wcjoZQAMDZ9Ft7BuRpaDuVOT3dYoPGVEaC9p/j1pJivTj7M1q4ZTh68EYHXx63yVcCsAp1f8kaKTp9zKz0rtzc2j+xAfE8rdlyQ3eF/9+d/37g4lfBRNotm82oQQAcDTwH0+lJ0jhEgTQqTl5+ef+8opfMaTh9qc19LqCQP9sxfXZHLrK1vJyi9l+Y7DbMwqJCQ4gJXpx0n74g2CD6yBjpFETP8zY4Z9R2ntUaKCe/H4pb8CMLIc3HtpfwAjg4GdcFFmojpB0jNoPIMiRiACK/gsNp/CsCQ6nTnKvg+fcCsfHeagZ1QouYXlPPfV/gaFiTURq9KAFI3Fn4LnMGAeCXq5jul0Ai4A1gghcoHRwHIhRKr1QlLKxVLKVCllamxsrB+rqDhbrHnR9EEoK7/MMJ3pA9G0lDjmTR9I+tHTrM7IZ85raRSXVQFa/rDKM2U80fkd7cJT/kBO1UkW71oMQGzljZSekbaCzlu25faEtwF/Y1YhBworCCqeRWhQKOuOrGbDmBsBGJn3CpS6T+gak9VbT8Sqa5YqV5uisfhT8GwDkoUQCUIIB3ADsFz/UEp5UkrZRUoZL6WMBzYDM5VX2/mGe1408yCk739j3fHzL1cPJjE2jKz8MrLySwFYty+fM+v+SeSZw1R3Gciison8acPDVNVW0UWOZ3tGDI+tSK83qLUnU5oZOyHjyeRlju15dMYE7h1+LwDPHF7ByaSpBDhPa3v2mLAKE2/3tdJefxNF0/GbO7WUsloI8WvgczR36pellN8LIR4B0qSUy71fQXE+YJcXzZqVwGoKS4wNZ+ldY1myKYcKZy1CCPZn7uOekA9BwtMBt/LvrW/Tsfs3RHeM5pnx81koD2sBkaEOt83h2mteMDs38VmpvVm3L5/VGfks2ZTDvZcOMI7rr9FhDiJCfsgrOz7gyJl0/t7zQv6aHQTbl3Ay5ae8nRfjlrjV+mw9uae3tZRFiubFr3E8UspPgE8sxx7yUHayP++taB58iTXRk3maTWHa34LF67OZM6EfD1V9QsfjZ9geMo4XDkcSkfwKNcCvL/odQ3v25JVb6/xSQh1BLPh0L6GOwHY7wFnjbPQYnsE9I1z7GgmjrPU3en/7YfalTyciKZPlh9cwfeiPGf/NuxS8czePn5pv7NYKGLFV72w9SPrR0/z20v622oxKq6M4G1TmAoXfsNN23NHMcz1Ld9H/+CdUCwf3lFxL9/4fclqUU106gMJjg2Cg+1ntdZCzBojaaRzmvG2e0D6bQFWnKhZ99yyPVObweocuJFbu5ZeRW3g+Y4xhyjRvmwAQHCh45daR9a7ZXjVPhX9QgqedcK5Tyeibkq3OyGfedE1yLFqb5Xa/2WMTCAsWXJ52MwDLw37MiY776Riwmw4inOsT7+O6EX1s69seBzlPZi6zgJ89Nt6n37PcWU1twXgGRH1JRvEenhk4hgU7P+IX1f8lbNIPDcFV7qymwllL76hQsgvKmD+jXgy4UQdv7UnPwQfC5zoq2g8qSWg7wV+eR54Wm82bks1K7V3P+w20WfKcsHV0K93DaUdXlnS6gA5xmmW25OCPCKiJMAYo5SnledG+sZ59evnnvsrm4pA7CQoIYsWpnSxxJBNeXcKEvH8bazyhjiAWr8+mZ1QIr98xisTYcI/X9BTLo09CFq7KZOGq/e36N1TYozSedoKv5ipfU/OD+yzcuqBtuyvo6eOw8hEAliXeRlb5KwhRS0z1FeSWDib96OlG17ct40nTa+yz0bedAMHMoT3Yt2Em35z6H6sSLuZnGdlccORdOPZr6Hahx11jgXpbWegmuaVph9ycS3TNd1xiDKnx0Wfd5hRtDyV42gm+mqsaSrJpN+jZDRy2u4J+MR8qT0LSpZT1j0PsqmJM9wnc2X8e9xXtpF+XMMNzrb2a13yhMc9G/21mj00A4L53d7DpOHTsDmVBYezsfi0XH30HPrmfous+ZMnXueiOCua2UO6sZuGqTMqd1dx76QDD/dosmPRzdM1Xd69viPM5sauiaSjBo3DDF8Fi3dDNuiMm1B8cT6WvpPN371IlOvDfiF9SWbMPgJSY/qQdKCa3sJzcwgMcKi73ecBSNIx1u4PVGfl06xlCGfD9kSJmH53B+tDPiTj4NbtW/JOFO3TPDk1T1R1FdJPp9gMlXicH9TXfhlHabftDrfEo3LBLOeNpvaVuUJNuO2JaN4H7z1fp1Cx3BTE6r+aRjRXszNO8pjoEdmBWam/mTEggPia0wch5lZ6lcejrRPr2E+OTulBwugaAqPBAhiT15aEzPwFgQs6zzJsQxdypyVQ4a137JGnZrGePjWfKgFg2ZBZ4zdHWlJRFKs1R+0MJHkWDWBe5s/JLufWVrYyIj3bbERM0U47mVJDDM19m8IvXt3Nq1VNEnTnI4aA+VKT+krlTk0mM6wiAI1CbOceEdyC3sLwuY7UHlNNB49AH9ZXpx1m4KpPhfSOJDQ8FIDCglmdvvJiC+JmsrRlCYGUJt5xeBEh25pW4riCM6zx13VAjrc7ZJAj1x+RBTUDOb5TgUTSIdUb62Ip0Vmfk89QXGZQ7q1myKdcwu+n2fRAsXJVJQe53/CroQwDuK7+ZzuFh3D21HwdOa6a5DoEdgDrh1pCZTaVnaRr6PkWT+ndlWK8uACR3CyE6zMHgXpH8sfp2nAEd6bB3Gd+tfpctOUWMS4wBcMsubhY+TRX+/pg8qAnI+Y1a41F4xJO3kRbbkU5yXCdj22qA2WPjgTpb/ZnKSn6881E6VFXxTfSVbD6SwuDKcn675rd8ffRrQoJCGNdzHOD7grlyOvCO3W9WVOY0JgsHi8rJLS8gtA+UOI9TU1sDUpInY/mq2x1cceSf/CPsNV4bNoOqoDAWrtrvljHC7FQwLSWuXqyWL3Xyx5qOWhc6v1GCR+GGPkhMS4kzBitw9zbSNxwrKnPy/eGTrpQtsp5QmFTyAQOq9lDm6MLm5Pvg2EHWn/47Rwu/o5OjEy9Me4GEiITm/optGjsPMV0T1RO1ThhwMQc7xpBZksmiXYsIcUwD4E1xFT1qP2ZIVQ6/lm+RNfQhduWV1NvNVf+dF63NYsGne7Wtyx1BTEuJY2X6cUPAeHI+8cfkQU1Azm+U4GlD+CMeQh+49BiN+JhQCsuchieT9V6P/PACY7AxH7+iRzkjsv8JwMrEeYwZ0p2kkgUcde4jukMMl8X8id6hg87+SyvcsNME9L/NgmHfqceZ88UcXtz5Ig+NGMCUAbHcfUky+777GxduvwmxdTG7q0ayOiOEIb0Ou7nG622sLttBjVub0QWRLnB8SetzvsTy+FLP8+W7tCRK8LQh/BEPYR6kQNN4Fq/LJsZlItE7lPVe+mJvubOGZ1dlMCni/wiureQTMYHf7HbQq/o2SqqP0COsB5dE/okXVp4kM2+HsabjqbOqTtw47DSBaMtvFx3mYHTYaOYMmcOiXYt4dOsDnD7wc0bnxnDnjCupCPw1IVsWcmXWIxyavIT12UVsySkyBIo1K7W2BiSpcNYypFckQL08cg39di0Ry9OUtuVLPVVcUsMowdOGMAsNX2zvdpgHrqeuG2rk25qWEmdEpG/OLjRyeJlnwQs+3cucCQk82v1rBhZ/R7GIYJ6cQqd+L1BSXUr/qP78dcw/+HRHOeOTgt2i3n1Nv69oGnYThaCTl9O3w/ccqNxARPwSRg+4BIDXO97I+NoPGHTqAJccfI4nc65xXUXYalR6qp2Fq/YagibUEehT+9NzulU4a70klz03NKVt+bK2pNafGkYJnvMMb7M0q+0dmj5YmyPeo8McLFqbZawTrM7IZ3S/426zaL2TdSzJ5NqixSBg5djbkYffQIoqxvUYx5OTnuSNrzW33rlTk5iQ3MUQkvo6gkq/f26wPsclm3JYuCqTkQk3Ul19AsL3ce+6u3j1ile4ZmQiq8r/zsDtPyPlyPs8fuE4VlRcwMyhPc46jY+1/ep55IB6W2mcazylB/JWB1/WltT6U8MowXMeYO4UvszSPA0CjTEtWO9jt05QbzOwcb2oXnwzgaKKxUnjee7I+yAk1yRfwx9H/5HggGBmpTqMukWHOXjmywy3VCxWVCf2D/Wfo3D9H8zQDnPZUf4E+eRy62e38vIVLzPrqukQ8UdY+TAzcv/GgtN/Y/mOKGaPjbdtQ+bre5v42LUrPY9cc08ufK2zwv8owXMeYO6svswsPQ3WjTEtWO9jvmbipHD7vXe+eoygE7t5uktvXqk5CMCYyJ/x5zH3I4Sodx2Aiqpat1dF8zB7bDy78kpYnZHP3KlJjA1awJbyJ9hZ8A23fnYrT098ka3VV/GzbisIP5bGX4NfYr989qwmPnafRYc5bCcczU1zaNZqvbIOJXho/Q3Cmv/KrsP78h3MWsszX2agZys2e6V52wPHztXaMI9kr4VNz1ErAlkS6gCquLTLvcyffJMhdOzqGRKsxTDrr4pzS1Z+KY+tSGf+jBS3JJ/RYQ5ur3qBO7/8JTvytzPnyzso2H8r6zrdwYtyNzMCt1AWtp7KobOBpk18GvoMGr+Pj7/67tn0K19R65V1KMFDnb3bk7mnpdE7he45ZtcJvG0aZhUmi9ZmGXZ1fdar461jWF2tpwyIZVpKHEu+TONnO+8kAEnFmPuoPfoOAPMm3VhvENGvoXtIzRza081Vt7VPAlorvrr5znktjaz8Mg4WpbH0rrFuv3NocCijQu4nrewPVIRlE9HvP2zKvoPXu97LXQULCPvqj4QljmZW6gCWph1iRHw0T32RweCeEdw1qfG51uzqbF7z8WWr83MxmHszbZ9N+1TrlXUowQPU7VcvvJY6l3iKODcP3N46madGbXfOrNTeFJZWkn70NHdfkszofjGmmIxqyp019eJ2zNeelhLHkF6HAcFH3x5i+Pq5BAQehd6jCZ58L7z1NgEE8frX9QcRa+yH9buoWWHT8NXNNyu/jOiwYLLyy9z20dH5ychkYAFpZ55i+4mtxPVfwo+u/hhWH4ZvXqPmndn8IfwZPttfagSkbswqJMY0OdK1YnMwqa911td8Kpy1HtuhmXMxmHszbZ9N+1TrlXX4TfAIIa4AFgKBwH+klI9bPv8tcAdQDeQDt0kpD/jr/mfD7LHxhvtnYzjXarh19uetkzXG20hPyrkhM4cJyV3cvNP0OA272aY5HkTLxbaf1/t9yfjA3ZQGRfJy9B8oX50BQE1NEBXOGm4e3YeVe07w+ffHmJYSR2JsuDFA2T1zNStsGo1x8zULBSvRYQ5+PSWFM9XPM/rN0ZTWFFJYUcC7jju4MXQjXYqzuKzg71T2f4i7p/Y3NB7dO1EPGtW1Yv2+dv3EU9u899IBxmJ/Q1rPuRjMvZm2Vfv0D0JK2XCphi4iRCCwD7gUyAO2ATdKKdNNZaYAW6SU5UKIXwCTpZTXN3Tt1NRUmZaWdtZ19EZTBYjeOeZNH3jWjd8XjedshZv5evrajnm9ZsqAWObPSKk3UzWfB7Bw1X7iY0L5Rc8srt93HzVScFPVH/ha9qdj9/cIjthJbWUXftrzefafKDUGoCkDYnnl1pFn9R0U5xbzGtDcDTdw4NQBftrjn7y4qpREcZiPO86no6ykbOoCKofdUa9dzJ2aVC99jm7KnjMhgRBHIN7aszeNqTHaVHtHCLFdSpna0vXwhL80npFAppQyG0AI8TZwNWAIHinlalP5zcBNfrp3k/G2kZkvnM3sxypoPEWc2605+dIBG7Kfg2T22AQjKNQcn2PnAqufN3dqEomxYTgLsrmq4mEAXnL8lM3VfYlMeIkaRy5BoiOnjs8kJCGI+TNSKCnfQUGpk7svSW70c1I0D3p7Wb+/gA2ZBUA6HWO0rStKq8oZlRDNlhz4feXPedbxT8K+ms9XhdEs3NwJoF5aHPPOpNsPFAPwRfpxcgvLgToN3tMaIHh2xTZrU2e77tLQ8zjXws16H3NAbYgj0C8TztaIvwRPT8CcnzwPGOWl/O3Ap54+FELMAeYA9OnTxx/1s0VvyN5ySXnTOryp+Q01XG8dzNdzrR2woevPSu3N+v35bMgsBITbNgZmTceKOdZi9th4fjioE45X7ya8+jTOxMvpOfInRH09jypRTPew7vx17DN8s7+jUf8rLujOgk/3si23iGF9o3z+normQ28vN4/uw9GTYdx9STJPfa9tWfH2jm+57aIfEhQg6NbzJioQhGx9jqsy5nFw5L/5uqgTM4f2JDE2HKizBGzOLmRIr0g2ZBYa60HxMaFcNribYX7ztAbozRV7WkqcEcBsrjv4b12wudYarfdxnxz65mBxPtLszgVCiJuAVGCSpzJSysXAYtBMbeeqLlZbrh2N9bIxn+et4XrrYL6ea+2ADV0/OszBszcOc5uN6mWiwxwkTgq3/S5umldNNdGrfwXVB9lb24vnOw1jw5ZfUy2qGdZ1GE9NfoouIV0Y0aPh72pdxFVCqOUwO31k5ZexLbeI4XHD2ZW/i7CeSzlUHc3GrN4cO3WG62+6n94Fe3Bkr2TW/t/zfPEfWb4jinsv7W9cS58UDekVYeyA+ucPd7Mhs5CQ4EBjvdAaPNqQK7ZdO7G2L39MaHyxZpyL+5idK0KasO58vuCvNZ4xwMNSystd7+cBSCkXWMpNA54DJkkpT/hy7eZY4zHjSfVt7DqL3XV8baQNlfX0udk+r88+m4LZlLd8h+a9pn/3Mx/eS8dvX+ZEcBQ/jx9PdvVOAGanzGbu8LkEBwT7fA87jz1/rJcpmo65bUWEBPJk2pO8vud1AEIqJnPiwDSm9O/OxD4OJq69gcSAo6yuuYgtI//Jgz8Y4tUM/LeP01m8PodRCVG8cFOqT33J2tZ9WVfVy4xLjCE1Pvqcmav8ucbrb9rLGs82IFkIkQAcBm4AfmIuIIS4GFgEXOGr0GkJrNqGdZ2lqTmdGtJi7OJtfK2jTt1mX1qMhl7WWldPjgzW2AWzKS/UEcidwZ/Q8duX2RfUkRtie1JVvRNZ04Eru/2G341o3JKdrknqAaj6TM8XF1rFucPa9h4Y+QCJkYn8dfNfqQhZQ9yAHG6b/ASDYxP53PkC3bffwhR20vXQ/4F8nSWbclm4aj/lzhpDAzJwBRJvySnm7je/8UkoWK83LSWOzdmF9fYIMmPWuDZmFfrVXGXuJ03RipR5WcMvgkdKWS2E+DXwOZo79ctSyu+FEI8AaVLK5cATQDiw1BXJflBKOdMf9z8bPO2O6CnDc1Ntv54aaVMcHDxda/6MFA4WpRkxGmAfEOrJdVsXNnrmaXO8zk87rIfP5pPuCOa2nn2o4hTxnZIYHjKXiKAejRYW1u+gZzj2xYVWcW7wNCEpPHoxz03+D49t+yOHSw9w38ZbGNH5Jv48+Xaqkt6m6q0fMvjECvjqMSqc1wJQ4aypd309O0V8TCgbswp9FArS7XVl+nHDCcabadicWd1OMDTVOlAvP2ED7dRuDUfvZw1t896W8dsaj5TyE+ATy7GHTH9P89e9/ImdhuMtw3NTPdk8paDRvcrGJcb4nBbek0aUGBvO0rvGGoNHcbmTzdmFjIiPdhOidt/BPEvUPdu0fVYgsXA1YZseJCc4iNt696Gs9gwj48ZwcYe5VFUHuwSmNDIQ+NKZ7L6DipFoWbxNSOZNH8iiS97gNyv/SGb5RtYVLeaa5Su4O/V2fjDrPwS/OxvWP0lnRxEwjRBHXQokXaDpWSrMGQ/Mv7WddgAwd2pyvW3VG2ojDeWAq9tdN93Wxd+TEK6Xn7AB7NZw9H62ZFNOo/pMW6LdZS7wdf93T5qPecA823Ub3assMTaMjVmFTOwfe9YNUP9OSzblsP1ACRsyC8jOL+NAUblhrvDkuj1/RgrO6t3GjqNL0w7xzer/8cvgJxGils8HXUnZ6V2M7TGWi4Lv5e+fZTFnQgJTBsRSUVXLwlV7KSytZP+J0iatM6nI7ubF3CYB20HVbAL9X1oh326fwdVjx5B2+r+UVB/h0c2P8g9HJx6fdA8T1zzDr50vUxsRxE1jLzXucc9b37Ahs9Bof4vWZtm2d08eXlMGxBpl/NVGNK0+3dDurXgL6DbXx4yn9S1zfXVtTLdytNcsHe1O8FjzsnlqyJ40n8ZuUQDumo1ZxdY7+Ij4aJ77an89u7Wvgs3bHie6G6uGZ0eSojInj61IN0wgIcEBJBRt4JUOTxEkq3ml+nLeK+kDgbsYGjuUy3r2YlvOKRCC1Rn5JHcNZ2R8FO+kHeJkRbWxzmS3htTeZnetFXP7BYw1N8DYI2ll+nH0LBVa2MEgZqX2plPIT1h1YBXzNszjtPM0nwdVMWjiX4ld90fuqVwM3/eHkT9nadohl/s+6O2vocmenXZgl95HpyltKzE23Gswc0NWAb0+5nvr61F1sVD244I592KoI/CsNm48X2mzgsdzY7TPy+apvLkBmgWI9TNP9zdnBtCDNPVGaxZumuvpYTfV21eHBOv6kNlNdebQHryz9SDpR08zc2hPj89L1750T6DEwtVcsWceQaKGb7vP4i85P8RR/BkdukB1TYBha3dWaztHAmzN1QIFQ4IDyMov475367a2VjnYWh+eBlerc4luBtYdAbLyS7nv3XTmzxjPggkLuH/t/Zw8U8br8lIu6lfCJdlPwCe/Y8P+E0y7bK7hHgyCrPxSw1PSijftwJtp62xi4qxanzcHH7v6uAtvTbCmdO9kpKLydr+GzPptmTYreKyaTZ2duYdtjjBr4zU3EPNmUXpHLCx1smRTrkevHLvOO7hHBDOGeM5PZla9deHhzZ7sKQBWt2/r3yHEEcSGzAJWpntekDUPQtGZ/0NumocQNRwffAdPFl/LzWM6sb4wmALg010FPDcjzvhuqfFaUOiwPpF8c7CEWcN7cai4wk3IKq+11od1cPWUqHNjViGp8dHGuos5w/Vvr9bc53OLilnxTSZwMa8Muo8pOU8xfv/f+eBYHlz4G0IcgSxctd+SDd19XdAXIWEX1e/N062hCY9V62tIAFifmfneUaEOr2s25omrnp3dV++4tkabFTxmzcaqqVi3GID69m1Pkf/lzmq2Hyhh8fpswHNQqd4g9ezPulZiNmVY1420xVRJubPGpbZnet0O2E1YuDqu+bq+ZGbQiQ5zcOfEfrDuSVj9GALY2usWni2ZxcbsQhzBgVwcH8aXh2D/8TOsTD9u2MmLy6t47esDzJmQwOWuqHT9GSqvtfMP/bdauEqLhdGQRj8yZ7h+ZX0RBEDemZ0MuiiQ4PJxbOpyLZ/sO82C4P/wo9Nv8ta6o+SOfoQpA2K5+5JkhvSKQO+f5omWdS0I7M2BdUgqnDV8uvsYh4orGNLrSD0XbqvFwirYrEKrIVdtHbO1weyQ461dL9mUy+qMfOJjQo01UXMd7a7fVs1vbVbwmDNOm1PD2KvI1Bvk7WYheofckFlgmKTMDdq8qFjn9hnjZs8dER/NrBc3GesuVtOCPjh7Ehbe4n2swtJsctOvo28AV09Tq6mCFffCt/+lFsEjVT/j1czLmDOhMxLoGLmHVYc+BuCCuL5uz3V8kjY4hTiCvObOao8zu/MJuxgVa7JP/ff+y9UXaCbkfVVcPGwsORVbyHNug6BtlJ78hMikGdy5/15ecDzHjUGrycn8PVcfu41t/WIMbzNzhnK7tSBwd/J5Z+tBRiVE0zcmlCMlZ6ioqmXx+hyjrJ0Lt3lSZ6dtmN2zAVuTtx3WSZ1v6zTa98otLCckOMA4z279t62bptus4DEPylbNwHpMdzuelhJna2IzYx7M9cHbnJvKvP5T7qymsNTJ3z7eY5gG9BljYmyY7QBsV1cz3hqk3exOXxgOdQS6tC6b9D+nj1H19s0EH95CdUAHfnnml3xRO4LxSTHcNTmJR9e8xpf5zyFELUM6XcXD468zBC3UT7PvqY7Ka611YzdxcZ9AaFrK8L7RJMaGu2JlIoAU/jDawfPb3mJ36Qryq/aRH/A0CUNHsrvfiwz6/LcklHzNF+EHcXZ7o56A09uSuV8VlTl5cU0m6UdP85erB7My/bghZEIdgWzILGB430jmTEgwNJ4QR4DHSY958gmi3n47umAbn9TF0Ea8xdpY++kzX+7zHDjrYvbYBCqcNcZ6a2JsuGG+t67/tvVJWpsVPGY85XjS0Wc9Q3odZlfeSY+JN/VrWU1G5sar504zmyt0Qh2Bbm6cdg26obp6a5DmgV3vCHMm9GPu1CTKnTXGjHBUQpQrJ1cpuzd9yg/2/YHgshMck1H8K+Zh4uNHMdclKD89uJSVBQsRAq5LupX5Y+9l8brseoLFvHbU1jtNW8X6u+mCqG632Lr1UevEZtm3oRwoGsKYpLFsL15BeNxaciq2ckf6LgKDr2Yp64ivysH57lV8nvQnFuzsw/r9+QzvG+0asOu0EMDNNK4Heq7acxwhhNvmhdFhDu6anGQIL91ct35/Ac/eeHG9iaaeBko3q5sX+HXBNrxvJFMGxHqNtak/iXIPdLUTgNFhdftgLd9x2Ng+Qq+XefLW1idp7ULwgPtiv+6ZYp316PbaxNgwj3ZeaxCZt9xU5kVaPeFfdJiDV24d6dM21ubOaG68vjXIuo6QllvMxqxCRiVEGR5o/1y1l8htz/CzyrcJELUc6jyMBWEP8ElODfMucDBlSA1T3h9BrawF4MzxK+kafzVCiAYFS1vvNG0Vs1lq0dosRsRHu8VoQd1EQ9fy50xIcHPZv7B7LJOT72Ly4Lk8t+NJ1h5eRVXcav4YmcpvD8cz/MRqfpDxIKVBU3gk82cM7xvFvOkD3RxrQDN59Y4KoUdkiJE5Xfea3JZbVC+WblpKHA8t283GLM1ctyGzwM0F2yxgrGZ1vU/PmdDPsEwA9WJtzNqZvpeV/gruga6etH7zOrHucu1totlWaTeCZ1aqe8QwUG/WU1TmNLxurB5g7q7LdQ3XamYzq+d20dNWF2g7dd7Oy81b0Kq3iG/A6IxbcooJdQTxu+GBXBH6CIOc+0FAWo+buD77cn49NYF5A4P40bDuTP1f3a4W4yLuok+XS9w80nzJN+fpmKLl8fa76IOmPuvXM0xb3a6hLpu1vt3BXZMSjetdGHw3n+V1p2O3/7GrJI0/devNh8MWEPTlw9zIai4PzyIwfjERyWPq7Uir96ekruFEhdpnsjbXdXN2odHOh/WJZEJyrO26i92kyZw30CxcADctz+qpqr++m6ZtKW6NgTLfx/y8zevEdY5E+88q+Pp8o80LHvMPrudv0mcbuseOeZYyf0aK7VYDnjzEzCo8pDcY7Ga+jlmd16hbN/K0PbS3tCb6MXOEtWbak1Q4a8k4UsjAzJdIOvgBHWQlJ0QM1T94nn4DpvJ73QMtsJQ/rP+tcb/B4Zfx2eZ4pgw4zeqMfK8eab7UTdE68GWt0G7fGz3Yef6MlHqTtRjXZMt8nTe3jORgTm/Ck/7OodOHODToCvolTIL3bic6fw+8MR1G3AFj5lHurDZCFJ66bqhhbtP7k10KHHNdndWaxjMhuYuxnmo1m1snTWZtxywEdGFiF/CtPxdzn9cdl8zWCsD4PlYLi6711AXmDjQCTw8WpbH45tQ2vctqmxc8VhObPtvQ09ToiQoBr8n7PC36mxuybsLz5uGiu2/OHNqT2WMT6gWA7sorcdOYrIOCXcyCuW5Z+aWs25fPqIRoIxfVU9cNJfrYRmqyfkdgcCZI2BQ2lbsKb2D4zhieGgA/HRPHh1lLWbxrMUVnigCIC43j2cse4YMux7zu/WNXD2/HFC2Pr2uFutava/Z1ZjUtx5meagnS65mno8McvHzrCB5bkU5JZH8yT+6hzFkGcRfCnNWwZgFs+ids+zcddi4jp/R6PqodjR7fo+9Ya5fRw7ylhl7X535St9eUnSerHeZJ2uqMfLcJob7gbw6ENj+XojInQ3pFMKRXpDFhtDOl6dkJ1u3Lp9BljTBrPboWNy0lzoiRqssl1zYnbG1e8Hgysc0c2sOIoNYb5rp9+UY56+yqoUV/vUxDkchm903d71+P39l+oMRWYzJrbebzo1Idbh5BSzblsjm7kC05RYxKiGLKgFhO7NvG6Zf+QnTRBgIBZ2QiX/T9HT2GTSf2vZ3G/T4qvJcDp7XnM6rbKB4d9yjdw7sb30E333nDU7R3W+w45zuN/V30PmLWeHS8ZYzWU9P8/IsIMk/CiuwVDO4ymIDgELj0EbhwFnw0l7DD23nW8U/+EL6K9SX3cP/2CEPI2WX0MMfz3Htp/wZzMOrmr+U7DrsFoNppd8XlTiCduy9J5rmv9nu0Yuj1mDs1ye3ediEXS9MOGRPdGFP99EmnbknQk/z6MtE7n2nzgsealE9Xba0bnAGkxke77MT26XR82brAUwZb8yIo1Lf9zh6bwOyx7kGX+uee0vTo2txbWw8ae9mPStCyCEyLLuDGyncJ7/ARFEFVYAjPnJnJ/s438+WWk4wryCArv4zxSTFc3KcT/8zNQUrBJdG/5R+X3UyAqMsuDMpk1t4xCyprjjN9lu8tK8XMhOvZenQbb+59k5POk9w79E8s++YYs1IHEH37l/Dtf2H1ArqV7mHW979gYMTFLCicDl1GUlxe5ZaFZFZqb9bt0/ba2X6gyOhDdlnmAdtwB6gLKbBqd7q2NLpfjGGe171AzeYvT2uxdtYR6xqVJ3O6nbbZFmnzggfcMzbrC+66ANEbn651mD1TdOzWdzzt52FeqDSnAfHkTGDtMFYnAqunnTU+af3+fCPwbkJiNH+/OB9n+XP0/X4rAJUymN09ZxFz+QPs/OIYiTHhzJ3alYoqLRXKkN4hPPLVO1pLqAkhIWRsPaGj38v8qlBYF8y9ZaU4ciSB0oOz6dz3TT7O/phv8g6y77trjLWQCucoOl3wLhceeoPhea9xYeW3vOn4lp2n+vHB9z8mmCHoud4eW5HOfZcNgC8y2JBZyJJNOcwemwDYZygwazVVNdr22+OTujAiPppbX9nK3Zcksy23yLa82WXcnO7HbP0APAoPO4rL6+/w295oF4IH6nuuaO7GdWY28+fgnp3XrKXoAqNOC3Hfz8NTLITVmUD34bdqR1ZBZT5vZXqdeU0vLyX0oIAbO27iJ/kbiFmRB0CZ7MBnwdN4ovQKrk8YhSO7gs1Ht7C16Bijkjpwqvoo8UOO8taJPKqDqgDoGBBdT+gqFJ6wC0mw8+KqG8xncrB0EMuO/YWjfEe/C8vYV34Fn20PQTq1rQaG9ZlJVuVwnk7YxpSTH3BReTYXOZ9kbkgEouxG/u+DkazO7sjBonKmDoozLBTmvad085g5HEEXBM/eOMwY9G9/dRu5heV8e7CEkgqtD5i1pWe+zGDhqkzmTEhg7tQkKpy1DOkVWc+JwC7fnNXl2hy8rQswdw1Met0/qK3RbgSPtWNY42PsTFi6NmENRgOMvWt6R4fyzJf7mDm0h6GGm7UW3WMGhLFYqsdF6OaCuVOT6mlAVu8586Lpgk/3EnbmGPEFa/nlwRWM6ZBOABKcUNohjoBRd/L77KF8nHmGkYkBrC56hiNV2wntewaAnaWmByMFNWd608MxlIUzfu7RDVqZ2hRWrCYlsxZgFUp1k6XedFgXyqeFfyHfmU2+81+EJ8Lo0PuoODmIcmc1JwnnlcBrmfqbRyjb+hrFa1+gV1Uu7HiRBbzIzzomsKJoJLHV1zJ3ahIAf/jfd2zJKWLHoRKKy6tsMxSYM5IsWptlmKdLKqo8OCFomkiIIwiQLF6fbdxP79dmJ4I7JyXW2xpBFy5zpyYxZ0IC6UdPM3tMPAB3X5JMVU2ty2LRvrSediN47HI2gb2bpW6P3ZyteXfpwWjmxpsYG87E/rFGwzar4XYbR63OyDfK6HERhWVOl526xLCNe/Keu3NMd8jbxM1n1vLjLp8Q+7Ur2C5QM6etDxzF1ogruP76m0mMi+CRkZXEblzJ5/lPUnrmBABJkUn0j+pPpCOWE4WRdKQbb208A7UhdEuKoaOINQIH739vJ1n5ZcaMTg8m9CWBoqJ94C1XoN1Ebt2+fFLjo+kU1I1jGbcQ2+9dKgJyAdhc/hTlR2/l9mFXEBESrDkvOEIJG38XlUNv5a0vVtA9811Sy9aQInJICc6Bb94hPySBj04P4KKYUXxPb4rLQ0mMDWP+jBRXtujAeoLQnMfQmulax2p6X7Ip1/WJMPq17vzgLrS0wO2qmhrDiUnfPl7PEg+SDZmFDOkVaeS8mzm0h9u923rsW7sRPDq+uFnqnilbcorqeaaY1esR8dGMjI+ql8YDXHvVl1YyPqmL8Zk1pU5RmZP9x0/X95qRtQQUZUJ2OhzbRVXu1wQc+YZAWU0IEAJUBXRkVdWF5MZMonbAdEpqw1i8Ppvqb3bRtcf3fJT1EQdPHzS+06TI3/DItJ/V61ydyOSL9ONsyCw0XDijw4IpKqsiMTYMfdaoC09ve90r2jdmV3/rRE6f+esa/uR+A1idcRdhiU8Q4NDWKEP7vEpR2Anuv+oXhlk5OsxBdHgHTsUMZd7WYEYmzmRgaBq9T+3l6uM76FqRw21BOXDyMx7sGEh2cCKbivuS+fkYLr/sSmYNizfij8wTJ7vgbsBYQ0qO68TiddnGWq1mgpZu31NfezWnv5o5tCcrdh01grVnDu1hZA6ZMyFBm3CWOl1ajjQ8Ao2wh3ZiXfCr4BFCXAEsBAKB/0gpH7d83gF4DRgOFALXSylz/VkHK9bZg95oPOVK07FqHmaz2ebsIrbkFDE+qYuRxuOz3UfZf6LU6HRLNuW45X4KdQQRFeo+QxSVp5gWeYRrh+Qxpeo7Kt/LoWj/Ln5UmU3kmjNGuWCgVgqyg5P46kwynVKmcemVs3jngwxWZxzj54NL2H96JaHx61h64hBoCg5dQrowo98MKosvZvGqCi6MPFTPeSHEEURuYbkp2LRu07rFN6cas8ZpKXEM6XVY7aej8IjVrdoueFv3JAVtz6yMsuvYePIF1xUknx/4mAMF5WxLuxyAmcND2XhkI99VbSAm5Wv2yHL2VAOh8GxCF2KIZXptDL+oKqLz0W9JqtpHUtA+yPwSMh8hXDiYUNOd0+EJXHQqmpzVI+gybBjLs+GqsUOJDu/o9h30yZezutYwdevfo8JZ69oORTJ/RgoHi7SYG11I6sHpehJgXaDoGRXM2dtjwuusG9ZdTb3tL9RW8JvgEUIEAs8DlwJ5wDYhxHIpZbqp2O1AsZQySQhxA/B/wPX+qoMdus1VzxrrLeZAx+yxpgsd3Tw3ZUAsW3I0E1xK904M7xuJkLXszj3GnuxDLKrIZnhcAGfyjjIz4DgXdw0g+dAG9mdnUby7lk6OU5wuyCPUWUiUPMNP9Zvu016OhHTkRwldub6ilvmdL2S7sxfP74ugY+I4jlV24JuDJYw8FcV1EZ2Jin+PiIAvefuI5hkUGAKyNpgBncbw2zE/YVT3UQQFBGkN3XHITRvTv8/cqUnGwunyHUe4+5JknNW1DO4ZQVSoexCr2k9H4Q3zZM0uDMC6rnrvpQOolcmszh7Lf9dVknphBou+f4K9ZasZkQrLC1/gn+8fcLtHz7BeHC7THGikDKRQ5PN6QD4rOkdy/bD5zHD0Y+tX6wkr3EVqUA5xMp+UgANQfoBRQUD6/yAdfgbUfB3ESUcspzt0I7Z7Hzp0juWp2HA+PeVk2qBBhEaWsHLVd+Q7g3ktLZ/4Xp0J6/U2awq70fH7X5Nzei+TBgwxrCG684A+idOdC3RTmy5wrSZKXSgXljp55kttILBqQW0NIaVsuJQvFxJiDPCwlPJy1/t5AFLKBaYyn7vKfC2ECAKOAbHSSyVSU1NlWlpa4ypTUcIzH9/Ky2WZjKyJIaAEOncMpHdUCCApKaskKkSTuSfLK4kICSJQALIWJBw/VU6Fs5qw4ABiwh2UlFZQVeUkNAjCgqG6qooAWaP9owbhSqTp9jx8qGYVgZTIMERIFE5HZ0oCwlgr9hifj+8xmb6dEkg/WkzPaAfr9x8jv6yc3jGCqoBjFFUdAqBXeC/CZTLfZvRkbI+x/OO6UbYeNnqn1+MadDMi1LmX6yY1gHnTB3oMZG2LnUHhG760A72NTRkQ65aU19qmAG59ZSurM/IZO0DwfeA8IzEtQGhQGF2CBlNT1p+MnB48OG0sd05KJCu/lPnLviWiSyZFgZ+zp/h7AIIDHPTpMBrn6ST25HbhiviePD4xhMCi/Wz79muOlmVQHnaKY9WnqKSKsoAAyoWgIkAQLCGstpaI2loia2oJQnI0KIiImlp+UFrGsk5hvN25EwCRNbWUBAYQXwPjqgM5IwQ11RBEABEdOuAIDAQRSC0BnDpTS+ewDgQFBCAQCBEACIQAQQAIwdFTlRw9WYkEunUOoaKqluLyai4eeBd3zbim0b+REGK7lDK10Sc2E/40tfUEDpne5wGjPJWRUlYLIU4CMUCBuZAQYg4wB6BPnz6NrsimvPW8XKbNPrYGFmp3ANCD74OBatffDsC6f5SbIlQOnSyfhwTh30dX4vrnzoYja9jAGgC+PQUEgyMSjtdg1Hlir4k8P/V5bTCIdR8MrEkNwT3ATXfZHpcYY2TmNS+GWtfAVBYCBfgWTGzn8WYNMtUFmO7tOX9GCgGOZXxz/BtqqaV/VH82ft+Rv3+m9eXxSV2M6ybGhjMyPo6Fq05xzyV/YVC3DN7Z/zp02kNWxToIWkd4EmSEdOW32fHsKzhCccfDiBB9jtvB9c833ohwHwRKArVYt9xAyA10dUZDBrsGGvOtqvFOMNDF9N4BhMENsQUeTji/aZXOBVLKxcBi0DSexp4/pvckghBUI7kmYhDRAWEcPVlJ98hQggMDQATgrJUcLKqguhYCAgKJ7xKGIygQhADXLKSyppas/HLyy2o4cqqKIX1jqBWBbMk9RQ0B9IzuRJeIML7OKWZCcix5ReXkFJYTGRpE547BHCwqJ75LGJendCPElQ9O26r3KAeLyhmVEENwoGBDZgG9o0OBak4Efmp8j6ja0Rwr6IyUgSAD6delMwkxnRnbrxeZR+GHQwZyUZzm3mlNaT8r1X0XSXP6Db1sVn4pK3YdYWNWIRP7xxqDSHuKJ1A0Hl+Cia2TFOs+VrNSexumuHnTB5pi4cJJiEgwzus9wsnG/dri/PC+kQCmLeu1oUEImDvuCnp0HEyvrqU8tvo9imszCI/II7/iBPkVrkVPAqg5043+kSl0D+nLyu9PcvVFCYzo251/rjzAkVNliIAzTBwYRq0oJf1YEcWnQomKLCaw8zYqqiu4of/PeCPjVQC6dYwjsKw7xSUBTO7Tlwt6huCsrsRZU4msrUbKGtIPF5ObfwqBRAhJvy6hBAVA5olSEmND6RMdgpSSA4WlHCoqp2OwYFC3TnQMDgAkiT2Hnv0P1grxp+A5DJhbYi/XMbsyeS5TWwSak4FfER078+3sXcZ7O9PAorVZvL/e5a21J9/WBLBobRYffafF1Jg3qbIukurXLi53akn+DpZx49QkdomTrN6bzzUJA7lzZKKxdvT3qXWR0gBLozSvm6e+yGBs1x+xvPAuena8kPRvf2h4zV3UO5K7TG7avxnvfYfSd9MOsfjmVK/pN1amH7dxB63DU3YGRfumIc3X29bnutazZFMuqzPyja0BvDmspMZHkRofXS/L8+yxCW7rRrowO3pwNFMG/IAnfjSEZ9dt4I3t3yFrQhnVcxAj+3Uz+u3QiLoQiUNHq+qZnudMSDC2KegZFUSNrKHSGcyeo6c4Uv01kyN/z7+/rWTKgFgemD4U0Prf9SPqvvczX2bwcXpdXrlLLkw2vod18zddEE+8sP5Y1Nbwp+DZBiQLIRLQBMwNwE8sZZYDs4GvgWuBr7yt75wN5q1zU3p0ZvG6bKAu4aWeMWDm0B62yfjM7tD6NrU6Vo1AT8cDwi2duZ5sUPdOMWectWY7mPXiJrLyy3AExbLlZ1soPQP/jT5gBKhdMrAr0WEOI5q6sLSSmPAO9Tr3tJQ4Fq3LMjLcmu/jLZGiXaev8/DZzcT+sWptR9EgnuLkwF3r0QOkrUkyrZhj4aC+Cc8aR6Q7AM2fkcL72w/TUfakpqySUQlRXNgjyq0u+rlmT9eoUIeRWssc26P1Hc2rc+2WkcybfjOzUnvTxVHXn/R1LXM2Az2VjzVeyC6ZsNXJwBpb1Jbwm+Bxrdn8GvgczZ36ZSnl90KIR4A0KeVy4CXgv0KITKAITTidE5amHTLcmZ3VNYYPv7ljzJs+kKhQ9x9WH5wPF1fw2mbNo8aaqsYa3W9Oh2EOWNPzrOlbXSd3DaeqRrpl9tXrqrtgzp+RQnllAO9v13JEbcgssGgkmttC+tHTbMjMqZeyY2X6cSMGx+4+nhIp2qG7VyfHdWrzcQUK/9BQnJxdmILurl9nQquzIphjgHR344bWlaalxBmTpjkTEpgyINaIywHqCTnd0/VgURpTB8UZ8TtQZ9azyyhi7j/mySzgJoBmj02wFSBWk6UumPXdXtuy96hf13iklJ8An1iOPWT6+wwwy5/39MSs1N4cLi5n3f4CBnbrzGubDzC633HKndWszshnfFKMbUTzkk05LFyV6fKAg77RobbloG4gnzMhgXGJMQzuGVGvkdZ5iaUbws6sPZkbrD7DMXud6VqZrr399tL+RkddmX6cwlInC1dpjXz22IR617I+E/NrQ+gp7YvKnMS4YqA87TOkUIDnfat07NZ+zLP/zdlaRL85C7zuFddQuzVfSxd+IY4gI1uIHjagaxR6mqtpKXHGxm81tceM3ImeMjHYfS993Jg7NYmZQ3uyK6/Edstwb89Cv77dTqttjVbpXOAPosMc9IwKJbewnKuHBhuzFD31RUr3CLe9bPTGWFGluXKO6ReD42AxT1x7kVvgqd0GbOVOLdNzanwUSzblUFxWZZjo9M3ePAVg6g12zoQE21nexP6xrEw/bmhvwYHCMJ8lTgo3/P7NqTzM+8nrnI0rtC/7DCkUQD0twNc2Z27z1q22m7p3UF1exkCjDrqpGurSXJU7q5k6sCs1tcfJLSwn1BFomKILSytZv7+AaSlxDdRBGK/LdxxhdUY+yV07MW/6QCMLti9rpVqWhARXVuy2S5sVPOCueq9MPw5o21yHOgLd9tAwq7dm23PW9jzW7jvBttwiCsucrk5xhHsv7Q/U326hwlnjimzWWJl+3M00oNu39R0WtfppDfaL9ONG0kK7WV5haSU7806SHNfJzR3VbisHu5mSp2y6nrAbNBqrMSnaN41J/WLeN+tsNWqroDJ7e+oTy1EJ0SR37cSQXpEALF6fYzgTmNPq7D9RyobMAh5a5n2dUx9X6tZ7IcQRwJ2TEo04JWsme0+olDnnOdaZunnQhbotad/ZepBRCdEM6NYJEMY+GWZ7bXxMKICx8ZQ5Vsa8nQLAsD6RBAcGUFjmNHbuNG8CZxZ6s8fGGzMvq13cvAf9XZOTuO/dHSxel01IcAChjiBjbcms4TRk/7aa5jx1dLvGb+3QKqBUoXM2ExXzud42WNQdeMyZ4H1pd9ZMHWanBl2z0iej2uSybsdTT+ucuseneS8fve66QwEIisqc3H1JMgeLyo14JW/PrDHP7XymzQoec9S+7p1mtbnqQkk3YwUFCDZmFfLhjsPkFpYzb/pApqXEGe97R4W4Np7KNbQes4vozKE9jAa8cNV+tuQUsf/4acNmrW8CpyEpd2qBZ3YzPXMKjlBXDJAunPTEndatE7xRt8eIu2nO08zKl8bfHmZmCt/wZaLi67l2A7K5P1gzwVvLm98DbumudKGgWynMTgK6o4N1R1HzOqeeDmjOa1qetoNF5WTll9X73uaYJYCs/DLDemKtl/lc/Xx9jamtTurarOCx7qkDWrJO60A9LSWOdfvyGdwzAoCNWYVG0kz9x9dNYHVItxmYvkOhbkLTUpxLth8oMRIOjk/q4uaZY8175mmRscJZawg1/Ti42619Qe+MunA0PwM74eLLoNEeZmYK37BrC75oxHZbxdsJMfOiuzUEwlre/B7qJmx63rNnvtzncgTQNBCzw4zddtTufV3L/5iVX0Z8TChPXHuRmzAxPw/zhBTcLR3melmfmZ5s1LzPT1ujzQoes2lJ2+o2xnBrNEf3r0w/bjgGANw8ui/ZBWVGglC9AekNQZ81mWdgZvVdb1j3Xjqgnoo/Ibku5YfZWcFTBw11BAH2cQ6+LuBadzQ1n2v9u7GoFDoKHbu24ItGbOcQYx60zWZt89qkOSDa3Kf0nYPNx/W/6/qHFjpY4ayup3U0pG1pGox2fq+oEOK7hBHfJazeOVatx06gmetlF5rhbeuW8502K3jMsxdzCvJFa7PcBmGzZ9rCVfuZMiCWDZkFrEzXslfr+3bYBV/qGol1XcjsjaMHhWnl6pLCmbNkA/V2SjRvfW3V0sx1aahze9rR1A61ZqPwJ75oxOYy5vZnHbS9tfOGvC6t5fWMB/qajnmAt3PCsbo4F5c7WbHrKBsyCw3vM7v7Wj1hrSY08+Tx7je/Mfbs0XY5FSqA9HzGPBMz76FuDQIzB7LZZTKwm9GZNRLd7mx1ZfYUFGbXKa0xQOY6mtFdsL/ae5yLekW5CTsrDcUfmFFrNgp/4k0jNpuv9AFW75+64wvUF15WAdWURXmrN6p5gNeFhXU92JytRA/4tmok1vvqbtVmT1jrbqwzh/bgsRXue/boZdsybV7wuCPcXovL3Ru/t7xmYK9pmIM87QQW2AeFeXL5PFxcTnxMKJP6d2VY36h61zLXf0tOMVtyipkyIJbicvvO2BhzmFqzUTQXVvOV1kbr+qenQFPAo2bTmLZuNfGZTV36uqw+obMKSbvJXJ0bdZ0grXA5D+mvejnzbqz6hHVcYgyDe2jrzO1ho8V2JXhmDu3BrrwSENTbLsCc6sbcoDyZtcwNaGL/WBJjw41dF63R/bq5zownW/Jrm7Xtqp/7ar9Hn//ZY+MNT7vIkGC3zAhgr634YkZTazaK5sJuMmaOhWnoXPNrU+9vfrXLNjKxf6yRScQqJK1hBea1ol15JcyfkcL3R04CWjyPmSG9Ikju2snYhkSfsGqWDM1Jqa1niG9XgkdfV9Ejo0fERxu7bdqtl5gblNZJcHPBtEvl4au5ypPnTmFpJelHT9fLs2YmOszBS7eMcIsjsJoIrYJGmdEUrQ1rILOniY+1LTd2gmQ36bJew04QmY95S2Oj56cblxiDEMJty2uzC7deVte06ltYNI1v+4GStq/1SClb9b/hw4dLf1FYWilfXJMpC0srpZRSvrgmU/Z9YIVx7Okv9sqnv8io9/ktL2+RT3+RYZRtzD3OtlxTMX+35rifQtEYrO3TH2U9tXG785vSHzydk3nitLzl5S0y88Rpo0zmidO2Zb3dt7C0Uv7031/Lvg+skE9/keFzvexAS8zc4uO3p3/tSuPxNsuxZiCwc8v0xQzQ0GzMlyhtf2CdwTUmh5byblOcaxpjLmuobEMhA3bnNzaFlPkc8/WLypxGJuzR/dxTZNmtFduND+b+NrxvlCvu8JzsFtNqaFeCxw49LY05QNOTWcq6w6fdOhDgddBuLpOXNwHoqwv2ua6jov3SUHS+XfhCQ+mdGgo9MDsR2KWvakrAdEPbQPiCub+ZN7dry7QbweNLYJh5rQWw3QbAbqakuzfrx7w15MbM9M6V5tFQHZR3m6I58DbBWbIpl4Wr9lPurOHeS/v7nN7J3E/Ma7R6uiq7dVxfLBngeRsD870b02etQtAc3tHWafOCx5sabl00tJrB7Nw29Ybmnv6izg20oUHbl4bVkOmgofMaavQN1aG9NH5Fy+K9r0i3V09lvbV5XRtJjA0z0lXZBXX6M3tHY6wF3sq2dXN3mxY85hmPvlmbnjFa944xuy1as1jbpd7wlM/JnDvtbIVEY7IN2J0HykSmaP146ytWk5OntRFPiTah/rYourlOX4/xFK/XWOwSEpv37fKEN8Hb1vtymxY8ZvtriCOIjVmFbMwqJMZDg7fTZjz96HYBoI2tm6+mg6IyJ898mYE5OE3tl6Noy/gygfN1fSUq1FHPYtGYPtKQ9qH3ZXNcoJ5yy9u55kmsddtvu0lvW6JNCx5zIysud5KWW2TE7NjhSZux4i0Nu3UBU99a17p/iLcOYPVAM8/qvOWuUiYyRXvC09oO2GtDTTVfNaR96Oa7uy9JZkivCMwmd0/ec54C08E+71tbo00LHvNAvDTtEBuzCo1oZF/Ps8NbGnbze30GZJ4JNWWvktUZ+YyMj8IRFGjMhOy24lYo2hMNeW56SgAKjRvUG9KS6hL+xtTLOODJimLNgmK9flvVdHT8IniEENHAO0A8kAtcJ6UstpQZCrwAdEZL0/xXKeU7/ri/J8yzCrsft6kzIOu1PL3q2QQ8JR5tzL10RwNdhTdntzbbqtv6oqRC4Qt22lBTTdENTRIbGlvsrCjW+vlrq5LzBn9EoQJ/Bx50/f0g8H82ZfoDya6/ewBHgciGrn02mQsainhuTPT0uaKpmQ4aE6WtULQ3WjpTR0v3Q9pJ5oKrgcmuv5cAa4AHLAJun+nvI0KIE0AsUOKnOtRD313U7MlmpjUsxvuq/nvL1mumNXwnhaKlMZu79Z1H/YUvVgXVD70T0HARn4iTUh51/X0M8LrwIIQYCTiALA+fzxFCpAkh0vLz85tcKX130cXrso0Nm8zog7evjVL3Pikqcza5TlZmpfa2jbZu6n0a+50UirbIrNTeRpZpu75/NuhCzdt1VT/0js8ajxBiJdDN5qM/mt9IKaUQwmOiISFEd+C/wGwpZa1dGSnlYmAxQGpqapOTFlkDRM+Wc+Fbb6e5tHUffoXiXOMpe7w/UNrM2SM0c+BZXkSIDGCylPKoS7CskVLW21BCCNEZzQz3Nynle75cOzU1VaalpZ11HZtCQ27T3sr6874KhaJ105Q+ey77uRBiu5Qy1a8X9SP+MrUtB2a7/p4NfGgtIIRwAB8Ar/kqdFoaq0rtTX32Rf1WKBRtk6b0//Y8ZvjLueBx4F0hxO3AAeA6ACFEKnCXlPIO17GJQIwQ4hbXebdIKXf4qQ5+pzEqdVNySXlCmdoUitZFQ/24Kea39myy84up7VzSkqY2f6DnfzPvONgQytSmUPgPf/SnpvTjlqS1m9radOaC1kBTZjUq9Y1C4T/8YUFoz9rJuUAJnnOMEiIKRcviD6Gh+rF/8ZdzgUKhULRKzjam5lzE77V3lOBpRagGrlC0Ptqz99m5QpnaWhHKm02haH34c31HOQ5pKMHTilALmApF6xuc/bm+oyaXGm3W1ObNbNVaTVoqv5NC0bZNW3a5GdsjbVbj8TazULMOhaL1YNVw2rLmr7zjNNqs4PHWeNtyw1YozjesE0E1OLd92qzg8dZ4VcNWKFoPaiLY/mizgkehUJwfqIlg+6PNOhcoFAqFonWiBI9CoVAomhUleBQKhULRrCjBo1AoFIpmRQkehUKhUDQrSvAoFAqFollRgkehUCgUzUqr3/paCJEPHGji6V2AAj9Wp62inpPvqGflG+o5+ca5ek59pZSx5+C6fqHVC56zQQiR1pr3HW8tqOfkO+pZ+YZ6Tr7RXp+TMrUpFAqFollRgkehUCgUzUpbFzyLW7oC5wnqOfmOela+oZ6Tb7TL59Sm13gUCoVC0fpo6xqPQqFQKFoZSvAoFAqFollpk4JHCHGFECJDCJEphHiwpevTGhBC5AohvhNC7BBCpLmORQshvhRC7He9RrmOCyHEs67nt0sIMaxla3/uEEK8LIQ4IYTYbTrW6OcihJjtKr9fCDG7Jb7LucTDc3pYCHHY1aZ2CCGuNH02z/WcMoQQl5uOt+m+KYToLYRYLYRIF0J8L4SY6zqu2pQZKWWb+gcEAllAP8AB7ARSWrpeLf0PyAW6WI79HXjQ9feDwP+5/r4S+BQQwGhgS0vX/xw+l4nAMGB3U58LEA1ku16jXH9HtfR3a4bn9DDwO5uyKa5+1wFIcPXHwPbQN4HuwDDX352Afa7nodqU6V9b1HhGAplSymwppRN4G7i6hevUWrkaWOL6ewnwQ9Px16TGZiBSCNG9Bep3zpFSrgOKLIcb+1wuB76UUhZJKYuBL4ErznnlmxEPz8kTVwNvSykrpZQ5QCZav2zzfVNKeVRK+Y3r79PAHqAnqk250RYFT0/gkOl9nutYe0cCXwghtgsh5riOxUkpj7r+PgbEuf5u78+wsc+lPT+vX7tMRC/r5iPUcwJACBEPXAxsQbUpN9qi4FHYM15KOQyYDvxKCDHR/KHU9HvlW29BPRevvAAkAkOBo8BTLVqbVoQQIhx4H/iNlPKU+TPVptqm4DkM9Da97+U61q6RUh52vZ4APkAzexzXTWiu1xOu4u39GTb2ubTL5yWlPC6lrJFS1gL/RmtT0M6fkxAiGE3ovCGl/J/rsGpTJtqi4NkGJAshEoQQDuAGYHkL16lFEUKECSE66X8DlwG70Z6L7i0zG/jQ9fdy4GaXx81o4KTJTNAeaOxz+Ry4TAgR5TI3XeY61qaxrPv9CK1NgfacbhBCdBBCJADJwFbaQd8UQgjgJWCPlPJp00eqTZlpae+Gc/EPzVNkH5oHzR9buj4t/Q/Ni2in69/3+jMBYoBVwH5gJRDtOi6A513P7zsgtaW/wzl8Nm+hmYmq0OzotzfluQC3oS2iZwK3tvT3aqbn9F/Xc9iFNoB2N5X/o+s5ZQDTTcfbdN8ExqOZ0XYBO1z/rlRtyv2fSpmjUCgUimalLZraFAqFQtGKUYJHoVAoFM2KEjwKhUKhaFaU4FEoFApFs6IEj0KhUCiaFSV4FAqFQtGsKMGjUCgUimbl/wG/tJxmejVSigAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "def eval_snr(x, x_est):\n", + " if np.array_equal(x, x_est):\n", + " return 0\n", + " num = np.sqrt(np.sum(np.abs(x) ** 2))\n", + " den = np.sqrt(np.sum(np.abs(x - x_est) ** 2))\n", + " return round(20*np.log10(num/den), 2)\n", + "\n", + "SNR_est = eval_snr(x, best_estimate)\n", + "SNR_data = eval_snr(x, y)\n", + "\n", + "plt.plot(np.real(y), \"o\", markersize=1)\n", + "plt.plot(np.real(x), linewidth=2)\n", + "plt.plot(np.real(best_estimate), linewidth=2)\n", + "plt.legend([\"data\", \"true\", \"fit\"])\n", + "\n", + "plt.title(\"Data SNR: {}dB, Reconstruction SNR: {}dB\".format(SNR_data, SNR_est), fontsize=16)\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.0" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/optimusprimal/__init__.py b/optimusprimal/__init__.py index 454edde..233b238 100644 --- a/optimusprimal/__init__.py +++ b/optimusprimal/__init__.py @@ -5,7 +5,9 @@ from . import linear_operators from . import map_uncertainty from . import primal_dual +from . import forward_backward from . import prox_operators +from . import ai_operators from . import Empty # create logger diff --git a/optimusprimal/ai_operators.py b/optimusprimal/ai_operators.py new file mode 100644 index 0000000..34cc315 --- /dev/null +++ b/optimusprimal/ai_operators.py @@ -0,0 +1,144 @@ +import abc +import numpy as np +import tensorflow as tf + + +class LearntPrior(metaclass=abc.ABCMeta): + """Base abstract class for learnt prior models""" + + @abc.abstractmethod + def __init__(self, tf_model): + """Constructor setting the hyper-parameters and domains of the model. + + Must be implemented by derived class (currently abstract). + + Args: + tf_model (KerasTensor): network pretrained as a prior model + """ + self.model = tf_model + + @abc.abstractmethod + def prox(self, x, gamma): + """Evaluates the l2-ball prox of x + + Args: + + x (np.ndarray): Array to evaluate proximal gradient + gamma (float): weighting of proximal gradient + """ + return self.model(x) + + @classmethod + def fun(self, x): + """Placeholder for loss of functional term + + Args: + + x (np.ndarray): Array to evaluate model loss of + + """ + return 0 + + @classmethod + def dir_op(self, x): + """Evaluates the forward sensing operator + + Args: + + x (np.ndarray): Array to transform + + Returns: + + Forward sensing operator applied to x + """ + return x + + @classmethod + def adj_op(self, x): + """Evaluates the forward adjoint sensing operator + + Args: + + x (np.ndarray): Array to adjoint transform + + Returns: + + Forward adjoint sensing operator applied to x + """ + return x + + +class PnpDenoiser(LearntPrior): + """This class integrates machine learning operators to PNP algorithms""" + + def __init__(self, tf_model, sigma): + """Initialises a pre-trained tensorflow model + + Args: + + tf_model (KerasTensor): network trained as a denoising prior + sigma (float): noise std of observed data + """ + + self.model = tf_model + + # Normalisation specific parameters + self.maxtmp = 0 + self.mintmp = 0 + self.scale_range = 1.0 + sigma / 2.0 + self.scale_shift = (1 - self.scale_range) / 2.0 + + def prox(self, x, gamma=1): + """Applies a keras model as a backward projection step + + Args: + + x (np.ndarray): Array to execute learnt backward denoising step + + Returns: + + Denoising plug & play model applied to input + """ + out = x.numpy() + out = self.__normalise(out) + out = self.__sigma_correction(out) + out = self.model(out) + out = self.__invert_sigma_correction(out) + return self.__invert_normalise(out) + + def __normalise(self, x): + """Maps tensor from [a,b] to [0,1] + + Args: + + x (np.ndarray): Array to normalise + """ + self.maxtmp, self.mintmp = x.max(), x.min() + return (x - self.mintmp) / (self.maxtmp - self.mintmp) + + def __invert_normalise(self, x): + """Maps tensor from [0,1] to [a,b] + + Args: + + x (np.ndarray): Array to invert normalise + """ + return x * (self.maxtmp - self.mintmp) + self.mintmp + + def __sigma_correction(self, x): + """Corrects normalisation [a,b] onto [0,1] for noise + + Args: + + x (np.ndarray): Array to apply sigma shifting + """ + return x * self.scale_range + self.scale_shift + + def __invert_sigma_correction(self, x): + """Invert corrects normalisation [0,1] onto [a,b] for noise + + Args: + + x (np.ndarray): Array to invert sigma shifting + """ + return (x - self.scale_shift) / self.scale_range diff --git a/optimusprimal/forward_backward.py b/optimusprimal/forward_backward.py new file mode 100644 index 0000000..d97e839 --- /dev/null +++ b/optimusprimal/forward_backward.py @@ -0,0 +1,103 @@ +import optimusprimal.Empty as Empty +import logging +import numpy as np +import time + +logger = logging.getLogger("Optimus Primal") + + +def FB(x_init, options=None, g=None, f=None, h=None, alpha=1, tau=1, viewer=None): + """Evaluates the base forward backward optimisation + + Note that currently this only supports real positive semi-definite + fields. + + Args: + + x_init (np.ndarray): First estimate solution + options (dict): Python dictionary of optimisation configuration parameters + g (Grad Class): Unconstrained data-fidelity class + f (Prox Class): Reality constraint + h (Prox/AI Class): Proximal or Learnt regularisation constraint + alpha (float): regularisation paremeter / step-size. + tau (float): custom weighting of proximal operator + viewer (function): Plotting function for real-time viewing (must accept: x, iteration) + """ + if f is None: + f = Empty.EmptyProx() + if g is None: + g = Empty.EmptyGrad() + if h is None: + h = Empty.EmptyProx() + + x = x_init + + if options is None: + options = {"tol": 1e-4, "iter": 500, "update_iter": 100, "record_iters": False} + + # algorithmic parameters + tol = options["tol"] + max_iter = options["iter"] + update_iter = options["update_iter"] + record_iters = options["record_iters"] + + # initialization + x = np.copy(x_init) + + logger.info("Running Base Forward Backward") + timing = np.zeros(max_iter) + criter = np.zeros(max_iter) + + # algorithm loop + for it in range(0, max_iter): + + t = time.time() + # forward step + x_old = np.copy(x) + x = x - alpha * g.grad(x) + x = f.prox(x, tau) + + # backward step + u = h.dir_op(x) + x = x + h.adj_op(h.prox(u, tau) - u) + + # time and criterion + if record_iters: + timing[it] = time.time() - t + criter[it] = f.fun(x) + g.fun(x) + h.fun(h.dir_op(x)) + + if np.allclose(x, 0): + x = x_old + logger.info("[Forward Backward] converged to 0 in %d iterations", it) + break + # stopping rule + if np.linalg.norm(x - x_old) < tol * np.linalg.norm(x_old) and it > 10: + logger.info("[Forward Backward] converged in %d iterations", it) + break + if update_iter >= 0: + if it % update_iter == 0: + logger.info( + "[Forward Backward] %d out of %d iterations, tol = %f", + it, + max_iter, + np.linalg.norm(x - x_old) / np.linalg.norm(x_old), + ) + if viewer is not None: + viewer(x, it) + logger.debug( + "[Forward Backward] %d out of %d iterations, tol = %f", + it, + max_iter, + np.linalg.norm(x - x_old) / np.linalg.norm(x_old), + ) + + criter = criter[0 : it + 1] + timing = np.cumsum(timing[0 : it + 1]) + solution = x + diagnostics = { + "max_iter": it, + "times": timing, + "Obj_vals": criter, + "x": x, + } + return solution, diagnostics diff --git a/optimusprimal/grad_operators.py b/optimusprimal/grad_operators.py index 7e10c79..a986b13 100644 --- a/optimusprimal/grad_operators.py +++ b/optimusprimal/grad_operators.py @@ -1,8 +1,46 @@ +import abc import numpy as np import optimusprimal.linear_operators as linear_operators +import tensorflow as tf -class l2_norm: +class Gradient(metaclass=abc.ABCMeta): + """Base abstract class for gradient classes""" + + @abc.abstractmethod + def __init__(self, data, Phi): + """Constructor setting the hyper-parameters and domains of the gradient. + + Must be implemented by derived class (currently abstract). + + Args: + data (np.ndarray): observed data + Phi (linear operator): sensing operator + """ + + @abc.abstractmethod + def grad(self, x, gamma): + """Evaluates the l2-ball prox of x + + Args: + + x (np.ndarray): Array to evaluate proximal gradient + gamma (float): weighting of proximal gradient + """ + return self.model(x) + + @abc.abstractmethod + def fun(self, x): + """Evaluates the loss of functional term + + Args: + + x (np.ndarray): Array to evaluate model loss of + + """ + + +class l2_norm(Gradient): """This class computes the gradient operator of the l2 norm function. f(x) = ||y - Phi x||^2/2/sigma^2 @@ -64,3 +102,87 @@ def fun(self, x): return np.sum(np.abs(self.data - self.Phi.dir_op(x)) ** 2.0) / ( 2 * self.sigma ** 2 ) + + +class l2_norm_tf(tf.keras.layers.Layer): + """This class computes the gradient operator of the l2 norm function in tensorflow. + + f(x) = ||y - Phi x||^2/2/sigma^2 + + When the input 'x' is a tensor. 'y' is a data tensor, `sigma` is a scalar uncertainty + """ + + def __init__(self, sigma, data, Phi, shape_x, shape_y): + """Initialises the l2_norm_tf class + + Args: + + sigma (double): Noise standard deviation + data (tf.tensor): Observed data + Phi (Linear operator): Sensing operator + + Raises: + + ValueError: Raised when noise std is not positive semi-definite + + """ + + if np.any(sigma <= 0): + raise ValueError("'sigma' must be positive") + self.sigma = sigma + self.data = data + self.beta = 1.0 / sigma ** 2 + self.Phi = Phi + self.input_spec = [ + tf.keras.layers.InputSpec(dtype=tf.float32, shape=shape_x), + tf.keras.layers.InputSpec(dtype=tf.complex64, shape=shape_y), + ] + self.depth = 1 + self.trainable = False + + def grad(self, x): + """Wraps the layer call for gradient of the l2_norm class + + Args: + + x (tf.tensor): Data estimate + + Returns: + + Gradient of the l2_norm expression + + """ + return self.__call__(x)[0] + + def __call__(self, x): + """Computes the gradient of the l2_norm class + + Args: + + x (tf.tensor): Data estimate + + Returns: + + Gradient of the l2_norm expression + + """ + tmp = tf.cast(x, tf.complex64) + tmp = self.Phi.dir_op(tmp) - self.data + return tf.cast(self.Phi.adj_op(tmp), tf.float32) + + def fun(self, x): + """Evaluates the l2_norm class + + Args: + + x (np.ndarray): Data estimate + + Returns: + + Computes the l2_norm loss + + """ + tmp = tf.cast(x, tf.complex64) + return np.sum(np.abs(self.data - self.Phi.dir_op(tmp)) ** 2.0) / ( + 2 * self.sigma ** 2 + ) diff --git a/optimusprimal/linear_operators.py b/optimusprimal/linear_operators.py index 3d26335..a572d12 100644 --- a/optimusprimal/linear_operators.py +++ b/optimusprimal/linear_operators.py @@ -1,3 +1,4 @@ +import abc import numpy as np import pywt import logging @@ -46,12 +47,52 @@ def power_method(op, x_init, tol=1e-3, iters=1000): return val_new, x_new -class identity: +class LinearOperator(metaclass=abc.ABCMeta): + """Base abstract class for general linear operators""" + + @abc.abstractmethod + def __init__(self): + """Constructor setting the hyper-parameters and domains of the operator + + Must be implemented by derived class (currently abstract). + """ + + @abc.abstractmethod + def dir_op(self, x): + """Evaluates the forward sensing operator + + Args: + + x (np.ndarray): Array to transform + + Returns: + + Forward sensing operator applied to x + """ + + @abc.abstractmethod + def adj_op(self, x): + """Evaluates the forward adjoint sensing operator + + Args: + + x (np.ndarray): Array to adjoint transform + + Returns: + + Forward adjoint sensing operator applied to x + """ + + +class identity(LinearOperator): """ Identity linear operator """ + def __init__(self): + """Initialises an identity operator class""" + def dir_op(self, x): """Computes the forward operator of the identity class @@ -71,7 +112,7 @@ def adj_op(self, x): return x -class projection: +class projection(LinearOperator): """ Projection wrapper for linear operator """ @@ -110,7 +151,7 @@ def adj_op(self, x): return z -class sum: +class sum(LinearOperator): """ Sum wrapper for abstract linear operator """ @@ -147,7 +188,7 @@ def adj_op(self, x): return z -class weights: +class weights(LinearOperator): """ weights wrapper for abstract linear operator """ @@ -244,7 +285,7 @@ def adj_op(self, x): return scipy.fft.idctn(x, norm="ortho") -class diag_matrix_operator: +class diag_matrix_operator(LinearOperator): """ Constructs a linear operator for coefficient wise multiplication W * x """ @@ -277,7 +318,7 @@ def adj_op(self, x): return np.conj(self.W) * x -class matrix_operator: +class matrix_operator(LinearOperator): """ Constructs a linear operator for matrix multiplication A * x """ @@ -311,7 +352,7 @@ def adj_op(self, x): return self.A_H @ x -class db_wavelets: +class db_wavelets(LinearOperator): """ Constructs a linear operator for abstract Daubechies Wavelets """ @@ -398,7 +439,7 @@ def adj_op(self, x): ) -class dictionary: +class dictionary(LinearOperator): """ Constructs class to permit sparsity averaging across a collection of wavelet dictionaries """ diff --git a/optimusprimal/prox_operators.py b/optimusprimal/prox_operators.py index fd0c925..b949533 100644 --- a/optimusprimal/prox_operators.py +++ b/optimusprimal/prox_operators.py @@ -1,8 +1,66 @@ +import abc import optimusprimal.linear_operators as linear_operators import numpy as np -class l2_ball: +class ProximalOperator(metaclass=abc.ABCMeta): + """Base abstract class for proximal functionals""" + + @abc.abstractmethod + def __init__(self): + """Constructor setting the hyper-parameters and domains of the proximal operator + + Must be implemented by derived class (currently abstract). + """ + + @abc.abstractmethod + def prox(self, x, gamma): + """Evaluates the l2-ball prox of x + + Args: + + x (np.ndarray): Array to evaluate proximal gradient + gamma (float): weighting of proximal gradient + """ + + @abc.abstractmethod + def fun(self, x): + """Placeholder for loss of functional term + + Args: + + x (np.ndarray): Array to evaluate model loss of + + """ + + @abc.abstractmethod + def dir_op(self, x): + """Evaluates the forward sensing operator + + Args: + + x (np.ndarray): Array to transform + + Returns: + + Forward sensing operator applied to x + """ + + @abc.abstractmethod + def adj_op(self, x): + """Evaluates the forward adjoint sensing operator + + Args: + + x (np.ndarray): Array to adjoint transform + + Returns: + + Forward adjoint sensing operator applied to x + """ + + +class l2_ball(ProximalOperator): """This class computes the proximity operator of the l2 ball. f(x) = (||Phi x - y|| < epsilon) ? 0. : infty @@ -93,7 +151,7 @@ def adj_op(self, x): return self.Phi.adj_op(x) -class l_inf_ball: +class l_inf_ball(ProximalOperator): """This class computes the proximity operator of the l_inf ball. f(x) = (||Phi x - y||_inf < epsilon) ? 0. : infty @@ -182,7 +240,7 @@ def adj_op(self, x): return self.Phi.adj_op(x) -class l1_norm: +class l1_norm(ProximalOperator): """This class computes the proximity operator of the l2 ball. f(x) = ||Psi x||_1 * gamma @@ -270,12 +328,12 @@ def adj_op(self, x): return self.Psi.adj_op(x) -class l2_square_norm: +class l2_square_norm(ProximalOperator): """This class computes the proximity operator of the l2 squared. f(x) = 0.5/sigma^2 * ||Psi x||_2^2 - When the input 'x' is an array. 0.5/sigma^2 is a regularization term. Psi is an operator. + When the input 'x' is an array. 0.5/sigma^2 is a regularisation term. Psi is an operator. """ def __init__(self, sigma, Psi=None): @@ -355,7 +413,7 @@ def adj_op(self, x): return self.Psi.adj_op(x) -class positive_prox: +class positive_prox(ProximalOperator): """This class computes the proximity operator of the indicator function for positivity. @@ -422,7 +480,7 @@ def adj_op(self, x): return x -class real_prox: +class real_prox(ProximalOperator): """This class computes the proximity operator of the indicator function for reality. @@ -489,7 +547,7 @@ def adj_op(self, x): return x -class zero_prox: +class zero_prox(ProximalOperator): """This class computes the proximity operator of the indicator function for zero. f(x) = (0 == x) ? 0. : infty @@ -564,7 +622,7 @@ def adj_op(self, x): return self.op.adj_op(x) -class poisson_loglike_ball: +class poisson_loglike_ball(ProximalOperator): """This class computes the proximity operator of the log of Poisson distribution f(x) = (1^t (x + b) - y^t log(x + b) < epsilon/2.) ? 0. : infty @@ -705,7 +763,7 @@ def adj_op(self, x): return self.Phi.adj_op(x) -class poisson_loglike: +class poisson_loglike(ProximalOperator): """This class computes the proximity operator of the log of Poisson distribution f(x) = 1^t (x + b) - y^t log(x + b) @@ -792,7 +850,7 @@ def adj_op(self, x): return self.Phi.adj_op(x) -class l21_norm: +class l21_norm(ProximalOperator): """This class computes the proximity operator of the l2 ball. f(x) = (||Phi x - y|| < epsilon) ? 0. : infty @@ -881,7 +939,7 @@ def adj_op(self, x): return self.Phi.adj_op(x) -class translate_prox: +class translate_prox(ProximalOperator): """ This class wraps an abstract proximal operator with an arbitrary translation diff --git a/requirements/requirements-core.txt b/requirements/requirements-core.txt index 118250a..4106e6b 100644 --- a/requirements/requirements-core.txt +++ b/requirements/requirements-core.txt @@ -2,6 +2,7 @@ numpy scipy PyWavelets +tensorflow # Formatting pacakges black From 2706e934f61c9406570705cee42cefc5805dec1d Mon Sep 17 00:00:00 2001 From: CosmoMatt Date: Thu, 17 Mar 2022 13:16:25 +0000 Subject: [PATCH 2/2] explicitly pull out theta acceleration in primal dual --- optimusprimal/primal_dual.py | 41 +++++++++++++++++++++++++++--------- 1 file changed, 31 insertions(+), 10 deletions(-) diff --git a/optimusprimal/primal_dual.py b/optimusprimal/primal_dual.py index 347d35f..d0ad249 100644 --- a/optimusprimal/primal_dual.py +++ b/optimusprimal/primal_dual.py @@ -6,7 +6,9 @@ logger = logging.getLogger("Optimus Primal") -def FBPD(x_init, options=None, g=None, f=None, h=None, p=None, r=None, viewer=None): +def FBPD( + x_init, options=None, g=None, f=None, h=None, p=None, r=None, viewer=None, Theta=1 +): """Evaluates the Primal dual forward backward optimization Args: @@ -14,10 +16,11 @@ def FBPD(x_init, options=None, g=None, f=None, h=None, p=None, r=None, viewer=No x_init (np.ndarray): First estimate solution options (dict): Python dictionary of optimisation configuration parameters g (Grad Class): Unconstrained data-fidelity class - f (Prox Class): Constrained data-fidelity class + p (Prox Class): Constrained data-fidelity class h (Prox Class): Proximal regularisation constraint - p (Prox Class): Positivity constraint + f (Prox Class): Positivity constraint r (Prox Class): Reality constraint + Theta (float): Primal-Dual acceleration balancing viewer (function): Plotting function for real-time viewing (must accept: x, iteration) """ if f is None: @@ -34,11 +37,22 @@ def FBPD(x_init, options=None, g=None, f=None, h=None, p=None, r=None, viewer=No y = h.dir_op(x) * 0.0 z = p.dir_op(x) * 0 w = r.dir_op(x) * 0 - return FBPD_warm_start(x_init, y, z, w, options, g, f, h, p, r, viewer) + return FBPD_warm_start(x_init, y, z, w, options, g, f, h, p, r, viewer, Theta) def FBPD_warm_start( - x_init, y, z, w, options=None, g=None, f=None, h=None, p=None, r=None, viewer=None + x_init, + y, + z, + w, + options=None, + g=None, + f=None, + h=None, + p=None, + r=None, + viewer=None, + Theta=1, ): """Evaluates the Primal dual forward backward optimization with warm-start @@ -50,10 +64,11 @@ def FBPD_warm_start( w (np.ndarray): First simulation from `r class' options (dict): Python dictionary of optimisation configuration parameters g (Grad Class): Unconstrained data-fidelity class - f (Prox Class): Constrained data-fidelity class + p (Prox Class): Constrained data-fidelity class h (Prox Class): Proximal regularisation constraint - p (Prox Class): Positivity constraint + f (Prox Class): Positivity constraint r (Prox Class): Reality constraint + Theta (float): Primal-Dual acceleration balancing viewer (function): Plotting function for real-time viewing (must accept: x, iteration) """ # default inputs @@ -85,11 +100,13 @@ def FBPD_warm_start( max_iter = options["iter"] update_iter = options["update_iter"] record_iters = options["record_iters"] + # step-sizes tau = 1 / (g.beta + 2) sigmah = (1 / tau - g.beta / 2) / (h.beta + p.beta + r.beta) sigmap = (1 / tau - g.beta / 2) / (h.beta + p.beta + r.beta) sigmar = (1 / tau - g.beta / 2) / (h.beta + p.beta + r.beta) + # initialization x = np.copy(x_init) @@ -105,14 +122,18 @@ def FBPD_warm_start( x_old = np.copy(x) x = x - tau * (g.grad(x) + h.adj_op(y) + p.adj_op(z) + r.adj_op(w)) x = f.prox(x, tau) + + # Primal-Dual acceleration step + x_accel = x + Theta * (x - x_old) + # dual forward-backward step - y = y + sigmah * h.dir_op(2 * x - x_old) + y = y + sigmah * h.dir_op(x_accel) y = y - sigmah * h.prox(y / sigmah, 1.0 / sigmah) - z = z + sigmap * p.dir_op(2 * x - x_old) + z = z + sigmap * p.dir_op(x_accel) z = z - sigmap * p.prox(z / sigmap, 1.0 / sigmap) - w = w + sigmar * r.dir_op(2 * x - x_old) + w = w + sigmar * r.dir_op(x_accel) w = w - sigmar * r.prox(w / sigmar, 1.0 / sigmar) # time and criterion if record_iters: