examples/seismic/tutorials/04_dask.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# 04 - Full waveform inversion with Devito and Dask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In this tutorial we show how [Devito](http://www.opesci.org/devito-public) and [scipy.optimize.minimize](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) are used with [Dask](https://dask.pydata.org/en/latest/#dask) to perform [full waveform inversion](https://www.slim.eos.ubc.ca/research/inversion) (FWI) on distributed memory parallel computers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## scipy.optimize.minimize \n",
    "\n",
    "In this tutorial we use [scipy.optimize.minimize](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) to solve the FWI gradient based minimization problem rather than the simple grdient decent algorithm in the previous tutorial.\n",
    "\n",
    "```python\n",
    "scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)\n",
    "```\n",
    "\n",
    "> Minimization of scalar function of one or more variables.\n",
    ">\n",
    "> In general, the optimization problems are of the form:\n",
    ">\n",
    "> minimize f(x) subject to\n",
    ">\n",
    "> g_i(x) >= 0,  i = 1,...,m\n",
    "> h_j(x)  = 0,  j = 1,...,p\n",
    "> where x is a vector of one or more variables. g_i(x) are the inequality constraints. h_j(x) are the equality constrains.\n",
    "\n",
    "[scipy.optimize.minimize](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) provides a wide variety of methods for solving minimization problems depending on the context. Here we are going to focus on using L-BFGS via [scipy.optimize.minimize(method=’L-BFGS-B’)](https://docs.scipy.org/doc/scipy/reference/optimize.minimize-lbfgsb.html#optimize-minimize-lbfgsb)\n",
    "\n",
    "```python\n",
    "scipy.optimize.minimize(fun, x0, args=(), method='L-BFGS-B', jac=None, bounds=None, tol=None, callback=None, options={'disp': None, 'maxls': 20, 'iprint': -1, 'gtol': 1e-05, 'eps': 1e-08, 'maxiter': 15000, 'ftol': 2.220446049250313e-09, 'maxcor': 10, 'maxfun': 15000})```\n",
    "\n",
    "The argument `fun` is a callable function that returns the misfit between the simulated and the observed data. If `jac` is a Boolean and is `True`, `fun` is assumed to return the gradient along with the objective function - as is our case when applying the adjoint-state method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is Dask?\n",
    "\n",
    "> [Dask](https://dask.pydata.org/en/latest/#dask) is a flexible parallel computing library for analytic computing.\n",
    ">\n",
    "> Dask is composed of two components:\n",
    ">\n",
    "> * Dynamic task scheduling optimized for computation...\n",
    "> * “Big Data” collections like parallel arrays, dataframes, and lists that extend common interfaces like NumPy, Pandas, or Python iterators to larger-than-memory or distributed environments. These parallel collections run on top of the dynamic task schedulers.\n",
    ">\n",
    "> Dask emphasizes the following virtues:\n",
    "> \n",
    "> * Familiar: Provides parallelized NumPy array and Pandas DataFrame objects\n",
    "> * Flexible: Provides a task scheduling interface for more custom workloads and integration with other projects.\n",
    "> * Native: Enables distributed computing in Pure Python with access to the PyData stack.\n",
    "> * Fast: Operates with low overhead, low latency, and minimal serialization necessary for fast numerical algorithms\n",
    "> * Scales up: Runs resiliently on clusters with 1000s of cores\n",
    "> * Scales down: Trivial to set up and run on a laptop in a single process\n",
    "> * Responsive: Designed with interactive computing in mind it provides rapid feedback and diagnostics to aid humans\n",
    "\n",
    "**We are going to use it here to parallelise the computation of the functional and gradient as this is the vast bulk of the computational expense of FWI and it is trivially parallel over data shots.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up (synthetic) data\n",
    "In a real world scenario we work with collected seismic data; for the tutorial we know what the actual solution is and we are using the workers to also generate the synthetic data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Set up inversion parameters.\n",
    "param = {'t0': 0.,\n",
    "         'tn': 1000.,              # Simulation last 1 second (1000 ms)\n",
    "         'f0': 0.010,              # Source peak frequency is 10Hz (0.010 kHz)\n",
    "         'nshots': 5,              # Number of shots to create gradient from\n",
    "         'm_bounds': (0.08, 0.25), # Set the min and max slowness\n",
    "         'shape': (101, 101),      # Number of grid points (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",
    "         'nbpml': 40}              # nbpml thickness.\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import scipy\n",
    "from scipy import signal, optimize\n",
    "\n",
    "from devito import Grid\n",
    "\n",
    "from distributed import Client, LocalCluster, wait\n",
    "\n",
    "import cloudpickle as pickle\n",
    "\n",
    "# Import acoustic solver, source and receiver modules.\n",
    "from examples.seismic import Model, demo_model\n",
    "from examples.seismic.acoustic import AcousticWaveSolver\n",
    "from examples.seismic import TimeAxis, PointSource, RickerSource, Receiver\n",
    "\n",
    "# Import convenience function for plotting results\n",
    "from examples.seismic import plot_image\n",
    "\n",
    "def get_true_model():\n",
    "    ''' Define the test phantom; in this case we are using\n",
    "    a simple circle so we can easily see what is going on.\n",
    "    '''\n",
    "    return demo_model('circle-isotropic', vp=3.0, vp_background=2.5, \n",
    "                      origin=param['origin'], shape=param['shape'],\n",
    "                      spacing=param['spacing'], nbpml=param['nbpml'])\n",
    "\n",
    "def get_initial_model():\n",
    "    '''The initial guess for the subsurface model.\n",
    "    '''\n",
    "    # Make sure both model are on the same grid\n",
    "    grid = get_true_model().grid\n",
    "    return demo_model('circle-isotropic', vp=2.5, vp_background=2.5, \n",
    "                      origin=param['origin'], shape=param['shape'],\n",
    "                      spacing=param['spacing'], nbpml=param['nbpml'],\n",
    "                      grid=grid)\n",
    "\n",
    "def wrap_model(x, astype=None):\n",
    "    '''Wrap a flat array as a subsurface model.\n",
    "    '''\n",
    "    model = get_initial_model()\n",
    "    if astype:\n",
    "        model.m.data[:] = x.astype(astype).reshape(model.m.data.shape)\n",
    "    else:\n",
    "        model.m.data[:] = x.reshape(model.m.data.shape)\n",
    "    return model\n",
    "\n",
    "def load_model(filename):\n",
    "    \"\"\" Returns the current model. This is used by the\n",
    "    worker to get the current model.\n",
    "    \"\"\"\n",
    "    pkl = pickle.load(open(filename, \"rb\"))\n",
    "    \n",
    "    return pkl['model']\n",
    "\n",
    "def dump_model(filename, model):\n",
    "    ''' Dump model to disk.\n",
    "    '''\n",
    "    pickle.dump({'model':model}, open(filename, \"wb\"))\n",
    "    \n",
    "def load_shot_data(shot_id, dt):\n",
    "    ''' Load shot data from disk, resampling to the model time step.\n",
    "    '''\n",
    "    pkl = pickle.load(open(\"shot_%d.p\"%shot_id, \"rb\"))\n",
    "    \n",
    "    return pkl['src'].resample(dt), pkl['rec'].resample(dt)\n",
    "\n",
    "def dump_shot_data(shot_id, src, rec):\n",
    "    ''' Dump shot data to disk.\n",
    "    '''\n",
    "    pickle.dump({'src':src, 'rec':rec}, open('shot_%d.p'%shot_id, \"wb\"))\n",
    "    \n",
    "def generate_shotdata_i(param):\n",
    "    \"\"\" Inversion crime alert! Here the worker is creating the\n",
    "        'observed' data using the real model. For a real case\n",
    "        the worker would be reading seismic data from disk.\n",
    "    \"\"\"\n",
    "    true_model = get_true_model()\n",
    "    shot_id = param['shot_id']\n",
    "    \n",
    "    # Time step from model grid spacing\n",
    "    dt = true_model.critical_dt\n",
    "\n",
    "    # Set up source data and geometry.\n",
    "\n",
    "    time_range = TimeAxis(start=param['t0'], stop=param['tn'], step=dt)\n",
    "    src = RickerSource(name='src', grid=true_model.grid, f0=param['f0'],\n",
    "                       time_range=time_range)\n",
    "\n",
    "    src.coordinates.data[0, :] = [30, param['shot_id']*1000./(param['nshots']-1)]\n",
    "    \n",
    "    # Number of receiver locations per shot.\n",
    "    nreceivers = 101\n",
    "\n",
    "    # Set up receiver data and geometry.\n",
    "    rec = Receiver(name='rec', grid=true_model.grid, time_range=time_range,\n",
    "                   npoint=nreceivers)\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.\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",
    "    dump_shot_data(shot_id, src, true_d)\n",
    "\n",
    "def generate_shotdata(param):\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",
    "        generate_shotdata_i(work[i])\n",
    "        \n",
    "    # Map worklist to cluster\n",
    "    futures = client.map(generate_shotdata_i, work)\n",
    "\n",
    "    # Wait for all futures\n",
    "    wait(futures)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "CustomCompiler: cache hit `/var/folders/qr/4hc1nw7d3hj55jbd1z7kkrvh0000gn/T/devito-jitcache-uid501/ff78393c8a1120844c3a982e7e2a89d3e5dc3820.c` [0.04 s]\n",
      "=========================================================================================\n",
      "section0<714,181,181> with OI=1.95 computed in 0.201 s [3.72 GFlops/s, 0.12 GPts/s]\n",
      "section1<714,1> with OI=4.60 computed in 0.000 s [0.54 GFlops/s, 0.01 GPts/s]\n",
      "section2<714,101> with OI=3.69 computed in 0.001 s [3.03 GFlops/s]\n",
      "=========================================================================================\n",
      "CustomCompiler: cache hit `/var/folders/qr/4hc1nw7d3hj55jbd1z7kkrvh0000gn/T/devito-jitcache-uid501/ff78393c8a1120844c3a982e7e2a89d3e5dc3820.c` [0.04 s]\n",
      "=========================================================================================\n",
      "section0<714,181,181> with OI=1.95 computed in 0.196 s [3.82 GFlops/s, 0.12 GPts/s]\n",
      "section1<714,1> with OI=4.60 computed in 0.000 s [0.68 GFlops/s, 0.01 GPts/s]\n",
      "section2<714,101> with OI=3.69 computed in 0.001 s [3.35 GFlops/s]\n",
      "=========================================================================================\n",
      "CustomCompiler: cache hit `/var/folders/qr/4hc1nw7d3hj55jbd1z7kkrvh0000gn/T/devito-jitcache-uid501/ff78393c8a1120844c3a982e7e2a89d3e5dc3820.c` [0.04 s]\n",
      "=========================================================================================\n",
      "section0<714,181,181> with OI=1.95 computed in 0.191 s [3.92 GFlops/s, 0.12 GPts/s]\n",
      "section1<714,1> with OI=4.60 computed in 0.000 s [0.94 GFlops/s, 0.01 GPts/s]\n",
      "section2<714,101> with OI=3.69 computed in 0.001 s [3.93 GFlops/s]\n",
      "=========================================================================================\n",
      "CustomCompiler: cache hit `/var/folders/qr/4hc1nw7d3hj55jbd1z7kkrvh0000gn/T/devito-jitcache-uid501/ff78393c8a1120844c3a982e7e2a89d3e5dc3820.c` [0.04 s]\n",
      "=========================================================================================\n",
      "section0<714,181,181> with OI=1.95 computed in 0.191 s [3.93 GFlops/s, 0.12 GPts/s]\n",
      "section1<714,1> with OI=4.60 computed in 0.000 s [0.85 GFlops/s, 0.01 GPts/s]\n",
      "section2<714,101> with OI=3.69 computed in 0.001 s [3.85 GFlops/s]\n",
      "=========================================================================================\n",
      "CustomCompiler: cache hit `/var/folders/qr/4hc1nw7d3hj55jbd1z7kkrvh0000gn/T/devito-jitcache-uid501/ff78393c8a1120844c3a982e7e2a89d3e5dc3820.c` [0.04 s]\n",
      "=========================================================================================\n",
      "section0<714,181,181> with OI=1.95 computed in 0.188 s [3.98 GFlops/s, 0.12 GPts/s]\n",
      "section1<714,1> with OI=4.60 computed in 0.000 s [0.94 GFlops/s, 0.01 GPts/s]\n",
      "section2<714,101> with OI=3.69 computed in 0.001 s [4.15 GFlops/s]\n",
      "=========================================================================================\n"
     ]
    }
   ],
   "source": [
    "#NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Start Dask cluster\n",
    "cluster = LocalCluster(n_workers=5, death_timeout=600)\n",
    "client = Client(cluster)\n",
    "\n",
    "# Generate shot data.\n",
    "generate_shotdata(param)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dask specifics\n",
    "\n",
    "Previously we defined a function to calculate the individual contribution to the functional and gradient for each shot, which was then used in a loop over all shots. However, when using distributed frameworks such as Dask we instead think in terms of creating a worklist which gets *mapped* onto the worker pool. The sum reduction is also performed in parallel. For now however we assume that the scipy.optimize.minimize itself is running on the *master* process; this is a reasonable simplification because the computational cost of calculating (f, g) far exceeds the other compute costs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we want to be able to use standard reduction operators such as sum on (f, g) we first define it as a type so that we can define the `__add__` (and `__rand__` method)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create operators for gradient based inversion\n",
    "To perform the inversion we are going to use [scipy.optimize.minimize(method=’L-BFGS-B’)](https://docs.scipy.org/doc/scipy/reference/optimize.minimize-lbfgsb.html#optimize-minimize-lbfgsb).\n",
    "\n",
    "First we define the functional, ```f```, and gradient, ```g```, operator (i.e. the function ```fun```) for a single shot of data. This is the work that is going to be performed by the worker on a unit of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "    # Load the current model and the shot data for this worker.\n",
    "    # Note, unlike the serial example the model is not passed in\n",
    "    # as an argument. Broadcasting large datasets is considered\n",
    "    # a programming anti-pattern and at the time of writing it\n",
    "    # it only worked relaiably with Dask master. Therefore, the\n",
    "    # the model is communicated via a file.\n",
    "    model0 = load_model(param['model'])\n",
    "    \n",
    "    dt = model0.critical_dt\n",
    "\n",
    "    src, rec = load_shot_data(param['shot_id'], dt)\n",
    "\n",
    "    # Set up solver.\n",
    "    solver = AcousticWaveSolver(model0, src, rec, space_order=4)\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 = Receiver(name='rec', grid=model0.grid,\n",
    "                        time_range=rec.time_range,\n",
    "                        coordinates=rec.coordinates.data)\n",
    "\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",
    "    grad = Function(name=\"grad\", grid=model0.grid)\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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the global functional-gradient operator. This does the following:\n",
    "* Maps the worklist (shots) to the workers so that the invidual contributions to (f, g) are computed.\n",
    "* Sum individual contributions to (f, g) and returns the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "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.p\"\n",
    "    dump_model(param['model'], wrap_model(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, retries=1)\n",
    "    \n",
    "    # Perform 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)"
   ]
  },
  {
   "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": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import optimize\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 box constraints and to monitor the true relative\n",
    "# solution error.\n",
    "relative_error = []\n",
    "def fwi_callbacks(x):\n",
    "    # Apply boundary constraint\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=0.1, maxiter=5):\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_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": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      fun: 215.09233345324492\n",
      " hess_inv: <32761x32761 LbfgsInvHessProduct with dtype=float64>\n",
      "      jac: array([5.97415130e-12, 2.72497060e-11, 6.92164034e-11, ...,\n",
      "       1.14447021e-10, 4.60593820e-11, 1.02192577e-11])\n",
      "  message: b'STOP: TOTAL NO. of ITERATIONS REACHED LIMIT'\n",
      "     nfev: 6\n",
      "      nit: 5\n",
      "   status: 1\n",
      "  success: False\n",
      "        x: array([0.16, 0.16, 0.16, ..., 0.16, 0.16, 0.16])\n"
     ]
    }
   ],
   "source": [
    "#NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "model0 = get_initial_model()\n",
    "\n",
    "# Baby steps\n",
    "result = fwi(model0, param)\n",
    "\n",
    "# Print out results of optimizer.\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUAAAAEKCAYAAABjU4ygAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztvXm4ZUV97/359Tndp+k+bTfQIsogKqiXGBMUNT6aqHG46E00en1vcJ6uaKK5QY25DnnRq28SnMMbSQhOqDEap6tEUdSI4w0oIsqkCeJAKwgNNnQDp7vPOb/7R1XtXbt2Va1aezhn7z71fZ7z1Kph1aq19jq1vvWbSlSVioqKirWIdas9gIqKiorVQp0AKyoq1izqBFhRUbFmUSfAioqKNYs6AVZUVKxZ1AmwoqJizWKsE6CInCQiPxSRq0XkVZH6ORH5Z1t/kYgcM87xVFRUVPgY2wQoIjPAmcDjgeOBp4nI8UGzFwC/UtVjgXcAbxrXeCoqKipCjJMBPhi4WlWvUdV9wEeAJwVtngS83x5/HHi0iMgYx1RRUVHRwewY+z4CuNbL7wAekmqjqosicgtwKLDTbyQipwCnALB+8wM55L5wkKlbf+heAI7hJ53213M4AHsW5wHQO2ZMxT7bYNGmS95FXNlyQ+o7ziwFZZpoG3O2KXXAkcSxn18X5GPntE1j1xsXVtsZKXX9kt+tzdjbtE09+9zvFdKZmaDNTJD6xws23W/TX31np6reGeBYEb29ZMzAdXC+qp5U2HzVMc4JMPYThq9ASRtU9WzgbAA5/ETlWRfD/UzdXZ7zHwCcyfM77U/nfwLw9Zt+B4D9l9/JVPzENnDT6y7vIu54j00XEvkFugjrFhvyPhYjZTHMRo7DdGMiX9I2labK2qLkPkufxbiQun6sPCxrypfWhQifecnvFf6280G5y2/zznHHP7Dp9Tb9iPzUNbkdeFHRoOH1sL2w6URgnBPgDuAoL38k8ItEmx0iMgtsBW7O9roO80PaH27XbebgLze/ttPkF9y19xz3w4c/zcbIcTjxpSZEv9+myTLsG7r/DGEaYjFyHL78Ybl/TjjxLQb51ASZQ+6tGecbFcMwk+cgE98wE+AgSD3P3AdrY5CGE58/AW4PynrWXgaSGca0Y5wywG8Dx4nIPURkA3AycG7Q5lzgOfb4qcCXtUZnqKiYKKzDSJxK/qYNY5vYrUzvpcD5GEnDe1X1ChF5A3Cxqp4LvAf4oIhcjWF+Jzd27D5H9ku75yd3BuCiIx/cabK0aG5r/64tveeGX8TY1zNkba7tniDvl7k0xSJd3r9euExuSmNoYoSxslQ+ZIh+WWopXPL2jOINyz2D2D237aOExY1y6dvE6mJtSsQd4TvblPePM7+TAOvT1VONsTJbVT0POC8oO807XgD+n3GOoaKiYjgcyEvg6bsv92s4BmWFtguLh3TbhGwuVELk5Fwh29mTyPtlG4O6kPHF5IdhXU5hQtAmxT5ycqEUA0zJBmNlbeSFbd6sVNtwzLG6FIZhaIMwwJLrDCLPa1J8xRhgKRNMXTtAZYAVFRVrFpUBVlRUrFlUBjhJUMwyI1wy+jZ94V2lliVtbOxy5zfZbIXLaOhfnqeWzTEzmFKTDL8sNdbUUjhWFi6FSeRL65pQouAoVU4Ms6zN1Q1j9tLGBjOl/MgtgcM0d51Unq4W+EDE9E2AFRUVK4rKACcJoUCixFwkbFuCYRhhE6uLtQlZYoz9NJnO5BC2iZnmxPIlbWLeCCmM641rYsVtjJvD+mExiFdHKfNrwwBj56TG2K56anGg3ldFRcWIUBngJKPkDkpYQZOcMObOFn6dU/K2EgbYZEQdazOI8bRDKSOMlbVhjcMYUedkdKXy0JLfehTubG2eW8nqYpQMsIRpZlC1wBUVFWsWVQkyiWjDIIZhSuF1fDeiVH+h9jRmkJ3SArs0dMGLtQn7z7GdQVjNKOSGKZlgzF2vaRxtmG0JAxxG9jeIoXcT84vJUpuYX44Bps4tYewe6hK4oqJizaIugScRg9hs5TBKGVUqbcMAQybY5pwSzXEJmtqG8stYWcgaY0wj5/Lm15docsfFAFP2k7n3pi3zy8kA29gBpvofRFZLZYAVFRVrGJUBThKWMawiFegAmu8q5vUwCqQYYY6ZhXW5CNSlDDB2To6NxsaeQ0o2COXMr0Rrmrtuk/Y3d3/DsOFB/mOankHsmTQFoGhjO9jm2VcGWFFRUdGFULXAFRUVaxQCrC+dKUa9qhozpm8CXMYsE8MozDnhc1ieyg+KJmVIzog6XKKG5SVL4JJlc7g0TF3fR+nLHFvat1nukSlrGk+piVMbQ+iS8TQZOcfKSo2bc3UlZjCDLJsz9yoCs3UCrKioWIsQgfUzze2mEdM3AS5hQl/lvoRNJgRt7rqELZQK1WPRllP5NozCpSWs0bUNhfqjMqJO3dcgzC8cW5vrtfndhsEozWBiSpBhtjEdxBUupgRpwwCnDOPcFa6iouIAgAisnyv7a+5LjhKRC0TkKhG5QkT+NNJmq4j8i4h8z7Z53jjuCw4EBhjb7Dn8Anb2QrA7bs5GPv2LwaNYtHu2lxjhpsrbOOQ7pBiaf9wm5FQTWwzHE5PnhRgkeECK4TaVte23qbypLjWOJpZaYlrSJKMr2dB+EFe4YeWvozUEXAReoaqXiMgW4Dsi8kVVvdJr8xLgSlX9fRG5M/BDEfmQqu4b2Sgspm8CrKioWFmMcAJU1euA6+zxbhG5CjgC8CdABbaIiGDoy82MSb0yfRPgIoYBbvPyIdxdOeY3vx+AdRt7PyDLixHJrmOCKSPmNkbGK+3CFZOZNcmfwr5jmuqwrkQ2l9vdLsSkvoUlrK6NTK6JzeVkgIMwwBIZ4Cg04r3YLiIXe/mzVfXsWEMROQY4AbgoqHoncC7wC2AL8Iequlw8ghaY1FevoqJiUiBAuRZ4p6qe2NilyDzwCeBUVb01qP7PwKXA7wL3Ar4oIl+PtBsa0zcBKr2ucCWO6xZ9jM/PL1hnn9Cmro1rWqlNWqysRJZVGnwh12+TFjj27FJtY2ijhS1Fm7e0pO0o5HklzLr0nNi5w2xw1CR7jCH3Do3YGVhE1mMmvw+p6icjTZ4HnK6qClwtIj8G7gt8a3SjMJi+CbCiomJlIUCBhreoKyPXew9wlaq+PdHsZ8Cjga+LyF2A+wDXjGYEvThwJsCY9tKxt0Xnyr2+t22Mze1JpGF9rJ823gijwCA2iuMQJefklMPY4+XezpT9WorJlvTVxqavRL7WVo43CAMcRAscw8oxwIcBzwIuE5FLbdlrgKMBVPUs4I3AOSJymb36/1TVnSMbgYcDZwKsqKgYD0arBf6G7THX5hfA40ZzxTzqBFhRUdGM6go3QYip8f2l6c6gLFwexZbATUve2C5tbcxeUmgSyA+LUlexWPu2ht65ukEMykO0UU6E9S2UZdnlbKlpS5u2o14CD6P8iJkvHcARUQ/Q26qoqBgZ6gQ4QXA/hvtiObaX26/XIcdGSkNMxRjgKJhfidC+6dfKmbA0lZewuRJlT1vFUBvkWFyTofcgpi0lrmk5tlVqwjIqBtiG+TmUMsARaYEnDdM3AVZUVKwsKgOcIChddzjwTF0y5wzyJSxxhRvGjWgQBthU1+Y+U+Ul8rySZ5JjFEloojyjNCxluG0MpHNGxk3Mrw2LG+TcNjLHsLyEsVcZYEVFRYWHdq5wU4WxToAichJwBubxvVtVTw/qXw78d8x36Ebg+ar602ynLiR+qOmNsZDwazkf5NtoxwbRYpb0X+JeVnq9QZBjsaNggH3QSIP9DYOM7UkWlKXuI/Z8UwypjTxvlAwwdr1SDXLuHIc2Mts1xgDHFhBVRGaAM4HHA8cDTxOR44Nm3wVOVNX7Ax8H3jyu8VRUVAwIpwQp+ZsyjHNefzBwtapeAyAiHwGehBf3S1Uv8NpfCDyzsddFDPu73uY7XyyfRbjbsrIjP1iqj9jXOuzCv27J2EaNUTPDJrllGy1wLiBFB/uDdDHIlyDl3+aXJQLYlpxaogVuksGVnNPGhnCQ0FZN/80xxp6zcnA4gBngOG/rCOBaL78DeEim/QuAz8UqROQU4BQAZo4e0fAqKiqKUCfAgRBT3UXVfCLyTOBE4BGxehtQ8WwAWXeicj2w4BjEzZEzDjOJ+2puS6QlX+sStLW1i5W1kTWWap9zbUvKmwK+RuF+4juCk3IMMMUGnZwvRuNc2UFB26B6ENvMEi3wMJ4Zw4S2GuS9bGP/Oog8e4oxztvaARzl5Y/ERHjtgYg8Bngt8AhV3TvG8VRUVAyCqgUeCN8GjhORewA/B04Gnu43EJETgH8ATlLVG4p6dQFRccFhHQO4U7fNkUGaYn7zdJHS6o3yi1viMTEKG7tBmGbueknGF7I9/yRXNowMMGwb0waHP1QgCyxBk1bYPx7ELm8YeV7J2EIM8xtXGeBooKqLIvJS4HzM9+O9qnqFiLwBuFhVzwXegpmGPmbiJPIzVX3iuMZUUVExAKor3GBQ1fOA84Ky07zjx4zz+hUVFSNAZYCTiENNcrjNPtyrOtamKYPnEqPVEiF0kyIhZ2KQqmtzzmyQzy1Z2441dk7fkjdUdADcHpSVLH1Ttio5G5YQQdTvsIuSJXHuty4NZFASnKBE3NLGNdKhackbU3SUhHKrE2BFRcWaRZ0AJwgzmJ1CnYLjZJMc8qqfd5rcvscoRhauPsQUOKNpP2gq9N596CY3H7QJDaVjCL+muf1EmvYgibn4hQFewzSs98cS5kvYQQeOtYXM7/Yg73cQnpMaUAyh/1ruNXWMz10nUIYMg5xSosSdrcltLud6FxuDj5xlUBslSKkZTNUCV1RUrElUBjhBmMPI+OzWy/On3gjAf8wc22ny/K3vBeD8Y/8zAAuLlgmGd+uzOsf4thk517p5w242bDSmiXMb9wEwM7vUN6Qlu7/w3oUNAOxbMCqzZZtnT7DnsH/swnqlGIV/Tsj0QsTkOE0MIisjCxmeS50JUmjy4h8PYvbi2NwgDDA8J9gBMIcSeW9TmzaG0G2u14Sc62KJLLrkPaha4IqKijWLygAnCOsxml+r/d1z/XYAnnKv/91p8hOOAWBh58GmINSKxr7WlgE65rdl224ADpoz+TksA6SfATr5yL45w/j2bjWfyztuM7LI2zduAmDZpmYM0juWQQyvQ4Ta4Rgajalv9453B2UuHzJCnwGmtMAlCNncpiAfa3tHkE8YTbcxjHbIGUKvlBa46X0YNnhFSUCLOgFWVFSsWVRXuAnCDL3hrX5gmNRXF07qlrmvWhg01WGeJJyMz6WO+W1gb0/ex5J9O2bsZ7TDEjf3tvO51fJiUNkGqS/9xiBfcm4Hjjnt9spuDdI7gjYlWuAUE8y5tR0UlDsm6A861PqG14tnk2WxYeTqSuR4bW1L29gBtnF3zNmJVjvAioqKigyEMjOwKcT0TYDrMD+G+1LtsKnP8lJfwpwW2LbZbzW3Hc3ujMk7VrcUkwEm4BjhBqtBXlzsriMW5u3x4sbeMTp2Gub9spDphfnYr5qU/TnvDsfyfJ7qynYH+ZwMMNQCp+CfE8r+ShBqfcOgC+6+BgiOEF7CPy6VBebalMh9h9ECt/EuKvGUqUvgioqKNYu6BK6oqFjTOEBnium7rfBr5Oj8rkjbkNZvDMr3ROpmzcHu2S0ALNml6tJmc9F9bOic4pbF3eXxTE8aYtYzol5nj5c32qXaYrBUC8fql4Vpm713+5Y64RootpxNGULHgiGULoFjSpDUmHIS/4TyIzSILrlMDk2GzyX7epQoP5pQMtZSg2j/uC6BKyoqKiKoS+AJgmK+VoOEC3KImQMEDNKZqeyxSpG91r1tbmM3ar9TbszMOCZoOl4KBrO02P/5dGY2y7N2MBsD160SBjgSI+pUxGb/OGyT6mNYDKIEGQAh2xnGFa4NA0yGZ4tulWPHmgjqEGOvpS5wgzDA6gpXUVGxJlEZ4ATBMUCH3Je3ScYTcwnaFeQXDDPbbwMa7J/vsp11lgG6gAmzgRG1Q4wBduDaztqfwrnIxeR7pYa0gzjVxwdn0zDYaIlsLjRmTvXt93tQUOfym4J87Jxwd7jMjZcy5mHNYMIQa47pOdYfCazRB/fuLNoLhYxwVGYwDrFnUifAioqKNYs6AU4QljHa29C42TcYDl3dUi5BOafw8GvZYWRdreKyldstzFvNsP2ir++EzhrA+raNHGqQ3cX6wmyFu+rF5Hkh20rt+Ablfma+djZkbWE+ZHn+ccgOt9g0YNIle0CPQgbov3udOvN81gUh1cL3Y2lx1js2zG/ZFSwW/Ks2vcMlwRBSl6la4IqKijWJygAnCIsYOd0eLw+9dxKGsw+/iLFw82Gb8NwYa+zIGi1jsYxwv/1677df+nUlsp4USlhIiTtW6GLXsYV09oeOAfosK7T/S4W5bxPs1CEmA0y1CeV7fl0gEwtXBDFmNp9Ic9raJq1vWA8wb16mcEWQlhV3885tshNct8MWA2uBnAywJBhC+P+TkgEeoFrgdas9gIqKigmHY4Alf01diRwlIheIyFUicoWI/GmkzStF5FL7d7mILInIIaO7oS6mjwHux2xytN3mDw9SH2EEi/ArFwtR75Bifjm7qw7Lsl/pWZMuh9q/nusE2r0S74Q2cqjwCx+2CdmwJ+Ps3Idjh8MGFIjlS9rkzhmGFYfPYluQj7HGFAMM5H3QZf7OdjQMtebsR5eW0gK2JWcvGla0sQPMBUQNMX4t8CLwClW9RES2AN8RkS+q6pWugaq+BXgLgIj8PvAyVb0516mIbAYWVLXVcqsywIqKijxGyABV9TpVvcQe7wauAo7InPI04MN9QxJZJyJPF5HPisgNwA+A6yyrfIuIHFdya3UCrKioaITOlP0B20XkYu/vlFSfInIMcAJwUaJ+E3AS8IlI9QXAvYBXA4er6lGqehjw28CFwOki8sym+5rOJfAOuktgt+fvdq9NuJSJCaghL0AOy1P1sTZ9AnK3/0eBY37q+pBe3oVmP/6vGhpUp0wiBgkMkMr7ZSVtS9FmCZ67bpNCo80SuM8kqbsC2xAsfTfMuX1lghuJrICTxvO593GY3zajBNF1sK88IOpOVT2xqZGIzGMmtlNV9dZEs98HvplY/j5GVfu0b7btJ4BPiEgu4gYwjRNgRUXFikIFFmdKF4t90so+2InpE8CHVPWTmaYnE1n+ArjJT0TuBexQ1b0i8kjg/sAHVHVXbIIMMX0T4BLGDGaHzW8PUuhnftuCfIwRhl/HJibotwmfYti21PUqdZ0QKeVOTgnS1O8gioacuUgbBpiKZBXmcwL/EqVViDaKkyYliFV0xUyeYntJl2KpzxWOeBori7m+hec4pFZJgIqwNFs6VfTvmeNDRAR4D3CVqr49024r8AigaRn7CeBEETnW9nsu8E/AE0pGO30TYEVFxYpjaWZkriAPA54FXCYil9qy1wBHA6jqWbbsycAXVPW2hv6WVXVRRJ4M/I2q/q2IfLd0MFM4AS4Dt8P11gXKMUF/p7jwa02Qj7GTnHwwlvfLUowvdKPz0cQKYwbeYZCHWMisFMLrtNnTInDlcrKtucC1C2DDjA0MEQSLdQiDx/pwZYudwLJmcPv2upBk3WC0HQNhV7YQGAjn3L9K5Z2xZ5+UH/bK+6Br8Nwnzwuyzgxmn3d/e/vuL2EmlXNrK3mHO4O1aYwBIskgv22hqt+gz3o92u4c4JyCLveLyNOA52BkhpCPttuDKZwAKyoqVhKKdD5KE4jnAS8G/lJVfywi9wD+sfTkKZwAl4E7YI9lgE4LHGOAoTuba1Oi0QrPDfPhcaxtqo9YmxJZWZNhd+wrnmBvm+aNO9umGbML3EF2N7hNnpvbJuJ1c3aP5IOCPHQZnttHebZgF72Q8Tm2sdduP7BvzrCh2+e6bnr7ttoyGwzhDusu5/Ld1JTvW+r6cu0Lgty6nQD7Qk/lEIS0KnF3dEwwZISdHQgXumPcHzLbYXZ4Kwl5lVlNKMK+CfOFE5Gzgc8BX1LV/+HKVfXHwOml/WR/aRE5EqOJ+W3gbhgn0MuBzwKfU9VmlU9FRcVUY5RL4BHivRgbwZeLyD7gC8DnVfV7bTpJToAi8j6MhfZngDcBN2C+E/e2F36tiLxKVb+W6eMk4AyM1OPdqhqdmUXkqcDHgAep6sVtbqDjyrXTK3NywZCBua9lTEYYY3gx5GSBgyDF/GIa1o68ybrWWTa3cd4wNCeTO2iuu7fvFrunr2NvLt+UAswHZa6PkBnOeZo/x/xcWSgDjCHcTGqvZRtuA6qQzQHcETC93TYMVpg6Zrh7Zkvn3Ns323M2b4lep8M8PdYYsrfFXJBbi8WGti78VSf0lScDTDK/PtdF/4KJNIemlYgb64RNgKp6Icbg+fUicijwOOAVInJ/4BLMZPjRpn5yDPBtqnp5pPxy4JMisgGruYlBRGaAM4HHYqakb4vIub7Pn223BfgfJKzBKyoqVhcTLgNEVW/C2At+GEBEHoghaY1IToCJyc+v3wdcnWnyYOBqVb3GDuojwJOAK4N2bwTeDPxZyYCN954Xsim2LaaTC4Z359o6m0HfXq5JGzsu5tfoldC15Vzv5HaW8TmmF7K3bd7DCMsOtum2hjTWr0tDmeCmvV3GuWHBSEWcSExSz8973upM3Oz/2L6Nxuh275xhRCHbgxjjm+/J77IC327+YO/c3rZ7Euzx9pnu9W6fsbJEK48MWWJMvrfYKQuYXkfWGMgcFzzlaBPjG0S7nbP1bJABhht9TQpEZBvwbOAYvDv05YI5NN6ViPweZpK6u20vpn+9U/ZEs3y+1svvAB4S9H0CcJSqfkZEkhOg9Se0PoVHNQ25oqJihDBKkA3NDVcH52GWwpdR4oYSoGRa/xvgKcBlqprZv68PMVufzvkisg54B/Dcpo5U9WzgbHPeiYrHBDpfLJ8BpuQZ4RfS1xyntiwM0SaUUFY7G6Qd/1PrUWBZ3pZtXZmcY3yOxYXs7lBu6smbsp0957g2rnx7cI5/7rbbbjFDvMUWOI9Nl3cmqr4cam9QVuAEIZYIrbfPYv2ceY83b7SdbHap5xK61ST77Wd491ZzsmN+jvF1890f2x3v5NBo25BF+mWOPXZY6Uxv6k8UiwEr7Mj4QsaX264yxQRjQX1T4a5y73TuHbUw+5BN7BJ4o6q+fNCTSybAa4HLW05+YBifT9eOBH7h5bcA9wO+YrxjOBw4V0Se2FoRUlFRMUZM7hIY+KCIvBCjrO3YYzXFD3Qouas/B84Tka8GF0j68Vl8GzjOGib+HGNO83Tv/FvwPHhF5CvAn9XJr6JisjChZjAO+zDBU19Ld4WpwD1LTi6ZAP8SQ7g3QrkgwPrnvRQ4H2MG815VvUJE3gBcrKrnlvbVA7c/Qd/eFgP05S8VQuPpkgAGqbo2e8Z2UqPs2GiXvFu2mhvyzVK28Sub9io0ustak96FX3bOcWXb7ZL3MG6w5XYJvNfUb77Bik9u8sbovqG3JFK3BO7aQXeXYa6sjSnGXJDfnEihswReb9NDtpoLH3KI0YDpoSbdeYh5wP5y9ia79HVL4Jvsd9iV/8q2vcmLsBEupXf1yE88U57Z/omio/RoMmr23+HSvX1ze/ym0MYg32KCJ8CXA8eq6s7GlhGUTICHqOrjBulcVc/DCCn9stMSbR85yDUqKirGiwlngFcAtze2SqBkAvySiDxOVb8w6EVGihkMawqFwLGvX7gHrkt3Rdomg5kG7XJfz/Ccop3JzA3MW+Y3v9mkobmKf5xidWEeumywww6XTP5Ov7TmNa5pmEKX6TlW6JQg7jmOmgGGTNA9twwD7Gxn7PKGxOG20LnzoWawdz6sS68OPcw8p23cxaa9SqUtthPfwNu5+4WG3WEAh0UvakrHcLtj7mIrQsaXU2iUBrINj320eYcjM4IiHeP0CcQScKmIXECviG40ZjDAS4A/F5G9mHjMpWYwFRUVBwAmnAF+yv4NhMYJUFW3hGU2qOHqYZY8s2gKqJkKFBk7JxdCKxX8oCjEei/z27a5l4XkTFqc6cphlt3dxdI2x/aiDPA2I9Db6PTwKeZXIgO8LUhjZjDhMw/NYWL/T+HzmwvSGAN0ZYcGY3abKIZjBg65zQx402E/BeCgzb0ufRsibC/F/JxBtGNId3gmWrNhgIQ2DDDc97pNmPsmJjgAA5zgCfByVf2OX2B3kitCY5xrq7Tw8+toEW6moqJi+rHITNHfKuBdIvLrLmNjA/5F6cklS+CjReTVqvrXIjKHCVpwSftxjhg5BrgtSEMNr4P/xQ2DmKYCiC5Gyord2roXDGV+oVFzV97XVW6FzG97YNQcGjlDhvm5vGNMvwzyUK79zTHApSAfg3tu7v/H9R/KAn1ZozveHORThtg+GbNj2WjTu9ztxp6+usEZui9Bx/UtEaDhjiAoA0RC4adWICVa4FFuXpVjgBFMsisc8FTg4yLyDODhGLe4YqVtyV09D/iQiLwaeBQmDNY7BhlpRUXF9GGSl8Cqeo2InIyRA14LPE5V72g4rYNcOKwHeNkzgH8Avgl8VUQe4DY3XjWEbOtwr+5Im6Y2Q4rddSgPHCTcfDJAqbHPXL/RCxtlj11AARdgIBeeqjSgwcFLXbnhRsfoXOpkfE1yvVhZyLJiMkBLetSWLZbsB2T7cyZ0EoYmi7HIQcLaJ+o22t/64I3mue2bcQEPuprPkOm538Uxvg1WftgTHNYxwA4TDML252TSKRlgDql7bWPbmgiGsHfCfIFF5DI811qM1HcGuEhEUNX7l/STDYcV5H8FHG/LFfjd8uFWVFRMKyZ0Cfx7o+gkFw7rUaO4wMixDsPOHOO7r01/02vj6sJACaHWLYYUo2gT+qrvi2tOnpntduI2D3KMIWSCm7Ih6huCm97ibYcasriU/C4mzxtE/pSo29/i+bnHljU1CJ/xbUHekbeZIA9dueEtvfktm81zu/2Qfva9xQqS9/SFBOtNN3i2gzMzDfQ3xwCbNjhyyIW4aoOG32cCl8A3qWrW/0tE5pvaJLXAIvLMnLmLiNxLRB7ePM6KiopphpMBlvytID4tIm8Tkd8RkY6BlIjcU0ReICLnUxAUNfetOBSqDSr9AAAgAElEQVRjYf0d4DvAjRjudSxmw+KdwKuGuYOKiorJxyQqQVT10SLyBOBFwMNE5GAMj/0hZs+i56jq9bk+IL8EPkNE3omR9T0MuD9mU6SrgGep6s+Gv40BMINRbjhFx8km2frU/nu95XK7FnbbL//EpjHBcm73t1IkzQ7MUsg3jHW7pW0I9s4Iyzd5bo7dPTj29rQJDXglt5RaCspLXNZipiSrjdSy3KWhwsY3oUmIBcSmm7baZe1M99k7UUT3d+nNu99vxnuQ7tjtGNcXrTMnWmij/HBoMnspWRonlSCT5woXizXQFtlHoqpLwBftX0VFxRrEJDLAUWHiVDuNmMVEEXykyT78GWZufguv7DT5LE8A4IMPfDYAP91lNSXOPjgWDMFhFMwvAd8wNnStcmyhyygW+9qF7LDLBK1CZa8VwEeMfovNRmIsbyZIc0gw6PW2PKYMWd/0/GJmMHNBXdgmxQhjdYESYsba7jhFFfTvd9z9fRbtcPr3QU7uiVyiaBtkb99wFTNC1AmwoqJiTWLSd4UbBtM3AbpwWMea7Iv4BwA+L939kJ+qRnbzrzwGgJ/OWwaY2iPER8hghtn5rQCOQYRsYbYjU+pnjSE77Nt71x9ziStaSX1LOGPmzuO042hkezGM+n8vFajBpjOLVlrnib3C3ylMu+1GLCgd8/vXd53p2xXurcD7VPWKQc5vvCvr//tf6d927g2pcyoqKg4sTPAS+AfA2SIyC7wP+LDdbqMIJdP6pzEmo9+hV5KyOljGyGp+YrLvs5vKna4/7DT5lJUB/tIGvOxzKxo1VuorPUpMkkZ3tTHm/+2ZUAs8mWQqiUneFlNV3w28W0Tug4lb8H0R+SbwLlW9oOn8kp/iSFUt2mW9oqLiwMOkywBFZAbjE3ZfjKrze8DLReRFqnpy7tySCfD/iMivq+plww91BFjE3OKFJvvlY4xL4BOfdHyniZNX3HjF0aZgh60ImeCov8SrxAQHWp60OaWNHaB9BmrTMBhCziVuIPmgQ+p+wlBbfll4PdtmadY4SPlyL3fcCX2fuOBQS8VcmKo279YgASIaA6JOJm0VkbcDvw98GfgrVf2WrXqTiPwwfaZBLhqMi7YwCzxPRK7BLIFdSPyiaAsVFRXTjwmWAV4O/IWqxjZGenDTyblpfSTRFkaOJYwdn/PusLZb11/qbQN6eG9dXzCEkuCcKY1xSUimRPnSYvclWprrDbrpGEW4yY4ve+kG5Zzrybs+Olsyzno+B+6S4/iAt2AlbYIhJBG7h/D/sik0GfSH2g8C1+6d6w1+6h+7Z++eufOQyC0RO797G0+NNt4bY8aEG0I/Q1Xf6xeIyL+q6qNLlCE5V7if2s4+qKrPCi7wQeBZ0RMrKioOKEyiDFBENgKbgO3WD9gFbrkTcLfSfkq+L78WXHgGeGDpBSoqKqYbRgs8cb7ALwJOxUx2fnDmW4EzSzvJyQBfDbwGOEhEbqU7w+4Dzm472pHBmcG4Za1Ld3ht3BI4jAidwyBLjiaH/E5ql7XeEjhc4obL2m6+++Kl2rg+9tpIxsx58QDDfUlcd6nl1yAf+pUyqckpMlL3GaaxtjaYktrURX32n/3ezl4gB/Wk/b9f95wka0qNPbZHzSiWwCUin1SeyVwCq+oZwBki8ieq+reD9pNbAv818Nci8teq+upBL1BRUTH9mLQJUER+V1W/DPxcRJ4S1qvqJ0v6Kfm+vMZe4OEYrfDXVXXgjYhHhnAnLf/r6VjhdpuGu8OFDMA/LhGip5BigAtmL4h9Cx6bm7eMwbI2xyg2BbuL+eGwwjIXkqkbKdqWb+2es/k2qxBxISPD1IWEiu28Fu7OFu6YFzPRSLCNEhOXvj1BEkqKnrrUfbnfemuQgpESQXfvYFu3e6v5ndxz9nd4c8eOHd7R2RXOsUWnkPIUXamXJrl3jNemULEWRUnghKYxeZhEGSAmJumXMSYwIRQY2QR4Jsbz9sM2/2IReayqvqTkAhUVFdONUdoBishRwAcwgqpl4Gy7nA3bPRL4G8xOUjtV9RE9Y1J9nU2fN8x4Su7qEcD9VFXtwN4PrJ5RtGK+cCH78CP/N31hY1/cYRhgkvn1pssLXbOK2/cYxrdha+/eILs7wU17g5xCf0im/rBYNj/n7T53iNkObn1qD5AwIGqbkEwx2Ocktu36QfYFDhngfJCHLtMLGZ5jdYcG+cO8c+/S2+a2w4zh884ZU7DLLhlu6nTSLQvZYbhbnL972tJSgwwwZLa54Byle4TE6oY0QRqxK9wi8ApVvUREtgDfEZEvquqVroGIbAP+DjhJVX8mIoelOhORvwLerKq7bP5g2/9flAwmuSeIhx8CR3v5o4Dvl3ReUVEx/XBL4JK/xr5Ur3Nb6qrqbkyE+SOCZk8HPumizqvqDZkuH+8mP9v2V2CDARSghAEeClwlIs7F5EHAv4nIufaCTyy92Mjgf9EGkW/kZIAhE0z1kRtXKCvrsNP1naYLlkF0fjlfRtWAViGXbL+H2Y2B1ze5jOWMccPnF8oIIc8oUwj7D2V/obwP+pnfAAzw1ruZ38MxvxtshQui4TPAm6xAeactC2WCoWYeugbQy6EhdPiuzdOP1LPPMcHUtg45d7pCV7sWS+DtInKxlz9bVaNWIyJyDHACcFFQdW9gvYh8BdgCnKGqH0hcb0ZE5lR1r+3zICi32Sm5q9NKO6uoqDjw0NIMZqeqntjUSETmgU8Ap6rqrUH1LMbW+NHAQRjCdaGq/nukq38E/lVE3ocRkD0feH/pYBsnQFX9qojcHThOVb9kZ9hZS19XHksYRhV+sXJsLvzSunQbXcRYYROa3OWyG9qst03iTDC2nFjquMuZC+8LXeFsfa/toLUR3GrS7fM3AXAnuwduhzm5Rcad6MLV3WRTx6acg1G4x7C5ITdYdyPNCGWAKebnMyU3ziYGaJnffo8B3rTVNHYMb6dldzfYxo7l3dQxI+i23cXBNu2VCUZlgJb5dTZFCtlcyPz892lPUNbEBFNlsfxspK4xIOrotMAish4z+X0oYa6yAzOR3gbcJiJfA34D6JsAVfXNIvJ9sNGP4Y2qen7pWBonQBF5IXAK5tW6F2Y/trMws3NFRcUawKjMYOxe4+8BrlLVtyeafRp4pw1yugF4CPCOTLffxTAKtcfFKFkCvwQTVeEiAFX9j5xWZuxQeplVjLmFDG97Io0xwJRMJMZkUm3abGXomOCiYRY3LVp2t8063c90fyLHLrryJued4GzSjEYyZr/mGMuuGUP1Dr2b2SFq+6GG3m0+zNoL3uwNzTE/VxYyvzCF/k2I2jDAVJCCUOMbKwuY34LN79rcy/YAftXR8m7vqdsZsLyYFjitDTa/QVZWttF56FhZcM4LI2R6ISMkqPfrBtn+IMMAl1k3Sle4h2HiCFwmIi6kyWuwilZVPUtVrxKRz2OUrcvAu1X18lhnIvLfgLcAX8F4q/2tiLxSVT9eMpiSCXCvqu4zEzfYWVlLOheRk4AzMIubd6vq6YkbeL3t83uq+vSSvisqKlYOo1oCq+o36LrV5tq9BTOxNeG1wIOcplhE7gx8CRjZBPhVEXE+wY8F/hj4l6aTbNCEM4HHYtb03xaRcwN7n+OAVwMPU9VfrSqzrKioiGISfYE9rAvMZG6izLwPKJsAXwW8AGP8/CLMTuzvLjjvwcDVqnoNgIh8BHgScKXX5oXAmdZ2p8nex2IJE/DhTr3FvkDZLW0PT6TbghTST6KNAWqINkviRfNR3O+WrFaAvrSt++Lts3HqQnescFnmC+8PtZshOwG/M/HYZtUuh86Z+oPvZvLb7tbdNHnbXnO8+Ra7PHZL4BIlSNMSOLbcC5UgKTc36Cx9b9tq3vXdc8FSP0h/5f3Y3bqDo21zrnC77YuWdIWLGD939oPeaAzUu3uD2BuftYTIf0/CpW/nnL7uu0iZvbRxn4sqQUYnAxwDPi8i59P1VPtDzBxVhBIt8LKIfAr4lKre2GJgRwDXevkdGGGmj3sD2E1MZoDXq+rnw45E5BSMIgZjh11RUbFymNyQ+Kr6ShH5rxjZomDsDv936fm5cFgCvA54qe1YRGQJ+NvCLTFj6/xQdjgLHAc8EqNd/rqI3M+37AawhpRnm3HdXw0DtF9n9/X0GaAjQE0M0FecNH01cyyulBXGjFbDfneZ+1leNHTnFt99bt4GSpi3bMSyHhcEwQnxt3WNatjCXXvKDrbpPLt78tuCFGDLnGmz5bDe9KAgGIPvrje317rwLRieM5uw2fYig7E061LD5roRmXsDRTjWBf0sLZWGrC5X566zx+Z9c6Ju+KtENG7L/Pyo3w5zG83zceHQluxD6RhKb7TneMEy2GPf69D8JYxqXhKhPKyPla2gGcyooaqfwJjVtEZuWj8VM6s+SFV/DCAi9wT+XkRepqo5tTQYxufTtSOBX0TaXKiq+4Ef201MjgO+3eIeKioqxghFeuwbJwEispu4MtbtWXSnSF0fchPgs4HHqupOV6Cq14jIM4EvkLfLATOJHSci9wB+DpyM8fHz8SngacA5IrIdsyS+pmTgHYLpWNx2ryrF/FwbxxZjdx/K/NqEw2pyn2vjRtcZT9d9br89vsUGUti90TCVDZZhbJq34bFmuuGwQra2xTI/F0LL5R0jdHlzvKenbV/YrRgDtIEYNszttbfc7LbXDQ7r9tsweSdfC819/OOQtYVmKSG789uG8rswEKovz3NsLcbwUpgJ6G83v6+nL8cM93ntl2fthLOn+/ubxjaNBVAoCVcWy/vnRlY6k7grnKpuaW7VjJy2ZL0/+XkXvhHfqTUBVV3ELJ/Pxzg8f1RVrxCRN4iI8x8+H7hJRK4ELgBeqao3xXusqKhYLSwxU/S3GhCRh4vI8+zxdku6ipCb1vcNWNeBqp5HoJFR1dO8YwVebv8KsQ68r3kfywOz2PZTVxcGRC0I5DnQDnIlIbWa2GFMJtP5wpvvz7Izop7dbFNTfbPflzW+XWc1kI4tztn8QXO9bG7OY3Mh03MMMBWGyxz3Mr+ZII2h6+IX7ojXG/LfZ4Aha0uxxU75Xm93PStXdQFqO2HKFgtYzqz9QSxbc25ujt05eZ8PVzczEzwDK/JzTHPWY4B77XFnc4NFyzncOxx7P0pXHDlZYJIBTqYMUEReB5wI3Ad4H8Zz5B8x4rtG5H7x37B7gfRdk3YesxUVFVMMRVhanswJEHgyJqKMC7H1CxtnsAjJCVBVJ/SOZ4FDuqwuZHv+cUrrm2NXYZrTrDUxvVSQ1VgfDjl5TdPmS9E+LFucdcEXetniLamx+mWdNM4mfeaywdY5tjPTwnAypVF1jG3JY2j7HWtzmtOFwJYuZDQl2tKca1rnGazvSZc3Wjn8vGHH/sZXIfNzzyKUiy7OzPRew0NXU7y+d8y53ytkcSUucjkZ4LKwd2FkrnCjxj5VVRFxAZs3N53gY7IkmxUVFRMHVWml/FlhfFRE/gHYZgO3PB94V+nJ0zcBCjAn/czPZ4ClzC9ned/GgyEVfiusT/UTG1vOziu1IVTsK55jQk3oY7YBm7Syx1ayzRJNeJgvsXkbJAx8E3P3bUtTIdWsNYKz5VvqkSP2ispTGnFXvuTJCGcC2eKyC6SwEARSiNmyhnUl4fNzKx1tp/1eSajqW62L7q0YOeBpqvrF0vOnbwKsqKhYUagKi/snawIUkXcC/6Sq/8dOeMWTno/pmwA3YNheStML/XZ+IQuK+VlmZCA9fcVkLqntNp18KNAcRuG+sO4LHxtH6j6c88aeII2VDSLrbEKMAZbmS+oGYa1hvoS5h+zO/w18v/FYH1ZGtuT9xh3bQbvfp5P1ldhGOszYd2e/63c2YIA5mW3KLtBHURBfYXlp4qaK/wDeJiJ3Bf4Z+LCqXtpwTh+KoyZUVFSsUSjmA13yt1JDUj1DVR+K2bXyZuB9InKViJwmIvcu7adOgBUVFXksizEbKPlbYajqT1X1Tap6AsbT7MkYx4siTByvbcR6zHK3JLBByiQitgQOlwIlgvFw6Ttvl7zWPGS9MwmZTa/hnNC8Y9bhQiTlYkY2LYX9UBLh8ji1XI6JAEa9PI7lS9qUmBGlfq/cDoCpJW+JyVP4vDa6IBZdFuRMYmbcUnhmBP9uqfvzj/uW50F5zHi6KXTbIL/9CsDuL3ISxtX20cBXgf9Vev70TYAVFRUrCxMQcKJgNb9PA/4L8C3gI8ApdiOlYkzfBDiLCWoQY34OfQEFgjRnVlEqIPePrYnCehuMwLlDhc7wMey1Y3KC8uWOwiTibp0SWOeYbcgKU2nnWXScrzCWBX7l7UF+f5CPoFM1G+T9+5sNytYnyv3X1blDBkzZNelj516bnDuZ30fOyDiVzno78iV+/5nZuKzMNzXpBl9o+BfNKXcWgnzuvU/lYSInQMw+Iv8E/Jmq3tzUOIXpmwArKipWFkrvd3ECoKqPGkU/0zcBrsN8wUPmFzP+LTHydGjaS9gxzvnum+BcwjZZN6jQDSwH55zv3Mg6JrOLAVPKGUI3MUH/OJQBdpifC6nmwmD5wXhc2R02bWKEPsL/mJDRzkbqQqZ3UJD6fTh3z029bRZt+a4gXNogcr0SV8k+BthlpMt2bAsBm8vJhB1c2w4rDDWsbcx7SgzKU3kwE2B/jIcDAtM3AVZUVKwsJnMJPBJM3wS4jl72534YX+4VarZSMp6YBq3D9FzeMKR1luU5tgdw0GbDjEIn9zD0kwsltDeyt6rTFHa0h4uBU39OK9vERnJ1nWcSsrtuQNT+Mpe/I0h9thfKB0PE5Hkp5rc+yPtBft21t0TqvPKFILQ8NMuE27jeZQ3oHRs0v/t++xvvD0JpxbAcMr6ULDAnpwwZYFjulxG09VEnwIqKijWLOgFOEITmUadsmkK5SWwrTZduNyfPb7MbAW3uDQoKfrDPXga42AnwaeU4kWCSfVq+Tlgn4mmsrA1zGQghm0uld9BFEwN0bXNaYNeHF/g2iZQGOSY3TKAkoEKK+bnLh9tX9vTjWOj6nnR5NralRUMfud8zFQQhDKHl91ES2bNOgBUVFWsadQKcMIRfO99ZPbTpC7/WMSv6BPOb39y7WdCcF+IoHea9d8tEB3+THSf76wb2TMj+SuR5g4S66iClgfWPwzbuQjE2VzqIGENcH9SFlL2EEY4AJRr4lHwt10+f10pEPgnpshhyocjCMea2goh5iTgsk/YQmXJM7wRYUVGxMqhL4IqKijWLOgFOEBwddyN38QCP8do4gfTOIO9+xJhbm10euAAGXaWH2y3N7oORsQhN7Z3a2eVsobszmduRjD0be8eYSv3jUmVIDH1DzC0z3XHT239HpCy3pooOhPSSu0CRke13QLRRgpT0k9ojpiTIQwlSS99UeaxNDHUCrKioWNOoE+CEYAnjxuV+kPuZT/Hjjji/0+SH3AeAn/7IpOwIFAxRNyJjkpByU1oK9qqFrhIk3NfW7UXb2bt2r92jdk93X9vlXXY/jZDp5aI7N+0F0oYBdgiaM9Z1hsQ5G5qSOFWONcaMpJsQmrDk2GkqUEJMmVOIYQyhc2h6bIOwvVTfflmJK1wJk60MsKKiYs1imbiU4wDA9E2A+4AddGR/hxxunPc/ufcpnSbPmXs/AL/cfhgACzsPMRU5tx9rhrJgWVpnr9vN7rKGzcVMX0JXNxfowO1n6/pkl2d3k2J8OQZYGsy0xPHfDaXTv2NOsT2l3Um7g7Yh24NugARXVxAyq+86IYvbFJT7/YfnOCTCZJUgNtRBZH+pa6eYX84oOdXHYqSsyRwmZgjdxADLtzGZKkzfBFhRUbHyqEvgCcE+4Cd0vsQ3n34EAEe/6medJrffYtjBwtWW+YXa4NgXt/OFN4V7rHZ2jw126kJf5bDstLzOvW1PIHvM7dbWRgvchvk1KWFDB3o29bftY1lhMASfAYbMbxSywDCfaxO80iVGxqlnVCIDLJEJpthbyPhaae8z5U395xhgrL8qA6yoqFizqBPgBGERw+gcq3unSW6+9Ihum2NsGu4PnAuGkAqIutE5rgf7sYZjiqVtApU2aXhzdW1c4po0kD3nhrK3lNtcTgsc2hKWMMGUZjdmo1gY/KCNvV6JtjSl/Y2FmkqFqUq51fkIf5eQ1fn3FbZJuYC2lQEewK5wdVvMioqKZiwW/jVARI4SkQvsHr5XiMifRto8UkRuEZFL7d9pI7yTHkwfA+zAhm6/3mY/fmi36sgg3W7TMNip/zUN2WLqy5hjgGE+Jx8Ky0pYY9M5JeGwmjSPMdaTZHwlWuBQBhj7L3F1qbD5MQYYjiEhxBpE+1sSDKEE4bXbbFbUNO4S2WZKqxxjgDkmONol8CLwClW9RES2AN8RkS+q6pVBu6+r6u+N7KoJTPEEWFFRsSIY4aZIqnodcJ093i0iVwFHAOEEuCKoS+CKioo8nB1gyR9sF5GLvb9TUt2KyDHACcBFkeqHisj3RORzIvJro7uZXoyVAYrIScAZwAzwblU9Pag/Gng/ZnE6A7xKVc9r7HgWWHRLH7cl6C+79TuMATQ7rRnK4bY8tRSGMnOAtkgtjXN1g7hjlQRDaLsci47NGRc75YiLaBwLSpcKNBejEk22OjETl5SCZgAD6NRv4S93w9soWRKnrt1mt7bUY4xdw0WYXgziDObEHuF7H1s2t1sC71TVE5saicg88AngVFW9Nai+BLi7qu4RkScAnwKOKx5BC4yNAYrIDHAm8HjgeOBpInJ80OwvgI+q6gnAycDfjWs8FRUVA0IxotySvwKIyHrM5PchVf1k3+VUb1XVPfb4PGC9iGwP240C42SADwauVtVrAETkI8CT6F3rK10v/K3ALxp7XQ/cBdhpGcCCYwIRZuG+zs6FLPyKxr7wqSdS8gUsCS1U2l8bo+YSs5fU2EoE8SmlCyEjhHRU55ABDmIYHXOFC5Ug8exQyitIM73cb52KGp1ikyVsv+TdSjHBnBlMkxJkRK5wIiLAe4CrVPXtiTaHA79UVRWRB2OI2k2xtsNinBPgEcC1Xn4H8JCgzeuBL4jIn2C8bh8T68jKEYwsYeboUY+zoqIih9FqgR8GPAu4TEQutWWvAY4GUNWzgKcCfyQiixheebKqtthBqhzjnAAlUhbexNOAc1T1bSLyUOCDInI/VV3uOUn1bOBsANl0orKdrvyuYy7isYPwaxmavaTkKX5ZG/lakwirbV1Jfdu2TaK5HAsJmWCKGQJpI+bQDCa2l3Bq0GHf0M/8ErK/YQygSxhgnythBk1tY0bNTe9fzzmO8WlvGj6bnLH28DLALFT1G/0D6mvzTjouDuPFOCfAHcBRXv5I+pe4LwBOAlDVfxORjRhVxQ1jHFdFRUUbjNAMZtIwzgnw28BxInIP4OcYJcfTgzY/Ax4NnCMi/wnz/bkx2+sMhtE5zW5OtlMiBwrPScq7gvLwOJYPEfvCtxlj6a81KtaTkj/l3ML65IMhI3Sd+UbNTf9dMS1wgkQMosUvMVxvw95KkZNJN7WJyfOSARMCJui3S7mAhqjhsNpBVRdF5KXA+Zhp672qeoWIvAG4WFXPBV4BvEtEXob5zjx3XGv9ioqKAXEA+wKP1Q7QqrDPC8pO846vxAhF22EWbx9fm8Zs+lKIya5C17OmwAOpfmL5ErelcTBBH00aT4eY612qbc7esO/ZhIzQhysLv31ZUVH8uqk0hpQ8NGRZsbJwQhhEZpvrI8UKc/sRd8YbyAI7/dv8Ru+5un5i7qEOdQlcUVGxZlEjQk8YFumO3DHBw716xwpTWt9cqKlUGgtqmpIbhvUxJOU1mXNSbXLXaetxEg0SG/SR8iKJlZVo0TsYIeMbhAmGyMkAmzxrfDTJ+lKMsASx+yplggCz0ntOSQCFAwjTOQFWVFSsHGpA1AlC+GO4L9Y2r2x7kIYsYBAGGDLBWD+DyMxGgdx1Su0a2/i05rwSmpjfMM+gjSx1EC1wDjG265cTqW96tsMwv9j1kmNJMEG/Tc4XuCpBKioq1iwqA6yoqFjTqBPghCNm3DlvKP+6eROleMbu9bu0aPbxXbYpAHYXuD5lR04JUhqZeZjABsO66zUtdVOmID6aAirkkOt3GIzClbDk3kM07bNRslQcxdI3htRY2iyFqxlMRUVFhYdqBjNBEMxXqsQUZLHXrGJmdrEnXVrsdrJk2WFnb1+7G1x/0AWvw7YMsMSNrg2bK3Xfi12vpH4QxcUw5j2DtB0FBmGC4bkxJckgCqfS65WMMWSEUSWSC5wQmMP4qDLAioqKNYtlioOdThumbwJ0wRDCXdx8BF/a5T0mYOeClfmtm83weVfX+RIG+wHHzA5KTUzamKfkzgnrBgnvVSKvLDVhKTFPaSpv26ZUhtpUHmvThgnmZLajlPHFDNVhsFBovsF5n5F0oo+6BK6oqFizOEBDlEzfBDgLPQFRY0ywT9Ziv3iL5vPZibYa04YFcsNW7mxtQhc1aWVjaGJzOTllKm3DUlNjbBO4oQ1bLEGTLLXNuWF5bnWRysfc50YhAxykr1T4slgQ1TWKui1mRUXFmsX0McAZjNtbuLVlTHuVYkgd+LKQxPWGlQP5+Zh8qInFtWFzOXe91HWGsVnMsbmmYASzieMcBtFUx9jcMLK+QQJRtAlpNQo0Ba7NPfvpmxGGwhq73YqKivY4cNXA0zcBOi1wjgGGX/YS74ph0MT4SgKwljDApiCtw9gqhmOPoQ2bS200P4wMMGdH2ebcYTxbcv2GfaWebVN4sUFRyvxaM+kD1xVk+ibAioqKFcaBawldJ8CKiooGVAY4OXCucLndrFKC6mEE/rn60iVwG4VGLvhCmyVwm3sPURpwIPbsm5bCba6TU2SM0t2sjcnTIH3kTGYgH9svBX9cpUvfnKlOFHUCrKioWLNQqhJkUiCYUeeE68M4tTcxpRIBcpPpSaysZBe6sE2TaUtsvE1GzG1MW0J2FysbhgEOYq4SjiPn4jdupDlHvvAAAApXSURBVBhfU1Rpv80g0aSbomSXMOkeVBlgRUXFmkVdAk8nSmVLubrcF7LUJS3HAFPMr41R8yDyvRAlzK+J3cXKRmEGMwpzlRhS5jDDkJ2YuU2ujY8cExzEBS5kiTmXzMoAKyoqKmKoDHC6MIxWMdU29uVNMbImVpera2PU3Ib5NTGynEY31LQPwwDboMTtbJh+U8xvFEzQxyD9DGMUnXr2bYK29qAywIqKijWL6go3ucgxlxAlX/gmu7I2mtySDZWamF/Mfa7Nx7hJI5jbFzZVl9PstgmCkEKb4Aul8sGSYAg5F7mmFUKurolRuvJYQI9UuLRh5LwLkbLqCldRUVGRQl0CTxbcyMPAqEAn0KkL9ljyFS1lfjkWNwgDbGJ+uTE6lAQnaGJ8OY+aNrZ9gzA+h2FYTYoJxtjXIDLAUcoFS7TAgzDo1HsRMr/qCdLB9E6AFRUVK4Q6AVZUVKxZVC3w5MEteQ83X6ZDjryhU7Vkd3/bvWsLAMs7N/eeGzNpGWQZm6prYwbTZNQcQxtXp7ZL30GMmocxbh62TXi9cRlNl2IQ98uS5XUsmnMKqfejTbCMHlQtcEVFxZpFXQJPDhwbtyM//O7XAvAQLuo0uXbmKAB+tO1YAG5xDNB9+XYFKfSztyZ2F6sbJExVSkGTwziNmXOhrYYRyOfQ1HYYs582jHAYRUfuOqX3V2J249DGUNq9n7HfMQy6sMYMoce2K5yIvFdEbhCRyxP1IiL/v4hcLSLfF5EHjGssFRUVw8AxwJK/6cI4GeA5wDuBDyTqHw8cZ/8eAvy9TfNYwjA3y94Otgf/i9d1mryR0wD4EYYBdr6AO216vU19BuiOS2WB/vEgJi1kyiAvV2syRykxZm7qK3fdNhi1rK8t2hhPj1p+2HSdNn2kMKrAr3sydQcwAxzbBKiqXxORYzJNngR8QFUVuFBEtonIXVX1unGNqaKiYhBUJcg4cARwrZffYcv6JkAROQU4xWb38tbusvqqk036mz1nPGukAx0S2+lyz0nHNI0Vpmu80zRWgPt0D687H16/vfC8abrHVZ0AJVKmsYaqejZwNoCIXKyqJ45zYKPENI13msYK0zXeaRormPG6Y1U9aTXHMk6MTQlSgB3AUV7+SOAXqzSWioqKNYjVnADPBZ5ttcG/BdxS5X8VFRUribEtgUXkw8Ajge0isgN4HbAeQFXPAs4DngBcDdwOPK+w67NHPtjxYprGO01jheka7zSNFaZvvANBjBK2oqKiYu1hNZfAFRUVFauKOgFWVFSsWUzsBCgiJ4nID62r3Ksi9XMi8s+2/qIGo+uxomCsLxeRK63L37+KyN1XY5zeeLLj9do9VURURFbNfKNkrCLy3+zzvUJE/mmlxxiMpeldOFpELhCR79r34QmrMU47luquqqoT9wfMAD8C7glsAL4HHB+0+WPgLHt8MvDPEzzWRwGb7PEfrdZYS8dr220BvgZcCJw4qWPFuFJ+FzjY5g+b5GeLUS78kT0+HvjJKo73d4AHAJcn6p8AfA5js/tbwEWrNdZx/U0qA3wwcLWqXqOq+4CPYFznfDwJeL89/jjwaBGJGVePG41jVdULVPV2m70QY/O4Wih5tgBvBN7McBs0DouSsb4QOFNVfwWgqjeweigZrwJ3ssdbWUXbV1X9GnBzpknHXVVVLwS2ichdV2Z0K4NJnQBTbnLRNqq6CNwCHLoio0uMwyI2Vh8vwHxVVwuN4xWRE4CjVPUzKzmwCEqe7b2Be4vIN0XkQhFZTa+FkvG+HnimNQ07D/iTlRnaQGj7bk8dJjUeYImbXLEr3ZhRPA4ReSZwIvCIsY4oj+x4RWQd8A7guSs1oAxKnu0sZhn8SAyz/rqI3E9Vd4UnrgBKxvs04BxVfZuIPBT4oB3v8viH1xqT8j82NkwqAyxxk+u0EZFZzHIiR+fHhSKXPhF5DPBa4ImquneFxhZD03i3APcDviIiP8HIfs5dJUVI6XvwaVXdr6o/Bn6ImRBXAyXjfQHwUQBV/TdMILLSQAMrjQPfXXW1hZAJ4esscA1wD7rC5F8L2ryEXiXIRyd4rCdghOPHTcOzDdp/hdVTgpQ825OA99vj7Zgl26ETPN7PAc+1x/8JM6HIKr4Px5BWgvwXepUg31qtcY7t/ld7AJkf5gnAv9uJ47W27A0YBgXmy/kxjCvdt4B7TvBYvwT8ErjU/p07yc82aLtqE2DhsxXg7cCVwGXAyZP8bDGa32/ayfFS4HGrONYPY8LP7cewvRcALwZe7D3bM+29XLaa78G4/qorXEVFxZrFpMoAKyoqKsaOOgFWVFSsWdQJsKKiYs2iToAVFRVrFnUCrKioWLOoE+ABBBE5SkR+LCKH2PzBNj+W6DMi8mIRebY9fq6I3M2re7eIHD+i6/yBiJxmj88RkacO2M+dReTzoxhTxYGBOgEeQFDVazEbzJ9ui04HzlbVn47pemepqtv4/rnA3by6/66qV47oUn8O/N2wnajqjcB1IvKw4YdUcSCgToAHHt4B/JaInAo8HHhb2EBEjhGRH4jI+22ct4+LyCZb92gbq+4yGy9uzpaf7sU0fKste72I/JllZCcCHxKRS0XkIBH5inOfE5Gn2f4uF5E3eePYIyJ/KSLfs4EM7hIZ672Bvarat9+siLzRMsJ1IvITEfkrEfk3EblYRB4gIueLyI9E5MXeaZ8CnjH44604kFAnwAMMqrofeCVmIjxVTVimGO6DYYf3B24F/lhENgLnAH+oqr+Oce36I7ukfjLGrev+wP8XXPPjwMXAM1T1N1X1Dldnl8VvAn4Xs3/9g0TkD2z1ZuBCVf0NTOzBF0bG+TDgkrBQRN4MHAY8T7uBBK5V1YcCX7f38VSMC9cbvFMvBn478Uwq1hjqBHhg4vEYF6f7Zdpcq6rftMf/iGGL9wF+rKr/bsvfjwmaeSsmLuC7ReQpmF38SvEg4CuqeqOasGUfsn0C7ANcyK3vYPxSQ9wVuDEo+3+Bbar6Iu11ZTrXppdhgnfutsveBRHZZutuwFuqV6xt1AnwAIOI/CbwWAzzeVkmgGXoA6nEwx9hJ64HA58A/gBoo0jIBand701gS8TDs92B8fv28W3ggU7Z48FF2Vn2jl3e9b3R9llRUSfAAwk2IvbfY5a+PwPeArw10fxoG48OTIy6bwA/AI4RkWNt+bOAr4rIPLBVVc8DTsUsZUPsxoTSCnER8AgR2S4iM/ZaX21xW1cBxwZln8coeD4rIrFr5nBvILoHRsXaQ50ADyy8EPiZqn7R5v8OuK+IxAKwXgU8R0S+DxwC/L2qLmA2qP+YiFyGYU5nYSa2z9i2XwVeFunvHOAspwRxhap6HfBq4AJMBJRLVPXTLe7pa8AJ4XYHqvox4F2YWIUHRc+M41HAZ1u0rziAUaPBrEGI2UHvM6qakxFODETkDOBfVPVLI+jra8CT1O4hUrG2URlgxTTgr4BNw3YiIncG3l4nvwqHygArKirWLCoDrKioWLOoE2BFRcWaRZ0AKyoq1izqBFhRUbFmUSfAioqKNYv/CwOdYJNow71vAAAAAElFTkSuQmCC\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": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASMAAAD8CAYAAAA8P8JjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztfX3QZUdZ5++ZmXzNkGHySUKSMmEriiyuQmURjbVLEVhDROJuiUtENqXZSm0VrhG1IJHdYq3SFUtLwJLCnTJIXNkENuCGYlFgYxBdJcsEWAkEhAVMJgnJDORjyPdkev84/czt93efp7vPvfd935PJ86t6q+/p06dP33PP2/3r51NSSggEAoHNxpbNHkAgEAgAMRkFAoGJICajQCAwCcRkFAgEJoGYjAKBwCQQk1EgEJgEYjIKBAKTwLpMRiJyoYh8WUS+KiJXrsc9AoHAkQVZtdGjiGwF8PcAXg5gL4BPA7gkpfTFld4oEAgcUdi2Dn2+CMBXU0pfAwARuQ7AxQDcyWiLSNoG4Jh8fEIuTz1u1ubBR4byQD5+nMonqQSAQ1QKlVvouKdN2dY6tupq031y2nB96jjX6qN2/zGwvvMqr1nPMfW0q7VZ5LszvO/X81srtjglAGzNJf9vHAL2p5ROAYALL7ww7d+/v2u8t9xyy0dTShd2NV4C6zEZnQHgjuJ4L4Af5EYicjmAy4Hh4Z0C4Nx87l/m8ornztp//LND+Zf5+Gu5vCuX91MJAA/n8tFc6g92LJVHFddw3dF0rOVWOi6xlY6frBzrRPkEnfOOy8/6sh1yjrkPayyHMB7e3p6/d+85Bo+xB17/PFarndemp23r/tZ34Wfu/eZlnb7Deh99T7fnckdxzTNz+Q+5PDArtQr79+/Hnj17nFGvhYic3NVwSazHZGQtHnOLQUppN4DdAHCsSNqO2YPVCeUjn521/yZdrz/GM6m+nBx0Mnool49SqT/SsZiBJyqenPi4fHG2Ul3tpeZr9CXXF1X7f4KOrTo93tI4tsao8P5BxqDsg/9p+ZmsArW+vAnGmky836nWdhVgNl9beFqTUPl/sCuX+3JZLtAzJAAH+we7AViPyWgvgLOK4zMxIzCBQGASSJgty9PAekxGnwZwroicA+BOAK8B8NO1CwTDdkhXhntyeXPRRs8dwFroNkpXinI1U+agsqjHcqk/ATMmYLaKeAyJWUnJqngr522BrNWVGcsWp76s460iM6XHqV15zttu6JjLa2rbzPWAxz4WYUJ8vvbse7Z0y0CfrceEuCzb6Ht2DB1vp+Py89Go4WnAjFJKB0Xk5wF8FMPv+e6U0hdWfZ9AILAMngaTEQCklD4C4CNjr9MV4Vu5fLw4xwxBz+lqZgmSGbwCWtcye+LVSu+rq0/JEpgt8crUI4Ng1sNMqTzHTMj7fuU9WNOixyw8tdjPGEbU27bH0G0RWY7HhGpC6UUYWYs1HjLOecoJS4DtMXRPllmOt/5snyaTUSAQmDpiMgoEApNBTEYulNKyTQUwTznZoHErHZd1vF2qbSGY4rKw26PWwIwqs52PZ7Nk3dertwTKLfW1dS238cwAymvK7XJ5/0UwRsXf2nqtl5q+Z2y927PHjbqWTVnN1ITfIT2u/a/YOISZUGIamNRkFAgENgqxTesCq0AteCtTz6rNLKqmAmUBryfYtupYyGh9L2YKNVbjjU3BpgRWH55wu/weVl/APCPr+Z3GoGVF7T0j65zX56rV9YWrxZpjrgfmlSCeQWPNXOQYp75k2z3vzoDVTEYi8m4ArwRwb0rp+XTuVwD8NoBTUkpV/5MIIRIIPC2hzKjnr4n3AJjzXRORszA4zN/e08lkmNEW57MHXonGrNI9xnOe+lyPmSGVn/Ucr3RcX/bHMoAe+QhjETMHLi0XEu7fw7LPntvw9/HGavXbIzNqocYOexlR+Rw9RlRjOZ7MiGVH4/0CV7dNSyl9UkTONk69DcAbAdzQ089kJqNAILCROIT1dAcRkVcBuDOl9H9F+mIdTGoy6lm9xsiKuK7HME6v8eQ8LA8qtSWPUR27ltRcSLxoAVvouKxTeFECVr0H9+5r3a8VDaBH7tNib4vIjHqwiHGnx4jK96PFiKwoEh5r4tLaWVjhRdaimxmdLCKli//u7OhuQkS2A3gzgH/RewNgYpNRIBDYKIzapu1PKZ03ovN/BOAcAMqKzgTwGRF5UUqJA3AcxmQmo0OYtx2yVldPVmTBY0I1mYoXSoPvx/v/8jPLkzhkycOY4Viq49XRkiN4cqUxLGAMa2pp3GoMs7dvqx+Wh3jnrf5WJSMCbO2n9x56GrMSLZZz9Ii2NVuy+nuxfqr9lNLnAZyqxyLyDQDnhTYtEAgYWJ02TUSuBfC3AL5HRPaKyGWLjGgSzChhWEk4lGyPVqjmFLoKtGQCtah8vEpqfWn3qufUwtuTH5RyJtakeAxpkZXGYo/cj2e9XQuu5qHGjFqMqGZ1vIgVdQ96GZHFsrwooWxFbbFgj+3UNIsbxYxSSpc0zp/d088kJqNAILDReHoEVwsEApNHuIOYOIRBgKtCXMv571iqswS7vfBU4d7YyjaeIBvwY1Nrm6PpGGgLu7WvcmunLgGtZAE1FThjzPbWi7lkxe7xjq37tVT7Yxxl+dj6jWvnynpLgO2V/L5YW65Wsgdri9wryC7Rjme03jE7x2ESk1EgENhoBDMy8SQGJuCptVE5t2pVrrfS8arMbKdErwC2rGNWxQaUD2EGjjbprbTW/XpDeNRW2ppKn+/jHdfu0xs6pOdarq+lbfLMRWrRGr0ojYuo6WtmHGPYb78rUUxGgUBg07G+7iCLYBKTkWCY8XX114wH1mqqK8MOOq4xI08dW1PD9l4zJla0tTovIvPSVVjHwEHdeHW27tfDbuCcY8bVk6hQ0RPYjtGjtm8xvR65ILe1ZEZskMttaoH0eo0dLTblqf97nIZtxDYtEAhMAjEZmdiCgQ09kI9rJvScP4pdI6z9Pa9mbKpfSx3NsgDPcdZqo6ixqVboE4u58DUcFpYxJoV1zVDOY59j3Gtq6NG4jUVNdtUKQ1NLD+6FA/FkPOXn3hTqVh3vBHpCiPi/QUxGgUBg0xHMyARnlLXQMnG37H5YI8XuGI9SCYxnRjWNnHdssbdW/2NsoWqhIzwZUU3m1uuQ2yO3q6Fl9bKI3IePa8zZY7AWc+ZEED2asRYjGhNCxHMTKcftHQ+IySgQCEwCoU0zIVi7qtaCmT9Kbfia8vEyA3qIjmuOq579SGsFXhZj7tNjd7Moao6Xi1hCM3o0cPwsWBtZ9mHJDK0+LJkRM6IemaL20wr5YgXSO5rOLcKmen7zWuDBAcGMAoHApiO2aYFAYBKIychEwlp6rZ9LF4gtRl0J3saVn7l8mI4t9XxPRgo+7t3C9Wzt+D5johr2RAFsHVvXjImftMw2jQXxy+TRqykcWtuzWp67XgNGy7Vjma3dmO2Ztc2cISajQCAwGYTXvoknMZvtNatGqab/di6ZBfSoYVlnwCvQ8ca5ReJKewLDRVb2RTAmiwaznFoeLs9Js8dZ2Xt+NWakv3tviI9aXS3kS0tgbb1THvNhVw42ygV8RlRzEO+Nd24J5tn5ei1CmxYIBCaB2Ka52IoZQ7H29zyH1wwJFbqaPDOXrXAM5TUtZlRblVvq5Z7+FmFKnvzHkuV4LMdanb2YzWOydShqTrzsvGvLOubbcz/luVo2l7HhQID+AGk9jrIeI7Jy5FmhakqU35+/q33t9CajyA4SCDwtsdLsIO8WkXtF5Nai7rdF5Esi8nci8qcisqvVz8LMSETOAvDHAE7DsADtTim9Q0ROBPA+AGcD+AaAn0op3VfrawuGVUIZDOcXB+blB6whq2sO1vZRM9pbJhSFgrVANZbFmiPFIqFFFD0yI15x2fGyFhRsmTAW/L3K58nZRjzUfuPecLBlP+sZIG2RECLl92dW42kYLQNhZUjHwMPKmNF7APw+hvlA8XEAV6WUDorIbwG4CsCbap0sw4wOAvjllNL3AngxgNeLyPMAXAngxpTSuQBuzMeBQGBSUAF2z18dKaVPYqZj0rqPpZR0tvsUhqyyVSzMjFJKdwO4O38+ICK3ATgDwMUAXpKbXQPgE2jMiNsAnAhAeRwzJGA282ug+gepD9YgAPOrhrfyjgkstgx6Zv6e+3n9eAzMWmm9EBRj2MAq7Ix6Ql949kc1eGFWLG1rixH1hAHpce1oObnW5EPed7bCyPD3moDM6Ocw7JaqWIkAW0TOBvACADcDeFaeqJBSultETnWuuRzA5UCNRgYCgfXBqMnoZBHZUxzvTint7rlQRN6cb/TeVtulJyMReQaADwD4xZTSgyLSdV3+MrsBYJdIOhHAKfmcatXKfbAyIl0BHqZ6tqouP/do3hQt1uEdl9f2ptyxznn3rdnwjLlfiwn1aNM8WVHN0lvhWVfXwCyqHWh+BmZXlmyl5fQ6RmbErKe8tqW57Alo59nV1ZzKfclx92S0P6V0Xm9jhYhcCuCVAC5IKaVW+6UmIxE5CsNE9N6U0gdz9T0icnpmRacDuHeZewQCgfXA+m7TRORCDOKZf55SerjVHlhOmyYArgZwW0rpd4tTHwJwKYC35vKGVl9bMQTY15VJZ/Z9RZu7qE4ZkfqqcXD68rOnOepBK4ljKYN4lM71yCC8sdXsSloMzGM/tXNj7GO8+5TotS+qacZYG9myOyqv8bRqpWyl5WfW8yy8QPls4V62GRPyxbNZq1mJs3zM1syubjISkWsxyIlPFpG9AN6CQXt2DICP593Sp1JK/67WzzLM6HwArwPweRH5XK77VQyT0PtF5DIAtwN49RL3CAQC64LVuYOklC4xqq8e288y2rS/xhAXzcIFi/YbCAQ2CuEoOwfN+q0Osntz+bmijQq12blVy5qT4RiBIcOjw5YTLtexQSZHnATmt6Y6tkN0XMIzTeByTOTAmkrac6Ltz1w6v33qCYHhbekWcV62BNheRo+ebB1em5o70TJGfd62kw0cy7r685qeO8gkJqNAILDRiMnIhLIiFU6r6P2qoo3W3ZHLb1G9wmIBXForHsNToXoB28rPbGagTEjtqcqY28yeSgF8iXJ189ThLUdMq67FlKx+l1Gx9wiuvdjXnPWkh13xcVnP373HabjFJGuhWHrD0tScsD0Xllqc7vUWYK8Kk5iMAoHAZiAmozk8jkF1r+r7N+TyTZfN2nwqy+b/LB/rrM8q8HIV2J7LHXTsrYglWoxI2U0p/1FGdCwdH0XHJZvaQnWefmOL8xmYZw4eQyo/e2yxlrurJ3gboxbczGurvym77yzimlMLbNd6bmOMHr2wKrXfrYdR9jIiixnVn1cEVwsEApNAbNNMHARwDwZnWWAmO/rzwlLha7lUjZu3L7bkI8qIVPO2g9pae2o2amTmoozIYhAto8TaiqjsSRnEIkaQHPi9ZEbMDj2G1BNCZIz2jN1AamyUf1OPIY1xzLXaec9tEaPHVkiW8j6LhOL1jDc5W3JZ12SUKVT7gUBgCuhxDNxATGIyOoSBaSgz+kYu/7Roo6sju4EwLG2UZ39zLNUD8w6ctbz1Hlq/sRVqlWUr+v2OpvMWWlqg7UVb/c476FhLK5B8j2Msj7FXe2alAWJW6tlR9TiW8nhqIXg9txqLbbdcb3oC23mw0iktIjOquj+pcd+EMInJKBAIbDAS+pz9NhCTmYy2Yt5B1nL1ZVkRr47WKuatIotoZ9gquGQQbCvkhTy1AsDpCsgrbW3xYtbmyYpKtuhpFj3ZkTWmHplHy1YIdL78zL9Pj8xtkd+yxYQsbaQnKxpjg1VjlHzMluOeB4BlgV1lYMGMAoHAZBAyo0AgsOkIZmRDsJZSclygGlgNbLlnqDCYaavlgsHUvZZdgvv06D1vxcr78RiYYntZakt4xnvap5XVtFXWIj16sITDnpq5ZhTYuz3r2Zrxc6u5Z3gq/po7TW+0zR5YIhxPlb+80SNiMgoEAhNAQmzTaqiZy3sCXa4vVwgWgLMQ0BLwtgz7etwLPCE3s5+yjvPL6zU1g31e+VqhRAB/5a6p7RcJfeHFHa8xvfX83+A410DbHYR/t/Kcp+Lnvmtj6Tn2BNheeBpgno2aSKj75mwCJjUZBQKBDUQwo3mwyYOuSDuKul6GUvbzELXhVYQN/co6lhvUVNIMj21YMghPvcwyI0sF7m35vUBc1rVWbGi+X8utwBqj5+Cp93mMjsvPXj4zL48af67BklF57h+1ECwt2ZEFLzMKsx7rOfKzYTcQS2akYzTf1RBgBwKBySCY0TwEwwrDoVat0BcKL+xmLS+Wx4zKvnXFZreIMaFqFS35TNlva4XtcZvg79cTIkXb1phgbzAwa0Xn0CE1LRCv9o/R8UPUrpYd1jNwrWnTFgnBMuZ34/fBE9mU1/S6gdSYkUmAVsiMROTdGPKj3ZtSen6uOxFDFtmzMXh4/VRK6b5aP8uE5Q0EAk9V6GTU89fGewBcSHVXArgxpXQugBvzcRWTYEZbMWS91JXPyve0vWgLtIOflecO0TldkWo5p1qhNaxZ3HN6ZSyyOvfYT3m2SeW1HPrECwo2xoanJ/SFJ/8p2YHHojzNUakt9cIB87tUCwfiadHGMCOWB5X3s4Lnl+iRvXmyo5r8zGVGK/JNSyl9Mqe4L3ExhlxqAHANgE9gSOroYhKTUSAQ2ASsrwD7WSmluwEgZ5c+tXXBJCajrRjCh3w1H+vqXa7onKpI0WMqwauyZwNTnmNNBwf6GhPMf4wVtSevKFmAhlHx2KGerwUHGxM6VtGyHVqkrfXsFa00UeVv78mZ9PuxY3D5uZWwYQybUli/ea9WcozMyMps3MXQxxk9niwie4rj3Sml3d1Xd2ISk1EgENgE9DOj/Sml80b2fo+InJ5Z0ekA7m1dEALsQODpCGVGPX+L4UMALs2fLwVwQ+uCSTCjowCcAuCBfKzxjHYVbXib5tHjEi2jwJpjIh/35MXia7wtnwUvJpGivI9uQb6dS31ui8SqVvDYaoaSVhurjx4sshpa4+FtLm/L1IDWMnD1DF2td8yLk92T+cPaUgG+qUb52RNc81a2534AVuoOIiLXYhBWnywiewG8BcBbAbxfRC4DcDuAV7f6mcRkFAgENgErMnpMKV3inLpgTD+TmIy2ATgJsxVeM4DsK9roSqcr0s5ceu4UgJ+dVWExJP59WGDNx4uowHtWdP1evLKX13srLPddy4jhZbX4m5SwClwksmZsNeEqG656DHNMOBU9rrl2MGtiJ9ha7jMP1v+5931qChZ+Jiy4rj0Lz9UHQLiDBAKBCSEmo3lswSATOiUfKzP6ZtGGjR4VNdmKroI1B0s+9s4xG7GYUe+qWctqYTEhHpcX97sVy7nsX+veNpiCYN5c8HPFVQcbpY+PpL/Kn/RVG0bwn+UFANayOzZU9IweLXeQFlOo5VrzHGPHsF9P5tbjsuKZLpT9tsxTSuh3tPo7jIhnFAgEJoNgRvPYgmFVUvaj+dPuL9rck0tepVTLVssKW96n7MNaGDzNF7tp6O9Yk8d42qxyrMfQuVY+rvJzK3SsPs+SZem5f5X+IH96by6/k8vH6BiYMaBH6biGbU45jOBX01UAgP8tv3n4Cs8xlt0/LNcfj1VZDtQePPmjxYwsdxbrvOWe1HJdsdyTPHmTBXaZMgP0RaqiQCAwCRyJAmwR2QpgD4A7U0qvFJFzAFyHgeB8BsDrUkpdFg26EulKflJx7s5cqhyJZ329xgok31pNai4JDCsAFsMLB2vJcLyAXhz4rWQ3zHiOb5TfnV5XXK1887/l8juNEphjRCl/e099BwCiPGMtI5qVzwAAnJ/+SXHRUPeI/A2AGRPiLMJ6XLrIPJjLx6itxz7Kzz1sg9FiKpZcy2NEXFrOw6330nK+ftQ4twYTkxmtwgL7CgC3Fce/BeBtOXTAfQAuW8E9AoHAKrHaECIrwVLMSETOBPBjAH4DwC+JiAB4KYCfzk2uAfCfALyr1s8hDLO4TtS6dh5ftFGWtK+4Bpjfh5cMguVIY+xVFoEnc/DYDjCf4bXFdoCZZfrxdKyytuP0Bo8+J38oLPFT5hBMN3h5fmx2SdMASLFmec6Nt+SLj87lMfn+x2ZXJcNz9bh87+N+YChP0rHmMj2w9isAM753wCn1a5ZsSus4BK6n7QL6baAs+U8vM7LYG8syayyC30NXm3aEbdPeDuCNmP1fnATg/pSSSjj3AjhjyXsEAoFV40gSYIuIhpm8RUReotVGU9OUV0QuB3A5MLCCA5hPYFiyAZ7E91G9tRKxJS1jkYWhR/7Dlrwe+wHmmc8uKp+Zy1J+pgxo51FUoeUtubz7a0NZUgiPKng0AZjPqukJSiwjHs/pjh8SMP+gPpDL/5jLzIgk06Cdhbp1pzrq5br7Dq255PDXLjW0LfZkae1aDMgLmF+e80io9Q637IoshlQLzLcGE5MZLcOMzgfwKhG5CMMrtRMDU9olItsyOzoTwF3WxTkeym4AOElkNb4HgUCgDxPcpi0swE4pXZVSOjOldDaA1wD4i5TSawHcBOAnc7Ou0AGBQGATcCQJsB28CcB1IvLrAD4L4OrWBU9ioMqecyMw28Ywbf021qJ8diw49rZr1ozsZWXtyRzRK5QG/O2Y7rjURWZnuR88hUpu/OVcqr5b9ypAe5vWE0y8S4CdS7YW5Xgd2zEDay70XBZkH95jfSuX5ffat/bcCd9eWz74xNquy268eNYKK/IiGzV6W7CaANtzd7EE2AweY83NxY2BfQRt0w4jpfQJDAG3kVL6GoAXraLfQCCwjpjYNm0SFtiq2ucwIZaLh64U7DpwwGjrqUNrZv6eer7lggHMMyE99lTxwIzUKCNi0rOVWY/VSMvX5vKOXN5PJTBvMcgMiVO+lnVjmJGCHzbTkFKazxae+sCeRaXSYaU2wOxB7qMy97EzHx9bsKlWgDQrBIeX6bfHZYWZUI8pQa8Au5ZxxtwRHEnatEAg8BTGBAXYk5iMdJK2nEIVXlAwXnjLlYjdQTikSO0+fOwxpFLk4cmGetT0THLk5PzhNDph1b0yl+pNzIzIUu23ZEUlM2I9ds9L7EU54wdYjo0fqp7jB/pddAzMHn4j1cfRBU04hQSO/DU52y4wzzK8HG/s7Gu1YWdbyxiXHzW/n7UAek3t1JEoMwoEAk8xBDOqg1edcnHWhfxBOuZtbznZtzQR1grCuea1rZdl1JIZ8ULuacqAEYzoNMzAF3laJrbqA+Yf4BifhFXEf2VLPEsd2UoW9kLqy7qPpzotXhD1UNn1nbVDUjZjhapltsFuSWzsWDJ1z+1kGePbWtA47xjAJCejSFUUCDxdscJURSLyBhH5gojcKiLXikgrBP0cJsWMdMXQlanc0qty5AC1rS3SHLeEw3HWwol6K0+PnRGby9TsjPSzKH1ifxAWOJWftfTkQFZqXlZD9kT48mKe1l5U9u1hAZ7VRy1fUoldxnlWnT5OpUVV8vPZkZ/b8flafZxWFtqWaxHLhUp2z1o0793t0fKOyQjs2hmtSJsmImcA+AUAz0spPSIi78dgCP2eMf1MajIKBAIbhNVv07YBOE5EnsAwj5tuYK0ONh2aT06ZkJrJlBO3yll0leJAW4pauIcexulZtrYsssuxeeFfrWD7x7G9DVsfM70qP6sQSmVGLa9NwGcmNRuiRXIG9aJ84DpOfugcJexPcqlOR2Ubfm6P0nEZQ4S0dTuyrM1L6lgOwbPS9mRI5eeeAH18v2Xg/kwrmoxSSneKyO9gSNb4CICPpZQ+NrafkBkFAk9HjEtvfbKI7Cn+Li+7EpETAFwM4BwAzwawQ0R+ZuyQJsGMAoHAJqCfGe1PKZ1XOf8yAF9PKe0DABH5IIAfxozDdmESk5E6yqrN3vNz+dKijVK4HKHnsC+obkxVq23tLtiEvweekWXN1N4TdntmAQBmLhAc38dL4lV+9pK61YTEq5QTsFB6EVgJ4VgIfRQd6/eseSlzAG1OG1tek8tjaZumTa1cfIxWBMhy2J5K39qSeWr6HsF1dRu4WneQ2wG8WES2Y9imXYAhLv4oTGIyCgQCG4wVCrBTSjeLyPUYEnAcxBCtY/fYfiYxGT2JQQarfpBvyOWP/8iszYN/PZTK+9SuT4Xe1Xi/xv1KLCMktPKwKxNqmQWsacQX1dLCMjNab7P+lgGhxZBWIXlVeILzmrUq01I+LutyuTVfe1S+j+We5OXCaw3VquuJZ70KuK/HCt+blNJbALxlmT4mMRkFAoENxgQtsCcxGR3CsMXXGFo//uz84a8+f7jNzh/9PgDAaVlhyOFGFt47Y/1/k6p5vueb0hMHwjMSXGTF8zwwx/S3Sja06P1athlW0HJ6xluJGS2DnneLH++GPcaYjAKBwKbjSI30uCy2YNBgaNbYj2cV2ctf+n2H2zx001CqWwh7LyyyAm0qI1LUGEkvJrbCTRq1+BwjUHO2ttpNDmppPCFMYjIKBAKbgGBG89iCwSL/q/n493L5+Zvm26p9kdokuTmhqP8aLDOcVXg88LVmXx5dG3PjVSy/PbTRk0lZbVdJCTwV1iKpXq0fm8pFfgJFTdvWsg2y6j1buaU1cCHADgQCk0DIjGxswxAVQn0Zb87ll4o26hOqk7n6hqrDrGVf5ClWaiuft5D2WNh6bbgst+rHeVHZ2eOy/IKegUrLCMZCz+rotaldq+eWYUieRlE1YlaYEy45dEj58Kntk4fWVvf4DLeGXCrvxjjIMtika4wGzn0dghkFAoFNR2zTAoHAJBCpimwIBjqrMYvUp3Ff0UYdZMtsCxasLAktFXvPds0ry9+T67ydQil0P85LK8rbCiswDrs69MaBLjHGj2YVqBlzsutLy5+mfJDeQ+ayjGdE51r5zID29qyWm89zsh2DZRKzLNTZBmISk1EgENhghADbx1bMVg5NelHmJFPhNqcCYwF2Lec4r1Y1VSof9wRR9BZjjYP8MJUAsCOfPJpDXnDusDLSo7bVh+E51VpOtsqm+IHVhN5eqJAe5sXew3w/ywG4lb73X+eyDJLOoT+9DClWDrkDa7vgxCmW7kDRigBai0qvj2CMGckiriPBjAKBwHQRzKgPumKUq4rnUfUXAAAc/klEQVQyoxZTKVcBjk/Gq1cPQ+JVi0U65cL+MNX1EJXDbDCv8sKNjqbS6ug5zqC5LNvwqlhbcplFLZI3zQvpYQVIU1a4g471JVBGpBH1yjotOZcc1xefH6KsIJxUxQohzo/L+7l6EqiwIeOYOWJhchPuIIFAYDIIZjQP1aaxnV/PrM+RSUs5ky60rVxTlrjE05rpamlpKrz+xuRnO2mf17KCf+rUW/IgluEoDWANVc2gUDGGGXmMqGRGGueVM6XoseaJ+2YuS5mRMqB9VN7j1AN4PGeSZTmkMlyV9dW0aZ6saDt86CPR/r2oweVnnjdq80iXzWvYGQUCgUkgJiMbWzCsJKy0KVckT1OlK4QumqXSycqVbqEnMNsYM5yWP6klwjlMPvLJU/KKvtViJSzMuD6Xr82lPgRd8kv2oedYa+dlll0zOBp0Da3wuRYzYhkRMyQd2z/k0pIZ7aNzxIge+c78JSxO8jS1wLwnziLaM5ZZsrzTEuWM0bixy4h7TWzTAoHApuNIY0YisgvAH2LILpQA/ByGKB/vA3A2gG8A+KmU0n21frZiWPx44S1XJDYj0VWMU9JbzMjTUvSEimhZ3D5snOu1VSo/MzHR76tatp0V+5jDS/s7c/nmXCpdVNN261o2qrHMxPmH6UkS36tFK6mEx4g+nUu1HVLWU2rGmOaQFu2+Q2tPl5drtw/mkrVpNVbsMaKaGRUnyK2Fwakl+rWOgXkLcvMdXrE7iDUXpJT+dkwfy4ZFeQeAP08pPRfA9wO4DcCVAG5MKZ0L4MZ8HAgEpoYnO//6YM0Fo7DwZCQiOwH8MwBXA0BK6fGU0v0Y0txek5tdA+AnFr1HIBBYJ4xLb11FZS4YhWW2ac/BIBb8IxH5fgC3ALgCwLNSSnfnQd0tIqe2OtoKYCfm6XFJX73tmbqOKKMvVaqlnSDgq+u3Gm08WDGJFI8ZdYBvOAnMb8vYK0S/70nFRSfmGOE7VTir2zB9KM+li60tHt+wtjdh3xfvBbV8cViQ7bl4AL6x4x/mkvdV5euet2VPkrqeH0F5CT9r7/2z3gl+txScBu9R45xnamLBe2d7ZM/NZDGrkxmZc0FKqeXXvgbLbNO2AXghgHellF6A4bfs3pKJyOUiskdE9jyyxCACgcACUAF23zbtZP1fzX+XU29LzQVlJ4tiL4C9KSUNzHh9HsA9InJ6ZkWnA7jXujiltBs5Be5pImk75tX15bSqi6EupMqIlBQoMypXLM97geV2lhyvtWhYMtyWdX0PM+IVXL9XaQupdbtyR7vuWVvueIYO6LuH8tT9xdW550cOrb2x0gSWqJd1tTCXwGLM6BjMoHXKjPLw+aE8/tja6uLUHMthOX35TnlRW1pCY8DPHsxBKGt+wGxSYr1TzN671fbosMToV+3vTymdVznvzQWjsDAzSil9E8AdIvI9ueoCAF8E8CEAl+a6SwHcsOg9AoHAOkG1aT1/ra78uWAUlrUz+vcA3isiR2OIf/azGCa494vIZQBuB/Dq3s50BtfV64HinLKOs3KpzEjFJLrAWgaMLS+GMfGDPWdbqx9deHjFHcOMxohWduby+HyD7fL3AIB/nH62uDr3fFwWrhynT1mtAR+l47LuIJU1bKOSv8kzqAT0Gz0oQ1oYZTOqcvdkPOU5lvuoHI/V3eXnFgu23g/PkJZ/45Ite8HOajHSjzLqLNTidJtjXb2dkTUXjMJSk1FK6XMALPp2wTL9BgKBDcAKJ6PKXNCNSVhgK2PkmGGlnOTZuWRGpPITa7/ZK+oo0XKqrTGjlgOuJYNg1qarfk3jwuIXlVscQ/XKnPbIHx2+VlnUq5OqqJgJMaco65gZ1RiSx4z0eBjdJ+Q/zN2FI8S26oF5JuLZaVrg382T6VjnPFZthYvxGFHN2JdlRp5cxZpX+D0pdxoRzygQCEwHR5I7yKrwJAZGoBb8GvXhmUUbZUSaP00ZEe+He1Jp8f7bWm14leSoqT0yI0XNhL+2Kpb11iLWkmfx2Mu6v5R/CwD4/cOeOjV/EI8R9TAjPh7W6V+T7wVgy8+838v7PQFf89XDMPX5sIKv5mht5UUrUQtLo/16qfEsNuXlKKglk7DY2WFEdpBAIDAVTIwYTWMyOoiBFSkjUo3Ic4o2bE/kMaJa9poe8xjPn7MnqD+vhp59k8WMdPX3tECWfKTHHobB4/4+OQHA/Pf6dEodvbXxUpE1x7Ux6zk29G5FyAXmn/0iUW49VuVZW1tta+D3wAtDYjFZz46px4nXGv8EnfanMRkFAoGNx8Tk19OYjA5iYEWqPWP5EDDTDPEKyIzIsrD1Vlhr5fDkBrzC1jQunj+RtRJxG2VEbHdUMiNLm2SNzRqj5Y9nje2sgtH0BImz+lg1ahFMel3grPRX3AeHMrbsfpYBMy79za0owa28ll4WqfKcNeZgRoFAYDIIZhQIBDYdhzC5TEXTmIxUgK20UelzKWxUWuqlptetS/mAWdDrBR+0XC08mu/R5fI+XriHmpCTvxe7NxzAPDjutzd2K6GsZ8zZA89UwVppPeG6tdXyhNu8zWXhPjAv+Gd3Gt6SAe1ssFa4DitzcTk2Rg/7YKfbWoo8fndrz9GKNjl2bBuJSUxGgUBgYxEyIwdPYjBV55W9BAcu8wzhelYIZg5W3GweS9WAzLnfk3RswXOS1FK/d7mKKSPaScccOtp6ni2m1xPwS1ETYDNb8tT1tRhuLRZcCvA5EoqO33OvAWZKipZxbHkfz3zDey9qBq4MS4DNbInfQxa2W/33GuNuNiYxGQUCgY3FBF3TpjEZHcKwsj2T6ssViYVtnuFfuQrwKqLMgROXWlloufQYQ/mD6hhb8pFDxjn+PswKelLSH++UJfNryZXGuLkoagaMnrsLM7+yjuU/bLahY7TYGzMxZhCWTJHP1dxOPHMRT5Zk5chrGamWLMdj9U/QcXkfTxZWYoLeINOYjAKBwMYjtmkGBGtnb8u1Q+FpbiyZDjMIZgrMLMprPK1ZT7jPQ1TymMsVqbUaWwGyPPlBj6Gfp3njPmqhURjWs2Cmx6yn5r7Dwew9tlGyRXYyZVkLM5ryGi75vlbCBv5e7HJRY229Tr1lv0/Qsedsa/Vv/W4hwA4EApNByIwc9MolPBmRtT/2ZCjMkMoV1gtf6zEyKw87a8a8rLHAvM0MX2PJqlpMBZXznn0Rsy0rsUHvb1Rer9+H++jJoqrX8jOwxtMbfMzS+HlatBoz0pKZS482kt+hnsBvnr2R5QDcIwsKZhQIBCaBmIwa0FXFknVwcHtetVhDVl5/OLUPHVt2ON7q7wVms+xIPG2QxYy8trXV2QO/XNqHtWqyjKXmhOqt8uykWY6Rrc5Zjlaz/G6F7+2BZbNT3r/87DEki8V51uE94WhrYylhaYRZNsQyo/L30nOezA2YpjZtmSSOgUDgKYwVZbcGAIjIVhH5rIh8eNHxTIoZBQKBjcE6bNOuAHAbZo4BozGJyUgpo9LT03L5rKKNUs5v55KTePckKPVU3z1ZaBWey0JZ15vlovzsqfh17D1xk7yYylY8nhqFL/ssP/MWpBaB0TPw6zH8qxmJemPsRXmNZxjZ477juf7U3Gr4N/Seo3UNC7LZ5MPKXNtKdL+qyUhEzgTwYwB+A8AvLdrPJCajQCCwsVixO8jbAbwRdnCEbkxiMmKjR419/bKi7iu5/HIuv5lLFjaWLKcl9GUVcnkNCyg9tmNlNeVz7MRpRaPksia89TLXMjMbI/xmNlILY8H9MZss4bmD1MJzeMJ0FoKXaK3yFgvxBNg1ITSPhQXMtewx3KYHbLLCz8Bi8NaztTCCGZ0sInuK490ppd0AICKvBHBvSukWEXlJf5fzmMRkFAgENhYjtWn7U0pettjzAbxKRC7CIPnYKSJ/klL6mbFjmsRkdBQGOZHGwNZMIC9//azN8e8cSs0gosGzePa3HCE9+QizkLKNxzY4RrWyHauOc8T3XMOq2x6n0JbMqGZg6IUwKc0BvAwVoHqr35b6vPy9vGBqnvyuJm/isVnsrcWIxphTeKyxZyvkyZKsOpYZ1ZhsjfmsSoCdUroKwFUAkJnRrywyEQETmYwCgcDGI9xBDBwF4NmYzeiHM8O/c9ZGmcMduVRmxKu+ZeCnbdTIkQNu1fbwHstg2VFZ58mIatdwqBDW9NVWaWZxnsGfBf5+HKIC8LPqjnF98BxMrXxwHkNqhd6oQcdqBXNjJuGFAynRCqpmoSXDWyT0r5aWFrkmO1oPC+yU0icAfGLR6ycxGQUCgY1HuIMYOApDrjS1K7orl79ZtOEg/Z6DrBWETFcNdgYdw4w87ZoVeMsr2eWjhBfszAqq7gWDt8JJMDxZg+feUI6JQ3zUVnrWXnnMyGIqHiOq2R9xiF8vv1gJT2bEsO7Ty4gsbVorHHEPo61pGPk3tSadKbqDTGIyCgQCG4twlK1gC2YZZLU8rTivbEnti1RmpPIfDi1b1nmhY2vMyAs27zGl8rNnpW1po7wUNJy+piZPsMbioWWnYjEjrvPCZfSMsScgfyvQ3BjBq2fBXPa7SLomD7VMr63UQbX++B2tpSpi5u9NOiHADgQCm44pMqOlvPZF5A0i8gURuVVErhWRY0XkHBG5WUS+IiLvExGOyBkIBDYZ6g6yKq/9VWBhZiQiZwD4BQDPSyk9IiLvB/AaABcBeFtK6ToR+QMAlwF4V62vQxjU3ZzlonSU1a2WZhBRA8kHcqkGhmW2Cd3CebGGe9DrZgDMU3SOV2xt07ytTo9wvccIkOu9OMy8Taxt07ytZIlW7B7Oo1bWear8RbZnNUNGK/vHomDTAev3YxOCnu1hLWZUWW9t09iFhHFEMSMMk9lxIrINgxnP3QBeCuD6fP4aAD+x5D0CgcCKodq0nr+NwsLMKKV0p4j8DoDbATwC4GMAbgFwf0rpYG62F8AZrb4OYggNojO1qpAtQZ+uwsqe2PmwVJuXkSIBf3VZheDSAs/0NRW4d00NLSZkqaE9ZlRzvPRMBzyBNuB/j5oAu+VMO8b4kQW/lvC4Jfgfwxy8d8piYsx2ehQA3jnL9YPNRKzvd0TJjETkBAAXAzgHgwH1DgCvMJom5/rLRWSPiOyZmr1DIPB0wJOdfxuFZbRpLwPw9ZTSPgAQkQ8C+GEAu0RkW2ZHZ2KmlV+DHIJgNwAcJ5K+hZncRwOoWbIVzy3D2hczs2qp+Mu6Fmryi9YPWLvWCzNhMRW+hpmfxXI8Y0RPLlTWcRse65i4z5YciOVILbcQC57shlXhZZ2iZ1H07s0MjB2Py3PsblJTwbdYvPW9+HfiuO3ANNNbLyMzuh3Ai0Vku4gIgAsAfBHATQB+Mre5FMANyw0xEAisB44YZpRSullErgfwGQxin89iYDr/E8B1IvLrue7qVl+qTVOWU1shPFZggeUEPVk/+X6t85amis8twqKYBdTG0zJCLL8nuwr0OHx6shovM0aJlvanFqp2FcaOzLJ6DA57GJIn66tpyPhcz3081snHliGt7gSs8LNTZEZLGT2mlN4C4C1U/TUAL1qm30AgsL5IsLdvm4lJWWCrXZCVN621olk2I/r5AWrr5ZcH5m1nWvtYa2XnY8+1xEIPLfZWSU8WYWnT+HnWNFXMKLmPHvbGsL4na+taDMla2T1NVY/lLT+3mh1SS3vWw4y882MYu9UX24p5DOyIYkaBQOCpiSmq9icxGQmG1UAfjtoQnVK04RAiXnjY0gLbC2rGQdDuxwyLhFhVtORZy/z4NU2VJyuyVr5lVkOWu9QY0RhGyf23nGprMkWW2/F9agzJG3MtbC/ffxFbIUaNlXp9lcfM+K3xx2QUCAQmg9imObBm9lJmpEH6lSHxisHW2+VnLyB+LamiJ6eorc6rQA+b8piQ54tkrfit0LSW5bBnz7RezIjbeoyp/Fyzz2LwO+Qxl7IPlr/wb8H3HxMorYae34nrajKjCK4WCAQmgSlu05Z1lA0EAk9RrMroUUTOEpGbROS2HFLoikXGMxlm1DKF50iO26mN5fSon1mQzbGoa9s0z0WhtlVg1LZenvq/5gLhbRV522QJmJnWewLf0gaFXQ5YebAM3a+97KsIfdEjNNbvp9/Ze0Yl+P1Y71W9xxFXUTPtUKzY6PEggF9OKX1GRI4HcIuIfDyl9MUxnUxmMgoEAhuLVW3TUkp3YwgfhJTSARG5DUO0jqfeZDRGmNZiAxab0pJZj5Xhw4tf3ZPl1DtXu8Yz8GOnyhJehlVGjxFdz4rOjKEVeqOGjZJT9Ah8+Rl7avSawa2XWWSVAfxqsNxFmBFZ32u9ZEYicjaAFwC4eey1k5iMAoHAxmKkNu1kEdlTHO/OUTfWQESeAeADAH4xpfQgn29hEpMRh521ViIvs6sXcqM8p/AMGEv5CMtMel0UrHNeXzWfoB7W46nW2TG4xhZbK/cYg8ll3Vs8tNhGD6urqfo9w8geVsNst8cwktEjs2nJilo534B2oLsO7E8pnVdrICJHYZiI3ptS+mB/1zNMYjIKBAIbi1Vu03IIoasB3JZS+t1F+5nEZHQIg8ZL3T/Y5UPbAPMrkGeEVqLlplFzNuRreY9eYzm8atXYhhc+g5lg2R87/LZysJXXeqEpxrCcGotrsaielZxZ75hre+o9Q0nGmIyyPdl8W+hxKfGCuY3pb4Uyo/MBvA7A50Xkc7nuV1NKHxnTySQmo0AgsLFYpWo/pfTXGFxMl8IkJqMtGNw81NVDZUdHUxtg3u1Dy2WC6lsrhKcB8+RAZRuP1Vgsp9VG+yzDnLAjpMeExgQ9U1h2TS0bqx5m5GlBS3iOvmMC1tdkiN41DGZgPY6rikVkR4qai46XYEDHUzKy3oy1U7PAnsRkFAgENhaHEL5pJrZgYENsXV2yAV75eOW1VqRVhMvgY0+rVn5u2TPVLL4Vnjyo/MxBtLxAYjWW0LKJKutatlYWm/LA9lSAHyp2TErimnxubB/esXUfDji3TOhaC16CUMsCvPcZBDMKBAKbjiMuBnYgEHjqIpiRAcFAb2tbEy//Vk9OrUUiLraiDbLbSFnnbcusLRnHbmLhNG/Fys8tNX0tH5e3DbWcX1tbVM/Z17pPbWvirdT6jMds17z7LfIPaG3TelnFGLlMT9xszoVmbc35f8H6zlMMITKJySgQCGwsIriaA8Ew89fyfXmrfM3oq7Vy14S1vcyoFEZzHG4WqivbOb6oO5rOeQaMliPkImErWiYLFuNjcwYv8mONlY4xXPTghfjouZ+ivLbX0djCIqzC++fvMVnQ8jGn3spcW2Os5fmpYBKTUSAQ2FiEAHsF8FaRnrzyipp6viUf6VHTsxyolqeN61g2sAgLqjEW73t5oVNq19RCozCWYUQ96GVgq/oHXOZ7tLZHNRbMBq1W3HMv7A0jmFEgENh0BDNaAq2VyNIoeTIiXuEtjVivZqzUcnFuN48Rldd4hos9ebh6GV+5EntB47wwu7V+F1lZl1mNme1Ysr4WE1o2POwqw8v2OHkztA2/f2VfzA6DGQUCgckitGmdYJN3wF+RxszuLNuw5CPMfDzNmLquHFPUMSPyHFl7bIZqaAXvX8QWqsamxsiGlsHYfnsCpdWcXj0mMob9LCI7YvbrOUf3XGsxdJYjWduxsDMKBAKTQExGDbAWqrTH4eDpXqiGWhhYZQg11sPMSEuP9ZSaMWVJnu2QZS3bG9jL+l6ebVDNFsqTDfVoGNcD5ardygZbCyniZbvlZ1XeYxENH2u1ONXTmD4YFiv1rq1lje2xwAZCgB0IBCaAYEaBQGAyCGbkYAtm25cTc3lacV5n8QO55JjY1pbEyxyr5UPUjq8HfHW9tU1ruXL0CKn5Bamp2nvV9Ks2YFwGPTF8akLnEtYWz9uWWVuy3uiQNXORWls+7hWMWwoV7oO3aaUAW9/J2lgT6vHbNwPN5yMi7xaRe0Xk1qLuRBH5uIh8JZcn5HoRkd8Tka+KyN+JyAvXc/CBQGAxqNFjz99GoYcZvQfA7wP446LuSgA3ppTeKiJX5uM3AXgFgHPz3w8CeFcuq1CbB53dz8pleeHeXP6/XD6QS2U1ypgewgwPUx0LrrW+XFXYhWOMAaPnyqFgoWcJz0F3jONqDzOqxfAu778sVqE29wTMY8ZYE1Iv46LSuqbHlKCVEw2YZ/6eo2z5HvbGh5+azKj5bqSUPgng21R9MYBr8udrAPxEUf/HacCnAOwSkdNXNdhAILAaqAC752+jsKjM6FkppbsBIKV0t4icmuvPAHBH0W5vrru71tkhDMxG2c2uXJ7/ilmbT/3ZUCozUlajOXS/lcsDs0sOf+bss8pqTqRjYCa38tT0tcwbvVlaa+YHHuuxnF1brhy1rLfLyIbW21jQC/y26pAerb7WgzGVYFMP/f1KluM5avP5Evrb6v/ICjLKbghWLcC2cicls6HI5QAuB9bPkzsQCNg4ktxB7hGR0zMrOh3Avbl+L2YiHwA4E8BdVgcppd0AdgOAiOy7C3joLmA/APyv3OayP1twdA0caDdp4WTksT4FEGNdH2zqWEe+wzrW79KKQ8BHHxrqe7Ah33PRyehDAC4F8NZc3lDU/7yIXIdB/vyAbudqSCmdIiJ7UkrnLTieDUWMdX0QY10fWGNNKV24WePx0JyMRORaAC8BcLKI7AXwFgyT0PtF5DIAtwN4dW7+EQAXAfgqBmXWz67DmAOBwBGI5mSUUrrEOXWB0TYBeP2ygwoEAk8/rDJW1LLYvdkDGIEY6/ogxro+eEqMVQYyEwgEApuLKTGjQCDwNMYkJiMRuVBEvpx92q7c7PGUEJGzROQmEblNRL4gIlfketM/b7MhIltF5LMi8uF8fI6I3JzH+T4RWSQx67pARHaJyPUi8qX8fH9ows/1Dfn3v1VErhWRY6fybI8U/9FNn4xEZCuAd2Lwa3segEtE5HmbO6o1OAjgl1NK3wvgxQBen8en/nnnArgxH08BVwC4rTj+LQBvy+O8D8BlmzIqG+8A8OcppecC+H4M457ccxWRMwD8AoDzUkrPx2Cn+xpM59m+BwCr6r3nWPqPXo7Bf3QaSClt6h+AHwLw0eL4KgBXbfa4KuO9AcDLAXwZwOm57nQAX57A2M7E8OK9FMCHMVjE7wewzXrWmzzWnQC+jiy3LOqn+FzVzelEDBroDwP40Sk9WwBnA7i19RwB/BcAl1jtNvtv05kRfH+2yUFEzgbwAgA3g/zzAJzqX7lheDuAN2LmdnQSgPtTSgfz8ZSe7XMA7APwR3lb+YcisgMTfK4ppTsB/A4Gm7q7MQSNuAXTfbaA/xwn+/82hcmo259tMyEizwDwAQC/mFJ6sNV+oyEirwRwb0rplrLaaDqVZ7sNwAsBvCul9AIMfp2bviWzkOUtFwM4B8CzMSSHeYXRdCrPtobJvhNTmIy6/dk2CyJyFIaJ6L0ppQ/m6ns0PAr5520WzgfwKhH5BoDrMGzV3o4hjIsat07p2e4FsDeldHM+vh7D5DS15woALwPw9ZTSvpTSEwA+COCHMd1nC/jPcbL/b1OYjD4N4NysmTgag2DwQ5s8psMQEQFwNYDbUkq/W5xS/zxgrX/epiCldFVK6cyU0tkYnuFfpJReC+AmAD+Zm236OBUppW8CuENEvidXXQDgi5jYc824HcCLRWR7fh90rJN8thnec/wQgH+TtWovRqf/6IZgs4VWWYh2EYC/xxCu6M2bPR4a249goLF/B+Bz+e8iDPKYGwF8JZcnbvZYizG/BMCH8+fnAPg/GPwF/zuAYzZ7fMU4fwDAnvxs/weAE6b6XAH8GoAvAbgVwH/FEPJqEs8WwLUYZFlPYGA+l3nPEcM27Z35f+3zGDSEm/58U0phgR0IBKaBKWzTAoFAICajQCAwDcRkFAgEJoGYjAKBwCQQk1EgEJgEYjIKBAKTQExGgUBgEojJKBAITAL/H6zFusE4wHFcAAAAAElFTkSuQmCC\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": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEaCAYAAADDgSq4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8FfW9//HXJwkQ1iCbQBIEAVEREBIUhbq0tu6IQFyqtqVWK1fa623v7a+Lv+tta2+XW/3ZVrGu19qqlLC41Vq1rVqxKAm7IAq4JOz7vmT5/P44E3tIEzhJzmROct7Px2Me5MzMmfM5JyHvzMx3MXdHREQkDBlRFyAiIq2XQkZEREKjkBERkdAoZEREJDQKGRERCY1CRkREQqOQERGR0ChkRAAz+7yZlZjZXjPbYGZ/NLNxUdcl0tIpZCTtmdk3gHuA/waOB/oB04EroqyrhpllRV2DSGMpZCStmVkO8APgVnef4+773L3C3Z9z9/8ws3Zmdo+ZrQ+We8ysXfDc88ys3My+aWabgzOgKcG2MWa20cwy417rSjNbGnydYWbfNrM1ZrbNzGaaWbdgW38zczO70cw+Bv4SrP+CmX0U7P9/zexDM7ugAcf7opl9bGZbzex7cXVlmtl3g+fuMbNSM8sPtp1sZi+b2XYzW2VmVzXDt0VaEYWMpLuzgGxgbj3bvweMAU4HRgBnALfHbe8N5AC5wI3AfWZ2nLvPB/YBn47b9/PAk8HXXwcmAOcCfYEdwH21Xvtc4BTgQjM7ldjZ1XVAn7jXrJHI8cYBQ4DPAP9pZqcE678BXAtcAnQBvgzsN7OOwMtBzb2Cfaab2dB6PiuRf+buWrSk7ULsl/bGo2xfA1wS9/hC4MPg6/OAA0BW3PbNwJjg6zuBR4OvOxMLnROCxyuBz8Q9rw9QAWQB/QEHTozb/p/AU3GPOwCHgQsacLy8uO1vA9cEX68CrqjjvV8N/K3WugeAO6L+vmlpOYuu9Uq62wb0MLMsd6+sY3tf4KO4xx8F6z55fq3n7Qc6BV8/CbxpZlOBicBCd6851gnAXDOrjntuFbF7QjXKatXxyWN3329m2+K2J3K8jfXUmU8sTGs7ATjTzHbGrcsCflvHviJ10uUySXd/Bw4Su9RUl/XEftnW6BesOyZ3X0EslC7myEtlEAuMi929a9yS7e7r4g8R9/UGIK/mgZm1B7o38Hj1KQMG1rP+tVrH7OTuUxM4pgigkJE05+67iF2Kus/MJphZBzNrY2YXm9nPgKeA282sp5n1CPb9XQNe4kli90vOAYrj1v8a+JGZnQAQHP9ordlmAZeb2dlm1hb4PmBNOF68h4EfmtlgixluZt2B54GTzOyG4DNpY2aj4+7liByTQkbSnrvfTezm9+3AFmJ/wU8DniZ2X6UEWAosAxYG6xL1FLF7N39x961x638BPAu8ZGZ7gPnAmUep8R3ga8AMYmc1e4jd/znUmOPVcjcwE3gJ2A08ArR39z3A54BriJ29bQR+CrRL8LgimLsmLRNpacysE7ATGOzuH0Rdj0h9dCYj0kKY2eXB5byOwM+JnVl9GG1VIkenkBFpOa4gdtlqPTCYWBNkXYqQlKbLZSIiEhqdyYiISGgUMiIiEpq07/Hfo0cP79+/f9RliIi0KKWlpVvdveex9kv7kOnfvz8lJSVRlyEi0qKY2UfH3kuXy0REJEQKGRERCY1CRkREQqOQERGR0ChkREQkNAoZEREJjUJGRERCo5AREZHQKGRERCQ0ChkREQmNQkZEREKjkBERkdAoZEREJDQKGRERCY1CRkREQqOQERGR0ChkREQkNAoZEREJjUJGRERCo5AREZHQKGRERCQ0ChkREQmNQkZEREKjkBERkdAoZEREJDRZURcQtQ+27uMLj77dqOd2yc7i8hF9+fTJvWiTqbwWEakt7UOmqtrZfaCiUc9dsX43zy/dQPeObblyZC5FhfkM6d05yRWKiLRc5u5R1xCpwsJCLykpadRzK6uqee29LRSXlPPKyk1UVjsj8nKYXJjP+BF9yWnfJsnVioikBjMrdffCY+6nkGl8yMTbtvcQTy9eT3FJGe9u3EPbrAwuGtqbosI8xg7sQUaGJaFaEZHUoJBJULJCpoa7s3zdbopLy3h60Tp2H6wkt2t7Jo3KZXJBPv26d0jaa4mIREUhk6Bkh0y8gxVVvLxiEzNLynhj9VbcYcyJ3SgqyOfiYb3p0Dbtb4mJSAulkElQmCETb/3OA8xZWE5xaTkfbdtPp3ZZXDa8D0WF+Yzq1xUzXU4TkZZDIZOg5gqZGu7O2x9sZ2ZJOS8s28CBiioG9uxIUWE+E0fm0qtLdrPVIiLSWAqZBDV3yMTbe6iSF5ZuYGZJGSUf7SAzwzj3pJ5cVZjHp08+nrZZ6nsjIqlJIZOgKEMm3toteykuLWfOwnI27T5Et45tmXB6LkWFeZzSp0vU5YmIHEEhk6BUCZkalVXV/O39rRSXlvHyik1UVDnDcnMoKsxj/Ii+dO3QNuoSRUQUMolKtZCJt33fYZ5ZvI7iknJWbNhN28wMPjf0eIoK8xk3qAeZ6nsjIhFRyCQolUMm3vJ1u5hVWs7Ti9exc38FfXKymTQqj8kFefTv0THq8kQkzShkEtRSQqbGocoqXlmxmeLSMl5/bwvVDmcM6EZRQR6XDOtDx3bqeyMi4VPIJKilhUy8jbsOMnthOcUlZXy4bT8d22Zy6fA+XFWYT8EJx6nvjYiERiGToJYcMjXcnZKPdlBcUsbzSzew/3AVJ/boyOTCPCaNyuN49b0RkSRTyCSoNYRMvH2HKnlh2QaKS8p5+8PtZBice1JPigrz+cwpvWiXlRl1iSLSCqR1yJhZR2A6cBh41d2fqG/f1hYy8T7Yuo9ZpWXMLl3Hxt0H6dqhzSd9b4b2zYm6PBFpwVIiZMzs34CvAA4sA6a4+8FGHOdR4DJgs7ufVmvbRcAvgEzgYXf/iZndAOx09+fM7PfufnV9x27NIVOjqtp5Y/VWZpaU8fI7mzhcVc3Qvl0oKsjjitNzOa6j+t6ISMNEHjJmlgu8AZzq7gfMbCbwgrs/FrdPL+CAu++JWzfI3VfXOtY5wF7g8fiQMbNM4D3gs0A5sAC4FrgC+KO7LzazJ9398/XVmQ4hE2/n/sM8s3g9xaVlLF8X63vz2VOPZ3JhHucM7qm+NyKSkERDJuz2rllAezOrADoA62ttPxeYamaXuPtBM7sJuBK4JH4nd3/dzPrXcfwzgNXuvhbAzGYQC5hyIA9YDGgAsDhdO7Tli2f354tn92fF+n/Me/OHZRvo3SWbiaNi00gPUN8bEUmC0ELG3deZ2c+Bj4EDwEvu/lKtfYrNbAAww8yKgS8TOytJVC5QFve4HDgT+CVwr5ldCjxX1xPN7HLg8kGDBjXg5VqXU/t24Y6+Q/nOxafw55WbKC4t59evrWH6q2sY3f84igrzuVR9b0SkCcK8XHYcMBu4GtgJFAOz3P13dew7g9jZy0B331LP8foDz9e6XFYEXOjuXwke3wCc4e5fS7TOdLtcdiybdh9kzsJ1FJeUsXbrPjq0zeSSYbG+N6P7q++NiMSkwuWyC4APakLDzOYAZwNHhIyZfQo4DZgL3AFMa8BrlAP5cY/z+OdLctIAx3fJZup5A7nl3BNZ+PEOikvKeW7JemaVltO/ewcmF+QxqSCPPjntoy5VRFqAMO9XfAyMMbMOFvvz9zPAyvgdzGwk8BCx+yhTgG5mdmcDXmMBMNjMBphZW+Aa4NmkVJ/mzIyCE7rxk0nDWXD7BdxVNILju2Tz85feY+xP/sIXHn2b55eu51BlVdSlikgKC7sJ8/eJXS6rBBYBX3H3Q3HbxwK73X1Z8LgN8CV3f6jWcZ4CzgN6AJuAO9z9kWDbJcA9xJowP+ruP2pIjbpc1jAfbdvH7NJyZpWWs37XQXLat2HC6X0pKsxnaN8uupwmkiYib8LcUihkGqeq2nlzzVZmlpTzp3c2criympN7d+aqwnwmjMylm/reiLRqCpkEKWSabtf+Cp5dso7i0nKWlu+iTaZxwSnHUxT0vcnKVCtykdZGIZMghUxyvbtxN8Ul5Ty9aB3b9h2mV+d2TByVR1FhHgN7doq6PBFJEoVMghQy4ThcWc1f3t3MrNIy/rpqC1XVTsEJx1FUkMelw/vQObtN1CWKSBMoZBKkkAnf5j0Hmbswdjlt9ea9tG+TycXDenNVYT5nDuimxgIiLZBCJkEKmebj7iwq2/lJ35u9hyrp160DRUHfm75d1fdGpKVQyCRIIRONA4erePGdDcxcUM7f127DDMYN6kFRYT6fO/V4stto3huRVKaQSZBCJnpl2/czK+h7s27nAbpkZ3FFMO/NsNwcXU4TSUEKmQQpZFJHdbXz97XbmFlSxovLN3Io6HszuSCPK0fm0r1Tu6hLFJGAQiZBCpnUtOtABc8vXc/MknKWlO0kK8P4zCm9KCrI57wh6nsjEjWFTIIUMqnvvU17KC4pY+6idWzde5iendsxcWTsctqgXp2jLk8kLSUlZIKZJye4++xkFpdKFDItR0VVNX99dzPFpeX89d3NVFY7I/t1paggn8tG9KGL+t6INJukncmY2d/c/VNJqyzFKGRapi17DvH0onXMLCnj/c17yW6TwcWn9aGoMI8xA7qToWmkRUKVzJC5HdgL/B7YV7Pe3Xc3tchUoJBp2dydJeW7KC4p49kl69lzsJL8bu2ZPCqfSQW55B3XIeoSRVqlZIZMWR2r3d37Nba4VKKQaT0OVlTxp3c2UlxSzrw1WwE4e2B3rirM58KhvdX3RiSJdOM/QQqZ1ql8x35ml66juLSM8h0H6JydxfgRsXlvRuSp741IUyXzTCYLuBk4J1j1KvCwu1c2tchUoJBp3aqrnfkfbKO4pJw/Lt/AwYpqTjq+E0UFsXlvenZW3xuRxkhmyDwAdAQeD1ZdDxx095ubXGUKUMikj90HK/jD0g3MLClj0cexvjfnn9yLooI8zj+5F23U90YkYckMmSXuPuJY61oqhUx6Wr15D8Ul5cxeuI6tew/Ro1NbrhyZS1FhPicdr743IseSzJBZBFzp7h8Gj/sDc919ZNPLjJ5CJr1VVFXz2qotFJeW8eeVsb43I/K7UlSQx+Uj+pLTXn1vROqSzJD5LPAI8B5gwCDgRnd/JRmFRk0hIzW27o31vZlVWs67G/fQLiuDi06LzXtz1onqeyMSL1k9/jOA0cBS4BRiIbPC3Q8kq9CoKWSkNndn2bpdFJeU88zidew+WElu1/ZMLshjckEe+d3U90YkmWcy8919TNIqSzEKGTmagxVVvLRiE8UlZbyxeivusb43RYV5XDS0D+3bqu+NpKdkhswPgRJ3fyZZxaUShYwkat3OA8wpLae4tJyPt++nc7ssLhvRl6LCPEbmd1XfG0kryQyZHUAOcAg4QOySmbt7t2QUGjWFjDRUdbXz9ofbKS4p54VlGzhQUcWgXp0oKsjjylG59OqcHXWJIqFL1j0ZA/KBdbW3uXtVkypMEQoZaYq9hyr5QzDvTelHO8jMMM4f0pPJBfl8+uRetM1S3xtpnZJ5JlPq7gVJqyzFKGQkWdZs2UtxSTlzFpazec8hundsy4Rg3puTe3eJujyRpEpmyNwPPOTuC5NVXCpRyEiyVVZV87f3tzKzpIxXVm6iosoZnpdDUWE+44f3JaeD+t5Iy5fMkFlGrPnyGmJD/dfckxmVjEKjppCRMG3fd/iTeW/e3biHtlkZXDi0N1cV5nH2wB5kqu+NtFDJDJmBda139zWNrC2lKGSkObg776zfTXFJGU8vXs+uAxX0zclmUtD35oTuHaMuUaRBkjrUv5mNAU5y98fNrDvQ0d0/TkKdkVPISHM7WFHFKys3UVxSzuvvb8EdzhzQjasK87l4WG86tM2KukSRY0r2zJhjgYHufpKZ5QK/d/dxySk1WgoZidKGXQeYs3AdxSVlfLhtP53aZXHZ8Ng00qP6Hae+N5Kykhkyi4GRwMKaQTHNbKm7D09KpRFTyEgqcHcWfLiD4pIy/rBsA/sPx/re/HTSMApOaBVd0qSVSTRkEmnEf8hjSeTBgTVwk0iSmRlnDOjG/xSNYMH3LuBnk4dzuLKaqx+Yz2PzPiDdZ7CVliuRkJljZvcBOWY2BXgJeDTcskTSV8d2WVxVmM9zXxvHeUN68l/PreBfZyxm/+FWMRmtpJljhoy7/xR4HngWGAH8yN3vCbswkXSX074ND95QyH9cOITnl65nwn3zWLtlb9RliTRIQq3LWjPdk5GW4I33t/L1GYs4XFnNz4tGcNFpvaMuSdJcMu/JiEjExg3uwXNfG8fAnh255Xel/PiPK6msqo66LJFjUsiItBC5Xdsz85azuO7Mfjzw2lpueORttuw5FHVZIkeVUMiYWVszGxR2MSJydO2yMvnRlcO4q2gECz/ewWW/+hulH22PuiyReh0zZMzsUmAZ8HLw+HQzmxt2YSJSv0kFecz9l7Fkt8lUM2dJaYmcyfwAOBPYCeDuiwGd1YhE7NS+XXh2mpo5S2pLJGQq3H1nrXX6k0kkBaiZs6S6REJmpZldBWSY2QAzuweYH3JdIpKgjAzj1vMH8fiXz2Tr3sOMv3ceLy7fEHVZIkBiITMNKACqgTnAQeC2MIsSkYY7spnzQn78gpo5S/QSGSBzmLsva6Z6mp06Y0prc6iyih88t4In3vqYMSd241fXjqJn53ZRlyWtTDI7Y043s3fM7A4zOzkJtYlIiOKbOS/6eKeaOUukEhm77FPA54A9wG/MbJGZfTv0ykSkSdTMWVJBQp0x3X2du98NfIlYn5kfhlmUiCSHmjlL1BLpjDnYzG4PJi97CFgA9Au9MhFJCjVzliglcibzJHAAGO/u49z9V+6u9pEiLYiaOUtUErknM9rd73L3j5ujIBEJj5o5S3OrN2TM7Kng30VmtjBuWWRmC5uvRBFJpiNGc359Ldc/8pZGc5bQ1NtPxszy3L3czAbWtd3d14RaWTNRPxlJZ7NLy/nu3GV07dCG6deNouCEblGXJC1Ek/vJuHt58OWN7r4mfgFuTFahIhIdNXOWsCVy4/+iOtZdmuxCRCQaauYsYTraPZmvmtkiYEitezLvAyubr0QRCZuaOUtYjnZP5jigO/BjIL6H/x5339wMtTUL3ZMROdIb72/l6zMWcbiymp8XDeei0/pEXZKkoGTck9nh7qvdvSi4D7ODWH+ZLDPrm8RaRSSFqJmzJFMiPf4vMbP3gHLgLaAM+EvYhYlIdGqaOV8/Rs2cpWkSufH/38BYYJW75xNrCPBqmEWJSPTaZWVy54Rh3H3VCBaXaTRnaZxEQqbS3bcQmxnT3P1lYFTIdYlIipg4Ko85U9XMWRonkZDZZWYdgTeAx83sLmKzZIpImqirmfO+Q2rmLMeWSMhMAA4Rm3L5VWAdcHmINYlICqrdzPnK6WrmLMeWyACZe9y90t0r3P0Rd787uHwmImlGozlLQx2tM+YOM9set+yI/7c5ixSR1PJJM+dendTMWY7qaGcyPYCecUuPWv+KSBrL7dqemV8do2bOclRH64xZVbMAo4Hrgq9zgNzmKlBEUpeaOcuxJNIZ83bgDuD2YFV7YrNliogA/9zM+X/VzFkCibQumwxcAuwDcPd1QJcwixKRlie+mfP31cxZAomEzCGP/UniAGbWIdySRKSlUjNnqS2RkJljZvcBOWY2BXgJeDTcskSkpVIzZ4mXSD+ZnwLPA88CI4Afufs9YRcmIi3buME9eF7NnNNe1tE2mlkm8IK7Xwj8sXlKEpHWom/QzPmHz6/ggdfXsqR8J7+6dhQ9O7eLujRpJkc9kwmaLB82M93oF5FGUTPn9JbIPZm9wBIze8DM7q5Zwi5MRFoXNXNOT0e9XBZ4JVhERJqkppnzN2cu5vvPrWDRxzv58cRhdGyXyK8iaYmO+Z1190eaoxARSQ81zZzvf20Nd720inc37ubX1xdwYs9OUZcmIUjkcpmISFKpmXP6UMiISGTUzLn1SzhkzExtDkUk6fpqNOdWLZEBMs8ws2XA+8HjEWb2q9ArE5G0oWbOrVciZzK/BC4DtgG4+xLg/DCLEpH0pGbOrU8iIZPh7h/VWlcVRjEiIhrNuXVJJGTKzOwMwM0s08xuA94LuS4RSWMazbn1SCRkpgLfAPoBm4AxwToRkdComXPrkMgozJvd/Rp37xEs17j71uYoTkSkdjPn/1Yz5xblmD3+zewhggnL4rn7zaFUJCJSS/xozg++vpalGs25xUjkctkrwJ+DZR7QC1AjdhFpVmrm3DIlcrns93HLb4CJwKnhlyYi8s/UzLllacywMgOAE5JdiIhIov7RzLmXmjmnuER6/O8ws+3BshN4Gfhu+KWJiNQv1sy54IhmzmvUzDnlHDVkzMyAEUDPYDnO3U9095nNUZyIyNHUbuZ8xb3z+OMyNXNOJceaftmBue5eFSy68CkiKSe+mfPUJ9TMOZUkck/mbTMbFXolIiJNED+a84MazTll1BsyZlbTh2YcsaBZZWYLzWyRmS1snvJERBKnZs6p52idMd8GRgETmqkWEZGkmDgqj1P6dOGW35Vy9QPz+d6lp/Cls/sTu80szelol8sMwN3X1LU0U30iIo1ySh81c04FRzuT6Wlm36hvo7vfHUI9IiJJU9PM+f7X1nDXS6tYuWE3v76hgIE9O0VdWto42plMJtAJ6FzPIiKS8uKbOW/bp2bOzc3qa5VsZgvdvdW3KissLPSSkpKoyxCRZrB+5wGmPrGQJWU7ufmcE/nWhUPIymzMwCdiZqXuXnis/Y55T0ZEpLWo3cz5uoffYuteNXMO09FC5jPNVoWISDOp3cz5invnsWL97qjLarXqDRl3V+NyEWm1Jo7Ko/iWs6iqdibd/6bu04REFyNFJG0Nz+vKs9PGMqR3Z6Y+sZB7XnmP6mqNnpVMChkRSWu9umQz4+YxTByVyz2vvM+0pxay/7D60ySLQkZE0l52m0zuKhrB9y45hReXb2Ty/X9n3c4DUZfVKihkREQAM+Omc07kkS+Npmz7fsb/6g1KPtSt6aZSyIiIxDl/SC/m3no2nbOzuPah+cxcUBZ1SS2aQkZEpJZBvTrzzK3jGHNid741eyk/eG6F5qdpJIWMiEgdcjq04X+/NJopY/vz6LwPmPLYAnbtr4i6rBZHISMiUo+szAzuuHwoP500jPlrtzFh+jxWb94bdVktikJGROQYrh7djydvGsPuAxVcOX0ef121OeqSWgyFjIhIAkb378azXxtH3nEduPGxBTz0+lrqG2BY/kEhIyKSoNyu7Zk99SwuHNqbH72wkn8vXsrBiqqoy0ppChkRkQbo0DaL+z4/itsuGMzsheVc+9B8Nu8+GHVZKUshIyLSQBkZxm0XnMT9143i3Q17GH/vPJaV74q6rJSkkBERaaSLh/Vh1tSzyMwwih54k+eWrI+6pJSjkBERaYKhfXN4ZtpYhuXm8LWnFvHzP63SSM5xFDIiIk3Uo1M7nvjKGK4uzOfev67mq78rZe8hjeQMChkRkaRom5XBTyYN478uP5W/vLuZSdPfpGz7/qjLipxCRkQkScyML40dwG+mnMGGXQcYf+8bzF+7LeqyIqWQERFJsnGDe/DMtHF069iW6x9+iyfe+ijqkiKjkBERCcGAHh2Ze+tYPjW4B9+bu5z/+/RyKtJwJGeFjIhISLpkt+HhL47mq+ecyG/nf8QXHnmbHfsOR11Ws1LIiIiEKDPD+M4lp3BX0QhKP9rBFffN471Ne6Iuq9koZEREmsGkgjxmfHUMByqquPK+ebyyYlPUJTULhYyISDMZ1e84np02lhN7duKm35Yw/dXVrX4kZ4WMiEgz6pPTnuJbzuKy4X352YuruO33i1v1SM5ZURcgIpJusttk8strTufk3p35+Uur+GDrPh68oZDeOdlRl5Z0OpMREYmAmXHr+YN48IZC1mzey/h732Bx2c6oy0o6hYyISIQ+e+rxzPmXsbRrk8FVD/yduYvKoy4pqRQyIiIRG9K7M8/cOo5R/bryb79fwo9fWElVKxnJWSEjIpICunVsy29vPJPrx/TjgdfXctPjJew5WBF1WU2mkBERSRFtMjO4c8IwfjjhNF5/bwtXTn+TD7fui7qsJlHIiIikmBvGnMDjN57Btr2HuOK+ecxbvTXqkhpNISMikoLOHtiDZ24dR+8u2Xzh0bd5bN4HLbLjpkJGRCRF9evegdn/cjbnD+nFfz23gu/OXcbhypY1krNCRkQkhXVql8WDNxQw7fxBPPV2Gdc//Bbb9h6KuqyEKWRERFJcRobx7xcO4ZfXjmRJ+U7G3zuPFet3R11WQhQyIiItxPgRfSm+5Syqqp3Jv36TF5dviLqkY1LIiIi0IMPzuvLstLGcdHxnbvndQn7xyvsp3SBAISMi0sL06pLNjJvHMHFULv/vlfeY9uQi9h+ujLqsOmkUZhGRFii7TSZ3FY3g5N6d+fEf3+XDbft48AuF5HZtH3VpR2iVZzJm1tHMfmNmD5nZdVHXIyISBjPj5nMG8ugXR/Pxtv1cce8blHy4PeqyjhBayJjZEDNbHLfsNrPbGnmsR81ss5ktr2PbRWa2ysxWm9m3g9UTgVnufhMwvglvQ0Qk5Z1/ci/m3no2ndplce1D85lZUhZ1SZ8ILWTcfZW7n+7upwMFwH5gbvw+ZtbLzDrXWjeojsM9BlxUe6WZZQL3ARcDpwLXmtmpQB5Q8ym33innREQCg3rFRnI+c0B3vjVrKT98fgWVVdF33Gyuy2WfAda4+0e11p8LPGNm2QBmdhPwy9pPdvfXgbrOAc8AVrv7Wnc/DMwArgDKiQUN1PMezexyM3tw165djXk/IiIpJ6dDGx6bMpopY/vzyBsfMOWxBezaH+1Izs0VMtcAT9Ve6e7FwIvAjODeyZeBqxpw3Fz+ccYCsXDJBeYAk8zsfuC5up7o7s+5+805OTkNeDkRkdSWlZnBHZcP5aeThjF/7TaunD6PNVv2RlZP6CFjZm2J3Rcprmu7u/8MOAjcD4x394Z8Glb3IX2fu09x96nu/kSDixYRaeGuHt2PJ28aw64DFUy4bx6vrtocSR3NcSZzMbDQ3Tf3JaRzAAAIUElEQVTVtdHMPgWcRux+zR0NPHY5kB/3OA9Y35giRURam9H9u/HMtLHkHdeBLz+2gIf/trbZO242R8hcSx2XygDMbCTwELH7KFOAbmZ2ZwOOvQAYbGYDgjOma4Bnm1iviEirkXdcB2ZPPYsLh/bmzj+s5D9mLeVQZfO1hwo1ZMysA/BZYvdI6tIBKHL3Ne5eDXwRqN04ADN7Cvg7MMTMys3sRgB3rwSmAX8CVgIz3f2d5L8TEZGWq0PbLO77/Chuu2Aws0rLufbB+Wzec7BZXttSecyb5lBYWOglJSVRlyEi0ixeWLaBb85cQtcObXjoC4Wcltu4xk9mVuruhcfar1X2+BcRkbpdMqwPs6aeRYYZ63YeCP31NHaZiEiaGdo3hz9/81yy22SG/lo6kxERSUPNETCgkBERkRApZEREJDQKGRERCY1CRkREQqOQERGR0ChkREQkNAoZEREJTdoPK2NmW6hjvLQE9QC2JrEckXg5gGbVC1e6f8ZNef8nuHvPY+2U9iHTFGZWksjYPSKNYWYPuvvNUdfRmqX7Z9wc71+Xy0RSV52zukpSpftnHPr715lME+hMRkTk6HQm0zQPRl2AiEgq05mMiIiERmcyIiISGoWMiIiERpOWibQyZtYRmA4cBl519yciLqnVSefPuKHvXWcySWRmHc3sN2b2kJldF3U9Eh0zyzezv5rZSjN7x8z+tQnHetTMNpvZ8jq2XWRmq8xstZl9O1g9EZjl7jcB4xv7uqnOzLLN7G0zWxJ8xt9vwrFa5GdsZplmtsjMnm/CMUJ97wqZY6jvG5DKP3iSEiqBb7r7KcAY4FYzOzV+BzPrZWada60bVMexHgMuqr3SzDKB+4CLgVOBa4PXyAPKgt2qmvg+Utkh4NPuPgI4HbjIzMbE75AGn/G/Aivr2pAq710hc2yPUesb0AJ+8CRi7r7B3RcGX+8h9osgt9Zu5wLPmFk2gJndBPyyjmO9Dmyv42XOAFa7+1p3PwzMAK4Ayon9LEIr/j/uMXuDh22CpXZz2Vb7GZtZHnAp8HA9u6TEe2+1P4DJUs83IGV/8CT1mFl/YCTwVvx6dy8GXgRmBJdXvwxc1YBD5/KPP2og9vOXC8wBJpnZ/bTyHu3B5aLFwGbgZXdPp8/4HuBbQHVdG1PlvevGf+PU9eGfSeyvhHvN7FJa+X9uSYyZdQJmA7e5++7a2939Z2Y2A7gfGBj3l3lCh69jnbv7PmBKowpuYdy9CjjdzLoCc83sNHdfXmufVvcZm9llwGZ3LzWz8+rbLxXeu/7abpx6P3x3n+LuU9OptYnUzczaEAuYJ9x9Tj37fAo4DZgL3NHAlygH8uMe5wHrG1Fqi+fuO4FXqfveQmv8jMcC483sQ2JXUj5tZr+rvVMqvHeFTOOk6g+epAgzM+ARYKW7313PPiOBh4hdap0CdDOzOxvwMguAwWY2wMzaAtcAzzat8pbDzHoGZzCYWXvgAuDdWvu0ys/Y3b/j7nnu3p9YTX9x9+vj90mV966QaZyU/MGTlDIWuIHYX5iLg+WSWvt0AIrcfY27VwNfpI65jczsKeDvwBAzKzezGwHcvRKYBvyJWMOCme7+TnhvKeX0Af5qZkuJ/Z982d1rN+VN5884Jd67xi47huAbcB6xCco2AXe4+yPBL4x7gEzgUXf/UXRVioikJoWMiIiERpfLREQkNAoZEREJjUJGRERCo5AREZHQKGRERCQ0ChkREQmNQkbSkpntDf7tb2afT/Kxv1vr8ZvJPH6ymdmXzOzeqOuQ1kkhI+muP9CgkAmmejiaI0LG3c9uYE0tSgKfh6QxhYyku58AnwqGffm3YOj4/zGzBWa21My+CmBm51lspssngWXBuqfNrNRiszLeHKz7CdA+ON4TwbqasyYLjr3czJaZ2dVxx37VzGaZ2btm9kQw9tkRgn1+arHZIN8LBj/8pzMRM3u+ZmReM9sbPKfUzF4xszOC46w1s/iJ9fLN7EWLTcR3R9yxrg9eb7GZPVATKMFxf2BmbwFnJeubIa2PhvqXdPdt4N/d/TKAICx2uftoM2sHzDOzl4J9zwBOc/cPgsdfdvftweCMC8xstrt/28ymufvpdbzWRGIzOI4gNkzRAjN7Pdg2EhhKbKDVecTGPnujjmNkufsZwbBGdxAbFPJoOhKbh/3/mNlc4E7gs8Qm2/sN/xhz7wxio/XuD+r6A7APuBoY6+4VZjYduA54PDjucnf/z2O8vqQ5hYzIkT4HDDezycHjHGAwcBh4Oy5gAL5uZlcGX+cH+207yrHHAU8Fc6BsMrPXgNHA7uDY5QAWm4SrP3WHTM2UAaXBPsdymNjEVRA7AzsUBMayWs9/2d23Ba8/J6i1EiggFjoA7YlNDgaxmV9nJ/D6kuYUMiJHMuBr7v6nI1bGLj/tq/X4AuAsd99vZq8C2Qkcuz6H4r6uov7/m4fq2KeSIy99x9dR4f8YoLC65vnuXm1m8a9RexBDD+r9jbt/p446DgZhKXJUuicj6W4P0Dnu8Z+AqRabcAwzO8nMOtbxvBxgRxAwJwNj4rZV1Dy/lteBq4P7Pj2Bc4C3k/AePiQ2O2SGmeUTu/TVUJ81s27Bpb8JxC7Z/RmYbGa9AILtJyShXkkjOpORdLcUqDSzJcBjwC+IXUZaGNx830Lsl25tLwK3BHOZrALmx217EFhqZgvd/bq49XOJ3SRfQuxM4VvuvjEIqaaYB3xA7HLYcmBhI47xBvBbYBDwpLuXAJjZ7cBLZpYBVAC3UsecJCL10VD/IiISGl0uExGR0ChkREQkNAoZEREJjUJGRERCo5AREZHQKGRERCQ0ChkREQmNQkZERELz/wGmm3N0WpNcKgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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.5"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}