{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 05 - FWI with total variation (TV) minimization as constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "There is a lot of research on regularization to improve the quality of the final result beyond the simple box constraints implemented in the previous tutorials. In this tutorial we look at how more advanced FWI techniques such as [total variation denoising]( https://en.wikipedia.org/wiki/Total_variation_denoising) applied as a constraint can be implemented using [Devito](http://www.opesci.org/devito-public) and [scipy.optimize.minimize](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) using [scikit-image](http://scikit-image.org/). This is a variant of the FWI with TV constrains algorithm described in [Peters and Herrmann 2017](https://doi.org/10.1190/tle36010094.1).\n",
    "\n",
    "[Dask](https://dask.pydata.org/en/latest/#dask) is also used here to speed up the examples.\n",
    "\n",
    "This tutorial uses the same synthetic datasets and model setup as the previous two tutorials, so check back if you get lost on parts of the code specific to Devito, SciPy.optimize or Dask. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up (synthetic) data\n",
    "We are going to set up the same synthetic test case as for the previous tutorial (refer back for details)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#NBVAL_IGNORE_OUTPUT\n",
    "from examples.seismic import Model, demo_model\n",
    "\n",
    "# Define the grid parameters\n",
    "def get_grid():\n",
    "    shape = (101, 101)    # Number of grid point (nx, nz)\n",
    "    spacing = (10., 10.)  # Grid spacing in m. The domain size is now 1km by 1km\n",
    "    origin = (0., 0.)     # Need origin to define relative source and receiver locations\n",
    "\n",
    "    return shape, spacing, origin\n",
    "\n",
    "# Define the test phantom; in this case we are using a simple circle\n",
    "# so we can easily see what is going on.\n",
    "def get_true_model():\n",
    "    shape, spacing, origin = get_grid()\n",
    "    return demo_model('circle-isotropic', vp=3.0, vp_background=2.5, \n",
    "                      origin=origin, shape=shape, spacing=spacing, nbpml=40)\n",
    "\n",
    "# The initial guess for the subsurface model.\n",
    "def get_initial_model():\n",
    "    shape, spacing, origin = get_grid()\n",
    "\n",
    "    return demo_model('circle-isotropic', vp=2.5, vp_background=2.5, \n",
    "                      origin=origin, shape=shape, spacing=spacing, nbpml=40)\n",
    "\n",
    "\n",
    "from examples.seismic.acoustic import AcousticWaveSolver\n",
    "from examples.seismic import TimeAxis, RickerSource, Receiver\n",
    "\n",
    "# This is used by the worker to get the current model. \n",
    "def get_current_model(param):\n",
    "    \"\"\" Returns the current model.\n",
    "    \"\"\"\n",
    "    model = get_initial_model()\n",
    "    model.m.data[:] = np.reshape(np.load(param['model']), model.m.data.shape)\n",
    "    return model\n",
    "\n",
    "# Inversion crime alert! Here the worker is creating the 'observed' data\n",
    "# using the real model. For a real case the worker would be reading\n",
    "# seismic data from disk.\n",
    "def get_data(param):\n",
    "    \"\"\" Returns source and receiver data for a single shot labeled 'shot_id'.\n",
    "    \"\"\"\n",
    "    true_model = get_true_model()\n",
    "    dt = true_model.critical_dt  # Time step from model grid spacing\n",
    "\n",
    "    # Set up source data and geometry.\n",
    "    src = RickerSource(name='src', grid=true_model.grid, f0=param['f0'],\n",
    "                       time_range=TimeAxis(start=param['t0'], stop=param['tn'], step=dt))\n",
    "    src.coordinates.data[0, :] = [30, param['shot_id']*1000./(param['nshots']-1)]\n",
    "\n",
    "    # Set up receiver data and geometry.\n",
    "    nreceivers = 101  # Number of receiver locations per shot\n",
    "    rec = Receiver(name='rec', grid=true_model.grid, npoint=nreceivers,\n",
    "                   time_range=src.time_range)\n",
    "    rec.coordinates.data[:, 1] = np.linspace(0, true_model.domain_size[0], num=nreceivers)\n",
    "    rec.coordinates.data[:, 0] = 980. # 20m from the right end\n",
    "\n",
    "    # Set up solver - using model_in so that we have the same dt,\n",
    "    # otherwise we should use pandas to resample the time series data. \n",
    "    solver = AcousticWaveSolver(true_model, src, rec, space_order=4)\n",
    "\n",
    "    # Generate synthetic receiver data from true model\n",
    "    true_d, _, _ = solver.forward(src=src, m=true_model.m)\n",
    "\n",
    "    return src, true_d, solver\n",
    "\n",
    "# Define a type to store the functional and gradient.\n",
    "class fg_pair:\n",
    "    def __init__(self, f, g):\n",
    "        self.f = f\n",
    "        self.g = g\n",
    "    \n",
    "    def __add__(self, other):\n",
    "        f = self.f + other.f\n",
    "        g = self.g + other.g\n",
    "        \n",
    "        return fg_pair(f, g)\n",
    "    \n",
    "    def __radd__(self, other):\n",
    "        if other == 0:\n",
    "            return self\n",
    "        else:\n",
    "            return self.__add__(other)\n",
    "\n",
    "from devito import Function\n",
    "\n",
    "# Create FWI gradient kernel for a single shot\n",
    "def fwi_gradient_i(param):\n",
    "    from devito import clear_cache\n",
    "\n",
    "    # Need to clear the workers cache.\n",
    "    clear_cache()\n",
    "\n",
    "    # Communicating the model via a file.\n",
    "    model0 = get_current_model(param)\n",
    "    src, rec, solver = get_data(param)\n",
    "    \n",
    "    # Create symbols to hold the gradient and the misfit between\n",
    "    # the 'measured' and simulated data.\n",
    "    grad = Function(name=\"grad\", grid=model0.grid)\n",
    "    residual = Receiver(name='rec', grid=model0.grid, time_range=rec.time_range,\n",
    "                        coordinates=rec.coordinates.data)\n",
    "    \n",
    "    # Compute simulated data and full forward wavefield u0\n",
    "    d, u0, _ = solver.forward(src=src, m=model0.m, save=True)\n",
    "        \n",
    "    # Compute the data misfit (residual) and objective function  \n",
    "    residual.data[:] = d.data[:] - rec.data[:]\n",
    "    f = .5*np.linalg.norm(residual.data.flatten())**2\n",
    "    \n",
    "    # Compute gradient using the adjoint-state method. Note, this\n",
    "    # backpropagates the data misfit through the model.\n",
    "    solver.gradient(rec=residual, u=u0, m=model0.m, grad=grad)\n",
    "    \n",
    "    # Copying here to avoid a (probably overzealous) destructor deleting\n",
    "    # the gradient before Dask has had a chance to communicate it.\n",
    "    g = np.array(grad.data[:])\n",
    "    \n",
    "    # return the objective functional and gradient.\n",
    "    return fg_pair(f, g)\n",
    "\n",
    "import numpy as np\n",
    "from distributed import LocalCluster, Client\n",
    "\n",
    "# Dumps the model to disk; workers will pick this up when they need it.\n",
    "def dump_model(param, model):\n",
    "    np.save(param['model'], model.astype(np.float32))\n",
    "\n",
    "def fwi_gradient(model, param):\n",
    "    # Dump a copy of the current model for the workers\n",
    "    # to pick up when they are ready.\n",
    "    param['model'] = \"model_0.npy\"\n",
    "    dump_model(param, model)\n",
    "\n",
    "    # Define work list\n",
    "    work = [dict(param) for i in range(param['nshots'])]\n",
    "    for i in  range(param['nshots']):\n",
    "        work[i]['shot_id'] = i\n",
    "        \n",
    "    # Distribute worklist to workers.\n",
    "    fgi = client.map(fwi_gradient_i, work)\n",
    "    \n",
    "    # Perform data reduction.\n",
    "    fg = client.submit(sum, fgi).result()\n",
    "    \n",
    "    # L-BFGS in scipy expects a flat array in 64-bit floats.\n",
    "    return fg.f, fg.g.flatten().astype(np.float64)\n",
    "\n",
    "# Start Dask cluster\n",
    "cluster = LocalCluster(n_workers=5, death_timeout=600)\n",
    "client = Client(cluster)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FWI with L-BFGS-B\n",
    "Equipped with a function to calculate the functional and gradient, we are finally ready to define the optimization function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import optimize\n",
    "from skimage.restoration import denoise_tv_chambolle as denoise\n",
    "\n",
    "# Define bounding box constraints on the solution.\n",
    "def apply_box_constraint(m):\n",
    "    # Maximum possible 'realistic' velocity is 3.5 km/sec\n",
    "    # Minimum possible 'realistic' velocity is 2 km/sec\n",
    "    return np.clip(m, 1/3.5**2, 1/2**2)\n",
    "\n",
    "# Many optimization methods in scipy.optimize.minimize accept a callback\n",
    "# function that can operate on the solution after every iteration. Here\n",
    "# we use this to apply constraints and to monitor the true relative\n",
    "# solution error.\n",
    "relative_error = []\n",
    "\n",
    "def fwi_tv_callbacks(x):\n",
    "    # Apply boundary constraint\n",
    "    x.data[:] = denoise(x.reshape(181, 181), weight=5.0e-3).flatten()\n",
    "    x.data[:] = apply_box_constraint(x)\n",
    "    \n",
    "    # Calculate true relative error\n",
    "    true_x = get_true_model().m.data.flatten()\n",
    "    relative_error.append(np.linalg.norm((x-true_x)/true_x))\n",
    "\n",
    "def fwi(model, param, ftol=1e-6, maxiter=20):\n",
    "    result = optimize.minimize(fwi_gradient,\n",
    "                               model.m.data.flatten().astype(np.float64),\n",
    "                               args=(param, ), method='L-BFGS-B', jac=True,\n",
    "                               callback=fwi_tv_callbacks,\n",
    "                               options={'ftol':ftol,\n",
    "                                        'maxiter':maxiter,\n",
    "                                        'disp':True})\n",
    "\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now apply our FWI function and have a look at the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      fun: 199.0333099294088\n",
      " hess_inv: <32761x32761 LbfgsInvHessProduct with dtype=float64>\n",
      "      jac: array([-7.24996760e-15, -4.09159536e-14, -1.42357399e-13, ...,\n",
      "       -1.36085370e-13, -3.84566307e-14, -6.74617866e-15])\n",
      "  message: b'STOP: TOTAL NO. of ITERATIONS EXCEEDS LIMIT'\n",
      "     nfev: 8\n",
      "      nit: 6\n",
      "   status: 1\n",
      "  success: False\n",
      "        x: array([0.16, 0.16, 0.16, ..., 0.16, 0.16, 0.16])\n"
     ]
    }
   ],
   "source": [
    "#NBVAL_SKIP\n",
    "\n",
    "# Change to the WARNING log level to reduce log output\n",
    "# as compared to the default DEBUG\n",
    "from devito import configuration\n",
    "configuration['log_level'] = 'WARNING'\n",
    "\n",
    "# Set up inversion parameters.\n",
    "param = {'t0': 0.,\n",
    "         'tn': 1000.,              # Simulation lasts 1 second (1000 ms)\n",
    "         'f0': 0.010,              # Source peak frequency is 10Hz (0.010 kHz)\n",
    "         'nshots': 8}             # Number of shots to create gradient from\n",
    "\n",
    "model0 = get_initial_model()\n",
    "\n",
    "# Apply FWI with TV.\n",
    "result = fwi(model0, param, ftol=1e-6, maxiter=5)\n",
    "\n",
    "# Print out results of optimizer.\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUYAAAEKCAYAAABuTfznAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztvXu8ZUV17/sdvXfTD7qhoZsGwcZGReNbDGL8YI74PMiNGnO8R4waJR5RY3KCGnN85KpX8/BtuPGB7QtJUKOCSgyKmCCIERQQBWlNCCAQWpoGgQb6we4e94+qWqtW7aqaNdeaa++9uuv7+ezPfNWcs+Zca9f61RijRomqUqlUKpU+i+a7ApVKpbLQqA1jpVKpBNSGsVKpVAJqw1ipVCoBtWGsVCqVgNowViqVSsBYG0YROV5EfiEi14rImyPHl4jIP9rjl4rI+nHWp1KpVEoYW8MoIlPAR4HnAI8EXiwijwyKvRL4tao+FPgw8N5x1adSqVRKGadiPAa4VlWvU9WdwBeB5wdlng98zq5/BXiGiMgY61SpVCqNTI/x2ocBN3nbNwNPSpVR1RkRuQtYDWzxC4nIycDJACze9zdZ/Rv9mrumfbd3QnjMocFyd+TY7mAZ7vcHCuWO5fbH6lLCXAxSqj9L3RG+S2m5P3Ys/E775yyK7Gu6fup/4tbLt6jqQQAPFdH7IqfG2ATnqerxhcUXLONsGHMfQ5syqOoGYAOAPOBo5RWXwRp7cIVd3uOd4I4tDS40Eyy3e8e2t1zGrpO6/gyzCcuG+2PkjnVFV9+IcX6zShj3uyp5vrDMdIv94T63XJrYnzuWui/M/q665Qfkl67IfcCrI6fGeGf/v2+iGefX92Zgnbf9QOCWRJmbRWQa2B+4I3tVwdQ6/LD9j6PpCzjKP41/ze3BsXA7vP8w9y05Z5jrpj75kmuVfGvmohFvQ1f1SX2HhvlPSn0//fWmpf/jn2o8Q4HgU/DD7P7l9ibGaWP8EXCkiBwhIvsAJwLnBGXOAV5u118I/KvWrBaVyoJiEbCs8G9PYWw/BNZm+MfAecAU8BlV/ZmIvAu4TFXPAT4N/L2IXItRiicW3yDXTSh9qtg5uV9XGFSFqbLbg+OhkozRRtF0oX5GUYZdqcq2jEuFtlF/qbIzwXZqX4zc8SalmFOMoXKMET5P5LsqwOLMJfZExqqQVfVc4Nxg39u99e3A/z3OOlQqldHYG7vSk/u8Xdh4cufnbDmO7Q1lw+NtGEUdtTm3jZooOWeYOnRBl/fLXWsc/zElzpdQBcYUY8q2mKtz2LuJlK2KsVKpVAKqYqxUKpWAqhgnAcV0dcI4QL/708VTNYX8xMqkutA5J0xTFzBmzE/RdTjQMM6XheJ06aJrnYv7yzldwvNL4wrbxCTmutK5WEe/zk11sDiv9N5Eza5TqVSyOMVY8td4LZF1InKBiGwUkZ+JyJ9GyuwvIv8kIj+xZU7yjr1cRP7D/r08PLcrJlcxhpH6JSEPpftzZXIhPk1KMXa/ppCemPLpIlC9hGGcL3PldCm9zyiOqK57IeG1ck6SlDJckdgf29fm+93wnjp8/Bngjap6hYisBC4XkfNV9RqvzOuAa1T1uSJyEPALETkT8/TvAI7GtASXi8g5qvrr7qpnqIqxUqlk6VIxquomVb3Crm8FNmJyJgwUA1bahDIrMDHOM8B/B85X1TtsY3g+MJZx2ZOnGHdj1FVOMTYFabehRDGG26X7c/fJMdc2t2GGNo7yzZrrYPfwWiW2RRLbsWPDDO9rCs8pCddpoxgzZcfllbb5V48CLg0OfQQzKu4WYCXwIlXdLSKxxDRho9oJk9cwViqVOaWl82WNiFzmbW+wSWAGEJEVwFnAKap6d3D4vwNXAk8HHgKcLyLfozDpTBdMZsPYZGNs84sZntOmXBeKsWnYWM7GmCvTJW0856XndHn/ru5X4mluOreNd7pE/TXZFnOKMUWbbE+0DtfZoqpH5wqIyGJMo3imqp4dKXIS8B6bN+FaEbke+A2MQjzOK/dA4LvlVSun2hgrlUoWP6FV01/jtYzd8NPARlX9UKLYjcAzbPmDgYcD12HyLjxbRA4QkQOAZ9t9nTN5ijGMYwyHNEH6128YL3R4za4VY5v7hcdS1xq3jW6hKdZhy47T0xzbN4y9sM05TT2lIZV3xwHexwIvA64SkSvtvrcChwOo6mnAu4HTReQqe/v/o6pbAETk3ZjMXQDvUtV8msIhmbyGsVKpzCldOl9U9WIa8sSr6i0YNRg79hngMx1VJ8nkNYwlI1/CfSVxjKV2tFhsW3ifpv2pOsTOiR0rqVsTc5WkogtGea6SmMsS22KTRzl2v7be6GEUY4ltvY1SHL9inAgmr2GsVCpzirD3DQmsDWOlUskiwOLSlmKhTWkxJJPbMOYmtkpNguXIdZWGSeqQOjfXRWvqxpWE6YTXKjmWKruQv9Bdd/vH6XSJDTQYJVi7qQs96mCG3P+RRQSma8NYqVQqfURg8dR812JumdyGsU24zjBPWaLOSp0sbRw3wzhfRqFUUc4HXdaljVovuU6JM6bU2RJTf00pxGJKcZQwqsxUv60U4x7CXva4lUqlLSKweMl812JumbyGMQzXceTCdYZRBG2ON6UXKwnXaWNj7CKJRBuluNBSibWhjVovuU64nVNy47AxNgVxl+C/i1ApxtLf7YVzG+xlj1upVFpTG8YJIFSM7gn8X7phkh6UBlrHbJmlCWr9c1MzDObqmKpTl17brpXjuJJJlCr7uU4IEbMxNtkUR5mmoI0tOjcYIkzKEjJ5LcVI7GWPW6lUWiNA9UpPALEkElu84yu8cjlyT9/mzYQewpin3N/vX79JZcbq1BQvmaNNHOM4YyBHVZIpJdjGu9+kvkYZ3uevNynEkuF9Oe93KTHF6L5v90SOOWpXulKpVAIEqF7pCSWWfmkYSm2NTXWAtHKM0WRzzB0bJvYxlTQgZn8q3S491rbsKHGoOe9+G7tdqVJsM/IlpSRjZbpUijHFmLMxVsVYqVQqAbVhrFQqlQjV+TJh5IZIdXn9NqS61LFuVhM5h02XQeGZfHxF3e7Uvi7CdUrMGU2OqRLaDO8rGdaX6maXhOKMswvtf6fCLnR1vgB73eNWKpXW1IZxDyFnbPa3Y7/Qqe0YoUrJJbYISTlb2qiWplCVUUm9x5y6HCcl92vzuYXbJaE3owRpd5GNOxfc3/R9H0UxduSVFpF1wBnAIZhZ4jeo6qlBmTcBL7Gb08AjgINU9Q4RuQHYCuwCZppmJByWPbNhrFQq3dGtYpwB3qiqV4jISuByETlfVa9xBVT1/cD7AUTkucDrg0mvnuYmxxoXk9cwuiGBjtwwp2GSSZQONcsRBmuPYv+MDYNrsjX6dWyyQ+ZU4EhKMTUPeu7k1MvNzJ2UulwXITgxBTeORBA5+2Rqu4SUUswpxhgdNoyqugnYZNe3ishG4DDgmsQpLwa+0M3dy6nzSlcqlTxuSGDJH6wRkcu8v5OTlxVZDxwFXJo4vhw4HjjL263At0Xk8ty1R2WsilFEjgdOxbyyT6nqe4LjbwD+F+Y37TbgD1X1l61uUmI3WRqUjSmfLn6hQ3LD/IYZYpYi54lta3csUYyz8NWhK3R/5sJNuEovDrazs262v3ybd5+yLbYJ0h4mlViX3uhcgHfuM26nGLeU2P1EZAWmwTtFVe9OFHsu8P2gG32sqt4iImuB80Xk56p6UXHtChmbYhSRKeCjwHOARwIvFpFHBsV+DBytqo8FvgK8b1z1qVQqQ+KcLyV/JZcTWYxpFM9U1bMzRU8k6EbbOadR1c3AV4Fjyh+knHEqxmOAa1X1OgAR+SLwfDxbgqpe4JW/BHhp8dVLhrKlvMS5+LvUGxnlTQ0Ta5nznKeUzT0002RrjO2bdUyDA746DPeFJ4dKEmbPWpwynvrlgnPaDN0sVYqxuNM2nuVS2+IwSrHk80rZFmO9q5yg79DGKCICfBrYqKofypTbH3gqXpsgIvsCi6xtcl/g2cC7uqnZIONsGA8DbvK2bwaelCn/SuCbsQPWlmDsCcsO76h6lUqliG690scCLwOuEpEr7b63AocDqOppdt8LgG+r6r3euQcDXzVtK9PA51X1W53VzGOcDWPMIBR1VYrIS4GjMb8Qs09S3QBsAJBVR2vjiIuUUgxtfTEbYxvlEZ7b5vgwtsWU9zk8J+fJLh0BEy3rPr5tQYFtXpn7g2XyYh7u/FAZ5j6MQrtjm3jQYTzMKxL7c+eUJIhIfaeGUfolXulxJPKIoKoXU2AsVtXTgdODfdcBj+umJnnG2TDeDKzzth8I3BIWEpFnAm8DnqqqO8ZYn0qlMgw1UW2n/Ag4UkSOAP4LY0j9fb+AiBwFfAI43hpTu6E0Vi+mjsahHEPF6pdNqb2YighVZJOCTO2LURSjmFKK90UuFNoYY7ZFR+h9Ds9ZHqmcW0/YGkNy76RNTGJKIXYRxzjqf2OTQoxFcJQEC9Qhgd2hqjMi8sfAeZjfm8+o6s9E5F3AZap6Dia6fQXwZWs3uFFVnzeuOlUqlSGoiWq7RVXPBc4N9r3dW3/mOO9fqVQ6oCrGCaKk+9Em0LXNsMFSct3yYbJ9h9cp6X43Ba6XGPNndYtdPO624Li/r6QLHd60jR0jfNiEPT/3TroY3ldyTmlAeexxh3GWjeJ8iX3/asNYqVQqAbVh3EMoTXYwzPC3XOTIKMInpRT9eqScL6FyXBopM0zA8KzwnKalv54K8I4RhuksC5a5QPLEpUocUk2qbwV9QqdLygnjn5MKGA8/x5LvR5dhOtmQrATVK12pVCoeVTHuoZTYGFPndBGuE/vVTwWfp1RhjPA5SgK8m64V3elUYGhbDLf99SbFGKvQ8si+1Dnh9cNhhcEpMftuF+E6uQDvlGJssvtCuZIrGQqb2l9K9UpXKpVKQFWME0D4IbWxteSUYsorXTIqLUXueCqhRalS8MvmPOqlXuneff1Rm2EAd6gUtwbH/TJNXmlf4aUCu0P8yjt16coujheN2fVKbYvD2BhjirGLQO42kRWppBHD2hhrw1ipVCoBdUjgBJFTWqEASf0qtlFnJbRRkCmV18b+0+SF9K/bqBTDAzDb6xwqxRIbY3jd6aAczPZGh/uJHL8vsi9CG690ShXmjrWZ0GoYO2/p8FZotinm/h9ydauKsVKpVAKE0eYtmkAmv2GM2U1S8WIlNsZwu+QNDWM7alK1w4zoiU2hkDo3tT2g5EKFmPJGlyhGgv0JL/IA4Yv17+OUYsnIGsoUY0oN5o4VKcbUxGCWmYIpG0p6B01xiyW9kdh3p3alK5VKJaB2pSuVSiXCXtZSTO7jlnQpRgnbGeXNtAngTYUF5bo9qefIdeearjFr+B/MzsYddpNj2bpHmR3QkbIz5O5j6x92S0cJ8I45X5qcLn63ebrh2WemB88ZpUsN43W+dNSVFpF1wBnAIcBuYIOqnhqUeRPwEq9mjwAOUtU7mmYe7Yo6r3SlUsnjutIlf83MAG9U1UcAvwW8Lpw9VFXfr6qPV9XHA28BLrSNYsnMo50wuYrR4Z4gNnjfUSJaRhE2TWEZwySeyNHkdBnGcRNVZTnFNnDyiLjrhhVPzR9TQJvEHSXOl1S4jlN7S2fPyrFoetfA9u6ZhOyayVS2yenSxvkSo6QF6HBIoKpuAjbZ9a0ishEzcd41iVNeTH8K1caZR7uiKsZKpZKnnWJcIyKXeX8nJy8rsh44Crg0cXw5cDxmDmqIzzx62JBPlWXyFKNSbiNpExybKjMKs5RjJGyjpxQD+1JulsCmBLVj+1RdiM22YH9J/FPueKrCixPLcB16iWpTar0k8DpnY5ylJq2KtarQqcOpQCXm2O1WhlGMORtjSjGWKMdUuE75d2qLqh7dVEhEVmAavFNU9e5EsecC31fVO7yahDTEQg3H5DWMlUplbuk4XEdEFmMaxTNV9exM0RPpd6OhcObRLpj8hrHkCdoMwG9SmSW2q1Rgb85L2VOOdmW7/XGMJZ0NbUm5tFZNdtDwvtHkDqmHjh0P7YElSSRSyjAnhd16MCQw5S3OpQNrmsVvoIx5nkVLdwJ9hTjV5IH22JWyNfo0fQ9LpilooxSb6M4rLcCngY2q+qFMuf0x88y/1NvdOPNoV0x+w1ipVMZLt4rxWOBlwFUicqXd91bgcABVPc3uewHwbVW9152Ymnm0s5p57JkNY6ltMacY25zT9BZjaiK0RTkV4couDZSjf5/SJZRP3nSPu89K7+QwGWwbw6srm0ry4FcynMpgv2Dp6nSgd85qu7Tpx4YZ3tekHAfOseo/sCU6pTidsS3O2M82VJVJL/XAycGyZI7oJm+0/+pLPtJuvdIXk5y5bKDc6cDpkf2zZh4dB3tmw1ipVLqjDgmcAFJe6RhtvNJhLGAbO00pOY9l6ti0Z4tLTXoVqqOY3cktVwXXCLnHn14gVHAH26VL+RW7SFOSiNhXLmVbDNORRa7pnnlVsFwTHG8ziiWadNY866KEUizxRs+UKMRe4ZZLGM3G2DTyZfJaCgBEZF9gu6qWhwswsY9bqVTmjAlqGEVkEcYp8xLgicAOYImI3Ibpgm9Q1f9ouk4N8K5UKo3oVNnfAuAC4CGYoYSHqOo6VV0L/DZwCfAeEXlp7gIwMb8DEXKOj6YhgSXOl6audWxfqTPGw3XNkob4gaQEQQhP2IVu03Vy57qu5yF2OWDMt/ebsV3p6f0YIOfsyZVJkfqcYtsph1PYTQ672tDvZq8KyszqSvvvfrDLnOpC+9tFYTlNNA0JLAk5GxFdBDsnJ1HtM1V1VoyYDRI/CzjLxlFmmdyGsVKpzAkqMDNV2rnc3VxkjLhGUUQeAtysqjtE5DjgscAZqnpnrOEMmbyGsc2QwKanKxnJlttO/TJ3oRxjxnwXwuOiHXqJDIJyMQUXOmZS6iIX6tM070nsnGEC8EOH0T3BdqxMSrXHArybVGXvufqJIRb3ArrbyzGnIhudL8OEj5Wk2xsRFWHXdOmXeWc3Nx2ds4CjReShmIDyc4DPAyeUnDx5DWOlUplzdk0tDANiC3bbgPAXAH+rqn8nIj8uPXlyG8aS0WLhdu4XdJjwiNS5YWKIIZSjI9oxccrR2RzDoYE55ZFMmGDsaYtXufldYKVdXzlllqu4E4DlNlxnpZ0LZpk3r/Rym2BiH4zaWpJQELu8MWYzdn2bDdbewT4AbLWB3W7/Vi/43K1v3WWXd5rl/XbZC4zPBTvPeheDKcScSgRYYvc1heXE7Irhvt52b+lsuV6hUQYcdIwiA5/XhHC/iLwYeDkmGQWUTTQETHLDWKlU5gRFej9eE8RJwGuAv1LV6+346n8oPXnyG8aSxAkhw/yylijFpHKMVCShPHJeTkeoQJxdMlSbAPtYpbN8hVFyofpzyzVsAWA1t/fOdeur+LUtc3v03OWeYnQqch+rFKfsy5hisG6+Atlpx5vdFyjDrdbYdycH2GXftXy7HRJ455TZt2W1cTXfvnr1QFm39NXmfbvMfXZu32egTs4GmAvaDvcl1WDkurvs96BnR3bfizZREiX2w6690kjvM1roiMgG4JvAd1T1f7v9qno9UDwNQrZhFJEHYoIlfxs4FJOI72rgn4Fvqur8uqAqlcrYmbCu9GcwyW3fICI7gW8D31LVn7S5SLJhFJHPYrLjfgN4L7AZY415mL3x20Tkzap6UeYaRRPXiMgLgS8DT1TVy7I19rMJE1nmaIpvLDk2TB6FgrqVpK8KExW4c/ZZYtRZaPuDvqpbbRXhwWwGYC23DmwfatPauf3+voN3mDL7bra/g66ISx96l1cpt+4cuqlJvXwvsRMj+9rl/nbp8kTY3BHbD+2fcvu+B9qqmGGKt/AAADbbbbd/M2sHtgFun7Kqct9BNbnDVmSntXH63cddwYe4a5dLDGE+gxKleL9TqKFtcZgUYm1sjB0oyElpGFX1Ekwg9ztFZDXwbOCNIvJY4ApMI/mlpuvk/mU/qKpXR/ZfDZwtIvtgUwXF8CaueRYmweSPROQcVb0mKLcS+N8k0ptXKpX5ZUJtjKjq7ZhEt18AEJHfxIi6RpINY6JR9I/vBK7NFCmduObdwPuAPyupcI82SnEURvFkj+Cdjtm5nDJcYuXYskAhOhvgWqsCAQ628m6dnSrDqUC37ZaHbbby70avEm52DbfP5Up2l3fmSD8xvVOMTvG46ofvxP8/c+rRKUaXXWz14HLpA/qnHHaoqe9hh5vlEw7fCMAd65baqq8bWG6iLzdvsetueatVlaEt8z4vZZqzsTmPuau/U44OP1YxqRS3W+doaiqC2L42URIUHGuB6UpPljtCRFYBfwCsx/vP8+2OORqfVkR+B9N4PciWF3N93S97YnzimicF1z4KWKeq3xCRZMNoJ9Qxk+osTYrUSqUyBozzZZ/mgguLczFd6qsYYjhOyc/A3wK/B1ylqm0mnslOXGOzYHwYeEXThVR1A7ABQFYdrQM2xhwpG1/s3KZf4DZjpUu807PshYNprPaxMXT7TPVHX7gYQacQnd3QKcVQDfrr67lhYHnQTXY4icszcp1dXu9VKlSKbmnF5e1WMd7hPcq2YBmmuo1NoOB0mfMb72eF1erAxuiZCekJQDcDiP29PPAII7EOfLB5sIce8Z8A3LBkfe/U6zHra+3JTjlu6ilIc6Pbe5K1b4ecsp5zpyZztrdZ8YqhUgxH9OQUY6zMOIh8/81gs4nrSi9V1TcMe3JJ83ITcHXLRhGaJ65ZCTwa+K6ZBoJDgHNE5HmNDphKpTKHdNeVFpF1wBmY//fdmDRgp0bKHYcRZYsxMw8+1e6/AdiKMdLMZGYk/HsReRXGedxTF96Mg1lKnvbPgXNF5MLgBsmJbCzZiWtU9S76eU4Qke8Cf1YbxUplYdFxuM4M8EZVvcI6Xi8XkfN9p6y1D34MOF5VbxSRtcE1nqaqWxrusxN4P/A2+j1VBR5cUsmShvGvMKJ/KZQbGlIT14jIu4DLVPWc0mtFCWse62KkUiWVhNG0CY9o26X2CNNYhV1oP/TmgF7ozWDX+QF2eXjQbfbX1991MwCLXdc5XIaOFujpe7VL559xPhf307vVOyXM7V0yR6B7LS53+DLbNT9w8+Dy4E3e+WH3PrHc91ZjXnrUkc5WAKsOHQxQX2n7tP3hjINOLhgM94F+F3rXlKn9riA0B7xA7qYudC5JxjBOlzHQVcOoqpuATXZ9q4hsxPgjfKfs7wNnq+qNttzmWRdq5g3AQwsa0CglDeOBqvrsYS4em7hGVd+eKHvcMPeoVCrjpaViXCMifq9vg/URzEJE1gNHMTtU72HAYtuLXAmcqqpn9KoD3xYRBT6RujbwM/CGZLWkpGH8jog8W1W/PexNxkL4i3mPtx4mUyhQbsX3G0UxzqQnR+sFawdK0alE6IfhOKXolqFjZb3nQXnQ5tvMSqgQXREnpJxi9KzA99vf6VttCI5TiG7ponS29U9JOl1yuI/DneucMfcH+7d5geSH2QsvD50SOxqWwGE7zBPsc4RThoMKsT+MMV1711D0wndiname8yWoY4ljZZhwnVwKudjxQhTpBb8XsCVj9+shIiswqcFOUdW7g8PTwG8Cz8B8HX4gIpeo6r8Dx6rqLbZ7fb6I/DwxyGQXcKWIXMCgCbCbcB3gdcCfi8gOzHe1NFynUqnsAXQ9JNBm0D4LOFNVz44UuRnTwN4L3CsiFwGPA/5d1Rh2VHWziHwVEy8daxi/Zv+GorFhVNWV4T6xbuQFxcCMbiNcp8mWUzLQPzcdQoAL03GKcXZITj+pgwvWTgVpr7MGwsPvuK1/g+uDZVOwtqfK7rDr7uc8XG4LtmG2bbGNYmxjPpu207A7y9/y1BDRqWDb4yDbzdh1hHl/7p/fhab4nth+gPfgsMHe9rRZDgwNdPbGlDJM2Rpz54zLtpi5bpcNo207Pg1szDhwvw58RESmMVL8ScCH7Yx/i6xtcl/McL93Ja5xtapeHtz7uYmys2jMV26dJf72Ilqk76lUKpPPDFNFfwUcC7wMeLqIXGn/ThCR14jIawBUdSPwLeCnwA8xeRauxvwOXiwiP7H7/1lVv5W4zydF5DFuw+Zm/IvS5y3RVoeLyFtU9W9EZAkm2cMVpTfoHDe1QajORrEbhuv+dpgENhdgG5YpkELOG718ytiJXSKIVYEHGvo2xrWBcuwt7/0VAOJHi7r1MPFDqBCtArvv3v6pznLtvM6h/dBt+49VamP0P66U5zoMBvc92a5Oy219l9nnEDes0D1fmKDCX7fLQ/Y1J+9Ya16WS3+2jf4c2y6g+74gwLsX+D01ewjn7jAB7TA2xpKeS1v8l19w3S6HBKrqxcQHf4Tl3o8Jt/H3XYfpUpfwQuArIvIS4CmY4YHFTuSSpz0JOFNE3gI8DZNu7MOlN6hUKpPNhKUdA0wjKiInYuyMNwHPVtVtDaf1yKUde4K3eSrwCeD7wIUi8gRVnR/VmJoMy6cpQW3MO930y5lLVBuqyJKpFCw9b7T1jIbeaJckFvrJZN0QwNXBcqlTg37UV5joIVCI4XKb9ywp9RfaEX3FF5ZNqUF/v1OCqVceeqv9fU7VbrN+R6cge8/p1KE/3sGlNdtvcHvN/uZl3blkdnJbt37PrCS6pgZ+zOMsUrbonP0wLDtP8YvgvNKTMVZaRK7CG3qMGVA6BVwqIqjqY0uuk007Fmz/Gnik3a/A08urW6lUJpUJy67zO11cJJd27Gld3KBzSqZPnZW8IbGdu054LLQ1xva1GKWwKJFk1qUSW8HgFATQV5Ph1AIr77XuTKeSfO9mqAzDuD5bN7XLGS8hRJNC7FrEuOtOJ7ZnImV7qtJ+FstT8Yue7bT3ngIVve9dZpTMqrWD7xf67959LqnRMn6i4Z4qbvpexGyMw4x4adFTacsEdaVvV9V7cgVEZEVTmaRXWkRemgvLEZGHiMhTmutZqVQmGWdjLPlbAHxdRD4oIv/NhvQAICIPFpFXish5FCSrzenj1ZjI8cuBy4HbMNroocBTgS3Am0d5gkqlsvCZJOeLqj5DRE4AXg0cKyIHYDT0LzBzVb1cVX/VdJ1cV/pUEfkIxpZ4LPBYjM17I/AyN8B7Xpih312IJYpo6kI64g6/AAAgAElEQVTnuhxN3ZDcOU3dHu9tO6eLM9q7eZhd1yxcgt+N22qPmW730oQjBZg990pozLddZ9eFTjlLFgKxuvW61GGm8FwoTNi9Dpbh+x3ct23g2PLA+RLOyzNAqgtdMrwvda0YTebA3LmJcJ0WQwLnnViOhrZkX6Gq7gLOt3+VSmUvZJIUY1dMjKuph3O+5FRfavB8SjnGyjZdK7avhaG8NxQQpxwHExks6SnIvmpxZfYJlrNUoS9aShNcLCBKqpYq45xIPeP4rmAJaRVtFeOSHdYhtqT/7t27dp/HPj2l7+bPTs9F3Vjpku9U0zXHTG0YK5VKxWNSZwkchdowNtGkIIekNwsgg8tQDU55UmefVFqswF7Yyla1pxK+i9x7CETe1IwJ25la0j/Qf+c2sXDvcxvcn1WMjnEM8xsjExbHCICIfAD4rKr+bJjzS2YJXAL8D2ZPQ5jKalGpVPYwJrAr/XNgg83Q81ngC3Y6lSJKfga+jgmBvRxy454qlcoA025kWhAOPFniayKnT1XVTwGfEpGHY/I9/FREvg98UlUvaDq/5CN6oKo2BkRWKpU9k0m1MYrIFPAb9m8L8BPgDSLyalU9MXduScP4byLyGFW9avSqdoAwWOuStGMl80qnrjdKWrPMOTNu8iTKltBPlursPT27Ty4p64SpE5+Rqh6+i+nMsWB717QZEOYnTui/83gy297nNdOiAcl9LxeQvXFCbYwfAp4L/Cvw16r6Q3vovSLyi6bzc9l1XJaKaeAkEbkO05V2UxsUZamoVCqTzwTaGK8G/kJVYxNiHdN0cu5noJMsFWPDjXiJKYKUysttN/0gltxnCHUZKsMdvZT5g0t/fWew7A1KcO/E/w5PB8e6UMJjpqRKsx7DPrOkpjaIvZPw3dj3eN8Sk4x2pzfaw71rl6h2ZzDFQbThCCc/a/q+zETKEjnWlhHV54QGeL9EVT/j7xCRf1HVZ5Q4YZJJJFT1l6r6S+Av3bq/b/R6VyqVScDZGLuY2kBE1onIBSKyUUR+JiJ/mih3nJ324GcicqG3/3gR+YWIXCsis3I1iMhSETkQM43rASJyoP1bDxxa+swlP86PCm48hZnasFKp7AUYr3RnY6VngDeq6hUishK4XETOV9VrXAERWQV8DDheVW+0U6W6tuejwLMwMwn+SETO8c/FJI84BdMI+sm077bnFpGzMb4FeCuwTETuph9zsBNITXI9N0zTrktY4nxp6uaUdKVbdFd32RnkZnpdaPPF2xbMN+JnkQ73uW7d9n1NeuqlwVwmQL+bHXa3g66m64r686o0vdq56oXH5nyZdSz1zt3zZuZ8CZfh/C7gv/PBYzuC2QMHnC8uXGda4nUr+Q636QY3DfcsuVakDl12pVV1E7DJrm8VkY3AYYDfuP0+cLZLVKOqLif9McC1du4XROSLwPP9c1X1VOBUEfkTVf27YeuZy67zN8DfiMjfqOpbhr1BpVKZfFo0jGtE5DJve4OqRoWU7d4eBVwaHHoYsFhEvgusBE5V1TMwDehNXrmbMVOr+td8uqr+K/BfIvJ74T0T81jPouRH/632Bk/BeKm/p6pDT2Q9Mi5cZxSlmDu3za/60mA5hDPGhUGESsQplJWeYrwTMxeJyyzd297XbC/d3yYl3p8+bt3NeeLUUaAgxS6nve//4mB0W6jcXEI0//FGSVsWKsJw27/PsmDftFsJZwV0yxXeyW6fm/NltVncvdrc0b1Xt4T0nC/bbE2cE2YmFq4TqtdUxnefcKbJEsYU4tMyjnGLqh7dVEhEVgBnAaeo6t3B4WmMue4ZmI/6ByJyCfHZBTXYfiomRCc2h7QCnTWMH8Ukp/2C3X6NiDxLVV9XcoNKpTLZdB3HKCKLMY3imQkFdzOmgb0XuFdELsJMm3ozsM4r90D6EwSbuqq+wy5PGqWOJU/7VODRqqoAIvI5YP6CvZ1iDBPU5lI3tYn/CLdTNquSMjmlaJXFzh029GaJkTpOOTq1ssxLO9ZPVGskjkuB5ZKorjzQKMalfjCCU4oH2mU4e15gg1zmPd+0m7PZbocKMTa7X6loyYn1UA26bV9Bun3OCrg8ZTcMZwSEnkLkYLPQtWa5ZcocuJ01drm6d8oWu/5rBmcQdAo/q6hS3482qd/CmSjnMAC8yyGBdrqUTwMbVfVDiWJfBz5ixznvg+kufxgz/vlIETkC+C/gRIw9Mnafvwbep6p32u0DME6fvyipZzJcx+MXwOHe9jrgpyUXr1Qqk0+X4TqY2QBeBjzdhuNcKSIniMhrROQ1AKq6EfgWpp35IfApVb1aVWeAPwbOw8wk8KVM9pznuEbRXvPXwAmlz1yipVYDG0XEDal5IqbPf4694fNKb9YJoY3R4c/5VRq83LWNsdHWGJpD+spxxxLzi+yUSD8p7Zpe2X4y23jC2n32NTk+jji0P6WFpKY9SCyXe97b/axK2WZtjaFydNu+XTFUe20coaENc1mw7Ftb+/sOdIrQLVcHS6sGByLYDh1c3movstlKyFvsgVudpPSOOTV5T6AYo11NN2PgtH2iptkkfcLeTvhC27zg1LVb0FVXWlUvJm4rDMu9H3h/ZH/ptAVTIrJEVXcAiMgyKI85Knnat5derFKp7HlM6MiXfwD+RUQ+i3G6/CHwudKTGxtGVb1QRB4EHKmq37Et77Sqbh22xmMhNynWMDbG1DVycYxD2BqdF3PnLhvHOGXk2dZIotowKWpIL1nqgf3jDzriNnsjBpfBvNIx9XKgXd92R7IIMKgYw2OLI2X8/TDbluiW+wXLA+lzsFW2i93OUBk6E70zAj3YO/kIs/jVOqMUb2A9ANfb5S08YGAJfXvjnYGNsTckcFfMG20/h6X2abuYVqLNucPcZ8xxjHOFqr5PRH4KPNPuereqnld6fmOTISKvAk7GfC8fgvEEnYZxpVcqlb2ASUw7BvwY8xusdr2YEi31OkzE+aUAqvofbojOvBLaWmKKcRhGGfnSuEz/dPdsjVNGOW4dsKjZMkGqqzCBgVv6HsSda8311i+5GYDF4UgYt4w8n0vI8EC7XLw5XtRXf/cFx9wTL2MQ/zWGNsVQIbrlYZ5klFAhHh4srSrkyGAJ/HLtQQDcYAvdZOXlDT3FeKhd+orR2BadYnQjXXINxiKrGHc72/JSa1obZTRLm15QGzI2y90s6nJI4JwgIv8TY6P8Lsam+Xci8iZV/UrJ+SWvd4eq7jRedrAu9NlehHjljgdOxQw8+5SqvifxAO+01/yJqkbd75VKZf6YtK408DbgiW44oYgcBHwH6KxhvFBE3JjpZwF/BPxT00klA75F5EjgLcCxqvrrBaFEK5XKAJNoYwQWeWOsAW6nLDwRKGsY3wy8EhPU/WqMq/xTBec1DvgGXgV81MYYETxIHDevdEnwdpfdjWG60gW4ZBK7Zky3a5cdk3fflAsD6X8hwxyNYY5A1/32h7K59dv3N86D9UffAMAhB9tI7yB0ZSCs5cbB5cGbB5e3326Wd3i+IDe2y4X0pHqNseF9vS60feTVLuTGdaFjITepLrRd3nGksbHc5A2YCLvQYdfZhebEhgT28zE6M8bghz3tzRK4y3WllzpPl+uONkarlM8HHvs+dhz8bf7lJq5h/JaInEd/xN6LKAvzAcq80rtF5GvA11T1thYVaxzwjRksjp2kZgp4p6p+K7yQiJyMcQDB0sPDw5VKZaxM3tQGqvomEfkfmIBywSSz+Grp+cmntUN33oGJNBe7axfwd4VTp5YM+J7GmMWPw3i7vycij/Yj1gFsdo4NALLf0TqgGGO/jqlfzmE+29RQQX+9cenST6XnHO4lH9hulIibn3iXl9VhZso5VwZTlPVUoQ0p8QOTnRpy6siFpKxbd9PA8vDHmN88uZE+bt39vLlRqVYxrraKcbU/BNFJRjeELTWvpG/LD1ODOYUYBmn7hhb3+xiE5fxq7f62yusGlrd4ctOtb7YXvDVQiH0Hi+fE6s23MzjnS4g/r/SUc7bZ5+uHLCWUY2xspVsOk1RiGPaQcB0AVT0LMya7NbnXfAqmtX2iql4PICIPBj4uIq9X1Q83XLtxwLctc4mq3g9cbyepORL4UYtnqFQqY0SRgR+JhYyIbCXuHHZzVe0XOTaLXMP4B8CzVHWL26Gq14nIS4FvYwZ15/gRzQO+vwa8GDhdRNZgutbXlVS86Jd0lF/XnFIM940QphMSpq3yE586u9Y2q0q2Thu715YpE0rikkr46sglmFiNkXcHcysAa63s620faJaHHripd+7Bj3dlzfLAzfZlO0uwVYy9RBXQH2IYKsZcuIkTUOHwvkA53n1oPzDIJXyYPYxv7cB+pwY3e3JzSxB64xJ3bAsSQozafZwOewgp5TgTCeMJv99Lg+1xhe1EmKRZAlV1dpzbEOS8NIv9RtG78W3EkymH5aIDvkXkXSLixlefB9wuItcAFwBvUtXb41esVCrzxS6miv4WEiLyFBE5ya6vsSKtiNzPwM4hj/WIDfhW1bd76wq8wf61Y5hhTsOkdh/JxlgU7gn0vdPOLhVLeNr3YA8e2x1sLxqwc5n15StsirIlRkG6ZLfhcjX938I1VhL2jq01x1atNdsH2P3Lej7ovkJdYqVib5hiMIwx5m0Pp24Ivey+l3hLkBrs9iAtWP+cAwauBbDtXqMQe3N7B+/esWTp7K/5VMZOXHIc+nbj3T0bZDBk0F9vUoqx72O4vRfOEigi7wCOBh4OfBaTvuwfMObBRnIN4+PsXC+z7kl340wqlcoCRxF27Z6shhF4AWbahCsAVPUWO/lWEcmGUVUn403kvNJtSP3alsQoNpWxKmJRiZqYmQ62+x9DTxn2PNhWaTg1YZe7vaSmu+37uWvG2Jzvmj4EgJvdbVaEy/7JK1YZ9bdiX5cgN1yaXG/LvWS6bn2fIAnGVPBB+TarnUF8plOOYXymr/p6++4yFd9+pz12j30nLg1dLHohkezj/hVu6J5Ru/5nsWTpoHu9RBmmyvTiG91xNy/DtOelbjvc1CdUiiMqR90t7Ng+WUMCgZ2qqiLiEmzv23SCz2RYVCuVyryhKrPMNxPAl0TkE8AqmwjnD4FPlp68ZzaMqV/K0AAwiqe56Ri08kY7svZD96u93SoLJ+7uCZZ+FGh4LFCXs2M9+y/pHrt+z1KTdOFXTlW6IuG2v970bkpi9sKln4w4rH9qREiYRNiv96pg23mHV5jCsYm9hlGOIc6WOcvWOO35NBOqdtbSn/JgXP/NOvt7udBR1Q/YIcx3Y+yMb1fV80vP3zMbxkql0hmqwsz93TSMIrIOOAM4BNiNGZFyalDmOMy8L9fbXWe7QSUicgOwFdgFzIQzEorIR4DPq+q/2YawuDH0mbyG0U1t0CZRZ0opDjWKJXJOSAtvtCOpFP39MwmleGdi6a9vSZRJKUpoTnQ6jIouoWSi+CZFFarB/gwR/WONz9e3qzm7oPNkD6MUG+nq+xguRx47Leze1VlTMYOZlOoK6wy5XETO95PLWL6nqr+TuMbTYqGElv8APigiDwD+EfiCql7ZtpLF2SYqlcpeimJ+oEv+mi6luklVnad4KybG+bDOqqp6qqo+GTO76R3AZ0Vko4i8XUQeVnqd2jBWKpU8uwW2T5f9wRoRucz7Ozl1WRFZjwmpuTRy+Mki8hMR+aaIPMrbr8C3ReTy3LVV9Zeq+l5VPQoz4u4FmEa4iMnrSjtKugmlITjQYqa/wusV0tiF9sN3Uk6JVNca+l3o1DLsavsOjp77IUwmti04PkxfzX9Z4UAqd8zNGh3OBuOdE3ad3XJNsMxVIeyO995rP3xm97RNM+aSe8wMdqVjXetOHBZN38OY4yucg7qEpjR+5R/xltDuF0NEVmASPJyiqmG89BXAg1T1HhE5ATN02OVgP9bGJK4FzheRn6vqRZHrLwaOxwxFfgZwIfD/lj5EVYyVSiWPnwO16a8A22idBZypqmfPup3q3ap6j10/F1hscymgqrfY5Wbgq5i8r/61nyUin8EkqDkZM/LuIar6IlX9WukjT55iDJ0vJU8wivobxYkQEA7dA29ekPCYU4ozXtBvaViLrxiaQnpmKcV+sDY2eUQ/S4SbGNINZ98W7IfZKjI1P2BMMabmCXTZJPyBCza12nZ7bEuQwqtN76AkFGapuYBTgbEhmymGUo6lSjEWrpMaRlgyf3UMZbhOQQSbzvDTwEZV/VCizCHArTZA+xiMgLvdBmkvUtWtdv3ZQJgC8a3A54E/U9U7GJLJaxgrlcrcosSDOofjWOBlwFUi4rzFb8Vm1FTV04AXAq8VkRnMr++JtpE8GPiqnX9qGhOWM5DYWlWf1kUlJ79hzP3idRHq0CZcZ1Zgeeb1hkoxtCmGoTn+9cPg5hLF2KguXYiRr/5Cm2Joa9wabPvrqf8kdzw2s/T9wTLH4sGlHfLYU75hEHrunaSUuK+SrL3R2Rp76SVswLdTkLNSjXmESSuKaFKKORvj6PZCg5JOONz2UqoX0zC3g6p+BPhIZP91wOO6qUmeyW8YK5XKeOmwKz0p7DkNYyxJQLjdtYc5pVad2nOB3k4h+EMEY15n/9xQvfjrKWUTC1gunVQpahPs0gsdElOFYaXcfZcF2/75bp+1P25PvL8B9Ue8TE55974XRqHutptOOfamovA+zzCNWZg2LmZznnW/Wfe3S6cO/VuU2EpDSj7K2jBWKpVKQG0YJxg//i6MRSTYbjO8r4Tklyac7CiT+Dyl5NrYC3OKcSx0Z5Fvf59Q4brtIOlr7p2UxoVC5PsxqBx3R1LLjTWOMeo5D5ZdRVbUhrFSqVQi1IZxQghjs5ZGjqV+KYeJfRyGNokummIUY/tKbIwpwufr2cZ8Vbs4so/I8XH/1+SmGIrFRUYosbvmbIwDI4L8a7kROOb+u2Op5lKe6pSdeeDcYBnaFv3vfeipDpexZ5iVdi5SZjd5O+UeyOQ2jJVKZW6oXelKpVIJqA3jBBEbRB8ea9Olnquuc1i2TVe6KRt3LDSlNJdib+ih/wLCIXouqUPo+Mh5sVIOmliAd7hcFpQd4gMbxvni3m/u0rM+NxeiFcvC3cG/Wcl3OhX8Hf6vpMwCEP++1IaxUqlUItSGccLI/XKmynSYGAKIDAVsKOevt3EAjCNcZ9Y7WBZZd0uXxKFk6F5TdoKYYgzvl1KQ4fkwOzQqUR2Y/b5ChTgd7I+Run7UCSiD273jbgBAdoRcnFiPKaxDU/iOT+7jqoqxUqlUAnYzOOhoL2DPbhibArxzx0pCfIZRiuG+NgkhmmyLORtjKiyjSDGGSnE/BomF+DQNH/RvHNoQw7Rjoa0xtS9CiY1xOtjOKcUUOVt37BjQkEth8Brhdm6QQlM6tbb/9QqMYYqbhcye3TBWKpVuqF3pCSFUa8MkhChxcubeUNPbazMsLaUYY3MpF6cSS9w7xizvNPRVWdNFYkquKf1YTOmFNsbU0j8/8SHk3n34nkKbYpsho7n7lMxO2ZYSxZgKCo8p2KbeDlQbY6VSqcyiNowTgPuQwprnPrg2XulRPNgpT2VOtbRJfTXMkMBSpRhVDqFiTKk//+UEyRx65DzYKRtjykvtrxfaGGPJFkLFWOKNHoWkrTG4v09pHGpsX0opxs7JdQr2wiGBdTKsSqXSzEzhXwMisk5ELrBzPf9MRP40UuY4EblLRK60f2/3jh0vIr8QkWtF5M2dPFuEyVOMIbn4qyb1l7MhtXkzpUpxmNRXOdvYMF7pVl0iZ290iq3EXhgbDdN041D9pZTjcvq09EbH3mOb70eXpGyPo3ZXQ2XYlFTCPydHt13pGeCNqnqFiKwELheR81X1mqDc91T1d/wdIjIFfBR4FmYWwB+JyDmRc0emKsZKpZLHTYZV8td0KdVNqnqFXd8KbAQOK6zJMcC1qnqdqu4Evgg8v82jlFIbxkqlksfFMZb8tUBE1gNHAZdGDj9ZRH4iIt8UkUfZfYcBN3llbqa8UW3FWLvSInI8cCowBXxKVd8THD8c+BywypZ5s51gO3NRmmvd1DXKGaG76EKH2zEHwDBDAlMhPiVd9qa6xp5/Jjzoh8vETob4LIBQ5nxJ3W+EJBKOkgDvkkuWhH6ljoX3zTljSr+HJc6Xkq507n7tutJrROQyb3uDqm4IC4nICuAs4BRVvTs4fAXwIFW9R0ROAL4GHEk8Il4j+0ZmbA1joT3gL4AvqerHReSRwLnA+nHVqVKpDIHSZkjgFlU9OldARBZjGsUzVfXsWbfzGkpVPVdEPiYiazDtyDqv6AOBW4pr1oJxKsaePQBARJw9wG8Ylf64r/0Z5SGHCXQtOafLMJ0SR0pOZTaF+OQcLEM5XxyhE8bhlNy2yL5wSGBKbfqETpiUM8Y/1jCsrkRFpxRkrGrDOOnCMk3pwfw6tXHMlPaUYnNR50KUOhwSKCICfBrYqKofSpQ5BLhVVVVEjsGY/G4H7gSOFJEjgP8CTgR+v5uaDTLOhjFmD3hSUOadwLdF5E+AfYFnxi4kIicDJwOw7PCu61mpVHJ065U+FngZcJWIXGn3vRU4HEBVTwNeCLxWRGYwv7onqqoCMyLyx8B5GNPbZ1T1Z53VzGOcDWOJPeDFwOmq+kEReTLw9yLyaFXdPXCSsVFsAJBVRw9eI/ZLXpp2qc2Qrxxhmab5n/31JhtjyTDCUQK8HTmF3LtGqBxjsVJhMHhYgZitMWVjDBVjLFXZEORUeerSpbbFNufG5oYOz2kK8PYJP44mW6O/nrtfhw2jql5Mg8xX1Y8AH0kcOxdjchsr42wYS+wBrwSOB1DVH4jIUmANsHmM9apUKm1w4Tp7EeNsGH9Esz3gRuAZwOki8gjM79dt2as6r3T46xjztI0yv24qcLzEftekBv31lELMDQksvV9qn0+oFIqGpbkf/Ji3uOnFDWNrDPf7dSikRK3nbI0pVTmMYmzjpW7zH5r6LFMK0l9vsjXWtGPdoKpRe4CIvAu4TFXPAd4IfFJEXo/5XXqFtSVUKpWFwl44VnqscYwxe4Cqvt1bvwZjjB2dmI2xSSmW2Gma4v/8fW1iEkeZpqCNYiwl900IbWFJBQmzPcrh71zJ0MDc9RtOTQnVmKhto/DD+6Q82CUe7fC6OVtj7rpNZUtsjCsS5/jUrnSlUqkE1AzeE0hXMYmlcWM5xdjkjR7GK51TmSVKsQtv4ige+1lqryH5Q1uavsEln3kbW22TB7tNDGR4PJYaLUUbW2ZJHGMYWxnSXbjORDD5DWOlUhkvNVHtBJKzBzm6eMphVF+Jh3mYpLNNSrHNKIlRKPJgN+wf9X5NKqmENjbG8D5tYiBTdY15pUujCHIM45WO1b86XyqVSiWgKsZKpVKJUBvGBU4q7VgsPVeKEoN8an+ua9tmjuhxhuuMQpuhZrmhbKlzu6bE+VFKiROryQmzPbOvTV27fK6SrrQL27kzcn4N16lUKpWAGq4zIcSCuddE9rUJM2kqWxLC0aQUSxLVltwvtT0KMRVYGgrTVWCyo81zNYXADMMwQeG5YPBhtkOHTPhOhnm+mBOoxPlSbYyVSqUSsJs2iWr3CCavYXQ2xlzQ6ighG6n9MQXXNnFs7tw2dsOmX+9hFFdOJY0SRJ2ia5U5jmuV2BpLPuNUXYaxI7aZDqEpqYR/vdj/kU/tSlcqlUrAXpbaZfIaxtArHbON5KYwKKWNza800HsYD3POxpi6VgkpNdHGxjgKJdcehw21zf1yHvqSIPDwOqMMJwzvl5sOoekaseGzzivdNBRxL6FOn1qpVCoBk6cYHeEvZomNcRSV0sUQvRKbVVM9mq7XlpyamG/aDH8b5rpN+3PRCyWfcZOqbKMYHSXTIYT3z5GLcRwDIrIOOAM4BOPW2aCqpybKPhG4BHiRqn7F7tsFXGWL3KiqzxtHPRfSv0GlUlmQdOqWngHeqKpXiMhK4HIROT+YVtlNv/xeTKJrn22q+viuKpNi8hvGmK2lC+9fG1U2TqVYYmNsQ0rxzNU3YRgVONfKMUfqvZV4pcNzRpmEyzGMwov9r2Rjf7sb+qKqm4BNdn2riGzEzCh6TVD0TzBzTz+xkxu3pNoYK5VKAy7Cu+SPNSJymfd3cuqqIrIeOAq4NNh/GPAC4LTIaUvtdS8Rkd8d8cGSTL5irFQqY6aVYtyiqkc3FRKRFRhFeIqq3h0c/lvg/6jqLpFZ01ocrqq3iMiDgX8VkatU9T9LK1fK5DWMLlwnF5g6SphOilHCZ0YJ1ymtS1tG6Z528V7H9c1r43AYx/1KzCbjym0YdquHyd0Y7Zp3m0VCRBZjGsUzVfXsSJGjgS/aRnENcIKIzKjq11T1FgBVvU5EvotRnLVhrFQqc43SlfNFTGv3aWCjqn4oejfVI7zypwPfUNWvicgBwH2qukNE1mAm0ntfJxULmLyGMTUkMGZQTi1jtA29KTmnjSOlTfhHSBfKca7PXUh0EdaVU2fDKMVR0o51lWCiR6dZJI4FXgZcJSJX2n1vBQ4HUNWYXdHxCOATIrIb4x95T+jN7oo95atdqVTGRqde6Yspmg+3V/4V3vq/AY/ppCINTGbD2GRjLB0SOEwozjCKMRwa2Obc3DltGCVMpqtyw1ISfD7fAeolNuFwe1zzFKVCeEq+P+NXjBPBZDaMlUplDtn7UnhPXsO4CPOL2LVXehQFN0wKsaZzS2gzjNAx1zbFLoOoc2XH8U3u6prD2JO78FyH/wdD9ziqYqxUKpWAvS9T7eQ1jM4rHaZJiinGURjG5teUqLbrqQ3Gwbi9n11dt61SHIdtddj7hGXGFdeY6jm1tlvXrnSlUqlEqF3phU0qjjFGG290F17iJm/0qF7p8NyQcdni5ivWsbX3NHPuqN/0cY6WiUUMjDJKxp2bis6oirGRyWsYK5XKHFMbxkqlUgmoXunJwO9Kj0KbLnVsbuim+aNzAd6jhOmMI1h7rrvY47r+XAV6j7NrDel5pUu61qlkEmBP6lwAAAhZSURBVCVd9uh3qnqlK5VKJaB2pSePmOJKZSMeJfSmjWIcJXg75zTo2qHQxCSmF5vrb3QXCShy1w2/S20UZOiEyZXJKsa9rys9tgzeIvIZEdksIlcnjouI/H8icq2I/FREnjCuulQqlVFwirHkb89gnL+vpwMfwcwIFuM5wJH270nAx+2ymRn6v3Cx+TPWeOXC8/yybeyF90TOKQ3sbhOsnbMfdjn3yaSnG1sIdfBpM19LF2nNQgXp70udE7tv7n+ix96nGMf29VLVi+ycDimeD5yhqgpcIiKrROQBdrKcSqWyYKjOl7nkMOAmb/tmu29Ww2gn1HGT6uzgL+Pd8wXKGmDLfFeikEmqK0xWfSeprgAP769uOg/euSZddIBJesYk89kwxpJVaqygqm4ANgCIyGUlk+0sFCapvpNUV5is+k5SXcHU162r6vHzWZf5YD6nT70ZWOdtPxC4ZZ7qUqlUKj3ms2E8B/gD653+LeCual+sVCoLgbF1pUXkC8BxmAm4bwbeASyG3oQ35wInANcC9wEnFV56Q+eVHS+TVN9JqitMVn0nqa4wefXtFDFO4UqlUqk45rMrXalUKguS2jBWKpVKwIJtGEXkeBH5hR0y+ObI8SUi8o/2+KUNweRjpaCubxCRa+zQx38RkQfNRz29+mTr65V7oYioiMxbmElJXUXkf9r3+zMR+fxc1zGoS9N34XARuUBEfmy/DyfMRz1tXeqw3RSquuD+gCngP4EHA/sAPwEeGZT5I+A0u34i8I8LuK5PA5bb9dfOV11L62vLrQQuAi4Bjl6odcUMKf0xcIDdXruQ3y3GqfFau/5I4IZ5rO9/A54AXJ04fgLwTUzM8W8Bl85XXef6b6EqxmOAa1X1OlXdCXwRM4TQ5/nA5+z6V4BniEgsaHzcNNZVVS9Q1fvs5iWYmM35ouTdArwbeB/jm6qphJK6vgr4qKr+GkBVN89xHX1K6qvAfnZ9f+YxdldVLwLuyBTpDdtV1UuAVSLygLmp3fyyUBvG1HDBaBlVnQHuAlbPSe0S9bDE6urzSsyv8HzRWF8ROQpYp6rfmMuKRSh5tw8DHiYi3xeRS0RkPkdplNT3ncBLbQjbucCfzE3VhqLtd3uPYaHlKHGUDBcsHlI4ZorrISIvBY4GnjrWGuXJ1ldEFgEfBl4xVxXKUPJupzHd6eMwSvx7IvJoVb1zzHWLUVLfFwOnq+oHReTJwN/b+u4ef/Vas1D+x+achaoYS4YL9sqIyDSmW5LrFoyLoqGNIvJM4G3A81R1xxzVLUZTfVcCjwa+KyI3YGxL58yTA6b0e/B1Vb1fVa8HfoFpKOeDkvq+EvgSgKr+AJM4rDRBw1yz9w7bnW8jZ8LoOw1cBxxB34j9qKDM6xh0vnxpAdf1KIxR/shJeLdB+e8yf86Xknd7PPA5u74G0/VbvYDr+03gFXb9EZiGRubx+7CetPPl/2LQ+fLD+arnnL+X+a5A5gM7Afh326C8ze57F0Zxgfml/TJmSOEPgQcv4Lp+B7gVuNL+nbOQ321Qdt4axsJ3K8CHgGuAq4ATF/K7xXiiv28bzSuBZ89jXb+ASfN3P0YdvhJ4DfAa791+1D7LVfP5PZjrvzoksFKpVAIWqo2xUqlU5o3aMFYqlUpAbRgrlUoloDaMlUqlElAbxkqlUgmoDeMehIisE5HrReRAu32A3R5LNh8ReY2I/IFdf4WIHOod+5SIPLKj+/yuiLzdrp8uIi8c8joHici3uqhTZc+mNox7EKp6E/Bx4D1213uADar6yzHd7zRVPcNuvgI41Dv2v1T1mo5u9efAx0a9iKreBmwSkWNHr1JlT6Y2jHseHwZ+S0ROAZ4CfDAsICLrReTnIvI5m2fvKyKy3B57hs0VeJXN17fE7n+Pl1PyA3bfO0Xkz6yCOxo4U0SuFJFlIvJdN4xQRF5sr3e1iLzXq8c9IvJXIvITmwDi4EhdHwbsUNVZ8xWLyLutglwkIjeIyF+LyA9E5DIReYKInCci/ykir/FO+xrwkuFfb2VvoDaMexiqej/wJkwDeYqa9FcxHo5Rk48F7gb+SESWAqcDL1LVx2CGuL3Wds1fgBne9ljgL4N7fgW4DHiJqj5eVbe5Y7Z7/V7g6cDjgSeKyO/aw/sCl6jq4zC5H18VqeexwBXhThF5H7AWOEn7CRhuUtUnA9+zz/FCzFC2d3mnXgb8duKdVCpAbRj3VJ6DGer16EyZm1T1+3b9HzDq8uHA9ar673b/5zDJTO/G5GX8lIj8HmZWx1KeCHxXVW9Tkx7uTHtNgJ2AS212OWbcbsgDgNuCff8PsEpVX62DQ7fOscurMElVt9ru83YRWWWPbcbr8lcqMWrDuIchIo8HnoVRSq/PJBYNx4Iq8TRT2AbtGOAs4HeBNg6MXPLg+72GbRfxNHjbMOPifX4E/KZzMnm4rEW7vXW37a691F6zUklSG8Y9CJvB/OOYLvSNwPuBDySKH27zAYLJEXgx8HNgvYg81O5/GXChiKwA9lfVc4FTMF3ikK2YlGUhlwJPFZE1IjJl73Vhi8faCDw02PctjGPpn0Ukds8cDwOic5xUKo7aMO5ZvAq4UVXPt9sfA35DRGKJcTcCLxeRnwIHAh9X1e3AScCXReQqjNI6DdPgfcOWvRB4feR6pwOnOeeL26mqm4C3ABdgMspcoapfb/FMFwFHhdNWqOqXgU9ickUui54Z52nAP7coX9kLqdl19kLEzKj4DVXN2SAXDCJyKvBPqvqdDq51EfB8tXPEVCoxqmKsTAJ/DSwf9SIichDwodooVpqoirFSqVQCqmKsVCqVgNowViqVSkBtGCuVSiWgNoyVSqUSUBvGSqVSCfj/AVj2/pReU6LdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#NBVAL_SKIP\n",
    "\n",
    "# Show what the update does to the model\n",
    "from examples.seismic import plot_image, plot_velocity\n",
    "\n",
    "model0.m.data[:] = result.x.astype(np.float32).reshape(model0.m.data.shape)\n",
    "model0.vp = np.sqrt(1. / model0.m.data[40:-40, 40:-40])\n",
    "plot_velocity(model0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASMAAAD8CAYAAAA8P8JjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnX/sp1V159+H+cHIAMIwoANDO9iytq4bF0MsarNrxO5SyoqbqJG2LmlJJruxK1qNQt2E3U2b6Nb4Y2PDOikWbAmISBfimqqhENNsZZ1RV5GR6qrAwMAwBQRBmF9n/3ieO9/7Pd9z7o/neT6fz52Z80q+uZ/n133uc5/ne++5555zLjEzHMdxFs1xiy6A4zgO4I2R4ziN4I2R4zhN4I2R4zhN4I2R4zhN4I2R4zhN4I2R4zhNMJPGiIguIqL7ieiHRHTVLO7hOM7RBU1t9EhEqwD8A4DfALALwDcAXMbM9016I8dxjipWzyDP1wD4ITP/CACI6GYAlwIwG6NVRLwaS2JaSA/F54hjATZSbd8hY792jYQy29a+1P4UNde0aEOfK9OQOpk3JfUqz8lta5BIU8eGfBch3Q/sZebTAeCiiy7ivXv3FuWzY8eOLzPzRRW3HsQsGqOzADwUbe8C8GvyJCLaCmAr0DU0mwCc1B9b16fPR+fLY4HQwOzr0/3Rsf0ifd44V7vmYJ/KRjBsr8FKrAZzVWZbo2b8fCh/isnB/CkzybekDhaN9gyyrg+K/QeN/Vp+8ntZK7aBpe9Mfnep+pP3Dt/9I8AD4Zy9e/di+/btiVyWIKKNRSeOZBaNkdZ4r+ggmHkbgG0AcDwRr8JSBYf0pOj8tViOfBnyeArZsMR5hd9a46bdP9VolDRCUyjtcnmkGquSRmFIgzVFY7OoBkt2RPHza1L7UKxGKO7o1op9VucYE8orO9blMIADhSWdD7NojHYBODva3gzgkRncx3GcwTCWjz0Wzywao28AOJeIzgHwMIB3APjtkgulzmitciwnbcyqN5USUo20MKRMY55Dlq1E+horPc2DWdmhyGfXJKIpkRK5lIji736NkabeiawnXbo/BiQjZj5ARH8A4Mvo6uwzzPy9qe/jOM4YjoHGCACY+UsAvlRzzUFMMw6Pewxr7C+lLE3JaJ1T0luOkSRmrWuxyp+SOqZ4LzX3yzGVpCbftfWc2jeVQ/tepBJaplI/lDpHfp8pZbte18dIY+Q4Tut4Y+Q4TjN4Y7QCRidmBkWbJuIeJ/aVKBnlEEtua/cpNRGQIr6GZosUl6OGIddodWPlk6rHeTswlj5rTZ2k6qJ0uKZh2RfJ7zTOXw7X5PAstqWzFNfjh6qHALwwOpcpaaIxchxn3vgwzeQQlgwNQxr3yHIatGbaVfbssieMpaFw75wUJQ3LSu5bUkatRx1KjSJ7VsaPNflPcU0uryHPEF9jTWRoFtfA8m9ASkCWRKQpsC3FdYrUt9kxTWNERJ8BcAmAPcz8SnHs/QD+FMDpzJz0P/EQIo5zTBIko5K/LNcDWOG7RkRno3OYf7AkkyYkI0bXiktfMa311/x3LKzeQ+7frxyTPaAmRclrLeO5knKlp2HLsK7VdCC5nrVG3zSEqf3zcki9ITBO0pPfh9xOGTBKCUhup9xBrDqJ95dJ19MN05j5a0S0RTn0cQAfAHB7ST5NNEaO48ybQ5ilOwgRvRnAw8z8f4nKYg000RjJ2TRNMrLG/LKniLdlb2j1GIeU39a5smwpr2xLLzO11JGTIIZEAJiXC8gYCTfFGKmnxPPeckatMWC00ng2LRUlIr5vXEehvHJ2eiXFktFGIopd/Lf1ju4qRHQCgA8B+FelNwAaaYwcx5k3VcO0vcx8fkXmvwTgHABBKtoM4JtE9BpmftS6qJnG6CBW9jZrxHHAdtPQelhpNyJ7WM1WKOcGknIPkb1mSmorpcYhsubawBCn2imZSneUc3YtsSGyjqVm0+R3Ybl8AHldUdiO66TEMTYuV4ws43JmN7XPzN8FcEbYJqKfADjfZ9Mcx1GYbjaNiG4C8PcAXk5Eu4joiiElakIyYnQtu4zamBpvD7HHsXRHmjNjTiLS7D2sGRZNvxSYsjdYtC3PEMbojDQL/YClLyxBvr/4vcmgeznnV+0blhJR0C9p+iEriF9KwrOswZcz6WzaZZnjW0ryaaIxchxn3hwbwdUcx2kedwdJYhmQxfuGTD1bimwtL0sRmlNcascspXvMFFPp1rUtKwRLnteacEjlMWRKP2fAqH2HgZpp+lz8Im1oFxgy3JRT/MsJBjXt0FRj5DjOvHDJyERTIu9XjpdIG6VokoPMP2dSoK3tlutpY0qNOWvILZW0CIZIsoES6dGqvxJFr7WdMt+wlhAKkpBm9LjWOFYT37oEaSJjO8x6Y+Q4zsKZrTvIEJpqjEKLHqZPj4+OWY6IkrG9SmoaXitHSrqyel5Nz1QjDZQ4S1rXzlJKmtqoM0fKAdjShmj7pQQbJAlt7byQvxUOxNof76sJmGaZo1iLScbH9inHlvBhmuM4TeCNURKpK4p7JE1aGpr/rLCcaFO6iVKXlRopR+6vka5yedcyizqv0b0FNCkqFzJWeyelzq4po0cr/OwQ9x3NmTe1fPtyvDFyHGfhuGRksgorXT2eiY6v71NtFsuitPcfEmIjFboh1+Nq4XRzYU5Srg9W+VPPVxP4TTKlbihl6VIaAiY+N5e/NjOWe29DwoCkJKMhIWQtpH4IWJKEgnpaz98bI8dxmsBn01QIem9rzZiNZVGWyVqwKyktWfZNMZZ+SUpZKYlIXltjmzRl/WnB8HKU2A1bs02pQGnSPkfOnAHlVtTS+TX+bUlENUslpQLASV2R/X/kkpHjOAvHh2mO4zSBN0ZVlKzumhrOBHLK0yHXDiG1Uq5lzJYyiCslNVTI7dfuPWWdaG5AOTeN0vyAMqdXK/Z6arWOUqNHTYFdMpy2sFyNNAV2SGNn3SW8MXIcpxnca38UqWlXQJ82t7ZrKJHSJFLy0hSJlntJSpFtxXeWaFKVlMQsZap2jXZM29awpJ2pFNhW/ilnZbmyh/ymUtKNlcp3rUmlQyQi63k0l5Xnxbn6fXw2zXGcJvBhWhGa7kDGxZYGkmFbC+lRY+BnSR2h56nRY8lzU72zvH/KUTYnEaWMIi2JKCUp5fRMJVixy4dM05ecK+takyDkVH7JCq+lRo418axT5Bx+NZeP/JQ+0GJj1HIwQMdxZsakq4N8hoj2ENG90b4/JaLvE9F3iOivieiUXD6DJSMiOhvAZwG8FF1HtI2ZP0lEGwB8DsAWAD8B8HZmfjKXn6aj0Fp7y4S/JADXFC1vyYqyskdKhT2xenCpO9IklRwpB1nLSVPr0S2JaIiuSBr2pd6JFeguUGPAqIVgtYwbrRmy+Jhl9Jhy9bDqsUY6lM8lt4Glb0gzvFzOZJLR9QA+ha49CHwVwNXMfICIPgLgagAfTGUy5v/zAID3MfOvArgAwLuI6BUArgJwJzOfC+DOfttxnKYICuySvzTM/DUAT4h9X2Hm0Np9Hd2qskkGS0bMvBvA7v73M0S0E8BZAC4F8Ib+tBsA3I1MixiQPUWJTYjURaRWh5X3mTWhzJakFB+zAnxpvadlm1SiG7MkolTo05qZt4CUBsfoimrOlfWX0q2EY9Y6ZnJ/fExKRLKOpIQUU/P9WZKe3NbsjHT7osBcdUa/j260lGQSBTYRbQFwHoB7ALykb6jAzLuJ6Azjmq0AtgKuuHKc+VPVGG0kou3R9jZm3lZyIRF9qL/RjblzRzdGRHQigC8AeA8zP01ERdf1D7MNANYQsdYTp4KrWTMhWkiPEktrSa730kKZyN5fSkSapCT3heeRerOU/ZRE9tKpkCWWDiS+ZkjIi1lKRJrknNMZaSFYrVmzmuWGLImypq5SIwKIfdYIQQuNm793cWO0l5nPLz05QESXA7gEwIXMzLnzRzVGRLQGXUN0IzPf1u9+jIg29VLRJgB7xtzDcZxZMNthGhFdhE498y+Z+bmSa8bMphGA6wDsZOaPRYfuAHA5gA/36e1D76Fh+RpJHQtgSyaBEjujnF9bSZgTS3+i7bPsZFIzR7kyaDM6paEw4t9DZtEsiS91TSklOsVUOA0p+eQkpdSxkhCy1qIL1na8z5pN0yQjaQ2u22lN1xgR0U3o9MQbiWgXgGvQzZ4dD+Cr/Wjp68z871P5jJGMXg/gnQC+S0Tf7vf9EbpG6BYiugLAgwDeNuIejuPMhOncQZj5MmX3dbX5jJlN+zt0cdE0Lhyar+M488IdZU3k8COljM7FLY5/l8aMjskN10oMCkum9q19oc/SDP7kPmu4pg2rpCGcNY2tuTEMWcVC1t8+47z4HOv9WIaucT7WyhglURtzw7b4tzURMIXLR0oxL6f2NTeXE4ru2p47SFONkeM488Ibo9FYkpDWW1quIiWKbHlM9vA1zpsS7Vr5HFLxO8QdRJvatxTYlpNo/Lum17cUrtZ58e+colyTgqUkJKfyw+oy8bp7UgLKKbIBW9Ff40w8JFhcToGtSZjp+3lj5DhOM3hjNIhc8CztPCtg2RB3ECkhyenTGoZM4WrPZa0hl9Jv5fQjNZJRypBRC4cak3LfserU0p/E97FCrmoGjLln1+oipyMa42qkmXFYQdQsnVg5HlzNcZwm8GGaCiM95gVsHZE8viazL95fItVYM3EpndGQ3tFyZ0j1ljndg9RvxPssqeB4cV58vTVbF9CMEC2JIaX3yaE9vzR6zEmA2j5LV6TVhRWcLoX1vlIzwlYIFKkTS4V8MWGf2nccpwXGzMLMgCYbI623KZ1F01wtagKwBUqD0JcEyq/pPQPy+Z5VjkksB9kSF4jjjf1xvpZkpPXOOfcPTd8l9R85aUq7VoYFsdLUsVQ4FcvWKufqoSHfoxYATtoTWbojLRxysgxyONIATTZGjuPMGEb5uHhONNUYWWEtANteZdaNe862RltqZ8jMiuz1Q0+Xmi2Rzy4lIS1sqtX7p2bTrEBhJXZHlo4t5B/PttVKkPuU39ZzSckv/l0jGdW+2xoJqcR+yppBrS6DS0aO4zSD64wcx1k4LhmlKVkRI5hpWVP72hT4kDLUmPfLY9YKGCXDEKkM1owrreFZLopj/Nta1UIrY25oohk/ynOtWOXafXLGjymj0VLnV22fNGsYsjJLTSTLXLxuYKXCWiqyU0PlbFm8MXIcZ+EwfJhWQqij2Fg95w6SCstRE5XRwpJqNAW2NcWvUaq4ju+bi0KZclHIuS+MWRyhxil0iLmDRAshEsJn5CJZArZEZBl5Annn3ZIyWqkWtdGKYJlad66ojWHYvjoLosnGyHGcOeCSUTmpIW2N0aPsgTRjQEluLXqI/fJ3Kg8NKb3Ja+JrNefZqRmyqkUNKYk2YMWITl0rg6hJaUdb6WOdce4YKbGkTnKB07R9cjv1DSfdhlyB7ThOM7hktBKCrhPRqHERyAVEq3EPscpRI0EEtHXTpK5I9s7abFou9K7W0wY1wTpxzAriJe8N2K4xmuOqpQ/RXB+sgGFyJim1IkZpCuQlotQ7tr7RWc+m1UhcSSaUjIjoM+jWR9vDzK/s921At4rsFgA/AfB2Zn4ylY8v5uo4xyKhMSr5y3M9gIvEvqsA3MnM5wK4s99O0oRkVELonSyzeOlQCNhBwOQ1JbNsJb3lkHAjMniWzF+bBZJ1IXvPMAu5TmzH18qA/yVo0kxMXPdyRjRc86zYHzsAPydSee4L4j5xOUpXh9VCbUgdYk143RrhIufsrUl8lkQpv7/4Gyt6pxP6pjHz1/ol7mMuRbeWGgDcAOBudIs6mhwxjZHjOBMzWwX2S5h5NwD0q0ufkbugmcYoFYpDIxdKJP5t2XVoywBJ/VJOJ1QToE2zOwq/pc4o9PT7xPH4WJAYcksIxc/3vEhLA8vHWNbv2jnWewrP9QKWkGV7pk+fFWnYH5dROsJadkZxHUmpsySErGU7JvWTKaHD0ulJSVc7J2djBhQG36szetxIRNuj7W3MvK346kKaaYwcx5kz5ZLRXmY+vzL3x4hoUy8VbQKwJ3eBK7Ad51gkSEYlf8O4A8Dl/e/LAdyeu6AJyUhO7WvkGnHNEM4awklXC80cYIjrSG6YpB0PQ4VQBjnlvk4c13iqT6VSOOyPh3ilSv3UNL1mRmFhDXlKDEClQl4q5k+KrllvpNLlo2bVE41aQ8/U92gN11IuJIHRUsSE7iBEdBM6ZfVGItoF4BoAHwZwCxFdAeBBAG/L5dNEY+Q4zgKYyOiRmS8zDl1Yk0+TjVFNzyslGa0nyimuUy4kY+JYB2rcQCy0FSqCU6hlLKitXSYVuVKCCHl/lTlTojLeRrSsDHKqPzY7yD1HIGXAGOpEPpcWQmSKtc4kJe4tloSUWvUk97+gTcIkcXcQx3GawRujPFoPkau3EsNCKYWk7jNE/5O7tiRwmQyMFnryeKp4g3E/K5BYkBYA4OQ+DfqW9/Nd/S+pafp0dJXU2hwQaWD1it+f50/02ycuS2+gdwBYmqaPSyCn9IMRpJSqNGlAmkqkjB5TQdQs5Hcot1O6N0tSt1Jt32RSnMczchynGVwyskn1EKWrqMZYjrElUpR1X81QEuKYta0he+V1mfOAlavCBsnnlD4NUs8GsT/e908O6xz/a5/+tE+fEikA/KxLuNf8WBZ9cSFXh40TRdqV5nJ+HQDgafrfhy8Jd3zCKEkoYZCcno5uJ11F5HvU3kWplFEys2gZMGproJXqjmaKL1XkOE4THI0KbCJaBWA7gIeZ+RIiOgfAzeg64W8CeCczF1k0pHoIWW9TjJ1TPZGUhMbMqpXY2Mj8pf4ntqkJv0/r0yDtnC7Sk4PI9GzkFnSgN4Tde1OXWuLGc0uXrFDWlEhGa/paXddnuL5PT3hk2UOc/MjSJSf3ItwvnNfveLxL9j2xbPOw5BS2gZVSlNQzaRJTqb4nZfdjhUTRrpXVZznB1uhKR9GYzmgKC+wrAeyMtj8C4ON96IAnAVwxwT0cx5mSaUOITMIoyYiINgP4LQB/AuAPiYgAvBHAb/en3ADgPwO4NpVPqBfZe01VD5YtUso2SUoqlu6oxFo3ZeFrLakTpJ+g7zkdS4TfL+3T04J58Zl9elaf3tinOyK3oCBW/KPYfkakcWwPaxrLMgCKf4cHCoY/UsEVK7Re3KfX9+mnumTtY1161sN92ktTPw9SHYD+lMPSUkiDpKQ9Vk7Q074/SycknVyDDkuzo8rNos1NZ3SUDdM+AeADWD5yeIqZw5zvLiz9aziO0wpHkwKbiEKYyR1E9IawWzlVNeUloq0AtgKdxBD75JXY8gwhJyFp987pjuLyWIsPWrYvQF4SCumZWOKwJHS2SMNJ/6FPv9GnsXLlCZFaokOsM5LdvtV1a1N+0lw8SEgniRRYqQS7pE/v69NHlqcvemjp0i39vtP753oUy9NQBZEwdfiRpRVVKkSK5eeYW2wxPmZZmmvfdqmUNEiH2pjOaIxk9HoAbyaii9H9T52MTlI6hYhW99LRZix9Qsvo46FsA4A1RNP4HjiOU0aDw7TBCmxmvpqZNzPzFgDvAPC3zPw7AO4C8Nb+tKLQAY7jLICjSYFt8EEANxPRHwP4FoDrhmakDZ/ktmxNhzjXpqJD5hTZKQdXK4517J4RhmVyej4op8MIbO3G6KIt4mBIN/fp/X0axiZBWQ3kFdZheJZaoiKQGk/L+ByWj0o8TJPz8yF9VZ+GytGsOfvf6/uh2y/16fpgYdCfpoVTSTlOy+3wWw7DrHAncSRLeY0c8k09cRNQh3BHqzsIM9+NLuA2mPlHAF4zRb6O48yQxoZpTVpga64RWsiJHDnpJrVGmKWolitJaOFHZF5SItLcM4IkJGfnV4UdvxhdJCWikIYbBK2tnMYHVkpEcikOraJTMS7ia7R9sqIsjW/820pfLfLU4qqI9KUPdOlxfR4pRW9uqh+wV/KQEpEWz9o6llJgz4yjaTbNcZwjmAYV2E01RlL/ozkZToEc95dgSUqazkjqiuT0/YboXKkjWiERBT1QPLd/pkhDxpalXzxNL/0krGhnJR6eKUpj1a5Sfls+N9/s06BDSuk8xH3P6CWk/dE1cmo9FfjNyNYU4mT1audMqSOK8yh2WToadUaO4xxhuGRkcxyWJIkxY2hN/xN6CLkGWUnPJNfBsvRO2jXBNjGockJgs1hnJI0bV1ler7E/SDAODBnlopGluucSC7yahd4tZAWGl6EtmCf1P+G5QoWG546fSz6rMa11+iP2Jc8tP1WVekslopLVYa1ZtalR7XcabIx8qSLHOVaZcKkiInovEX2PiO4lopuIyArNZdKEZFSzVJG0CbHsgrRjMq8UIZ/9YjsltcmOPUhgwQMi+IEq5jFYH3r9DSLVHEqlK0XQFVndclxYOYtl6XS0jzC3lGyJT0LJOlFSzAifdRBZwnPHXq+hfoJ4E6TF05bvXxtFZNvws+WnhmzDKamVeq3lhqw0/p1TwaW+T8u+rpoJZ9OI6CwA7wbwCmb+ORHdgs4Q+vqafJpojBzHmTPTD9NWA3gREe1Hp5lQ3cByGSwcQtcLSWkk7jTlLJZl/5MKXDYG6eyqCQ7SrshaPic2Oj4s8EhpR6bh4jjDYE9kSUKaY2vNmjfWPusjLlG+yXNT5s1BWpRKllv79E3RtSeI1KrHkw9fgZN/tvyQfE/BqbbENknaDGlCqSVVTxFSeZCj7ESNETM/TEQfRbdY488BfIWZv1Kbj+uMHOdYpG55641EtD362xpnRUSnArgUwDnoDE7WE9Hv1hapCcnIcZwFUC4Z7WXm8xPH3wTgx8z8OAAQ0W0AXgfgr2qK01RjJEXNWL92SJyTW8VDOybRxELLHaTGDCAM08IowxquAcDaNWKntRSqFkUxKGtzITJTQ6+WpndzL1UOR+P5mtyicTLSJIBV/Qta/8LyS6RlgfadWHr/VDDMIbajQ5CqDNNRdjp3kAcBXEBEJ6Abpl2ILi5+FU01Ro7jzIkJFdjMfA8R3YrORv4Aumgd22rzaaYxiltvLTLFGIOwmnAjliRkCRKp9bhk6JCwfXx8suzRrSVQtaVQra521lZ0U5ISZa3K1z4QawlZuV+RptZkJKOalWZTK8o2x4TfBzNfA+CaMXk00xg5jjNHGrTAbrIxsvRCJWgOgwGpdyrJJxU8TWKtBiI752X3l4oJy1k0vii1rG3rDKlYi5RYKl+CZvuxanlRhqyJZ+l/Ujajzby2ZgrS0WRj5DjOjDlaIz1OzZgGe4qVZqfMx5kzE/b2Y1YNbh5GXQydOdBkY+Q4zhxwycimxjymBqk+sI6njk2hT1BnWKxIW6mLLOXXkdA915Q5d25JcLVUjOH+txUhpcah2tqvhbSpyX9muALbcZwmcJ2RzUGsNLTVrFetbW2yJCcJlVhg51aJ1cpkhSbVlq/JRuXSYlGEsb71IFahY1qWomS5LQkp1nnI0ChW/Nc4INsL+iWpCCyyCDKV6w+kvuGF01iBmmmMHMeZIz5McxynCXypIp3QSJesmmDFcdFGHXIUY23H9nc5m8OaYIa5sMwAsK8/eW3Y+axIU3GsLZeHEoNJiGOpCIzzovRFhTT+ZwpjX6v+ZHxwAAf3py+xFjTRihyKJEfQKSPcWQkmxSNwl4wcx1k4rsC2iWN/10zppxTLpR4WcU9iOUmmvApkGaUCVHbKz2CJw6Gaw1pnJxknx+Ehw7GXiHMsyUhbgaMmrJ6UlrRuP94f55+rwNS6aZan8b/t08eia61KlmvIRTGwreXltKXjZBGlq88LYlv7lq0IKZKSMDWp77A4P5eMHMdZOC4Z6QSdkbRO1zpNq/5KbOdyKoj4t9Uph1QTLGRvKNfjCj1wHMXicOjmvmt90RP9DhkMLJaMpLQUzgnSQbhBKEC8Jn3owi3xUIu3UrokhSaBWalUtsS/ZRqeJ8ReSYmYTxlpX6/7frZ0iZSMpM4oNbUv/ZuPF+eWCB0lzrUzw91BHMdpBpeMdLT14taJ40B+NlKTpuR2QVSJFR13icojYBk7StUOsDKg11l9D75KSgNxZchjYW0wa7H4eCbOkj6sOKox1hItGlZlW6JnvE+GjJVL8v6jSIGlteNk+tjy7ceXrgjC0goJqURnJNVaJfHsLGnpebE9VVDBJG5n5DhOE3hjVIbUuYwlN6sWqzqs2bS1YjslFEgpTkpGmmpFbp/ZL4FHKfFNpmERtpQ/gyUBpZQd1vK6h8R2jKWYCxWprU4g9WNBEgrP9co+fUSkAPCQkfbnPNqX9dHokiBYBclIzqZpSAHPEixLdEcyztyYULUpCd3Eh2mO4yyco00yIqJTAPw5uj6LAfw+gPsBfA7AFgA/AfB2Zn6yJD8ZdkFruKeovxpTl5zDbKpHKpHwLMEkCCFn7urStXEmllPtmSLV1trJiYdhW7P4tuyLZN7x79wSQtryuiHd0Kfn9umP+jRIREH6iX/3Ke9afmqQiGLTpCAZSV2RfCep7yNW5ZViRTeRxvAzZWJ3EK0tYOa/r8lj7IqynwTwN8z8KwBeBWAngKsA3MnM5wK4s992HKc1Dhb+laG1BVUMboyI6GQA/wLAdQDAzPuY+Sl0y9ze0J92A4C3DL2H4zgzom556ySJtqCKMcO0l6GbKf0LInoVgB0ArgTwEmbe3RdqNxGdUZLZcZi9j6blOpLyYrAU1yWteI3hW7A/k4aSwRDvpXuXrjlNGvjJqe5f71M57AGW5rZf3Kc/FXkpDqVFVoCA7nFsrewahmdx2TaI9PQ+/VafPtynu8Q2cHh49nRv1CiHZaFqfhpdIodn4R2k3ptmq5m7RiKHZ+G+2vdojaQmCUc1nc5IbQuY+dn0ZcsZM0xbDeDVAK5l5vPQfcbFQzIi2kpE24loe2NKfcc5+gkK7LJh2sbwv9r/bRW5jWoL4kyGsgvALma+p9++tS/AY0S0qZeKNgHYo13MzNvQL4G7lohXwV4YdSy5UMopP9KcsWOqNbfKH++X0/+Wn2dsrHd679Lx0h936WlBUnigT4Oi9/PKxUFEkBZ/ltdoXDgZW0MSV2xuzfsgmcWSUTDeDJLRf+tTqYXut5+OyhgOhUfNGTQCKwU9KzxN/LhrUYc2yFWUAAAWFElEQVRmVWGFG5lV8E1TACr/B9vLzOcnjlttQRWDJSNmfhTAQ0T08n7XhQDuA3AHgMv7fZcDuH3oPRzHmRFhNq3kL5eV3RZUMdbO6D8CuJGI1qLrj38PXQN3CxFdAeBBAG8rzUyGZUj58Vk2e6mFSi0JqcSesCZ0g9UTaWElLOdgKRHFtoFSqNjQfzCnP7A8PfWccJMw1w8cliGe7EtjhdrQJKNQWCvWS4lrh9QZvRhLhAf65T7tH/7Z55ZtHk5jbxBL9WVN16ewJKTcPo2US9NasS1tSuPfVgz4kiBuKtPbGWltQRWjGiNm/jYATXy7cEy+juPMgQkbo0RbUEwzFtjaeDEen08Z7SCl98lFPh0yrk+tbGKFkUjFLZN2hNKP9LDk1He1p9CS30RQx/xTfmf349R+mu4XgmwRpu2iWBuHf/dyBhtvg+KaDKU7UaRB/OlKuYfuOnxFkGpCSeQMmBXyIyrZCj1QSZy5mneaOze8Py0yrnQhyX1rM8XjGTmO0wxHkzvIlMSNdMkYPSU5yH2lM2TaudasWck6bTndUXyOtDmxPD60a6SuLfS8KV/U9fSXAID/dNhiP8gZQf6IJSMxnUYH+u0DWM5q5XeQiJZLSp+mzhY2jo8WSmCFgU2ZO1l+uVKK1NRaktzsa4wULqTTa3w/6VUT7v+82K9Jb5O7iPjqII7jtEJjgtGR0xjlfDQDmv6nNDC/tU8rxxA0fVBuAdQgLcSr0FqrGFkG0pq+KfTKt9Frl20HCeJu5hXlH8K/IQKwcsFXLYBZ6XOE9PjoWCi3jNYbS4US+S5TM7EW1myX9h1Z0rYM6p/6Hq0IwLXfZYNO+0dOY+Q4zrQ0pr9uozEKrfQQqaNE/5OzGdJm8qaYRQukLMst37SS5Y3krJKUCjQ9iaUDk1LWP+slGsCWVEo+ZmuGqCT0imWdHvKMo4+EmUSrTKnvQ56TKqOFJbmnrPvltrZAhPxGp5JmXDJyHKcZXDJyHGfhHEJzKxW11RiVLDohKVFG54zNUtPzNeR6mhIFdkitpeOBlcOzYEZ4UiYFloZy0o9VGyIErPeRMuYMyGe2FPbxb2uKXxo9xuFAApb5hlz3DsgbttYMY1LDM7lPfpfWFH98jmbiMRaXjBzHWTiuM0pwCMOWf7eUgPFvea6139qXIuWgaEkOmgFj2LdPpCkJwlp0VsYtiyWj8NuaAteU3rK+rPeUei6pIJfB44CV68pZkkpIY8koSJIhf2nsKOs1PkdKbylF9hSrwFrfrnSc1faFVMbNjik1vvXGyHGchdOga1rbjVHc6lsrdpbof6wp/BpJbEgvYq38UaIzstL4emtxVqkP0hajtdJgSBjrVkolI82xOaXTk8h3LENuyLLGhqBBIgqSXkk9SvcMSdhfskqsteKHlp+1yK6mt7Pqfr84XqvrbNAbpO3GyHGc2eHDtAq01l5KBVIPVLJEWIk+KDezV7JSqNxO9c7yHKvHjctUGhI3NcOYk3aGIqUKWZ+afsSa9bQWpY3LLOtPSlVa3Qdpynr2IbO7JVjfo3w+YKUEpNWblX8KV2A7jtMMrjNSCK20ZRsCLNcPAHaPnjL3H8IYW5OcjU382woZklsdaCxW/ZUsNHCc2B5yP+3e1nbufKBuVnKV2GfpI7Xns46l3pOcFcxJSPFvmcq8NBeSVFlcMnIcpwm8McogZyJS4Vlr9EA1EtIUL8hy+EzpLWSPbumOgLx+x1qsQMvPqvMah9IaSmamLFIOzjKomWXprekULTQJXR6T1AT+lzoxzQ5NWo7L50jVZ2o2sMXZtKn1lo7jHCFMtLo1AICIVhHRt4joi0PL05Rk5DjOfJjBMO1KADuxFM2lmqYaI1k5KQO1EsVrSdzqoWVL7ZfTypZrhLbPGtLFPZR0ZrWMK1MKczkstFwjgHLxuWQ4KOM3abGdtKGphjaUtIZnqfjS8v5yhY+4jFOs4CGHT9K1I36/8v3IYVtwp0lJMDnThbEQ0WYAvwXgTwD84dB8mmqMHMeZDxO7g3wCwAew3A2ymiYaI0LXK6R6RsvYMWXgB3FMUtK7yRdWooC1JKKUZGRJLJqkIutAlk2uNBI7h1oSQ01d1Kx6kpP04rKlXClitDKHOpEOuXJ/SsqzjFW1SJnSvGEIlruStoKJNH60vgEt/wkcZTcS0fZoexszbwMAIroEwB5m3kFEbyjPciVNNEaO48yXytm0vcxsrRb7egBvJqKL0bkNnkxEf8XMv1tbpiYbI80FIlASv1qeOySOteUKkAohYYUDKZnat1b6kA6ngL2KhaWnSblAWDq4+DmlI25uWju+3pL4tNVBpESZi2etSSwyGFkshcbnaawT21o55DS89W3VGITKPLT7WWm4TyxhxqumWEylwGbmqwFcDQC9ZPT+IQ0R0Ghj5DjO7HF3kJFYLgOpQGmSMaFDrJmr+FhOIirRGclUkwIsqSa1GomcnbGM6OLeVT6XJQ1oOr4hoVGsHjv1rq0QrlI3JiWlFKnVXIassZZDC3drrRZsPS+wUp+l/R/MwgKbme8GcPfQ64+4xshxnGlwdxADTVrRwmbI7ZQOqXRF2RJJyZKINKkjN0MWSwOWQ6e1Lntc3lK3EC3UqsxfOotqAeBKZnDkvaXeSj5vyuXHCnJf4lAarpUS0qyYcnZNk+6thQXC88bhe63wuTEtuoM00xg5jjM/3FF2AnLrlWtL0YwJIGaFiCixk7EkIs3uRx6TuomaoHGpEKjhPjWOnlL/UxNmxLK90mb6Su2MNMlQ2uNICVaTjHJhiaegRIqU9kCaZCSlwVSo2lLFtCuwHcdZOC1KRqM6AyJ6LxF9j4juJaKbiGgdEZ1DRPcQ0Q+I6HNEpEVgcBxngQR3kKm89qdgsGRERGcBeDeAVzDzz4noFgDvAHAxgI8z881E9D8AXAHg2lx+mrI6FuGt4VmJAnuI8WNueKYpsK2hVskKFdZ6WGvFdvy71IgzNeQKpIweLQO/mmGaFaM6vo+m1NbQnt+KDWQNnYGV0UOt+6SYcminGT3K57LcQrSVWXLGj0eVZISuMXsREa1Gt0rMbgBvBHBrf/wGAG8ZeQ/HcSYmzKaV/M2LwZIRMz9MRB8F8CCAnwP4CoAdAJ5i5gP9absAnDW6lLAV1iWrw+Za3Jp1sUqceS1pR4tnbcW6Dq4JmutD6fpv2hRv7uOS4TPifZYRXcosw5IoNQktFekzvq9UpMe/w3dh5RUrsoMEMWQl49LtmNL7aJMV8r3INfKewxIlYUWOKp0REZ0K4FIA5wA4E90Kyb+pnMrG9VuJaDsRbW9Nq+84xwIHC//mxZjZtDcB+DEzPw4ARHQbgNcBOIWIVvfS0WYAj2gX9yEItgHAaiI+iJVT03FFhJVCcxJRSYxjywlW3jPezoUFSZ1jGTbG+3IhUoZIfBpSWtL0c/F52v1KpAGrHlPuNLVB1TS3CSsch5a3JSXWuBhZ31hq1sZy30mVxXpv2tS+vMZyB2lNCBijM3oQwAVEdAIREYALAdwH4C4Ab+3PuRzA7eOK6DjOLDhqJCNmvoeIbgXwTQAHAHwLnaTzvwDcTER/3O+7rjZvrUW3QifUGKzJfLWKLtVxaJKRZdxo6Y7ifSXuLbWkjAjDfaTzq6YzkoaSNXqSUglJ21ez0oYsv9QdaSFkA1IHVTNbadXF1P/ElvFjeM44/Inl8hPTomQ0yuiRma8BcI3Y/SMArxmTr+M4s4Wx3BOgBZq0wNZ6Z9kj1PRicjYoFSAN4pycRJRy7ZCSUCpsRu65Uu4gVtlTvbQlEWm6FVkG2eOWTP/mQrHEv2skooClW7EWL9D2lTxHTnIdMzOX2p+z9YqfM0hJzyjHYo4qychxnCOTFqf2m26MUlbHJRKR7E2snqCkd7YcMEuWAbIss+OyWT1uytpY5iEZ8rFpluU5e59UOawyaPtzYVqscsTnyHctr42lBLnPClFbEpZGHk+Rcoy1kJKf1Blp4VS0dxnwxshxnGbwYZpBbGcUSCnYxswyWTM8QF4isvRC8T4pEclrU7OEU4Q9kZT0gCm9mnwv1izNENeBkhnNIVgSjBbsXt5XvseS8C3Wcc0WSr7bIRKS1KvGdk3rxDGtjj24muM4TdDiMG0W8aQcxzkCmMrokYjOJqK7iGhnH1LoyiHlaUoykg8+lRiZU6Km1vuqmdovWScNSIv9UxgWSuJrrSn21JDBGi7lokbWkptgSE35W2YblumCdkwaSqYccnPT8iVDrnC/IcO11Cq00olWi3I5sdHjAQDvY+ZvEtFJAHYQ0VeZ+b6aTJpqjBzHmR9TDdOYeTe68EFg5meIaCe6aB1HXmMkW+lQSSdF+6zVTKdcAw2ol4hKQmBYBmtAfpq+pJfMoUkDuQ+xptec9Vi/xCSjNG62JvHlDCXjd5yTiFJuPLIM8h1Y0p2Wn/ymUiulPIOVzEpnRERbAJwH4J7aa5tojBzHmS+Vs2kbiWh7tL2tj7qxDCI6EcAXALyHmZ+uLVMzjVGsLEuFRagZk1stvzV9r+3LSUSa0aMmNQHpFT7HSHySMcaOQySxIW4bs8q/5v6WUWzq+5DmJjUO20OCq1nIMsfPJZ1nrTqvkH73MvP5qROIaA26huhGZr6tPOslmmmMHMeZH1MO0/oQQtcB2MnMHxuaT1ONUWjBtTXQUlJFjBbMvDS4PmC7feRCi2jnBGqcXsfoiCwH2ZKZsdLjtaQC/Utyz17i4FxDzk1Iu0+NtCv3D3m31jdsnQekXUViJpQ4Xw/gnQC+S0Tf7vf9ETN/qSaTphojx3Hmw5RT+8z8dwBobD5NNUap8becTSvpZYYsNyT35Zw3Uw6lgRpd0SxYpKWtdW+tLuT7qgnPWmprpeUhZ7FSNjyyLJZrzJAAbVN9C3KGbQ6S0SQ01Rg5jjMfDsF900xWwQ6pCdh6l1QIh5wYmuoZSnVFKUfPnDOlRip0iEUueFzq3mN6x1T95nr5ktm73DlDyq6Vy7pPSvrN/SOnnmsWtmSa7ZpLRo7jNM9RFwPbcZwjF5eMEsjhWckaaEMoiWdUo7iW1wZKDBlLh2Mlyk05XNPIGYKWUPIRjzGiHEKuTEOGrtpEilwPrSTaZe5YTR2VuvPE52r5txhCpKnGyHGc+eDB1QwIXe+TmpIsVeymVvqQx1IrVJQqrrVrc/G6S5woS8gpXmsY0kuWKLAX1fuO0YeU1KcVmsRyFxlLrUFofG9XYDuO0yyuwDYgdGPx4NgnU0DXI8XbJRJRzsUjdW6JE6VldpBaUWIKPcKUU+BTfaCtfOhTuJBo0k0u4NzY/Eup0R0NuXaeNNEYOY4zX1wyMiB041qpK0o5yg7pTYY4u5ZIU7kyppxhh+iKLKaUgGbda5Y4qS5K71Sz9llAllGTkKaQgEpWhy1tZFwychxn4fhsWoLjsFIy0uyMrN5F0x1ZDrH7xHbJ2mcy1fQ0cjYwp0PSKJklLO1hS3rIMYHYSmhZXzGmDDmXDk2Ctv75S96nlIhSbj0p6T3gdkaO4zSBN0YGYTZN6opKdCspa2q5T/YY2nJDpSFEtEDopc68MUN0RWMUj2MslWdxP40aK2N5zRTlGFK/lm3S2Pq0wuBIS/CUrdwUVvfzoInGyHGc+eKSkeM4zeCSUYIhQ5ZwjaYcHDJNn1sNxFJSx/usKf0SReUYB9Oanm5WcYxKGeLMK9FW+qhhlpJBKm6S9ewl7zx8YymnaPkNa++LsdJ1ZdFk/z+I6DNEtIeI7o32bSCirxLRD/r01H4/EdF/J6IfEtF3iOjVsyy84zjDCEaPJX/zokQyuh7ApwB8Ntp3FYA7mfnDRHRVv/1BAL8J4Nz+79cAXNunVWg9liaJaOdqUfmk64ZM47XIayWiuDWf0tl1il57EW4grd4vRcmKtTlqVkGxVo9J2f1YEygpVxVphmLl2QpZyYiZvwbgCbH7UgA39L9vAPCWaP9nuePrAE4hok1TFdZxnGkICuySv3kxVGf0EmbeDQDMvJuIzuj3nwXgoei8Xf2+3bkMD2KlFPJsdPyUPrVWjNCkHCndPC+2U5KRNYVfM20vSY3zZ/nSx0ohrfWgpeud1TJGgrCkHS2fEnMUeWytsZ26Rlv1OKYl6RSYXoGtrZ3E6olEWwFsBeYXDdBxnI6jyR3kMSLa1EtFmwDs6ffvAnB2dN5mAI9oGTDzNgDbAICIHn8AePYBYO/A8sybjfCyzgIv62wIZf3FsOMQ8OVnu/0lzOU5hzZGdwC4HMCH+/T2aP8fENHN6BTXPw3DuRTMfDoRbWfm8weWZ654WWeDl3U2aGVl5osWVR6LbGNERDcBeAOAjUS0C8A16BqhW4joCgAPAnhbf/qXAFwM4IcAngPwezMos+M4RyHZxoiZLzMOXaicywDeNbZQjuMce8xjqfdSti26ABV4WWeDl3U2HBFlpU6YcRzHWSwtSUaO4xzDNNEYEdFFRHR/79N21aLLE0NEZxPRXUS0k4i+R0RX9vtV/7xFQ0SriOhbRPTFfvscIrqnL+fniChlLzdXiOgUIrqViL7f1+9rG67X9/bv/14iuomI1rVSt0eL/+jCGyMiWgXgz9D5tb0CwGVE9IrFlmoZBwC8j5l/FcAFAN7Vly/4550L4M5+uwWuBLAz2v4IgI/35XwSwBULKZXOJwH8DTP/CoBXoSt3c/VKRGcBeDeA85n5lejsdN+Bdur2egByqt6qx9h/dCs6/9E2YOaF/gF4LYAvR9tXA7h60eVKlPd2AL8B4H4Am/p9mwDc30DZNqP78N4I4IvoLOL3Alit1fWCy3oygB+j11tG+1us1+DmtAHdDPQXAfzrluoWwBYA9+bqEcCnAVymnbfov4VLRrD92ZqDiLYAOA/APRD+eQDOsK+cG58A8AEsuR2dBuApZj7Qb7dUty8D8DiAv+iHlX9OROvRYL0y88MAPorOpm43gJ8C2IF26xaw67HZ/7cWGqNif7ZFQkQnAvgCgPcw89OLLo+EiC4BsIeZd8S7lVNbqdvVAF4N4FpmPg+dX/TCh2Qavb7lUgDnADgTwHp0wx1JK3WbotlvooXGqNifbVEQ0Rp0DdGNzHxbv/uxEB5F+OctitcDeDMR/QTAzeiGap9AF8YlGLe2VLe7AOxi5nv67VvRNU6t1SsAvAnAj5n5cWbeD+A2AK9Du3UL2PXY7P9bC43RNwCc289MrEWnGLxjwWU6DBERgOsA7GTmj0WHgn8esNw/byEw89XMvJmZt6Crw79l5t8BcBeAt/anLbycAWZ+FMBDRPTyfteFAO5DY/Xa8yCAC4johP57CGVtsm57rHq8A8C/62fVLkCh/+hcWLTSqleiXQzgHwD8PwAfWnR5RNl+HZ0Y+x0A3+7/Lkanj7kTwA/6dMOiyxqV+Q0Avtj/fhmA/4POX/DzAI5fdPmicv5zANv7uv2fAE5ttV4B/BcA3wdwL4C/BHB8K3UL4CZ0uqz96CSfK6x6RDdM+7P+f+276GYIF16/zOwW2I7jtEELwzTHcRxvjBzHaQNvjBzHaQJvjBzHaQJvjBzHaQJvjBzHaQJvjBzHaQJvjBzHaYL/D50BFDPhTnzOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#NBVAL_SKIP\n",
    "\n",
    "# Plot percentage error\n",
    "plot_image(100*np.abs(model0.vp-get_true_model().vp.data)/get_true_model().vp.data, vmax=15, cmap=\"hot\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEaCAYAAADDgSq4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8FfXZ///XlYWENewECAICIhB2RMC1dQNbFWVz62K93brZ2m9729a7/mpr79q79VZb96W1rYrgWr3LZt1F0bCHVUA0LGHft0By/f6YAWOawAnJyeSc834+HudBzpw5M9eJMe/MzGeuj7k7IiIi8ZAWdQEiIpK8FDIiIhI3ChkREYkbhYyIiMSNQkZEROJGISMiInGjkBERkbhRyIgAZnalmRWY2W4zW29mU8zs9KjrEkl0ChlJeWZ2C3AP8BugHXAC8ABwSZR1HWZmGVHXIHK8FDKS0swsB7gD+I67v+Due9z9oLu/4u4/NrMsM7vHzNaFj3vMLCt879lmtsbMfmRmG8MjoGvC14aZWbGZpZfb16VmtiD8Os3MbjWzlWa2xcwmmVnL8LUuZuZmdq2ZfQa8Hi7/upl9Gq7/X2a22szOrcb2vmFmn5nZZjP7ebm60s3sZ+F7d5nZbDPrFL52spnNMLOtZrbMzMbXwX8WSSIKGUl1w4Fs4MUqXv85MAwYAPQHhgK3lXs9F8gBOgLXAvebWQt3/wDYA3y53LpXAk+HX38fGA2cBXQAtgH3V9j3WUAv4AIz601wdHUV0L7cPg+LZXunAz2Bc4BfmFmvcPktwBXAhUAz4FvAXjNrDMwIa24brvOAmfWp4nsl8u/cXQ89UvZB8Eu7+CivrwQuLPf8AmB1+PXZwD4go9zrG4Fh4de/Bp4Iv25KEDqdw+dLgHPKva89cBDIALoADpxY7vVfAM+Ue94IKAHOrcb28sq9/iFwefj1MuCSSj77BOCdCsseBm6P+r+bHonz0LleSXVbgNZmluHuhyp5vQPwabnnn4bLjry/wvv2Ak3Cr58GZprZTcBlwBx3P7ytzsCLZlZW7r2lBNeEDiuqUMeR5+6+18y2lHs9lu0VV1FnJ4IwragzcKqZbS+3LAP4WyXrilRKp8sk1b0P7Cc41VSZdQS/bA87IVx2TO6+mCCURvHFU2UQBMYod29e7pHt7mvLb6Lc1+uBvMNPzKwh0Kqa26tKEdCtiuVvVdhmE3e/KYZtigAKGUlx7r6D4FTU/WY22swamVmmmY0ys98BzwC3mVkbM2sdrvv3auziaYLrJWcCk8stfwi408w6A4TbP9potueAi8xshJk1AH4JWA22V95jwK/MrIcF+plZK+BV4CQz+1r4Pck0s1PKXcsROSaFjKQ8d7+b4OL3bcAmgr/gvwu8RHBdpQBYACwE5oTLYvUMwbWb1919c7nl9wL/AKab2S7gA+DUo9S4CPgeMJHgqGYXwfWfA8ezvQruBiYB04GdwONAQ3ffBZwPXE5w9FYM3AVkxbhdEcxdk5aJJBozawJsB3q4+ydR1yNSFR3JiCQIM7soPJ3XGPg9wZHV6mirEjk6hYxI4riE4LTVOqAHwRBknYqQek2ny0REJG50JCMiInGjkBERkbhJ+Tv+W7du7V26dIm6DBGRhDJ79uzN7t7mWOulfMh06dKFgoKCqMsQEUkoZvbpsddK0pAJh3g+QNBA8E13fyrikkREUlJcr8mY2Q/NbJGZFZrZM2aWfZzbeSKcr6OwktdGhvNcrDCzW8PFlwHPuft1wMU1+AgiIlIDcQsZM+tI0LNpiLvnA+kE7SnKr9PWzJpWWNa9ks39BRhZyT7SCebMGAX0Bq4I593I4/OOtaU1+yQiInK84j26LANoGE4f24h/7157FvDy4SMcM7sOuK/iRtz9bWBrJdsfCqxw91XuXkLQ1+kSYA2fd6zVCDoRkYjE7Rdw2GL898BnBA39drj79ArrTAamAhPN7CqCGfmqM71rR74458aacNkLwBgzexB4pbI3hi06HtmxY0c1diciItURz9NlLQiOKroSTLjU2Myurrieu/+OYD6PB4GL3X13dXZTyTL3YJ72a9z9pqou+nswh/v1OTk51didiIhURzxPJZ0LfOLum9z9IMHRxYiKK5nZGUA+wRzrt1dzH2sIZvU7LI8YJ5QSEZH4i2fIfAYMC7vGGnAOwTzkR5jZQOBRgiOea4CWZladuTo+AnqYWddwIqfLCebUEBGReiCe12RmEczmN4egJXka8EiF1RoB49x9pbuXAd/gi/OpA2BmzxBMk9vTzNaY2bXhPg4RTC41jSDAJoWTO4mISD2Q8l2YhwwZ4rrjX0SkesxstrsPOdZ6Gt4rIiJxo5AREZG4ScreZdWxc99Bpi8qPq73Ns3OZHDnFjTIUFaLiFQm5UPm0617uf5vs4/7/c2yMzi3VztG5udy5kltyM5Mr8XqREQSW8qHTPe2TXj2e6cf13vX79jPtEXFzFi8gRfmrqVRg3S+1LMtI/Nz+dLJbWmSlfLfXhFJcRpdVgujyw6WljFr1VamFK5n2qINbN59gAYZaZzZozUj89tzbq+2NG/UoJYqFhGJXqyjyxQytTyEubTMmfPZNqYsLGbaomLWbt9HRpoxvFsrRubncl7vdrRtelwzHoiI1BsKmRjF8z4Zd2fh2h1MKSxmamExn2zegxmc0rklF+TnMjI/l47NG8Zl3yIi8aSQiVFd3Yzp7izfsJuphcVMKVzP0uJdAPTLy2Fkfi6j8tvTtXXjuNchIlIbFDIxiuqO/9Wb9zB1UTFTCouZX7QdgJ7tmgaB0zeXnu2aErR8ExGpfxQyMaoPbWXWbd/HtDBwPlq9FXfo0qoRI/PbMzI/l/55OQocEalXFDIxqg8hU96mXQeYsXgDUwrX8/7KLRwqczrkZAfXcPrkMqRLS9LTFDgiEi2FTIzqW8iUt2PvQV5bsoGpi4p5a/kmSg6V0bpJA87rncuo/FyGd2tFZrq6DYhI3VPIxKg+h0x5ew4c4o1lG5laWMwbSzeyp6Q06DbQux2j8ttzRo/W6jYgInVGIROjRAmZ8vYfLOWdjzcztbCY15ZsYMe+g0G3gZPbMio/l7N7qtuAiMRXrCGj30QJKDsznfN6t+O83u04WFrG+yu3MHVRMdMXFfN/C9aH3QbaMCo/l3N7tSOnUWbUJYtIitKRTAIeyVSltMyZ/em2oL1NYTHrduz/QreB83vn0qZpVtRlikgS0OmyGCVTyJTn7ixYc7jbwHpWb9kbdBvo0pKRfYJuAx3UbUBEjpNCJkbJGjLluTvLNuw60k/tcLeB/nk5jMxvz6j8XLqo24CIVINCJkapEDIVrdq0m6mLiplWWMz8NTsAODm36ZH2Nie1a6KbP0XkqBQyMUrFkClv7fZ9TAsbeH70adBtoGvrxowMb/7sp24DIlIJhUyMUj1kytu06wDTFweBU7HbwKj89gzu3ELdBkQEUMjETCFTue17S3htyUamFq7n7Y83h90Gsji/TztG5ecy7ER1GxBJZQqZGClkjm33gUO8sTTsNrBsI3tLSslpmMm5vYLAOV3dBkRSjkImRgqZ6tl/sJS3l29iamExM5ZsYNf+QzQ+0m2gPWf3bENjdRsQSXq641/iIjsznfP75HJ+n1xKDpXx/qotQeAsLubVBevJykjjzJOCbgPnnKxuAyKpTkcyOpKpFaVlTsHqrUwpDO7FWR92GxjRvTUj++Ryfp92tG6ibgMiyUKny2KkkKl9ZWXOgrU7mFK4nqmFxXy6ZS9ph7sN5AfdBtrnqNuASCJTyMRIIRNf7s7S4l1H2tss37AbgP6dmjMqP5gXp3MrdRsQSTQKmRgpZOrWyk27mRqeUltQrtvAqPz2jOqbS4+26jYgkggUMjFSyERnzba9TFu0gamF6yn4dBvucOLhbgP5ufTtqG4DIvWVQiZGCpn6YePO/UxfvCHoNrBqC6VlTsfmDbmgTy6j+uYy+IQWpKnbgEi9oZCJkUKm/tm2p4TXlgSB887HmykpLaNN0yzOD6eaPvXEluo2IBIxhUyMFDL12679B3lj2SamFq7njaWb2HewlOaNvthtICtD3QZE6ppCJkYKmcSxr6SUtz8Oug28FnYbaJKVEXYbyOWsk9RtQKSu6I5/SToNG6RzQZ9cLgi7DcxcuZlpi4qZvmgDr8xfR1ZGGmed1IaR+bmc06sdOQ3VbUAkajqS0ZFMwjtUWsZHq7cxbVEwTUHxzv1kphsjurVmZH4u5/duRyt1GxCpVTpdFiOFTHIpK3PmrdnOtMJiphQW89nWoNvA0K4tGdknlwvUbUCkVihkYqSQSV7uzpL1u5hauJ4phcV8vDHoNjDgSLeB9pzQqlHEVYokJoVMjBQyqWPFxt1HTqktXBt0G+jVvtmR9jbd1W1AJGYKmRgpZFJT0da9RwJn9mdht4E2jRmVn8vIPu3J79hMgSNyFAqZGClkZOPO/UHgLCrmg1Vbj3QbGBke4QxStwGRf6OQiZFCRsrbuqeE1xZvYOqiYt4t123ggj5ht4GuLclQtwERhUysFDJSlZ37D/LG0o1MLSzmzWWfdxs4r1c7RvXN5bTu6jYgqUshEyOFjMRiX0kpby0P2tv8a8lGdh0Iug18+XC3gZ5taNRA9zZL6tAd/yK1qGGD9CNTEJQcKuO9lZuZurCY6YuL+cf8dTTNzuCBqwZxRo82UZcqUq/oSEZHMlIDh0rL+HD1Vu54ZTErN+3mD+MHcHH/DlGXJRJ3sR7J6AqmSA1kpKcxoltrnr1hOAM7teD7z8zlz+99EnVZIvWGQkakFuQ0zOSv1w7l/N7t+OUri/mfaUtJ9bMEIqCQEak12ZnpPHDVIK4Y2on731jJrc8v5FBpWdRliURKF/5FalFGehq/ubQvbZpkcd/rK9iyp4Q/XjGQhg001FlSk45kRGqZmXHL+T2545I+/GvpBr72+Cx27D0YdVkikThqyJhZupmNqatiRJLJ14d34U9XDGLBmh2Me3gmxTv2R12SSJ07asi4eynwgzqqRSTpfKVfe/5yzSms276fMQ/OZEU43YBIqojldNk0M/uBmbU3s2aHH3GvTCRJjOjemonXD+PAoVLGPTSTuZ9ti7okkToTS8jcAPwI+BBYFD4K41mUSLLJ75jD8zeNoGl2Jlc+Oos3l22MuiSROnHMkHH3TpU8TqiL4kSSSedWjXn+phGc2KYx//FkAS/MWRN1SSJxd8yQMbMMM/u2mU0MHzeamYY+ixyHNk2zmHj9MIZ2bcktk+bz6Nuroi5JJK5iOV12PzACeCJ8jAAeiGdRIsmsaXYmf77mFC7sm8ud/1zCb/65hLIydQeQ5BTLEckwd+9f7vl0M5sfr4JEUkFWRjp/vGIQrRov4pG3V7F59wHuGtOPTE2IJkkmlpApM7Mu7r4awMy6AOqVIVJD6WnGHZf0oU3TLO6esZxte0q4/6pBmpdGkkosfzb9BHjbzF4zs38BbwE/jm9ZIqnBzPj+OT34zaV9eWv5Jq56bBbb9pREXZZIrTnqn0xmlgbsBHoCvQADFrv7vjqoTSRlXHnqCbRs3IDvT5zL2Idm8tdrT6Vj84ZRlyVSY8e6478MuNfd97n7HHefrYARiY+R+bn87VtD2bjrAGMemMnyDbuiLkmkxmI5XTbDzC6JeyUiwqkntmLSDcMpdWfsgzMpWL016pJEaiSWkPku8KKZ7TOzrWa2zcz0ky8SJ73aN+OFm0bQqkkWVz02i9cWb4i6JJHjdqwuzAb0BzKBJkAboHX4r4jESaeWjXjuxuH0zG3KDX+fzaSCoqhLEjkux7om48CL7l5a8VFH9YmkrFZNsnj6umGM6NaKnzy3gAfeXKEpnSXhxHK67EMzGxT3SkTk3zTJyuDxb5zCxf078Lupy7jj1cXqDiAJJZa7vk4HrjOzlcAegmHM7u4KHpE60CAjjXsmDKBVkwb8+b3VbNldwu/H9adBhroDSP0XS8iMjnsVInJUaWnGL77amzZNs/jd1GVs21vCg1cPpkmWugNI/RZLq/+VBBf6Twu/3g5ownKROmZmfPvs7vxubD9mrtzClY9+wJbdB6IuS+SoYmn1fxtwO3BbuCgbeDqeRYlI1cYP6cTDVw9mWfEuxj70PkVb90ZdkkiVYjmpOxa4kOB6DO6+FtD0yyIROrd3O576j1PZuqeEMQ/OZMn6nVGXJFKpWELmQDiU2QHMrFF8SxKRWAzp0pLJNw4nzYzxD7/PB6u2RF2SyL+JJWReMLP7gRwzuwaYTjB5mYhE7KR2TXn+2yNo2zSLrz/xIVMLi6MuSeQLYrnwfxfwKvAPgrv/73T3e+JdmIjEpmPzhjx34wh6t2/Gt5+azdOzPou6JJEjYhr/6O5TgClxrkVEjlOLxg14+rpT+fZTc/jZiwvZvPsA3/tyd4LOUCLR0d1cIkmiUYMMHv36EC4b1JG7Zyzn9n8solTdASRiupNLJIlkpqfxh3H9adMki4ffXsWW3SXcPaE/WRnpUZcmKSqmkDGzBsAJ7r4izvWISA2ZGT+9sBetm2Rx5z+XsG1vCQ9/bTBNszOjLk1SUCw3Y34FWAjMCJ8PMLMX412YiNTMdWeeyP9O6M+Hn2xlwsMfsHHX/qhLkhQUyzWZO4BTCdrJ4O7zgO7xLEpEaselA/N49BtD+GTzHsY++D6fbtkTdUmSYmIJmYPuvr3CMl1NFEkQX+rZlqevO5Vd+w8y5sGZFK7dEXVJkkJiCZklZjYeSDOzrmZ2D/BBnOsSkVo08IQWTL5xBFkZ6Vz+yAfMXLE56pIkRcQSMt8FBgNlwAvAfuAH8SxKRGpf97ZNeP6mEXRs3pBv/vkjXl2wLuqSJAXEEjInuvt/uvvA8HGru6vtq0gCys3JZtINw+mXl8P3npnLX99fHXVJkuRiCZkHzGyRmd1uZifHvSIRiaucRpn8/T9O5ZyT2/GLlxfxh+nLCHrgitS+WHqXnQGcD+wCnjSzuWZ2a9wrE5G4yc5M56GrBzF+SB5/fH0FP3txIYdKy6IuS5JQTG1l3H2tu98NfJPgnplfxbOomjKzxmb2pJk9amZXRV2PSH2UkZ7GXWP68Z0vdeOZD4v49lNz2H+wNOqyJMnEcjNmDzO7zczmAY8CHwEnxPC+nmY2r9xjp5kd14ABM3vCzDaaWWElr400s2VmtqLcEdZlwHPufh1w8fHsUyQVmBk/vuBkbr+oNzOWbODrT3zIjn2aXV1qTyxHMk8D+4CL3f10d/+ju68/1pvcfZm7D3D3AQSj0/YCX+gUYGZtzaxphWWV3ej5F2BkxYVmlg7cD4wCegNXmFlvIA8oClfTn2Yix3DNaV259/KBzP1sGxMefp8NO9UdQGpHLNdkTnH3P7h7TSapOAdY6e6fVlh+FvCymWUDmNl1wH2V1PA2sLWS7Q4FVrj7KncvASYClwBrCIIGqviMZnaRmT2yY4duTBMBuLh/B/78zaEUbd3LZQ/MZNWm3VGXJEmgypAxs2fCf+ea2Zxyj7lmNqea+7kceKbiQnefDEwFJobXTr4FjK/Gdjvy+RELBOHSkeB+njFm9iDwSmVvdPdX3P36nJycauxOJLmd3qM1z1w/jP0HSxn70PvML6rY7EOkeo7WhfnH4b9ja7KDsIPzxcBPK3vd3X9nZhOBB4Fu7l6dP58qm5HJ3X0PcE21ixUR+uU157mbRvC1x2dxxaMf8NDVgznzpDZRlyUJqsojGXdfE355rbuvLP8Arq3GPkYBc9x9Q2UvmtkZQD7B9Zrbq7FdCI5cOpV7ngfoNmaRGuraujEv3DSCE1o24lt/+YiX562NuiRJULFc+P+3C+7AV6qxjyuo5FQZgJkNJBixdgnBkUdLM/t1Nbb9EdAj7KnWgOC03D+q8X4RqULbZtlMunE4gzu34OaJ83ji3U+iLkkS0NGuydxgZnOBnhWuyXwMLIll42bWCDiP4BpJZRoB48IjpDLgG0DFwQGHrw+9H9ayxsyuBXD3QwS91aaFNU1y90Wx1CYix9YsO5MnvzWUkX1yuePVxdw1dam6A0i1WFU/MGbWAmgF/DdQ/g7/Xe6+sQ5qqxNDhgzxgoKCqMsQqddKy5z/ermQp2d9xtjBefz2sr5kpMd0L7ckKTOb7e5DjrVelRf+3X0bsA0YF26wJZANZJhZB3fXtQ+RFJGeZtw5Op82TbK4918fs21PCX+6chANG6RHXZrUc7Hc8X+hmS0nuMg+i2DI8OvxLkxE6hcz44fnncSvRufz+rKNXP34LLbvLYm6LKnnYjne/Q1wGrDM3TsRDAR4M55FiUj99bVhnbn/ykEsXLODcQ+9z/od+6IuSeqxWELmkLtvIpgZ09x9BjAoznWJSD12Yd/2PPmtoRTv2M+YB2ayYuOuqEuSeiqWkNlhZo2Bd4G/mtkfCGbJFJEUNrxbKybeMIySUmfsQ+8z+9NtUZck9VAsITMaOEAw5fKbwFrgojjWJCIJok+HHF64aQQ5DTO56rEPeH1ppfdcSwqLpUHmLnc/5O4H3f1xd787PH0mIsIJrRrx3I0j6N62Cdf9dTbPz15z7DdJyjjazZjbzGxruce28v/WZZEiUr+1aZrFM9cNY9iJLfnR5Pk8/NbKqEuSeuJoDTJb11kVIpLwmmZn8sQ3T+GWSfP57ylL2bz7AD8d1Yu0tMr62EqqONrNmEcm+zKzYcBJ7v7X8KbMJkBN5pcRkSSUlZHOHy8fSJsmWTz6zids3l3C78b2I1PdAVLW0Y5kADCz2wjuk+kG/BVoSDBb5unxLU1EElFamnH7Rb1p0zSL/5m2jK17SnjgqkE0zjrmrxtJQrH8eTEWuBDYA+Dua4Fm8SxKRBKbmfGdL3Xnt5f15Z2PN3HlY7PYukfdAVJRLCFzwIMumg5HOiuLiBzT5UNP4KGrB7N0/U7GPjSTNdv2Rl2S1LFYQuYFM7sfyDGza4DpwBPxLUtEksX5fXL527WnsmnXAcY8OJNlxeoOkEpiuU/mLuBVgsnA+gN3uvs98S5MRJLH0K4tmXzjcADGPTSTj1brLohUcdSQMbN0M5vm7lPc/Yfu/gN3n1JXxYlI8jg5txnP3zSC1k2zuPqxWUwuKKLkkDpUJbujhkw4jLnEzHShX0RqLK9F0B2gV/tm/Pi5BQz/73/x61cXs3yDTqElq1jGFO4G5pvZdMIRZgDufkvcqhKRpNWycQOeu3E473y8mUkFRTz5/moee/cT+ndqzoQhnfhq//Y0y86MukypJVVOv3xkBbNrK1vu7o/HpaI6pumXRaK1ZfcBXpy7lkkFRSzfsJvszDQuzG/P+FM6cWrXlpipY0B9FOv0y8cMmWSnkBGpH9yd+Wt2MKmgiFfmrWPXgUN0btWIcYPzGDM4j/Y5DaMuUcpRyMRIISNS/+wrKWVK4XomFRTxwaqtpBmceVIbxg/pxLm92tEgQ21qoqaQiZFCRqR++3TLHiYXrOG52Wso3rmflo0bMHpAR8afksfJuRqTFJVaDxkzy3L3AzWurJ5RyIgkhtIy552PNzG5YA3TFxdzsNTpl5fD+CGduKh/B3IaarBAXaq1kDGzocDjQI67n2Bm/YH/cPfv1U6p0VLIiCSerXtKeCkcLLC0eBdZGWlc2Lc944bkMaxrK00vUAdqM2Q+ACYAL7n7wHBZobvn10qlEVPIiCQud2fh2mCwwMvz1rFr/yFOaPn5YIEOzTVYIF5qM2Q+dPehZja3XMjMd/f+tVRrpBQyIslh/8FSphYWM6mgiJkrt2AGZ/Row4QhnTi3d1uyMtKjLjGpxBoysdyMWRSeMnMzSwe+ByyvaYEiIrUpOzOd0QM7MnpgR4q27mVyQRHPzV7Dd56eQ/NGmYwe0JEJp3SiV3sNFqhLsRzJtAXuA84NF70GfNfdN8e5tjqhIxmR5FVa5ry3YjPPFhQxY9EGSkrL6Nsxh/FD8rh4QEcNFqgBDWGOkUJGJDVs21PCy/PW8mzBGpas30lWRhoj83MZP6QTw0/UYIHqqs1rMo8STlhWnrtff/zl1R8KGZHU4u4sWreTSQVFvDR3LTv3HyKvRUPGDe7E2CF5dNRggZjUZshMKPc0G7gUKNIQZhFJdPsPljJtUTGTC9bw7orNmMHp3VszfkgnzuvdjuxMDRaoStxOl5lZGjDD3c853uLqE4WMiAAUbd3Lc7ODzgJrt+8jp2Emlw7syLghefTpkBN1efVOPEOmGzDN3bsfb3H1iUJGRMorK3PeW7mZSQVrmLaomJJDZfTp0IwJp3Tikv4dyWmkwQJQu6fLtvH5NZk0YCtwq7tPqnGV9YBCRkSqsn1vCS/PW8ekgiIWrdtJg4w0LuiTy4QhnRjRLbUHC9RKyFgwkUMnYG24qMyTbDiaQkZEYlG4dgeTC4p4ad46duw7SMfmDRk7OI9xQ/LIa9Eo6vLqXG0eycx298G1Vlk9o5ARkerYf7CUGYs3MKmgiHdXBLcLntatNeOG5HFBn9yUGSxQm3f8f2hmg9x9Ti3UJSKS0LIz07mofwcu6t+BNdv28vzstUyeXcTNE+fRLDuDMYPz+M+RJ6dM2BxLlSFjZhnufgg4HbjOzFYCewAD3N0H1VGNIiL1Ul6LRtx8bg++9+XuvL9qC89+VMSf31vNoVLnV6OToodwjR3tSOZDYBAwuo5qERFJSGlpxmndW3Na99a0bZrFY+9+wlknteHc3u2iLi1yR5vD1ADcfWVljzqqT0Qkofx4ZE96t2/GT55fwMad+6MuJ3JHO5JpY2a3VPWiu98dh3pERBJaVkY6910xkK/+8R1+NHk+T14zNKWHOh/tSCYdaAI0reIhIiKV6N62Cb/4ah/e+XgzT7z3SdTlROpoRzLr3f2OOqtERCSJXDG0E28t38hdU5cy7MRW5HdMzdY0x7wmIyIi1Wdm/PayfrRs3ICbJ85lX0lp1CVF4mghkxQNMEVEotKicQPuHj+AVZv38Kv/Wxx1OZGoMmTcfWtdFiIikoxO696a6888kadnfcbUwuKoy6lzRzuSERGRWvCj83rSt2MOt76wgOIdqTWsWSEjIhJnDTLSuPfyARw4WMYtk+ZRVpZUfYaPSiEjIlIHTmzThP/v4t7MXLmFR95ZFXU5dUYhIyJSR8YP6cT9fO0PAAALtklEQVSFfXP5/bRlLFizPepy6oRCRkSkjpgZ/31pP9o0zeLmifPYc+BQ1CXFnUJGRKQO5TTK5H8nDGD1lj3c8UryD2tWyIiI1LFhJ7bi22d349mCIv65cH3U5cSVQkZEJAI/OPck+ndqzq3PL2Dd9n1RlxM3ChkRkQhkpqdx74QBlJY5P3x2HqVJOqxZISMiEpEurRvzy0vymfXJVh56Kzmn6VLIiIhEaMygjny1X3v+d8Zy5hUl37BmhYyISITMjDsv7Uu7ZtncPHEuu5NsWLNCRkQkYjkNM7nn8gEUbd3L7S8virqcWqWQERGpB07p0pLvfqk7z89Zwz/mr4u6nFqjkBERqSe+f04PBp3QnJ+/uJA12/ZGXU6tUMiIiNQTGelp3Hv5QNzhh8/O41BpWdQl1ZhCRkSkHunUshG/Hp3PR6u38cCbiT+sWSEjIlLPjB7YkdEDOnDvvz5m9qfboi6nRhQyIiL10B2j82mfEwxr3rn/YNTlHDeFjIhIPdQsO5N7Lx/I+h37+cVLhVGXc9wUMiIi9dTgzi34/pd78NK8dbw0d23U5RwXhYyISD32nS91Y0jnFtz2UiGfbUm8Yc0KGRGReiwjPY17Lh+AGfzg2bkJN6xZISMiUs/ltWjEnZf2Zc5n27nv9RVRl1MtChkRkQRwcf8OjBmUx59e/5iPVm+NupyYKWRERBLELy/pQ6eWjfjBxHns2JcYw5oVMiIiCaJJVgb3TBhA8c79/PzFhbjX/9k0FTIiIglk4AktuOW8k3h1wXqen1P/hzUrZEREEsyNZ3Xj1K4tuf3lQlZv3hN1OUelkBERSTDpacb/ThhAeppx87PzOFiPhzUrZEREElCH5g357Zh+zC/azj2vLY+6nCopZEREEtSFfdszfkgeD7y5kg9WbYm6nEopZEREEtjtF/WhS6vG/PDZeWzfWxJ1Of9GISMiksAaZ2Vw7+UD2LTrAD+rh8OaFTIiIgmuX15zfnR+T/65sJjJBWuiLucLFDIiIknghjNPZES3Vtz+j0Ws2rQ76nKOUMiIiCSBtDTj7vEDyMpM4+aJ8yg5VD+GNStkRESSRG5ONr+9rB8L1+7g7hn1Y1izQkZEJImMzM/liqEn8PDbK5m5YnPU5ShkRESSzX99tRddWzfmh5PmsW1PtMOaFTIiIkmmUYMM7rt8IFv3lPCfzy+IdFizQkZEJAnld8zhJxeczPTFG3jmw6LI6lDIiIgkqWtP78oZPVpzx6uLWLExmmHNChkRkSSVlmb8YVx/GjXI4PvPzOXAodK6r6HO9ygiInWmbbNs7hrTj8Xrd/L7acvqfP8KGRGRJHde73Z8bVhnHn3nE95evqlO962QERFJAT//Si96tG3CjybPZ8vuA3W2X4WMiEgKyM5M574rBrJj78E6HdaskBERSRG92jfj1lEn89qSjfz9g0/rZJ8KGRGRFHLNaV0466Q2/Pr/lrB8w664708hIyKSQsyM34/rT6vGDVhaHP+QyYj7HkREpF5p0zSL1//f2WRnpsd9XzqSERFJQXURMKCQERGROFLIiIhI3ChkREQkbhQyIiISNwoZERGJG4WMiIjEjUJGRETixqKc+7k+MLNNwPE28WkNbK7FciRx5QA7oi4iwSTr9yyRPldNau3s7m2OtVLKh0xNmFmBuw+Jug6Jnpk94u7XR11HIknW71kifa66qFWny0RqxytRF5CAkvV7lkifK+616kimBnQkIyJydDqSqZlHoi5ARKQ+05GMiIjEjY5kREQkbhQyIiISN5q0TCRCZtYYeAAoAd5096ciLqleS+bvV7J+Nh3J1CIza2xmT5rZo2Z2VdT1SGzMrJOZvWFmS8xskZndXINtPWFmG82ssJLXRprZMjNbYWa3hosvA55z9+uAi493v3XJzLLN7EMzmx9+v35Zg23Vy++XmaWb2Vwze7UG26iXn62uKWSOoaoflFT6IUkBh4AfuXsvYBjwHTPrXX4FM2trZk0rLOteybb+AoysuNDM0oH7gVFAb+CKcB95QFG4WmkNP0ddOQB82d37AwOAkWY2rPwKSfD9uhlYUtkLSfDZ6pRC5tj+QoUflFT7IUl27r7e3eeEX+8i+OXSscJqZwEvm1k2gJldB9xXybbeBrZWspuhwAp3X+XuJcBE4BJgDcHPDSTI/48e2B0+zQwfFYepJuz3y8zygK8Aj1WxSsJ+tigkzQeJlyp+UFLqhySVmFkXYCAwq/xyd58MTAUmhqdCvwWMr8amO/L5HyAQ/Kx0BF4AxpjZgyTQneLh6aR5wEZghrsn0/frHuAnQFllLyb4Z6tzuvB/fCr7ITmV4K+ZP5nZV0iiH5JUYWZNgOeBH7j7zoqvu/vvzGwi8CDQrdxf8zFtvpJl7u57gGuOq+AIuXspMMDMmgMvmlm+uxdWWCfhvl9m9lVgo7vPNrOzq1ovET9bVPTX9vGp8ofE3a9x95uSZWRIqjCzTIKAecrdX6hinTOAfOBF4PZq7mIN0Knc8zxg3XGUWq+4+3bgTSq/9pCI36/TgIvNbDXBGYovm9nfK66UoJ8tEgqZ45NSPyTJzswMeBxY4u53V7HOQOBRgtOi1wAtzezX1djNR0APM+tqZg2Ay4F/1KzyaJhZm/AIBjNrCJwLLK2wTkJ+v9z9p+6e5+5dwn2+7u5Xl18nUT9bVBQyxyelfkhSwGnA1wj+ap0XPi6ssE4jYJy7r3T3MuAbVDIPkZk9A7wP9DSzNWZ2LYC7HwK+C0wjGFgwyd0Xxe8jxVV74A0zW0Dw/8IMd6841DeZv1/J/NlqnXqXHUP4g3I2wQRlG4Db3f3x8JfQPUA68IS73xldlSIi9ZNCRkRE4kany0REJG4UMiIiEjcKGRERiRuFjIiIxI1CRkRE4kYhIyIicaOQkZRkZrvDf7uY2ZW1vO2fVXg+sza3X9vM7Jtm9qeo65DkpJCRVNcFqFbIhFM9HM0XQsbdR1SzpoQSw/dDUphCRlLdb4EzwlYyPwxb2P+PmX1kZgvM7AYAMzvbgtkznwYWhsteMrPZFswOeX247LdAw3B7T4XLDh81WbjtQjNbaGYTym37TTN7zsyWmtlTYT+1LwjXucuCWSmXh00a/+1IxMxePdxB2Mx2h++ZbWavmdnQcDurzKz8xHqdzGyqBRPx3V5uW1eH+5tnZg8fDpRwu3eY2SxgeG39x5Dko1b/kupuBf6fu38VIAyLHe5+ipllAe+Z2fRw3aFAvrt/Ej7/lrtvDZtEfmRmz7v7rWb2XXcfUMm+LiOYSbI/QZuij8zs7fC1gUAfgkar7xH0U3u3km1kuPvQsK3R7QTNKY+mMcF88f9pZi8CvwbOI5hs70k+77k3lKCr8N6wrv8D9gATgNPc/aCZPQBcBfw13G6hu//iGPuXFKeQEfmi84F+ZjY2fJ4D9ABKgA/LBQzA983s0vDrTuF6W46y7dOBZ8K5WDaY2VvAKcDOcNtrACyYDKwLlYfM4WkIZofrHEsJwQRbEByBHQgDY2GF989w9y3h/l8Iaz0EDCYIHYCGBJOUQTDz6/Mx7F9SnEJG5IsM+J67T/vCwuD0054Kz88Fhrv7XjN7E8iOYdtVOVDu61Kq/n/zQCXrHOKLp77L13HQP29QWHb4/e5eZmbl91GxiaGH9T7p7j+tpI79YViKHJWuyUiq2wU0Lfd8GnCTBZOYYWYnmVnjSt6XA2wLA+ZkYFi51w4efn8FbwMTwus+bYAzgQ9r4TOsJpilMs3MOhGc+qqu88ysZXjqbzTBKbt/AWPNrC1A+HrnWqhXUoiOZCTVLQAOmdl84C/AvQSnkeaEF983EfzSrWgqcGM4p8oy4INyrz0CLDCzOe5+VbnlLxJcJJ9PcKTwE3cvDkOqJt4DPiE4HVYIzDmObbwL/A3oDjzt7gUAZnYbMN3M0oCDwHeoZO4Ukaqo1b+IiMSNTpeJiEjcKGRERCRuFDIiIhI3ChkREYkbhYyIiMSNQkZEROJGISMiInGjkBERkbj5/wG8MjbkF7HFoQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "distributed.core - WARNING - Event loop was unresponsive in Nanny for 89.32s.  This is often caused by long-running GIL-holding functions or moving large chunks of data. This can cause timeouts and instability.\n",
      "distributed.core - WARNING - Event loop was unresponsive in Nanny for 93.64s.  This is often caused by long-running GIL-holding functions or moving large chunks of data. This can cause timeouts and instability.\n",
      "distributed.core - WARNING - Event loop was unresponsive in Nanny for 93.65s.  This is often caused by long-running GIL-holding functions or moving large chunks of data. This can cause timeouts and instability.\n",
      "distributed.core - WARNING - Event loop was unresponsive in Nanny for 93.65s.  This is often caused by long-running GIL-holding functions or moving large chunks of data. This can cause timeouts and instability.\n",
      "distributed.core - WARNING - Event loop was unresponsive in Nanny for 93.65s.  This is often caused by long-running GIL-holding functions or moving large chunks of data. This can cause timeouts and instability.\n",
      "distributed.core - WARNING - Event loop was unresponsive in Scheduler for 93.65s.  This is often caused by long-running GIL-holding functions or moving large chunks of data. This can cause timeouts and instability.\n",
      "distributed.core - WARNING - Event loop was unresponsive in Nanny for 4.41s.  This is often caused by long-running GIL-holding functions or moving large chunks of data. This can cause timeouts and instability.\n",
      "distributed.core - WARNING - Event loop was unresponsive in Nanny for 4.50s.  This is often caused by long-running GIL-holding functions or moving large chunks of data. This can cause timeouts and instability.\n",
      "distributed.core - WARNING - Event loop was unresponsive in Nanny for 4.52s.  This is often caused by long-running GIL-holding functions or moving large chunks of data. This can cause timeouts and instability.\n",
      "distributed.core - WARNING - Event loop was unresponsive in Nanny for 4.57s.  This is often caused by long-running GIL-holding functions or moving large chunks of data. This can cause timeouts and instability.\n",
      "distributed.core - WARNING - Event loop was unresponsive in Nanny for 4.57s.  This is often caused by long-running GIL-holding functions or moving large chunks of data. This can cause timeouts and instability.\n",
      "distributed.core - WARNING - Event loop was unresponsive in Scheduler for 4.57s.  This is often caused by long-running GIL-holding functions or moving large chunks of data. This can cause timeouts and instability.\n",
      "distributed.core - WARNING - Event loop was unresponsive in Nanny for 4.56s.  This is often caused by long-running GIL-holding functions or moving large chunks of data. This can cause timeouts and instability.\n"
     ]
    }
   ],
   "source": [
    "#NBVAL_SKIP\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Plot objective function decrease\n",
    "plt.figure()\n",
    "plt.loglog(relative_error)\n",
    "plt.xlabel('Iteration number')\n",
    "plt.ylabel('True relative error')\n",
    "plt.title('Convergence')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<sup>This notebook is part of the tutorial \"Optimised Symbolic Finite Difference Computation with Devito\" presented at the Intel® HPC Developer Conference 2017.</sup>"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}